--- /dev/null
+/*=========================================================================
+
+ 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 TInputImage, class TOutputImage, class TtNorm=Function::Minimum<
+ typename TOutputImage::PixelType,
+ typename TOutputImage::PixelType,
+ typename TOutputImage::PixelType> >
+class ITK_EXPORT RelativePositionPropImageFilter :
+ public ImageToImageFilter< TInputImage, TOutputImage >
+{
+public:
+ /** Standard class typedefs. */
+ typedef RelativePositionPropImageFilter Self;
+ typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
+ typedef SmartPointer<Self> Pointer;
+ typedef SmartPointer<const Self> 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<typename TInputImage::IndexType , ImageDimension>
+ CorrespondanceMapType;
+ typedef float TabulationPixelType;
+ typedef typename itk::Image<TabulationPixelType , ImageDimension> TabulationImageType;
+
+
+ /** Image typedef support. */
+ typedef TInputImage InputImageType;
+ typedef TOutputImage OutputImageType;
+
+ typedef TtNorm TNormType;
+
+ typedef itk::Vector<double, ImageDimension> 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<OutputPixelType>));
+
+ /** 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
--- /dev/null
+/*=========================================================================
+
+ 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 <class TInputImage, class TOutputImage,class TtNorm>
+ void
+ RelativePositionPropImageFilter<TInputImage,TOutputImage,TtNorm>
+ ::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<b?a:b);
+ }
+
+ template< class TInputImage, class TOutputImage ,class TtNorm>
+ 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<InputImageType> InputConstIteratorWithIndexType;
+ typedef itk::ImageRegionIteratorWithIndex<OutputImageType> 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;i<ImageDimension;i++)
+ {
+ totalSize *= size[i];
+ }
+ const typename InputImageType::OffsetValueType * offsetTable = input->GetOffsetTable();
+
+
+
+ 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<pow((double)2,(int)ImageDimension);i++)
+ {
+ int orient[ImageDimension];
+ int init[ImageDimension];
+ typename InputImageType::OffsetType tempOffset;
+ tempOffset.Fill(0);
+
+ it.GoToBegin();
+ outputIt.GoToBegin();
+ for(int j=0;j<ImageDimension;j++)
+ {
+ orient[j] = ((int)(i/pow((double)2,j)))%2;
+ tempOffset[j] = (orient[j])*(size[j]-1);
+ }
+ it += tempOffset;
+
+ outputIt.SetIndex(it.GetIndex());
+
+ for(int j=0;j<ImageDimension;j++)
+ {
+ if(orient[j]==0)
+ init[j] = -size[j]+1;
+ else
+ init[j] = size[j]-1;
+ }
+
+ // DD(totalSize);
+ for(unsigned int p=1;p<=totalSize;p++)
+ {
+ // DD(p);
+ if ((m_VerboseProgress && (p % totalSize/100) == 0)) {
+ DD(p);
+ }
+ max = -100;
+ indmax = -1;
+ for (unsigned i = 0; i < it.Size(); i++)
+ {
+ if( it.GetPixel(i,in)!=nullIndex && in )
+ {
+ for(int j=0;j<ImageDimension;j++)
+ {
+ ind1[j] = it.GetPixel(i)[j];
+ vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
+ }
+
+ tmp = static_cast<double>(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;j<ImageDimension;j++)
+ {
+ if(p%offsetTable[j]==0)
+ {
+ if(orient[j]==0)
+ {
+ tempOffset[j] = 1;
+ for(int k=0;k<j;k++)
+ tempOffset[k] = init[k];
+ }
+ else
+ {
+ tempOffset[j] = -1;
+ for(int k=0;k<j;k++)
+ tempOffset[k] = init[k];
+ }
+ }
+ }
+ it += tempOffset;
+ outputIt.SetIndex(it.GetIndex());
+ }
+ }
+ }
+
+ // ==================================================================================================
+ else //if fast compute in two pass
+ {
+ // std::cout<<"pass 1"<<std::endl;
+
+ // DD(it.Size());
+
+ //pass 1
+ long i=0;
+ long nb = input->GetLargestPossibleRegion().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;j<ImageDimension;j++)
+ vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
+ // DD(vecttemp);
+
+ tmp = vecttemp.GetNorm();
+ if(tmp!=0)
+ tmp = acos((m_DirectionVector*vecttemp)/tmp);
+ //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
+ tmp=std::max(0.,1.-tmp/m_K1);
+ tmp = min( input->GetPixel(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"<<std::endl;
+ //pass2
+ it.GoToEnd();
+ --it;
+
+ i=0;
+ nb = this->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;j<ImageDimension;j++)
+ vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
+
+ tmp = vecttemp.GetNorm();
+
+ if(tmp!=0)
+ tmp = acos((m_DirectionVector*vecttemp)/tmp);
+ //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
+ tmp=std::max(0.,1. - tmp/m_K1);
+ tmp = min( input->GetPixel(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;i<ImageDimension;i++)
+ start[i]=0;
+
+ typename TabulationImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
+
+ for(register int i=0;i<ImageDimension;i++)
+ size[i]*=2;
+
+ typename TabulationImageType::RegionType region;
+
+ region.SetSize(size);
+ region.SetIndex(start);
+
+ m_AngleTabulation->SetRegions(region);
+
+ m_AngleTabulation->Allocate();
+ m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion());
+ typedef typename itk::ImageRegionIteratorWithIndex<TabulationImageType>
+ 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<ImageDimension;i++)
+ center[i]-=1;
+
+ VectorType vecttemp;
+
+ for(inputIt.GoToBegin();!inputIt.IsAtEnd();++inputIt)
+ {
+ typename TabulationImageType::IndexType idx = inputIt.GetIndex();
+ typename TabulationImageType::IndexType diff = idx - center;
+ for(int i=0;i<ImageDimension;i++)
+ vecttemp[i]=diff[i];
+ double tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
+ if(tmp==0)
+ {//std::cout<<"tem"<<std::endl;
+ inputIt.Set(0);}
+ else
+ {
+ double cos_angle = acos((m_DirectionVector*vecttemp)/tmp);
+ inputIt.Set(cos_angle);
+ }
+ }
+
+ typedef itk::ImageFileWriter<TabulationImageType> WT;
+ typename WT::Pointer wt = WT::New();
+ wt->SetFileName("testangle.nii");
+ wt->SetInput(m_AngleTabulation);
+ wt->Write();
+ std::cout<<"end compute angle"<<std::endl;
+
+ return m_AngleTabulation;
+ }
+
+ template< class TInputImage, class TOutputImage, class TtNorm>
+ 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
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff