}
/** Compute the Jacobian Matrix of the transformation at one point */
+#if ITK_VERSION_MAJOR >= 4
+ virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
+ virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
+ }
+#else
virtual const JacobianType& GetJacobian(const InputPointType &point ) const;
+#endif
/** Return the number of parameters that completely define the Transfom */
virtual unsigned int GetNumberOfParameters(void) const;
// VD Add MultipleBSplineDeformableTransform as friend to facilitate wrapping
friend class MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>;
+#if ITK_VERSION_MAJOR >= 4
+ mutable JacobianType m_SharedDataBSplineJacobian;
+#endif
}; //class BSplineDeformableTransform
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+#if ITK_VERSION_MAJOR >= 4
+ ::BSplineDeformableTransform():Superclass(0)
+#else
::BSplineDeformableTransform():Superclass(OutputDimension,0)
+#endif
{
unsigned int i;
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<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
+#else
JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
+#endif
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
+#if ITK_VERSION_MAJOR >= 4
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ void
+ BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::ComputeJacobianWithRespectToParameters(const InputPointType & point, JacobianType & jacobian) const
+ {
+ if (m_NeedResetJacobian)
+ ResetJacobian();
+
+ //========================================================
+ // For each dimension, copy the weight to the support region
+ //========================================================
+
+ // 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 );
+
+ //Reset the iterators
+ 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
+ while ( ! (m_Iterator[0]).IsAtEnd() )
+ {
+ //copy weight to jacobian image
+ for ( j = 0; j < OutputDimension; j++ )
+ {
+ m_ZeroVector[j]=*weights;
+ (m_Iterator[j]).Set( m_ZeroVector);
+ m_ZeroVector[j]=itk::NumericTraits<JacobianValueType>::Zero;
+ ++(m_Iterator[j]);
+ }
+ // go to next coefficient in the support region
+ weights++;
+ }
+ m_NeedResetJacobian = true;
+
+ // Return the result
+ jacobian = m_SharedDataBSplineJacobian;
+
+ }
+#else
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
const
typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
return this->m_Jacobian;
}
+#endif
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint , jacobian);
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
+#endif
const RealType fixedValue = ti.Value();
if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
+#endif
const RealType fixedValue = ti.Value();
return OutputCovariantVectorType();
}
+#if ITK_VERSION_MAJOR >= 4
+ virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( << "DeformationFieldTransform doesn't declare ComputeJacobianWithRespectToParameters" );
+ }
+ virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( << "DeformationFieldTransform doesn't declare ComputeJacobianWithRespectToPosition" );
+ }
+#else
virtual const JacobianType& GetJacobian(const InputPointType &point ) const
{
itkExceptionMacro( << "DeformationFieldTransform doesn't declare GetJacobian" );
return this->m_Jacobian;
}
+#endif
protected:
DeformationFieldTransform();
// Constructor
template<class TScalarType, unsigned int InputDimension, unsigned int OutputDimension, unsigned int SpaceDimension>
DeformationFieldTransform<TScalarType, InputDimension, OutputDimension, SpaceDimension>
+#if ITK_VERSION_MAJOR >= 4
+ ::DeformationFieldTransform():Superclass(1)
+#else
::DeformationFieldTransform():Superclass(OutputDimension,1)
+#endif
{
m_DeformationField=NULL;
m_Interpolator=DefaultInterpolatorType::New();
}
/** Compute the Jacobian Matrix of the transformation at one point */
+#if ITK_VERSION_MAJOR >= 4
+ virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
+ virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
+ }
+#else
virtual const JacobianType& GetJacobian(const InputPointType &point ) const;
+#endif
/** Return the number of parameters that completely define the Transfom */
virtual unsigned int GetNumberOfParameters(void) const;
std::vector<ParametersType> m_parameters;
mutable std::vector<CoefficientImagePointer> m_CoefficientImages;
mutable int m_LastJacobian;
+#if ITK_VERSION_MAJOR >= 4
+ mutable JacobianType m_SharedDataBSplineJacobian;
+#endif
void InitJacobian();
// FIXME it seems not used
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+#if ITK_VERSION_MAJOR >= 4
+ ::MultipleBSplineDeformableTransform() : Superclass(0)
+#else
::MultipleBSplineDeformableTransform() : Superclass(OutputDimension, 0)
+#endif
{
m_nLabels = 1;
m_labels = 0;
return m_trans[lidx]->DeformablyTransformPoint(inputPoint);
}
+#if ITK_VERSION_MAJOR >= 4
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::ComputeJacobianWithRespectToParameters (const InputPointType &point, JacobianType &jacobian) const
+ {
+ if (m_LastJacobian != -1)
+ m_trans[m_LastJacobian]->ResetJacobian();
+
+ int lidx = 0;
+ if (m_labels)
+ lidx = m_labelInterpolator->Evaluate(point) - 1;
+ if (lidx == -1)
+ {
+ m_LastJacobian = lidx;
+ jacobian = this->m_SharedDataBSplineJacobian;
+ return;
+ }
+
+ JacobianType unused;
+ m_trans[lidx]->ComputeJacobianWithRespectToParameters(point, unused);
+ m_LastJacobian = lidx;
+
+ jacobian = this->m_SharedDataBSplineJacobian;
+ }
+
+#else
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::JacobianType &
MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::GetJacobian( const InputPointType & point ) const
{
- /*
- for (unsigned i = 0; i < m_nLabels; ++i)
- m_trans[i]->ResetJacobian();
- */
if (m_LastJacobian != -1)
m_trans[m_LastJacobian]->ResetJacobian();
return this->m_Jacobian;
}
+#endif
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
inline void
MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>::InitJacobian()
{
unsigned numberOfParameters = this->GetNumberOfParameters();
+#if ITK_VERSION_MAJOR >= 4
+ this->m_SharedDataBSplineJacobian.set_size(OutputDimension, numberOfParameters);
+ JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
+#else
this->m_Jacobian.set_size(OutputDimension, numberOfParameters);
JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
+#endif
memset(jacobianDataPointer, 0, numberOfParameters * sizeof (JacobianPixelType));
unsigned tot = 0;
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
const RealType fixedValue = ti.Get();
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
+#endif
// Get the gradient by NearestNeighboorInterpolation:
// which is equivalent to round up the point components.
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
const RealType fixedValue = ti.Get();
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
+#endif
// Get the gradient by NearestNeighboorInterpolation:
// which is equivalent to round up the point components.
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
const RealType fixedValue = ti.Get();
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
+#endif
// Get the gradient by NearestNeighboorInterpolation:
// which is equivalent to round up the point components.
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
const RealType fixedValue = ti.Get();
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
+#endif
// Get the gradient by NearestNeighboorInterpolation:
// which is equivalent to round up the point components.
return OutputCovariantVectorType();
}
+#if ITK_VERSION_MAJOR >= 4
+ virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( << "PointListTransform doesn't declare ComputeJacobianWithRespectToParameters" );
+ }
+ virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( << "PointListTransform doesn't declare ComputeJacobianWithRespectToPosition" );
+ }
+#else
virtual const JacobianType& GetJacobian(const InputPointType &point ) const
{
itkExceptionMacro( << "PointListTransform doesn't declare GetJacobian" );
return this->m_Jacobian;
}
+#endif
protected:
PointListTransform();
// Constructor
template<class TScalarType, unsigned int NDimensions, unsigned int NOutputDimensions>
PointListTransform<TScalarType, NDimensions, NOutputDimensions>
+#if ITK_VERSION_MAJOR >= 4
+ ::PointListTransform():Superclass(1)
+#else
::PointListTransform():Superclass(NOutputDimensions,1)
+#endif
{
m_PointLists.resize(0);
m_PointList.resize(1);
}
/** Compute the Jacobian Matrix of the transformation at one point */
+#if ITK_VERSION_MAJOR >= 4
+ virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
+ virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
+ {
+ itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
+ }
+#else
virtual const JacobianType& GetJacobian(const InputPointType &point ) const;
+#endif
/** Return the number of parameters that completely define the Transfom */
virtual unsigned int GetNumberOfParameters(void) const;
// JV Shape
unsigned int m_TransformShape;
+#if ITK_VERSION_MAJOR >= 4
+ mutable JacobianType m_SharedDataBSplineJacobian;
+#endif
}; //class ShapedBLUTSpatioTemporalDeformableTransform
// Constructor with default arguments
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+#if ITK_VERSION_MAJOR >= 4
+ ::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(0)
+#else
::ShapedBLUTSpatioTemporalDeformableTransform():Superclass(OutputDimension,0)
+#endif
{
unsigned int i;
//=====================================
//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
// JV four rows of three comps of parameters
+#if ITK_VERSION_MAJOR >= 4
+ JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_SharedDataBSplineJacobian.data_block());
+#else
JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(this->m_Jacobian.data_block());
+#endif
memset(jacobianDataPointer, 0, OutputDimension*numberOfPixels*sizeof(JacobianPixelType));
for (unsigned int j=0; j<OutputDimension; j++)
// 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>
+#if ITK_VERSION_MAJOR >= 4
+ void
+ ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::ComputeJacobianWithRespectToParameters( const InputPointType & point, JacobianType & jacobian) const
+#else
const
typename ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::JacobianType &
ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
::GetJacobian( const InputPointType & point ) const
+#endif
{
//========================================================
if(m_Mask && !(m_Mask->IsInside(point) ) )
{
// Outside: no (deformable) displacement
+#if ITK_VERSION_MAJOR >= 4
+ jacobian = m_SharedDataBSplineJacobian;
+ return;
+#else
return this->m_Jacobian;
+#endif
}
// Get index
// we assume zero displacement and return the input point
if ( !this->InsideValidRegion( m_Index ) )
{
+#if ITK_VERSION_MAJOR >= 4
+ jacobian = m_SharedDataBSplineJacobian;
+ return;
+#else
return this->m_Jacobian;
+#endif
}
// Compute interpolation weights
}
// Return the result
+#if ITK_VERSION_MAJOR >= 4
+ jacobian = m_SharedDataBSplineJacobian;
+#else
return this->m_Jacobian;
+#endif
}
// Compute the transform Jacobian.
typedef typename TransformType::JacobianType JacobianType;
- const JacobianType& jacobian =
- this->m_Transform->GetJacobian(
- m_FixedImageSamples[sampleNumber].FixedImagePointValue );
+#if ITK_VERSION_MAJOR >= 4
+ JacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( m_FixedImageSamples[sampleNumber].FixedImagePointValue, jacobian );
+#else
+ const JacobianType & jacobian =
+ this->m_Transform->GetJacobian( m_FixedImageSamples[sampleNumber].FixedImagePointValue );
+#endif
for ( unsigned int mu = 0; mu < m_NumberOfParameters; mu++ ) {
double innerProduct = 0.0;
weights = m_BSplineTransformWeightsArray[sampleNumber];
indices = m_BSplineTransformIndicesArray[sampleNumber];
} else {
+#if ITK_VERSION_MAJOR >= 4
+ m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition(
+ m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
+#else
m_BSplineTransform->GetJacobian(
m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
+#endif
}
for( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
-
+#endif
const RealType fixedValue = ti.Value();
this->m_NumberOfPixelsCounted++;
if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
+#if ITK_VERSION_MAJOR >= 4
+ TransformJacobianType jacobian;
+ this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
+#else
const TransformJacobianType & jacobian =
this->m_Transform->GetJacobian( inputPoint );
-
+#endif
const RealType fixedValue = ti.Value();
this->m_NumberOfPixelsCounted++;