[Insight-users] Re: region growing from local maxima
Gaetan Lehmann
gaetan.lehmann at jouy.inra.fr
Thu Oct 6 04:21:16 EDT 2005
Hi,
I'm not sure this filter is as useful as a regional maxima filter, because
it will miss all maxima which are plateaus. But that's my point of view :-)
And this filter works with real types...
Some comments about the filter:
- The docstring for the threshold is the docstring for the radius
- with that filter, the user can choose the radius of the neighborhood,
but not the neighborhood. I don't know if it's useful to have a radius
greater than one, but I'm sure it's useful to be able to choose the
connectivity used (4-connected or 8-connected in 2D).
If be able to have a radius greater than 1 is useful, you should give the
user the ability to give a neighborhood; if it's not useful, a
SetFullyConnected() method should be used.
- the user should be able to choose the output image type
- You should give to the user the ability to choose the background and the
foreground values. The default foreground should be
NumericTraits<OutputPixelType>::max(), and the default background value
NumericTraits<OutputPixelType>::NonpositiveMin()
- you should cast your data in the PrintSelf method to avoid garbage if
pixel type is unsigned char:
os << indent << "Threshold: " << static_cast<typename
NumericTraits<InputPixelType>::PrintType>(m_Threshold) << std::endl;
Things are often better after a discussion, so I hope it can help to make
the filter better.
Regards,
Gaetan
On Wed, 05 Oct 2005 12:48:58 +0200, Bryn <blloyd at bwh.harvard.edu> wrote:
> Hi,
>
> I have written some code and posted it a while ago on the InsightJournal:
>
> http://www.insightsoftwareconsortium.org/InsightJournal/
>
>
> Since then I have rewritten it, inheriting from the ImageToImageFilter,
> permitting multithreading.
>
> Here is the code:
>
> itkLocalMaximumImageFilter.h:
> /*=========================================================================
>
>
> authors: Bryn A. Lloyd, Simon K. Warfield, Computational Radiology
> Laborotory (CRL), Brigham and Womans
> date: 06/30/2005
> Acknowledgements:
> This investigation was supported in part by NSF ITR 0426558 and
> NIH grants R21 MH67054 and P41 RR13218.
>
> =========================================================================*/
> #ifndef __itkLocalMaximumImageFilter_h
> #define __itkLocalMaximumImageFilter_h
>
> #include "itkImageToImageFilter.h"
> #include "itkImage.h"
>
> #include "itkImageRegionConstIteratorWithIndex.h"
> #include "itkImageRegionIteratorWithIndex.h"
>
>
>
>
> namespace itk
> {
>
>
> /** \class LocalMaximumImageFilter
> * \brief Calculates the Local Maximum values and puts them in a binary
> image and a point-set/ mesh.
> *
> * This filter finds the local maxima of an image. The result is
> returned as a PointSet: GetLocalMaximaPointSet( void ).
> *
> * Additionally a binary image, with 0 for background, and 1 for the
> points can be retrieved: GetOutput()
> *
> * This filter assumes a pixel value must be larger than all of it's
> neighbors to be a maximum.
> * The neighborhood size can be set by setting the Radius. The second
> parameter is a threshold.
> * Only maxima above the threshold are selected.
> *
> * TODO: allow user to set which information should be saved in the
> point-set (i.e. index or physical point, data)
> *
> *
> * \sa Image
> * \sa Neighborhood
> * \sa NeighborhoodOperator
> * \sa NeighborhoodIterator
> *
> * \ingroup IntensityImageFilters
> */
> template <class TInputImage, class TOutputMesh>
> class ITK_EXPORT LocalMaximumImageFilter :
> public ImageToImageFilter<TInputImage, TInputImage>
> {
> public:
>
>
> /** Run-time type information (and related methods). */
> itkTypeMacro(LocalMaximumImageFilter, ImageToImageFilter);
>
>
> /** Extract dimension from input and output image. */
> itkStaticConstMacro(InputImageDimension, unsigned int,
> TInputImage::ImageDimension);
>
> /** Some typedefs associated with the input images. */
> typedef TInputImage
> InputImageType;
> typedef typename InputImageType::Pointer
> InputImagePointer;
> typedef typename InputImageType::ConstPointer
> InputImageConstPointer;
> typedef typename InputImageType::RegionType InputImageRegionType;
> typedef typename InputImageType::RegionType
> OutputImageRegionType;
> typedef ImageRegionConstIteratorWithIndex<InputImageType>
> InputImageIterator;
> typedef ImageRegionIteratorWithIndex<InputImageType>
> OutputImageIterator;
>
> /** Some typedefs associated with the output mesh. */
> typedef TOutputMesh OutputMeshType;
> typedef typename OutputMeshType::PointType PointType;
> typedef typename OutputMeshType::Pointer OutputMeshPointer;
> typedef typename OutputMeshType::ConstPointer
> OutputMeshConstPointer;
> typedef typename OutputMeshType::PointsContainer PointsContainer;
> typedef typename PointsContainer::Pointer
> PointsContainerPointer;
> typedef typename OutputMeshType::PointDataContainer PointDataContainer;
> typedef typename PointDataContainer::Pointer
> PointDataContainerPointer;
>
>
> /** Image typedef support. */
> typedef typename InputImageType::PixelType InputPixelType;
> typedef typename InputImageType::SizeType InputSizeType;
> typedef typename InputImageType::IndexType InputIndexType;
>
> /** Standard class typedefs. */
> typedef LocalMaximumImageFilter Self;
> typedef ImageToImageFilter< InputImageType, InputImageType> Superclass;
> typedef SmartPointer<Self> Pointer;
> typedef SmartPointer<const Self> ConstPointer;
>
>
> /** Method for creation through the object factory. */
> itkNewMacro(Self);
>
> /** Set the radius of the neighborhood used to compute the mean. */
> itkSetMacro(Radius, InputSizeType);
> /** Get the radius of the neighborhood used to compute the mean */
> itkGetConstReferenceMacro(Radius, InputSizeType);
>
> /** Set the radius of the neighborhood used to compute the mean. */
> itkSetMacro(Threshold, InputPixelType);
> /** Get the radius of the neighborhood used to compute the mean */
> itkGetConstMacro(Threshold, InputPixelType);
>
>
> /** get the points as point set */
> OutputMeshType * GetLocalMaximaPointSet( void );
>
> virtual void GenerateInputRequestedRegion()
> throw(InvalidRequestedRegionError);
> protected:
> LocalMaximumImageFilter();
> virtual ~LocalMaximumImageFilter() {}
> void PrintSelf(std::ostream& os, Indent indent) const;
>
> /** LocalMaximumImageFilter can be implemented as a multithreaded
> filter.
> * Therefore, this implementation provides a ThreadedGenerateData()
> * routine which is called for each processing thread. The output
> * image data is allocated automatically by the superclass prior to
> * calling ThreadedGenerateData(). ThreadedGenerateData can only
> * write to the portion of the output image specified by the
> * parameter "outputRegionForThread"
> *
> * \sa ImageToImageFilter::ThreadedGenerateData(),
> * ImageToImageFilter::GenerateData() */
> void ThreadedGenerateData(const OutputImageRegionType&
> outputRegionForThread,
> int threadId );
>
> void AfterThreadedGenerateData ();
>
>
> private:
> LocalMaximumImageFilter(const Self&); //purposely not implemented
> void operator=(const Self&); //purposely not implemented
>
> InputSizeType m_Radius;
> InputPixelType m_Threshold;
>
> OutputMeshPointer m_PointSet;
> };
>
> } // end namespace itk
>
> #ifndef ITK_MANUAL_INSTANTIATION
> #include "itkLocalMaximumImageFilter.txx"
> #endif
>
> #endif
>
>
>
>
>
>
>
> itkLocalMaximumImageFilter.txx:
>
>
> #ifndef _itkLocalMaximumImageFilter_txx
> #define _itkLocalMaximumImageFilter_txx
> #include "itkLocalMaximumImageFilter.h"
>
> #include "itkConstNeighborhoodIterator.h"
> #include "itkNeighborhoodInnerProduct.h"
>
> #include "itkImageRegionIterator.h"
> #include "itkNeighborhoodAlgorithm.h"
> #include "itkZeroFluxNeumannBoundaryCondition.h"
> #include "itkConstantBoundaryCondition.h"
>
> #include "itkOffset.h"
> #include "itkProgressReporter.h"
>
> #include "itkNumericTraits.h"
>
>
> namespace itk
> {
>
> template <class TInputImage, class TOutputMesh>
> LocalMaximumImageFilter< TInputImage, TOutputMesh>
> ::LocalMaximumImageFilter()
> {
>
> // Modify superclass default values, can be overridden by subclasses
> this->SetNumberOfRequiredInputs(1);
> /*
> PointDataContainerPointer pointData = PointDataContainer::New();
> OutputMeshPointer mesh = OutputMeshType::New();
> mesh->SetPointData( pointData.GetPointer() );
> */ this->SetNumberOfRequiredOutputs( 1 );
>
> this->m_Radius.Fill(1);
> this->m_Threshold = 0.0;
> }
>
>
>
>
> template <class TInputImage, class TOutputMesh>
> typename LocalMaximumImageFilter< TInputImage,
> TOutputMesh>::OutputMeshType *
> LocalMaximumImageFilter< TInputImage, TOutputMesh>
> ::GetLocalMaximaPointSet()
> {
> return this->m_PointSet;
> }
>
>
>
> template <class TInputImage, class TOutputMesh>
> void
> LocalMaximumImageFilter< TInputImage, TOutputMesh>
> ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
> {
> // call the superclass' implementation of this method
> Superclass::GenerateInputRequestedRegion();
> // get pointers to the input and output
> typename Superclass::InputImagePointer inputPtr =
> const_cast< TInputImage * >( this->GetInput() );
> typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
> if ( !inputPtr || !outputPtr )
> {
> return;
> }
>
> // get a copy of the input requested region (should equal the output
> // requested region)
> typename TInputImage::RegionType inputRequestedRegion;
> inputRequestedRegion = inputPtr->GetRequestedRegion();
>
> // pad the input requested region by the operator radius
> inputRequestedRegion.PadByRadius( m_Radius );
>
> // crop the input requested region at the input's largest possible
> region
> if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
> {
> inputPtr->SetRequestedRegion( inputRequestedRegion );
> return;
> }
> else
> {
> // Couldn't crop the region (requested region is outside the largest
> // possible region). Throw an exception.
>
> // store what we tried to request (prior to trying to crop)
> inputPtr->SetRequestedRegion( inputRequestedRegion );
> // build an exception
> InvalidRequestedRegionError e(__FILE__, __LINE__);
> OStringStream msg;
> msg << static_cast<const char *>(this->GetNameOfClass())
> << "::GenerateInputRequestedRegion()";
> e.SetLocation(msg.str().c_str());
> e.SetDescription("Requested region is (at least partially) outside
> the largest possible region.");
> e.SetDataObject(inputPtr);
> throw e;
> }
> }
>
>
>
> template <class TInputImage, class TOutputMesh>
> void
> LocalMaximumImageFilter< TInputImage, TOutputMesh>
> ::ThreadedGenerateData(const OutputImageRegionType&
> outputRegionForThread,
> int threadId )
> {
> unsigned int i;
> ConstantBoundaryCondition<InputImageType> cbc;
> cbc.SetConstant( NumericTraits<InputPixelType>::NonpositiveMin() );
>
> ConstNeighborhoodIterator<InputImageType> bit;
> ImageRegionIterator<InputImageType> it;
>
> InputImageConstPointer input = this->GetInput(0);
> InputImagePointer outputImage = this->GetOutput();
> //Set background value for binary local maxima image
> //outputImage->FillBuffer( 0 );
> // Find the data-set boundary "faces"
> typename
> NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType
> faceList;
> NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
> faceList = bC(input, outputRegionForThread, m_Radius);
>
> typename
> NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
> fit;
>
> // support progress methods/callbacks
> ProgressReporter progress(this, threadId,
> outputRegionForThread.GetNumberOfPixels());
>
> // Process each of the boundary faces. These are N-d regions which
> border
> // the edge of the buffer.
> for (fit=faceList.begin(); fit != faceList.end(); ++fit)
> {
> bit = ConstNeighborhoodIterator<InputImageType>(m_Radius,
> input, *fit);
> unsigned int neighborhoodSize = bit.Size();
> it = ImageRegionIterator<InputImageType>(outputImage, *fit);
> bit.OverrideBoundaryCondition(&cbc);
> bit.GoToBegin();
>
> while ( ! bit.IsAtEnd() )
> {
> bool isMaximum = true;
> InputPixelType centerValue = bit.GetCenterPixel();
> //NumericTraits<InputRealType>::NonpositiveMin();
> for (i = 0; i < neighborhoodSize; ++i)
> {
> InputPixelType tmp = bit.GetPixel(i);
> // if we find a neighbor with a larger value than the
> center, tthe center is not a maximum...
> if (tmp > centerValue) {
> isMaximum = false;
> break;
> }
> }
>
> if (isMaximum & (centerValue>m_Threshold))
> {
> it.Set( 1 );
> }
> ++bit;
> ++it;
> progress.CompletedPixel();
> }
> }
> }
>
>
> template <class TInputImage, class TOutputMesh>
> void
> LocalMaximumImageFilter< TInputImage, TOutputMesh>
> ::AfterThreadedGenerateData()
> {
>
> InputImageConstPointer input = this->GetInput(0);
> InputImagePointer outputImage = this->GetOutput();
> m_PointSet = OutputMeshType::New();
> PointsContainerPointer points = PointsContainer::New();
> PointDataContainerPointer data = PointDataContainer::New();
> ImageRegionConstIterator<InputImageType> it(input,
> input->GetBufferedRegion());
> ImageRegionConstIterator<InputImageType> ot(outputImage,
> outputImage->GetBufferedRegion());
> for (it.GoToBegin(), ot.GoToBegin(); !it.IsAtEnd(); ++it, ++ot) {
> if (ot.Value() > 0 ) {
> InputIndexType index = ot.GetIndex();
> PointType point;
> input->TransformIndexToPhysicalPoint( index , point );
> data->push_back( it.Value() );
> points->push_back( point );
> }
> }
> m_PointSet->SetPoints( points );
> m_PointSet->SetPointData( data );
> }
>
>
>
>
>
> template <class TInputImage, class TOutputMesh>
> void
> LocalMaximumImageFilter< TInputImage, TOutputMesh>
> ::PrintSelf(
> std::ostream& os,
> Indent indent) const
> {
> Superclass::PrintSelf( os, indent );
> os << indent << "Radius: " << m_Radius << std::endl;
> os << indent << "Threshold: " << m_Threshold << std::endl;
>
> }
>
> } // end namespace itk
>
> #endif
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
--
Gaëtan Lehmann
Biologie du Développement et de la Reproduction
INRA de Jouy-en-Josas (France)
tel: +33 1 34 65 29 66 fax: 01 34 65 29 09
http://voxel.jouy.inra.fr
More information about the Insight-users
mailing list