1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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
21 /* =================================================
22 * @file clitkConditionalBinaryDilateImageFilter.txx
28 ===================================================*/
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"
45 template <class TInputImage, class TOutputImage, class TKernel>
46 ConditionalBinaryDilateImageFilter<TInputImage, TOutputImage, TKernel>
47 ::ConditionalBinaryDilateImageFilter()
49 this->m_BoundaryToForeground = false;
53 template< class TInputImage, class TOutputImage, class TKernel>
55 ConditionalBinaryDilateImageFilter< TInputImage, TOutputImage, TKernel>
58 this->AllocateOutputs();
62 // Retrieve input and output pointers
63 typename OutputImageType::Pointer output = this->GetOutput();
64 typename InputImageType::ConstPointer input = this->GetInput();
66 // Get values from superclass
67 InputPixelType foregroundValue = this->GetForegroundValue();
68 InputPixelType backgroundValue = this->GetBackgroundValue();
69 KernelType kernel = this->GetKernel();
72 typename TInputImage::RegionType inputRegion = input->GetBufferedRegion();
73 typename TOutputImage::RegionType outputRegion = output->GetBufferedRegion();
75 // compute the size of the temp image. It is needed to create the progress
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)
87 padBy[i] = (padBy[i]>kernel.GetRadius(i) ? padBy[i] : kernel.GetRadius(i));
89 tmpRequestedRegion.PadByRadius( padBy );
90 tmpRequestedRegion.Crop( paddedInputRegion );
92 typename TInputImage::RegionType requiredInputRegion
93 = input->GetBufferedRegion();
94 requiredInputRegion.Crop( tmpRequestedRegion );
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() );
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 );
111 for( inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt )
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
120 { outIt.Set( static_cast<OutputPixelType>( value ) ); }
121 progress.CompletedPixel();
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();
131 // Define regions of temp image
132 tmpImage->SetRegions( tmpRequestedRegion );
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();
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;
154 if( this->m_BoundaryToForeground )
155 { tmpImage->FillBuffer( onTag ); }
157 { tmpImage->FillBuffer( backgroundTag ); }
159 // Iterators on input and tmp image
161 itk::ImageRegionConstIterator<TInputImage> iRegIt( input, requiredInputRegion );
162 // iterator on tmp image
163 itk::ImageRegionIterator<TempImageType> tmpRegIt( tmpImage, requiredInputRegion );
165 for( iRegIt.GoToBegin(), tmpRegIt.GoToBegin();
167 ++iRegIt, ++tmpRegIt )
169 OutputPixelType pxl = iRegIt.Get();
170 // JV == foregroundValue
171 if( pxl == foregroundValue )
172 { tmpRegIt.Set( onTag ); }
174 else if ( pxl == backgroundValue )
175 { tmpRegIt.Set( bgTag ); }
178 // by default if it is not foreground, consider
180 tmpRegIt.Set( backgroundTag );
182 progress.CompletedPixel();
187 // Border tracking and encoding
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 );
198 itk::ConstNeighborhoodIterator<TempImageType>
199 oNeighbIt( radius, tmpImage, tmpRequestedRegion );
201 // Define boundaries conditions
202 itk::ConstantBoundaryCondition<TempImageType> cbc;
203 cbc.SetConstant( backgroundTag );
204 oNeighbIt.OverrideBoundaryCondition(&cbc);
206 unsigned int neighborhoodSize = oNeighbIt.Size();
207 unsigned int centerPixelCode = neighborhoodSize / 2;
209 std::queue<IndexType> propagQueue;
211 // Neighborhood iterators used to track the surface.
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);
225 itk::ConstNeighborhoodIterator<TempImageType>
226 nnit( radius, tmpImage, tmpRequestedRegion );
227 nnit.OverrideBoundaryCondition(&cbc);
230 for( tmpRegIndexIt.GoToBegin(), oNeighbIt.GoToBegin();
231 !tmpRegIndexIt.IsAtEnd();
232 ++tmpRegIndexIt, ++oNeighbIt )
234 // Test current pixel: it is active ( on ) or not?
235 if( tmpRegIndexIt.Get() == onTag )
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
241 // Test current pixel: it is a border pixel or an inner pixel?
242 bool bIsOnContour = false;
244 for (i = 0; i < neighborhoodSize; ++i)
246 // If at least one neighbour pixel is off the center pixel
247 // belongs to contour
249 if( oNeighbIt.GetPixel( i ) == bgTag )
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
262 // mark pixel value as a border pixel
263 tmpRegIndexIt.Set( borderTag );
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
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();
277 IndexType idx = tmpRegIndexIt.GetIndex() + *itIdx;
278 if( outputRegion.IsInside( idx ) )
280 // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
281 if (input->GetPixel(idx)==backgroundValue)
282 output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
287 propagQueue.push( tmpRegIndexIt.GetIndex() );
289 // now find all the border pixels
290 while ( !propagQueue.empty() )
292 // Extract pixel index from queue
293 IndexType currentIndex = propagQueue.front();
296 nit += currentIndex - nit.GetIndex();
298 for (i = 0; i < neighborhoodSize; ++i)
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
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 )
309 // Check if it is an inner or border neighbour pixel
310 // Get index of current neighbour pixel
311 IndexType neighbIndex = nit.GetIndex( i );
313 // Force location of neighbour iterator
314 nnit += neighbIndex - nnit.GetIndex();
316 bool bIsOnBorder = false;
318 for( j = 0; j < neighborhoodSize; ++j)
320 // If at least one neighbour pixel is off the center
321 // pixel belongs to border
322 if( nnit.GetPixel(j) == bgTag )
332 // neighbour pixel is a border pixel
335 nit.SetPixel( i, borderTag, status );
337 // check whether we could set the pixel. can only set
338 // the pixel if it is within the tmpimage
342 propagQueue.push( neighbIndex );
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();
352 IndexType idx = neighbIndex + *itIndex;
353 if( outputRegion.IsInside( idx ) )
355 // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
356 if (input->GetPixel(idx)==backgroundValue)
357 output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
364 // neighbour pixel is an inner pixel
366 nit.SetPixel( i, innerTag, status );
369 progress.CompletedPixel();
370 } // if( nit.GetPixel( i ) == onTag )
372 } // for (i = 0; i < neighborhoodSize; ++i)
374 } // while ( !propagQueue.empty() )
376 } // if( bIsOnCountour )
378 { tmpRegIndexIt.Set( innerTag ); }
380 } // if( tmpRegIndexIt.Get() == onTag )
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
390 // Deallocate tmpImage
391 tmpImage->Initialize();
394 // traverse structure of border and SE CCs, and paint output image
396 // Let's consider the the set of the ON elements of the input image as X.
398 // Let's consider the structuring element as B = {B0, B1, ..., Bn},
399 // where Bi denotes a connected component of B.
401 // Let's consider bi, i in [0,n], an arbitrary point of Bi.
404 // We use hence the next property in order to compute minkoswki
405 // addition ( which will be written (+) ):
407 // X (+) B = ( Xb0 UNION Xb1 UNION ... Xbn ) UNION ( BORDER(X) (+) B ),
409 // where Xbi is the set X translated with respect to vector bi :
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 )
414 // Define boundaries conditions
415 itk::ConstantBoundaryCondition<TOutputImage> obc;
416 obc.SetConstant( static_cast<OutputPixelType> ( backgroundValue ) );
418 itk::NeighborhoodIterator<OutputImageType>
419 onit( kernel.GetRadius(), output, outputRegion );
420 onit.OverrideBoundaryCondition(&obc);
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();
432 // iterator on output image
433 itk::ImageRegionIteratorWithIndex<OutputImageType>
434 ouRegIndexIt(output, outputRegion );
435 ouRegIndexIt.GoToBegin();
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;
447 // Pad the input region by the kernel
448 inputRegionForThread.PadByRadius( kernel.GetRadius() );
449 inputRegionForThread.Crop( input->GetBufferedRegion() );
451 if( this->m_BoundaryToForeground )
453 while( !ouRegIndexIt.IsAtEnd() )
455 // Retrieve index of current output pixel
456 IndexType currentIndex = ouRegIndexIt.GetIndex();
457 for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
460 IndexType translatedIndex = currentIndex - *vecIt;
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 )
470 || ( (input->GetPixel( translatedIndex ) == foregroundValue )
472 && (input->GetPixel( currentIndex ) == backgroundValue ) ) )
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.
481 progress.CompletedPixel();
486 while( !ouRegIndexIt.IsAtEnd() )
488 IndexType currentIndex = ouRegIndexIt.GetIndex();
489 for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
491 IndexType translatedIndex = currentIndex - *vecIt;
493 if( inputRegionForThread.IsInside( translatedIndex )
494 // JV == foregroundValue
495 && ( (input->GetPixel( translatedIndex ) == foregroundValue )
497 && (input->GetPixel( currentIndex ) == backgroundValue ) ) )
500 ouRegIndexIt.Set( static_cast<OutputPixelType> ( foregroundValue ) );
506 progress.CompletedPixel();
513 * Standard "PrintSelf" method
515 template <class TInputImage, class TOutput, class TKernel>
517 ConditionalBinaryDilateImageFilter<TInputImage, TOutput, TKernel>
518 ::PrintSelf( std::ostream& os, itk::Indent indent) const
520 Superclass::PrintSelf( os, indent );
521 os << indent << "Dilate Value: " << static_cast<typename itk::NumericTraits<InputPixelType>::PrintType>( this->GetForegroundValue() ) << std::endl;
524 } // end namespace itk