X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkBSplineDeformableTransform.txx;h=db47430b5a82ab339d4103558a29d923bce7d6c8;hb=d64e6b39ec8d4bb7b99f9d1a531e9dba3121e2bc;hp=4290191679d182a02d1f66f1b8bf19c01f912f9c;hpb=657652a78c2e2717a6f77e027049173442ca29f0;p=clitk.git

diff --git a/registration/clitkBSplineDeformableTransform.txx b/registration/clitkBSplineDeformableTransform.txx
index 4290191..db47430 100644
--- a/registration/clitkBSplineDeformableTransform.txx
+++ b/registration/clitkBSplineDeformableTransform.txx
@@ -3,7 +3,7 @@
 
   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 __clitkBSplineDeformableTransform_txx
 #define __clitkBSplineDeformableTransform_tx
 #include "clitkBSplineDeformableTransform.h"
@@ -86,14 +86,7 @@ namespace clitk
     for (i=0; i <InputDimension;i++)
       {
 	m_Offset[i] = m_SplineOrders[i] / 2;
-	if ( m_SplineOrders[i] % 2 ) 
-	  {
-	    m_SplineOrderOdd[i] = true;
-	  }
-	else
-	  {
-	    m_SplineOrderOdd[i] = false;
-	  }
+        m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
       }
     m_ValidRegion = m_GridRegion;
 	
@@ -172,7 +165,6 @@ namespace clitk
     m_SupportRegion.SetIndex(m_SupportIndex);
     for (  i = 0; i < OutputDimension; i++ ) 
       m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
-  
   }
     
 
@@ -204,28 +196,20 @@ namespace clitk
   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
   ::SetSplineOrders(const SizeType & splineOrders)
   {
-    if(m_SplineOrders!=splineOrders)
+    if (m_SplineOrders!=splineOrders)
       {
 	m_SplineOrders=splineOrders;
 	
 	//update the interpolation function
 	m_VectorInterpolator->SetSplineOrders(m_SplineOrders);
 	
-	//update the varaibles for computing interpolation
+	//update the variables for computing interpolation
 	for (unsigned int i=0; i <InputDimension;i++)
-	  {
-	    m_SupportSize[i] = m_SplineOrders[i]+1;
-	    m_Offset[i] = m_SplineOrders[i] / 2;
-	    
-	    if ( m_SplineOrders[i] % 2 ) 
-	      {
-		m_SplineOrderOdd[i] = true;
-	      }
-	    else
-	      {
-		m_SplineOrderOdd[i] = false;
-	      }
-	  }
+        {
+          m_SupportSize[i] = m_SplineOrders[i]+1;
+          m_Offset[i] = m_SplineOrders[i] / 2;
+          m_SplineOrderOdd[i] = m_SplineOrders[i] % 2;
+        }
 	this->Modified();
       }
   }
@@ -271,7 +255,6 @@ namespace clitk
   BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
   ::GetNumberOfParameters(void) const
   {
-
     // The number of parameters equal OutputDimension * number of
     // of pixels in the grid region.
     return ( static_cast<unsigned int>( OutputDimension ) *
@@ -307,7 +290,7 @@ namespace clitk
 	m_WrappedImage->SetRegions( m_GridRegion );
 	for (unsigned int j=0; j <OutputDimension;j++)
 	  m_JacobianImage[j]->SetRegions( m_GridRegion );
-	
+
 	// Set the valid region
 	// If the grid spans the interval [start, last].
 	// The valid interval for evaluation is [start+offset, last-offset]
@@ -448,7 +431,6 @@ namespace clitk
       }
   }
 
-
   // Set the parameters
   template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
   void
@@ -600,6 +582,7 @@ namespace clitk
     JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
     memset(jacobianDataPointer, 0,  OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
     m_LastJacobianIndex = m_ValidRegion.GetIndex();
+    m_NeedResetJacobian = false;
 
     for (unsigned int j=0; j<OutputDimension; j++)
       {
@@ -935,33 +918,9 @@ namespace clitk
     // 	itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" );
     //       }
 
-  
-    //========================================================
-    // Zero all components of jacobian
-    //========================================================
-    // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
-    // NOTE: for efficiency, we only need to zero out the coefficients
-    // that got fill last time this method was called.
+    if (m_NeedResetJacobian)
+      ResetJacobian();
 
-    unsigned int j=0;
-      
-    //Define the  region for each jacobian image 
-    m_SupportRegion.SetIndex( m_LastJacobianIndex );
-    
-    //Initialize the iterators
-    for ( j = 0; j < OutputDimension; j++ ) 
-      m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
-    
-    //Set the previously-set to zero 
-    while ( ! (m_Iterator[0]).IsAtEnd() )
-      {
-	for ( j = 0; j < OutputDimension; j++ )
-	  {
-	    m_Iterator[j].Set( m_ZeroVector );
-	    ++(m_Iterator[j]);
-	  }
-      }
-    
     //========================================================
     // For each dimension, copy the weight to the support region
     //========================================================
@@ -989,6 +948,7 @@ namespace clitk
     m_SupportRegion.SetIndex( m_LastJacobianIndex );
 
     //Reset the iterators
+    unsigned int j = 0;
     for ( j = 0; j < OutputDimension; j++ ) 
       m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion);
 
@@ -1006,6 +966,7 @@ namespace clitk
 	// go to next coefficient in the support region
 	weights++;
       }
+    m_NeedResetJacobian = true;
 
     // Return the result
     return this->m_Jacobian;
@@ -1037,7 +998,48 @@ namespace clitk
       }
   }
 
+  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+  unsigned
+  BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
+  ::SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim)
+  {
+    unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
+    m_JacobianImage[dim]->GetPixelContainer()->SetImportPointer(jacobianDataPointer, numberOfPixels);
+    return numberOfPixels;
+  }
 
+  template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+  void
+  BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
+  ::ResetJacobian() const
+  {
+    //========================================================
+    // Zero all components of jacobian
+    //========================================================
+    // JV not thread safe (m_LastJacobianIndex), instantiate N transforms
+    // NOTE: for efficiency, we only need to zero out the coefficients
+    // that got fill last time this method was called.
+
+    unsigned int j = 0;
+
+    //Define the  region for each jacobian image
+    m_SupportRegion.SetIndex(m_LastJacobianIndex);
+
+    //Initialize the iterators
+    for (j = 0; j < OutputDimension; j++)
+      m_Iterator[j] = IteratorType(m_JacobianImage[j], m_SupportRegion);
+
+    //Set the previously-set to zero
+    while (!(m_Iterator[0]).IsAtEnd())
+      {
+	for (j = 0; j < OutputDimension; j++)
+	  {
+	    m_Iterator[j].Set(m_ZeroVector);
+	    ++(m_Iterator[j]);
+	  }
+      }
+    m_NeedResetJacobian = false;
+  }
 } // namespace
 
 #endif