X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkBSplineDeformableTransform.txx;h=60a2136d7160c9bbfce3431e376bd2d9232dd08c;hb=3504c907c279d7f51f0e1d81eead7243135794a5;hp=4290191679d182a02d1f66f1b8bf19c01f912f9c;hpb=657652a78c2e2717a6f77e027049173442ca29f0;p=clitk.git diff --git a/registration/clitkBSplineDeformableTransform.txx b/registration/clitkBSplineDeformableTransform.txx index 4290191..60a2136 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,9 +14,9 @@ - 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 +#define __clitkBSplineDeformableTransform_txx #include "clitkBSplineDeformableTransform.h" //itk @@ -31,8 +31,7 @@ namespace clitk // Constructor with default arguments template - BSplineDeformableTransform - ::BSplineDeformableTransform():Superclass(OutputDimension,0) + BSplineDeformableTransform::BSplineDeformableTransform():Superclass(0) { unsigned int i; @@ -86,14 +85,7 @@ namespace clitk for (i=0; i ::Zero; - } @@ -204,28 +195,20 @@ namespace clitk BSplineDeformableTransform ::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 Modified(); } } @@ -267,11 +250,10 @@ namespace clitk // Get the number of parameters template - unsigned int + typename BSplineDeformableTransform::NumberOfParametersType BSplineDeformableTransform ::GetNumberOfParameters(void) const { - // The number of parameters equal OutputDimension * number of // of pixels in the grid region. return ( static_cast( OutputDimension ) * @@ -307,7 +289,7 @@ namespace clitk m_WrappedImage->SetRegions( m_GridRegion ); for (unsigned int j=0; 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 +430,6 @@ namespace clitk } } - // Set the parameters template void @@ -594,12 +575,14 @@ namespace clitk m_CoefficientImage = m_WrappedImage; //JV Wrap jacobian into OutputDimension X Vectorial images - this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); + this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); // Use memset to set the memory - JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_Jacobian.data_block()); + JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_SharedDataBSplineJacobian.data_block()); + memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType)); m_LastJacobianIndex = m_ValidRegion.GetIndex(); + m_NeedResetJacobian = false; for (unsigned int j=0; j - const - typename BSplineDeformableTransform - ::JacobianType & + void BSplineDeformableTransform - ::GetJacobian( const InputPointType & point ) const + ::ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const { - // Can only compute Jacobian if parameters are set via - // SetParameters or SetParametersByValue - // if( m_InputParametersPointer == NULL ) - // { - // 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 //======================================================== // Check if inside mask - if(m_Mask && !(m_Mask->IsInside(point) ) ) + if (m_Mask && !(m_Mask->IsInside(point) ) ) { // Outside: no (deformable) displacement - return this->m_Jacobian; - } + jacobian = m_SharedDataBSplineJacobian; + return; + } - //Get index + //Get index this->TransformPointToContinuousIndex( point, m_Index ); // NOTE: if the support region does not lie totally within the grid // we assume zero displacement and return the input point if ( !this->InsideValidRegion( m_Index ) ) { - return this->m_Jacobian; + jacobian = m_SharedDataBSplineJacobian; + return; } //Compute interpolation weights - const WeightsDataType *weights=NULL; + const WeightsDataType *weights = NULL; m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex); m_SupportRegion.SetIndex( m_LastJacobianIndex ); //Reset the iterators - for ( j = 0; j < OutputDimension; j++ ) + unsigned int j = 0; + for ( j = 0; j < OutputDimension; j++ ) m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion); // For each dimension, copy the weight to the support region @@ -1006,13 +959,15 @@ namespace clitk // go to next coefficient in the support region weights++; } + m_NeedResetJacobian = true; // Return the result - return this->m_Jacobian; + jacobian = m_SharedDataBSplineJacobian; } + template void BSplineDeformableTransform @@ -1037,7 +992,48 @@ namespace clitk } } + template + unsigned + BSplineDeformableTransform + ::SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim) + { + unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels(); + m_JacobianImage[dim]->GetPixelContainer()->SetImportPointer(jacobianDataPointer, numberOfPixels); + return numberOfPixels; + } + + template + void + BSplineDeformableTransform + ::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