itk::SparseKernelTransform< TScalarType, NDimensions > Class Template Reference

#include <itkSparseKernelTransform.h>

Inheritance diagram for itk::SparseKernelTransform< TScalarType, NDimensions >:

Collaboration diagram for itk::SparseKernelTransform< TScalarType, NDimensions >:

Public Types
typedef SparseKernelTransform	Self
typedef Transform< TScalarType, NDimensions, NDimensions >	Superclass
typedef SmartPointer< Self >	Pointer
typedef SmartPointer< const Self >	ConstPointer
typedef Superclass::ScalarType	ScalarType
typedef Superclass::ParametersType	ParametersType
typedef Superclass::JacobianType	JacobianType
typedef Superclass::InputPointType	InputPointType
typedef Superclass::OutputPointType	OutputPointType
typedef Superclass::InputVectorType	InputVectorType
typedef Superclass::OutputVectorType	OutputVectorType
typedef DefaultStaticMeshTraits< TScalarType, NDimensions, NDimensions, TScalarType, TScalarType >	PointSetTraitsType
typedef PointSet< InputPointType, NDimensions, PointSetTraitsType >	PointSetType
typedef PointSetType::Pointer	PointSetPointer
typedef PointSetType::PointsContainer	PointsContainer
typedef PointSetType::PointsContainerIterator	PointsIterator
typedef PointSetType::PointsContainerConstIterator	PointsConstIterator
typedef itk::VectorContainer< unsigned long, InputVectorType >	VectorSetType
typedef VectorSetType::Pointer	VectorSetPointer
typedef Eigen::Matrix< TScalarType, NDimensions, NDimensions >	IMatrixType
typedef Eigen::Triplet< TScalarType >	TripletType
typedef Eigen::Matrix< TScalarType, NDimensions, NDimensions >	GMatrixType
typedef Eigen::SparseMatrix< TScalarType >	LMatrixType
typedef Eigen::SparseMatrix< TScalarType >	KMatrixType
typedef Eigen::SparseMatrix< TScalarType >	PMatrixType
typedef Eigen::Matrix< TScalarType, Eigen::Dynamic, Eigen::Dynamic >	YMatrixType
typedef Eigen::Matrix< TScalarType, Eigen::Dynamic, Eigen::Dynamic >	WMatrixType
typedef Eigen::Matrix< TScalarType, Eigen::Dynamic, Eigen::Dynamic >	DMatrixType
typedef Eigen::Matrix< TScalarType, NDimensions, NDimensions >	AMatrixType
typedef Eigen::Matrix< TScalarType, NDimensions, 1 >	BMatrixType
typedef Eigen::Matrix< TScalarType, 1, NDimensions >	RowMatrixType
typedef Eigen::Matrix< TScalarType, NDimensions, 1 >	ColumnMatrixType

Public Member Functions
	itkTypeMacro (SparseKernelTransform, Transform)
	itkNewMacro (Self)
	itkStaticConstMacro (SpaceDimension, unsigned int, NDimensions)
	itkGetObjectMacro (SourceLandmarks, PointSetType)
virtual void	SetSourceLandmarks (PointSetType *)
	itkGetObjectMacro (TargetLandmarks, PointSetType)
virtual void	SetTargetLandmarks (PointSetType *)
	itkGetObjectMacro (Displacements, VectorSetType)
void	ComputeWMatrix (void) const
virtual OutputPointType	TransformPoint (const InputPointType &thisPoint) const
virtual void	SetIdentity ()
virtual void	SetParameters (const ParametersType &)
virtual void	SetFixedParameters (const ParametersType &)
virtual void	UpdateParameters (void) const
virtual const ParametersType &	GetParameters (void) const
virtual const ParametersType &	GetFixedParameters (void) const
virtual void	ComputeJacobianWithRespectToParameters (const InputPointType &in, JacobianType &jacobian) const
virtual void	SetStiffness (double stiffness)
	itkGetMacro (Stiffness, double)

Public Attributes
PointSetPointer	m_SourceLandmarks
PointSetPointer	m_TargetLandmarks

Protected Member Functions
void	PrintSelf (std::ostream &os, Indent indent) const
virtual const GMatrixType &	ComputeG (const InputVectorType &landmarkVector) const
virtual const GMatrixType &	ComputeReflexiveG (PointsIterator) const
virtual void	ComputeDeformationContribution (const InputPointType &inputPoint, OutputPointType &result) const
void	ComputeK () const
void	ComputeL () const
void	ComputeP () const
void	ComputeY () const
void	ComputeD () const
void	ReorganizeW (void) const

Protected Attributes
double	m_Stiffness
VectorSetPointer	m_Displacements
LMatrixType	m_LMatrix
LMatrixType	m_LMatrixInverse
KMatrixType	m_KMatrix
PMatrixType	m_PMatrix
YMatrixType	m_YMatrix
WMatrixType	m_WMatrix
DMatrixType	m_DMatrix
AMatrixType	m_AMatrix
BMatrixType	m_BVector
GMatrixType	m_GMatrix
bool	m_WMatrixComputed
bool	m_LMatrixComputed
bool	m_LInverseComputed
IMatrixType	m_I

