From 23e1206b38d3b26f254fd01ebb5e0f0b5bb193fa Mon Sep 17 00:00:00 2001 From: dsarrut Date: Wed, 30 Jun 2010 06:00:16 +0000 Subject: [PATCH] segmentation utils (jef) --- itk/clitkConditionalBinaryDilateImageFilter.h | 115 ++++ ...litkConditionalBinaryDilateImageFilter.txx | 509 ++++++++++++++++++ itk/clitkDecomposeAndReconstructImageFilter.h | 147 +++++ ...litkDecomposeAndReconstructImageFilter.txx | 81 +++ itk/clitkDecomposeThroughErosionImageFilter.h | 142 +++++ ...litkDecomposeThroughErosionImageFilter.txx | 233 ++++++++ ...nControlledThresholdConnectedImageFilter.h | 182 +++++++ ...ontrolledThresholdConnectedImageFilter.txx | 341 ++++++++++++ 8 files changed, 1750 insertions(+) create mode 100644 itk/clitkConditionalBinaryDilateImageFilter.h create mode 100644 itk/clitkConditionalBinaryDilateImageFilter.txx create mode 100644 itk/clitkDecomposeAndReconstructImageFilter.h create mode 100644 itk/clitkDecomposeAndReconstructImageFilter.txx create mode 100644 itk/clitkDecomposeThroughErosionImageFilter.h create mode 100644 itk/clitkDecomposeThroughErosionImageFilter.txx create mode 100644 itk/clitkExplosionControlledThresholdConnectedImageFilter.h create mode 100644 itk/clitkExplosionControlledThresholdConnectedImageFilter.txx diff --git a/itk/clitkConditionalBinaryDilateImageFilter.h b/itk/clitkConditionalBinaryDilateImageFilter.h new file mode 100644 index 0000000..ff7926c --- /dev/null +++ b/itk/clitkConditionalBinaryDilateImageFilter.h @@ -0,0 +1,115 @@ +#ifndef clitkConditionalBinaryDilateImageFilter_h +#define clitkConditionalBinaryDilateImageFilter_h + +/* ================================================= + * @file clitkConditionalBinaryDilateImageFilter.h + * @author + * @date + * + * @brief + * + ===================================================*/ + +// clitk include +#include "clitkIO.h" +#include "clitkCommon.h" + +// itk include +#include +#include +#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 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 Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer 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::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 diff --git a/itk/clitkConditionalBinaryDilateImageFilter.txx b/itk/clitkConditionalBinaryDilateImageFilter.txx new file mode 100644 index 0000000..d82468c --- /dev/null +++ b/itk/clitkConditionalBinaryDilateImageFilter.txx @@ -0,0 +1,509 @@ +#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 + ConditionalBinaryDilateImageFilter + ::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 outIt( output, outputRegion ); + itk::ImageRegionConstIterator 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 ( backgroundValue ) ); } + // keep all of the original background values intact + else + { outIt.Set( static_cast( 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 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 iRegIt( input, requiredInputRegion ); + // iterator on tmp image + itk::ImageRegionIterator 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 + tmpRegIndexIt( tmpImage, tmpRequestedRegion ); + + itk::ConstNeighborhoodIterator + oNeighbIt( radius, tmpImage, tmpRequestedRegion ); + + // Define boundaries conditions + itk::ConstantBoundaryCondition cbc; + cbc.SetConstant( backgroundTag ); + oNeighbIt.OverrideBoundaryCondition(&cbc); + + unsigned int neighborhoodSize = oNeighbIt.Size(); + unsigned int centerPixelCode = neighborhoodSize / 2; + + std::queue 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 + nit( radius, tmpImage, tmpRequestedRegion ); + nit.OverrideBoundaryCondition(&cbc); + nit.GoToBegin(); + + itk::ConstNeighborhoodIterator + 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 ( foregroundValue ) ); } + if (input->GetPixel(idx)==backgroundValue) + output->SetPixel( idx, static_cast ( 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 ( foregroundValue ) ); } + if (input->GetPixel(idx)==backgroundValue) + output->SetPixel( idx, static_cast ( 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 obc; + obc.SetConstant( static_cast ( backgroundValue ) ); + + itk::NeighborhoodIterator + 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 + 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 ( 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 ( foregroundValue ) ); + break; + } + } + + ++ouRegIndexIt; + progress.CompletedPixel(); + } + } + } + + + /** + * Standard "PrintSelf" method + */ + template + void + ConditionalBinaryDilateImageFilter + ::PrintSelf( std::ostream& os, itk::Indent indent) const + { + Superclass::PrintSelf( os, indent ); + os << indent << "Dilate Value: " << static_cast::PrintType>( this->GetForegroundValue() ) << std::endl; + } + +} // end namespace itk + +#endif diff --git a/itk/clitkDecomposeAndReconstructImageFilter.h b/itk/clitkDecomposeAndReconstructImageFilter.h new file mode 100644 index 0000000..43cbecb --- /dev/null +++ b/itk/clitkDecomposeAndReconstructImageFilter.h @@ -0,0 +1,147 @@ +#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 ITK_EXPORT DecomposeAndReconstructImageFilter : + public itk::ImageToImageFilter + { + public: + //---------------------------------------- + // ITK + //---------------------------------------- + typedef DecomposeAndReconstructImageFilter Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer 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 + + diff --git a/itk/clitkDecomposeAndReconstructImageFilter.txx b/itk/clitkDecomposeAndReconstructImageFilter.txx new file mode 100644 index 0000000..757a2ad --- /dev/null +++ b/itk/clitkDecomposeAndReconstructImageFilter.txx @@ -0,0 +1,81 @@ +#ifndef clitkDecomposeAndReconstructImageFilter_txx +#define clitkDecomposeAndReconstructImageFilter_txx + +/* ================================================= + * @file clitkDecomposeAndReconstructImageFilter.txx + * @author + * @date + * + * @brief + * + ===================================================*/ + + +namespace clitk +{ + + //------------------------------------------------------------------- + // Update with the number of dimensions + //------------------------------------------------------------------- + template + DecomposeAndReconstructImageFilter::DecomposeAndReconstructImageFilter() + { + m_Verbose=false; + for (unsigned int i=0; i + void + DecomposeAndReconstructImageFilter::GenerateData() + { + + + // Internal type + typedef itk::Image InternalImageType; + + // Filters used + typedef clitk::DecomposeThroughErosionImageFilter DecomposeThroughErosionImageFilterType; + typedef clitk::ReconstructThroughDilationImageFilter 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 diff --git a/itk/clitkDecomposeThroughErosionImageFilter.h b/itk/clitkDecomposeThroughErosionImageFilter.h new file mode 100644 index 0000000..3c9ceac --- /dev/null +++ b/itk/clitkDecomposeThroughErosionImageFilter.h @@ -0,0 +1,142 @@ +#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 ITK_EXPORT DecomposeThroughErosionImageFilter : + public itk::ImageToImageFilter + { + public: + //---------------------------------------- + // ITK + //---------------------------------------- + typedef DecomposeThroughErosionImageFilter Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer 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 + + diff --git a/itk/clitkDecomposeThroughErosionImageFilter.txx b/itk/clitkDecomposeThroughErosionImageFilter.txx new file mode 100644 index 0000000..f5f33b0 --- /dev/null +++ b/itk/clitkDecomposeThroughErosionImageFilter.txx @@ -0,0 +1,233 @@ +#ifndef clitkDecomposeThroughErosionImageFilter_txx +#define clitkDecomposeThroughErosionImageFilter_txx + +/* ================================================= + * @file clitkDecomposeThroughErosionImageFilter.txx + * @author + * @date + * + * @brief + * + ===================================================*/ + + +namespace clitk +{ + + //------------------------------------------------------------------- + // Update with the number of dimensions + //------------------------------------------------------------------- + template + DecomposeThroughErosionImageFilter::DecomposeThroughErosionImageFilter() + { + m_Verbose=false; + m_Lower =1; + m_Upper=1; + m_Inside=1; + m_Outside=0; + m_ErosionPaddingValue=static_cast(-1); + for (unsigned int i=0; i + void + DecomposeThroughErosionImageFilter::GenerateData() + { + + //--------------------------------- + // Typedefs + //--------------------------------- + + // Internal type + typedef itk::Image InternalImageType; + + // Filters used + typedef itk::BinaryThresholdImageFilter InputBinaryThresholdImageFilter; + typedef itk::BinaryBallStructuringElement KernelType; + typedef itk::BinaryErodeImageFilter BinaryErodeImageFilterType; + typedef itk::BinaryThresholdImageFilter BinaryThresholdImageFilterType; + typedef itk::StatisticsImageFilter StatisticsImageFilterType; + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typedef itk::RelabelComponentImageFilter RelabelImageFilterType; + typedef clitk::SetBackgroundImageFilter 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..."<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..."<Update(); + + //--------------------------------- + // Count the initial labels + //--------------------------------- + typename StatisticsImageFilterType::Pointer inputStatisticsImageFilter=StatisticsImageFilterType::New(); + inputStatisticsImageFilter->SetInput(inputConnectFilter->GetOutput()); + if(m_Verbose) std::cout<<"Counting the initial labels..."<Update(); + unsigned int initialNumberOfLabels= inputStatisticsImageFilter->GetMaximum(); + if(m_Verbose) std::cout<<"The input contained "<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 "<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..."<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..."<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 "<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..."<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..."<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 "< 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 diff --git a/itk/clitkExplosionControlledThresholdConnectedImageFilter.h b/itk/clitkExplosionControlledThresholdConnectedImageFilter.h new file mode 100644 index 0000000..e325f37 --- /dev/null +++ b/itk/clitkExplosionControlledThresholdConnectedImageFilter.h @@ -0,0 +1,182 @@ +#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 ITK_EXPORT ExplosionControlledThresholdConnectedImageFilter: + public itk::ImageToImageFilter +{ +public: + /** Standard class typedefs. */ + typedef ExplosionControlledThresholdConnectedImageFilter Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer 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)); + itkConceptMacro(OutputEqualityComparableCheck, + (itk::Concept::EqualityComparable)); + itkConceptMacro(SameDimensionCheck, + (itk::Concept::SameDimension)); + itkConceptMacro(InputOStreamWritableCheck, + (itk::Concept::OStreamWritable)); + itkConceptMacro(OutputOStreamWritableCheck, + (itk::Concept::OStreamWritable)); + /** 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 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 diff --git a/itk/clitkExplosionControlledThresholdConnectedImageFilter.txx b/itk/clitkExplosionControlledThresholdConnectedImageFilter.txx new file mode 100644 index 0000000..0674414 --- /dev/null +++ b/itk/clitkExplosionControlledThresholdConnectedImageFilter.txx @@ -0,0 +1,341 @@ +#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 + ExplosionControlledThresholdConnectedImageFilter + ::ExplosionControlledThresholdConnectedImageFilter() + { + m_Lower = itk::NumericTraits::NonpositiveMin(); + m_Upper = itk::NumericTraits::max(); + m_FinalLower = itk::NumericTraits::NonpositiveMin(); + m_FinalUpper = itk::NumericTraits::max(); + m_AdaptLowerBorder=false; + m_AdaptUpperBorder=false; + m_MaximumUpperThreshold=itk::NumericTraits::max(); + m_MinimumLowerThreshold=itk::NumericTraits::NonpositiveMin(); + m_Multiplier=2.0; + m_ReplaceValue = itk::NumericTraits::One; + m_MinimumThresholdStepSize=itk::NumericTraits::One; + m_ThresholdStepSize=64; + m_Verbose=false; + m_FullyConnected=false; + m_MinimumSize=5000; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExplosionControlledThresholdConnectedImageFilter + ::ClearSeeds() + { + if( this->m_Seeds.size() > 0 ) + { + this->m_Seeds.clear(); + this->Modified(); + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExplosionControlledThresholdConnectedImageFilter + ::SetSeed(const IndexType & seed) + { + this->ClearSeeds(); + this->AddSeed ( seed ); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExplosionControlledThresholdConnectedImageFilter + ::AddSeed ( const IndexType & seed ) + { + this->m_Seeds.push_back ( seed ); + this->Modified(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + /** + * Standard PrintSelf method. + */ + template + void + ExplosionControlledThresholdConnectedImageFilter + ::PrintSelf(std::ostream& os, itk::Indent indent) const + { + this->Superclass::PrintSelf(os, indent); + os << indent << "Upper: " + << static_cast::PrintType>(m_Upper) + << std::endl; + os << indent << "Lower: " + << static_cast::PrintType>(m_Lower) + << std::endl; + os << indent << "ReplaceValue: " + << static_cast(m_Multiplier) + << std::endl; + os << indent << "ReplaceValue: " + << static_cast::PrintType>(m_ReplaceValue) + << std::endl; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExplosionControlledThresholdConnectedImageFilter + ::GenerateInputRequestedRegion() + { + Superclass::GenerateInputRequestedRegion(); + if ( this->GetInput() ) + { + InputImagePointer image = + const_cast< InputImageType * >( this->GetInput() ); + image->SetRequestedRegionToLargestPossibleRegion(); + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExplosionControlledThresholdConnectedImageFilter + ::EnlargeOutputRequestedRegion(itk::DataObject *output) + { + Superclass::EnlargeOutputRequestedRegion(output); + output->SetRequestedRegionToLargestPossibleRegion(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ExplosionControlledThresholdConnectedImageFilter + ::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::Zero ); + + + //-------------------------------------------------- + // Get initial region size + //-------------------------------------------------- + typedef itk::BinaryThresholdImageFunction FunctionType; + //typedef itk::ShapedFloodFilledImageFunctionConditionalIterator IteratorType; + typedef itk::FloodFilledImageFunctionConditionalIterator 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_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 "< m_Multiplier*(double) previousSize) + { + // Raise lower threshold + if (m_Verbose) std::cout<<"Explosion detected, adapting lower threshold ..."<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 "< 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="<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 "< m_Multiplier* (double) previousSize) + { + // Lower upper threshold + if (m_Verbose) std::cout<<"Explosion detected, adapting upper threshold ..."<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 "< 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="<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 -- 2.47.1