X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkBSplineDeformableTransform.txx;h=9d7cc6b8e2427d195544854ce540d75390a0789b;hb=29d00baffa433b3f4a60c69116b840c97e470ec9;hp=4290191679d182a02d1f66f1b8bf19c01f912f9c;hpb=657652a78c2e2717a6f77e027049173442ca29f0;p=clitk.git diff --git a/registration/clitkBSplineDeformableTransform.txx b/registration/clitkBSplineDeformableTransform.txx index 4290191..9d7cc6b 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" @@ -32,7 +32,11 @@ namespace clitk // Constructor with default arguments template BSplineDeformableTransform +#if ITK_VERSION_MAJOR >= 4 + ::BSplineDeformableTransform():Superclass(0) +#else ::BSplineDeformableTransform():Superclass(OutputDimension,0) +#endif { unsigned int i; @@ -86,14 +90,7 @@ namespace clitk for (i=0; i ::Zero; - } @@ -204,28 +200,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 +255,14 @@ namespace clitk // Get the number of parameters template +#if ITK_VERSION_MAJOR >= 4 + typename BSplineDeformableTransform::NumberOfParametersType +#else unsigned int +#endif BSplineDeformableTransform ::GetNumberOfParameters(void) const { - // The number of parameters equal OutputDimension * number of // of pixels in the grid region. return ( static_cast( OutputDimension ) * @@ -307,7 +298,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 +439,6 @@ namespace clitk } } - // Set the parameters template void @@ -594,12 +584,21 @@ namespace clitk m_CoefficientImage = m_WrappedImage; //JV Wrap jacobian into OutputDimension X Vectorial images +#if ITK_VERSION_MAJOR >= 4 + this->m_SharedDataBSplineJacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); +#else this->m_Jacobian.set_size( OutputDimension, this->GetNumberOfParameters() ); +#endif // Use memset to set the memory +#if ITK_VERSION_MAJOR >= 4 + JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_SharedDataBSplineJacobian.data_block()); +#else JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_Jacobian.data_block()); +#endif memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType)); m_LastJacobianIndex = m_ValidRegion.GetIndex(); + m_NeedResetJacobian = false; for (unsigned int j=0; j= 4 template - 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" ); - // } + if (m_NeedResetJacobian) + ResetJacobian(); - //======================================================== - // Zero all components of jacobian + // For each dimension, copy the weight to the support region //======================================================== - // 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 + // Check if inside mask + if (m_Mask && !(m_Mask->IsInside(point) ) ) + { + // Outside: no (deformable) displacement + jacobian = m_SharedDataBSplineJacobian; + return; + } + + //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 ) ) + { + jacobian = m_SharedDataBSplineJacobian; + return; + } + + //Compute interpolation weights + const WeightsDataType *weights = NULL; + m_VectorInterpolator->EvaluateWeightsAtContinuousIndex( m_Index, &weights, m_LastJacobianIndex); m_SupportRegion.SetIndex( m_LastJacobianIndex ); - - //Initialize the iterators - for ( j = 0; j < OutputDimension; j++ ) + + //Reset the iterators + unsigned int j = 0; + for ( j = 0; j < OutputDimension; j++ ) m_Iterator[j] = IteratorType( m_JacobianImage[j], m_SupportRegion); - - //Set the previously-set to zero + + // For each dimension, copy the weight to the support region while ( ! (m_Iterator[0]).IsAtEnd() ) { + //copy weight to jacobian image for ( j = 0; j < OutputDimension; j++ ) { - m_Iterator[j].Set( m_ZeroVector ); + m_ZeroVector[j]=*weights; + (m_Iterator[j]).Set( m_ZeroVector); + m_ZeroVector[j]=itk::NumericTraits::Zero; ++(m_Iterator[j]); } + // go to next coefficient in the support region + weights++; } - + m_NeedResetJacobian = true; + + // Return the result + jacobian = m_SharedDataBSplineJacobian; + + } +#else + template + const + typename BSplineDeformableTransform + ::JacobianType & + BSplineDeformableTransform + ::GetJacobian( const InputPointType & point ) const + { + // Can only compute Jacobian if parameters are set via + // SetParameters or SetParametersByValue + // if( m_InputParametersPointer == NULL ) + // { + // itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" ); + // } + + if (m_NeedResetJacobian) + ResetJacobian(); + //======================================================== // For each dimension, copy the weight to the support region //======================================================== @@ -989,6 +1027,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,11 +1045,13 @@ namespace clitk // go to next coefficient in the support region weights++; } + m_NeedResetJacobian = true; // Return the result return this->m_Jacobian; } +#endif template @@ -1037,7 +1078,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