[Insight-users] Re: region growing from local maxima
Gaetan Lehmann
gaetan.lehmann at jouy.inra.fr
Mon Oct 10 07:49:30 EDT 2005
Done now
However, I have some comments to do about review submission to the insight
journal:
+ The proposed title are not so good for code review
+ I have take care to change < and > to < and >, but the system
translate them to < and > which breaks the layout I tried to set
+ Only the first word of the title is saved
+ when review is saved, affiliation and expertise field are loose
+ there is no way to preview the review
and also, related to the insight journal: I no more receive new submission
email
Hope it can help to improve the journal :-)
Gaetan
On Thu, 06 Oct 2005 17:00:31 +0200, Miller, James V (Research)
<millerjv at crd.ge.com> wrote:
> Gaetan,
>
> Did you put these comments as a review to the submission to the Insight
> Software
> Journal?
>
> Jim
>
> -----Original Message-----
> From: insight-users-bounces+millerjv=crd.ge.com at itk.org
> [mailto:insight-users-bounces+millerjv=crd.ge.com at itk.org]On Behalf Of
> Gaetan Lehmann
> Sent: Thursday, October 06, 2005 4:21 AM
> To: Bryn
> Cc: insight-users at itk.org
> Subject: Re: [Insight-users] Re: region growing from local maxima
>
>
>
> 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