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
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>
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;
37 /** Method for creation through the object factory. */
40 /** Run-time type information (and related methods) */
41 itkTypeMacro( HelperClass1, ImageToImageFilter );
43 /** Constants for the image dimensions */
44 itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
48 typedef typename OutputImageType::PixelType PixelType;
49 typedef itk::Image<double, ImageDimension > WeightsImageType;
50 typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
52 //===================================================================================
54 void SetWeights(const typename WeightsImageType::Pointer input) {
58 void SetMutexImage(const typename MutexImageType::Pointer input) {
65 typename WeightsImageType::Pointer GetWeights() {
69 /** Typedef to describe the output image region type. */
70 typedef typename OutputImageType::RegionType OutputImageRegionType;
76 //the actual processing
77 void BeforeThreadedGenerateData();
78 #if ITK_VERSION_MAJOR >= 4
79 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId );
81 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
85 typename WeightsImageType::Pointer m_Weights;
86 typename MutexImageType::Pointer m_MutexImage;
93 //=========================================================================================================================
94 //Member functions of the helper class 1
95 //=========================================================================================================================
98 //=========================================================================================================================
100 template<class InputImageType, class OutputImageType >
101 HelperClass1<InputImageType, OutputImageType>::HelperClass1()
106 //=========================================================================================================================
107 //Before threaded data
108 template<class InputImageType, class OutputImageType >
109 void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
111 //std::cout << "HelperClass1::BeforeThreadedGenerateData - IN" << std::endl;
112 //Since we will add, put to zero!
113 this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
114 this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
117 //=========================================================================================================================
118 //update the output for the outputRegionForThread
119 template<class InputImageType, class OutputImageType>
120 #if ITK_VERSION_MAJOR >= 4
121 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
123 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
126 // std::cout << "HelperClass1::ThreadedGenerateData - IN " << threadId << std::endl;
127 //Get pointer to the input
128 typename InputImageType::ConstPointer inputPtr = this->GetInput();
130 //Get pointer to the output
131 typename OutputImageType::Pointer outputPtr = this->GetOutput();
132 //typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
134 //Iterator over input
135 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
137 //define them over the outputRegionForThread
138 InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
141 typename InputImageType::IndexType index;
142 itk::ContinuousIndex<double,ImageDimension> contIndex, inContIndex;
143 typename InputImageType::PointType ipoint;
144 typename OutputImageType::PointType opoint;
145 typedef typename OutputImageType::PixelType DisplacementType;
146 DisplacementType displacement;
149 typename OutputImageType::SizeType size = outputPtr->GetLargestPossibleRegion().GetSize();
151 //define some temp variables
152 signed long baseIndex[ImageDimension];
153 double distance[ImageDimension];
154 unsigned int dim, counter, upper;
155 double totalOverlap,overlap;
156 typename OutputImageType::IndexType neighIndex;
158 //Find the number of neighbors
159 unsigned int neighbors = 1 << ImageDimension;
161 //==================================================================================================
162 //Loop over the output region and add the intensities from the input to the output and the weight to the weights
163 //==================================================================================================
164 while( !inputIt.IsAtEnd() ) {
165 // get the input image index
166 index = inputIt.GetIndex();
167 inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
169 // get the required displacement
170 displacement = inputIt.Get();
172 // compute the required output image point
173 for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
175 // Update the output and the weights
176 if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) {
177 for(dim = 0; dim < ImageDimension; dim++) {
178 // The following block is equivalent to the following line without
179 // having to call floor. (Only for positive inputs, we already now that is in the image)
180 // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
182 baseIndex[dim] = (long) contIndex[dim];
183 distance[dim] = contIndex[dim] - double( baseIndex[dim] );
186 //Add contribution for each neighbor
187 totalOverlap = itk::NumericTraits<double>::Zero;
188 for( counter = 0; counter < neighbors ; counter++ ) {
189 overlap = 1.0; // fraction overlap
190 upper = counter; // each bit indicates upper/lower neighbour
192 // get neighbor index and overlap fraction
193 for( dim = 0; dim < ImageDimension; dim++ ) {
195 neighIndex[dim] = baseIndex[dim] + 1;
196 overlap *= distance[dim];
198 neighIndex[dim] = baseIndex[dim];
199 overlap *= 1.0 - distance[dim];
203 if (neighIndex[dim] >= size[dim])
204 neighIndex[dim] = size[dim] - 1;
208 //Set neighbor value only if overlap is not zero
209 if( (overlap>0.0)) // &&
210 // (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
211 // (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
212 // (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
213 // (neighIndex[0]>=0) &&
214 // (neighIndex[1]>=0) &&
215 // (neighIndex[2]>=0) )
217 //what to store? the original displacement vector?
218 if (! m_ThreadSafe) {
219 //Set the pixel and weight at neighIndex
220 outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
221 m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
225 //Entering critilal section: shared memory
226 m_MutexImage->GetPixel(neighIndex).Lock();
228 //Set the pixel and weight at neighIndex
229 outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
230 m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
233 m_MutexImage->GetPixel(neighIndex).Unlock();
236 //Add to total overlap
237 totalOverlap += overlap;
240 if( totalOverlap == 1.0 ) {
250 // std::cout << "HelperClass1::ThreadedGenerateData - OUT " << threadId << std::endl;
255 //=========================================================================================================================
256 //helper class 2 to allow a threaded execution of normalisation by the weights
257 //=========================================================================================================================
258 template<class InputImageType, class OutputImageType> class HelperClass2 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
262 /** Standard class typedefs. */
263 typedef HelperClass2 Self;
264 typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
265 typedef itk::SmartPointer<Self> Pointer;
266 typedef itk::SmartPointer<const Self> ConstPointer;
268 /** Method for creation through the object factory. */
271 /** Run-time type information (and related methods) */
272 itkTypeMacro( HelperClass2, ImageToImageFilter );
274 /** Constants for the image dimensions */
275 itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
278 typedef typename OutputImageType::PixelType PixelType;
279 typedef itk::Image<double,ImageDimension> WeightsImageType;
282 void SetWeights(const typename WeightsImageType::Pointer input) {
286 void SetEdgePaddingValue(PixelType value) {
287 m_EdgePaddingValue = value;
291 /** Typedef to describe the output image region type. */
292 typedef typename OutputImageType::RegionType OutputImageRegionType;
299 //the actual processing
300 #if ITK_VERSION_MAJOR >= 4
301 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId );
303 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
307 typename WeightsImageType::Pointer m_Weights;
308 PixelType m_EdgePaddingValue;
314 //=========================================================================================================================
315 //Member functions of the helper class 2
316 //=========================================================================================================================
319 //=========================================================================================================================
321 template<class InputImageType, class OutputImageType > HelperClass2<InputImageType, OutputImageType>::HelperClass2()
324 for(unsigned int i=0;i <PixelType::Dimension; i++) zero[i] = 0.0;
325 m_EdgePaddingValue=zero;
326 //m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
330 //=========================================================================================================================
331 //update the output for the outputRegionForThread
332 #if ITK_VERSION_MAJOR >= 4
333 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
335 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
338 // std::cout << "HelperClass2::ThreadedGenerateData - IN " << threadId << std::endl;
340 //Get pointer to the input
341 typename InputImageType::ConstPointer inputPtr = this->GetInput();
343 //Get pointer to the output
344 typename OutputImageType::Pointer outputPtr = this->GetOutput();
346 //Iterators over input, weigths and output
347 typedef itk::ImageRegionConstIterator<InputImageType> InputImageIteratorType;
348 typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
349 typedef itk::ImageRegionIterator<WeightsImageType> WeightsImageIteratorType;
351 //define them over the outputRegionForThread
352 OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
353 InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
354 WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread);
357 //==================================================================================================
358 //loop over the output and normalize the input, remove holes
359 PixelType neighValue;
360 double zero = itk::NumericTraits<double>::Zero;
361 while (!outputIt.IsAtEnd()) {
362 //the weight is not zero
363 if (weightsIt.Get() != zero) {
364 //divide by the weight
365 outputIt.Set(static_cast<PixelType>(inputIt.Get()/weightsIt.Get()));
368 //copy the value of the neighbour that was just processed
370 if(!outputIt.IsAtBegin()) {
374 //Neighbour cannot have zero weight because it should be filled already
375 neighValue=outputIt.Get();
377 outputIt.Set(neighValue);
380 //DD("is at begin, setting edgepadding value");
381 outputIt.Set(m_EdgePaddingValue);
390 // std::cout << "HelperClass2::ThreadedGenerateData - OUT " << threadId << std::endl;
395 }//end nameless namespace
402 //=========================================================================================================================
403 // The rest is the InvertVFFilter
404 //=========================================================================================================================
406 //=========================================================================================================================
408 template <class InputImageType, class OutputImageType>
409 InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
412 //m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
414 for(unsigned int i=0;i <PixelType::Dimension; i++) zero[i] = 0.0;
415 m_EdgePaddingValue=zero; //no other reasonable value?
421 //=========================================================================================================================
423 template <class InputImageType, class OutputImageType> void InvertVFFilter<InputImageType, OutputImageType>::GenerateData()
425 // std::cout << "InvertVFFilter::GenerateData - IN" << std::endl;
427 //Get the properties of the input
428 typename InputImageType::ConstPointer inputPtr=this->GetInput();
429 typename WeightsImageType::RegionType region = inputPtr->GetLargestPossibleRegion();
431 //Allocate the weights
432 typename WeightsImageType::Pointer weights=WeightsImageType::New();
433 weights->SetOrigin(inputPtr->GetOrigin());
434 weights->SetRegions(region);
436 weights->SetSpacing(inputPtr->GetSpacing());
438 //===========================================================================
439 //Inversion is divided in in two loops, for each we will call a threaded helper class
440 //1. add contribution of input to output and update weights
441 //2. normalize the output by the weight and remove holes
442 //===========================================================================
445 //===========================================================================
446 //1. add contribution of input to output and update weights
448 //Define an internal image type
450 typedef itk::Image<itk::Vector<double,ImageDimension>, ImageDimension > InternalImageType;
452 //Call threaded helper class 1
453 typedef HelperClass1<InputImageType, InternalImageType > HelperClass1Type;
454 typename HelperClass1Type::Pointer helper1=HelperClass1Type::New();
457 if(m_NumberOfThreadsIsGiven)helper1->SetNumberOfThreads(m_NumberOfThreads);
458 helper1->SetInput(inputPtr);
459 helper1->SetWeights(weights);
463 //Allocate the mutex image
464 typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New();
465 mutex->SetRegions(region);
467 mutex->SetSpacing(inputPtr->GetSpacing());
468 helper1->SetMutexImage(mutex);
469 if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <<std::endl;
470 } else if(m_Verbose)std::cout <<"Inverting using a thread-unsafe algorithm" <<std::endl;
472 //Execute helper class
476 typename InternalImageType::Pointer temp= helper1->GetOutput();
477 weights=helper1->GetWeights();
480 //===========================================================================
481 //2. Normalize the output by the weights and remove holes
482 //Call threaded helper class
483 typedef HelperClass2<InternalImageType, OutputImageType> HelperClass2Type;
484 typename HelperClass2Type::Pointer helper2=HelperClass2Type::New();
486 //Set temporary output as input
487 if(m_NumberOfThreadsIsGiven)helper2->SetNumberOfThreads(m_NumberOfThreads);
488 helper2->SetInput(temp);
489 helper2->SetWeights(weights);
490 helper2->SetEdgePaddingValue(m_EdgePaddingValue);
492 //Execute helper class
493 if (m_Verbose) std::cout << "Normalizing the output VF..."<<std::endl;
497 this->SetNthOutput(0, helper2->GetOutput());
499 //std::cout << "InvertVFFilter::GenerateData - OUT" << std::endl;