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 ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
55 // call the superclass' implementation of this method
56 Superclass::GenerateInputRequestedRegion ();
58 // We need all the input.
59 typename InputImageType::Pointer input =
60 const_cast < InputImageType * >(this->GetInput ());
62 input->SetRequestedRegion (input->GetLargestPossibleRegion ());
66 * Generate all the output data
68 template < class TInputImage, class TOutputImage,class TtNorm>
69 void RelativePositionPropImageFilter < TInputImage,
70 TOutputImage ,TtNorm>::EnlargeOutputRequestedRegion (DataObject * output)
72 Superclass::EnlargeOutputRequestedRegion (output);
73 output->SetRequestedRegionToLargestPossibleRegion ();
77 inline float min(float a,float b)
82 template< class TInputImage, class TOutputImage ,class TtNorm>
84 RelativePositionPropImageFilter< TInputImage, TOutputImage ,TtNorm>
85 //::GenerateThreadedData(const typename TOutputImage::RegionType& outputRegionForThread, int threadId)
89 // std::cout <<"GenerateData" << std::endl;
91 this->AllocateOutputs();
93 typename InputImageType::ConstPointer input = this->GetInput();
95 typename CorrespondanceMapType::Pointer m_CorrespondanceMap = InitCorrespondanceMap();
97 typename InputImageType::IndexType nullIndex;
100 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputConstIteratorWithIndexType;
101 typedef itk::ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
103 typedef itk::NeighborhoodIterator< CorrespondanceMapType > NeighborhoodIteratorType;
106 typename NeighborhoodIteratorType::RadiusType radius;
107 radius.Fill(m_Radius);
108 NeighborhoodIteratorType it( radius,
110 m_CorrespondanceMap->GetLargestPossibleRegion());
111 OutputIteratorType outputIt( this->GetOutput(),
112 this->GetOutput()->GetLargestPossibleRegion() );
114 typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
115 unsigned int totalSize = 1;
116 for(int i=0;i<ImageDimension;i++)
118 totalSize *= size[i];
120 const typename InputImageType::OffsetValueType * offsetTable = input->GetOffsetTable();
128 typename InputImageType::IndexType ind1;
130 //unsigned int centeroffset = (int)(it.Size() / 2);
135 for(int i=0;i<pow((double)2,(int)ImageDimension);i++)
137 int orient[ImageDimension];
138 int init[ImageDimension];
139 typename InputImageType::OffsetType tempOffset;
143 outputIt.GoToBegin();
144 for(int j=0;j<ImageDimension;j++)
146 orient[j] = ((int)(i/pow((double)2,j)))%2;
147 tempOffset[j] = (orient[j])*(size[j]-1);
151 outputIt.SetIndex(it.GetIndex());
153 for(int j=0;j<ImageDimension;j++)
156 init[j] = -size[j]+1;
162 for(unsigned int p=1;p<=totalSize;p++)
165 if ((m_VerboseProgress && (p % totalSize/100) == 0)) {
170 for (unsigned i = 0; i < it.Size(); i++)
172 if( it.GetPixel(i,in)!=nullIndex && in )
174 for(int j=0;j<ImageDimension;j++)
176 ind1[j] = it.GetPixel(i)[j];
177 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
180 tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
183 if((m_DirectionVector*vecttemp)/tmp>=1)
186 if((m_DirectionVector*vecttemp)/tmp<=-1)
189 tmp = acos((m_DirectionVector*vecttemp)/tmp);
191 tmp=std::max(0.,1.-(tmp/m_K1));
192 // tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
194 tmp = m_TNorm( input->GetPixel(ind1), tmp);
205 it.SetCenterPixel(it.GetPixel(indmax));
208 else outputIt.Set( 0 );
216 for(int j=0;j<ImageDimension;j++)
218 if(p%offsetTable[j]==0)
224 tempOffset[k] = init[k];
230 tempOffset[k] = init[k];
235 outputIt.SetIndex(it.GetIndex());
240 // ==================================================================================================
241 else //if fast compute in two pass
243 // std::cout<<"pass 1"<<std::endl;
249 long nb = input->GetLargestPossibleRegion().GetNumberOfPixels( );
250 for (it.GoToBegin(), outputIt.GoToBegin(); ! it.IsAtEnd(); ++it, ++outputIt)
252 if (m_VerboseProgress && (i%(nb/10)) == 0) {
254 std::cout << i << " / " << nb << std::endl;
260 for (unsigned i = 0; i < it.Size(); i++) {
261 if( it.GetPixel(i,in)!=nullIndex && in )
264 // DS : can be precomputed ??? no
265 for(int j=0;j<ImageDimension;j++)
266 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
269 tmp = vecttemp.GetNorm();
271 tmp = acos((m_DirectionVector*vecttemp)/tmp);
272 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
273 tmp=std::max(0.,1.-tmp/m_K1);
274 tmp = min( input->GetPixel(it.GetPixel(i)), tmp);
285 it.SetCenterPixel(it.GetPixel(indmax));
293 // std::cout<<"pass 2"<<std::endl;
299 nb = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels( );
301 for ( outputIt.GoToReverseBegin(); ! outputIt.IsAtReverseEnd(); --it, --outputIt)
303 if (m_VerboseProgress && (i%(nb/10)) == 0) {
305 std::cout << i << " / " << nb << std::endl;
312 for (unsigned i = 0; i < it.Size(); i++)
313 if( it.GetPixel(i,in)!=nullIndex && in )
315 for(int j=0;j<ImageDimension;j++)
316 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
318 tmp = vecttemp.GetNorm();
321 tmp = acos((m_DirectionVector*vecttemp)/tmp);
322 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
323 tmp=std::max(0.,1. - tmp/m_K1);
324 tmp = min( input->GetPixel(it.GetPixel(i)),tmp);
335 it.SetCenterPixel(it.GetPixel(indmax));
345 template< class TInputImage, class TOutputImage, class TtNorm >
347 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
348 PrintSelf(std::ostream& os, Indent indent) const
350 Superclass::PrintSelf(os,indent);
352 os << indent << "First Angle: " << m_Alpha1 << std::endl;
353 os << indent << "Second Angle: " << m_Alpha2 << std::endl;
354 os << indent << "K1: " << m_K1 << std::endl;
357 template< class TInputImage, class TOutputImage, class TtNorm >
358 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::TabulationImageType::Pointer
359 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
360 ComputeAngleTabulation(){
362 typename TabulationImageType::Pointer m_AngleTabulation;
363 m_AngleTabulation = TabulationImageType::New();
365 typename TabulationImageType::IndexType start;
367 for(register int i=0;i<ImageDimension;i++)
370 typename TabulationImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
372 for(register int i=0;i<ImageDimension;i++)
375 typename TabulationImageType::RegionType region;
377 region.SetSize(size);
378 region.SetIndex(start);
380 m_AngleTabulation->SetRegions(region);
382 m_AngleTabulation->Allocate();
383 m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion());
384 typedef typename itk::ImageRegionIteratorWithIndex<TabulationImageType>
386 InputIteratorType inputIt = InputIteratorType(m_AngleTabulation, m_AngleTabulation->GetRequestedRegion());
388 typename TabulationImageType::IndexType requestedIndex =
389 m_AngleTabulation->GetRequestedRegion().GetIndex();
391 typename TabulationImageType::SizeType center = this->GetInput()->GetLargestPossibleRegion().GetSize();
392 for(register int i=0;i<ImageDimension;i++)
397 for(inputIt.GoToBegin();!inputIt.IsAtEnd();++inputIt)
399 typename TabulationImageType::IndexType idx = inputIt.GetIndex();
400 typename TabulationImageType::IndexType diff = idx - center;
401 for(int i=0;i<ImageDimension;i++)
403 double tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
405 {//std::cout<<"tem"<<std::endl;
409 double cos_angle = acos((m_DirectionVector*vecttemp)/tmp);
410 inputIt.Set(cos_angle);
414 typedef itk::ImageFileWriter<TabulationImageType> WT;
415 typename WT::Pointer wt = WT::New();
416 wt->SetFileName("testangle.nii");
417 wt->SetInput(m_AngleTabulation);
419 std::cout<<"end compute angle"<<std::endl;
421 return m_AngleTabulation;
424 template< class TInputImage, class TOutputImage, class TtNorm>
425 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::CorrespondanceMapType::Pointer
426 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
427 InitCorrespondanceMap()
429 typename CorrespondanceMapType::Pointer m_CorrespondanceMap;
430 m_CorrespondanceMap = CorrespondanceMapType::New();
431 m_CorrespondanceMap->SetRegions(this->GetInput()->GetLargestPossibleRegion());
432 m_CorrespondanceMap->Allocate();
434 typename InputImageType::IndexType nullIndex;
437 typedef itk::ImageRegionIterator< CorrespondanceMapType > CorrIteratorType;
438 CorrIteratorType it = CorrIteratorType(m_CorrespondanceMap,
439 m_CorrespondanceMap->GetLargestPossibleRegion());
441 typedef itk::ImageRegionConstIteratorWithIndex< InputImageType>
443 InputIteratorType inputIt = InputIteratorType(this->GetInput(),
444 this->GetInput()->GetLargestPossibleRegion());
446 for(it.GoToBegin(),inputIt.GoToBegin();!it.IsAtEnd();++it,++inputIt)
448 it.Set(inputIt.GetIndex());
455 return m_CorrespondanceMap;
466 } // end namespace itk