lib/fpa/Common/SliceBySliceRandomWalker.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5 #ifndef __fpa__Common__SliceBySliceRandomWalker__hxx__
   6 #define __fpa__Common__SliceBySliceRandomWalker__hxx__
   7
   8 #include <vector>
   9 #include <itkExtractImageFilter.h>
  10 #include <itkImageRegionConstIterator.h>
  11 #include <itkImageRegionIterator.h>
  12 #include <itkJoinSeriesImageFilter.h>
  13 #include <itkMinimumMaximumImageCalculator.h>
  14
  15 // -------------------------------------------------------------------------
  16 template< class _TImage, class _TLabels, class _TScalarImage >
  17 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
  18 SliceBySliceRandomWalker( )
  19   : Superclass( )
  20 {
  21   ivqITKInputConfigureMacro( InputLabels, TLabels );
  22   ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
  23 }
  24
  25 // -------------------------------------------------------------------------
  26 template< class _TImage, class _TLabels, class _TScalarImage >
  27 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
  28 ~SliceBySliceRandomWalker( )
  29 {
  30 }
  31
  32 // -------------------------------------------------------------------------
  33 template< class _TImage, class _TLabels, class _TScalarImage >
  34 void
  35 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
  36 GenerateData( )
  37 {
  38   typedef itk::Image< TPixel, TImage::ImageDimension - 1 > _TSliceImage;
  39   typedef itk::Image< TLabel, TImage::ImageDimension - 1 > _TSliceLabels;
  40   typedef itk::Image< TScalar, TImage::ImageDimension - 1 > _TSliceScalarImage;
  41
  42   // Preare I/O objects
  43   const TImage* input = this->GetInput( );
  44   const TLabels* labels = this->GetInputLabels( );
  45   const TScalarImage* vesselness = this->GetInputVesselness( );
  46   this->AllocateOutputs( );
  47
  48   // Some values
  49   typename TImage::RegionType region = input->GetRequestedRegion( );
  50   typename TImage::SizeType regionIndex = region.GetIndex( );
  51   typename TImage::SizeType sliceSize = region.GetSize( );
  52   int numSlices = sliceSize[ 2 ];
  53   sliceSize[ 2 ] = 0;
  54
  55   // Composite image
  56   typename TScalarImage::Pointer composite = TScalarImage::New( );
  57   composite->CopyInformation( vesselness );
  58   composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
  59   composite->Allocate( );
  60   this->_Composite( composite, labels, vesselness );
  61
  62   // Extract slices
  63   std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
  64   std::vector< typename _TSliceLabels::Pointer > binaryTree( numSlices );
  65   std::vector< typename _TSliceScalarImage::Pointer > fusion( numSlices );
  66   for( int i = 0; i < numSlices; ++i )
  67   {
  68     // Prepare extraction region
  69     typename TImage::IndexType sliceIndex = regionIndex;
  70     sliceIndex[ 2 ] += i;
  71     typename TImage::RegionType sliceRegion( sliceIndex, sliceSize );
  72
  73     // Extract regions
  74     typename _TSliceImage::Pointer input_slice;
  75     typename _TSliceLabels::Pointer labels_slice;
  76     typename _TSliceScalarImage::Pointer composite_slice;
  77     this->_Slice( input_slice, input, sliceRegion );
  78     this->_Slice( labels_slice, labels, sliceRegion );
  79     this->_Slice( composite_slice, composite, sliceRegion );
  80
  81     // Update vectors with each image
  82     data3D[ i ] = input_slice;
  83     binaryTree[ i ] = labels_slice;
  84     fusion[ i ] = composite_slice;
  85
  86   } // rof
  87
  88   // Random walker slice-by-slice
  89   this->_GoDown( binaryTree, data3D, fusion, numSlices );
  90   this->_GoUp( binaryTree, data3D, fusion, numSlices );
  91
  92   // Join results
  93   typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
  94   typename _TJoin::Pointer join = _TJoin::New( );
  95   for( int i = 0; i < numSlices; ++i )
  96     join->SetInput( i, binaryTree[ i ] );
  97   join->Update( );
  98   this->GetOutput( )->Graft( join->GetOutput( ) );
  99 }
 100
 101 // -------------------------------------------------------------------------
 102 template< class _TImage, class _TLabels, class _TScalarImage >
 103 void
 104 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 105 _Composite(
 106   typename TScalarImage::Pointer& composite,
 107   const TLabels* labels,
 108   const TScalarImage* vesselness
 109   )
 110 {
 111   // Min-Max vesselness values
 112   typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
 113   typename _TMinMax::Pointer minMax = _TMinMax::New( );
 114   minMax->SetImage( vesselness );
 115   minMax->Compute( );
 116   TScalar maxVess = minMax->GetMaximum( );
 117   typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
 118
 119   // Composite image
 120   typename TScalarImage::Pointer composite = TScalarImage::New( );
 121   composite->CopyInformation( vesselness );
 122   composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
 123   composite->Allocate( );
 124   itk::ImageRegionConstIterator< TScalarImage > v( vesselness, region );
 125   itk::ImageRegionConstIterator< TLabels > l( labels, region );
 126   itk::ImageRegionIterator< TScalarImage > c( composite, region );
 127   v.GoToBegin( );
 128   l.GoToBegin( );
 129   c.GoToBegin( );
 130   for( ; !v.IsAtEnd( ); ++v, ++l, ++c )
 131     c.Set( ( l.Get( ) != 0 )? maxVess: v.Get( ) );
 132 }
 133
 134 // -------------------------------------------------------------------------
 135 template< class _TImage, class _TLabels, class _TScalarImage >
 136 template< class _TSlicePtr, class _TInput >
 137 void
 138 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 139 _Slice(
 140   _TSlicePtr& slice,
 141   const _TInput* input,
 142   typename _TInput::RegionType region
 143   )
 144 {
 145   typedef typename _TSlicePTr::ObjectType _TSlice;
 146   typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
 147
 148   typename _TExtract::Pointer extract = _TExtract::New( );
 149   extract->SetInput( input );
 150   extract->SetExtractionRegion( region );
 151   extract->SetDirectionCollapseToIdentity( );
 152   extract->Update( );
 153   slice = extract->GetOutput( );
 154   slice->DisconnectPipeline( );
 155 }
 156
 157 // -------------------------------------------------------------------------
 158 template< class _TImage, class _TLabels, class _TScalarImage >
 159 template< class _TBinaryTree, class _TData3D, class _TFusion >
 160 void
 161 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
 162 _GoDown(
 163   _TBinaryTree& binaryTree,
 164   const _TData3D& data3D,
 165   const _TFusion& fusion,
 166   int numSlices
 167   )
 168 {
 169   typedef typename _TBinaryTree::value_type     _TSliceBinaryPtr;
 170   typedef typename _TData3D::value_type         _TSliceData3DPtr;
 171   typedef typename _TFusion::value_type         _TSliceFusionPtr;
 172   typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
 173   typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
 174   typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
 175
 176   int z = -1;
 177   while( z < numSlices - 2 )
 178   {
 179     z++;
 180     _TSliceBinaryPtr tmp = binaryTree[ z ];
 181     _TSliceBinaryPtr next = binaryTree[ z  + 1];
 182     _TSliceData3DPtr data = data3D[ z + 1 ];
 183     _TSliceFusionPtr vess = fusion[ z + 1 ];
 184
 185     // Create seeds
 186     _TSliceBinaryPtr seeds = _TSliceBinary::New( );
 187     seeds->CopyInformation( tmp );
 188     seeds->SetRequestedRegion( tmp->GetRequestedRegion( ) );
 189     seeds->Allocate( );
 190     seeds->FillBuffer( 0 );
 191
 192     /* TODO
 193
 194
 195        SliceType::RegionType      sliceRegion = tmp->GetLargestPossibleRegion();
 196        SliceType::SizeType                size = sliceRegion.GetSize();
 197        SliceType::SpacingType     spacing = tmp->GetSpacing();
 198        SliceType::PointType       origin = tmp->GetOrigin();
 199        SliceType::IndexType       start = sliceRegion.GetIndex();
 200
 201        SliceType::Pointer seedsObject = SliceType::New();
 202        seedsObject->SetRegions(sliceRegion);
 203        seedsObject->SetOrigin(origin);
 204        seedsObject->SetSpacing(spacing);
 205        seedsObject->Allocate();
 206        seedsObject->FillBuffer(0);
 207
 208        SliceType::Pointer seedsBackground = SliceType::New();
 209        seedsBackground->SetRegions(sliceRegion);
 210        seedsBackground->SetOrigin(origin);
 211        seedsBackground->SetSpacing(spacing);
 212        seedsBackground->Allocate();
 213        seedsBackground->FillBuffer(0);
 214
 215        typedef itk::ImageRegionIteratorWithIndex<SliceType> SliceIteratorType;
 216        SliceIteratorType itSliceMori(tmp, tmp->GetLargestPossibleRegion());
 217
 218        int len = 0;
 219        int numO = 0;
 220        int numB = 0;
 221
 222        itSliceMori.GoToBegin();
 223
 224        while (!itSliceMori.IsAtEnd())
 225        {
 226        if (itSliceMori.Get() != 0)
 227        len++;
 228
 229        typename TImage::PixelType vessVal = vess->GetPixel(itSliceMori.GetIndex());
 230        typename TImage::PixelType imgVal = data->GetPixel(itSliceMori.GetIndex());
 231        typename TImage::PixelType bw = tmp->GetPixel(itSliceMori.GetIndex());
 232
 233        if (bw != 0 && vessVal > this->m_VessTh * this->m_vessMax / 100)
 234        {
 235        seedsObject->SetPixel(itSliceMori.GetIndex(), 1);
 236        numO++;
 237        }
 238
 239        itSliceMori++;
 240        }
 241
 242        if (len == 0)
 243        {
 244        std::cout << "Slice " << z << ": - brak pikseli zakwalifikowanych do drzewa" << std::endl;
 245        continue;
 246        }
 247
 248        SliceIteratorType itSliceVess(vess, vess->GetLargestPossibleRegion());
 249        itSliceVess.GoToBegin();
 250
 251        while (!itSliceVess.IsAtEnd())
 252        {
 253        if (itSliceVess.Get() < this->m_VessTh * this->m_vessMax / 100)
 254        {
 255        seedsBackground->SetPixel(itSliceVess.GetIndex(), 1);
 256        numB++;
 257        }
 258        itSliceVess++;
 259
 260        }
 261
 262        std::cout << "slice: " << z << "seeds obj: " << numO << " seeds bkg " << numB << std::endl;
 263        // przygotowania do Random Walkera
 264
 265        typename TImage::PixelType  *bufferImg = 0;
 266        typename TImage::PixelType  *bufferSeedsObj = 0;
 267        typename TImage::PixelType  *bufferSeedsBkg = 0;
 268
 269        int dims[2], dims2[2], dims3[2];
 270        double spac[2], orig[2], spac2[2], orig2[2], spac3[2], orig3[2];
 271
 272        image2buffer(bufferImg, dims, spac, orig, 2, data);
 273        image2buffer(bufferSeedsObj, dims2, spac2, orig2, 2, seedsObject);
 274        image2buffer(bufferSeedsBkg, dims3, spac3, orig3, 2, seedsBackground);
 275
 276        // TO DO HERE
 277        int width = dims[0];
 278        int height = dims[1];
 279
 280        PixelType *mask = new PixelType[width*height];
 281        double *probs = new double[2 * width*height];
 282
 283        int i, j, no_seeds = 0;
 284
 285        for (i = 0; i < width*height; i++)
 286
 287        if (bufferSeedsObj[i] == 1 || bufferSeedsBkg[i] == 1)
 288        no_seeds++;
 289
 290        int *seed_indexes = new int[no_seeds];
 291        int *seed_labels = new int[no_seeds];
 292
 293        for (j = 0, i = 0; i < width*height; i++)
 294        {
 295        if (bufferSeedsObj[i] == 1)
 296        {
 297        seed_indexes[j] = i + 1;//=i+1 (i)
 298        seed_labels[j] = this->m_replaceValue;
 299        j++;
 300        }
 301        else if (bufferSeedsBkg[i] == 1)
 302        {
 303        seed_indexes[j] = i + 1;//=i+1 (i)
 304        seed_labels[j] = 0;
 305        j++;
 306        }
 307        }
 308
 309        call_walker(mask, probs, bufferImg, height, width, seed_indexes, no_seeds, seed_labels, 2, this->m_beta);
 310
 311        delete[] bufferSeedsObj;
 312        delete[] bufferSeedsBkg;
 313        delete[] seed_indexes;
 314        delete[] seed_labels;
 315
 316        SliceType::Pointer img2 = buffer2image(mask, dims, spac, orig, 2);
 317
 318        SliceType::PixelType *probs_img = new PixelType[width*height];
 319        for (i = 0; i < width*height; i++)
 320        probs_img[i] = floor(255 * probs[i]);
 321
 322        SliceType::Pointer img3 = buffer2image(probs_img, dims, spac, orig, 2);
 323
 324        typedef itk::MaximumImageFilter<SliceType> MaxFilterType;
 325        MaxFilterType::Pointer maxFilter = MaxFilterType::New();
 326        maxFilter->SetInput(0, img2);
 327        maxFilter->SetInput(1, next);
 328        maxFilter->Update();
 329
 330        SliceType::Pointer sth = maxFilter->GetOutput();
 331        ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion());
 332
 333        delete[] probs;
 334     */
 335   }
 336 }
 337
 338 // -------------------------------------------------------------------------
 339 template< class _TImage, class _TLabels, class _TScalarImage >
 340 template< class _TBinaryTree, class _TData3D, class _TFusion >
 341 void
 342 _GoUp(
 343   _TBinaryTree& binaryTree,
 344   const _TData3D& data3D,
 345   const _TFusion& fusion,
 346   int numSlices
 347   )
 348 {
 349 }
 350
 351 #endif // __fpa__Common__SliceBySliceRandomWalker__hxx__
 352 // eof - $RCSfile$