Detailed Description

template<class TScalarType, unsigned int NDimensions>
class itk::SparseKernelTransform< TScalarType, NDimensions >

Intended to be a base class for elastic body spline and thin plate spline. This is implemented in as straightforward a manner as possible from the IEEE TMI paper by Davis, Khotanzad, Flamig, and Harms, Vol. 16, No. 3 June 1997. Notation closely follows their paper, so if you have it in front of you, this code will make a lot more sense.

SparseKernelTransform: Provides support for defining source and target landmarks Defines a number of data types used in the computations Defines the mathematical framework used to compute all splines, so that subclasses need only provide a kernel specific to that spline

This formulation allows the stiffness of the spline to be adjusted, allowing the spline to vary from interpolating the landmarks to approximating the landmarks. This part of the formulation is based on the short paper by R. Sprengel, K. Rohr, H. Stiehl. "Thin-Plate Spline Approximation for Image Registration". In 18th International Conference of the IEEE Engineering in Medicine and Biology Society. 1996.

This class was modified to support its use in the ITK registration framework by Rupert Brooks, McGill Centre for Intelligent Machines, Montreal, Canada March 2007. See the Insight Journal Paper by Brooks, R., Arbel, T. "Improvements to the itk::KernelTransform and its subclasses."

Definition at line 76 of file itkSparseKernelTransform.h.

Member Typedef Documentation

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType, NDimensions, NDimensions> itk::SparseKernelTransform< TScalarType, NDimensions >::AMatrixType

'A' matrix typedef. Rotational part of the Affine component

Definition at line 249 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType,NDimensions,1> itk::SparseKernelTransform< TScalarType, NDimensions >::BMatrixType

'B' matrix typedef. Translational part of the Affine component

Definition at line 253 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType,NDimensions,1> itk::SparseKernelTransform< TScalarType, NDimensions >::ColumnMatrixType

Column matrix typedef.

Definition at line 261 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType, Eigen::Dynamic, Eigen::Dynamic> itk::SparseKernelTransform< TScalarType, NDimensions >::DMatrixType

'D' matrix typedef. Deformation component

Definition at line 245 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType, NDimensions, NDimensions> itk::SparseKernelTransform< TScalarType, NDimensions >::GMatrixType

'G' matrix typedef.

Definition at line 220 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType, NDimensions, NDimensions> itk::SparseKernelTransform< TScalarType, NDimensions >::IMatrixType

'I' (identity) matrix typedef.

Definition at line 155 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Superclass::InputPointType itk::SparseKernelTransform< TScalarType, NDimensions >::InputPointType

Standard coordinate point type for this class.

Definition at line 105 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Superclass::InputVectorType itk::SparseKernelTransform< TScalarType, NDimensions >::InputVectorType

Standard vector type for this class.

Definition at line 109 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Superclass::JacobianType itk::SparseKernelTransform< TScalarType, NDimensions >::JacobianType

Jacobian type.

Definition at line 102 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::SparseMatrix<TScalarType> itk::SparseKernelTransform< TScalarType, NDimensions >::KMatrixType

'K' matrix typedef.

Definition at line 228 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::SparseMatrix<TScalarType> itk::SparseKernelTransform< TScalarType, NDimensions >::LMatrixType

'L' matrix typedef.

Definition at line 224 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Superclass::ParametersType itk::SparseKernelTransform< TScalarType, NDimensions >::ParametersType

Parameters type.

Definition at line 99 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::SparseMatrix<TScalarType> itk::SparseKernelTransform< TScalarType, NDimensions >::PMatrixType

'P' matrix typedef.

Definition at line 232 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef DefaultStaticMeshTraits<TScalarType, NDimensions, NDimensions, TScalarType, TScalarType> itk::SparseKernelTransform< TScalarType, NDimensions >::PointSetTraitsType

PointList typedef. This type is used for maintaining lists of points, specifically, the source and target landmark lists.

Definition at line 118 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType,1,NDimensions> itk::SparseKernelTransform< TScalarType, NDimensions >::RowMatrixType

Row matrix typedef.

Definition at line 257 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Superclass::ScalarType itk::SparseKernelTransform< TScalarType, NDimensions >::ScalarType

Scalar type.

Definition at line 96 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef SparseKernelTransform itk::SparseKernelTransform< TScalarType, NDimensions >::Self

Standard class typedefs.

Definition at line 81 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Triplet<TScalarType> itk::SparseKernelTransform< TScalarType, NDimensions >::TripletType

triplets used to fill sparse matrices.

Definition at line 217 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef itk::VectorContainer<unsigned long,InputVectorType> itk::SparseKernelTransform< TScalarType, NDimensions >::VectorSetType

VectorSet typedef.

Definition at line 126 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType, Eigen::Dynamic, Eigen::Dynamic> itk::SparseKernelTransform< TScalarType, NDimensions >::WMatrixType

'W' matrix typedef.

Definition at line 241 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

typedef Eigen::Matrix<TScalarType, Eigen::Dynamic, Eigen::Dynamic> itk::SparseKernelTransform< TScalarType, NDimensions >::YMatrixType

'Y' matrix typedef.

