X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkBSplineDeformableTransformInitializer.txx;fp=registration%2FclitkBSplineDeformableTransformInitializer.txx;h=0a93ab0432ae814347454e5b47121e016e8bbedd;hb=657652a78c2e2717a6f77e027049173442ca29f0;hp=0000000000000000000000000000000000000000;hpb=8abe9edbdc63a33f8fbb1f321073b762e08be9aa;p=clitk.git

diff --git a/registration/clitkBSplineDeformableTransformInitializer.txx b/registration/clitkBSplineDeformableTransformInitializer.txx
new file mode 100644
index 0000000..0a93ab0
--- /dev/null
+++ b/registration/clitkBSplineDeformableTransformInitializer.txx
@@ -0,0 +1,244 @@
+/*=========================================================================
+  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://oncora1.lyon.fnclcc.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 __clitkBSplineDeformableTransformInitializer_txx
+#define __clitkBSplineDeformableTransformInitializer_txx
+#include "clitkBSplineDeformableTransformInitializer.h"
+
+namespace clitk
+{
+
+
+  template < class TTransform, class TImage >
+  BSplineDeformableTransformInitializer<TTransform, TImage >
+  ::BSplineDeformableTransformInitializer()
+  {
+    this->m_NumberOfControlPointsInsideTheImage.Fill( 5 );
+    this->m_ControlPointSpacing.Fill(64.);
+    this->m_ChosenSpacing.Fill(64.);
+    m_NumberOfControlPointsIsGiven=false;
+    m_ControlPointSpacingIsGiven=false;
+    m_TransformRegionIsGiven=false;
+    m_SamplingFactorIsGiven=false;
+    m_SplineOrders.Fill(3);
+    m_Parameters=NULL;
+  }
+
+  template < class TTransform, class TImage >
+  void
+  BSplineDeformableTransformInitializer<TTransform, TImage >
+  ::InitializeTransform()
+  {
+    // Sanity check:
+    // The image is required for the region and spacing
+    if( ! this->m_Image )
+      {
+	itkExceptionMacro( "Reference Image has not been set" );
+	return;
+      }
+    
+    // A pointer to the transform is required
+    if( ! this->m_Transform )
+      {
+	itkExceptionMacro( "Transform has not been set" );
+	return;
+      }
+
+    // If the image come from a filter, then update that filter.
+    if( this->m_Image->GetSource() )
+      {
+	this->m_Image->GetSource()->Update();
+      }
+
+
+    //--------------------------------------
+    // Spacing & Size on image
+    //--------------------------------------
+    SpacingType fixedImageSpacing = m_Image->GetSpacing();
+    SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize();
+    typename RegionType::SizeType   gridBorderSize;
+    typename RegionType::SizeType   totalGridSize;
+
+    // Only spacing given: adjust if necessary
+    if (m_ControlPointSpacingIsGiven  && !m_NumberOfControlPointsIsGiven)
+      {
+	for(unsigned int r=0; r<InputDimension; r++) 
+	  {
+	    // JV
+	    m_ChosenSpacing[r]= m_ControlPointSpacing[r];
+	    m_ControlPointSpacing[r]= ( round(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
+	    m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
+	    if (  ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
+		 == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
+	      {
+		m_NumberOfControlPointsInsideTheImage[r] +=1;
+		DD("Adding control point");
+	      }
+	  }
+      }
+  
+    // Only number of CP given: adjust if necessary
+    // JV -1 ?
+    else if (!m_ControlPointSpacingIsGiven  && m_NumberOfControlPointsIsGiven)
+      {
+	for(unsigned int r=0; r<InputDimension; r++) 
+	  {
+	    m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r])  / 
+						      (double)(m_NumberOfControlPointsInsideTheImage[r]) );
+	    m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1)  / 
+								  (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
+	  }
+      }
+  
+    // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
+    else 
+      {
+	for(unsigned int r=0; r<InputDimension; r++) 
+	  {
+	    m_ChosenSpacing[r]= m_ControlPointSpacing[r];
+	    if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r]) 
+	      {
+		std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
+			 <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
+			 <<"mm) for dimension "<<r<<"!" <<std::endl
+			 <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+	      }
+	    if (  fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) ) 
+	      {
+		std::cout<<"WARNING: Specified control point spacing for dimension "<<r
+			 <<" does not allow exact representation of BLUT FFD!"<<std::endl
+			 <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
+			 <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;  
+	      }
+	  }
+      }
+
+    if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
+    if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl; 
+    if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl; 
+    
+
+    //--------------------------------------
+    // Border size 
+    //--------------------------------------
+    for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
+    totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
+    if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
+    RegionType gridRegion;
+    gridRegion.SetSize(totalGridSize);
+ 
+  
+    //--------------------------------------
+    // Origin
+    //--------------------------------------
+    typedef typename TransformType::OriginType OriginType;
+    OriginType fixedImageOrigin, gridOrigin;
+
+    // In case of non-zero index
+    m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
+    typename ImageType::DirectionType gridDirection = m_Image->GetDirection();      
+    SizeType orderShift;
+
+    // Shift depends on order
+    for(unsigned int r=0; r<InputDimension; r++)
+      {
+	orderShift[r] = m_SplineOrders[r] / 2;
+	gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
+      } 
+    if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
+
+
+    //--------------------------------------
+    // LUT sampling factors
+    //--------------------------------------
+    for (unsigned int i=0; i< InputDimension; i++)
+      {
+	if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
+				<< m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
+	if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
+	if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<m_SamplingFactors[i]<<"..."<<std::endl;
+      }
+
+
+    //--------------------------------------
+    // Set
+    //--------------------------------------
+    this->m_Transform->SetSplineOrders(m_SplineOrders);
+    this->m_Transform->SetGridRegion( gridRegion );
+    this->m_Transform->SetGridOrigin( gridOrigin );
+    this->m_Transform->SetGridSpacing( m_ControlPointSpacing );
+    this->m_Transform->SetGridDirection( gridDirection );
+    this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors);
+
+  }
+
+  template < class TTransform, class TImage >
+  void
+  BSplineDeformableTransformInitializer<TTransform, TImage >
+  ::SetInitialParameters( const std::string& s, ParametersType& params) 
+  {
+    //---------------------------------------
+    // Read Initial parameters
+    //---------------------------------------
+    typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
+    typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
+    coeffReader->SetFileName(s);
+    if(m_Verbose) std::cout<<"Reading initial coefficients from file: "<<s<<"..."<<std::endl;
+    coeffReader->Update();
+    typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
+    this->SetInitialParameters(coefficientImage, params);
+  }
+  
+  template < class TTransform, class TImage >
+  void
+  BSplineDeformableTransformInitializer<TTransform, TImage >
+  ::SetInitialParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params) 
+  {
+    // Keep a reference 
+    m_Parameters=&params;
+
+    // Resample
+    typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
+    typename ResamplerType::Pointer resampler=ResamplerType::New();
+    resampler->SetInput(coefficientImage);	
+    resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage()); 
+    if(m_Verbose) std::cout<<"Resampling initial coefficients..."<<std::endl;
+    resampler->Update();
+
+    // Copy parameters into the existing array, I know its crappy
+    typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
+    Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
+    it.GoToBegin();
+    unsigned int i=0;
+    while(! it.IsAtEnd())
+      {
+	for (unsigned int j=0; j<InputDimension;j++)
+	  params[i+j]=it.Get()[j];
+
+	++it;
+	i+=InputDimension;
+      }
+
+    // JV pass the reference to the array !!
+    this->m_Transform->SetParameters(params);
+  }
+
+
+  
+}  // namespace itk
+
+#endif