- 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
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::BSplineDeformableTransform():Superclass(OutputDimension,0)
+ BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::BSplineDeformableTransform():Superclass(0)
{
unsigned int i;
// Get the number of parameters
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- unsigned int
+ typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::GetNumberOfParameters(void) const
{
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<JacobianPixelType *>(this->m_Jacobian.data_block());
+ JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
+
memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
m_LastJacobianIndex = m_ValidRegion.GetIndex();
m_NeedResetJacobian = false;
// JV weights are identical as for transformpoint, could be done simultaneously in metric!!!!
// Compute the Jacobian in one position
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
- const
- typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::JacobianType &
+ void
BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
- ::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();
//========================================================
// 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
unsigned int j = 0;
- for ( j = 0; j < OutputDimension; j++ )
+ 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
m_NeedResetJacobian = true;
// Return the result
- return this->m_Jacobian;
+ jacobian = m_SharedDataBSplineJacobian;
}
+
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
void
BSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>