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 #ifndef __clitkComposeVFFilter_txx
19 #define __clitkComposeVFFilter_txx
20 #include "clitkComposeVFFilter.h"
26 //=========================================================================================================================
28 template <class InputImageType, class OutputImageType>
29 ComposeVFFilter<InputImageType, OutputImageType>::ComposeVFFilter()
31 for (unsigned int i=0; i<ImageDimension; i++) m_EdgePaddingValue[i]=0.0; //no other reasonable value?
36 //=========================================================================================================================
37 //update the output for the outputRegionForThread
38 template<class InputImageType, class OutputImageType>
39 void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
42 //Get pointer to the output
43 typename OutputImageType::Pointer outputPtr = this->GetOutput();
45 //Iterator over output
46 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> Input1ImageIteratorType;
47 typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
49 //define the output and input1 iterator over the outputRegionForThread
50 OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
51 Input1ImageIteratorType input1It(m_Input1, outputRegionForThread);
54 typename InputImageType::IndexType index;
55 itk::ContinuousIndex<double,ImageDimension> contIndex;
56 typename InputImageType::PointType point1;
57 typename InputImageType::PointType point2;
58 typedef typename OutputImageType::PixelType DisplacementType;
59 DisplacementType displacement;
60 DisplacementType totalDisplacement;
63 //define some temp variables outside the loop
64 signed long baseIndex[ImageDimension];
66 double distance[ImageDimension];
72 typename OutputImageType::IndexType neighIndex;
75 //Find the number of neighbors
76 unsigned int neighbors = 1 << ImageDimension;
78 //==================================================================================================
79 //1. loop over the region of the deformationfield and compose the displacements
80 while( !input1It.IsAtEnd() )
82 // get the input image index
83 index = input1It.GetIndex();
84 m_Input1->TransformIndexToPhysicalPoint( index, point1 );
86 // get the displacement of input 1
87 displacement = input1It.Get();
89 // compute the required image point in input2
90 for(unsigned int j = 0; j < ImageDimension; j++ ) point2[j] = point1[j] + (double)displacement[j];
92 //JV TODO add something for the borders
93 //true if inside image
94 if(m_Input2->TransformPhysicalPointToContinuousIndex(point2, contIndex ) )
97 for(dim = 0; dim < ImageDimension; dim++)
99 // The following block is equivalent to the following line without
100 // having to call floor. (Only for positive inputs, we already now that is in the image)
101 // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
103 baseIndex[dim] = (long) contIndex[dim];
104 distance[dim] = contIndex[dim] - double( baseIndex[dim] );
106 //Reset the accumulator
107 totalDisplacement[dim]=0.0;
110 //Add contribution for each neighbor
111 totalOverlap = itk::NumericTraits<double>::Zero;
112 for( counter = 0; counter < neighbors ; counter++ )
114 overlap = 1.0; // fraction overlap
115 upper = counter; // each bit indicates upper/lower neighbour
118 // get neighbor index and overlap fraction
119 bool neighbIndexSupZero = 1;
120 for( dim = 0; dim < ImageDimension; dim++ )
124 neighIndex[dim] = baseIndex[dim] + 1;
125 overlap *= distance[dim];
129 neighIndex[dim] = baseIndex[dim];
130 overlap *= 1.0 - distance[dim];
132 if (neighIndex[dim] < 0)
133 neighbIndexSupZero = 0;
137 //JV shouldn't we verify that the index is not over the upper border instead of zero?
138 // Set neighbor value only if overlap is not zero and index is still in image
139 if ( overlap>0.0 && neighbIndexSupZero )
141 //what to store? the weighted displacement vector of Input2?
142 totalDisplacement+=m_Input2->GetPixel(neighIndex)*overlap;
143 //JV don't we loose more by adding and verifyning each time, it will never be 1.0 before the end no?
144 totalOverlap += overlap;
147 if( totalOverlap == 1.0 )
153 //add the displacement of input1 and the interpolated displacement of input2 to the output
154 outputIt.Set(static_cast<PixelType>(input1It.Get())+totalDisplacement);
156 //displacement from input1 goes outside input2
157 else outputIt.Set(m_EdgePaddingValue);