From f5fbb8b3c473ceced4864212387ebe26929c4696 Mon Sep 17 00:00:00 2001 From: dsarrut Date: Wed, 30 Jun 2010 05:58:56 +0000 Subject: [PATCH] add Bloch's filter (relative position) --- itk/RelativePositionPropImageFilter.h | 216 ++++++++++++ itk/RelativePositionPropImageFilter.txx | 451 ++++++++++++++++++++++++ 2 files changed, 667 insertions(+) create mode 100644 itk/RelativePositionPropImageFilter.h create mode 100644 itk/RelativePositionPropImageFilter.txx diff --git a/itk/RelativePositionPropImageFilter.h b/itk/RelativePositionPropImageFilter.h new file mode 100644 index 0000000..e7f7da5 --- /dev/null +++ b/itk/RelativePositionPropImageFilter.h @@ -0,0 +1,216 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: RelativePositionPropImageFilter.h,v $ + Language: C++ + Date: $Date: 2010/06/30 05:58:56 $ + Version: $Revision: 1.1 $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __RelativePositionPropImageFilter_h +#define __RelativePositionPropImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkImage.h" +#include "itkConceptChecking.h" + +#include "itkPointSet.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkMinimumImageFilter.h" + +namespace itk +{ + + +/** \class RelativePositionPropImageFilter + * \brief Compute the fuzzy subset of an image which satisfies some directional relative position. + * \author Jamal Atif and Olivier Nempont + * + * This filter computes a fuzzy subset of an image satisfying a particular directionnal relative position from an object (crisp or fuzzy). + * + * + * \par INPUT / OUTPUT + * This filter takes a crisp or a fuzzy object as input. + * In fuzzy case, the values have to be defined between 0 and 1. + * + * The result is a fuzzy subset which values are defined between + * 0 if the relation isn't fulfilled in this point to 1 is the relation is + * fully satisfied. + * WARNING: the output image type as to be decimal. + * + * \par PARAMETERS + * \par + * The Alpha1 and Alpha2 parameters are used to specify the direction. + * Alpha1 is the angle in 'xy' plane from 'x' unit vector. + * Alpha2 is used in 3D to specify the angle with 'xy' plane + * + * \par + * K is an opening parameter. Higher value enlarge the support of the result. + * By default it is fixed at PI/2 + * + * \par REFERENCE + * Fuzzy Relative Position Between Objects in Image Processing: A Morphological Approach + * Isabelle Bloch + * IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 21, NO. 7, JULY 1999 + * + * This filter is implemented using the propagation algorithm + */ + +template > +class ITK_EXPORT RelativePositionPropImageFilter : + public ImageToImageFilter< TInputImage, TOutputImage > +{ +public: + /** Standard class typedefs. */ + typedef RelativePositionPropImageFilter Self; + typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Extract some information from the image types. Dimensionality + * of the two images is assumed to be the same. */ + typedef typename TOutputImage::PixelType OutputPixelType; + typedef typename TOutputImage::InternalPixelType OutputInternalPixelType; + typedef typename TInputImage::PixelType InputPixelType; + typedef typename TInputImage::InternalPixelType InputInternalPixelType; + + /** Extract some information from the image types. Dimensionality + * of the two images is assumed to be the same. */ + itkStaticConstMacro(ImageDimension, unsigned int, + TOutputImage::ImageDimension); + + + typedef typename itk::Image + CorrespondanceMapType; + typedef float TabulationPixelType; + typedef typename itk::Image TabulationImageType; + + + /** Image typedef support. */ + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + typedef TtNorm TNormType; + + typedef itk::Vector VectorType; + + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RelativePositionPropImageFilter, ImageToImageFilter); + + /** The output pixel type must be signed. */ + itkConceptMacro(SignedOutputPixelType, (Concept::Signed)); + + /** Standard get/set macros for filter parameters. */ + + + itkSetMacro(Alpha1, double); + itkGetMacro(Alpha1, double); + itkSetMacro(Alpha2, double); + itkGetMacro(Alpha2, double); + + itkSetMacro(K1, double); + itkGetMacro(K1, double); +// itkSetMacro(K2, double); +// itkGetMacro(K2, double); + + itkSetMacro(Radius, double); + itkGetMacro(Radius, double); + + itkSetMacro(VerboseProgress, bool); + itkGetMacro(VerboseProgress, bool); + itkBooleanMacro(VerboseProgress); + + itkSetMacro(Fast, bool); + itkGetMacro(Fast, bool); + itkBooleanMacro(Fast); + + void computeDirection() + { + switch(ImageDimension) + { +// case 2: // DS comment to avoid warning +// m_DirectionVector[0]=cos(m_Alpha1); +// m_DirectionVector[1]=sin(m_Alpha1); +// break; + case 3: + m_DirectionVector[0]=cos(m_Alpha1)*cos(m_Alpha2); + m_DirectionVector[1]=cos(m_Alpha2)*sin(m_Alpha1); + m_DirectionVector[2]=sin(m_Alpha2); + break; + } + } + + + virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError); + void EnlargeOutputRequestedRegion (DataObject * output); + +protected: + RelativePositionPropImageFilter() + { + m_Alpha1 = 0; + m_Alpha2 = 0; + m_K1 = vcl_acos(-1.0)/2; + // m_K2 = 3.1417/2; + m_Radius = 2; // DS + m_Fast = true; // DS + m_VerboseProgress = false; + } + virtual ~RelativePositionPropImageFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + //void GenerateThreadedData(const typename TOutputImage::RegionType& outputRegionForThread, int threadId); + void GenerateData(); + +private: + RelativePositionPropImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + + /** The angles*/ + double m_Alpha1; + double m_Alpha2; + double m_K1; + // double m_K2; + + unsigned int m_Radius; + TNormType m_TNorm; + bool m_VerboseProgress; + + VectorType m_DirectionVector; + + /** + * 2 pass instead of 2^NDimension. Warning this may cause some artifacts + */ + bool m_Fast; + + //allocation et initialisation de la carte de correspondance + typename CorrespondanceMapType::Pointer InitCorrespondanceMap(); + + //compute the tabulation map + typename TabulationImageType::Pointer ComputeAngleTabulation(); + + +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "RelativePositionPropImageFilter.txx" +#endif + +#endif diff --git a/itk/RelativePositionPropImageFilter.txx b/itk/RelativePositionPropImageFilter.txx new file mode 100644 index 0000000..56a27d4 --- /dev/null +++ b/itk/RelativePositionPropImageFilter.txx @@ -0,0 +1,451 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: RelativePositionPropImageFilter.txx,v $ + Language: C++ + Date: $Date: 2010/06/30 05:58:56 $ + Version: $Revision: 1.1 $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + + =========================================================================*/ +#ifndef _RelativePositionPropImageFilter_txx +#define _RelativePositionPropImageFilter_txx +#include "RelativePositionPropImageFilter.h" + +#include "itkNeighborhoodOperatorImageFilter.h" +#include "itkDerivativeOperator.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkProgressAccumulator.h" +#include "itkImageFileWriter.h" + +#include "itkSimpleContourExtractorImageFilter.h" +#include "itkUnaryFunctorImageFilter.h" + +namespace itk +{ + + template + void + RelativePositionPropImageFilter + ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError) + { + // call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion (); + + // We need all the input. + typename InputImageType::Pointer input = + const_cast < InputImageType * >(this->GetInput ()); + + input->SetRequestedRegion (input->GetLargestPossibleRegion ()); + } + + /** + * Generate all the output data + */ + template < class TInputImage, class TOutputImage,class TtNorm> + void RelativePositionPropImageFilter < TInputImage, + TOutputImage ,TtNorm>::EnlargeOutputRequestedRegion (DataObject * output) + { + Superclass::EnlargeOutputRequestedRegion (output); + output->SetRequestedRegionToLargestPossibleRegion (); + } + + + inline float min(float a,float b) + { + return (a + void + RelativePositionPropImageFilter< TInputImage, TOutputImage ,TtNorm> + //::GenerateThreadedData(const typename TOutputImage::RegionType& outputRegionForThread, int threadId) + ::GenerateData() + { + + // std::cout <<"GenerateData" << std::endl; + + this->AllocateOutputs(); + computeDirection(); + typename InputImageType::ConstPointer input = this->GetInput(); + + typename CorrespondanceMapType::Pointer m_CorrespondanceMap = InitCorrespondanceMap(); + + typename InputImageType::IndexType nullIndex; + nullIndex.Fill(-1); + + typedef itk::ImageRegionConstIteratorWithIndex InputConstIteratorWithIndexType; + typedef itk::ImageRegionIteratorWithIndex OutputIteratorType; + + typedef itk::NeighborhoodIterator< CorrespondanceMapType > NeighborhoodIteratorType; + + + typename NeighborhoodIteratorType::RadiusType radius; + radius.Fill(m_Radius); + NeighborhoodIteratorType it( radius, + m_CorrespondanceMap, + m_CorrespondanceMap->GetLargestPossibleRegion()); + OutputIteratorType outputIt( this->GetOutput(), + this->GetOutput()->GetLargestPossibleRegion() ); + + typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize(); + unsigned int totalSize = 1; + for(int i=0;iGetOffsetTable(); + + + + VectorType vecttemp; + OutputPixelType max; + int indmax; + OutputPixelType tmp; + typename InputImageType::IndexType ind1; + bool in=false; + //unsigned int centeroffset = (int)(it.Size() / 2); + + if(!m_Fast) + { + DD(m_Fast); + for(int i=0;i(sqrt(vecttemp*vecttemp)) ; + if(tmp!=0) + { + if((m_DirectionVector*vecttemp)/tmp>=1) + tmp = 0; + else + if((m_DirectionVector*vecttemp)/tmp<=-1) + tmp = 4.*atan(1.0); + else + tmp = acos((m_DirectionVector*vecttemp)/tmp); + } + tmp=std::max(0.,1.-(tmp/m_K1)); + // tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K; + + tmp = m_TNorm( input->GetPixel(ind1), tmp); + if ( tmp > max) + { + max = tmp; + indmax = i; + } + } + } + + if(indmax!=-1 ) + { + it.SetCenterPixel(it.GetPixel(indmax)); + if( max>0) + outputIt.Set( max ); + else outputIt.Set( 0 ); + } + else + { + outputIt.Set( 0 ); + } + + tempOffset.Fill(0); + for(int j=0;jGetLargestPossibleRegion().GetNumberOfPixels( ); + for (it.GoToBegin(), outputIt.GoToBegin(); ! it.IsAtEnd(); ++it, ++outputIt) + { + if (m_VerboseProgress && (i%(nb/10)) == 0) { + //DD(i); + std::cout << i << " / " << nb << std::endl; + } + i++; + + max = -100; + indmax = -1; + for (unsigned i = 0; i < it.Size(); i++) { + if( it.GetPixel(i,in)!=nullIndex && in ) + { + + // DS : can be precomputed ??? no + for(int j=0;jGetPixel(it.GetPixel(i)), tmp); + if ( tmp > max) + { + max = tmp; + indmax = i; + } + } + } + + if(max>0) + { + it.SetCenterPixel(it.GetPixel(indmax)); + outputIt.Set( max ); + } + else + outputIt.Set( 0 ); + + } + + // std::cout<<"pass 2"<GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels( ); + + for ( outputIt.GoToReverseBegin(); ! outputIt.IsAtReverseEnd(); --it, --outputIt) + { + if (m_VerboseProgress && (i%(nb/10)) == 0) { + //DD(i); + std::cout << i << " / " << nb << std::endl; + } + i++; + + max = -100; + indmax = -1; + + for (unsigned i = 0; i < it.Size(); i++) + if( it.GetPixel(i,in)!=nullIndex && in ) + { + for(int j=0;jGetPixel(it.GetPixel(i)),tmp); + + if ( tmp > max) + { + max = tmp; + indmax = i; + } + } + + if(max>0) + { + it.SetCenterPixel(it.GetPixel(indmax)); + outputIt.Set( max ); + } + else + outputIt.Set( 0 ); + } + } + + } + + template< class TInputImage, class TOutputImage, class TtNorm > + void + RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >:: + PrintSelf(std::ostream& os, Indent indent) const + { + Superclass::PrintSelf(os,indent); + + os << indent << "First Angle: " << m_Alpha1 << std::endl; + os << indent << "Second Angle: " << m_Alpha2 << std::endl; + os << indent << "K1: " << m_K1 << std::endl; + } + + template< class TInputImage, class TOutputImage, class TtNorm > + typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::TabulationImageType::Pointer + RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >:: + ComputeAngleTabulation(){ + computeDirection(); + typename TabulationImageType::Pointer m_AngleTabulation; + m_AngleTabulation = TabulationImageType::New(); + + typename TabulationImageType::IndexType start; + + for(register int i=0;iGetInput()->GetLargestPossibleRegion().GetSize(); + + for(register int i=0;iSetRegions(region); + + m_AngleTabulation->Allocate(); + m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion()); + typedef typename itk::ImageRegionIteratorWithIndex + InputIteratorType; + InputIteratorType inputIt = InputIteratorType(m_AngleTabulation, m_AngleTabulation->GetRequestedRegion()); + + typename TabulationImageType::IndexType requestedIndex = + m_AngleTabulation->GetRequestedRegion().GetIndex(); + + typename TabulationImageType::SizeType center = this->GetInput()->GetLargestPossibleRegion().GetSize(); + for(register int i=0;i(sqrt(vecttemp*vecttemp)) ; + if(tmp==0) + {//std::cout<<"tem"< WT; + typename WT::Pointer wt = WT::New(); + wt->SetFileName("testangle.nii"); + wt->SetInput(m_AngleTabulation); + wt->Write(); + std::cout<<"end compute angle"< + typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::CorrespondanceMapType::Pointer + RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >:: + InitCorrespondanceMap() + { + typename CorrespondanceMapType::Pointer m_CorrespondanceMap; + m_CorrespondanceMap = CorrespondanceMapType::New(); + m_CorrespondanceMap->SetRegions(this->GetInput()->GetLargestPossibleRegion()); + m_CorrespondanceMap->Allocate(); + + typename InputImageType::IndexType nullIndex; + nullIndex.Fill(-1); + + typedef itk::ImageRegionIterator< CorrespondanceMapType > CorrIteratorType; + CorrIteratorType it = CorrIteratorType(m_CorrespondanceMap, + m_CorrespondanceMap->GetLargestPossibleRegion()); + + typedef itk::ImageRegionConstIteratorWithIndex< InputImageType> + InputIteratorType; + InputIteratorType inputIt = InputIteratorType(this->GetInput(), + this->GetInput()->GetLargestPossibleRegion()); + + for(it.GoToBegin(),inputIt.GoToBegin();!it.IsAtEnd();++it,++inputIt) + if(inputIt.Get()>0) + it.Set(inputIt.GetIndex()); + else + it.Set(nullIndex); + + + + + return m_CorrespondanceMap; + } + + + + + + + + + +} // end namespace itk + +#endif -- 2.47.1