X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkForwardWarpImageFilter.txx;h=f5e98536516cb5e1657d577c50b1c8bff54d1687;hb=595624eb4297e747630105d45017de69733caef2;hp=38836be401d2bba9afef8d84c850d299eb20c657;hpb=0b7c9b1e1215634b02cbd38d4e4ba101d6111ba8;p=clitk.git diff --git a/itk/clitkForwardWarpImageFilter.txx b/itk/clitkForwardWarpImageFilter.txx index 38836be..f5e9853 100644 --- a/itk/clitkForwardWarpImageFilter.txx +++ b/itk/clitkForwardWarpImageFilter.txx @@ -1,9 +1,9 @@ /*========================================================================= Program: vv http://www.creatis.insa-lyon.fr/rio/vv - Authors belong to: + Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -======================================================================-====*/ +===========================================================================**/ #ifndef __clitkForwardWarpImageFilter_txx #define __clitkForwardWarpImageFilter_txx #include "clitkForwardWarpImageFilter.h" @@ -23,382 +23,364 @@ // Put the helper classes in an anonymous namespace so that it is not // exposed to the user -namespace -{//nameless namespace - - //========================================================================================================================= - //helper class 1 to allow a threaded execution: add contributions of input to output and update weights - //========================================================================================================================= - template class HelperClass1 : public itk::ImageToImageFilter - { - - public: - /** Standard class typedefs. */ - typedef HelperClass1 Self; - typedef itk::ImageToImageFilter Superclass; - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer 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 WeightsImageType; - typedef itk::Image 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 HelperClass1 : public itk::ImageToImageFilter +{ + +public: + /** Standard class typedefs. */ + typedef HelperClass1 Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer 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 WeightsImageType; + typedef itk::Image 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 +HelperClass1::HelperClass1() +{ + m_ThreadSafe=false; +} + + +//========================================================================================================================= +//Before threaded data +template +void HelperClass1::BeforeThreadedGenerateData() +{ + //Since we will add, put to zero! + this->GetOutput()->FillBuffer(itk::NumericTraits::Zero); + this->GetWeights()->FillBuffer(itk::NumericTraits::Zero); +} + + +//========================================================================================================================= +//update the output for the outputRegionForThread +template +void HelperClass1::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 InputImageIteratorType; + typedef itk::ImageRegionIterator DeformationFieldIteratorType; + + //define them over the outputRegionForThread + InputImageIteratorType inputIt(inputPtr, outputRegionForThread); + DeformationFieldIteratorType fieldIt(m_DeformationField,outputRegionForThread); + + //Initialize + typename InputImageType::IndexType index; + itk::ContinuousIndex 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; iTransformIndexToPhysicalPoint( 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::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(neighIndex[0])(neighIndex[1])(neighIndex[2])=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(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(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; + } + } } - - //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 - HelperClass1::HelperClass1() - { - m_ThreadSafe=false; + + ++fieldIt; + ++inputIt; } - //========================================================================================================================= - //Before threaded data - template - void HelperClass1::BeforeThreadedGenerateData() - { - //Since we will add, put to zero! - this->GetOutput()->FillBuffer(itk::NumericTraits::Zero); - this->GetWeights()->FillBuffer(itk::NumericTraits::Zero); +} + + + +//========================================================================================================================= +//helper class 2 to allow a threaded execution of normalisation by the weights +//========================================================================================================================= +template +class HelperClass2 : public itk::ImageToImageFilter +{ + +public: + /** Standard class typedefs. */ + typedef HelperClass2 Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods) */ + itkTypeMacro( HelperClass2, ImageToImageFilter ); + + /** Constants for the image dimensions */ + itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension); + + //Typedefs + typedef typename InputImageType::PixelType InputPixelType; + typedef typename OutputImageType::PixelType OutputPixelType; + typedef itk::Image WeightsImageType; + + + //Set methods + void SetWeights(const typename WeightsImageType::Pointer input) { + m_Weights = input; + this->Modified(); } - - - //========================================================================================================================= - //update the output for the outputRegionForThread - template - void HelperClass1::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 InputImageIteratorType; - typedef itk::ImageRegionIterator DeformationFieldIteratorType; - - //define them over the outputRegionForThread - InputImageIteratorType inputIt(inputPtr, outputRegionForThread); - DeformationFieldIteratorType fieldIt(m_DeformationField,outputRegionForThread); - - //Initialize - typename InputImageType::IndexType index; - itk::ContinuousIndex 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]; - 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::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 < 3; 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(neighIndex[0])(neighIndex[1])(neighIndex[2])=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(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(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; - } - } - } - - ++fieldIt; - ++inputIt; - } - - + void SetEdgePaddingValue(OutputPixelType value) { + m_EdgePaddingValue = value; + this->Modified(); } + /** Typedef to describe the output image region type. */ + typedef typename OutputImageType::RegionType OutputImageRegionType; + +protected: + HelperClass2(); + ~HelperClass2() {}; + + //the actual processing + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); - //========================================================================================================================= - //helper class 2 to allow a threaded execution of normalisation by the weights - //========================================================================================================================= - template - class HelperClass2 : public itk::ImageToImageFilter - { - - public: - /** Standard class typedefs. */ - typedef HelperClass2 Self; - typedef itk::ImageToImageFilter Superclass; - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods) */ - itkTypeMacro( HelperClass2, ImageToImageFilter ); - - /** Constants for the image dimensions */ - itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension); - - //Typedefs - typedef typename InputImageType::PixelType InputPixelType; - typedef typename OutputImageType::PixelType OutputPixelType; - typedef itk::Image WeightsImageType; - - - //Set methods - void SetWeights(const typename WeightsImageType::Pointer input) - { - m_Weights = input; - this->Modified(); + + //member data + typename WeightsImageType::Pointer m_Weights; + OutputPixelType m_EdgePaddingValue; +} ; + + + +//========================================================================================================================= +//Member functions of the helper class 2 +//========================================================================================================================= + + +//========================================================================================================================= +//Empty constructor +template +HelperClass2::HelperClass2() +{ + m_EdgePaddingValue=static_cast(0.0); +} + + +//========================================================================================================================= +//update the output for the outputRegionForThread +template void +HelperClass2::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(); + + //Iterators over input, weigths and output + typedef itk::ImageRegionConstIterator InputImageIteratorType; + typedef itk::ImageRegionIterator OutputImageIteratorType; + typedef itk::ImageRegionIterator WeightsImageIteratorType; + + //define them over the outputRegionForThread + OutputImageIteratorType outputIt(outputPtr, outputRegionForThread); + InputImageIteratorType inputIt(inputPtr, outputRegionForThread); + WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread); + + + //================================================================================================== + //loop over the output and normalize the input, remove holes + OutputPixelType neighValue; + double zero = itk::NumericTraits::Zero; + while (!outputIt.IsAtEnd()) { + //the weight is not zero + if (weightsIt.Get() != zero) { + //divide by the weight + outputIt.Set(static_cast(inputIt.Get()/weightsIt.Get())); } - void SetEdgePaddingValue(OutputPixelType value) - { - m_EdgePaddingValue = value; - this->Modified(); + + //copy the value of the neighbour that was just processed + else { + if(!outputIt.IsAtBegin()) { + //go back + --outputIt; + neighValue=outputIt.Get(); + ++outputIt; + outputIt.Set(neighValue); + } else { + //DD("is at begin, setting edgepadding value"); + outputIt.Set(m_EdgePaddingValue); + } } - - /** Typedef to describe the output image region type. */ - typedef typename OutputImageType::RegionType OutputImageRegionType; - - protected: - HelperClass2(); - ~HelperClass2(){}; - - - //the actual processing - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); - - - //member data - typename WeightsImageType::Pointer m_Weights; - OutputPixelType m_EdgePaddingValue; - } ; - - - - //========================================================================================================================= - //Member functions of the helper class 2 - //========================================================================================================================= - - - //========================================================================================================================= - //Empty constructor - template - HelperClass2::HelperClass2() - { - m_EdgePaddingValue=static_cast(0.0); - } - - - //========================================================================================================================= - //update the output for the outputRegionForThread - template void - HelperClass2::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(); - - //Iterators over input, weigths and output - typedef itk::ImageRegionConstIterator InputImageIteratorType; - typedef itk::ImageRegionIterator OutputImageIteratorType; - typedef itk::ImageRegionIterator WeightsImageIteratorType; - - //define them over the outputRegionForThread - OutputImageIteratorType outputIt(outputPtr, outputRegionForThread); - InputImageIteratorType inputIt(inputPtr, outputRegionForThread); - WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread); - - - //================================================================================================== - //loop over the output and normalize the input, remove holes - OutputPixelType neighValue; - double zero = itk::NumericTraits::Zero; - while (!outputIt.IsAtEnd()) - { - //the weight is not zero - if (weightsIt.Get() != zero) - { - //divide by the weight - outputIt.Set(static_cast(inputIt.Get()/weightsIt.Get())); - } - - //copy the value of the neighbour that was just processed - else - { - if(!outputIt.IsAtBegin()) - { - //go back - --outputIt; - neighValue=outputIt.Get(); - ++outputIt; - outputIt.Set(neighValue); - } - else{ - //DD("is at begin, setting edgepadding value"); - outputIt.Set(m_EdgePaddingValue); - } - } - ++weightsIt; - ++outputIt; - ++inputIt; - - }//end while - }//end member - + ++weightsIt; + ++outputIt; + ++inputIt; + + }//end while +}//end member + }//end nameless namespace @@ -407,109 +389,107 @@ namespace namespace clitk { - //========================================================================================================================= - // The rest is the ForwardWarpImageFilter - //========================================================================================================================= - - //========================================================================================================================= - //constructor - template - ForwardWarpImageFilter::ForwardWarpImageFilter() - { - // mIsUpdated=false; - m_NumberOfThreadsIsGiven=false; - m_EdgePaddingValue=static_cast(0.0); - m_ThreadSafe=false; - m_Verbose=false; - } +//========================================================================================================================= +// The rest is the ForwardWarpImageFilter +//========================================================================================================================= +//========================================================================================================================= +//constructor +template +ForwardWarpImageFilter::ForwardWarpImageFilter() +{ + // mIsUpdated=false; + m_NumberOfThreadsIsGiven=false; + m_EdgePaddingValue=static_cast(0.0); + m_ThreadSafe=false; + m_Verbose=false; +} + + +//========================================================================================================================= +//Update +template +void ForwardWarpImageFilter::GenerateData() +{ + + //Get the properties of the input + typename InputImageType::ConstPointer inputPtr=this->GetInput(); + typename WeightsImageType::RegionType region; + typename WeightsImageType::RegionType::SizeType size=inputPtr->GetLargestPossibleRegion().GetSize(); + region.SetSize(size); + typename OutputImageType::IndexType start; + for (unsigned int i =0; i< ImageDimension ; i ++)start[i]=0; + region.SetIndex(start); + + //Allocate the weights + typename WeightsImageType::Pointer weights=ForwardWarpImageFilter::WeightsImageType::New(); + weights->SetRegions(region); + weights->Allocate(); + weights->SetSpacing(inputPtr->GetSpacing()); + + + //=========================================================================== + //warp is divided in in two loops, for each we call a threaded helper class + //1. Add contribution of input to output and update weights + //2. Normalize the output by the weight and remove holes + //=========================================================================== + + //=========================================================================== + //1. Add contribution of input to output and update weights + + //Define an internal image type in double precision + typedef itk::Image InternalImageType; + + //Call threaded helper class 1 + typedef HelperClass1 HelperClass1Type; + typename HelperClass1Type::Pointer helper1=HelperClass1Type::New(); + + //Set input + if(m_NumberOfThreadsIsGiven)helper1->SetNumberOfThreads(m_NumberOfThreads); + helper1->SetInput(inputPtr); + helper1->SetDeformationField(m_DeformationField); + helper1->SetWeights(weights); + + //Threadsafe? + if(m_ThreadSafe) { + //Allocate the mutex image + typename MutexImageType::Pointer mutex=ForwardWarpImageFilter::MutexImageType::New(); + mutex->SetRegions(region); + mutex->Allocate(); + mutex->SetSpacing(inputPtr->GetSpacing()); + helper1->SetMutexImage(mutex); + if (m_Verbose) std::cout <<"Forwarp warping using a thread-safe algorithm" <Update(); + + //Get the output + typename InternalImageType::Pointer temp= helper1->GetOutput(); + + //For clarity + weights=helper1->GetWeights(); + + + //=========================================================================== + //2. Normalize the output by the weights and remove holes + //Call threaded helper class + typedef HelperClass2 HelperClass2Type; + typename HelperClass2Type::Pointer helper2=HelperClass2Type::New(); + + //Set temporary output as input + if(m_NumberOfThreadsIsGiven)helper2->SetNumberOfThreads(m_NumberOfThreads); + helper2->SetInput(temp); + helper2->SetWeights(weights); + helper2->SetEdgePaddingValue(m_EdgePaddingValue); + + //Execute helper class + helper2->Update(); + + //Set the output + this->SetNthOutput(0, helper2->GetOutput()); +} - //========================================================================================================================= - //Update - template - void ForwardWarpImageFilter::GenerateData() - { - - //Get the properties of the input - typename InputImageType::ConstPointer inputPtr=this->GetInput(); - typename WeightsImageType::RegionType region; - typename WeightsImageType::RegionType::SizeType size=inputPtr->GetLargestPossibleRegion().GetSize(); - region.SetSize(size); - typename OutputImageType::IndexType start; - for (unsigned int i =0; i< ImageDimension ;i ++)start[i]=0; - region.SetIndex(start); - - //Allocate the weights - typename WeightsImageType::Pointer weights=ForwardWarpImageFilter::WeightsImageType::New(); - weights->SetRegions(region); - weights->Allocate(); - weights->SetSpacing(inputPtr->GetSpacing()); - - - //=========================================================================== - //warp is divided in in two loops, for each we call a threaded helper class - //1. Add contribution of input to output and update weights - //2. Normalize the output by the weight and remove holes - //=========================================================================== - - //=========================================================================== - //1. Add contribution of input to output and update weights - - //Define an internal image type in double precision - typedef itk::Image InternalImageType; - - //Call threaded helper class 1 - typedef HelperClass1 HelperClass1Type; - typename HelperClass1Type::Pointer helper1=HelperClass1Type::New(); - - //Set input - if(m_NumberOfThreadsIsGiven)helper1->SetNumberOfThreads(m_NumberOfThreads); - helper1->SetInput(inputPtr); - helper1->SetDeformationField(m_DeformationField); - helper1->SetWeights(weights); - - //Threadsafe? - if(m_ThreadSafe) - { - //Allocate the mutex image - typename MutexImageType::Pointer mutex=ForwardWarpImageFilter::MutexImageType::New(); - mutex->SetRegions(region); - mutex->Allocate(); - mutex->SetSpacing(inputPtr->GetSpacing()); - helper1->SetMutexImage(mutex); - if (m_Verbose) std::cout <<"Forwarp warping using a thread-safe algorithm" <Update(); - - //Get the output - typename InternalImageType::Pointer temp= helper1->GetOutput(); - - //For clarity - weights=helper1->GetWeights(); - - - //=========================================================================== - //2. Normalize the output by the weights and remove holes - //Call threaded helper class - typedef HelperClass2 HelperClass2Type; - typename HelperClass2Type::Pointer helper2=HelperClass2Type::New(); - - //Set temporary output as input - if(m_NumberOfThreadsIsGiven)helper2->SetNumberOfThreads(m_NumberOfThreads); - helper2->SetInput(temp); - helper2->SetWeights(weights); - helper2->SetEdgePaddingValue(m_EdgePaddingValue); - - //Execute helper class - helper2->Update(); - - //Set the output - this->SetNthOutput(0, helper2->GetOutput()); - } - } #endif