Definition at line 237 of file itkSparseKernelTransform.h.

Member Function Documentation

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeD ( void  ) const

protected

Compute displacements .

Definition at line 199 of file itkSparseKernelTransform.hxx.

{
    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();

    PointsIterator sp  = m_SourceLandmarks->GetPoints()->Begin();
    PointsIterator tp  = m_TargetLandmarks->GetPoints()->Begin();
    PointsIterator end = m_SourceLandmarks->GetPoints()->End();

    m_Displacements->Reserve( numberOfLandmarks );
    typename VectorSetType::Iterator vt = m_Displacements->Begin();

    while( sp != end )
    {
        vt->Value() = tp->Value() - sp->Value();
        vt++;
        sp++;
        tp++;
    }
    //    std::cout<<" Computed displacements "<<m_Displacements<<std::endl;
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeDeformationContribution ( const InputPointType &  thisPoint, OutputPointType &  result  ) const

protectedvirtual

Compute the contribution of the landmarks weighted by the kernel funcion to the global deformation of the space

Default implementation of the the method. This can be overloaded in transforms whose kernel produce diagonal G matrices.

Reimplemented in itk::CompactlySupportedRBFSparseKernelTransform< TScalarType, NDimensions >.

Definition at line 170 of file itkSparseKernelTransform.hxx.

{

    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();

    PointsIterator sp  = m_SourceLandmarks->GetPoints()->Begin();

    for(unsigned int lnd=0; lnd < numberOfLandmarks; lnd++ )
    {
        const GMatrixType & Gmatrix = ComputeG( thisPoint - sp->Value() );
        for(unsigned int dim=0; dim < NDimensions; dim++ )
        {
            for(unsigned int odim=0; odim < NDimensions; odim++ )
            {
                result[ odim ] += Gmatrix(dim, odim ) * m_DMatrix(dim,lnd);
            }
        }
        ++sp;
    }

}

template<class TScalarType , unsigned int NDimensions>

const SparseKernelTransform< TScalarType, NDimensions >::GMatrixType & itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeG ( const InputVectorType &  landmarkVector ) const

protectedvirtual

Compute G(x) This is essentially the kernel of the transform. By overriding this method, we can obtain (among others): Elastic body spline Thin plate spline Volume spline

Reimplemented in itk::CompactlySupportedRBFSparseKernelTransform< TScalarType, NDimensions >.

Definition at line 133 of file itkSparseKernelTransform.hxx.

{
    //
    // Should an Exception be thrown here  ?
    //
    itkWarningMacro(<< "ComputeG() should be reimplemented in the subclass !!");
    return m_GMatrix;
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeK ( void  ) const

protected

Compute K matrix.

Definition at line 417 of file itkSparseKernelTransform.hxx.

{
    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
    GMatrixType G;
    //std::vector<TripletType> tripletList;

    m_KMatrix = KMatrixType( NDimensions * numberOfLandmarks,
                             NDimensions * numberOfLandmarks );
    //m_KMatrix.set_size( NDimensions * numberOfLandmarks,
    //                    NDimensions * numberOfLandmarks );

    //m_KMatrix.fill( 0.0 );

    PointsIterator p1  = m_SourceLandmarks->GetPoints()->Begin();
    PointsIterator end = m_SourceLandmarks->GetPoints()->End();

    // K matrix is symmetric, so only evaluate the upper triangle and
    // store the values in both the upper and lower triangle
    unsigned int i = 0;
    while( p1 != end )
    {
        PointsIterator p2 = p1; // start at the diagonal element
        unsigned int j = i;

        // Compute the block diagonal element, i.e. kernel for pi->pi
        //G = ComputeReflexiveG(p1);

        // force to compute the basis on the diagonal
        const InputVectorType s = p1.Value() - p1.Value();
        G = ComputeG(s); // the basis

        //m_KMatrix.update(G, i*NDimensions, i*NDimensions);
        for(unsigned int d = 0; d < NDimensions; d++)
        {
            if(G(d,d) != 0 )
                m_KMatrix.insert(i*NDimensions+d ,i*NDimensions+d) = G(d,d) + m_Stiffness; // this is as a regularizer
            //tripletList.push_back(TripletType(i*NDimensions+d ,i*NDimensions+d, m_GMatrix(d,d)));
        }

        p2++;
        j++;

        // Compute the upper (and copy into lower) triangular part of K
        // only save the lower part, don't need it
        while( p2 != end )
        {
            const InputVectorType s = p1.Value() - p2.Value();
            G = ComputeG(s); // the basis

            // write value in upper and lower triangle of matrix
            for(unsigned int ii = 0 ; ii < NDimensions; ii++)
                for(unsigned int jj = 0 ; jj < NDimensions; jj++)
                {
                    if (G(ii,jj) != 0)
                    {
                        m_KMatrix.insert(i*NDimensions+ii ,j*NDimensions+jj) = G(ii,jj); // upper
                        m_KMatrix.insert(j*NDimensions+ii ,i*NDimensions+jj) = G(ii,jj); // lower
                    }

                    //tripletList.push_back(TripletType(i*NDimensions+ii ,j*NDimensions+jj, G(ii,jj)));
                    //tripletList.push_back(TripletType(j*NDimensions+ii ,i*NDimensions+jj, G(ii,jj)));
                }

            // m_KMatrix.update(G, i*NDimensions, j*NDimensions);
            // m_KMatrix.update(G, j*NDimensions, i*NDimensions);
            p2++;
            j++;
        }
        p1++;
        i++;
    }
    //std::cout<<"K matrix: "<<std::endl;
    //std::cout<<m_KMatrix<<std::endl;

    //    // shireen: for debugging let's make sure that the L matrix is sparse (based on gaussian basis)
    //    std::ofstream ofs;
    //    ofs.open("Ksparse.csv");
    //    for (int k = 0; k < m_KMatrix.outerSize(); ++k) // column index
    //        for (typename KMatrixType::InnerIterator it(m_KMatrix,k); it; ++it)
    //            ofs << it.row()  << ", " << it.col()  << ", " << it.value() << std::endl;
    //    ofs.close();

    std::cout << "Kmatrix - nnz = "  << m_KMatrix.nonZeros() << std::endl;
    m_KMatrix.makeCompressed();

    //m_KMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeL ( void  ) const

protected

Compute L matrix.

postponed till needing the jacobian for this class

Definition at line 351 of file itkSparseKernelTransform.hxx.

{
    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
    //vnl_matrix<TScalarType> O2(NDimensions*(NDimensions+1),
    //                           NDimensions*(NDimensions+1), 0);

    this->ComputeP();
    this->ComputeK();

    //m_LMatrix.set_size( NDimensions*(numberOfLandmarks+NDimensions+1),
    //                    NDimensions*(numberOfLandmarks+NDimensions+1) );
    //m_LMatrix.fill( 0.0 );
    m_LMatrix = LMatrixType( NDimensions*(numberOfLandmarks+NDimensions+1),
                             NDimensions*(numberOfLandmarks+NDimensions+1) );

    //std::vector<TripletType> tripletList;

    // putting KMATRIX
    //m_LMatrix.update( m_KMatrix, 0, 0 );
    //it.value();
    //it.row(); // row index
    //it.col(); // col index (here it is equal to k)
    //it.index(); // inner index, here it is equal to it.row()
    for (int k = 0; k < m_KMatrix.outerSize(); ++k) // column index
        for (typename KMatrixType::InnerIterator it(m_KMatrix,k); it; ++it)
            m_LMatrix.insert(it.row(), it.col()) = it.value() ;
    //tripletList.push_back( TripletType ( it.row(), it.col(), it.value() ) );

    // putting PMATRIX - will keep the lower only ??
    //m_LMatrix.update( m_PMatrix, 0, m_KMatrix.columns() );
    //m_LMatrix.update( m_PMatrix.transpose(), m_KMatrix.rows(), 0);
    for (int p = 0; p < m_PMatrix.outerSize(); ++p) // column index
        for (typename PMatrixType::InnerIterator it(m_PMatrix,p); it; ++it)
        {
            // fill P -> upper
            m_LMatrix.insert(it.row(), NDimensions*numberOfLandmarks + it.col()) = it.value() ;
            //tripletList.push_back( TripletType ( it.row(), NDimensions*numberOfLandmarks + it.col(), it.value() ) );

            // fill P.transpose -> lower
            m_LMatrix.insert(NDimensions*numberOfLandmarks + it.col(),  it.row()) = it.value() ;
            //tripletList.push_back( TripletType ( NDimensions*numberOfLandmarks + it.col(),  it.row(), it.value() ) );
        }
    //m_LMatrix.update( O2, m_KMatrix.rows(), m_KMatrix.columns());

    //    // shireen: for debugging let's make sure that the L matrix is sparse (based on gaussian basis)
    //    std::ofstream ofs;
    //    ofs.open("Lsparse.csv");
    //    for (int k = 0; k < m_LMatrix.outerSize(); ++k) // column index
    //        for (typename LMatrixType::InnerIterator it(m_LMatrix,k); it; ++it)
    //            ofs << it.row()  << ", " << it.col()  << ", " << it.value() << std::endl;
    //    ofs.close();


    std::cout << "Lmatrix - nnz = "  << m_LMatrix.nonZeros() << std::endl;
    m_LMatrix.makeCompressed();

    //m_LMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
    m_LMatrixComputed=1;
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeP ( ) const

protected

Compute P matrix.

Definition at line 512 of file itkSparseKernelTransform.hxx.

{
    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();

    //IMatrixType I = IMatrixType::Identity();
    //IMatrixType temp;
    InputPointType p;

    //I.set_identity();

    //std::vector<TripletType> tripletList;

    //m_PMatrix.set_size( NDimensions*numberOfLandmarks,
    //                    NDimensions*(NDimensions+1) );
    //m_PMatrix.fill( 0.0 );
    m_PMatrix = PMatrixType(NDimensions*numberOfLandmarks,
                            NDimensions*(NDimensions+1) );
    for (unsigned int i = 0; i < numberOfLandmarks; i++)
    {
        m_SourceLandmarks->GetPoint(i, &p);
        for (unsigned int j = 0; j < NDimensions; j++)
        {
            //temp = I * p[j];
            for(unsigned int d = 0 ; d < NDimensions; d++)
            {
                m_PMatrix.insert(i*NDimensions + d, j*NDimensions + d) = p[j];
                //tripletList.push_back( TripletType( i*NDimensions + d, j*NDimensions + d, p[j] ) );
            }
            //m_PMatrix.update(temp, i*NDimensions, j*NDimensions);
        }

        for(unsigned int d = 0 ; d < NDimensions; d++)
        {
            m_PMatrix.insert(i*NDimensions + d, NDimensions*NDimensions + d) = 1;
            //tripletList.push_back( TripletType( i*NDimensions + d, NDimensions*NDimensions + d, 1 ) );
        }
        //m_PMatrix.update(I, i*NDimensions, NDimensions*NDimensions);
    }

    std::cout << "Pmatrix - nnz = "  << m_PMatrix.nonZeros() << std::endl;
    m_PMatrix.makeCompressed();
    //m_PMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
}

template<class TScalarType , unsigned int NDimensions>

const SparseKernelTransform< TScalarType, NDimensions >::GMatrixType & itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeReflexiveG ( PointsIterator  ) const

protectedvirtual

Compute a G(x) for a point to itself (i.e. for the block diagonal elements of the matrix K. Parameter indicates for which landmark the reflexive G is to be computed. The default implementation for the reflexive contribution is a diagonal matrix where the diagonal elements are the stiffness of the spline.

Definition at line 148 of file itkSparseKernelTransform.hxx.

{
    m_GMatrix = GMatrixType::Zero();
    for(unsigned d = 0; d < NDimensions; d++)
        m_GMatrix(d,d) = m_Stiffness;

    //m_GMatrix.fill( NumericTraits< TScalarType >::Zero );
    //m_GMatrix.fill_diagonal( m_Stiffness );

    return m_GMatrix;
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeWMatrix

(

void

)

const

Compute W matrix.

Definition at line 225 of file itkSparseKernelTransform.hxx.

{
    itk::TimeProbe clock;

    //    std::cout<<"Computing W matrix"<<std::endl;

    //typedef vnl_svd<TScalarType>  SVDSolverType;

    if(!m_LMatrixComputed) {
        this->ComputeL();
    }
    this->ComputeY();

    //SVDSolverType svd( m_LMatrix, 1e-8 );
    //m_WMatrix = svd.solve( m_YMatrix );

    clock.Start();

    //Eigen::BiCGSTAB<LMatrixType>  solver;
    Eigen::BiCGSTAB<LMatrixType, Eigen::IncompleteLUT<double> >  solver;
    solver.preconditioner().setDroptol(1e-10);
    solver.preconditioner().setFillfactor(1000);

    //    Eigen::SparseLU<LMatrixType> solver;
    solver.compute(m_LMatrix);

    if(solver.info()!= Eigen::Success) {
        // decomposition failed
        throw std::runtime_error("LMatrix failed to decompose ...!");
    }

    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
    m_WMatrix = WMatrixType::Zero(NDimensions*(numberOfLandmarks+NDimensions+1), 1);
    m_WMatrix = solver.solve(m_YMatrix);

    //    m_WMatrix = WMatrixType::Random(NDimensions*(numberOfLandmarks+NDimensions+1), 1);
    //    solver.setMaxIterations(1);
    //    int i = 0;
    //    do {
    //        m_WMatrix = solver.solveWithGuess(m_YMatrix,m_WMatrix);
    //        std::cout << "#iteration: " << i << " " << "estimated error: " << solver.error() << std::endl;
    //        ++i;
    //    } while (solver.info()!= Eigen::Success && i<100);

    std::cout << "#iterations: " << solver.iterations() << std::endl;
    std::cout << "estimated error: " << solver.error() << std::endl;

    std::cout  << solver.error() << std::endl;
    if(solver.info() != Eigen::Success) {
        // solving failed
        throw std::runtime_error("solving sparse system failed ...!");
    }

    clock.Stop();
    std::cout << "Computing Wmatrix:" << std::endl;
    std::cout << "Mean: " << clock.GetMean() << std::endl;
    std::cout << "Total: " << clock.GetTotal() << std::endl;

    this->ReorganizeW();
    m_WMatrixComputed=true;
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ComputeY ( void  ) const

protected

Compute Y matrix.

Definition at line 563 of file itkSparseKernelTransform.hxx.

{
    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();

    this->ComputeD();

    typename VectorSetType::ConstIterator displacement =
            m_Displacements->Begin();

    //m_YMatrix.set_size( NDimensions*(numberOfLandmarks+NDimensions+1), 1);
    //m_YMatrix.fill( 0.0 );
    m_YMatrix = YMatrixType::Zero(NDimensions*(numberOfLandmarks+NDimensions+1), 1);
    
    for (unsigned int i = 0; i < numberOfLandmarks; i++)
    {
        for (unsigned int j = 0; j < NDimensions; j++)
        {
            m_YMatrix(i*NDimensions+j, 0) = displacement.Value()[j];
            //m_YMatrix.put(i*NDimensions+j, 0, displacement.Value()[j]);
        }
        displacement++;
    }

    for (unsigned int i = 0; i < NDimensions*(NDimensions+1); i++)
    {
        m_YMatrix(numberOfLandmarks*NDimensions+i, 0) = 0;
        //m_YMatrix.put(numberOfLandmarks*NDimensions+i, 0, 0);
    }
}

template<class TScalarType , unsigned int NDimensions>

const SparseKernelTransform< TScalarType, NDimensions >::ParametersType & itk::SparseKernelTransform< TScalarType, NDimensions >::GetFixedParameters ( void  ) const

virtual

Get Transform Fixed Parameters - Gets the Target Landmarks

Definition at line 893 of file itkSparseKernelTransform.hxx.

{
    this->m_FixedParameters = ParametersType( m_SourceLandmarks->GetNumberOfPoints() * NDimensions );

    PointsIterator itr = m_SourceLandmarks->GetPoints()->Begin();
    PointsIterator end = m_SourceLandmarks->GetPoints()->End();

    unsigned int pcounter = 0;
    while( itr != end )
    {
        InputPointType  landmark = itr.Value();
        for(unsigned int dim=0; dim<NDimensions; dim++)
        {
            this->m_FixedParameters[ pcounter ] = landmark[ dim ];
            pcounter++;
        }
        itr++;
    }

    return this->m_FixedParameters;

}

template<class TScalarType , unsigned int NDimensions>

const SparseKernelTransform< TScalarType, NDimensions >::ParametersType & itk::SparseKernelTransform< TScalarType, NDimensions >::GetParameters ( void  ) const

virtual

Get the Transformation Parameters - Gets the Source Landmarks

Definition at line 879 of file itkSparseKernelTransform.hxx.

{
    this->UpdateParameters();
    return this->m_Parameters;

}

template<class TScalarType, unsigned int NDimensions>

itk::SparseKernelTransform< TScalarType, NDimensions >::itkGetObjectMacro	(	SourceLandmarks	,
		PointSetType
	)

Get the source landmarks list, which we will denote .

template<class TScalarType, unsigned int NDimensions>

itk::SparseKernelTransform< TScalarType, NDimensions >::itkGetObjectMacro	(	TargetLandmarks	,
		PointSetType
	)

Get the target landmarks list, which we will denote .

template<class TScalarType, unsigned int NDimensions>

itk::SparseKernelTransform< TScalarType, NDimensions >::itkGetObjectMacro	(	Displacements	,
		VectorSetType
	)

Get the displacements list, which we will denote , where .

template<class TScalarType, unsigned int NDimensions>

itk::SparseKernelTransform< TScalarType, NDimensions >::itkNewMacro

(

Self

)

New macro for creation of through a Smart Pointer

template<class TScalarType, unsigned int NDimensions>

itk::SparseKernelTransform< TScalarType, NDimensions >::itkStaticConstMacro	(	SpaceDimension	,
		unsigned	int,
		NDimensions
	)

Dimension of the domain space.

template<class TScalarType, unsigned int NDimensions>

itk::SparseKernelTransform< TScalarType, NDimensions >::itkTypeMacro	(	SparseKernelTransform< TScalarType, NDimensions >	,
		Transform
	)

Run-time type information (and related methods).

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::ReorganizeW ( void  ) const

protected

Reorganize the components of W into D (deformable), A (rotation part of affine) and B (translational part of affine ) components.

Warning

This method release the memory of the W Matrix

Definition at line 600 of file itkSparseKernelTransform.hxx.

{
    unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();

    // The deformable (non-affine) part of the registration goes here
    m_DMatrix = DMatrixType::Zero(NDimensions,numberOfLandmarks);
    //m_DMatrix.set_size(NDimensions,numberOfLandmarks);

    unsigned int ci = 0;
    for(unsigned int lnd=0; lnd < numberOfLandmarks; lnd++ )
    {
        for(unsigned int dim=0; dim < NDimensions; dim++ )
        {
            //std::cout << m_WMatrix(ci,0) << std::endl;
            m_DMatrix(dim,lnd) = m_WMatrix(ci++,0);
        }
    }

    // This matrix holds the rotational part of the Affine component
    m_AMatrix = AMatrixType::Zero(NDimensions, NDimensions);
    for(unsigned int j=0; j < NDimensions; j++ )
    {
        for(unsigned int i=0; i < NDimensions; i++ )
        {
            m_AMatrix(i,j) = m_WMatrix(ci++,0);
        }
    }

    // This vector holds the translational part of the Affine component
    m_BVector = BMatrixType::Zero(NDimensions,1);
    for(unsigned int k=0; k < NDimensions; k++ )
    {
        m_BVector(k) = m_WMatrix(ci++,0);
    }

    // release WMatrix memory by assigning a small one.
    m_WMatrix = WMatrixType(1,1);

    m_WMatrixComputed=1;

    //    std::ofstream ofs;
    //    ofs.open("D.csv");
    //    for (int k = 0; k < m_DMatrix.outerSize(); ++k) // column index
    //        for (typename DMatrixType::InnerIterator it(m_DMatrix,k); it; ++it)
    //            ofs << it.row()  << ", " << it.col()  << ", " << it.value() << std::endl;
    //    ofs.close();
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::SetFixedParameters ( const ParametersType &  parameters )

virtual

Set Transform Fixed Parameters: To support the transform file writer this function was added to set the target landmarks similar to the SetParameters function setting the source landmarks

Definition at line 807 of file itkSparseKernelTransform.hxx.

{
    typename PointsContainer::Pointer landmarks = PointsContainer::New();
    const unsigned int numberOfLandmarks =  parameters.Size() / NDimensions;

    landmarks->Reserve( numberOfLandmarks );

    PointsIterator itr = landmarks->Begin();
    PointsIterator end = landmarks->End();

    InputPointType  landMark;

    unsigned int pcounter = 0;
    while( itr != end )
    {
        for(unsigned int dim=0; dim<NDimensions; dim++)
        {
            landMark[ dim ] = parameters[ pcounter ];
            pcounter++;
        }
        itr.Value() = landMark;
        itr++;
    }

    //  m_TargetLandmarks->SetPoints( landmarks );
    m_SourceLandmarks->SetPoints( landmarks );

    // these are invalidated when the source lms change
    m_WMatrixComputed=false;
    m_LMatrixComputed=false;
    m_LInverseComputed=false;

    // you must recompute L and Linv - this does not require the targ lms
    //this->ComputeLInverse();
    this->ComputeL();

}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::SetIdentity ( )

virtual

Compute the Jacobian Matrix of the transformation at one point Set the Transformation Parameters to be an identity transform

Definition at line 745 of file itkSparseKernelTransform.hxx.

{
    this->SetParameters(this->GetFixedParameters());
}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::SetParameters ( const ParametersType &  parameters )

virtual

Set the Transformation Parameters and update the internal transformation. The parameters represent the source landmarks. Each landmark point is represented by NDimensions doubles. All the landmarks are concatenated to form one flat Array<double>.

Definition at line 759 of file itkSparseKernelTransform.hxx.

{
    //    std::cout<<"Setting parameters to "<<parameters<<std::endl;
    typename PointsContainer::Pointer landmarks = PointsContainer::New();
    const unsigned int numberOfLandmarks =  parameters.Size() / NDimensions;
    landmarks->Reserve( numberOfLandmarks );

    PointsIterator itr = landmarks->Begin();
    PointsIterator end = landmarks->End();

    InputPointType  landMark;

    unsigned int pcounter = 0;
    while( itr != end )
    {
        for(unsigned int dim=0; dim<NDimensions; dim++)
        {
            landMark[ dim ] = parameters[ pcounter ];
            pcounter++;
        }
        itr.Value() = landMark;
        itr++;
    }

    // m_SourceLandmarks->SetPoints( landmarks );
    m_TargetLandmarks->SetPoints( landmarks );

    // W MUST be recomputed if the target lms are set
    this->ComputeWMatrix();

    //  if(!m_LInverseComputed) {
    //  this->ComputeLInverse();
    //  }

    // Modified is always called since we just have a pointer to the
    // parameters and cannot know if the parameters have changed.
    this->Modified();

}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::SetSourceLandmarks ( PointSetType *  landmarks )

virtual

Set the source landmarks list.

Definition at line 79 of file itkSparseKernelTransform.hxx.

{
    itkDebugMacro("setting SourceLandmarks to " << landmarks );
    if (this->m_SourceLandmarks != landmarks)
    {
        this->m_SourceLandmarks = landmarks;
        this->UpdateParameters();
        this->Modified();

        // these are invalidated when the source lms change
        m_WMatrixComputed  = false;
        m_LMatrixComputed  = false;
        m_LInverseComputed = false;

        // you must recompute L and Linv - this does not require the targ lms
        // Linverse is only needed ofr Jacobian computation, I will defer this in case Jacobian is needed
        // we will assume by default that this transform is used only for warping, if it is a part of optimization
        // the GetJacobian will feel that the inverse is not computed and will compute it
        //this->ComputeLInverse();
        this->ComputeL();

    }
}

template<class TScalarType, unsigned int NDimensions>

virtual void itk::SparseKernelTransform< TScalarType, NDimensions >::SetStiffness ( double  stiffness )

inlinevirtual

Stiffness of the spline. A stiffness of zero results in the standard interpolating spline. A non-zero stiffness allows the spline to approximate rather than interpolate the landmarks. Stiffness values are usually rather small, typically in the range of 0.001 to 0.1. The approximating spline formulation is based on the short paper by R. Sprengel, K. Rohr, H. Stiehl. "Thin-Plate Spline Approximation for Image Registration". In 18th International Conference of the IEEE Engineering in Medicine and Biology Society. 1996.

Definition at line 199 of file itkSparseKernelTransform.h.

    {m_Stiffness=(stiffness>0)?stiffness:0.0;
        m_LMatrixComputed=false;
        m_LInverseComputed=false;
        m_WMatrixComputed=false;
    }

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::SetTargetLandmarks ( PointSetType *  landmarks )

virtual

Set the target landmarks list.

Definition at line 110 of file itkSparseKernelTransform.hxx.

{
    itkDebugMacro("setting TargetLandmarks to " << landmarks );
    if (this->m_TargetLandmarks != landmarks)
    {
        this->m_TargetLandmarks = landmarks;
        // this is invalidated when the target lms change
        m_WMatrixComputed=false;
        this->ComputeWMatrix();
        this->UpdateParameters();
        this->Modified();
    }

}

template<class TScalarType , unsigned int NDimensions>

SparseKernelTransform< TScalarType, NDimensions >::OutputPointType itk::SparseKernelTransform< TScalarType, NDimensions >::TransformPoint ( const InputPointType &  thisPoint ) const

virtual

Compute L matrix inverse. Compute the position of point in the new space

Definition at line 656 of file itkSparseKernelTransform.hxx.

{

    OutputPointType result;

    typedef typename OutputPointType::ValueType ValueType;

    result.Fill( NumericTraits< ValueType >::Zero );

    this->ComputeDeformationContribution( thisPoint, result );

    // Add the rotational part of the Affine component
    for(unsigned int j=0; j < NDimensions; j++ )
    {
        for(unsigned int i=0; i < NDimensions; i++ )
        {
            result[i] += m_AMatrix(i,j) * thisPoint[j];
        }
    }



    // This vector holds the translational part of the Affine component
    for(unsigned int k=0; k < NDimensions; k++ )
    {
        result[k] += m_BVector(k) + thisPoint[k];
    }

    return result;

}

template<class TScalarType , unsigned int NDimensions>

void itk::SparseKernelTransform< TScalarType, NDimensions >::UpdateParameters ( void  ) const

virtual

Update the Parameters array from the landmarks corrdinates.

Definition at line 851 of file itkSparseKernelTransform.hxx.

{
    this->m_Parameters = ParametersType( m_TargetLandmarks->GetNumberOfPoints() * NDimensions );

    PointsIterator itr = m_TargetLandmarks->GetPoints()->Begin();
    PointsIterator end = m_TargetLandmarks->GetPoints()->End();

    unsigned int pcounter = 0;
    while( itr != end )
    {
        InputPointType  landmark = itr.Value();
        for(unsigned int dim=0; dim<NDimensions; dim++)
        {
            this->m_Parameters[ pcounter ] = landmark[ dim ];
            pcounter++;
        }
        itr++;
    }
}

Member Data Documentation

template<class TScalarType, unsigned int NDimensions>

AMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_AMatrix

mutableprotected

Rotatinoal/Shearing part of the Affine component of the Transformation

Definition at line 348 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

BMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_BVector

mutableprotected

Translational part of the Affine component of the Transformation

Definition at line 351 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

VectorSetPointer itk::SparseKernelTransform< TScalarType, NDimensions >::m_Displacements

protected

The list of displacements. d[i] = q[i] - p[i];

Definition at line 320 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

DMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_DMatrix

mutableprotected

The Deformation matrix. This is an auxiliary matrix that will hold the Deformation (non-affine) part of the transform. Those are the coefficients that will multiply the Kernel function

Definition at line 345 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

GMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_GMatrix

mutableprotected

The G matrix. It is made mutable because m_GMatrix was made an ivar only to avoid copying the matrix at return time

Definition at line 356 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

IMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_I

protected

Identity matrix.

Definition at line 366 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

KMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_KMatrix

mutableprotected

The K matrix.

Definition at line 329 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

bool itk::SparseKernelTransform< TScalarType, NDimensions >::m_LInverseComputed

mutableprotected

Has the L inverse matrix been computed?

Definition at line 363 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

LMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_LMatrix

mutableprotected

The L matrix.

Definition at line 323 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

bool itk::SparseKernelTransform< TScalarType, NDimensions >::m_LMatrixComputed

mutableprotected

Has the L matrix been computed?

Definition at line 361 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

LMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_LMatrixInverse

mutableprotected

The inverse of L, which we also cache.

Definition at line 326 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

PMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_PMatrix

mutableprotected

The P matrix.

Definition at line 332 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

PointSetPointer itk::SparseKernelTransform< TScalarType, NDimensions >::m_SourceLandmarks

The list of source landmarks, denoted 'p'.

Definition at line 265 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

double itk::SparseKernelTransform< TScalarType, NDimensions >::m_Stiffness

protected

Stiffness parameter

Definition at line 316 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

PointSetPointer itk::SparseKernelTransform< TScalarType, NDimensions >::m_TargetLandmarks

The list of target landmarks, denoted 'q'.

Definition at line 268 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

WMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_WMatrix

mutableprotected

The W matrix.

Definition at line 338 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

bool itk::SparseKernelTransform< TScalarType, NDimensions >::m_WMatrixComputed

mutableprotected

Has the W matrix been computed?

Definition at line 359 of file itkSparseKernelTransform.h.

template<class TScalarType, unsigned int NDimensions>

YMatrixType itk::SparseKernelTransform< TScalarType, NDimensions >::m_YMatrix

mutableprotected

The Y matrix.

Definition at line 335 of file itkSparseKernelTransform.h.

---

The documentation for this class was generated from the following files:

• C:/Users/Brig/Documents/ShapeWorksStudio/src/Surfworks/itkSparseKernelTransform.h
• C:/Users/Brig/Documents/ShapeWorksStudio/src/Surfworks/itkSparseKernelTransform.hxx