From 42f68068ef3488c8ca2d921b6835d34fe6494d5f Mon Sep 17 00:00:00 2001 From: Romulo Pinho Date: Tue, 10 Jan 2012 16:36:21 +0100 Subject: [PATCH] invert Vf on given image grid with resampler - linear interpolation (itk resampler's default) - undoes previous version --- itk/clitkInvertVFFilter.h | 2 - itk/clitkInvertVFFilter.txx | 101 +++------------------------ tools/clitkInvertVFGenericFilter.txx | 26 +++++-- 3 files changed, 30 insertions(+), 99 deletions(-) diff --git a/itk/clitkInvertVFFilter.h b/itk/clitkInvertVFFilter.h index 586d872..8475241 100644 --- a/itk/clitkInvertVFFilter.h +++ b/itk/clitkInvertVFFilter.h @@ -79,8 +79,6 @@ namespace clitk InvertVFFilter(); ~InvertVFFilter() {}; void GenerateData( ); - void GenerateOutputInformation(); - void GenerateInputRequestedRegion(); bool m_Verbose; bool m_NumberOfThreadsIsGiven; diff --git a/itk/clitkInvertVFFilter.txx b/itk/clitkInvertVFFilter.txx index ae6be89..684b355 100644 --- a/itk/clitkInvertVFFilter.txx +++ b/itk/clitkInvertVFFilter.txx @@ -74,8 +74,6 @@ protected: ~HelperClass1() {}; //the actual processing - void GenerateInputRequestedRegion(); - void GenerateOutputInformation(); void BeforeThreadedGenerateData(); void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); @@ -101,33 +99,6 @@ HelperClass1::HelperClass1() m_ThreadSafe=false; } -//========================================================================================================================= -//Set input image params -template -void HelperClass1::GenerateInputRequestedRegion() -{ - //std::cout << "HelperClass1::GenerateInputRequestedRegion - IN" << std::endl; - - typename InputImageType::Pointer inputPtr = const_cast< InputImageType * >(this->GetInput()); - inputPtr->SetRequestedRegion(inputPtr->GetLargestPossibleRegion()); - - //std::cout << "HelperClass1::GenerateInputRequestedRegion - OUT" << std::endl; -} - -//========================================================================================================================= -//Set output image params -template -void HelperClass1::GenerateOutputInformation() -{ - //std::cout << "HelperClass1::GenerateOutputInformation - IN" << std::endl; - Superclass::GenerateOutputInformation(); - - // it's the weight image that determines the size of the output - typename OutputImageType::Pointer outputPtr = this->GetOutput(); - outputPtr->SetLargestPossibleRegion( m_Weights->GetLargestPossibleRegion() ); - outputPtr->SetSpacing(m_Weights->GetSpacing()); -} - //========================================================================================================================= //Before threaded data template @@ -150,18 +121,13 @@ void HelperClass1::ThreadedGenerateData(const O //Get pointer to the output typename OutputImageType::Pointer outputPtr = this->GetOutput(); -/* outputPtr->SetLargestPossibleRegion( m_Weights->GetLargestPossibleRegion() ); - outputPtr->SetSpacing(m_Weights->GetSpacing());*/ - typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize(); - //std::cout << outputPtr->GetSpacing() << size << std::endl; + //typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize(); //Iterator over input - //typedef itk::ImageRegionConstIteratorWithIndex InputImageIteratorType; - typedef itk::ImageRegionConstIteratorWithIndex InputImageIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex InputImageIteratorType; //define them over the outputRegionForThread - //InputImageIteratorType outputIt(inputPtr, outputRegionForThread); - InputImageIteratorType outputIt(outputPtr, outputPtr->GetLargestPossibleRegion()); + InputImageIteratorType inputIt(inputPtr, outputRegionForThread); //Initialize typename InputImageType::IndexType index; @@ -170,14 +136,14 @@ void HelperClass1::ThreadedGenerateData(const O typename OutputImageType::PointType opoint; typedef typename OutputImageType::PixelType DisplacementType; DisplacementType displacement; - outputIt.GoToBegin(); + 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, outIndex; + typename OutputImageType::IndexType neighIndex; //Find the number of neighbors unsigned int neighbors = 1 << ImageDimension; @@ -185,17 +151,13 @@ void HelperClass1::ThreadedGenerateData(const O //================================================================================================== //Loop over the output region and add the intensities from the input to the output and the weight to the weights //================================================================================================== - while( !outputIt.IsAtEnd() ) { + while( !inputIt.IsAtEnd() ) { // get the input image index - outIndex = outputIt.GetIndex(); - outputPtr->TransformIndexToPhysicalPoint( outIndex,opoint ); - for(unsigned int j = 0; j < ImageDimension; j++ ) ipoint[j] = opoint[j]; - inputPtr->TransformPhysicalPointToIndex(ipoint, index); + index = inputIt.GetIndex(); inputPtr->TransformIndexToPhysicalPoint( index,ipoint ); // get the required displacement - //displacement = outputIt.Get(); - displacement = inputPtr->GetPixel(index); + displacement = inputIt.Get(); // compute the required output image point for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j]; @@ -270,7 +232,7 @@ void HelperClass1::ThreadedGenerateData(const O } } - ++outputIt; + ++inputIt; } //std::cout << "HelperClass1::ThreadedGenerateData - OUT" << std::endl; @@ -429,41 +391,6 @@ InvertVFFilter::InvertVFFilter() m_Verbose=false; } -//========================================================================================================================= -//Set input image params -template -void InvertVFFilter::GenerateInputRequestedRegion() -{ - //std::cout << "InvertVFFilter::GenerateInputRequestedRegion - IN" << std::endl; - // call the superclass' implementation of this method - // Superclass::GenerateInputRequestedRegion(); - - typename InputImageType::Pointer inputPtr = const_cast< InputImageType * >(this->GetInput()); - inputPtr->SetRequestedRegion(inputPtr->GetLargestPossibleRegion()); - - //std::cout << "InvertVFFilter::GenerateInputRequestedRegion - OUT" << std::endl; -} - -//========================================================================================================================= -//Set output image params -template -void InvertVFFilter::GenerateOutputInformation() -{ - //std::cout << "InvertVFFilter::GenerateOutputInformation - IN" << std::endl; - Superclass::GenerateOutputInformation(); - - typename InputImageType::ConstPointer inputPtr=this->GetInput(); - typename WeightsImageType::RegionType region = inputPtr->GetLargestPossibleRegion(); - region.SetSize(m_OutputSize); - - typename OutputImageType::Pointer outputPtr = this->GetOutput(); - outputPtr->SetRegions( region ); - outputPtr->SetSpacing(m_OutputSpacing); - outputPtr->SetOrigin(inputPtr->GetOrigin()); - - //std::cout << "InvertVFFilter::GenerateOutputInformation - OUT" << std::endl; -} - //========================================================================================================================= //Update template void InvertVFFilter::GenerateData() @@ -473,16 +400,13 @@ template void InvertVFFilterGetInput(); typename WeightsImageType::RegionType region = inputPtr->GetLargestPossibleRegion(); - //region.SetSize(size); - region.SetSize(m_OutputSize); //Allocate the weights typename WeightsImageType::Pointer weights=WeightsImageType::New(); weights->SetOrigin(inputPtr->GetOrigin()); weights->SetRegions(region); - weights->SetSpacing(m_OutputSpacing); weights->Allocate(); - //weights->SetSpacing(inputPtr->GetSpacing()); + weights->SetSpacing(inputPtr->GetSpacing()); //=========================================================================== //Inversion is divided in in two loops, for each we will call a threaded helper class @@ -513,8 +437,7 @@ template void InvertVFFilterSetRegions(region); mutex->Allocate(); - //mutex->SetSpacing(inputPtr->GetSpacing()); - mutex->SetSpacing(m_OutputSpacing); + mutex->SetSpacing(inputPtr->GetSpacing()); helper1->SetMutexImage(mutex); if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" < void InvertVFFilterGetSpacing() << temp->GetLargestPossibleRegion().GetSize() << std::endl; - //std::cout << weights->GetSpacing() << weights->GetLargestPossibleRegion().GetSize() << std::endl; helper2->SetInput(temp); helper2->SetWeights(weights); helper2->SetEdgePaddingValue(m_EdgePaddingValue); diff --git a/tools/clitkInvertVFGenericFilter.txx b/tools/clitkInvertVFGenericFilter.txx index 2ddd46c..558b341 100644 --- a/tools/clitkInvertVFGenericFilter.txx +++ b/tools/clitkInvertVFGenericFilter.txx @@ -18,6 +18,8 @@ #ifndef clitkInvertVFGenericFilter_txx #define clitkInvertVFGenericFilter_txx +#include "itkVectorResampleImageFilter.h" + /* ================================================= * @file clitkInvertVFGenericFilter.txx * @author @@ -117,7 +119,7 @@ InvertVFGenericFilter::UpdateWithDimAndPixelType() typename InputReaderType::Pointer reader = InputReaderType::New(); reader->SetFileName( m_InputFileName); reader->Update(); - typename InputImageType::Pointer input= reader->GetOutput(); + typename InputImageType::Pointer input = reader->GetOutput(); // Filter typename OutputImageType::Pointer output; @@ -128,19 +130,29 @@ InvertVFGenericFilter::UpdateWithDimAndPixelType() // Create the InvertVFFilter typedef clitk::InvertVFFilter FilterType; typename FilterType::Pointer filter =FilterType::New(); - filter->SetInput(input); - typename FilterType::SpacingType spacing = input->GetSpacing(); - typename FilterType::SizeType size = input->GetLargestPossibleRegion().GetSize(); if (m_ArgsInfo.like_given) { + typename FilterType::SpacingType spacing; + typename FilterType::SizeType size; itk::ImageIOBase::Pointer header = readImageHeader(m_ArgsInfo.like_arg); for(unsigned int i=0; iGetDimensions(i); spacing[i] = header->GetSpacing(i); } + + typedef itk::VectorResampleImageFilter ResampleFilterType; + typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); + resampler->SetInput(input); + resampler->SetOutputOrigin(input->GetOrigin()); + resampler->SetOutputDirection(input->GetDirection()); + resampler->SetOutputSpacing(spacing); + resampler->SetSize(size); + resampler->Update(); + spacing = resampler->GetOutput()->GetSpacing(); + size = resampler->GetOutput()->GetLargestPossibleRegion().GetSize(); + filter->SetInput(resampler->GetOutput()); } - std::cout << spacing << size << std::endl; - filter->SetOutputSpacing(spacing); - filter->SetOutputSize(size); + else + filter->SetInput(input); filter->SetVerbose(m_Verbose); if (m_ArgsInfo.threads_given) filter->SetNumberOfThreads(m_ArgsInfo.threads_arg); -- 2.45.1