1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================*/
18 /*=========================================================================
20 Program: Insight Segmentation & Registration Toolkit
21 Module: $RCSfile: RelativePositionPropImageFilter.txx,v $
23 Date: $Date: 2010/07/12 06:57:25 $
24 Version: $Revision: 1.2 $
26 Copyright (c) Insight Software Consortium. All rights reserved.
27 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
29 This software is distributed WITHOUT ANY WARRANTY; without even
30 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
31 PURPOSE. See the above copyright notices for more information.
33 =========================================================================*/
34 #ifndef _RelativePositionPropImageFilter_txx
35 #define _RelativePositionPropImageFilter_txx
37 #include "RelativePositionPropImageFilter.h"
39 #include "itkNeighborhoodOperatorImageFilter.h"
40 #include "itkDerivativeOperator.h"
41 #include "itkZeroFluxNeumannBoundaryCondition.h"
42 #include "itkProgressAccumulator.h"
43 #include "itkImageFileWriter.h"
44 #include "itkSimpleContourExtractorImageFilter.h"
45 #include "itkUnaryFunctorImageFilter.h"
50 template <class TInputImage, class TOutputImage,class TtNorm>
52 RelativePositionPropImageFilter<TInputImage,TOutputImage,TtNorm>
53 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
54 ::GenerateInputRequestedRegion()
56 ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
59 // call the superclass' implementation of this method
60 Superclass::GenerateInputRequestedRegion ();
62 // We need all the input.
63 typename InputImageType::Pointer input =
64 const_cast < InputImageType * >(this->GetInput ());
66 input->SetRequestedRegion (input->GetLargestPossibleRegion ());
70 * Generate all the output data
72 template < class TInputImage, class TOutputImage,class TtNorm>
73 void RelativePositionPropImageFilter < TInputImage,
74 TOutputImage ,TtNorm>::EnlargeOutputRequestedRegion (DataObject * output)
76 Superclass::EnlargeOutputRequestedRegion (output);
77 output->SetRequestedRegionToLargestPossibleRegion ();
81 inline float min(float a,float b)
86 template< class TInputImage, class TOutputImage ,class TtNorm>
88 RelativePositionPropImageFilter< TInputImage, TOutputImage ,TtNorm>
89 //::GenerateThreadedData(const typename TOutputImage::RegionType& outputRegionForThread, int threadId)
93 this->AllocateOutputs();
95 typename InputImageType::ConstPointer input = this->GetInput();
97 typename CorrespondanceMapType::Pointer m_CorrespondanceMap = InitCorrespondanceMap();
99 typename InputImageType::IndexType nullIndex;
102 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputConstIteratorWithIndexType;
103 typedef itk::ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
105 typedef itk::NeighborhoodIterator< CorrespondanceMapType > NeighborhoodIteratorType;
108 typename NeighborhoodIteratorType::RadiusType radius;
109 radius.Fill(m_Radius);
110 NeighborhoodIteratorType it( radius,
112 m_CorrespondanceMap->GetLargestPossibleRegion());
113 OutputIteratorType outputIt( this->GetOutput(),
114 this->GetOutput()->GetLargestPossibleRegion() );
116 typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
117 unsigned int totalSize = 1;
118 for(int i=0;i<ImageDimension;i++)
120 totalSize *= size[i];
122 const typename InputImageType::OffsetValueType * offsetTable = input->GetOffsetTable();
130 typename InputImageType::IndexType ind1;
132 //unsigned int centeroffset = (int)(it.Size() / 2);
136 for(int i=0;i<pow((double)2,(int)ImageDimension);i++)
138 int orient[ImageDimension];
139 int init[ImageDimension];
140 typename InputImageType::OffsetType tempOffset;
144 outputIt.GoToBegin();
145 for(int j=0;j<ImageDimension;j++)
147 orient[j] = ((int)(i/pow((double)2,j)))%2;
148 tempOffset[j] = (orient[j])*(size[j]-1);
152 outputIt.SetIndex(it.GetIndex());
154 for(int j=0;j<ImageDimension;j++)
157 init[j] = -size[j]+1;
163 for(unsigned int p=1;p<=totalSize;p++)
166 if ((m_VerboseProgress && (p % totalSize/100) == 0)) {
171 for (unsigned i = 0; i < it.Size(); i++)
173 if( it.GetPixel(i,in)!=nullIndex && in )
175 for(int j=0;j<ImageDimension;j++)
177 ind1[j] = it.GetPixel(i)[j];
178 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
181 tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
184 if((m_DirectionVector*vecttemp)/tmp>=1)
187 if((m_DirectionVector*vecttemp)/tmp<=-1)
190 tmp = acos((m_DirectionVector*vecttemp)/tmp);
192 tmp=std::max(0.,1.-(tmp/m_K1));
193 // tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
195 tmp = m_TNorm( input->GetPixel(ind1), tmp);
206 it.SetCenterPixel(it.GetPixel(indmax));
209 else outputIt.Set( 0 );
217 for(int j=0;j<ImageDimension;j++)
219 if(p%offsetTable[j]==0)
225 tempOffset[k] = init[k];
231 tempOffset[k] = init[k];
236 outputIt.SetIndex(it.GetIndex());
241 // ==================================================================================================
242 else //if fast compute in two pass
244 // std::cout<<"pass 1"<<std::endl;
250 long nb = input->GetLargestPossibleRegion().GetNumberOfPixels( );
251 for (it.GoToBegin(), outputIt.GoToBegin(); ! it.IsAtEnd(); ++it, ++outputIt)
253 if (m_VerboseProgress && (i%(nb/10)) == 0) {
255 std::cout << i << " / " << nb << std::endl;
261 for (unsigned i = 0; i < it.Size(); i++) {
262 if( it.GetPixel(i,in)!=nullIndex && in )
265 // DS : can be precomputed ??? no
266 for(int j=0;j<ImageDimension;j++)
267 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
270 tmp = vecttemp.GetNorm();
272 tmp = acos((m_DirectionVector*vecttemp)/tmp);
273 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
274 tmp=std::max(0.,1.-tmp/m_K1);
275 tmp = min( input->GetPixel(it.GetPixel(i)), tmp);
286 it.SetCenterPixel(it.GetPixel(indmax));
294 // std::cout<<"pass 2"<<std::endl;
300 nb = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels( );
302 for ( outputIt.GoToReverseBegin(); ! outputIt.IsAtReverseEnd(); --it, --outputIt)
304 if (m_VerboseProgress && (i%(nb/10)) == 0) {
306 std::cout << i << " / " << nb << std::endl;
313 for (unsigned i = 0; i < it.Size(); i++)
314 if( it.GetPixel(i,in)!=nullIndex && in )
316 for(int j=0;j<ImageDimension;j++)
317 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
319 tmp = vecttemp.GetNorm();
322 tmp = acos((m_DirectionVector*vecttemp)/tmp);
323 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
324 tmp=std::max(0.,1. - tmp/m_K1);
325 tmp = min( input->GetPixel(it.GetPixel(i)),tmp);
336 it.SetCenterPixel(it.GetPixel(indmax));
346 template< class TInputImage, class TOutputImage, class TtNorm >
348 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
349 PrintSelf(std::ostream& os, Indent indent) const
351 Superclass::PrintSelf(os,indent);
353 os << indent << "First Angle: " << m_Alpha1 << std::endl;
354 os << indent << "Second Angle: " << m_Alpha2 << std::endl;
355 os << indent << "K1: " << m_K1 << std::endl;
358 template< class TInputImage, class TOutputImage, class TtNorm >
359 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::TabulationImageType::Pointer
360 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
361 ComputeAngleTabulation(){
363 typename TabulationImageType::Pointer m_AngleTabulation;
364 m_AngleTabulation = TabulationImageType::New();
366 typename TabulationImageType::IndexType start;
368 for(register int i=0;i<ImageDimension;i++)
371 typename TabulationImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
373 for(register int i=0;i<ImageDimension;i++)
376 typename TabulationImageType::RegionType region;
378 region.SetSize(size);
379 region.SetIndex(start);
381 m_AngleTabulation->SetRegions(region);
383 m_AngleTabulation->Allocate();
384 m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion());
385 typedef typename itk::ImageRegionIteratorWithIndex<TabulationImageType>
387 InputIteratorType inputIt = InputIteratorType(m_AngleTabulation, m_AngleTabulation->GetRequestedRegion());
389 typename TabulationImageType::IndexType requestedIndex =
390 m_AngleTabulation->GetRequestedRegion().GetIndex();
392 typename TabulationImageType::SizeType center = this->GetInput()->GetLargestPossibleRegion().GetSize();
393 for(register int i=0;i<ImageDimension;i++)
398 for(inputIt.GoToBegin();!inputIt.IsAtEnd();++inputIt)
400 typename TabulationImageType::IndexType idx = inputIt.GetIndex();
401 typename TabulationImageType::IndexType diff = idx - center;
402 for(int i=0;i<ImageDimension;i++)
404 double tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
406 {//std::cout<<"tem"<<std::endl;
410 double cos_angle = acos((m_DirectionVector*vecttemp)/tmp);
411 inputIt.Set(cos_angle);
415 typedef itk::ImageFileWriter<TabulationImageType> WT;
416 typename WT::Pointer wt = WT::New();
417 wt->SetFileName("testangle.nii");
418 wt->SetInput(m_AngleTabulation);
420 std::cout<<"end compute angle"<<std::endl;
422 return m_AngleTabulation;
425 template< class TInputImage, class TOutputImage, class TtNorm>
426 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::CorrespondanceMapType::Pointer
427 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
428 InitCorrespondanceMap()
430 typename CorrespondanceMapType::Pointer m_CorrespondanceMap;
431 m_CorrespondanceMap = CorrespondanceMapType::New();
432 m_CorrespondanceMap->SetRegions(this->GetInput()->GetLargestPossibleRegion());
433 m_CorrespondanceMap->Allocate();
435 typename InputImageType::IndexType nullIndex;
438 typedef itk::ImageRegionIterator< CorrespondanceMapType > CorrIteratorType;
439 CorrIteratorType it = CorrIteratorType(m_CorrespondanceMap,
440 m_CorrespondanceMap->GetLargestPossibleRegion());
442 typedef itk::ImageRegionConstIteratorWithIndex< InputImageType>
444 InputIteratorType inputIt = InputIteratorType(this->GetInput(),
445 this->GetInput()->GetLargestPossibleRegion());
447 for(it.GoToBegin(),inputIt.GoToBegin();!it.IsAtEnd();++it,++inputIt)
449 it.Set(inputIt.GetIndex());
456 return m_CorrespondanceMap;
467 } // end namespace itk