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 __clitkInvertVFFilter_txx
19 #define __clitkInvertVFFilter_txx
23 //=========================================================================================================================
24 //helper class 1 to allow a threaded execution: add contributions of input to output and update weights
25 //=========================================================================================================================
26 template<class InputImageType, class OutputImageType> class ITK_EXPORT HelperClass1 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
30 /** Standard class typedefs. */
31 typedef HelperClass1 Self;
32 typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
33 typedef itk::SmartPointer<Self> Pointer;
34 typedef itk::SmartPointer<const Self> ConstPointer;
36 /** Method for creation through the object factory. */
39 /** Run-time type information (and related methods) */
40 itkTypeMacro( HelperClass1, ImageToImageFilter );
42 /** Constants for the image dimensions */
43 itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
47 typedef typename OutputImageType::PixelType PixelType;
48 typedef itk::Image<double, ImageDimension > WeightsImageType;
49 typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
51 //===================================================================================
53 void SetWeights(const typename WeightsImageType::Pointer input) {
57 void SetMutexImage(const typename MutexImageType::Pointer input) {
64 typename WeightsImageType::Pointer GetWeights() {
68 /** Typedef to describe the output image region type. */
69 typedef typename OutputImageType::RegionType OutputImageRegionType;
75 //the actual processing
76 void BeforeThreadedGenerateData();
77 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
80 typename WeightsImageType::Pointer m_Weights;
81 typename MutexImageType::Pointer m_MutexImage;
88 //=========================================================================================================================
89 //Member functions of the helper class 1
90 //=========================================================================================================================
93 //=========================================================================================================================
95 template<class InputImageType, class OutputImageType >
96 HelperClass1<InputImageType, OutputImageType>::HelperClass1()
101 //=========================================================================================================================
102 //Before threaded data
103 template<class InputImageType, class OutputImageType >
104 void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
106 //Since we will add, put to zero!
107 this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
108 this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
111 //=========================================================================================================================
112 //update the output for the outputRegionForThread
113 template<class InputImageType, class OutputImageType>
114 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
117 //Get pointer to the input
118 typename InputImageType::ConstPointer inputPtr = this->GetInput();
120 //Get pointer to the output
121 typename OutputImageType::Pointer outputPtr = this->GetOutput();
122 //typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
124 //Iterator over input
125 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
127 //define them over the outputRegionForThread
128 InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
131 typename InputImageType::IndexType index;
132 itk::ContinuousIndex<double,ImageDimension> contIndex;
133 typename InputImageType::PointType ipoint;
134 typename OutputImageType::PointType opoint;
135 typedef typename OutputImageType::PixelType DisplacementType;
136 DisplacementType displacement;
139 //define some temp variables
140 signed long baseIndex[ImageDimension];
141 double distance[ImageDimension];
142 unsigned int dim, counter, upper;
143 double totalOverlap,overlap;
144 typename OutputImageType::IndexType neighIndex;
146 //Find the number of neighbors
147 unsigned int neighbors = 1 << ImageDimension;
149 //==================================================================================================
150 //Loop over the region and add the intensities to the output and the weight to the weights
151 //==================================================================================================
152 while( !inputIt.IsAtEnd() ) {
153 // get the input image index
154 index = inputIt.GetIndex();
155 inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
157 // get the required displacement
158 displacement = inputIt.Get();
160 // compute the required output image point
161 for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
163 // Update the output and the weights
164 if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) {
165 for(dim = 0; dim < ImageDimension; dim++) {
166 // The following block is equivalent to the following line without
167 // having to call floor. (Only for positive inputs, we already now that is in the image)
168 // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
170 baseIndex[dim] = (long) contIndex[dim];
171 distance[dim] = contIndex[dim] - double( baseIndex[dim] );
174 //Add contribution for each neighbor
175 totalOverlap = itk::NumericTraits<double>::Zero;
176 for( counter = 0; counter < neighbors ; counter++ ) {
177 overlap = 1.0; // fraction overlap
178 upper = counter; // each bit indicates upper/lower neighbour
180 // get neighbor index and overlap fraction
181 for( dim = 0; dim < ImageDimension; dim++ ) {
183 neighIndex[dim] = baseIndex[dim] + 1;
184 overlap *= distance[dim];
186 neighIndex[dim] = baseIndex[dim];
187 overlap *= 1.0 - distance[dim];
194 //Set neighbor value only if overlap is not zero
195 if( (overlap>0.0)) // &&
196 // (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
197 // (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
198 // (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
199 // (neighIndex[0]>=0) &&
200 // (neighIndex[1]>=0) &&
201 // (neighIndex[2]>=0) )
203 //what to store? the original displacement vector?
204 if (! m_ThreadSafe) {
205 //Set the pixel and weight at neighIndex
206 outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
207 m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
211 //Entering critilal section: shared memory
212 m_MutexImage->GetPixel(neighIndex).Lock();
214 //Set the pixel and weight at neighIndex
215 outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
216 m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
219 m_MutexImage->GetPixel(neighIndex).Unlock();
222 //Add to total overlap
223 totalOverlap += overlap;
226 if( totalOverlap == 1.0 ) {
240 //=========================================================================================================================
241 //helper class 2 to allow a threaded execution of normalisation by the weights
242 //=========================================================================================================================
243 template<class InputImageType, class OutputImageType> class HelperClass2 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
247 /** Standard class typedefs. */
248 typedef HelperClass2 Self;
249 typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
250 typedef itk::SmartPointer<Self> Pointer;
251 typedef itk::SmartPointer<const Self> ConstPointer;
253 /** Method for creation through the object factory. */
256 /** Run-time type information (and related methods) */
257 itkTypeMacro( HelperClass2, ImageToImageFilter );
259 /** Constants for the image dimensions */
260 itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
263 typedef typename OutputImageType::PixelType PixelType;
264 typedef itk::Image<double,ImageDimension> WeightsImageType;
267 void SetWeights(const typename WeightsImageType::Pointer input) {
271 void SetEdgePaddingValue(PixelType value) {
272 m_EdgePaddingValue = value;
276 /** Typedef to describe the output image region type. */
277 typedef typename OutputImageType::RegionType OutputImageRegionType;
284 //the actual processing
285 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
289 typename WeightsImageType::Pointer m_Weights;
290 PixelType m_EdgePaddingValue;
296 //=========================================================================================================================
297 //Member functions of the helper class 2
298 //=========================================================================================================================
301 //=========================================================================================================================
303 template<class InputImageType, class OutputImageType > HelperClass2<InputImageType, OutputImageType>::HelperClass2()
305 m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
309 //=========================================================================================================================
310 //update the output for the outputRegionForThread
311 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
314 //Get pointer to the input
315 typename InputImageType::ConstPointer inputPtr = this->GetInput();
317 //Get pointer to the output
318 typename OutputImageType::Pointer outputPtr = this->GetOutput();
320 //Iterators over input, weigths and output
321 typedef itk::ImageRegionConstIterator<InputImageType> InputImageIteratorType;
322 typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
323 typedef itk::ImageRegionIterator<WeightsImageType> WeightsImageIteratorType;
325 //define them over the outputRegionForThread
326 OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
327 InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
328 WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread);
331 //==================================================================================================
332 //loop over the output and normalize the input, remove holes
333 PixelType neighValue;
334 double zero = itk::NumericTraits<double>::Zero;
335 while (!outputIt.IsAtEnd()) {
336 //the weight is not zero
337 if (weightsIt.Get() != zero) {
338 //divide by the weight
339 outputIt.Set(static_cast<PixelType>(inputIt.Get()/weightsIt.Get()));
342 //copy the value of the neighbour that was just processed
344 if(!outputIt.IsAtBegin()) {
348 //Neighbour cannot have zero weight because it should be filled already
349 neighValue=outputIt.Get();
351 outputIt.Set(neighValue);
354 //DD("is at begin, setting edgepadding value");
355 outputIt.Set(m_EdgePaddingValue);
366 }//end nameless namespace
373 //=========================================================================================================================
374 // The rest is the InvertVFFilter
375 //=========================================================================================================================
377 //=========================================================================================================================
379 template <class InputImageType, class OutputImageType>
380 InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
382 m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
388 //=========================================================================================================================
390 template <class InputImageType, class OutputImageType> void InvertVFFilter<InputImageType, OutputImageType>::GenerateData()
393 //Get the properties of the input
394 typename InputImageType::ConstPointer inputPtr=this->GetInput();
395 typename WeightsImageType::RegionType region;
396 typename WeightsImageType::RegionType::SizeType size=inputPtr->GetLargestPossibleRegion().GetSize();
397 region.SetSize(size);
398 typename OutputImageType::IndexType start;
399 for (unsigned int i=0; i< ImageDimension; i++) start[i]=0;
400 region.SetIndex(start);
403 //Allocate the weights
404 typename WeightsImageType::Pointer weights=WeightsImageType::New();
405 weights->SetRegions(region);
407 weights->SetSpacing(inputPtr->GetSpacing());
409 //===========================================================================
410 //Inversion is divided in in two loops, for each we will call a threaded helper class
411 //1. add contribution of input to output and update weights
412 //2. normalize the output by the weight and remove holes
413 //===========================================================================
416 //===========================================================================
417 //1. add contribution of input to output and update weights
419 //Define an internal image type
421 typedef itk::Image<itk::Vector<double,ImageDimension>, ImageDimension > InternalImageType;
423 //Call threaded helper class 1
424 typedef HelperClass1<InputImageType, InternalImageType > HelperClass1Type;
425 typename HelperClass1Type::Pointer helper1=HelperClass1Type::New();
428 if(m_NumberOfThreadsIsGiven)helper1->SetNumberOfThreads(m_NumberOfThreads);
429 helper1->SetInput(inputPtr);
430 helper1->SetWeights(weights);
434 //Allocate the mutex image
435 typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New();
436 mutex->SetRegions(region);
438 mutex->SetSpacing(inputPtr->GetSpacing());
439 helper1->SetMutexImage(mutex);
440 if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <<std::endl;
441 } else if(m_Verbose)std::cout <<"Inverting using a thread-unsafe algorithm" <<std::endl;
443 //Execute helper class
447 typename InternalImageType::Pointer temp= helper1->GetOutput();
448 weights=helper1->GetWeights();
451 //===========================================================================
452 //2. Normalize the output by the weights and remove holes
453 //Call threaded helper class
454 typedef HelperClass2<InternalImageType, OutputImageType> HelperClass2Type;
455 typename HelperClass2Type::Pointer helper2=HelperClass2Type::New();
457 //Set temporary output as input
458 helper2->SetInput(temp);
459 helper2->SetWeights(weights);
460 helper2->SetEdgePaddingValue(m_EdgePaddingValue);
462 //Execute helper class
463 if (m_Verbose) std::cout << "Normalizing the output VF..."<<std::endl;
467 this->SetNthOutput(0, helper2->GetOutput());