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 this->AllocateOutputs();
91 typename InputImageType::ConstPointer input = this->GetInput();
93 typename CorrespondanceMapType::Pointer m_CorrespondanceMap = InitCorrespondanceMap();
95 typename InputImageType::IndexType nullIndex;
98 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputConstIteratorWithIndexType;
99 typedef itk::ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
101 typedef itk::NeighborhoodIterator< CorrespondanceMapType > NeighborhoodIteratorType;
104 typename NeighborhoodIteratorType::RadiusType radius;
105 radius.Fill(m_Radius);
106 NeighborhoodIteratorType it( radius,
108 m_CorrespondanceMap->GetLargestPossibleRegion());
109 OutputIteratorType outputIt( this->GetOutput(),
110 this->GetOutput()->GetLargestPossibleRegion() );
112 typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
113 unsigned int totalSize = 1;
114 for(int i=0;i<ImageDimension;i++)
116 totalSize *= size[i];
118 const typename InputImageType::OffsetValueType * offsetTable = input->GetOffsetTable();
126 typename InputImageType::IndexType ind1;
128 //unsigned int centeroffset = (int)(it.Size() / 2);
132 for(int i=0;i<pow((double)2,(int)ImageDimension);i++)
134 int orient[ImageDimension];
135 int init[ImageDimension];
136 typename InputImageType::OffsetType tempOffset;
140 outputIt.GoToBegin();
141 for(int j=0;j<ImageDimension;j++)
143 orient[j] = ((int)(i/pow((double)2,j)))%2;
144 tempOffset[j] = (orient[j])*(size[j]-1);
148 outputIt.SetIndex(it.GetIndex());
150 for(int j=0;j<ImageDimension;j++)
153 init[j] = -size[j]+1;
159 for(unsigned int p=1;p<=totalSize;p++)
162 if ((m_VerboseProgress && (p % totalSize/100) == 0)) {
167 for (unsigned i = 0; i < it.Size(); i++)
169 if( it.GetPixel(i,in)!=nullIndex && in )
171 for(int j=0;j<ImageDimension;j++)
173 ind1[j] = it.GetPixel(i)[j];
174 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
177 tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
180 if((m_DirectionVector*vecttemp)/tmp>=1)
183 if((m_DirectionVector*vecttemp)/tmp<=-1)
186 tmp = acos((m_DirectionVector*vecttemp)/tmp);
188 tmp=std::max(0.,1.-(tmp/m_K1));
189 // tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
191 tmp = m_TNorm( input->GetPixel(ind1), tmp);
202 it.SetCenterPixel(it.GetPixel(indmax));
205 else outputIt.Set( 0 );
213 for(int j=0;j<ImageDimension;j++)
215 if(p%offsetTable[j]==0)
221 tempOffset[k] = init[k];
227 tempOffset[k] = init[k];
232 outputIt.SetIndex(it.GetIndex());
237 // ==================================================================================================
238 else //if fast compute in two pass
240 // std::cout<<"pass 1"<<std::endl;
246 long nb = input->GetLargestPossibleRegion().GetNumberOfPixels( );
247 for (it.GoToBegin(), outputIt.GoToBegin(); ! it.IsAtEnd(); ++it, ++outputIt)
249 if (m_VerboseProgress && (i%(nb/10)) == 0) {
251 std::cout << i << " / " << nb << std::endl;
257 for (unsigned i = 0; i < it.Size(); i++) {
258 if( it.GetPixel(i,in)!=nullIndex && in )
261 // DS : can be precomputed ??? no
262 for(int j=0;j<ImageDimension;j++)
263 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
266 tmp = vecttemp.GetNorm();
268 tmp = acos((m_DirectionVector*vecttemp)/tmp);
269 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
270 tmp=std::max(0.,1.-tmp/m_K1);
271 tmp = min( input->GetPixel(it.GetPixel(i)), tmp);
282 it.SetCenterPixel(it.GetPixel(indmax));
290 // std::cout<<"pass 2"<<std::endl;
296 nb = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels( );
298 for ( outputIt.GoToReverseBegin(); ! outputIt.IsAtReverseEnd(); --it, --outputIt)
300 if (m_VerboseProgress && (i%(nb/10)) == 0) {
302 std::cout << i << " / " << nb << std::endl;
309 for (unsigned i = 0; i < it.Size(); i++)
310 if( it.GetPixel(i,in)!=nullIndex && in )
312 for(int j=0;j<ImageDimension;j++)
313 vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
315 tmp = vecttemp.GetNorm();
318 tmp = acos((m_DirectionVector*vecttemp)/tmp);
319 //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
320 tmp=std::max(0.,1. - tmp/m_K1);
321 tmp = min( input->GetPixel(it.GetPixel(i)),tmp);
332 it.SetCenterPixel(it.GetPixel(indmax));
342 template< class TInputImage, class TOutputImage, class TtNorm >
344 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
345 PrintSelf(std::ostream& os, Indent indent) const
347 Superclass::PrintSelf(os,indent);
349 os << indent << "First Angle: " << m_Alpha1 << std::endl;
350 os << indent << "Second Angle: " << m_Alpha2 << std::endl;
351 os << indent << "K1: " << m_K1 << std::endl;
354 template< class TInputImage, class TOutputImage, class TtNorm >
355 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::TabulationImageType::Pointer
356 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
357 ComputeAngleTabulation(){
359 typename TabulationImageType::Pointer m_AngleTabulation;
360 m_AngleTabulation = TabulationImageType::New();
362 typename TabulationImageType::IndexType start;
364 for(register int i=0;i<ImageDimension;i++)
367 typename TabulationImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
369 for(register int i=0;i<ImageDimension;i++)
372 typename TabulationImageType::RegionType region;
374 region.SetSize(size);
375 region.SetIndex(start);
377 m_AngleTabulation->SetRegions(region);
379 m_AngleTabulation->Allocate();
380 m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion());
381 typedef typename itk::ImageRegionIteratorWithIndex<TabulationImageType>
383 InputIteratorType inputIt = InputIteratorType(m_AngleTabulation, m_AngleTabulation->GetRequestedRegion());
385 typename TabulationImageType::IndexType requestedIndex =
386 m_AngleTabulation->GetRequestedRegion().GetIndex();
388 typename TabulationImageType::SizeType center = this->GetInput()->GetLargestPossibleRegion().GetSize();
389 for(register int i=0;i<ImageDimension;i++)
394 for(inputIt.GoToBegin();!inputIt.IsAtEnd();++inputIt)
396 typename TabulationImageType::IndexType idx = inputIt.GetIndex();
397 typename TabulationImageType::IndexType diff = idx - center;
398 for(int i=0;i<ImageDimension;i++)
400 double tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
402 {//std::cout<<"tem"<<std::endl;
406 double cos_angle = acos((m_DirectionVector*vecttemp)/tmp);
407 inputIt.Set(cos_angle);
411 typedef itk::ImageFileWriter<TabulationImageType> WT;
412 typename WT::Pointer wt = WT::New();
413 wt->SetFileName("testangle.nii");
414 wt->SetInput(m_AngleTabulation);
416 std::cout<<"end compute angle"<<std::endl;
418 return m_AngleTabulation;
421 template< class TInputImage, class TOutputImage, class TtNorm>
422 typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::CorrespondanceMapType::Pointer
423 RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
424 InitCorrespondanceMap()
426 typename CorrespondanceMapType::Pointer m_CorrespondanceMap;
427 m_CorrespondanceMap = CorrespondanceMapType::New();
428 m_CorrespondanceMap->SetRegions(this->GetInput()->GetLargestPossibleRegion());
429 m_CorrespondanceMap->Allocate();
431 typename InputImageType::IndexType nullIndex;
434 typedef itk::ImageRegionIterator< CorrespondanceMapType > CorrIteratorType;
435 CorrIteratorType it = CorrIteratorType(m_CorrespondanceMap,
436 m_CorrespondanceMap->GetLargestPossibleRegion());
438 typedef itk::ImageRegionConstIteratorWithIndex< InputImageType>
440 InputIteratorType inputIt = InputIteratorType(this->GetInput(),
441 this->GetInput()->GetLargestPossibleRegion());
443 for(it.GoToBegin(),inputIt.GoToBegin();!it.IsAtEnd();++it,++inputIt)
445 it.Set(inputIt.GetIndex());
452 return m_CorrespondanceMap;
463 } // end namespace itk