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