[Insight-users] Re: region growing from local maxima
Bryn
blloyd at bwh.harvard.edu
Wed Oct 5 06:48:58 EDT 2005
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
More information about the Insight-users
mailing list