-namespace
-{//nameless namespace
-
- //=========================================================================================================================
- //helper class 1 to allow a threaded execution: add contributions of input to output and update weights
- //=========================================================================================================================
- template<class InputImageType, class OutputImageType, class DeformationFieldType> class HelperClass1 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
- {
-
- public:
- /** Standard class typedefs. */
- typedef HelperClass1 Self;
- typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- /** Method for creation through the object factory. */
- itkNewMacro(Self);
-
- /** Run-time type information (and related methods) */
- itkTypeMacro( HelperClass1, ImageToImageFilter );
-
- /** Constants for the image dimensions */
- itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
-
-
- //Typedefs
- typedef typename OutputImageType::PixelType OutputPixelType;
- typedef itk::Image<double, ImageDimension > WeightsImageType;
- typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
- //===================================================================================
- //Set methods
- void SetWeights(const typename WeightsImageType::Pointer input)
- {
- m_Weights = input;
- this->Modified();
- }
- void SetDeformationField(const typename DeformationFieldType::Pointer input)
- {
- m_DeformationField=input;
- this->Modified();
- }
- void SetMutexImage(const typename MutexImageType::Pointer input)
- {
- m_MutexImage=input;
- this->Modified();
- m_ThreadSafe=true;
+namespace
+{
+//nameless namespace
+
+//=========================================================================================================================
+//helper class 1 to allow a threaded execution: add contributions of input to output and update weights
+//=========================================================================================================================
+template<class InputImageType, class OutputImageType, class DeformationFieldType> class HelperClass1 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
+{
+
+public:
+ /** Standard class typedefs. */
+ typedef HelperClass1 Self;
+ typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods) */
+ itkTypeMacro( HelperClass1, ImageToImageFilter );
+
+ /** Constants for the image dimensions */
+ itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
+
+
+ //Typedefs
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef itk::Image<double, ImageDimension > WeightsImageType;
+ typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
+ //===================================================================================
+ //Set methods
+ void SetWeights(const typename WeightsImageType::Pointer input) {
+ m_Weights = input;
+ this->Modified();
+ }
+ void SetDeformationField(const typename DeformationFieldType::Pointer input) {
+ m_DeformationField=input;
+ this->Modified();
+ }
+ void SetMutexImage(const typename MutexImageType::Pointer input) {
+ m_MutexImage=input;
+ this->Modified();
+ m_ThreadSafe=true;
+ }
+
+ //Get methods
+ typename WeightsImageType::Pointer GetWeights() {
+ return m_Weights;
+ }
+
+ /** Typedef to describe the output image region type. */
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+protected:
+ HelperClass1();
+ ~HelperClass1() {};
+
+ //the actual processing
+ void BeforeThreadedGenerateData();
+ void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
+
+ //member data
+ typename itk::Image< double, ImageDimension>::Pointer m_Weights;
+ typename DeformationFieldType::Pointer m_DeformationField;
+ typename MutexImageType::Pointer m_MutexImage;
+ bool m_ThreadSafe;
+
+};
+
+
+
+//=========================================================================================================================
+//Member functions of the helper class 1
+//=========================================================================================================================
+
+
+//=========================================================================================================================
+//Empty constructor
+template<class InputImageType, class OutputImageType, class DeformationFieldType >
+HelperClass1<InputImageType, OutputImageType, DeformationFieldType>::HelperClass1()
+{
+ m_ThreadSafe=false;
+}
+
+
+//=========================================================================================================================
+//Before threaded data
+template<class InputImageType, class OutputImageType, class DeformationFieldType >
+void HelperClass1<InputImageType, OutputImageType, DeformationFieldType>::BeforeThreadedGenerateData()
+{
+ //Since we will add, put to zero!
+ this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
+ this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
+}
+
+
+//=========================================================================================================================
+//update the output for the outputRegionForThread
+template<class InputImageType, class OutputImageType, class DeformationFieldType >
+void HelperClass1<InputImageType, OutputImageType, DeformationFieldType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+
+ //Get pointer to the input
+ typename InputImageType::ConstPointer inputPtr = this->GetInput();
+
+ //Get pointer to the output
+ typename OutputImageType::Pointer outputPtr = this->GetOutput();
+ //typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
+
+ //Iterators over input and deformation field
+ typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
+ typedef itk::ImageRegionIterator<DeformationFieldType> DeformationFieldIteratorType;
+
+ //define them over the outputRegionForThread
+ InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
+ DeformationFieldIteratorType fieldIt(m_DeformationField,outputRegionForThread);
+
+ //Initialize
+ typename InputImageType::IndexType index;
+ itk::ContinuousIndex<double,ImageDimension> contIndex;
+ typename InputImageType::PointType point;
+ typedef typename DeformationFieldType::PixelType DisplacementType;
+ DisplacementType displacement;
+ fieldIt.GoToBegin();
+ inputIt.GoToBegin();
+
+ //define some temp variables
+ signed long baseIndex[ImageDimension];
+ double distance[ImageDimension];
+ for(unsigned int i=0; i<ImageDimension; i++) distance[i] = 0.0; // to avoid warning
+ unsigned int dim, counter, upper;
+ double overlap, totalOverlap;
+ typename OutputImageType::IndexType neighIndex;
+
+ //Find the number of neighbors
+ unsigned int neighbors = 1 << ImageDimension;
+
+
+ //==================================================================================================
+ //Loop over the region and add the intensities to the output and the weight to the weights
+ //==================================================================================================
+ while( !inputIt.IsAtEnd() ) {
+
+ // get the input image index
+ index = inputIt.GetIndex();
+ inputPtr->TransformIndexToPhysicalPoint( index, point );
+
+ // get the required displacement
+ displacement = fieldIt.Get();
+
+ // compute the required output image point
+ for(unsigned int j = 0; j < ImageDimension; j++ ) point[j] += displacement[j];
+
+
+ // Update the output and the weights
+ if(outputPtr->TransformPhysicalPointToContinuousIndex(point, contIndex ) ) {
+ for(dim = 0; dim < ImageDimension; dim++) {
+ // The following block is equivalent to the following line without
+ // having to call floor. For positive inputs!!!
+ // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
+ baseIndex[dim] = (long) contIndex[dim];
+ distance[dim] = contIndex[dim] - double( baseIndex[dim] );
+ }
+
+ //Add contribution for each neighbor
+ totalOverlap = itk::NumericTraits<double>::Zero;
+ for( counter = 0; counter < neighbors ; counter++ ) {
+ overlap = 1.0; // fraction overlap
+ upper = counter; // each bit indicates upper/lower neighbour
+
+ // get neighbor index and overlap fraction
+ for( dim = 0; dim < ImageDimension; dim++ ) {
+ if ( upper & 1 ) {
+ neighIndex[dim] = baseIndex[dim] + 1;
+ overlap *= distance[dim];
+ } else {
+ neighIndex[dim] = baseIndex[dim];
+ overlap *= 1.0 - distance[dim];
+ }
+ upper >>= 1;
+ }
+
+ //Set neighbor value only if overlap is not zero
+ if( (overlap>0.0)) // &&
+ // (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
+ // (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
+ // (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
+ // (neighIndex[0]>=0) &&
+ // (neighIndex[1]>=0) &&
+ // (neighIndex[2]>=0) )
+ {
+
+ if (! m_ThreadSafe) {
+ //Set the pixel and weight at neighIndex
+ outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) + overlap * static_cast<OutputPixelType>(inputIt.Get()));
+ m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
+
+ } else {
+ //Entering critilal section: shared memory
+ m_MutexImage->GetPixel(neighIndex).Lock();
+
+ //Set the pixel and weight at neighIndex
+ outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) + overlap * static_cast<OutputPixelType>(inputIt.Get()));
+ m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
+
+ //Unlock
+ m_MutexImage->GetPixel(neighIndex).Unlock();
+
+ }
+ //Add to total overlap
+ totalOverlap += overlap;
+ }
+
+ //check for totaloverlap: not very likely
+ if( totalOverlap == 1.0 ) {
+ // finished
+ break;
+ }
+ }