itk/clitkConditionalBinaryErodeImageFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef clitkConditionalBinaryErodeImageFilter_txx
  19 #define clitkConditionalBinaryErodeImageFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkConditionalBinaryErodeImageFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "itkConstNeighborhoodIterator.h"
  31 #include "itkNeighborhoodIterator.h"
  32 #include "itkImageRegionIteratorWithIndex.h"
  33 #include "itkNeighborhoodInnerProduct.h"
  34 #include "itkImageRegionIterator.h"
  35 #include "itkImageRegionConstIterator.h"
  36 #include "itkNeighborhoodAlgorithm.h"
  37 #include "itkConstantBoundaryCondition.h"
  38 #include "itkOffset.h"
  39 #include "itkProgressReporter.h"
  40 #include "itkBinaryErodeImageFilter.h"
  41
  42 namespace clitk
  43 {
  44
  45   template <class TInputImage, class TOutputImage, class TKernel>
  46   ConditionalBinaryErodeImageFilter<TInputImage, TOutputImage, TKernel>
  47   ::ConditionalBinaryErodeImageFilter()
  48   {
  49     this->m_BoundaryToForeground = true;
  50   }
  51
  52
  53   template< class TInputImage, class TOutputImage, class TKernel>
  54   void
  55   ConditionalBinaryErodeImageFilter< TInputImage, TOutputImage, TKernel>
  56   ::GenerateData()
  57   {
  58     this->AllocateOutputs();
  59
  60     unsigned int i,j;
  61
  62     // Retrieve input and output pointers
  63     typename OutputImageType::Pointer output = this->GetOutput();
  64     typename InputImageType::ConstPointer input  = this->GetInput();
  65
  66     // Get values from superclass
  67     InputPixelType foregroundValue = this->GetForegroundValue();
  68     InputPixelType backgroundValue = this->GetBackgroundValue();
  69     KernelType kernel = this->GetKernel();
  70     InputSizeType radius;
  71     radius.Fill(1);
  72     typename TInputImage::RegionType inputRegion = input->GetBufferedRegion();
  73     typename TOutputImage::RegionType outputRegion = output->GetBufferedRegion();
  74
  75     // compute the size of the temp image. It is needed to create the progress
  76     // reporter.
  77     // The tmp image needs to be large enough to support:
  78     //   1. The size of the structuring element
  79     //   2. The size of the connectivity element (typically one)
  80     typename TInputImage::RegionType tmpRequestedRegion = outputRegion;
  81     typename TInputImage::RegionType paddedInputRegion
  82       = input->GetBufferedRegion();
  83     paddedInputRegion.PadByRadius( radius ); // to support boundary values
  84     InputSizeType padBy = radius;
  85     for (i=0; i < KernelDimension; ++i)
  86       {
  87         padBy[i] = (padBy[i]>kernel.GetRadius(i) ? padBy[i] : kernel.GetRadius(i));
  88       }
  89     tmpRequestedRegion.PadByRadius( padBy );
  90     tmpRequestedRegion.Crop( paddedInputRegion );
  91
  92     typename TInputImage::RegionType requiredInputRegion
  93       = input->GetBufferedRegion();
  94     requiredInputRegion.Crop( tmpRequestedRegion );
  95
  96     // Support progress methods/callbacks
  97     // Setup a progress reporter.  We have 4 stages to the algorithm so
  98     // pretend we have 4 times the number of pixels
  99     itk::ProgressReporter progress(this, 0,
 100                               outputRegion.GetNumberOfPixels() * 3
 101                               + tmpRequestedRegion.GetNumberOfPixels()
 102                               + requiredInputRegion.GetNumberOfPixels() );
 103
 104     // Allocate and reset output. We copy the input to the output,
 105     // except for pixels with DilateValue.  These pixels are initially
 106     // replaced with BackgroundValue and potentially replaced later with
 107     // DilateValue as the Minkowski sums are performed.
 108     itk::ImageRegionIterator<OutputImageType> outIt( output, outputRegion );
 109     //ImageRegionConstIterator<InputImageType> inIt( input, outputRegion );
 110
 111     for( outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt )
 112       {
 113         outIt.Set( foregroundValue );
 114         progress.CompletedPixel();
 115       }
 116
 117
 118     // Create the temp image for surface encoding
 119     // The temp image size is equal to the output requested region for thread
 120     // padded by max( connectivity neighborhood radius, SE kernel radius ).
 121     typedef itk::Image<unsigned char, TInputImage::ImageDimension> TempImageType;
 122     typename TempImageType::Pointer tmpImage = TempImageType::New();
 123
 124     // Define regions of temp image
 125     tmpImage->SetRegions( tmpRequestedRegion );
 126
 127     // Allocation.
 128     // Pay attention to the fact that here, the output is still not
 129     // allocated (so no extra memory needed for tmp image, if you
 130     // consider that you reserve som memory space for output)
 131     tmpImage->Allocate();
 132
 133     // First Stage
 134     // Copy the input image to the tmp image.
 135     // Tag the tmp Image.
 136     //     zero means background
 137     //     one means pixel on but not treated
 138     //     two means border pixel
 139     //     three means inner pixel
 140     static const unsigned char backgroundTag  = 0;
 141     static const unsigned char onTag          = 1;
 142     static const unsigned char borderTag      = 2;
 143     static const unsigned char innerTag       = 3;
 144
 145     if( !this->m_BoundaryToForeground )
 146       { tmpImage->FillBuffer( onTag ); }
 147     else
 148       { tmpImage->FillBuffer( backgroundTag ); }
 149
 150     // Iterators on input and tmp image
 151     // iterator on input
 152     itk::ImageRegionConstIterator<TInputImage> iRegIt( input, requiredInputRegion );
 153     // iterator on tmp image
 154     itk::ImageRegionIterator<TempImageType> tmpRegIt( tmpImage, requiredInputRegion );
 155
 156     for( iRegIt.GoToBegin(), tmpRegIt.GoToBegin();
 157          !tmpRegIt.IsAtEnd();
 158          ++iRegIt, ++tmpRegIt )
 159       {
 160         OutputPixelType pxl = iRegIt.Get();
 161
 162         //JV pxl != foregroundValue
 163         if( pxl == backgroundValue )
 164           { tmpRegIt.Set( onTag ); }
 165         else
 166           {
 167             // by default if it is not foreground, consider
 168             // it as background
 169             tmpRegIt.Set( backgroundTag );
 170           }
 171         progress.CompletedPixel();
 172       }
 173
 174     // Second stage
 175     // Border tracking and encoding
 176
 177     // Need to record index, use an iterator with index
 178     // Define iterators that will traverse the OUTPUT requested region
 179     // for thread and not the padded one. The tmp image has been padded
 180     // because in that way we will take care carefully at boundary
 181     // pixels of output requested region.  Take care means that we will
 182     // check if a boundary pixel is or not a border pixel.
 183     itk::ImageRegionIteratorWithIndex<TempImageType>
 184       tmpRegIndexIt( tmpImage, tmpRequestedRegion );
 185
 186     itk::ConstNeighborhoodIterator<TempImageType>
 187       oNeighbIt( radius, tmpImage, tmpRequestedRegion );
 188
 189     // Define boundaries conditions
 190     itk::ConstantBoundaryCondition<TempImageType> cbc;
 191     cbc.SetConstant( backgroundTag );
 192     oNeighbIt.OverrideBoundaryCondition(&cbc);
 193
 194     unsigned int neighborhoodSize = oNeighbIt.Size();
 195     unsigned int centerPixelCode = neighborhoodSize / 2;
 196
 197     std::queue<IndexType> propagQueue;
 198
 199     // Neighborhood iterators used to track the surface.
 200     //
 201     // Note the region specified for the first neighborhood iterator is
 202     // the requested region for the tmp image not the output image. This
 203     // is necessary because the NeighborhoodIterator relies on the
 204     // specified region to determine if you will ever query a boundary
 205     // condition pixel.  Since we call SetLocation on the neighbor of a
 206     // specified pixel, we have to set the region for the interator to
 207     // include any pixel we may set our location to.
 208     itk::NeighborhoodIterator<TempImageType>
 209       nit( radius, tmpImage, tmpRequestedRegion );
 210     nit.OverrideBoundaryCondition(&cbc);
 211     nit.GoToBegin();
 212
 213     itk::ConstNeighborhoodIterator<TempImageType>
 214       nnit( radius, tmpImage, tmpRequestedRegion );
 215     nnit.OverrideBoundaryCondition(&cbc);
 216     nnit.GoToBegin();
 217
 218     for( tmpRegIndexIt.GoToBegin(), oNeighbIt.GoToBegin();
 219          !tmpRegIndexIt.IsAtEnd();
 220          ++tmpRegIndexIt, ++oNeighbIt )
 221       {
 222
 223         unsigned char tmpValue = tmpRegIndexIt.Get();
 224
 225         // Test current pixel: it is active ( on ) or not?
 226         if( tmpValue == onTag )
 227           {
 228             // The current pixel has not been treated previously.  That
 229             // means that we do not know that it is an inner pixel of a
 230             // border pixel.
 231
 232             // Test current pixel: it is a border pixel or an inner pixel?
 233             bool bIsOnContour = false;
 234
 235             for (i = 0; i < neighborhoodSize; ++i)
 236               {
 237                 // If at least one neighbour pixel is off the center pixel
 238                 // belongs to contour
 239                 if( oNeighbIt.GetPixel( i ) == backgroundTag )
 240                   {
 241                     bIsOnContour = true;
 242                     break;
 243                   }
 244               }
 245
 246             if( bIsOnContour )
 247               {
 248                 // center pixel is a border pixel and due to the parsing, it is also
 249                 // a pixel which belongs to a new border connected component
 250                 // Now we will parse this border thanks to a burn procedure
 251
 252                 // mark pixel value as a border pixel
 253                 tmpRegIndexIt.Set( borderTag );
 254
 255                 // add it to border container.
 256                 // its code is center pixel code because it is the first pixel
 257                 // of the connected component border
 258
 259                 // paint the structuring element
 260                 typename NeighborIndexContainer::const_iterator itIdx;
 261                 NeighborIndexContainer& idxDifferenceSet
 262                   = this->GetDifferenceSet( centerPixelCode );
 263                 for( itIdx = idxDifferenceSet.begin();
 264                      itIdx != idxDifferenceSet.end();
 265                      ++itIdx )
 266                   {
 267                     IndexType idx = tmpRegIndexIt.GetIndex() + *itIdx;
 268                     if( outputRegion.IsInside( idx ) )
 269                       { output->SetPixel( idx, backgroundValue ); }
 270                   }
 271
 272                 // add it to queue
 273                 propagQueue.push( tmpRegIndexIt.GetIndex() );
 274
 275                 // now find all the border pixels
 276                 while ( !propagQueue.empty() )
 277                   {
 278                     // Extract pixel index from queue
 279                     IndexType currentIndex = propagQueue.front();
 280                     propagQueue.pop();
 281
 282                     nit += currentIndex - nit.GetIndex();
 283
 284                     for (i = 0; i < neighborhoodSize; ++i)
 285                       {
 286                         // If pixel has not been already treated and it is a pixel
 287                         // on, test if it is an inner pixel or a border pixel
 288
 289                         // Remark: all the pixels outside the image are set to
 290                         // backgroundTag thanks to boundary conditions. That means that if
 291                         // we enter in the next if-statement we are sure that the
 292                         // current neighbour pixel is in the image
 293                         if( nit.GetPixel( i ) == onTag )
 294                           {
 295                             // Check if it is an inner or border neighbour pixel
 296                             // Get index of current neighbour pixel
 297                             IndexType neighbIndex = nit.GetIndex( i );
 298
 299                             // Force location of neighbour iterator
 300                             nnit += neighbIndex - nnit.GetIndex();
 301
 302                             bool bIsOnBorder = false;
 303
 304                             for( j = 0; j < neighborhoodSize; ++j)
 305                               {
 306                                 // If at least one neighbour pixel is off the center
 307                                 // pixel belongs to border
 308                                 if( nnit.GetPixel(j) == backgroundTag )
 309                                   {
 310                                     bIsOnBorder = true;
 311                                     break;
 312                                   }
 313                               }
 314
 315
 316                             if( bIsOnBorder )
 317                               {
 318                                 // neighbour pixel is a border pixel
 319                                 // mark it
 320                                 bool status;
 321                                 nit.SetPixel( i, borderTag, status );
 322
 323                                 // check whether we could set the pixel.  can only set
 324                                 // the pixel if it is within the tmpimage
 325                                 if (status)
 326                                   {
 327                                     // add it to queue
 328                                     propagQueue.push( neighbIndex );
 329
 330                                     // paint the structuring element
 331                                     typename NeighborIndexContainer::const_iterator itIndex;
 332                                     NeighborIndexContainer& indexDifferenceSet
 333                                       = this->GetDifferenceSet( i );
 334                                     for( itIndex = indexDifferenceSet.begin();
 335                                          itIndex != indexDifferenceSet.end();
 336                                          ++itIndex )
 337                                       {
 338                                         IndexType idx = neighbIndex + *itIndex;
 339                                         if( outputRegion.IsInside( idx ) )
 340                                           { output->SetPixel( idx, backgroundValue ); }
 341                                       }
 342                                   }
 343                               }
 344                             else
 345                               {
 346                                 // neighbour pixel is an inner pixel
 347                                 bool status;
 348                                 nit.SetPixel( i, innerTag, status );
 349                               }
 350
 351                             progress.CompletedPixel();
 352                           } // if( nit.GetPixel( i ) == onTag )
 353
 354                       } // for (i = 0; i < neighborhoodSize; ++i)
 355
 356                   } // while ( !propagQueue.empty() )
 357
 358               } // if( bIsOnCountour )
 359             else
 360               { tmpRegIndexIt.Set( innerTag ); }
 361
 362             progress.CompletedPixel();
 363
 364           } // if( tmpRegIndexIt.Get() == onTag )
 365         else if( tmpValue == backgroundTag )
 366           { progress.CompletedPixel(); }
 367         // Here, the pixel is a background pixel ( value at 0 ) or an
 368         // already treated pixel:
 369         //     2 for border pixel, 3 for inner pixel
 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( 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                 //JV  != foregroundValue
 452                 if( !inputRegionForThread.IsInside( translatedIndex )
 453                     || input->GetPixel( translatedIndex ) == backgroundValue )
 454                   {
 455                     ouRegIndexIt.Set( backgroundValue );
 456                     break; // Do not need to examine other offsets because at least one
 457                     // input pixel has been translated on current output pixel.
 458                   }
 459               }
 460
 461             ++ouRegIndexIt;
 462             progress.CompletedPixel();
 463           }
 464       }
 465     else
 466       {
 467         while( !ouRegIndexIt.IsAtEnd() )
 468           {
 469             IndexType currentIndex  = ouRegIndexIt.GetIndex();
 470             for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
 471               {
 472                 IndexType translatedIndex = currentIndex - *vecIt;
 473
 474                 if( inputRegionForThread.IsInside( translatedIndex )
 475                     //JV != foregroundValue
 476                     && input->GetPixel( translatedIndex ) == backgroundValue )
 477                   {
 478                     ouRegIndexIt.Set( backgroundValue );
 479                     break;
 480                   }
 481               }
 482
 483             ++ouRegIndexIt;
 484             progress.CompletedPixel();
 485           }
 486       }
 487
 488     // now, we must to restore the background values
 489     itk::ImageRegionConstIterator<InputImageType> inIt( input, outputRegion );
 490
 491     for( inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt )
 492       {
 493         InputPixelType inValue = inIt.Get();
 494         OutputPixelType outValue = outIt.Get();
 495         // JV != foregroundValue
 496         if ( outValue == backgroundValue && inValue == backgroundValue )
 497           { outIt.Set( static_cast<OutputPixelType>( inValue ) ); }
 498         progress.CompletedPixel();
 499       }
 500
 501
 502   }
 503
 504
 505   /**
 506    * Standard "PrintSelf" method
 507    */
 508   template <class TInputImage, class TOutput, class TKernel>
 509   void
 510   ConditionalBinaryErodeImageFilter<TInputImage, TOutput, TKernel>
 511   ::PrintSelf( std::ostream& os, itk::Indent indent) const
 512   {
 513     Superclass::PrintSelf( os, indent );
 514     os << indent << "Dilate Value: " << static_cast<typename itk::NumericTraits<InputPixelType>::PrintType>( this->GetForegroundValue() ) << std::endl;
 515   }
 516
 517 }//end clitk
 518
 519 #endif //#define clitkConditionalBinaryErodeImageFilter_txx