1 /*=========================================================================
3 Program: Insight Segmentation & Registration Toolkit
4 Module: $RCSfile: RelativePositionPropImageFilter.txx,v $
6 Date: $Date: 2010/07/12 06:57:25 $
7 Version: $Revision: 1.2 $
9 Copyright (c) Insight Software Consortium. All rights reserved.
10 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
12 This software is distributed WITHOUT ANY WARRANTY; without even
13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14 PURPOSE. See the above copyright notices for more information.
16 =========================================================================*/
17 #ifndef _RelativePositionPropImageFilter_txx
18 #define _RelativePositionPropImageFilter_txx
20 #include "RelativePositionPropImageFilter.h"
22 #include "itkNeighborhoodOperatorImageFilter.h"
23 #include "itkDerivativeOperator.h"
24 #include "itkZeroFluxNeumannBoundaryCondition.h"
25 #include "itkProgressAccumulator.h"
26 #include "itkImageFileWriter.h"
27 #include "itkSimpleContourExtractorImageFilter.h"
28 #include "itkUnaryFunctorImageFilter.h"
33 template <class TInputImage, class TOutputImage,class TtNorm>
35 RelativePositionPropImageFilter<TInputImage,TOutputImage,TtNorm>
36 ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
38 // call the superclass' implementation of this method
39 Superclass::GenerateInputRequestedRegion ();
41 // We need all the input.
42 typename InputImageType::Pointer input =
43 const_cast < InputImageType * >(this->GetInput ());
45 input->SetRequestedRegion (input->GetLargestPossibleRegion ());
49 * Generate all the output data
51 template < class TInputImage, class TOutputImage,class TtNorm>
52 void RelativePositionPropImageFilter < TInputImage,
53 TOutputImage ,TtNorm>::EnlargeOutputRequestedRegion (DataObject * output)
55 Superclass::EnlargeOutputRequestedRegion (output);
56 output->SetRequestedRegionToLargestPossibleRegion ();
60 inline float min(float a,float b)
65 template< class TInputImage, class TOutputImage ,class TtNorm>
67 RelativePositionPropImageFilter< TInputImage, TOutputImage ,TtNorm>
68 //::GenerateThreadedData(const typename TOutputImage::RegionType& outputRegionForThread, int threadId)
72 // std::cout <<"GenerateData" << std::endl;
74 this->AllocateOutputs();
76 typename InputImageType::ConstPointer input = this->GetInput();
78 typename CorrespondanceMapType::Pointer m_CorrespondanceMap = InitCorrespondanceMap();
80 typename InputImageType::IndexType nullIndex;
83 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputConstIteratorWithIndexType;
84 typedef itk::ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
86 typedef itk::NeighborhoodIterator< CorrespondanceMapType > NeighborhoodIteratorType;
89 typename NeighborhoodIteratorType::RadiusType radius;
90 radius.Fill(m_Radius);
91 NeighborhoodIteratorType it( radius,
93 m_CorrespondanceMap->GetLargestPossibleRegion());
94 OutputIteratorType outputIt( this->GetOutput(),
95 this->GetOutput()->GetLargestPossibleRegion() );
97 typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
98 unsigned int totalSize = 1;
99 for(int i=0;i<ImageDimension;i++)
101 totalSize *= size[i];
103 const typename InputImageType::OffsetValueType * offsetTable = input->GetOffsetTable();
111 typename InputImageType::IndexType ind1;
113 //unsigned int centeroffset = (int)(it.Size() / 2);
118 for(int i=0;i<pow((double)2,(int)ImageDimension);i++)
120 int orient[ImageDimension];
121 int init[ImageDimension];
122 typename InputImageType::OffsetType tempOffset;
126 outputIt.GoToBegin();
127 for(int j=0;j<ImageDimension;j++)
129 orient[j] = ((int)(i/pow((double)2,j)))%2;
130 tempOffset[j] = (orient[j])*(size[j]-1);
134 outputIt.SetIndex(it.GetIndex());
136 for(int j=0;j<ImageDimension;j++)
139 init[j] = -size[j]+1;
145 for(unsigned int p=1;p<=totalSize;p++)
148 if ((m_VerboseProgress && (p % totalSize/100) == 0)) {
153 for (unsigned i = 0; i < it.Size(); i++)
155 if( it.GetPixel(i,in)!=nullIndex && in )
157 for(int j=0;j<ImageDimension;j++)
159 ind1[j] = it.GetPixel(i)[j];
160 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
163 tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
166 if((m_DirectionVector*vecttemp)/tmp>=1)
169 if((m_DirectionVector*vecttemp)/tmp<=-1)
172 tmp = acos((m_DirectionVector*vecttemp)/tmp);
174 tmp=std::max(0.,1.-(tmp/m_K1));
175 // tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
177 tmp = m_TNorm( input->GetPixel(ind1), tmp);
188 it.SetCenterPixel(it.GetPixel(indmax));
191 else outputIt.Set( 0 );
199 for(int j=0;j<ImageDimension;j++)
201 if(p%offsetTable[j]==0)
207 tempOffset[k] = init[k];
213 tempOffset[k] = init[k];
218 outputIt.SetIndex(it.GetIndex());
223 // ==================================================================================================
224 else //if fast compute in two pass
226 // std::cout<<"pass 1"<<std::endl;
232 long nb = input->GetLargestPossibleRegion().GetNumberOfPixels( );
233 for (it.GoToBegin(), outputIt.GoToBegin(); ! it.IsAtEnd(); ++it, ++outputIt)
235 if (m_VerboseProgress && (i%(nb/10)) == 0) {
237 std::cout << i << " / " << nb << std::endl;
243 for (unsigned i = 0; i < it.Size(); i++) {
244 if( it.GetPixel(i,in)!=nullIndex && in )
247 // DS : can be precomputed ??? no
248 for(int j=0;j<ImageDimension;j++)
249 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
252 tmp = vecttemp.GetNorm();
254 tmp = acos((m_DirectionVector*vecttemp)/tmp);
255 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
256 tmp=std::max(0.,1.-tmp/m_K1);
257 tmp = min( input->GetPixel(it.GetPixel(i)), tmp);
268 it.SetCenterPixel(it.GetPixel(indmax));
276 // std::cout<<"pass 2"<<std::endl;
282 nb = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels( );
284 for ( outputIt.GoToReverseBegin(); ! outputIt.IsAtReverseEnd(); --it, --outputIt)
286 if (m_VerboseProgress && (i%(nb/10)) == 0) {
288 std::cout << i << " / " << nb << std::endl;
295 for (unsigned i = 0; i < it.Size(); i++)
296 if( it.GetPixel(i,in)!=nullIndex && in )
298 for(int j=0;j<ImageDimension;j++)
299 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
301 tmp = vecttemp.GetNorm();
304 tmp = acos((m_DirectionVector*vecttemp)/tmp);
305 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
306 tmp=std::max(0.,1. - tmp/m_K1);
307 tmp = min( input->GetPixel(it.GetPixel(i)),tmp);
318 it.SetCenterPixel(it.GetPixel(indmax));
328 template< class TInputImage, class TOutputImage, class TtNorm >
330 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
331 PrintSelf(std::ostream& os, Indent indent) const
333 Superclass::PrintSelf(os,indent);
335 os << indent << "First Angle: " << m_Alpha1 << std::endl;
336 os << indent << "Second Angle: " << m_Alpha2 << std::endl;
337 os << indent << "K1: " << m_K1 << std::endl;
340 template< class TInputImage, class TOutputImage, class TtNorm >
341 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::TabulationImageType::Pointer
342 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
343 ComputeAngleTabulation(){
345 typename TabulationImageType::Pointer m_AngleTabulation;
346 m_AngleTabulation = TabulationImageType::New();
348 typename TabulationImageType::IndexType start;
350 for(register int i=0;i<ImageDimension;i++)
353 typename TabulationImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
355 for(register int i=0;i<ImageDimension;i++)
358 typename TabulationImageType::RegionType region;
360 region.SetSize(size);
361 region.SetIndex(start);
363 m_AngleTabulation->SetRegions(region);
365 m_AngleTabulation->Allocate();
366 m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion());
367 typedef typename itk::ImageRegionIteratorWithIndex<TabulationImageType>
369 InputIteratorType inputIt = InputIteratorType(m_AngleTabulation, m_AngleTabulation->GetRequestedRegion());
371 typename TabulationImageType::IndexType requestedIndex =
372 m_AngleTabulation->GetRequestedRegion().GetIndex();
374 typename TabulationImageType::SizeType center = this->GetInput()->GetLargestPossibleRegion().GetSize();
375 for(register int i=0;i<ImageDimension;i++)
380 for(inputIt.GoToBegin();!inputIt.IsAtEnd();++inputIt)
382 typename TabulationImageType::IndexType idx = inputIt.GetIndex();
383 typename TabulationImageType::IndexType diff = idx - center;
384 for(int i=0;i<ImageDimension;i++)
386 double tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
388 {//std::cout<<"tem"<<std::endl;
392 double cos_angle = acos((m_DirectionVector*vecttemp)/tmp);
393 inputIt.Set(cos_angle);
397 typedef itk::ImageFileWriter<TabulationImageType> WT;
398 typename WT::Pointer wt = WT::New();
399 wt->SetFileName("testangle.nii");
400 wt->SetInput(m_AngleTabulation);
402 std::cout<<"end compute angle"<<std::endl;
404 return m_AngleTabulation;
407 template< class TInputImage, class TOutputImage, class TtNorm>
408 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::CorrespondanceMapType::Pointer
409 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
410 InitCorrespondanceMap()
412 typename CorrespondanceMapType::Pointer m_CorrespondanceMap;
413 m_CorrespondanceMap = CorrespondanceMapType::New();
414 m_CorrespondanceMap->SetRegions(this->GetInput()->GetLargestPossibleRegion());
415 m_CorrespondanceMap->Allocate();
417 typename InputImageType::IndexType nullIndex;
420 typedef itk::ImageRegionIterator< CorrespondanceMapType > CorrIteratorType;
421 CorrIteratorType it = CorrIteratorType(m_CorrespondanceMap,
422 m_CorrespondanceMap->GetLargestPossibleRegion());
424 typedef itk::ImageRegionConstIteratorWithIndex< InputImageType>
426 InputIteratorType inputIt = InputIteratorType(this->GetInput(),
427 this->GetInput()->GetLargestPossibleRegion());
429 for(it.GoToBegin(),inputIt.GoToBegin();!it.IsAtEnd();++it,++inputIt)
431 it.Set(inputIt.GetIndex());
438 return m_CorrespondanceMap;
449 } // end namespace itk