--- /dev/null
+#ifndef clitkConditionalBinaryDilateImageFilter_h
+#define clitkConditionalBinaryDilateImageFilter_h
+
+/* =================================================
+ * @file clitkConditionalBinaryDilateImageFilter.h
+ * @author
+ * @date
+ *
+ * @brief
+ *
+ ===================================================*/
+
+// clitk include
+#include "clitkIO.h"
+#include "clitkCommon.h"
+
+// itk include
+#include <vector>
+#include <queue>
+#include "itkBinaryMorphologyImageFilter.h"
+#include "itkImage.h"
+#include "itkNumericTraits.h"
+#include "itkNeighborhoodIterator.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhood.h"
+#include "itkImageBoundaryCondition.h"
+#include "itkImageRegionIterator.h"
+#include "itkConceptChecking.h"
+
+namespace clitk
+{
+
+ template <class TInputImage, class TOutputImage, class TKernel>
+ class ITK_EXPORT ConditionalBinaryDilateImageFilter :
+ public itk::BinaryMorphologyImageFilter< TInputImage, TOutputImage, TKernel >
+ {
+ public:
+ /** Extract dimension from input and output image. */
+ itkStaticConstMacro(InputImageDimension, unsigned int,
+ TInputImage::ImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int,
+ TOutputImage::ImageDimension);
+
+ /** Extract the dimension of the kernel */
+ itkStaticConstMacro(KernelDimension, unsigned int,
+ TKernel::NeighborhoodDimension);
+
+ /** Convenient typedefs for simplifying declarations. */
+ typedef TInputImage InputImageType;
+ typedef TOutputImage OutputImageType;
+ typedef TKernel KernelType;
+
+
+ /** Standard class typedefs. */
+ typedef ConditionalBinaryDilateImageFilter Self;
+ typedef itk::BinaryMorphologyImageFilter<InputImageType, OutputImageType,
+ KernelType> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(ConditionalBinaryDilateImageFilter, BinaryMorphologyImageFilter);
+
+ /** Kernel (structuring element) iterator. */
+ typedef typename KernelType::ConstIterator KernelIteratorType;
+
+ /** Image typedef support. */
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef typename itk::NumericTraits<InputPixelType>::RealType InputRealType;
+ typedef typename InputImageType::OffsetType OffsetType;
+ typedef typename InputImageType::IndexType IndexType;
+
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+ typedef typename InputImageType::SizeType InputSizeType;
+
+ /** Set the value in the image to consider as "foreground". Defaults to
+ * maximum value of PixelType. This is an alias to the
+ * ForegroundValue in the superclass. */
+ void SetDilateValue(const InputPixelType& value)
+ { this->SetForegroundValue( value ); }
+
+ /** Get the value in the image considered as "foreground". Defaults to
+ * maximum value of PixelType. This is an alias to the
+ * ForegroundValue in the superclass. */
+ InputPixelType GetDilateValue() const
+ { return this->GetForegroundValue(); }
+
+ protected:
+ ConditionalBinaryDilateImageFilter();
+ virtual ~ConditionalBinaryDilateImageFilter(){}
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+ void GenerateData();
+
+ // type inherited from the superclass
+ typedef typename Superclass::NeighborIndexContainer NeighborIndexContainer;
+
+ private:
+ ConditionalBinaryDilateImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ };
+
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkConditionalBinaryDilateImageFilter.txx"
+#endif
+
+#endif
--- /dev/null
+#ifndef clitkConditionalBinaryDilateImageFilter_txx
+#define clitkConditionalBinaryDilateImageFilter_txx
+
+/* =================================================
+ * @file clitkConditionalBinaryDilateImageFilter.txx
+ * @author
+ * @date
+ *
+ * @brief
+ *
+ ===================================================*/
+
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhoodIterator.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkNeighborhoodInnerProduct.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkConstantBoundaryCondition.h"
+#include "itkOffset.h"
+#include "itkProgressReporter.h"
+#include "itkBinaryDilateImageFilter.h"
+
+namespace clitk
+{
+
+ template <class TInputImage, class TOutputImage, class TKernel>
+ ConditionalBinaryDilateImageFilter<TInputImage, TOutputImage, TKernel>
+ ::ConditionalBinaryDilateImageFilter()
+ {
+ this->m_BoundaryToForeground = false;
+ }
+
+
+ template< class TInputImage, class TOutputImage, class TKernel>
+ void
+ ConditionalBinaryDilateImageFilter< TInputImage, TOutputImage, TKernel>
+ ::GenerateData()
+ {
+ this->AllocateOutputs();
+
+ unsigned int i,j;
+
+ // Retrieve input and output pointers
+ typename OutputImageType::Pointer output = this->GetOutput();
+ typename InputImageType::ConstPointer input = this->GetInput();
+
+ // Get values from superclass
+ InputPixelType foregroundValue = this->GetForegroundValue();
+ InputPixelType backgroundValue = this->GetBackgroundValue();
+ KernelType kernel = this->GetKernel();
+ InputSizeType radius;
+ radius.Fill(1);
+ typename TInputImage::RegionType inputRegion = input->GetBufferedRegion();
+ typename TOutputImage::RegionType outputRegion = output->GetBufferedRegion();
+
+ // compute the size of the temp image. It is needed to create the progress
+ // reporter.
+ // The tmp image needs to be large enough to support:
+ // 1. The size of the structuring element
+ // 2. The size of the connectivity element (typically one)
+ typename TInputImage::RegionType tmpRequestedRegion = outputRegion;
+ typename TInputImage::RegionType paddedInputRegion
+ = input->GetBufferedRegion();
+ paddedInputRegion.PadByRadius( radius ); // to support boundary values
+ InputSizeType padBy = radius;
+ for (i=0; i < KernelDimension; ++i)
+ {
+ padBy[i] = (padBy[i]>kernel.GetRadius(i) ? padBy[i] : kernel.GetRadius(i));
+ }
+ tmpRequestedRegion.PadByRadius( padBy );
+ tmpRequestedRegion.Crop( paddedInputRegion );
+
+ typename TInputImage::RegionType requiredInputRegion
+ = input->GetBufferedRegion();
+ requiredInputRegion.Crop( tmpRequestedRegion );
+
+ // Support progress methods/callbacks
+ // Setup a progress reporter. We have 4 stages to the algorithm so
+ // pretend we have 4 times the number of pixels
+ itk::ProgressReporter progress(this, 0,
+ outputRegion.GetNumberOfPixels() * 2
+ + tmpRequestedRegion.GetNumberOfPixels()
+ + requiredInputRegion.GetNumberOfPixels() );
+
+ // Allocate and reset output. We copy the input to the output,
+ // except for pixels with DilateValue. These pixels are initially
+ // replaced with BackgroundValue and potentially replaced later with
+ // DilateValue as the Minkowski sums are performed.
+ itk::ImageRegionIterator<OutputImageType> outIt( output, outputRegion );
+ itk::ImageRegionConstIterator<InputImageType> inIt( input, outputRegion );
+
+ for( inIt.GoToBegin(), outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt, ++inIt )
+ {
+ InputPixelType value = inIt.Get();
+ // replace foreground pixels with the default background
+ //JV == foregroundValue
+ if (value == backgroundValue)
+ { outIt.Set( static_cast<OutputPixelType> ( backgroundValue ) ); }
+ // keep all of the original background values intact
+ else
+ { outIt.Set( static_cast<OutputPixelType>( value ) ); }
+ progress.CompletedPixel();
+ }
+
+
+ // Create the temp image for surface encoding
+ // The temp image size is equal to the output requested region for thread
+ // padded by max( connectivity neighborhood radius, SE kernel radius ).
+ typedef itk::Image<unsigned char, TInputImage::ImageDimension> TempImageType;
+ typename TempImageType::Pointer tmpImage = TempImageType::New();
+
+ // Define regions of temp image
+ tmpImage->SetRegions( tmpRequestedRegion );
+
+ // Allocation.
+ // Pay attention to the fact that here, the output is still not
+ // allocated (so no extra memory needed for tmp image, if you
+ // consider that you reserve som memory space for output)
+ tmpImage->Allocate();
+
+ // First Stage
+ // Copy the input image to the tmp image.
+ // Tag the tmp Image.
+ // zero means background
+ // one means pixel on but not treated
+ // two means border pixel
+ // three means inner pixel
+ static const unsigned char backgroundTag = 0;
+ static const unsigned char onTag = 1;
+ static const unsigned char borderTag = 2;
+ static const unsigned char innerTag = 3;
+ // JV add the tag in which can be dilated
+ static const unsigned char bgTag = 4;
+
+ if( this->m_BoundaryToForeground )
+ { tmpImage->FillBuffer( onTag ); }
+ else
+ { tmpImage->FillBuffer( backgroundTag ); }
+
+ // Iterators on input and tmp image
+ // iterator on input
+ itk::ImageRegionConstIterator<TInputImage> iRegIt( input, requiredInputRegion );
+ // iterator on tmp image
+ itk::ImageRegionIterator<TempImageType> tmpRegIt( tmpImage, requiredInputRegion );
+
+ for( iRegIt.GoToBegin(), tmpRegIt.GoToBegin();
+ !tmpRegIt.IsAtEnd();
+ ++iRegIt, ++tmpRegIt )
+ {
+ OutputPixelType pxl = iRegIt.Get();
+ // JV == foregroundValue
+ if( pxl == foregroundValue )
+ { tmpRegIt.Set( onTag ); }
+ // JV add the bgTag
+ else if ( pxl == backgroundValue )
+ { tmpRegIt.Set( bgTag ); }
+ else
+ {
+ // by default if it is not foreground, consider
+ // it as background
+ tmpRegIt.Set( backgroundTag );
+ }
+ progress.CompletedPixel();
+ }
+
+
+ // Second stage
+ // Border tracking and encoding
+
+ // Need to record index, use an iterator with index
+ // Define iterators that will traverse the OUTPUT requested region
+ // for thread and not the padded one. The tmp image has been padded
+ // because in that way we will take care carefully at boundary
+ // pixels of output requested region. Take care means that we will
+ // check if a boundary pixel is or not a border pixel.
+ itk::ImageRegionIteratorWithIndex<TempImageType>
+ tmpRegIndexIt( tmpImage, tmpRequestedRegion );
+
+ itk::ConstNeighborhoodIterator<TempImageType>
+ oNeighbIt( radius, tmpImage, tmpRequestedRegion );
+
+ // Define boundaries conditions
+ itk::ConstantBoundaryCondition<TempImageType> cbc;
+ cbc.SetConstant( backgroundTag );
+ oNeighbIt.OverrideBoundaryCondition(&cbc);
+
+ unsigned int neighborhoodSize = oNeighbIt.Size();
+ unsigned int centerPixelCode = neighborhoodSize / 2;
+
+ std::queue<IndexType> propagQueue;
+
+ // Neighborhood iterators used to track the surface.
+ //
+ // Note the region specified for the first neighborhood iterator is
+ // the requested region for the tmp image not the output image. This
+ // is necessary because the NeighborhoodIterator relies on the
+ // specified region to determine if you will ever query a boundary
+ // condition pixel. Since we call SetLocation on the neighbor of a
+ // specified pixel, we have to set the region for the interator to
+ // include any pixel we may set our location to.
+ itk::NeighborhoodIterator<TempImageType>
+ nit( radius, tmpImage, tmpRequestedRegion );
+ nit.OverrideBoundaryCondition(&cbc);
+ nit.GoToBegin();
+
+ itk::ConstNeighborhoodIterator<TempImageType>
+ nnit( radius, tmpImage, tmpRequestedRegion );
+ nnit.OverrideBoundaryCondition(&cbc);
+ nnit.GoToBegin();
+
+ for( tmpRegIndexIt.GoToBegin(), oNeighbIt.GoToBegin();
+ !tmpRegIndexIt.IsAtEnd();
+ ++tmpRegIndexIt, ++oNeighbIt )
+ {
+ // Test current pixel: it is active ( on ) or not?
+ if( tmpRegIndexIt.Get() == onTag )
+ {
+ // The current pixel has not been treated previously. That
+ // means that we do not know that it is an inner pixel of a
+ // border pixel.
+
+ // Test current pixel: it is a border pixel or an inner pixel?
+ bool bIsOnContour = false;
+
+ for (i = 0; i < neighborhoodSize; ++i)
+ {
+ // If at least one neighbour pixel is off the center pixel
+ // belongs to contour
+ //JV verify bg
+ if( oNeighbIt.GetPixel( i ) == bgTag )
+ {
+ bIsOnContour = true;
+ break;
+ }
+ }
+
+ if( bIsOnContour )
+ {
+ // center pixel is a border pixel and due to the parsing, it is also
+ // a pixel which belongs to a new border connected component
+ // Now we will parse this border thanks to a burn procedure
+
+ // mark pixel value as a border pixel
+ tmpRegIndexIt.Set( borderTag );
+
+ // add it to border container.
+ // its code is center pixel code because it is the first pixel
+ // of the connected component border
+
+ // paint the structuring element
+ typename NeighborIndexContainer::const_iterator itIdx;
+ NeighborIndexContainer& idxDifferenceSet
+ = this->GetDifferenceSet( centerPixelCode );
+ for( itIdx = idxDifferenceSet.begin();
+ itIdx != idxDifferenceSet.end();
+ ++itIdx )
+ {
+ IndexType idx = tmpRegIndexIt.GetIndex() + *itIdx;
+ if( outputRegion.IsInside( idx ) )
+ {
+ // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
+ if (input->GetPixel(idx)==backgroundValue)
+ output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
+ }
+ }
+
+ // add it to queue
+ propagQueue.push( tmpRegIndexIt.GetIndex() );
+
+ // now find all the border pixels
+ while ( !propagQueue.empty() )
+ {
+ // Extract pixel index from queue
+ IndexType currentIndex = propagQueue.front();
+ propagQueue.pop();
+
+ nit += currentIndex - nit.GetIndex();
+
+ for (i = 0; i < neighborhoodSize; ++i)
+ {
+ // If pixel has not been already treated and it is a pixel
+ // on, test if it is an inner pixel or a border pixel
+
+ // Remark: all the pixels outside the image are set to
+ // backgroundTag thanks to boundary conditions. That means that if
+ // we enter in the next if-statement we are sure that the
+ // current neighbour pixel is in the image
+ if( nit.GetPixel( i ) == onTag )
+ {
+ // Check if it is an inner or border neighbour pixel
+ // Get index of current neighbour pixel
+ IndexType neighbIndex = nit.GetIndex( i );
+
+ // Force location of neighbour iterator
+ nnit += neighbIndex - nnit.GetIndex();
+
+ bool bIsOnBorder = false;
+
+ for( j = 0; j < neighborhoodSize; ++j)
+ {
+ // If at least one neighbour pixel is off the center
+ // pixel belongs to border
+ if( nnit.GetPixel(j) == bgTag )
+ {
+ bIsOnBorder = true;
+ break;
+ }
+ }
+
+
+ if( bIsOnBorder )
+ {
+ // neighbour pixel is a border pixel
+ // mark it
+ bool status;
+ nit.SetPixel( i, borderTag, status );
+
+ // check whether we could set the pixel. can only set
+ // the pixel if it is within the tmpimage
+ if (status)
+ {
+ // add it to queue
+ propagQueue.push( neighbIndex );
+
+ // paint the structuring element
+ typename NeighborIndexContainer::const_iterator itIndex;
+ NeighborIndexContainer& indexDifferenceSet
+ = this->GetDifferenceSet( i );
+ for( itIndex = indexDifferenceSet.begin();
+ itIndex != indexDifferenceSet.end();
+ ++itIndex )
+ {
+ IndexType idx = neighbIndex + *itIndex;
+ if( outputRegion.IsInside( idx ) )
+ {
+ // JV output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) ); }
+ if (input->GetPixel(idx)==backgroundValue)
+ output->SetPixel( idx, static_cast<OutputPixelType> ( foregroundValue ) );
+ }
+ }
+ }
+ }
+ else
+ {
+ // neighbour pixel is an inner pixel
+ bool status;
+ nit.SetPixel( i, innerTag, status );
+ }
+
+ progress.CompletedPixel();
+ } // if( nit.GetPixel( i ) == onTag )
+
+ } // for (i = 0; i < neighborhoodSize; ++i)
+
+ } // while ( !propagQueue.empty() )
+
+ } // if( bIsOnCountour )
+ else
+ { tmpRegIndexIt.Set( innerTag ); }
+
+ } // if( tmpRegIndexIt.Get() == onTag )
+ else
+ { progress.CompletedPixel(); }
+ // Here, the pixel is a background pixel ( value at 0 ) or an
+ // already treated pixel:
+ // 2 for border pixel, 3 for inner pixel
+
+ }
+
+
+ // Deallocate tmpImage
+ tmpImage->Initialize();
+
+ // Third Stage
+ // traverse structure of border and SE CCs, and paint output image
+
+ // Let's consider the the set of the ON elements of the input image as X.
+ //
+ // Let's consider the structuring element as B = {B0, B1, ..., Bn},
+ // where Bi denotes a connected component of B.
+ //
+ // Let's consider bi, i in [0,n], an arbitrary point of Bi.
+ //
+
+ // We use hence the next property in order to compute minkoswki
+ // addition ( which will be written (+) ):
+ //
+ // X (+) B = ( Xb0 UNION Xb1 UNION ... Xbn ) UNION ( BORDER(X) (+) B ),
+ //
+ // where Xbi is the set X translated with respect to vector bi :
+ //
+ // Xbi = { x + bi, x belongs to X } where BORDER(X) is the extracted
+ // border of X ( 8 connectivity in 2D, 26 in 3D )
+
+ // Define boundaries conditions
+ itk::ConstantBoundaryCondition<TOutputImage> obc;
+ obc.SetConstant( static_cast<OutputPixelType> ( backgroundValue ) );
+
+ itk::NeighborhoodIterator<OutputImageType>
+ onit( kernel.GetRadius(), output, outputRegion );
+ onit.OverrideBoundaryCondition(&obc);
+ onit.GoToBegin();
+
+
+ // Paint input image translated with respect to the SE CCs vectors
+ // --> "( Xb0 UNION Xb1 UNION ... Xbn )"
+ typename Superclass::ComponentVectorConstIterator vecIt;
+ typename Superclass::ComponentVectorConstIterator vecBeginIt
+ = this->KernelCCVectorBegin();
+ typename Superclass::ComponentVectorConstIterator vecEndIt
+ = this->KernelCCVectorEnd();
+
+ // iterator on output image
+ itk::ImageRegionIteratorWithIndex<OutputImageType>
+ ouRegIndexIt(output, outputRegion );
+ ouRegIndexIt.GoToBegin();
+
+ // InputRegionForThread is the output region for thread padded by
+ // kerne lradius We must traverse this padded region because some
+ // border pixel in the added band ( the padded band is the region
+ // added after padding ) may be responsible to the painting of some
+ // pixel in the non padded region. This happens typically when a
+ // non centered SE is used, a kind of shift is done on the "on"
+ // pixels of image. Consequently some pixels in the added band can
+ // appear in the current region for thread due to shift effect.
+ typename InputImageType::RegionType inputRegionForThread = outputRegion;
+
+ // Pad the input region by the kernel
+ inputRegionForThread.PadByRadius( kernel.GetRadius() );
+ inputRegionForThread.Crop( input->GetBufferedRegion() );
+
+ if( this->m_BoundaryToForeground )
+ {
+ while( !ouRegIndexIt.IsAtEnd() )
+ {
+ // Retrieve index of current output pixel
+ IndexType currentIndex = ouRegIndexIt.GetIndex();
+ for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
+ {
+ // Translate
+ IndexType translatedIndex = currentIndex - *vecIt;
+
+ // translated index now is an index in input image in the
+ // output requested region padded. Theoretically, this translated
+ // index must be inside the padded region.
+ // If the pixel in the input image at the translated index
+ // has a value equal to the dilate one, this means
+ // that the output pixel at currentIndex will be on in the output.
+ if( !inputRegionForThread.IsInside( translatedIndex )
+ // JV == foreground
+ || ( (input->GetPixel( translatedIndex ) == foregroundValue )
+ // JV add bg
+ && (input->GetPixel( currentIndex ) == backgroundValue ) ) )
+ {
+ ouRegIndexIt.Set( static_cast<OutputPixelType> ( foregroundValue ) );
+ break; // Do not need to examine other offsets because at least one
+ // input pixel has been translated on current output pixel.
+ }
+ }
+
+ ++ouRegIndexIt;
+ progress.CompletedPixel();
+ }
+ }
+ else
+ {
+ while( !ouRegIndexIt.IsAtEnd() )
+ {
+ IndexType currentIndex = ouRegIndexIt.GetIndex();
+ for( vecIt = vecBeginIt; vecIt != vecEndIt; ++vecIt )
+ {
+ IndexType translatedIndex = currentIndex - *vecIt;
+
+ if( inputRegionForThread.IsInside( translatedIndex )
+ // JV == foregroundValue
+ && ( (input->GetPixel( translatedIndex ) == foregroundValue )
+ // JV add bg
+ && (input->GetPixel( currentIndex ) == backgroundValue ) ) )
+
+ {
+ ouRegIndexIt.Set( static_cast<OutputPixelType> ( foregroundValue ) );
+ break;
+ }
+ }
+
+ ++ouRegIndexIt;
+ progress.CompletedPixel();
+ }
+ }
+ }
+
+
+ /**
+ * Standard "PrintSelf" method
+ */
+ template <class TInputImage, class TOutput, class TKernel>
+ void
+ ConditionalBinaryDilateImageFilter<TInputImage, TOutput, TKernel>
+ ::PrintSelf( std::ostream& os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf( os, indent );
+ os << indent << "Dilate Value: " << static_cast<typename itk::NumericTraits<InputPixelType>::PrintType>( this->GetForegroundValue() ) << std::endl;
+ }
+
+} // end namespace itk
+
+#endif
--- /dev/null
+#ifndef clitkDecomposeAndReconstructImageFilter_h
+#define clitkDecomposeAndReconstructImageFilter_h
+
+/* =================================================
+ * @file clitkDecomposeAndReconstructImageFilter.h
+ * @author
+ * @date
+ *
+ * @brief
+ *
+ ===================================================*/
+
+
+// clitk include
+#include "clitkIO.h"
+#include "clitkCommon.h"
+#include "clitkDecomposeThroughErosionImageFilter.h"
+#include "clitkReconstructThroughDilationImageFilter.h"
+
+//itk include
+#include "itkImageToImageFilter.h"
+#include "itkRelabelComponentImageFilter.h"
+
+
+namespace clitk
+{
+
+ template <class InputImageType, class OutputImageType>
+ class ITK_EXPORT DecomposeAndReconstructImageFilter :
+ public itk::ImageToImageFilter<InputImageType, OutputImageType>
+ {
+ public:
+ //----------------------------------------
+ // ITK
+ //----------------------------------------
+ typedef DecomposeAndReconstructImageFilter Self;
+ typedef itk::ImageToImageFilter<InputImageType, OutputImageType> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ // Method for creation through the object factory
+ itkNewMacro(Self);
+
+ // Run-time type information (and related methods)
+ itkTypeMacro( DecomposeAndReconstructImageFilter, ImageToImageFilter );
+
+ /** Dimension of the domain space. */
+ itkStaticConstMacro(InputImageDimension, unsigned int, Superclass::InputImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int, Superclass::OutputImageDimension);
+
+ //----------------------------------------
+ // Typedefs
+ //----------------------------------------
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+ typedef int InternalPixelType;
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef typename InputImageType::SizeType SizeType;
+
+
+ //----------------------------------------
+ // Set & Get
+ //----------------------------------------
+ itkBooleanMacro(Verbose);
+ itkSetMacro( Verbose, bool);
+ itkGetConstReferenceMacro( Verbose, bool);
+ itkBooleanMacro(FullyConnected);
+ itkSetMacro( FullyConnected, bool);
+ itkGetConstReferenceMacro( FullyConnected, bool);
+ void SetRadius ( const SizeType& s){ m_Radius=s; this->Modified();}
+ SizeType GetRadius(void){return m_Radius;}
+ itkSetMacro( MaximumNumberOfLabels, unsigned int);
+ itkGetConstMacro( MaximumNumberOfLabels, unsigned int);
+ itkSetMacro( BackgroundValue, InternalPixelType);
+ itkGetConstMacro( BackgroundValue, InternalPixelType);
+ itkSetMacro( ForegroundValue, InternalPixelType);
+ itkGetConstMacro( ForegroundValue, InternalPixelType);
+ itkSetMacro( NumberOfNewLabels, unsigned int);
+ itkGetConstMacro( NumberOfNewLabels, unsigned int);
+ itkSetMacro( MinimumObjectSize, unsigned int);
+ itkGetConstMacro( MinimumObjectSize, unsigned int);
+ itkSetMacro( MinimumNumberOfIterations, unsigned int);
+ itkGetConstMacro( MinimumNumberOfIterations, unsigned int);
+ // // Convenience macro's: Built-in
+ // itkBooleanMacro (flag); //FlagOn FlagOff
+ // itkGetMacro(name, type);
+ // itkSetMacro(name, type);
+ // itkSetConstMacro( name, type);
+ // itkGetConstMacro( name, type);
+ // itkSetConstReferenceMacro(name, type);
+ // itkGetConstReferenceMacro(name, type);
+ // // Convenience macro's: Smartpointers
+ // itkSetObjectMacro(name, type);
+ // itkGetObjectMacro(name, type);
+ // itkSetConstObjectMacro(name, type);
+ // itkGetConstObjectMacro(name, type);
+ // itkSetConstReferenceObjectMacro(name, type);
+ // itkSetConstReference(name, type);
+
+
+ protected:
+
+ //----------------------------------------
+ // Constructor & Destructor
+ //----------------------------------------
+ DecomposeAndReconstructImageFilter();
+ ~DecomposeAndReconstructImageFilter() {};
+
+ //----------------------------------------
+ // Update
+ //----------------------------------------
+ // Generate Data
+ void GenerateData(void);
+
+ // // Threaded Generate Data
+ // void BeforeThreadedGenerateData(void );
+ // void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId );
+ // void AfterThreadedGenerateData(void );
+ // // Override defaults
+ // virtual void GenerateInputRequestedRegion();
+ // virtual void GenerateOutputInformation (void);
+ // virtual void EnlargeOutputRequestedRegion(DataObject *data);
+ void AllocateOutputs(){;}
+ //----------------------------------------
+ // Data members
+ //----------------------------------------
+ bool m_Verbose;
+ SizeType m_Radius;
+ unsigned int m_NumberOfNewLabels;
+ bool m_FullyConnected;
+ InputPixelType m_BackgroundValue;
+ InputPixelType m_ForegroundValue;
+ unsigned int m_MaximumNumberOfLabels;
+ unsigned int m_MinimumObjectSize;
+ unsigned int m_MinimumNumberOfIterations;
+ };
+
+
+} // end namespace clitk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkDecomposeAndReconstructImageFilter.txx"
+#endif
+
+#endif // #define clitkDecomposeAndReconstructImageFilter_h
+
+
--- /dev/null
+#ifndef clitkDecomposeAndReconstructImageFilter_txx
+#define clitkDecomposeAndReconstructImageFilter_txx
+
+/* =================================================
+ * @file clitkDecomposeAndReconstructImageFilter.txx
+ * @author
+ * @date
+ *
+ * @brief
+ *
+ ===================================================*/
+
+
+namespace clitk
+{
+
+ //-------------------------------------------------------------------
+ // Update with the number of dimensions
+ //-------------------------------------------------------------------
+ template<class InputImageType, class OutputImageType>
+ DecomposeAndReconstructImageFilter<InputImageType, OutputImageType>::DecomposeAndReconstructImageFilter()
+ {
+ m_Verbose=false;
+ for (unsigned int i=0; i<InputImageDimension; i++)
+ m_Radius[i]=1;
+ m_NumberOfNewLabels=1;
+ m_FullyConnected=true;
+ m_BackgroundValue=0;
+ m_ForegroundValue=1;
+ m_MaximumNumberOfLabels=10;
+ m_MinimumObjectSize=10;
+ m_MinimumNumberOfIterations=1;
+
+ }
+
+
+ //-------------------------------------------------------------------
+ // Update with the number of dimensions and the pixeltype
+ //-------------------------------------------------------------------
+ template <class InputImageType, class OutputImageType>
+ void
+ DecomposeAndReconstructImageFilter<InputImageType, OutputImageType>::GenerateData()
+ {
+
+
+ // Internal type
+ typedef itk::Image<InternalPixelType, InputImageDimension> InternalImageType;
+
+ // Filters used
+ typedef clitk::DecomposeThroughErosionImageFilter<InputImageType, InternalImageType> DecomposeThroughErosionImageFilterType;
+ typedef clitk::ReconstructThroughDilationImageFilter<InternalImageType, InputImageType> ReconstructThroughDilationImageFilterType;
+
+ // Erode
+ typename DecomposeThroughErosionImageFilterType::Pointer erosionFilter=DecomposeThroughErosionImageFilterType::New();
+ erosionFilter->SetInput(this->GetInput());
+ erosionFilter->SetVerbose(m_Verbose);
+ erosionFilter->SetFullyConnected(m_FullyConnected);
+ erosionFilter->SetRadius(m_Radius);
+ erosionFilter->SetNumberOfNewLabels(m_NumberOfNewLabels);
+ erosionFilter->SetMinimumObjectSize(m_MinimumObjectSize);
+ erosionFilter->SetMinimumNumberOfIterations(m_MinimumNumberOfIterations);
+ erosionFilter->Update();
+
+ // Reconstruct
+ typename ReconstructThroughDilationImageFilterType::Pointer reconstructionFilter =ReconstructThroughDilationImageFilterType::New();
+ reconstructionFilter->SetInput(erosionFilter->GetOutput());
+ reconstructionFilter->SetVerbose(m_Verbose);
+ reconstructionFilter->SetRadius(m_Radius);
+ reconstructionFilter->SetMaximumNumberOfLabels(m_MaximumNumberOfLabels);
+ reconstructionFilter->SetBackgroundValue(m_BackgroundValue);
+ reconstructionFilter->SetForegroundValue(m_ForegroundValue);
+ reconstructionFilter->Update();
+
+ // Output
+ this->SetNthOutput(0,reconstructionFilter->GetOutput());
+ }
+
+
+}//end clitk
+
+#endif //#define clitkDecomposeAndReconstructImageFilter_txx
--- /dev/null
+#ifndef clitkDecomposeThroughErosionImageFilter_h
+#define clitkDecomposeThroughErosionImageFilter_h
+
+/* =================================================
+ * @file clitkDecomposeThroughErosionImageFilter.h
+ * @author
+ * @date
+ *
+ * @brief
+ *
+ ===================================================*/
+
+
+// clitk include
+#include "clitkIO.h"
+#include "clitkCommon.h"
+#include "clitkSetBackgroundImageFilter.h"
+
+//itk include
+#include "itkImageToImageFilter.h"
+#include "itkBinaryThresholdImageFilter.h"
+#include "itkBinaryErodeImageFilter.h"
+#include "itkStatisticsImageFilter.h"
+#include "itkConnectedComponentImageFilter.h"
+#include "itkCastImageFilter.h"
+#include "itkBinaryBallStructuringElement.h"
+#include "itkRelabelComponentImageFilter.h"
+
+namespace clitk
+{
+
+ template <class InputImageType, class OutputImageType>
+ class ITK_EXPORT DecomposeThroughErosionImageFilter :
+ public itk::ImageToImageFilter<InputImageType, OutputImageType>
+ {
+ public:
+ //----------------------------------------
+ // ITK
+ //----------------------------------------
+ typedef DecomposeThroughErosionImageFilter Self;
+ typedef itk::ImageToImageFilter<InputImageType, OutputImageType> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ // Method for creation through the object factory
+ itkNewMacro(Self);
+
+ // Run-time type information (and related methods)
+ itkTypeMacro( DecomposeThroughErosionImageFilter, ImageToImageFilter );
+
+ /** Dimension of the domain space. */
+ itkStaticConstMacro(InputImageDimension, unsigned int, Superclass::InputImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int, Superclass::OutputImageDimension);
+
+ //----------------------------------------
+ // Typedefs
+ //----------------------------------------
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+ typedef int InternalPixelType;
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef typename InputImageType::SizeType SizeType;
+
+ //----------------------------------------
+ // Set & Get
+ //----------------------------------------
+ itkBooleanMacro(Verbose);
+ itkSetMacro( Verbose, bool);
+ itkGetConstReferenceMacro( Verbose, bool);
+ itkBooleanMacro(FullyConnected);
+ itkSetMacro( FullyConnected, bool);
+ itkGetConstReferenceMacro( FullyConnected, bool);
+ itkSetMacro( Lower, InputPixelType);
+ itkGetConstMacro( Lower, InputPixelType);
+ itkSetMacro( Upper, InputPixelType);
+ itkGetConstMacro( Upper, InputPixelType);
+ itkSetMacro( Inside, InternalPixelType);
+ itkGetConstMacro( Inside, InternalPixelType);
+ itkSetMacro( Outside, InternalPixelType);
+ itkGetConstMacro( Outside, InternalPixelType);
+ itkSetMacro( ErosionPaddingValue, OutputPixelType);
+ itkGetConstMacro( ErosionPaddingValue, OutputPixelType);
+ void SetRadius ( const SizeType& s){ m_Radius=s; this->Modified();}
+ SizeType GetRadius(void){return m_Radius;}
+ itkSetMacro( NumberOfNewLabels, unsigned int);
+ itkGetConstMacro( NumberOfNewLabels, unsigned int);
+ itkSetMacro( MinimumObjectSize, unsigned int);
+ itkGetConstMacro( MinimumObjectSize, unsigned int);
+ itkSetMacro( MinimumNumberOfIterations, unsigned int);
+ itkGetConstMacro( MinimumNumberOfIterations, unsigned int);
+
+ protected:
+
+ //----------------------------------------
+ // Constructor & Destructor
+ //----------------------------------------
+ DecomposeThroughErosionImageFilter();
+ ~DecomposeThroughErosionImageFilter() {};
+
+ //----------------------------------------
+ // Update
+ //----------------------------------------
+ // Generate Data
+ void GenerateData(void);
+ void AllocateOutput(){;}
+
+ // // Threaded Generate Data
+ // void BeforeThreadedGenerateData(void );
+ // void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId );
+ // void AfterThreadedGenerateData(void );
+ // // Override defaults
+ // virtual void GenerateInputRequestedRegion();
+ // virtual void GenerateOutputInformation (void);
+ // virtual void EnlargeOutputRequestedRegion(DataObject *data);
+ // void AllocateOutputs();
+ //----------------------------------------
+ // Data members
+ //----------------------------------------
+ bool m_Verbose;
+ bool m_FullyConnected;
+ InputPixelType m_Lower;
+ InputPixelType m_Upper;
+ OutputPixelType m_ErosionPaddingValue;
+ InputPixelType m_Inside;
+ InputPixelType m_Outside;
+ SizeType m_Radius;
+ unsigned int m_NumberOfNewLabels;
+ unsigned int m_MinimumObjectSize;
+ unsigned int m_MinimumNumberOfIterations;
+
+ };
+
+
+} // end namespace clitk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkDecomposeThroughErosionImageFilter.txx"
+#endif
+
+#endif // #define clitkDecomposeThroughErosionImageFilter_h
+
+
--- /dev/null
+#ifndef clitkDecomposeThroughErosionImageFilter_txx
+#define clitkDecomposeThroughErosionImageFilter_txx
+
+/* =================================================
+ * @file clitkDecomposeThroughErosionImageFilter.txx
+ * @author
+ * @date
+ *
+ * @brief
+ *
+ ===================================================*/
+
+
+namespace clitk
+{
+
+ //-------------------------------------------------------------------
+ // Update with the number of dimensions
+ //-------------------------------------------------------------------
+ template<class InputImageType, class OutputImageType>
+ DecomposeThroughErosionImageFilter<InputImageType, OutputImageType>::DecomposeThroughErosionImageFilter()
+ {
+ m_Verbose=false;
+ m_Lower =1;
+ m_Upper=1;
+ m_Inside=1;
+ m_Outside=0;
+ m_ErosionPaddingValue=static_cast<OutputPixelType>(-1);
+ for (unsigned int i=0; i<InputImageDimension; i++)
+ m_Radius[i]=1;
+ m_NumberOfNewLabels=1;
+ m_FullyConnected=true;
+ m_MinimumObjectSize=10;
+ m_MinimumNumberOfIterations=1;
+ }
+
+
+ //-------------------------------------------------------------------
+ // Update with the number of dimensions and the pixeltype
+ //-------------------------------------------------------------------
+ template<class InputImageType, class OutputImageType>
+ void
+ DecomposeThroughErosionImageFilter<InputImageType, OutputImageType>::GenerateData()
+ {
+
+ //---------------------------------
+ // Typedefs
+ //---------------------------------
+
+ // Internal type
+ typedef itk::Image<InternalPixelType, InputImageDimension> InternalImageType;
+
+ // Filters used
+ typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinaryThresholdImageFilter;
+ typedef itk::BinaryBallStructuringElement<InputPixelType,InputImageDimension > KernelType;
+ typedef itk::BinaryErodeImageFilter<InternalImageType, InternalImageType , KernelType> BinaryErodeImageFilterType;
+ typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> BinaryThresholdImageFilterType;
+ typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
+ typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
+ typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelImageFilterType;
+ typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType> SetBackgroundImageFilterType;
+
+ //---------------------------------
+ // Binarize input
+ //---------------------------------
+ typename InputBinaryThresholdImageFilter::Pointer inputBinarizer=InputBinaryThresholdImageFilter::New();
+ inputBinarizer->SetInput(this->GetInput());
+ inputBinarizer->SetLowerThreshold(m_Lower);
+ inputBinarizer->SetUpperThreshold(m_Upper);
+ inputBinarizer ->SetInsideValue(m_Inside);
+ inputBinarizer ->SetOutsideValue(m_Outside);
+ if(m_Verbose) std::cout<<"Binarizing the input..."<<std::endl;
+ inputBinarizer->Update();
+
+ //---------------------------------
+ // Label the input
+ //---------------------------------
+ typename ConnectFilterType::Pointer inputConnectFilter=ConnectFilterType::New();
+ inputConnectFilter->SetInput(inputBinarizer->GetOutput());
+ inputConnectFilter->SetBackgroundValue(0);
+ inputConnectFilter->SetFullyConnected(m_FullyConnected);
+ if(m_Verbose) std::cout<<"Labelling the connected components..."<<std::endl;
+ // inputConnectFilter->Update();
+
+ //---------------------------------
+ // Count the initial labels
+ //---------------------------------
+ typename StatisticsImageFilterType::Pointer inputStatisticsImageFilter=StatisticsImageFilterType::New();
+ inputStatisticsImageFilter->SetInput(inputConnectFilter->GetOutput());
+ if(m_Verbose) std::cout<<"Counting the initial labels..."<<std::endl;
+ inputStatisticsImageFilter->Update();
+ unsigned int initialNumberOfLabels= inputStatisticsImageFilter->GetMaximum();
+ if(m_Verbose) std::cout<<"The input contained "<<initialNumberOfLabels<<" disctictive label(s)..."<<std::endl;
+ if(m_Verbose) std::cout<<"Performing erosions till at least "<<initialNumberOfLabels+m_NumberOfNewLabels<<" distinctive labels are counted..."<<std::endl;
+
+ //---------------------------------
+ // Structuring element
+ //---------------------------------
+ KernelType structuringElement;
+ structuringElement.SetRadius(m_Radius);
+ structuringElement.CreateStructuringElement();
+
+ //---------------------------------
+ // Repeat while not decomposed
+ //---------------------------------
+ typename InternalImageType::Pointer current=inputBinarizer->GetOutput();
+ typename InternalImageType::Pointer output=inputBinarizer->GetOutput();
+ unsigned int iteration=0;
+ unsigned int max =initialNumberOfLabels;
+
+ while ( (iteration < m_MinimumNumberOfIterations) || ( (max< initialNumberOfLabels + m_NumberOfNewLabels ) && (iteration<100 ) ) )
+ {
+
+
+ if(m_Verbose) std::cout<<"Eroding image (iteration "<<iteration<<")..."<<std::endl;
+
+ //---------------------------------
+ // Erode
+ //---------------------------------
+ typename BinaryErodeImageFilterType::Pointer erosionFilter=BinaryErodeImageFilterType::New();
+ erosionFilter->SetInput (current);
+ erosionFilter->SetForegroundValue (1);
+ erosionFilter->SetBackgroundValue (-1);
+ erosionFilter->SetBoundaryToForeground(false);
+ erosionFilter->SetKernel(structuringElement);
+ erosionFilter->Update();
+ current=erosionFilter->GetOutput();
+
+ //---------------------------------
+ // Binarize (remove -1)
+ //---------------------------------
+ typename BinaryThresholdImageFilterType::Pointer binarizer=BinaryThresholdImageFilterType::New();
+ binarizer->SetInput(erosionFilter->GetOutput());
+ binarizer->SetLowerThreshold(1);
+ binarizer->SetUpperThreshold(1);
+ binarizer ->SetInsideValue(1);
+ binarizer ->SetOutsideValue(0);
+ if(m_Verbose) std::cout<<"Binarizing the eroded image..."<<std::endl;
+ //binarizer->Update();
+
+
+ //---------------------------------
+ // ReLabel the connected components
+ //---------------------------------
+ typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
+ connectFilter->SetInput(binarizer->GetOutput());
+ connectFilter->SetBackgroundValue(0);
+ connectFilter->SetFullyConnected(m_FullyConnected);
+ if(m_Verbose) std::cout<<"Labelling the connected components..."<<std::endl;
+ //connectFilter->Update();
+
+ //---------------------------------
+ // Sort
+ //---------------------------------
+ typename RelabelImageFilterType::Pointer relabelFilter=RelabelImageFilterType::New();
+ relabelFilter->SetInput(connectFilter->GetOutput());
+ relabelFilter->SetMinimumObjectSize(m_MinimumObjectSize);
+ //relabelFilter->Update();
+
+
+ //---------------------------------
+ // Count the labels
+ //---------------------------------
+ typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
+ statisticsImageFilter->SetInput(relabelFilter->GetOutput());
+ statisticsImageFilter->Update();
+ max= statisticsImageFilter->GetMaximum();
+ if(m_Verbose) std::cout<<"Counted "<<max<<" label (s) larger then "<<m_MinimumObjectSize<<" voxels..."<<std::endl;
+ output=statisticsImageFilter->GetOutput();
+
+ // Next iteration
+ iteration++;
+ }
+
+
+ //---------------------------------
+ // Binarize current (remove -1)
+ //---------------------------------
+ typename BinaryThresholdImageFilterType::Pointer binarizer=BinaryThresholdImageFilterType::New();
+ binarizer->SetInput(current);
+ binarizer->SetLowerThreshold(1);
+ binarizer->SetUpperThreshold(1);
+ binarizer ->SetInsideValue(1);
+ binarizer ->SetOutsideValue(0);
+ if(m_Verbose) std::cout<<"Binarizing the eroded image..."<<std::endl;
+ //binarizer->Update();
+
+ //---------------------------------
+ // ReLabel the connected components
+ //---------------------------------
+ typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
+ connectFilter->SetInput(binarizer->GetOutput());
+ connectFilter->SetBackgroundValue(0);
+ connectFilter->SetFullyConnected(m_FullyConnected);
+ if(m_Verbose) std::cout<<"Labelling the connected components..."<<std::endl;
+ connectFilter->Update();
+
+ //---------------------------------
+ // Sort
+ //---------------------------------
+ typename RelabelImageFilterType::Pointer relabelFilter=RelabelImageFilterType::New();
+ relabelFilter->SetInput(connectFilter->GetOutput());
+ //relabelFilter->SetMinimumObjectSize(m_MinimumObjectSize); // Preserve all intensities
+ //relabelFilter->Update();
+
+ //---------------------------------
+ // Set -1 to padding value
+ //---------------------------------
+ typename SetBackgroundImageFilterType::Pointer setBackgroundFilter =SetBackgroundImageFilterType::New();
+ setBackgroundFilter->SetInput(relabelFilter->GetOutput());
+ setBackgroundFilter->SetInput2(current);
+ setBackgroundFilter->SetMaskValue(-1);
+ setBackgroundFilter->SetOutsideValue(m_ErosionPaddingValue);
+ if(m_Verbose) std::cout<<"Setting the eroded region to "<<m_ErosionPaddingValue<<"..."<<std::endl;
+
+ //---------------------------------
+ // Cast
+ //---------------------------------
+ typedef itk::CastImageFilter<InternalImageType, OutputImageType> CastImageFilterType;
+ typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
+ caster->SetInput(setBackgroundFilter->GetOutput());
+ caster->Update();
+
+ //---------------------------------
+ // SetOutput
+ //---------------------------------
+ this->SetNthOutput(0, caster->GetOutput());
+ }
+
+
+}//end clitk
+
+#endif //#define clitkDecomposeThroughErosionImageFilter_txx
--- /dev/null
+#ifndef __clitkExplosionControlledThresholdConnectedImageFilter_h
+#define __clitkExplosionControlledThresholdConnectedImageFilter_h
+
+#include "itkImage.h"
+#include "itkImageToImageFilter.h"
+
+namespace clitk {
+
+/** \class ExplosionControlledThresholdConnectedImageFilter
+ * \brief Label pixels that are connected to a seed and lie within a neighborhood
+ *
+ * ExplosionControlledThresholdConnectedImageFilter labels pixels with ReplaceValue that
+ * are connected to an initial Seed AND whose neighbors all lie within a
+ * Lower and Upper threshold range.
+ *
+ * \ingroup RegionGrowingSegmentation
+ */
+template <class TInputImage, class TOutputImage>
+class ITK_EXPORT ExplosionControlledThresholdConnectedImageFilter:
+ public itk::ImageToImageFilter<TInputImage,TOutputImage>
+{
+public:
+ /** Standard class typedefs. */
+ typedef ExplosionControlledThresholdConnectedImageFilter Self;
+ typedef itk::ImageToImageFilter<TInputImage,TOutputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(ExplosionControlledThresholdConnectedImageFilter,
+ ImageToImageFilter);
+
+ typedef TInputImage InputImageType;
+ typedef typename InputImageType::Pointer InputImagePointer;
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename InputImageType::PixelType InputImagePixelType;
+ typedef typename InputImageType::IndexType IndexType;
+ typedef typename InputImageType::SizeType InputImageSizeType;
+
+ typedef TOutputImage OutputImageType;
+ typedef typename OutputImageType::Pointer OutputImagePointer;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+ typedef typename OutputImageType::PixelType OutputImagePixelType;
+
+ void PrintSelf ( std::ostream& os, itk::Indent indent ) const;
+
+ /** Clear the seeds */
+ void ClearSeeds();
+
+ /** Set seed point. */
+ void SetSeed(const IndexType & seed);
+
+ /** Add a seed point */
+ void AddSeed ( const IndexType & seed );
+
+ // Set/Get the lower threshold. The default is 0.
+ itkSetMacro(Lower, InputImagePixelType);
+ itkGetConstMacro(Lower, InputImagePixelType);
+
+ // Set/Get the upper threshold. The default is the largest possible value for the InputPixelType.
+ itkSetMacro(Upper, InputImagePixelType);
+ itkGetConstMacro(Upper, InputImagePixelType);
+
+ /** Set/Get value to replace thresholded pixels. Pixels that lie *
+ * within Lower and Upper (inclusive) will be replaced with this
+ * value. The default is 1. */
+ itkSetMacro(ReplaceValue, OutputImagePixelType);
+ itkGetConstMacro(ReplaceValue, OutputImagePixelType);
+
+ // /** Set the radius of the neighborhood used for a mask. */
+ // itkSetMacro(Radius, InputImageSizeType);
+
+ // /** Get the radius of the neighborhood used to compute the median */
+ // itkGetConstReferenceMacro(Radius, InputImageSizeType);
+
+ /** ImageDimension constants */
+ itkStaticConstMacro(InputImageDimension, unsigned int,
+ TInputImage::ImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int,
+ TOutputImage::ImageDimension);
+
+#ifdef ITK_USE_CONCEPT_CHECKING
+ /** Begin concept checking */
+ itkConceptMacro(InputEqualityComparableCheck,
+ (itk::Concept::EqualityComparable<InputImagePixelType>));
+ itkConceptMacro(OutputEqualityComparableCheck,
+ (itk::Concept::EqualityComparable<OutputImagePixelType>));
+ itkConceptMacro(SameDimensionCheck,
+ (itk::Concept::SameDimension<InputImageDimension, OutputImageDimension>));
+ itkConceptMacro(InputOStreamWritableCheck,
+ (itk::Concept::OStreamWritable<InputImagePixelType>));
+ itkConceptMacro(OutputOStreamWritableCheck,
+ (itk::Concept::OStreamWritable<OutputImagePixelType>));
+ /** End concept checking */
+#endif
+
+ // JV
+ itkBooleanMacro(Verbose);
+ itkSetMacro( Verbose, bool);
+ itkGetConstReferenceMacro( Verbose, bool);
+
+ itkBooleanMacro(FullyConnected);
+ itkSetMacro( FullyConnected, bool);
+ itkGetConstReferenceMacro( FullyConnected, bool);
+
+ itkSetMacro( FinalLower, InputImagePixelType);
+ itkGetConstMacro( FinalLower, InputImagePixelType);
+
+ itkSetMacro( FinalUpper, InputImagePixelType);
+ itkGetConstMacro( FinalUpper, InputImagePixelType);
+
+ itkBooleanMacro(AdaptLowerBorder);
+ itkSetMacro( AdaptLowerBorder, bool);
+ itkGetConstReferenceMacro( AdaptLowerBorder, bool);
+
+ itkBooleanMacro(AdaptUpperBorder);
+ itkSetMacro( AdaptUpperBorder, bool);
+ itkGetConstReferenceMacro( AdaptUpperBorder, bool);
+
+ itkSetMacro( Multiplier, double);
+ itkGetConstMacro( Multiplier, double);
+
+ itkSetMacro( MaximumUpperThreshold, InputImagePixelType);
+ itkGetConstMacro( MaximumUpperThreshold, InputImagePixelType);
+
+ itkSetMacro( MinimumLowerThreshold, InputImagePixelType);
+ itkGetConstMacro( MinimumLowerThreshold, InputImagePixelType);
+
+ itkSetMacro(ThresholdStepSize, InputImagePixelType);
+ itkGetConstMacro( ThresholdStepSize, InputImagePixelType);
+
+ itkSetMacro( MinimumThresholdStepSize, InputImagePixelType);
+ itkGetConstMacro( MinimumThresholdStepSize, InputImagePixelType);
+
+ itkSetMacro( MinimumSize, unsigned int);
+ itkGetConstMacro( MinimumSize,unsigned int);
+
+protected:
+ ExplosionControlledThresholdConnectedImageFilter();
+ ~ExplosionControlledThresholdConnectedImageFilter(){};
+ std::vector<IndexType> m_Seeds;
+ InputImagePixelType m_Lower;
+ InputImagePixelType m_Upper;
+ OutputImagePixelType m_ReplaceValue;
+
+ // JV
+ bool m_Verbose;
+ bool m_FullyConnected;
+ bool m_AdaptLowerBorder;
+ bool m_AdaptUpperBorder;
+ InputImagePixelType m_FinalLower;
+ InputImagePixelType m_FinalUpper;
+ double m_Multiplier;
+ InputImagePixelType m_MaximumUpperThreshold;
+ InputImagePixelType m_MinimumLowerThreshold;
+ InputImagePixelType m_MinimumThresholdStepSize;
+ InputImagePixelType m_ThresholdStepSize;
+ unsigned int m_MinimumSize;
+
+ // Override since the filter needs all the data for the algorithm
+ void GenerateInputRequestedRegion();
+
+ // Override since the filter produces the entire dataset
+ void EnlargeOutputRequestedRegion(itk::DataObject *output);
+ void GenerateData();
+
+private:
+ ExplosionControlledThresholdConnectedImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+};
+
+} // end namespace clitk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkExplosionControlledThresholdConnectedImageFilter.txx"
+#endif
+
+#endif
--- /dev/null
+#ifndef __clitkExplosionControlledThresholdConnectedImageFilter_txx
+#define __clitkExplosionControlledThresholdConnectedImageFilter_txx
+
+#include "clitkExplosionControlledThresholdConnectedImageFilter.h"
+#include "itkBinaryThresholdImageFunction.h"
+#include "itkFloodFilledImageFunctionConditionalIterator.h"
+#include "itkShapedFloodFilledImageFunctionConditionalIterator.h"
+#include "itkProgressReporter.h"
+
+namespace clitk
+{
+
+ //--------------------------------------------------------------------
+ /**
+ * Constructor
+ */
+ template <class TInputImage, class TOutputImage>
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
+ ::ExplosionControlledThresholdConnectedImageFilter()
+ {
+ m_Lower = itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
+ m_Upper = itk::NumericTraits<InputImagePixelType>::max();
+ m_FinalLower = itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
+ m_FinalUpper = itk::NumericTraits<InputImagePixelType>::max();
+ m_AdaptLowerBorder=false;
+ m_AdaptUpperBorder=false;
+ m_MaximumUpperThreshold=itk::NumericTraits<InputImagePixelType>::max();
+ m_MinimumLowerThreshold=itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
+ m_Multiplier=2.0;
+ m_ReplaceValue = itk::NumericTraits<OutputImagePixelType>::One;
+ m_MinimumThresholdStepSize=itk::NumericTraits<InputImagePixelType>::One;
+ m_ThresholdStepSize=64;
+ m_Verbose=false;
+ m_FullyConnected=false;
+ m_MinimumSize=5000;
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
+ ::ClearSeeds()
+ {
+ if( this->m_Seeds.size() > 0 )
+ {
+ this->m_Seeds.clear();
+ this->Modified();
+ }
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
+ ::SetSeed(const IndexType & seed)
+ {
+ this->ClearSeeds();
+ this->AddSeed ( seed );
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
+ ::AddSeed ( const IndexType & seed )
+ {
+ this->m_Seeds.push_back ( seed );
+ this->Modified();
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ /**
+ * Standard PrintSelf method.
+ */
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
+ ::PrintSelf(std::ostream& os, itk::Indent indent) const
+ {
+ this->Superclass::PrintSelf(os, indent);
+ os << indent << "Upper: "
+ << static_cast<typename itk::NumericTraits<InputImagePixelType>::PrintType>(m_Upper)
+ << std::endl;
+ os << indent << "Lower: "
+ << static_cast<typename itk::NumericTraits<InputImagePixelType>::PrintType>(m_Lower)
+ << std::endl;
+ os << indent << "ReplaceValue: "
+ << static_cast<double>(m_Multiplier)
+ << std::endl;
+ os << indent << "ReplaceValue: "
+ << static_cast<typename itk::NumericTraits<OutputImagePixelType>::PrintType>(m_ReplaceValue)
+ << std::endl;
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
+ ::GenerateInputRequestedRegion()
+ {
+ Superclass::GenerateInputRequestedRegion();
+ if ( this->GetInput() )
+ {
+ InputImagePointer image =
+ const_cast< InputImageType * >( this->GetInput() );
+ image->SetRequestedRegionToLargestPossibleRegion();
+ }
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
+ ::EnlargeOutputRequestedRegion(itk::DataObject *output)
+ {
+ Superclass::EnlargeOutputRequestedRegion(output);
+ output->SetRequestedRegionToLargestPossibleRegion();
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template <class TInputImage, class TOutputImage>
+ void
+ ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
+ ::GenerateData()
+ {
+ typename Superclass::InputImageConstPointer inputImage = this->GetInput();
+ typename Superclass::OutputImagePointer outputImage = this->GetOutput();
+
+ // Zero the output
+ outputImage->SetBufferedRegion( outputImage->GetRequestedRegion() );
+ outputImage->Allocate();
+ outputImage->FillBuffer ( itk::NumericTraits<OutputImagePixelType>::Zero );
+
+
+ //--------------------------------------------------
+ // Get initial region size
+ //--------------------------------------------------
+ typedef itk::BinaryThresholdImageFunction<InputImageType> FunctionType;
+ //typedef itk::ShapedFloodFilledImageFunctionConditionalIterator<OutputImageType, FunctionType> IteratorType;
+ typedef itk::FloodFilledImageFunctionConditionalIterator<OutputImageType, FunctionType> IteratorType;
+ typename FunctionType::Pointer function = FunctionType::New();
+ function->SetInputImage ( inputImage );
+ function->ThresholdBetween ( m_Lower, m_Upper );
+ IteratorType it = IteratorType ( outputImage, function, m_Seeds );
+ //it.SetFullyConnected(m_FullyConnected);
+ unsigned int initialSize=0;
+
+ it.GoToBegin();
+ while( !it.IsAtEnd())
+ {
+ ++initialSize;
+ ++it;
+ }
+ if (m_Verbose) std::cout<<"Initial region size using thresholds ["<<m_Lower<<", "<<m_Upper<<"], is "<<initialSize<<"..."<<std::endl;
+
+
+ //--------------------------------------------------
+ // Decrease lower threshold
+ //--------------------------------------------------
+ m_FinalLower=m_Lower;
+ if (m_AdaptLowerBorder)
+ {
+ // Init
+ InputImagePixelType currentLower=m_Lower;
+ unsigned int currentSize=initialSize;
+ unsigned int previousSize=initialSize;
+
+ // Lower the threshold till explosion
+ while ( (previousSize<m_MinimumSize) || ( (currentLower> m_MinimumLowerThreshold) && ( (unsigned int) currentSize <= (unsigned int )m_Multiplier* previousSize) ) )
+ {
+ // Update threshold
+ currentLower-=m_ThresholdStepSize;
+ function->ThresholdBetween(currentLower, m_Upper);
+
+ // Get current region size
+ previousSize=currentSize;
+ currentSize=0;
+ IteratorType it = IteratorType ( outputImage, function, m_Seeds );
+ it.GoToBegin();
+ while( !it.IsAtEnd())
+ {
+ ++currentSize;
+ ++it;
+ }
+ if (m_Verbose) std::cout<<"Decreasing lower threshold to "<<currentLower<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
+ }
+
+
+ // Explosion occured, increase lower theshold
+ if ( (double)currentSize > m_Multiplier*(double) previousSize)
+ {
+ // Raise lower threshold
+ if (m_Verbose) std::cout<<"Explosion detected, adapting lower threshold ..."<<std::endl;
+ currentLower+=m_ThresholdStepSize;
+ InputImagePixelType currentStepSize=m_ThresholdStepSize;
+
+ while (currentStepSize>m_MinimumThresholdStepSize)
+ {
+ currentStepSize/=2;
+ currentLower-=currentStepSize;
+
+ // Get current region size
+ currentSize=0;
+ IteratorType it = IteratorType ( outputImage, function, m_Seeds );
+ it.GoToBegin();
+ while( !it.IsAtEnd())
+ {
+ ++currentSize;
+ ++it;
+ }
+
+ if (m_Verbose) std::cout<<"Adapting lower threshold to "<<currentLower<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
+
+ // Explosion: go back
+ if ((double)currentSize > m_Multiplier* (double) previousSize)
+ currentLower+=currentStepSize;
+
+ // No explosion: adapt previous
+ else previousSize=currentSize;
+ }
+
+ // Converged
+ m_FinalLower=currentLower;
+ if (m_Verbose) std::cout<<"Final lower threshold (precision="<<m_MinimumThresholdStepSize<<") is set to "<<m_FinalLower<<"..."<<std::endl;
+ }
+ else {
+ if (m_Verbose) std::cout<<"No explosion before minimum lower threshold reached!"<<std::endl;
+ }
+
+ } // end update lower
+
+
+ //--------------------------------------------------
+ // Increase upper threshold
+ //--------------------------------------------------
+ m_FinalUpper=m_Upper;
+ if (m_AdaptUpperBorder)
+ {
+ // Init
+ InputImagePixelType currentUpper=m_Upper;
+ unsigned int currentSize=initialSize;
+ unsigned int previousSize=initialSize;
+
+ // Upper the threshold till explosion
+ while ((previousSize<m_MinimumSize) || ( (currentUpper< m_MaximumUpperThreshold) && ((double)currentSize <= m_Multiplier* (double) previousSize) ) )
+ {
+ // Update threshold
+ currentUpper+=m_ThresholdStepSize;
+ function->ThresholdBetween(m_Lower,currentUpper);
+
+ // Get current region size
+ previousSize=currentSize;
+ currentSize=0;
+ IteratorType it = IteratorType ( outputImage, function, m_Seeds );
+ it.GoToBegin();
+ while( !it.IsAtEnd())
+ {
+ ++currentSize;
+ ++it;
+ }
+ if (m_Verbose) std::cout<<"Increasing upper threshold to "<<currentUpper<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
+
+ }
+
+ // Explosion occured, decrease upper theshold
+ if ((double)currentSize > m_Multiplier* (double) previousSize)
+ {
+ // Lower upper threshold
+ if (m_Verbose) std::cout<<"Explosion detected, adapting upper threshold ..."<<std::endl;
+ currentUpper-=m_ThresholdStepSize;
+ InputImagePixelType currentStepSize=m_ThresholdStepSize;
+
+ while (currentStepSize>m_MinimumThresholdStepSize)
+ {
+ currentStepSize/=2;
+ currentUpper+=currentStepSize;
+
+ // Get current region size
+ currentSize=0;
+ function->ThresholdBetween(m_Lower,currentUpper);
+ IteratorType it = IteratorType ( outputImage, function, m_Seeds );
+ it.GoToBegin();
+ while( !it.IsAtEnd())
+ {
+ ++currentSize;
+ ++it;
+ }
+ if (m_Verbose) std::cout<<"Adapting upper threshold to "<<currentUpper<<", size is "<<currentSize <<"(previously "<<previousSize<<") ..."<<std::endl;
+
+ // Explosion: go back
+ if ((double)currentSize > m_Multiplier* (double) previousSize)
+ currentUpper-=currentStepSize;
+
+ // No explosion: adapt previous
+ else previousSize=currentSize;
+ }
+
+ // Converged
+ m_FinalUpper=currentUpper;
+ if (m_Verbose) std::cout<<"Final upper threshold (precision="<<m_MinimumThresholdStepSize<<") is set to "<<m_FinalUpper<<"..."<<std::endl;
+
+ }
+ else { if (m_Verbose) std::cout<<"No explosion before maximum upper threshold reached!"<<std::endl; }
+
+ } // end update upper
+
+
+
+ //--------------------------------------------------
+ // Update the output with final thresholds
+ //--------------------------------------------------
+ function->ThresholdBetween(m_FinalLower, m_FinalUpper);
+ IteratorType it2 = IteratorType ( outputImage, function, m_Seeds );
+ it2.GoToBegin();
+ while( !it2.IsAtEnd())
+ {
+ it2.Set(m_ReplaceValue);
+ ++it2;
+ }
+
+ }
+ //--------------------------------------------------------------------
+
+
+} // end namespace clitk
+
+#endif
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff