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 #if ITK_VERSION_MAJOR >= 4
39 template<class InputImageType, class OutputImageType>
40 void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
42 template<class InputImageType, class OutputImageType>
43 void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
47 //Get pointer to the output
48 typename OutputImageType::Pointer outputPtr = this->GetOutput();
50 //Iterator over output
51 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> Input1ImageIteratorType;
52 typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
54 //define the output and input1 iterator over the outputRegionForThread
55 OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
56 Input1ImageIteratorType input1It(m_Input1, outputRegionForThread);
59 typename InputImageType::IndexType index;
60 itk::ContinuousIndex<double,ImageDimension> contIndex;
61 typename InputImageType::PointType point1;
62 typename InputImageType::PointType point2;
63 typedef typename OutputImageType::PixelType DisplacementType;
64 DisplacementType displacement;
65 DisplacementType totalDisplacement;
68 //define some temp variables outside the loop
69 signed long baseIndex[ImageDimension];
71 double distance[ImageDimension];
77 typename OutputImageType::IndexType neighIndex;
80 //Find the number of neighbors
81 unsigned int neighbors = 1 << ImageDimension;
83 //==================================================================================================
84 //1. loop over the region of the deformationfield and compose the displacements
85 while( !input1It.IsAtEnd() )
87 // get the input image index
88 index = input1It.GetIndex();
89 m_Input1->TransformIndexToPhysicalPoint( index, point1 );
91 // get the displacement of input 1
92 displacement = input1It.Get();
94 // compute the required image point in input2
95 for(unsigned int j = 0; j < ImageDimension; j++ ) point2[j] = point1[j] + (double)displacement[j];
97 //JV TODO add something for the borders
98 //true if inside image
99 if(m_Input2->TransformPhysicalPointToContinuousIndex(point2, contIndex ) )
102 for(dim = 0; dim < ImageDimension; dim++)
104 // The following block is equivalent to the following line without
105 // having to call floor. (Only for positive inputs, we already now that is in the image)
106 // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
108 baseIndex[dim] = (long) contIndex[dim];
109 distance[dim] = contIndex[dim] - double( baseIndex[dim] );
111 //Reset the accumulator
112 totalDisplacement[dim]=0.0;
115 //Add contribution for each neighbor
116 totalOverlap = itk::NumericTraits<double>::Zero;
117 for( counter = 0; counter < neighbors ; counter++ )
119 overlap = 1.0; // fraction overlap
120 upper = counter; // each bit indicates upper/lower neighbour
123 // get neighbor index and overlap fraction
124 bool neighbIndexSupZero = 1;
125 for( dim = 0; dim < ImageDimension; dim++ )
129 neighIndex[dim] = baseIndex[dim] + 1;
130 overlap *= distance[dim];
134 neighIndex[dim] = baseIndex[dim];
135 overlap *= 1.0 - distance[dim];
137 if (neighIndex[dim] < 0)
138 neighbIndexSupZero = 0;
142 //JV shouldn't we verify that the index is not over the upper border instead of zero?
143 // Set neighbor value only if overlap is not zero and index is still in image
144 if ( overlap>0.0 && neighbIndexSupZero )
146 //what to store? the weighted displacement vector of Input2?
147 totalDisplacement+=m_Input2->GetPixel(neighIndex)*overlap;
148 //JV don't we loose more by adding and verifyning each time, it will never be 1.0 before the end no?
149 totalOverlap += overlap;
152 if( totalOverlap == 1.0 )
158 //add the displacement of input1 and the interpolated displacement of input2 to the output
159 outputIt.Set(static_cast<PixelType>(input1It.Get())+totalDisplacement);
161 //displacement from input1 goes outside input2
162 else outputIt.Set(m_EdgePaddingValue);