X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkInvertVFFilter.txx;h=6bda75eb0d2b65ba4a21b6f438862201c6fffd14;hb=a6957c4825e83c61b977ec316dd841878617ffbd;hp=a326c030db79e897c915e472a916aff582dc6e38;hpb=0083c3fb2c66812489631c7551709d121de51625;p=clitk.git diff --git a/itk/clitkInvertVFFilter.txx b/itk/clitkInvertVFFilter.txx index a326c03..6bda75e 100644 --- a/itk/clitkInvertVFFilter.txx +++ b/itk/clitkInvertVFFilter.txx @@ -1,478 +1,475 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.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 + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ #ifndef __clitkInvertVFFilter_txx #define __clitkInvertVFFilter_txx namespace { - //========================================================================================================================= - //helper class 1 to allow a threaded execution: add contributions of input to output and update weights - //========================================================================================================================= - template class ITK_EXPORT 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 PixelType; - typedef itk::Image WeightsImageType; - typedef itk::Image MutexImageType; - - //=================================================================================== - //Set methods - void SetWeights(const typename WeightsImageType::Pointer input) - { - m_Weights = 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 ); +//========================================================================================================================= +//helper class 1 to allow a threaded execution: add contributions of input to output and update weights +//========================================================================================================================= +template class ITK_EXPORT HelperClass1 : public itk::ImageToImageFilter +{ - //member data - typename WeightsImageType::Pointer m_Weights; - typename MutexImageType::Pointer m_MutexImage; - bool m_ThreadSafe; +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); - //========================================================================================================================= - //Member functions of the helper class 1 - //========================================================================================================================= + //Typedefs + typedef typename OutputImageType::PixelType PixelType; + typedef itk::Image WeightsImageType; + typedef itk::Image MutexImageType; - //========================================================================================================================= - //Empty constructor - template - HelperClass1::HelperClass1() - { - m_ThreadSafe=false; + //=================================================================================== + //Set methods + void SetWeights(const typename WeightsImageType::Pointer input) { + m_Weights = input; + this->Modified(); + } + void SetMutexImage(const typename MutexImageType::Pointer input) { + m_MutexImage=input; + this->Modified(); + m_ThreadSafe=true; } - //========================================================================================================================= - //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); + //Get methods + typename WeightsImageType::Pointer GetWeights() { + return m_Weights; } - //========================================================================================================================= - //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(); - - //Iterator over input - typedef itk::ImageRegionConstIteratorWithIndex InputImageIteratorType; - - //define them over the outputRegionForThread - InputImageIteratorType inputIt(inputPtr, outputRegionForThread); - - //Initialize - typename InputImageType::IndexType index; - itk::ContinuousIndex contIndex; - typename InputImageType::PointType ipoint; - typename OutputImageType::PointType opoint; - typedef typename OutputImageType::PixelType DisplacementType; - DisplacementType displacement; - inputIt.GoToBegin(); - - //define some temp variables - signed long baseIndex[ImageDimension]; - double distance[ImageDimension]; - unsigned int dim, counter, upper; - double totalOverlap,overlap; - 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,ipoint ); - - // get the required displacement - displacement = inputIt.Get(); - - // compute the required output image point - for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j]; - - // Update the output and the weights - if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) - { - for(dim = 0; dim < ImageDimension; dim++) - { - // The following block is equivalent to the following line without - // having to call floor. (Only for positive inputs, we already now that is in the image) - // 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) ) - { - //what to store? the original displacement vector? - if (! m_ThreadSafe) - { - //Set the pixel and weight at neighIndex - outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap)); - 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) - (displacement*overlap)); - m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap); - - //Unlock - m_MutexImage->GetPixel(neighIndex).Unlock(); - - } - //Add to total overlap - totalOverlap += overlap; - } - - if( totalOverlap == 1.0 ) - { - // finished - break; - } - } - } - - ++inputIt; + /** 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, itk::ThreadIdType threadId ); + + //member data + typename WeightsImageType::Pointer m_Weights; + 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, itk::ThreadIdType 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(); + + //Iterator over input + typedef itk::ImageRegionConstIteratorWithIndex InputImageIteratorType; + + //define them over the outputRegionForThread + InputImageIteratorType inputIt(inputPtr, outputRegionForThread); + + //Initialize + typename InputImageType::IndexType index; + itk::ContinuousIndex contIndex; + typename InputImageType::PointType ipoint; + typename OutputImageType::PointType opoint; + typedef typename OutputImageType::PixelType DisplacementType; + DisplacementType displacement; + inputIt.GoToBegin(); + + //define some temp variables + signed long baseIndex[ImageDimension]; + double distance[ImageDimension]; + unsigned int dim, counter, upper; + double totalOverlap,overlap; + 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,ipoint ); + + // get the required displacement + displacement = inputIt.Get(); + + // compute the required output image point + for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j]; + + // Update the output and the weights + if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) { + for(dim = 0; dim < ImageDimension; dim++) { + // The following block is equivalent to the following line without + // having to call floor. (Only for positive inputs, we already now that is in the image) + // 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) ) + { + //what to store? the original displacement vector? + if (! m_ThreadSafe) { + //Set the pixel and weight at neighIndex + outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap)); + 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) - (displacement*overlap)); + m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap); + + //Unlock + m_MutexImage->GetPixel(neighIndex).Unlock(); + + } + //Add to total overlap + totalOverlap += overlap; + } + + if( totalOverlap == 1.0 ) { + // finished + break; + } + } + } + + ++inputIt; } +} - //========================================================================================================================= - //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 OutputImageType::PixelType PixelType; - typedef itk::Image WeightsImageType; - - //Set methods - void SetWeights(const typename WeightsImageType::Pointer input) - { - m_Weights = input; - this->Modified(); - } - void SetEdgePaddingValue(PixelType 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; - //member data - typename WeightsImageType::Pointer m_Weights; - PixelType m_EdgePaddingValue; + /** 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 OutputImageType::PixelType PixelType; + typedef itk::Image WeightsImageType; - //========================================================================================================================= - //Member functions of the helper class 2 - //========================================================================================================================= - - - //========================================================================================================================= - //Empty constructor - template HelperClass2::HelperClass2() - { - m_EdgePaddingValue=itk::NumericTraits::Zero; + //Set methods + void SetWeights(const typename WeightsImageType::Pointer input) { + m_Weights = input; + this->Modified(); + } + void SetEdgePaddingValue(PixelType value) { + m_EdgePaddingValue = value; + this->Modified(); } - - - //========================================================================================================================= - //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 - PixelType 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; - - //Neighbour cannot have zero weight because it should be filled already - neighValue=outputIt.Get(); - ++outputIt; - outputIt.Set(neighValue); - //DD("hole filled"); - } - else{ - //DD("is at begin, setting edgepadding value"); - outputIt.Set(m_EdgePaddingValue); - } - } - ++weightsIt; - ++outputIt; - ++inputIt; - - }//end while - }//end member + /** Typedef to describe the output image region type. */ + typedef typename OutputImageType::RegionType OutputImageRegionType; -}//end nameless namespace +protected: + HelperClass2(); + ~HelperClass2() {}; + //the actual processing + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ); -namespace clitk + + //member data + typename WeightsImageType::Pointer m_Weights; + PixelType m_EdgePaddingValue; + +} ; + + + +//========================================================================================================================= +//Member functions of the helper class 2 +//========================================================================================================================= + + +//========================================================================================================================= +//Empty constructor +template HelperClass2::HelperClass2() +{ + m_EdgePaddingValue=itk::NumericTraits::Zero; +} + + +//========================================================================================================================= +//update the output for the outputRegionForThread +template void HelperClass2::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ) { - //========================================================================================================================= - // The rest is the InvertVFFilter - //========================================================================================================================= - - //========================================================================================================================= - //constructor - template - InvertVFFilter::InvertVFFilter() - { - m_EdgePaddingValue=itk::NumericTraits::Zero; //no other reasonable value? - m_ThreadSafe=false; - m_Verbose=false; - } + //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 + PixelType 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())); + } - //========================================================================================================================= - //Update - template void InvertVFFilter::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); - PixelType zero = itk::NumericTraits::Zero; - - - //Allocate the weights - typename WeightsImageType::Pointer weights=WeightsImageType::New(); - weights->SetRegions(region); - weights->Allocate(); - weights->SetSpacing(inputPtr->GetSpacing()); - - //=========================================================================== - //Inversion is divided in in two loops, for each we will 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 - - typedef itk::Image, ImageDimension > 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->SetWeights(weights); - - //Threadsafe? - if(m_ThreadSafe) - { - //Allocate the mutex image - typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New(); - mutex->SetRegions(region); - mutex->Allocate(); - mutex->SetSpacing(inputPtr->GetSpacing()); - helper1->SetMutexImage(mutex); - if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <Update(); + }//end while +}//end member - //Get the output - typename InternalImageType::Pointer temp= helper1->GetOutput(); - weights=helper1->GetWeights(); +}//end nameless namespace + + + +namespace clitk +{ + +//========================================================================================================================= +// The rest is the InvertVFFilter +//========================================================================================================================= + +//========================================================================================================================= +//constructor +template +InvertVFFilter::InvertVFFilter() +{ + m_EdgePaddingValue=itk::NumericTraits::Zero; //no other reasonable value? + m_ThreadSafe=false; + m_Verbose=false; +} + + +//========================================================================================================================= +//Update +template void InvertVFFilter::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); + PixelType zero = itk::NumericTraits::Zero; + + + //Allocate the weights + typename WeightsImageType::Pointer weights=WeightsImageType::New(); + weights->SetRegions(region); + weights->Allocate(); + weights->SetSpacing(inputPtr->GetSpacing()); + + //=========================================================================== + //Inversion is divided in in two loops, for each we will 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 + + typedef itk::Image, ImageDimension > 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->SetWeights(weights); + + //Threadsafe? + if(m_ThreadSafe) { + //Allocate the mutex image + typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New(); + mutex->SetRegions(region); + mutex->Allocate(); + mutex->SetSpacing(inputPtr->GetSpacing()); + helper1->SetMutexImage(mutex); + if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <Update(); + + //Get the output + typename InternalImageType::Pointer temp= helper1->GetOutput(); + 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 + helper2->SetInput(temp); + helper2->SetWeights(weights); + helper2->SetEdgePaddingValue(m_EdgePaddingValue); + + //Execute helper class + if (m_Verbose) std::cout << "Normalizing the output VF..."<Update(); + + //Set the output + this->SetNthOutput(0, helper2->GetOutput()); +} - //=========================================================================== - //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 - helper2->SetInput(temp); - helper2->SetWeights(weights); - helper2->SetEdgePaddingValue(m_EdgePaddingValue); - //Execute helper class - if (m_Verbose) std::cout << "Normalizing the output VF..."<Update(); - //Set the output - this->SetNthOutput(0, helper2->GetOutput()); - } - - - } #endif