invert VF on given image grid

- with --like option

diff --git a/itk/clitkInvertVFFilter.h b/itk/clitkInvertVFFilter.h
index 2afceb6992919a5ffd7a74506a2313421ebc4325..586d872fd619ce4351ace1dbe4979fd760d42d41 100644 (file)
--- a/itk/clitkInvertVFFilter.h
+++ b/itk/clitkInvertVFFilter.h
@@ -71,15 +71,21 @@ namespace clitk
       m_NumberOfThreads=r;
     }
     itkSetMacro(ThreadSafe, bool);
+    itkSetMacro(OutputSpacing, SpacingType);
+    itkSetMacro(OutputSize, SizeType);
 
   
   protected:
     InvertVFFilter();
     ~InvertVFFilter() {};
     void GenerateData( );
-   
+    void GenerateOutputInformation();
+    void GenerateInputRequestedRegion();
+    
     bool m_Verbose;
     bool m_NumberOfThreadsIsGiven;
+    SpacingType m_OutputSpacing;
+    SizeType m_OutputSize;
     unsigned int m_NumberOfThreads;
     PixelType m_EdgePaddingValue;
     bool m_ThreadSafe;

diff --git a/itk/clitkInvertVFFilter.txx b/itk/clitkInvertVFFilter.txx
index 0c443ca8e4812764371fda4520cb3f0cbe73b021..ae6be89f8e85c424dbbc9c7a2c39bbc5f0947da3 100644 (file)
--- a/itk/clitkInvertVFFilter.txx
+++ b/itk/clitkInvertVFFilter.txx
@@ -17,6 +17,7 @@
 ===========================================================================**/
 #ifndef __clitkInvertVFFilter_txx
 #define __clitkInvertVFFilter_txx
+
 namespace
 {
 
@@ -73,6 +74,8 @@ protected:
   ~HelperClass1() {};
 
   //the actual processing
+  void GenerateInputRequestedRegion();
+  void GenerateOutputInformation();  
   void BeforeThreadedGenerateData();
   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId );
 
@@ -98,11 +101,39 @@ HelperClass1<InputImageType, OutputImageType>::HelperClass1()
   m_ThreadSafe=false;
 }
 
+//=========================================================================================================================
+//Set input image params
+template <class InputImageType, class OutputImageType>
+void HelperClass1<InputImageType, OutputImageType>::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<class InputImageType, class OutputImageType >
+void HelperClass1<InputImageType, OutputImageType>::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<class InputImageType, class OutputImageType >
 void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
 {
+  //std::cout << "HelperClass1::BeforeThreadedGenerateData - IN" << std::endl;
   //Since we will add, put to zero!
   this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
   this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
@@ -113,49 +144,58 @@ void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
 template<class InputImageType, class OutputImageType>
 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
 {
-
+  //std::cout << "HelperClass1::ThreadedGenerateData - IN" << std::endl;
   //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();
+/*  outputPtr->SetLargestPossibleRegion( m_Weights->GetLargestPossibleRegion() );
+  outputPtr->SetSpacing(m_Weights->GetSpacing());*/
+  typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
+  //std::cout << outputPtr->GetSpacing() << size << std::endl;
 
   //Iterator over input
-  typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
+  //typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
+  typedef itk::ImageRegionConstIteratorWithIndex<OutputImageType> InputImageIteratorType;
 
   //define them over the outputRegionForThread
-  InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
+  //InputImageIteratorType outputIt(inputPtr, outputRegionForThread);
+  InputImageIteratorType outputIt(outputPtr, outputPtr->GetLargestPossibleRegion());
 
   //Initialize
   typename InputImageType::IndexType index;
-  itk::ContinuousIndex<double,ImageDimension> contIndex;
+  itk::ContinuousIndex<double,ImageDimension> contIndex, inContIndex;
   typename InputImageType::PointType ipoint;
   typename OutputImageType::PointType opoint;
   typedef typename OutputImageType::PixelType DisplacementType;
   DisplacementType displacement;
-  inputIt.GoToBegin();
+  outputIt.GoToBegin();
 
   //define some temp variables
   signed long baseIndex[ImageDimension];
   double distance[ImageDimension];
   unsigned int dim, counter, upper;
   double totalOverlap,overlap;
-  typename OutputImageType::IndexType neighIndex;
+  typename OutputImageType::IndexType neighIndex, outIndex;
 
   //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
+  //Loop over the output region and add the intensities from the input to the output and the weight to the weights
   //==================================================================================================
-  while( !inputIt.IsAtEnd() ) {
+  while( !outputIt.IsAtEnd() ) {
     // get the input image index
-    index = inputIt.GetIndex();
+    outIndex = outputIt.GetIndex();
+    outputPtr->TransformIndexToPhysicalPoint( outIndex,opoint );
+    for(unsigned int j = 0; j < ImageDimension; j++ ) ipoint[j] = opoint[j];
+    inputPtr->TransformPhysicalPointToIndex(ipoint, index);
     inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
 
     // get the required displacement
-    displacement = inputIt.Get();
+    //displacement = outputIt.Get();
+    displacement = inputPtr->GetPixel(index);
 
     // compute the required output image point
     for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
@@ -230,9 +270,10 @@ void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const O
       }
     }
 
-    ++inputIt;
+    ++outputIt;
   }
 
+  //std::cout << "HelperClass1::ThreadedGenerateData - OUT" << std::endl;
 }
 
 
