itk/clitkConditionalBinaryDilateImageFilter.txx

   1 #ifndef clitkConditionalBinaryDilateImageFilter_txx
   2 #define clitkConditionalBinaryDilateImageFilter_txx
   3
   4 /* =================================================
   5  * @file   clitkConditionalBinaryDilateImageFilter.txx
   6  * @author
   7  * @date
   8  *
   9  * @brief
  10  *
  11  ===================================================*/
  12
  13 #include "itkConstNeighborhoodIterator.h"
  14 #include "itkNeighborhoodIterator.h"
  15 #include "itkImageRegionIteratorWithIndex.h"
  16 #include "itkNeighborhoodInnerProduct.h"
  17 #include "itkImageRegionIterator.h"
  18 #include "itkImageRegionConstIterator.h"
  19 #include "itkNeighborhoodAlgorithm.h"
  20 #include "itkConstantBoundaryCondition.h"
  21 #include "itkOffset.h"
  22 #include "itkProgressReporter.h"
  23 #include "itkBinaryDilateImageFilter.h"
  24
  25 namespace clitk
  26 {
  27
  28   template <class TInputImage, class TOutputImage, class TKernel>
  29   ConditionalBinaryDilateImageFilter<TInputImage, TOutputImage, TKernel>
  30   ::ConditionalBinaryDilateImageFilter()
  31   {
  32     this->m_BoundaryToForeground = false;
  33   }
  34
  35
  36   template< class TInputImage, class TOutputImage, class TKernel>
  37   void
  38   ConditionalBinaryDilateImageFilter< TInputImage, TOutputImage, TKernel>
  39   ::GenerateData()
  40   {
  41     this->AllocateOutputs();
  42
  43     unsigned int i,j;
  44
  45     // Retrieve input and output pointers
  46     typename OutputImageType::Pointer output = this->GetOutput();
  47     typename InputImageType::ConstPointer input  = this->GetInput();
  48
  49     // Get values from superclass
  50     InputPixelType foregroundValue = this->GetForegroundValue();
  51     InputPixelType backgroundValue = this->GetBackgroundValue();
  52     KernelType kernel = this->GetKernel();
  53     InputSizeType radius;
  54     radius.Fill(1);
  55     typename TInputImage::RegionType inputRegion = input->GetBufferedRegion();
  56     typename TOutputImage::RegionType outputRegion = output->GetBufferedRegion();
  57
  58     // compute the size of the temp image. It is needed to create the progress
  59     // reporter.
  60     // The tmp image needs to be large enough to support:
  61     //   1. The size of the structuring element
  62     //   2. The size of the connectivity element (typically one)
  63     typename TInputImage::RegionType tmpRequestedRegion = outputRegion;
  64     typename TInputImage::RegionType paddedInputRegion
  65       = input->GetBufferedRegion();
  66     paddedInputRegion.PadByRadius( radius ); // to support boundary values
  67     InputSizeType padBy = radius;
  68     for (i=0; i < KernelDimension; ++i)
  69       {
  70         padBy[i] = (padBy[i]>kernel.GetRadius(i) ? padBy[i] : kernel.GetRadius(i));
  71       }
  72     tmpRequestedRegion.PadByRadius( padBy );
  73     tmpRequestedRegion.Crop( paddedInputRegion );
  74
  75     typename TInputImage::RegionType requiredInputRegion
  76       = input->GetBufferedRegion();
  77     requiredInputRegion.Crop( tmpRequestedRegion );
  78
  79     // Support progress methods/callbacks
  80     // Setup a progress reporter.  We have 4 stages to the algorithm so
  81     // pretend we have 4 times the number of pixels
  82     itk::ProgressReporter progress(this, 0,
  83                                    outputRegion.GetNumberOfPixels() * 2
  84                                    + tmpRequestedRegion.GetNumberOfPixels()
  85                                    + requiredInputRegion.GetNumberOfPixels() );
  86
  87     // Allocate and reset output. We copy the input to the output,
  88     // except for pixels with DilateValue.  These pixels are initially
  89     // replaced with BackgroundValue and potentially replaced later with
  90     // DilateValue as the Minkowski sums are performed.
  91     itk::ImageRegionIterator<OutputImageType> outIt( output, outputRegion );
  92     itk::ImageRegionConstIterator<InputImageType> inIt( input, outputRegion );
  93
  94     for( inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt )
  95       {
  96         InputPixelType value = inIt.Get();
  97         // replace foreground pixels with the default background
  98         //JV == foregroundValue
  99         if (value == backgroundValue)
 100           { outIt.Set( static_cast<OutputPixelType> ( backgroundValue ) ); }
 101         // keep all of the original background values intact
 102         else
 103           { outIt.Set( static_cast<OutputPixelType>( value ) ); }
 104         progress.CompletedPixel();
 105       }
 106
 107
 108     // Create the temp image for surface encoding
 109     // The temp image size is equal to the output requested region for thread
 110     // padded by max( connectivity neighborhood radius, SE kernel radius ).
 111     typedef itk::Image<unsigned char, TInputImage::ImageDimension> TempImageType;
 112     typename TempImageType::Pointer tmpImage = TempImageType::New();
 113
 114     // Define regions of temp image
 115     tmpImage->SetRegions( tmpRequestedRegion );
 116
 117     // Allocation.
 118     // Pay attention to the fact that here, the output is still not
 119     // allocated (so no extra memory needed for tmp image, if you
 120     // consider that you reserve som memory space for output)
 121     tmpImage->Allocate();
 122
 123     // First Stage
 124     // Copy the input image to the tmp image.
 125     // Tag the tmp Image.
 126     //     zero means background
 127     //     one means pixel on but not treated
 128     //     two means border pixel
 129     //     three means inner pixel
 130     static const unsigned char backgroundTag  = 0;
 131     static const unsigned char onTag          = 1;
 132     static const unsigned char borderTag      = 2;
 133     static const unsigned char innerTag       = 3;
 134     // JV add the tag in which can be dilated
 135     static const unsigned char bgTag       = 4;
 136
 137     if( this->m_BoundaryToForeground )
 138       { tmpImage->FillBuffer( onTag ); }
 139     else
 140       { tmpImage->FillBuffer( backgroundTag ); }
 141
 142     // Iterators on input and tmp image
 143     // iterator on input
 144     itk::ImageRegionConstIterator<TInputImage> iRegIt( input, requiredInputRegion );
 145     // iterator on tmp image
 146     itk::ImageRegionIterator<TempImageType> tmpRegIt( tmpImage, requiredInputRegion );
 147
 148     for( iRegIt.GoToBegin(), tmpRegIt.GoToBegin();
 149          !tmpRegIt.IsAtEnd();
 150          ++iRegIt, ++tmpRegIt )
 151       {
 152         OutputPixelType pxl = iRegIt.Get();
 153         // JV == foregroundValue
 154         if( pxl == foregroundValue )
 155           { tmpRegIt.Set( onTag ); }
 156         // JV add the bgTag
 157         else if ( pxl == backgroundValue )
 158           { tmpRegIt.Set( bgTag ); }
 159         else
 160           {
 161             // by default if it is not foreground, consider
 162             // it as background
 163             tmpRegIt.Set( backgroundTag );
 164           }
 165         progress.CompletedPixel();
 166       }
 167
 168
 169     // Second stage
 170     // Border tracking and encoding
 171
 172     // Need to record index, use an iterator with index
 173     // Define iterators that will traverse the OUTPUT requested region
 174     // for thread and not the padded one. The tmp image has been padded
 175     // because in that way we will take care carefully at boundary
 176     // pixels of output requested region.  Take care means that we will
 177     // check if a boundary pixel is or not a border pixel.
 178     itk::ImageRegionIteratorWithIndex<TempImageType>
 179       tmpRegIndexIt( tmpImage, tmpRequestedRegion );
 180
 181     itk::ConstNeighborhoodIterator<TempImageType>
 182       oNeighbIt( radius, tmpImage, tmpRequestedRegion );
 183
 184     // Define boundaries conditions
 185     itk::ConstantBoundaryCondition<TempImageType> cbc;
 186     cbc.SetConstant( backgroundTag );
 187     oNeighbIt.OverrideBoundaryCondition(&cbc);
 188
 189     unsigned int neighborhoodSize = oNeighbIt.Size();
 190     unsigned int centerPixelCode = neighborhoodSize / 2;
 191
 192     std::queue<IndexType> propagQueue;
 193
 194     // Neighborhood iterators used to track the surface.
 195     //
 196     // Note the region specified for the first neighborhood iterator is
 197     // the requested region for the tmp image not the output image. This
 198     // is necessary because the NeighborhoodIterator relies on the
 199     // specified region to determine if you will ever query a boundary
 200     // condition pixel.  Since we call SetLocation on the neighbor of a
 201     // specified pixel, we have to set the region for the interator to
 202     // include any pixel we may set our location to.
 203     itk::NeighborhoodIterator<TempImageType>
 204       nit( radius, tmpImage, tmpRequestedRegion );
 205     nit.OverrideBoundaryCondition(&cbc);
 206     nit.GoToBegin();
 207
 208     itk::ConstNeighborhoodIterator<TempImageType>
 209       nnit( radius, tmpImage, tmpRequestedRegion );
 210     nnit.OverrideBoundaryCondition(&cbc);
 211     nnit.GoToBegin();
 212
 213     for( tmpRegIndexIt.GoToBegin(), oNeighbIt.GoToBegin();
 214          !tmpRegIndexIt.IsAtEnd();
 215          ++tmpRegIndexIt, ++oNeighbIt )
 216       {
 217         // Test current pixel: it is active ( on ) or not?
 218         if( tmpRegIndexIt.Get() == onTag )
 219           {
 220             // The current pixel has not been treated previously.  That
 221             // means that we do not know that it is an inner pixel of a
 222             // border pixel.
 223
 224             // Test current pixel: it is a border pixel or an inner pixel?
 225             bool bIsOnContour = false;
 226
 227             for (i = 0; i < neighborhoodSize; ++i)
 228               {
 229                 // If at least one neighbour pixel is off the center pixel
 230                 // belongs to contour
 231                 //JV verify bg
 232                 if( oNeighbIt.GetPixel( i ) == bgTag )
 233                   {
 234                     bIsOnContour = true;
 235                     break;
 236                   }
 237               }
 238
 239             if( bIsOnContour )
 240               {
 241                 // center pixel is a border pixel and due to the parsing, it is also
 242                 // a pixel which belongs to a new border connected component
 243                 // Now we will parse this border thanks to a burn procedure
 244
 245                 // mark pixel value as a border pixel
 246                 tmpRegIndexIt.Set( borderTag );
 247
 248                 // add it to border container.
 249                 // its code is center pixel code because it is the first pixel
 250                 // of the connected component border
 251
 252                 // paint the structuring element
 253                 typename NeighborIndexContainer::const_iterator itIdx;
 254                 NeighborIndexContainer& idxDifferenceSet
 255                   = this->GetDifferenceSet( centerPixelCode );
 256                 for( itIdx = idxDifferenceSet.begin();
 257                      itIdx != idxDifferenceSet.end();
 258                      ++itIdx )
 259                   {
 260                     IndexType idx = tmpRegIndexIt.GetIndex() + *itIdx;
 261                     if( outputRegion.IsInside( idx ) )
 262                       {
 263                         // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
 264                         if (input->GetPixel(idx)==backgroundValue)
 265                           output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
 266                       }
 267                   }
 268
 269                 // add it to queue
 270                 propagQueue.push( tmpRegIndexIt.GetIndex() );
 271
 272                 // now find all the border pixels
 273                 while ( !propagQueue.empty() )
 274                   {
 275                     // Extract pixel index from queue
 276                     IndexType currentIndex = propagQueue.front();
 277                     propagQueue.pop();
 278
 279                     nit += currentIndex - nit.GetIndex();
 280
 281                     for (i = 0; i < neighborhoodSize; ++i)
 282                       {
 283                         // If pixel has not been already treated and it is a pixel
 284                         // on, test if it is an inner pixel or a border pixel
 285
 286                         // Remark: all the pixels outside the image are set to
 287                         // backgroundTag thanks to boundary conditions. That means that if
 288                         // we enter in the next if-statement we are sure that the
 289                         // current neighbour pixel is in the image
 290                         if( nit.GetPixel( i ) == onTag )
 291                           {
 292                             // Check if it is an inner or border neighbour pixel
 293                             // Get index of current neighbour pixel
 294                             IndexType neighbIndex = nit.GetIndex( i );
 295
 296                             // Force location of neighbour iterator
 297                             nnit += neighbIndex - nnit.GetIndex();
 298
 299                             bool bIsOnBorder = false;
 300
 301                             for( j = 0; j < neighborhoodSize; ++j)
 302                               {
 303                                 // If at least one neighbour pixel is off the center
 304                                 // pixel belongs to border
 305                                 if( nnit.GetPixel(j) == bgTag )
 306                                   {
 307                                     bIsOnBorder = true;
 308                                     break;
 309                                   }
 310                               }
 311
 312
 313                             if( bIsOnBorder )
 314                               {
 315                                 // neighbour pixel is a border pixel
 316                                 // mark it
 317                                 bool status;
 318                                 nit.SetPixel( i, borderTag, status );
 319
 320                                 // check whether we could set the pixel.  can only set
 321                                 // the pixel if it is within the tmpimage
 322                                 if (status)
 323                                   {
 324                                     // add it to queue
 325                                     propagQueue.push( neighbIndex );
 326
 327                                     // paint the structuring element
 328                                     typename NeighborIndexContainer::const_iterator itIndex;
 329                                     NeighborIndexContainer& indexDifferenceSet
 330                                       = this->GetDifferenceSet( i );
 331                                     for( itIndex = indexDifferenceSet.begin();
 332                                          itIndex != indexDifferenceSet.end();
 333                                          ++itIndex )
 334                                       {
 335                                         IndexType idx = neighbIndex + *itIndex;
 336                                         if( outputRegion.IsInside( idx ) )
 337                                           {
 338                                             // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
 339                                             if (input->GetPixel(idx)==backgroundValue)
 340                                               output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
 341                                           }
 342                                       }
 343                                   }
 344                               }
 345                             else
 346                               {
 347                                 // neighbour pixel is an inner pixel
 348                                 bool status;
 349                                 nit.SetPixel( i, innerTag, status );
 350                               }
 351
 352                             progress.CompletedPixel();
 353                           } // if( nit.GetPixel( i ) == onTag )
 354
 355                       } // for (i = 0; i < neighborhoodSize; ++i)
 356
 357                   } // while ( !propagQueue.empty() )
 358
 359               } // if( bIsOnCountour )
 360             else
 361               { tmpRegIndexIt.Set( innerTag ); }
 362
 363           } // if( tmpRegIndexIt.Get() == onTag )
 364         else
 365           { progress.CompletedPixel(); }
 366         // Here, the pixel is a background pixel ( value at 0 ) or an
 367         // already treated pixel:
 368         //     2 for border pixel, 3 for inner pixel
 369
 370       }
 371
 372
 373     // Deallocate tmpImage
 374     tmpImage->Initialize();
 375
 376     // Third Stage
 377     // traverse structure of border and SE CCs, and paint output image
 378
 379     // Let's consider the the set of the ON elements of the input image as X.
 380     //
 381     // Let's consider the structuring element as B = {B0, B1, ..., Bn},
 382     // where Bi denotes a connected component of B.
 383     //
 384     // Let's consider bi, i in [0,n], an arbitrary point of Bi.
 385     //
 386
 387     // We use hence the next property in order to compute minkoswki
 388     // addition ( which will be written (+) ):
 389     //
 390     // X (+) B = ( Xb0 UNION Xb1 UNION ... Xbn ) UNION ( BORDER(X) (+) B ),
 391     //
 392     // where Xbi is the set X translated with respect to vector bi :
 393     //
 394     // Xbi = { x + bi, x belongs to X } where BORDER(X) is the extracted
 395     // border of X ( 8 connectivity in 2D, 26 in 3D )
 396
 397     // Define boundaries conditions
 398     itk::ConstantBoundaryCondition<TOutputImage> obc;
 399     obc.SetConstant( static_cast<OutputPixelType> ( backgroundValue ) );
 400
 401     itk::NeighborhoodIterator<OutputImageType>
 402       onit( kernel.GetRadius(), output, outputRegion );
 403     onit.OverrideBoundaryCondition(&obc);
 404     onit.GoToBegin();
 405
 406
 407     // Paint input image translated with respect to the SE CCs vectors
 408     // --> "( Xb0 UNION Xb1 UNION ... Xbn )"
 409     typename Superclass::ComponentVectorConstIterator vecIt;
 410     typename Superclass::ComponentVectorConstIterator vecBeginIt
 411       = this->KernelCCVectorBegin();
 412     typename Superclass::ComponentVectorConstIterator vecEndIt
 413       = this->KernelCCVectorEnd();
 414
 415     // iterator on output image
 416     itk::ImageRegionIteratorWithIndex<OutputImageType>
 417       ouRegIndexIt(output, outputRegion );
 418     ouRegIndexIt.GoToBegin();
 419
 420     // InputRegionForThread is the output region for thread padded by
 421     // kerne lradius We must traverse this padded region because some
 422     // border pixel in the added band ( the padded band is the region
 423     // added after padding ) may be responsible to the painting of some
 424     // pixel in the non padded region.  This happens typically when a
 425     // non centered SE is used, a kind of shift is done on the "on"
 426     // pixels of image. Consequently some pixels in the added band can
 427     // appear in the current region for thread due to shift effect.
 428     typename InputImageType::RegionType inputRegionForThread = outputRegion;
 429
 430     // Pad the input region by the kernel
 431     inputRegionForThread.PadByRadius( kernel.GetRadius() );
 432     inputRegionForThread.Crop( input->GetBufferedRegion() );
 433
 434     if( this->m_BoundaryToForeground )
 435       {
 436         while( !ouRegIndexIt.IsAtEnd() )
 437           {
 438             // Retrieve index of current output pixel
 439             IndexType currentIndex  = ouRegIndexIt.GetIndex();
 440             for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
 441               {
 442                 // Translate
 443                 IndexType translatedIndex = currentIndex - *vecIt;
 444
 445                 // translated index now is an index in input image in the
 446                 // output requested region padded. Theoretically, this translated
 447                 // index must be inside the padded region.
 448                 // If the pixel in the input image at the translated index
 449                 // has a value equal to the dilate one, this means
 450                 // that the output pixel at currentIndex will be on in the output.
 451                 if( !inputRegionForThread.IsInside( translatedIndex )
 452                     // JV == foreground
 453                     ||   (  (input->GetPixel( translatedIndex ) == foregroundValue )
 454                             // JV add bg
 455                            &&  (input->GetPixel( currentIndex ) == backgroundValue ) ) )
 456                   {
 457                     ouRegIndexIt.Set( static_cast<OutputPixelType> ( foregroundValue ) );
 458                     break; // Do not need to examine other offsets because at least one
 459                     // input pixel has been translated on current output pixel.
 460                   }
 461               }
 462
 463             ++ouRegIndexIt;
 464             progress.CompletedPixel();
 465           }
 466       }
 467     else
 468       {
 469         while( !ouRegIndexIt.IsAtEnd() )
 470           {
 471             IndexType currentIndex  = ouRegIndexIt.GetIndex();
 472             for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
 473               {
 474                 IndexType translatedIndex = currentIndex - *vecIt;
 475
 476                 if( inputRegionForThread.IsInside( translatedIndex )
 477                     // JV == foregroundValue
 478                     &&  (  (input->GetPixel( translatedIndex ) == foregroundValue )
 479                            // JV add bg
 480                            &&  (input->GetPixel( currentIndex ) == backgroundValue ) ) )
 481
 482                   {
 483                     ouRegIndexIt.Set( static_cast<OutputPixelType> ( foregroundValue ) );
 484                     break;
 485                   }
 486               }
 487
 488             ++ouRegIndexIt;
 489             progress.CompletedPixel();
 490           }
 491       }
 492   }
 493
 494
 495   /**
 496    * Standard "PrintSelf" method
 497    */
 498   template <class TInputImage, class TOutput, class TKernel>
 499   void
 500   ConditionalBinaryDilateImageFilter<TInputImage, TOutput, TKernel>
 501   ::PrintSelf( std::ostream& os, itk::Indent indent) const
 502   {
 503     Superclass::PrintSelf( os, indent );
 504     os << indent << "Dilate Value: " << static_cast<typename itk::NumericTraits<InputPixelType>::PrintType>( this->GetForegroundValue() ) << std::endl;
 505   }
 506
 507 } // end namespace itk
 508
 509 #endif