1 #ifndef clitkConditionalBinaryDilateImageFilter_txx
2 #define clitkConditionalBinaryDilateImageFilter_txx
4 /* =================================================
5 * @file clitkConditionalBinaryDilateImageFilter.txx
11 ===================================================*/
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"
28 template <class TInputImage, class TOutputImage, class TKernel>
29 ConditionalBinaryDilateImageFilter<TInputImage, TOutputImage, TKernel>
30 ::ConditionalBinaryDilateImageFilter()
32 this->m_BoundaryToForeground = false;
36 template< class TInputImage, class TOutputImage, class TKernel>
38 ConditionalBinaryDilateImageFilter< TInputImage, TOutputImage, TKernel>
41 this->AllocateOutputs();
45 // Retrieve input and output pointers
46 typename OutputImageType::Pointer output = this->GetOutput();
47 typename InputImageType::ConstPointer input = this->GetInput();
49 // Get values from superclass
50 InputPixelType foregroundValue = this->GetForegroundValue();
51 InputPixelType backgroundValue = this->GetBackgroundValue();
52 KernelType kernel = this->GetKernel();
55 typename TInputImage::RegionType inputRegion = input->GetBufferedRegion();
56 typename TOutputImage::RegionType outputRegion = output->GetBufferedRegion();
58 // compute the size of the temp image. It is needed to create the progress
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)
70 padBy[i] = (padBy[i]>kernel.GetRadius(i) ? padBy[i] : kernel.GetRadius(i));
72 tmpRequestedRegion.PadByRadius( padBy );
73 tmpRequestedRegion.Crop( paddedInputRegion );
75 typename TInputImage::RegionType requiredInputRegion
76 = input->GetBufferedRegion();
77 requiredInputRegion.Crop( tmpRequestedRegion );
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() );
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 );
94 for( inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt )
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
103 { outIt.Set( static_cast<OutputPixelType>( value ) ); }
104 progress.CompletedPixel();
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();
114 // Define regions of temp image
115 tmpImage->SetRegions( tmpRequestedRegion );
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();
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;
137 if( this->m_BoundaryToForeground )
138 { tmpImage->FillBuffer( onTag ); }
140 { tmpImage->FillBuffer( backgroundTag ); }
142 // Iterators on input and tmp image
144 itk::ImageRegionConstIterator<TInputImage> iRegIt( input, requiredInputRegion );
145 // iterator on tmp image
146 itk::ImageRegionIterator<TempImageType> tmpRegIt( tmpImage, requiredInputRegion );
148 for( iRegIt.GoToBegin(), tmpRegIt.GoToBegin();
150 ++iRegIt, ++tmpRegIt )
152 OutputPixelType pxl = iRegIt.Get();
153 // JV == foregroundValue
154 if( pxl == foregroundValue )
155 { tmpRegIt.Set( onTag ); }
157 else if ( pxl == backgroundValue )
158 { tmpRegIt.Set( bgTag ); }
161 // by default if it is not foreground, consider
163 tmpRegIt.Set( backgroundTag );
165 progress.CompletedPixel();
170 // Border tracking and encoding
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 );
181 itk::ConstNeighborhoodIterator<TempImageType>
182 oNeighbIt( radius, tmpImage, tmpRequestedRegion );
184 // Define boundaries conditions
185 itk::ConstantBoundaryCondition<TempImageType> cbc;
186 cbc.SetConstant( backgroundTag );
187 oNeighbIt.OverrideBoundaryCondition(&cbc);
189 unsigned int neighborhoodSize = oNeighbIt.Size();
190 unsigned int centerPixelCode = neighborhoodSize / 2;
192 std::queue<IndexType> propagQueue;
194 // Neighborhood iterators used to track the surface.
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);
208 itk::ConstNeighborhoodIterator<TempImageType>
209 nnit( radius, tmpImage, tmpRequestedRegion );
210 nnit.OverrideBoundaryCondition(&cbc);
213 for( tmpRegIndexIt.GoToBegin(), oNeighbIt.GoToBegin();
214 !tmpRegIndexIt.IsAtEnd();
215 ++tmpRegIndexIt, ++oNeighbIt )
217 // Test current pixel: it is active ( on ) or not?
218 if( tmpRegIndexIt.Get() == onTag )
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
224 // Test current pixel: it is a border pixel or an inner pixel?
225 bool bIsOnContour = false;
227 for (i = 0; i < neighborhoodSize; ++i)
229 // If at least one neighbour pixel is off the center pixel
230 // belongs to contour
232 if( oNeighbIt.GetPixel( i ) == bgTag )
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
245 // mark pixel value as a border pixel
246 tmpRegIndexIt.Set( borderTag );
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
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();
260 IndexType idx = tmpRegIndexIt.GetIndex() + *itIdx;
261 if( outputRegion.IsInside( idx ) )
263 // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
264 if (input->GetPixel(idx)==backgroundValue)
265 output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
270 propagQueue.push( tmpRegIndexIt.GetIndex() );
272 // now find all the border pixels
273 while ( !propagQueue.empty() )
275 // Extract pixel index from queue
276 IndexType currentIndex = propagQueue.front();
279 nit += currentIndex - nit.GetIndex();
281 for (i = 0; i < neighborhoodSize; ++i)
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
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 )
292 // Check if it is an inner or border neighbour pixel
293 // Get index of current neighbour pixel
294 IndexType neighbIndex = nit.GetIndex( i );
296 // Force location of neighbour iterator
297 nnit += neighbIndex - nnit.GetIndex();
299 bool bIsOnBorder = false;
301 for( j = 0; j < neighborhoodSize; ++j)
303 // If at least one neighbour pixel is off the center
304 // pixel belongs to border
305 if( nnit.GetPixel(j) == bgTag )
315 // neighbour pixel is a border pixel
318 nit.SetPixel( i, borderTag, status );
320 // check whether we could set the pixel. can only set
321 // the pixel if it is within the tmpimage
325 propagQueue.push( neighbIndex );
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();
335 IndexType idx = neighbIndex + *itIndex;
336 if( outputRegion.IsInside( idx ) )
338 // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
339 if (input->GetPixel(idx)==backgroundValue)
340 output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
347 // neighbour pixel is an inner pixel
349 nit.SetPixel( i, innerTag, status );
352 progress.CompletedPixel();
353 } // if( nit.GetPixel( i ) == onTag )
355 } // for (i = 0; i < neighborhoodSize; ++i)
357 } // while ( !propagQueue.empty() )
359 } // if( bIsOnCountour )
361 { tmpRegIndexIt.Set( innerTag ); }
363 } // if( tmpRegIndexIt.Get() == onTag )
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
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( static_cast<OutputPixelType> ( 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 if( !inputRegionForThread.IsInside( translatedIndex )
453 || ( (input->GetPixel( translatedIndex ) == foregroundValue )
455 && (input->GetPixel( currentIndex ) == backgroundValue ) ) )
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.
464 progress.CompletedPixel();
469 while( !ouRegIndexIt.IsAtEnd() )
471 IndexType currentIndex = ouRegIndexIt.GetIndex();
472 for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
474 IndexType translatedIndex = currentIndex - *vecIt;
476 if( inputRegionForThread.IsInside( translatedIndex )
477 // JV == foregroundValue
478 && ( (input->GetPixel( translatedIndex ) == foregroundValue )
480 && (input->GetPixel( currentIndex ) == backgroundValue ) ) )
483 ouRegIndexIt.Set( static_cast<OutputPixelType> ( foregroundValue ) );
489 progress.CompletedPixel();
496 * Standard "PrintSelf" method
498 template <class TInputImage, class TOutput, class TKernel>
500 ConditionalBinaryDilateImageFilter<TInputImage, TOutput, TKernel>
501 ::PrintSelf( std::ostream& os, itk::Indent indent) const
503 Superclass::PrintSelf( os, indent );
504 os << indent << "Dilate Value: " << static_cast<typename itk::NumericTraits<InputPixelType>::PrintType>( this->GetForegroundValue() ) << std::endl;
507 } // end namespace itk