@@ -310,7 +351,8 @@ template<class InputImageType, class OutputImageType > HelperClass2<InputImageTy
 //update the output for the outputRegionForThread
 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
 {
-
+  //std::cout << "HelperClass2::ThreadedGenerateData - IN" << std::endl;
+  
   //Get pointer to the input
   typename InputImageType::ConstPointer inputPtr = this->GetInput();
 
@@ -360,6 +402,9 @@ template<class InputImageType, class OutputImageType > void HelperClass2<InputIm
     ++inputIt;
 
   }//end while
+  
+  //std::cout << "HelperClass2::ThreadedGenerateData - OUT" << std::endl;
+  
 }//end member
 
 
@@ -384,27 +429,60 @@ InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
   m_Verbose=false;
 }
 
+//=========================================================================================================================
+//Set input image params
+template <class InputImageType, class OutputImageType>
+void InvertVFFilter<InputImageType, OutputImageType>::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<class InputImageType, class OutputImageType >
+void InvertVFFilter<InputImageType, OutputImageType>::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 <class InputImageType, class OutputImageType> void InvertVFFilter<InputImageType, OutputImageType>::GenerateData()
 {
+  //std::cout << "InvertVFFilter::GenerateData - IN" << std::endl;
 
   //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);
-
+  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
@@ -435,7 +513,8 @@ template <class InputImageType, class OutputImageType> void InvertVFFilter<Input
     typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New();
     mutex->SetRegions(region);
     mutex->Allocate();
-    mutex->SetSpacing(inputPtr->GetSpacing());
+    //mutex->SetSpacing(inputPtr->GetSpacing());
+    mutex->SetSpacing(m_OutputSpacing);
     helper1->SetMutexImage(mutex);
     if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <<std::endl;
   } else  if(m_Verbose)std::cout <<"Inverting using a thread-unsafe algorithm" <<std::endl;
@@ -455,6 +534,8 @@ template <class InputImageType, class OutputImageType> void InvertVFFilter<Input
   typename HelperClass2Type::Pointer helper2=HelperClass2Type::New();
 
   //Set temporary output as input
+  //std::cout << temp->GetSpacing() << 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/clitkInvertVF.ggo b/tools/clitkInvertVF.ggo
index 54c1e21ab101d3b9c5f32787803812a9df476ba2..19c451a478bba22e6a86b204eb4a6627c17f961e 100644 (file)
--- a/tools/clitkInvertVF.ggo
+++ b/tools/clitkInvertVF.ggo
@@ -15,3 +15,4 @@ option "type"                 t       "Type of filter: 0=clitk (fast forward splat using line
 option "threadSafe"                    -       "Clitk: use thread safe algorithm"                              flag            off
 option "pad"                   p       "Clitk: edge padding value (1 or N number of values!, defautls to zeros)"               double          multiple        no
 option "sampling"              s       "Itk: subsampling factor"                                       int             no      default="20"
+option "like"   l "Image whose grid (spacing and size) will be used for output"         string    no

diff --git a/tools/clitkInvertVFGenericFilter.txx b/tools/clitkInvertVFGenericFilter.txx
index 9835c2e739f64d3edd1af378e7955efc43059ba7..2ddd46c5b4d6cc51422efd9fb17e5f6d97a4a0a7 100644 (file)
--- a/tools/clitkInvertVFGenericFilter.txx
+++ b/tools/clitkInvertVFGenericFilter.txx
@@ -129,6 +129,19 @@ InvertVFGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
     typedef clitk::InvertVFFilter<InputImageType,OutputImageType> 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) {
+      itk::ImageIOBase::Pointer header = readImageHeader(m_ArgsInfo.like_arg);
+      for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
+        size[i] = header->GetDimensions(i);
+        spacing[i] = header->GetSpacing(i);
+      }
+    }
+    std::cout << spacing << size << std::endl;
+    filter->SetOutputSpacing(spacing);
+    filter->SetOutputSize(size);
+
     filter->SetVerbose(m_Verbose);
     if (m_ArgsInfo.threads_given) filter->SetNumberOfThreads(m_ArgsInfo.threads_arg);
     if (m_ArgsInfo.pad_given) {