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
- 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"
// 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;
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;
m_SupportRegion.SetIndex(m_SupportIndex);
for ( i = 0; i < OutputDimension; i++ )
m_ZeroVector[i]=itk::NumericTraits<JacobianValueType>::Zero;
-
}
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();
}
}
// Get the number of parameters
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+#if ITK_VERSION_MAJOR >= 4
+ typename BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::NumberOfParametersType
+#else
unsigned int
+#endif
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 ) *
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]
}
}
-
// Set the parameters
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
void
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;
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
+#if ITK_VERSION_MAJOR >= 4
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();
-
//========================================================
- // 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<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>
+ ::JacobianType &
+ BSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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
//========================================================
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);
// go to next coefficient in the support region
weights++;
}
+ m_NeedResetJacobian = true;
// Return the result
return this->m_Jacobian;
}
+#endif
template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
}
}
+ 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
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff