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 clitkConditionalBinaryErodeImageFilter_txx
19 #define clitkConditionalBinaryErodeImageFilter_txx
21 /* =================================================
22 * @file clitkConditionalBinaryErodeImageFilter.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 "itkBinaryErodeImageFilter.h"
45 template <class TInputImage, class TOutputImage, class TKernel>
46 ConditionalBinaryErodeImageFilter<TInputImage, TOutputImage, TKernel>
47 ::ConditionalBinaryErodeImageFilter()
49 this->m_BoundaryToForeground = true;
53 template< class TInputImage, class TOutputImage, class TKernel>
55 ConditionalBinaryErodeImageFilter< 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() * 3
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 //ImageRegionConstIterator<InputImageType> inIt( input, outputRegion );
111 for( outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt )
113 outIt.Set( foregroundValue );
114 progress.CompletedPixel();
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();
124 // Define regions of temp image
125 tmpImage->SetRegions( tmpRequestedRegion );
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();
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;
145 if( !this->m_BoundaryToForeground )
146 { tmpImage->FillBuffer( onTag ); }
148 { tmpImage->FillBuffer( backgroundTag ); }
150 // Iterators on input and tmp image
152 itk::ImageRegionConstIterator<TInputImage> iRegIt( input, requiredInputRegion );
153 // iterator on tmp image
154 itk::ImageRegionIterator<TempImageType> tmpRegIt( tmpImage, requiredInputRegion );
156 for( iRegIt.GoToBegin(), tmpRegIt.GoToBegin();
158 ++iRegIt, ++tmpRegIt )
160 OutputPixelType pxl = iRegIt.Get();
162 //JV pxl != foregroundValue
163 if( pxl == backgroundValue )
164 { tmpRegIt.Set( onTag ); }
167 // by default if it is not foreground, consider
169 tmpRegIt.Set( backgroundTag );
171 progress.CompletedPixel();
175 // Border tracking and encoding
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 );
186 itk::ConstNeighborhoodIterator<TempImageType>
187 oNeighbIt( radius, tmpImage, tmpRequestedRegion );
189 // Define boundaries conditions
190 itk::ConstantBoundaryCondition<TempImageType> cbc;
191 cbc.SetConstant( backgroundTag );
192 oNeighbIt.OverrideBoundaryCondition(&cbc);
194 unsigned int neighborhoodSize = oNeighbIt.Size();
195 unsigned int centerPixelCode = neighborhoodSize / 2;
197 std::queue<IndexType> propagQueue;
199 // Neighborhood iterators used to track the surface.
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);
213 itk::ConstNeighborhoodIterator<TempImageType>
214 nnit( radius, tmpImage, tmpRequestedRegion );
215 nnit.OverrideBoundaryCondition(&cbc);
218 for( tmpRegIndexIt.GoToBegin(), oNeighbIt.GoToBegin();
219 !tmpRegIndexIt.IsAtEnd();
220 ++tmpRegIndexIt, ++oNeighbIt )
223 unsigned char tmpValue = tmpRegIndexIt.Get();
225 // Test current pixel: it is active ( on ) or not?
226 if( tmpValue == onTag )
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
232 // Test current pixel: it is a border pixel or an inner pixel?
233 bool bIsOnContour = false;
235 for (i = 0; i < neighborhoodSize; ++i)
237 // If at least one neighbour pixel is off the center pixel
238 // belongs to contour
239 if( oNeighbIt.GetPixel( i ) == backgroundTag )
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
252 // mark pixel value as a border pixel
253 tmpRegIndexIt.Set( borderTag );
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
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();
267 IndexType idx = tmpRegIndexIt.GetIndex() + *itIdx;
268 if( outputRegion.IsInside( idx ) )
269 { output->SetPixel( idx, backgroundValue ); }
273 propagQueue.push( tmpRegIndexIt.GetIndex() );
275 // now find all the border pixels
276 while ( !propagQueue.empty() )
278 // Extract pixel index from queue
279 IndexType currentIndex = propagQueue.front();
282 nit += currentIndex - nit.GetIndex();
284 for (i = 0; i < neighborhoodSize; ++i)
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
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 )
295 // Check if it is an inner or border neighbour pixel
296 // Get index of current neighbour pixel
297 IndexType neighbIndex = nit.GetIndex( i );
299 // Force location of neighbour iterator
300 nnit += neighbIndex - nnit.GetIndex();
302 bool bIsOnBorder = false;
304 for( j = 0; j < neighborhoodSize; ++j)
306 // If at least one neighbour pixel is off the center
307 // pixel belongs to border
308 if( nnit.GetPixel(j) == backgroundTag )
318 // neighbour pixel is a border pixel
321 nit.SetPixel( i, borderTag, status );
323 // check whether we could set the pixel. can only set
324 // the pixel if it is within the tmpimage
328 propagQueue.push( neighbIndex );
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();
338 IndexType idx = neighbIndex + *itIndex;
339 if( outputRegion.IsInside( idx ) )
340 { output->SetPixel( idx, backgroundValue ); }
346 // neighbour pixel is an inner pixel
348 nit.SetPixel( i, innerTag, status );
351 progress.CompletedPixel();
352 } // if( nit.GetPixel( i ) == onTag )
354 } // for (i = 0; i < neighborhoodSize; ++i)
356 } // while ( !propagQueue.empty() )
358 } // if( bIsOnCountour )
360 { tmpRegIndexIt.Set( innerTag ); }
362 progress.CompletedPixel();
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
373 // Deallocate tmpImage
374 tmpImage->Initialize();
377 // traverse structure of border and SE CCs, and paint output image
379 // Let's consider the the set of the ON elements of the input image as X.
381 // Let's consider the structuring element as B = {B0, B1, ..., Bn},
382 // where Bi denotes a connected component of B.
384 // Let's consider bi, i in [0,n], an arbitrary point of Bi.
387 // We use hence the next property in order to compute minkoswki
388 // addition ( which will be written (+) ):
390 // X (+) B = ( Xb0 UNION Xb1 UNION ... Xbn ) UNION ( BORDER(X) (+) B ),
392 // where Xbi is the set X translated with respect to vector bi :
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 )
397 // Define boundaries conditions
398 itk::ConstantBoundaryCondition<TOutputImage> obc;
399 obc.SetConstant( backgroundValue );
401 itk::NeighborhoodIterator<OutputImageType>
402 onit( kernel.GetRadius(), output, outputRegion );
403 onit.OverrideBoundaryCondition(&obc);
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();
415 // iterator on output image
416 itk::ImageRegionIteratorWithIndex<OutputImageType>
417 ouRegIndexIt(output, outputRegion );
418 ouRegIndexIt.GoToBegin();
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;
430 // Pad the input region by the kernel
431 inputRegionForThread.PadByRadius( kernel.GetRadius() );
432 inputRegionForThread.Crop( input->GetBufferedRegion() );
434 if( !this->m_BoundaryToForeground )
436 while( !ouRegIndexIt.IsAtEnd() )
438 // Retrieve index of current output pixel
439 IndexType currentIndex = ouRegIndexIt.GetIndex();
440 for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
443 IndexType translatedIndex = currentIndex - *vecIt;
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 )
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.
462 progress.CompletedPixel();
467 while( !ouRegIndexIt.IsAtEnd() )
469 IndexType currentIndex = ouRegIndexIt.GetIndex();
470 for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
472 IndexType translatedIndex = currentIndex - *vecIt;
474 if( inputRegionForThread.IsInside( translatedIndex )
475 //JV != foregroundValue
476 && input->GetPixel( translatedIndex ) == backgroundValue )
478 ouRegIndexIt.Set( backgroundValue );
484 progress.CompletedPixel();
488 // now, we must to restore the background values
489 itk::ImageRegionConstIterator<InputImageType> inIt( input, outputRegion );
491 for( inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt )
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();
506 * Standard "PrintSelf" method
508 template <class TInputImage, class TOutput, class TKernel>
510 ConditionalBinaryErodeImageFilter<TInputImage, TOutput, TKernel>
511 ::PrintSelf( std::ostream& os, itk::Indent indent) const
513 Superclass::PrintSelf( os, indent );
514 os << indent << "Dilate Value: " << static_cast<typename itk::NumericTraits<InputPixelType>::PrintType>( this->GetForegroundValue() ) << std::endl;
519 #endif //#define clitkConditionalBinaryErodeImageFilter_txx