]> Creatis software - clitk.git/blob - itk/clitkInvertVFFilter.txx
invert Vf on given image grid with resampler
[clitk.git] / itk / clitkInvertVFFilter.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
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
8
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.
12
13   It is distributed under dual licence
14
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
20
21 namespace
22 {
23
24 //=========================================================================================================================
25 //helper class 1 to allow a threaded execution: add contributions of input to output and update weights
26 //=========================================================================================================================
27 template<class InputImageType, class OutputImageType> class ITK_EXPORT HelperClass1 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
28 {
29
30 public:
31   /** Standard class typedefs. */
32   typedef HelperClass1  Self;
33   typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
34   typedef itk::SmartPointer<Self>         Pointer;
35   typedef itk::SmartPointer<const Self>   ConstPointer;
36
37   /** Method for creation through the object factory. */
38   itkNewMacro(Self);
39
40   /** Run-time type information (and related methods) */
41   itkTypeMacro( HelperClass1, ImageToImageFilter );
42
43   /** Constants for the image dimensions */
44   itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
45
46
47   //Typedefs
48   typedef typename OutputImageType::PixelType        PixelType;
49   typedef itk::Image<double, ImageDimension > WeightsImageType;
50   typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
51
52   //===================================================================================
53   //Set methods
54   void SetWeights(const typename WeightsImageType::Pointer input) {
55     m_Weights = input;
56     this->Modified();
57   }
58   void SetMutexImage(const typename MutexImageType::Pointer input) {
59     m_MutexImage=input;
60     this->Modified();
61     m_ThreadSafe=true;
62   }
63
64   //Get methods
65   typename  WeightsImageType::Pointer GetWeights() {
66     return m_Weights;
67   }
68
69   /** Typedef to describe the output image region type. */
70   typedef typename OutputImageType::RegionType OutputImageRegionType;
71
72 protected:
73   HelperClass1();
74   ~HelperClass1() {};
75
76   //the actual processing
77   void BeforeThreadedGenerateData();
78   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
79
80   //member data
81   typename  WeightsImageType::Pointer m_Weights;
82   typename MutexImageType::Pointer m_MutexImage;
83   bool m_ThreadSafe;
84
85 };
86
87
88
89 //=========================================================================================================================
90 //Member functions of the helper class 1
91 //=========================================================================================================================
92
93
94 //=========================================================================================================================
95 //Empty constructor
96 template<class InputImageType, class OutputImageType >
97 HelperClass1<InputImageType, OutputImageType>::HelperClass1()
98 {
99   m_ThreadSafe=false;
100 }
101
102 //=========================================================================================================================
103 //Before threaded data
104 template<class InputImageType, class OutputImageType >
105 void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
106 {
107   //std::cout << "HelperClass1::BeforeThreadedGenerateData - IN" << std::endl;
108   //Since we will add, put to zero!
109   this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
110   this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
111 }
112
113 //=========================================================================================================================
114 //update the output for the outputRegionForThread
115 template<class InputImageType, class OutputImageType>
116 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
117 {
118   //std::cout << "HelperClass1::ThreadedGenerateData - IN" << std::endl;
119   //Get pointer to the input
120   typename InputImageType::ConstPointer inputPtr = this->GetInput();
121
122   //Get pointer to the output
123   typename OutputImageType::Pointer outputPtr = this->GetOutput();
124   //typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
125
126   //Iterator over input
127   typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
128
129   //define them over the outputRegionForThread
130   InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
131
132   //Initialize
133   typename InputImageType::IndexType index;
134   itk::ContinuousIndex<double,ImageDimension> contIndex, inContIndex;
135   typename InputImageType::PointType ipoint;
136   typename OutputImageType::PointType opoint;
137   typedef typename OutputImageType::PixelType DisplacementType;
138   DisplacementType displacement;
139   inputIt.GoToBegin();
140
141   //define some temp variables
142   signed long baseIndex[ImageDimension];
143   double distance[ImageDimension];
144   unsigned int dim, counter, upper;
145   double totalOverlap,overlap;
146   typename OutputImageType::IndexType neighIndex;
147
148   //Find the number of neighbors
149   unsigned int neighbors =  1 << ImageDimension;
150
151   //==================================================================================================
152   //Loop over the output region and add the intensities from the input to the output and the weight to the weights
153   //==================================================================================================
154   while( !inputIt.IsAtEnd() ) {
155     // get the input image index
156     index = inputIt.GetIndex();
157     inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
158
159     // get the required displacement
160     displacement = inputIt.Get();
161
162     // compute the required output image point
163     for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
164
165     // Update the output and the weights
166     if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) {
167       for(dim = 0; dim < ImageDimension; dim++) {
168         // The following  block is equivalent to the following line without
169         // having to call floor. (Only for positive inputs, we already now that is in the image)
170         // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
171
172         baseIndex[dim] = (long) contIndex[dim];
173         distance[dim] = contIndex[dim] - double( baseIndex[dim] );
174       }
175
176       //Add contribution for each neighbor
177       totalOverlap = itk::NumericTraits<double>::Zero;
178       for( counter = 0; counter < neighbors ; counter++ ) {
179         overlap = 1.0;          // fraction overlap
180         upper = counter;  // each bit indicates upper/lower neighbour
181
182         // get neighbor index and overlap fraction
183         for( dim = 0; dim < ImageDimension; dim++ ) {
184           if ( upper & 1 ) {
185             neighIndex[dim] = baseIndex[dim] + 1;
186             overlap *= distance[dim];
187           } else {
188             neighIndex[dim] = baseIndex[dim];
189             overlap *= 1.0 - distance[dim];
190           }
191           upper >>= 1;
192         }
193
194
195
196         //Set neighbor value only if overlap is not zero
197         if( (overlap>0.0)) // &&
198           //                    (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
199           //                    (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
200           //                    (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
201           //                    (neighIndex[0]>=0) &&
202           //                    (neighIndex[1]>=0) &&
203           //                    (neighIndex[2]>=0) )
204         {
205           //what to store? the original displacement vector?
206           if (! m_ThreadSafe) {
207             //Set the pixel and weight at neighIndex
208             outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
209             m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
210           }
211
212           else {
213             //Entering critilal section: shared memory
214             m_MutexImage->GetPixel(neighIndex).Lock();
215
216             //Set the pixel and weight at neighIndex
217             outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
218             m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
219
220             //Unlock
221             m_MutexImage->GetPixel(neighIndex).Unlock();
222
223           }
224           //Add to total overlap
225           totalOverlap += overlap;
226         }
227
228         if( totalOverlap == 1.0 ) {
229           // finished
230           break;
231         }
232       }
233     }
234
235     ++inputIt;
236   }
237
238   //std::cout << "HelperClass1::ThreadedGenerateData - OUT" << std::endl;
239 }
240
241
242
243 //=========================================================================================================================
244 //helper class 2 to allow a threaded execution of normalisation by the weights
245 //=========================================================================================================================
246 template<class InputImageType, class OutputImageType> class HelperClass2 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
247 {
248
249 public:
250   /** Standard class typedefs. */
251   typedef HelperClass2  Self;
252   typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
253   typedef itk::SmartPointer<Self>         Pointer;
254   typedef itk::SmartPointer<const Self>   ConstPointer;
255
256   /** Method for creation through the object factory. */
257   itkNewMacro(Self);
258
259   /** Run-time type information (and related methods) */
260   itkTypeMacro( HelperClass2, ImageToImageFilter );
261
262   /** Constants for the image dimensions */
263   itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
264
265   //Typedefs
266   typedef typename OutputImageType::PixelType        PixelType;
267   typedef itk::Image<double,ImageDimension> WeightsImageType;
268
269   //Set methods
270   void SetWeights(const typename WeightsImageType::Pointer input) {
271     m_Weights = input;
272     this->Modified();
273   }
274   void SetEdgePaddingValue(PixelType value) {
275     m_EdgePaddingValue = value;
276     this->Modified();
277   }
278
279   /** Typedef to describe the output image region type. */
280   typedef typename OutputImageType::RegionType OutputImageRegionType;
281
282 protected:
283   HelperClass2();
284   ~HelperClass2() {};
285
286
287   //the actual processing
288   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
289
290
291   //member data
292   typename     WeightsImageType::Pointer m_Weights;
293   PixelType m_EdgePaddingValue;
294
295 } ;
296
297
298
299 //=========================================================================================================================
300 //Member functions of the helper class 2
301 //=========================================================================================================================
302
303
304 //=========================================================================================================================
305 //Empty constructor
306 template<class InputImageType, class OutputImageType > HelperClass2<InputImageType, OutputImageType>::HelperClass2()
307 {
308   m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
309 }
310
311
312 //=========================================================================================================================
313 //update the output for the outputRegionForThread
314 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
315 {
316   //std::cout << "HelperClass2::ThreadedGenerateData - IN" << std::endl;
317   
318   //Get pointer to the input
319   typename InputImageType::ConstPointer inputPtr = this->GetInput();
320
321   //Get pointer to the output
322   typename OutputImageType::Pointer outputPtr = this->GetOutput();
323
324   //Iterators over input, weigths  and output
325   typedef itk::ImageRegionConstIterator<InputImageType> InputImageIteratorType;
326   typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
327   typedef itk::ImageRegionIterator<WeightsImageType> WeightsImageIteratorType;
328
329   //define them over the outputRegionForThread
330   OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
331   InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
332   WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread);
333
334
335   //==================================================================================================
336   //loop over the output and normalize the input, remove holes
337   PixelType neighValue;
338   double  zero = itk::NumericTraits<double>::Zero;
339   while (!outputIt.IsAtEnd()) {
340     //the weight is not zero
341     if (weightsIt.Get() != zero) {
342       //divide by the weight
343       outputIt.Set(static_cast<PixelType>(inputIt.Get()/weightsIt.Get()));
344     }
345
346     //copy the value of the  neighbour that was just processed
347     else {
348       if(!outputIt.IsAtBegin()) {
349         //go back
350         --outputIt;
351
352         //Neighbour cannot have zero weight because it should be filled already
353         neighValue=outputIt.Get();
354         ++outputIt;
355         outputIt.Set(neighValue);
356         //DD("hole filled");
357       } else {
358         //DD("is at begin, setting edgepadding value");
359         outputIt.Set(m_EdgePaddingValue);
360       }
361     }
362     ++weightsIt;
363     ++outputIt;
364     ++inputIt;
365
366   }//end while
367   
368   //std::cout << "HelperClass2::ThreadedGenerateData - OUT" << std::endl;
369   
370 }//end member
371
372
373 }//end nameless namespace
374
375
376
377 namespace clitk
378 {
379
380 //=========================================================================================================================
381 // The rest is the InvertVFFilter
382 //=========================================================================================================================
383
384 //=========================================================================================================================
385 //constructor
386 template <class InputImageType, class OutputImageType>
387 InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
388 {
389   m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
390   m_ThreadSafe=false;
391   m_Verbose=false;
392 }
393
394 //=========================================================================================================================
395 //Update
396 template <class InputImageType, class OutputImageType> void InvertVFFilter<InputImageType, OutputImageType>::GenerateData()
397 {
398   //std::cout << "InvertVFFilter::GenerateData - IN" << std::endl;
399
400   //Get the properties of the input
401   typename InputImageType::ConstPointer inputPtr=this->GetInput();
402   typename WeightsImageType::RegionType region = inputPtr->GetLargestPossibleRegion();
403
404   //Allocate the weights
405   typename WeightsImageType::Pointer weights=WeightsImageType::New();
406   weights->SetOrigin(inputPtr->GetOrigin());
407   weights->SetRegions(region);
408   weights->Allocate();
409   weights->SetSpacing(inputPtr->GetSpacing());
410
411   //===========================================================================
412   //Inversion is divided in in two loops, for each we will call a threaded helper class
413   //1. add contribution of input to output and update weights
414   //2. normalize the output by the weight and remove holes
415   //===========================================================================
416
417
418   //===========================================================================
419   //1. add contribution of input to output and update weights
420
421   //Define an internal image type
422
423   typedef itk::Image<itk::Vector<double,ImageDimension>, ImageDimension > InternalImageType;
424
425   //Call threaded helper class 1
426   typedef HelperClass1<InputImageType, InternalImageType > HelperClass1Type;
427   typename HelperClass1Type::Pointer helper1=HelperClass1Type::New();
428
429   //Set input
430   if(m_NumberOfThreadsIsGiven)helper1->SetNumberOfThreads(m_NumberOfThreads);
431   helper1->SetInput(inputPtr);
432   helper1->SetWeights(weights);
433
434   //Threadsafe?
435   if(m_ThreadSafe) {
436     //Allocate the mutex image
437     typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New();
438     mutex->SetRegions(region);
439     mutex->Allocate();
440     mutex->SetSpacing(inputPtr->GetSpacing());
441     helper1->SetMutexImage(mutex);
442     if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <<std::endl;
443   } else  if(m_Verbose)std::cout <<"Inverting using a thread-unsafe algorithm" <<std::endl;
444
445   //Execute helper class
446   helper1->Update();
447
448   //Get the output
449   typename InternalImageType::Pointer temp= helper1->GetOutput();
450   weights=helper1->GetWeights();
451
452
453   //===========================================================================
454   //2. Normalize the output by the weights and remove holes
455   //Call threaded helper class
456   typedef HelperClass2<InternalImageType, OutputImageType> HelperClass2Type;
457   typename HelperClass2Type::Pointer helper2=HelperClass2Type::New();
458
459   //Set temporary output as input
460   helper2->SetInput(temp);
461   helper2->SetWeights(weights);
462   helper2->SetEdgePaddingValue(m_EdgePaddingValue);
463
464   //Execute helper class
465   if (m_Verbose) std::cout << "Normalizing the output VF..."<<std::endl;
466   helper2->Update();
467
468   //Set the output
469   this->SetNthOutput(0, helper2->GetOutput());
470 }
471
472
473
474 }
475
476 #endif