1 #ifndef __clitkInvertVFFilter_txx
2 #define __clitkInvertVFFilter_txx
5 // Put the helper classes in an anonymous namespace so that it is not
11 //=========================================================================================================================
12 //helper class 1 to allow a threaded execution: add contributions of input to output and update weights
13 //=========================================================================================================================
14 template<class InputImageType, class OutputImageType> class ITK_EXPORT HelperClass1 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
18 /** Standard class typedefs. */
19 typedef HelperClass1 Self;
20 typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
21 typedef itk::SmartPointer<Self> Pointer;
22 typedef itk::SmartPointer<const Self> ConstPointer;
24 /** Method for creation through the object factory. */
27 /** Run-time type information (and related methods) */
28 itkTypeMacro( HelperClass1, ImageToImageFilter );
30 /** Constants for the image dimensions */
31 itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
35 typedef typename OutputImageType::PixelType PixelType;
36 typedef itk::Image<double, ImageDimension > WeightsImageType;
37 typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
39 //===================================================================================
41 void SetWeights(const typename WeightsImageType::Pointer input)
46 void SetMutexImage(const typename MutexImageType::Pointer input)
54 typename WeightsImageType::Pointer GetWeights(){return m_Weights;}
56 /** Typedef to describe the output image region type. */
57 typedef typename OutputImageType::RegionType OutputImageRegionType;
63 //the actual processing
64 void BeforeThreadedGenerateData();
65 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
68 typename WeightsImageType::Pointer m_Weights;
69 typename MutexImageType::Pointer m_MutexImage;
76 //=========================================================================================================================
77 //Member functions of the helper class 1
78 //=========================================================================================================================
81 //=========================================================================================================================
83 template<class InputImageType, class OutputImageType >
84 HelperClass1<InputImageType, OutputImageType>::HelperClass1()
89 //=========================================================================================================================
90 //Before threaded data
91 template<class InputImageType, class OutputImageType >
92 void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
94 //Since we will add, put to zero!
95 this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
96 this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
99 //=========================================================================================================================
100 //update the output for the outputRegionForThread
101 template<class InputImageType, class OutputImageType>
102 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
105 //Get pointer to the input
106 typename InputImageType::ConstPointer inputPtr = this->GetInput();
108 //Get pointer to the output
109 typename OutputImageType::Pointer outputPtr = this->GetOutput();
110 typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
112 //Iterator over input
113 typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
115 //define them over the outputRegionForThread
116 InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
119 typename InputImageType::IndexType index;
120 itk::ContinuousIndex<double,ImageDimension> contIndex;
121 typename InputImageType::PointType ipoint;
122 typename OutputImageType::PointType opoint;
123 typedef typename OutputImageType::PixelType DisplacementType;
124 DisplacementType displacement;
127 //define some temp variables
128 signed long baseIndex[ImageDimension];
129 double distance[ImageDimension];
130 unsigned int dim, counter, upper;
131 double totalOverlap,overlap;
132 typename OutputImageType::IndexType neighIndex;
134 //Find the number of neighbors
135 unsigned int neighbors = 1 << ImageDimension;
137 //==================================================================================================
138 //Loop over the region and add the intensities to the output and the weight to the weights
139 //==================================================================================================
140 while( !inputIt.IsAtEnd() )
142 // get the input image index
143 index = inputIt.GetIndex();
144 inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
146 // get the required displacement
147 displacement = inputIt.Get();
149 // compute the required output image point
150 for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
152 // Update the output and the weights
153 if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) )
155 for(dim = 0; dim < ImageDimension; dim++)
157 // The following block is equivalent to the following line without
158 // having to call floor. (Only for positive inputs, we already now that is in the image)
159 // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
161 baseIndex[dim] = (long) contIndex[dim];
162 distance[dim] = contIndex[dim] - double( baseIndex[dim] );
165 //Add contribution for each neighbor
166 totalOverlap = itk::NumericTraits<double>::Zero;
167 for( counter = 0; counter < neighbors ; counter++ )
169 overlap = 1.0; // fraction overlap
170 upper = counter; // each bit indicates upper/lower neighbour
172 // get neighbor index and overlap fraction
173 for( dim = 0; dim < 3; dim++ )
177 neighIndex[dim] = baseIndex[dim] + 1;
178 overlap *= distance[dim];
182 neighIndex[dim] = baseIndex[dim];
183 overlap *= 1.0 - distance[dim];
190 //Set neighbor value only if overlap is not zero
191 if( (overlap>0.0)) // &&
192 // (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
193 // (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
194 // (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
195 // (neighIndex[0]>=0) &&
196 // (neighIndex[1]>=0) &&
197 // (neighIndex[2]>=0) )
199 //what to store? the original displacement vector?
202 //Set the pixel and weight at neighIndex
203 outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
204 m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
209 //Entering critilal section: shared memory
210 m_MutexImage->GetPixel(neighIndex).Lock();
212 //Set the pixel and weight at neighIndex
213 outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
214 m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
217 m_MutexImage->GetPixel(neighIndex).Unlock();
220 //Add to total overlap
221 totalOverlap += overlap;
224 if( totalOverlap == 1.0 )
239 //=========================================================================================================================
240 //helper class 2 to allow a threaded execution of normalisation by the weights
241 //=========================================================================================================================
242 template<class InputImageType, class OutputImageType> class HelperClass2 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
246 /** Standard class typedefs. */
247 typedef HelperClass2 Self;
248 typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
249 typedef itk::SmartPointer<Self> Pointer;
250 typedef itk::SmartPointer<const Self> ConstPointer;
252 /** Method for creation through the object factory. */
255 /** Run-time type information (and related methods) */
256 itkTypeMacro( HelperClass2, ImageToImageFilter );
258 /** Constants for the image dimensions */
259 itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
262 typedef typename OutputImageType::PixelType PixelType;
263 typedef itk::Image<double,ImageDimension> WeightsImageType;
266 void SetWeights(const typename WeightsImageType::Pointer input)
271 void SetEdgePaddingValue(PixelType value)
273 m_EdgePaddingValue = value;
277 /** Typedef to describe the output image region type. */
278 typedef typename OutputImageType::RegionType OutputImageRegionType;
285 //the actual processing
286 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
290 typename WeightsImageType::Pointer m_Weights;
291 PixelType m_EdgePaddingValue;
297 //=========================================================================================================================
298 //Member functions of the helper class 2
299 //=========================================================================================================================
302 //=========================================================================================================================
304 template<class InputImageType, class OutputImageType > HelperClass2<InputImageType, OutputImageType>::HelperClass2()
306 m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
310 //=========================================================================================================================
311 //update the output for the outputRegionForThread
312 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
315 //Get pointer to the input
316 typename InputImageType::ConstPointer inputPtr = this->GetInput();
318 //Get pointer to the output
319 typename OutputImageType::Pointer outputPtr = this->GetOutput();
321 //Iterators over input, weigths and output
322 typedef itk::ImageRegionConstIterator<InputImageType> InputImageIteratorType;
323 typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
324 typedef itk::ImageRegionIterator<WeightsImageType> WeightsImageIteratorType;
326 //define them over the outputRegionForThread
327 OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
328 InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
329 WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread);
332 //==================================================================================================
333 //loop over the output and normalize the input, remove holes
334 PixelType neighValue;
335 double zero = itk::NumericTraits<double>::Zero;
336 while (!outputIt.IsAtEnd())
338 //the weight is not zero
339 if (weightsIt.Get() != zero)
341 //divide by the weight
342 outputIt.Set(static_cast<PixelType>(inputIt.Get()/weightsIt.Get()));
345 //copy the value of the neighbour that was just processed
348 if(!outputIt.IsAtBegin())
353 //Neighbour cannot have zero weight because it should be filled already
354 neighValue=outputIt.Get();
356 outputIt.Set(neighValue);
360 //DD("is at begin, setting edgepadding value");
361 outputIt.Set(m_EdgePaddingValue);
372 }//end nameless namespace
379 //=========================================================================================================================
380 // The rest is the InvertVFFilter
381 //=========================================================================================================================
383 //=========================================================================================================================
385 template <class InputImageType, class OutputImageType>
386 InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
388 m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
394 //=========================================================================================================================
396 template <class InputImageType, class OutputImageType> void InvertVFFilter<InputImageType, OutputImageType>::GenerateData()
399 //Get the properties of the input
400 typename InputImageType::ConstPointer inputPtr=this->GetInput();
401 typename WeightsImageType::RegionType region;
402 typename WeightsImageType::RegionType::SizeType size=inputPtr->GetLargestPossibleRegion().GetSize();
403 region.SetSize(size);
404 typename OutputImageType::IndexType start;
405 for (unsigned int i=0; i< ImageDimension; i++) start[i]=0;
406 region.SetIndex(start);
407 PixelType zero = itk::NumericTraits<double>::Zero;
410 //Allocate the weights
411 typename WeightsImageType::Pointer weights=WeightsImageType::New();
412 weights->SetRegions(region);
414 weights->SetSpacing(inputPtr->GetSpacing());
416 //===========================================================================
417 //Inversion is divided in in two loops, for each we will call a threaded helper class
418 //1. add contribution of input to output and update weights
419 //2. normalize the output by the weight and remove holes
420 //===========================================================================
423 //===========================================================================
424 //1. add contribution of input to output and update weights
426 //Define an internal image type
428 typedef itk::Image<itk::Vector<double,ImageDimension>, ImageDimension > InternalImageType;
430 //Call threaded helper class 1
431 typedef HelperClass1<InputImageType, InternalImageType > HelperClass1Type;
432 typename HelperClass1Type::Pointer helper1=HelperClass1Type::New();
435 if(m_NumberOfThreadsIsGiven)helper1->SetNumberOfThreads(m_NumberOfThreads);
436 helper1->SetInput(inputPtr);
437 helper1->SetWeights(weights);
442 //Allocate the mutex image
443 typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New();
444 mutex->SetRegions(region);
446 mutex->SetSpacing(inputPtr->GetSpacing());
447 helper1->SetMutexImage(mutex);
448 if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <<std::endl;
450 else if(m_Verbose)std::cout <<"Inverting using a thread-unsafe algorithm" <<std::endl;
452 //Execute helper class
456 typename InternalImageType::Pointer temp= helper1->GetOutput();
457 weights=helper1->GetWeights();
460 //===========================================================================
461 //2. Normalize the output by the weights and remove holes
462 //Call threaded helper class
463 typedef HelperClass2<InternalImageType, OutputImageType> HelperClass2Type;
464 typename HelperClass2Type::Pointer helper2=HelperClass2Type::New();
466 //Set temporary output as input
467 helper2->SetInput(temp);
468 helper2->SetWeights(weights);
469 helper2->SetEdgePaddingValue(m_EdgePaddingValue);
471 //Execute helper class
472 if (m_Verbose) std::cout << "Normalizing the output VF..."<<std::endl;
476 this->SetNthOutput(0, helper2->GetOutput());