itk::PSMRBFCorrespondenceInterpolator< VDimension > Class Template Reference

Inheritance diagram for itk::PSMRBFCorrespondenceInterpolator< VDimension >:

Collaboration diagram for itk::PSMRBFCorrespondenceInterpolator< VDimension >:

Public Types
typedef PSMRBFCorrespondenceInterpolator	Self
typedef Point< double, VDimension >	PointType
typedef FunctionBase< PointType, PointType >	Superclass
typedef SmartPointer< Self >	Pointer
typedef SmartPointer< const Self >	ConstPointer

Public Member Functions
	itkStaticConstMacro (Dimension, unsigned int, VDimension)
	itkNewMacro (Self)
	itkTypeMacro (PSMRBFCorrespondenceInterpolator, FunctionBase)
void	SetPointSetA (const std::vector< PointType > &v)
const std::vector< PointType > &	GetPointSetA () const
void	SetPointSetB (const std::vector< PointType > &v)
const std::vector< PointType > &	GetPointSetB () const
void	Initialize ()
	itkGetMacro (Initialized, bool)
virtual PointType	Evaluate (const PointType &pt) const

Protected Member Functions
void	PrintSelf (std::ostream &os, Indent indent) const

Detailed Description

template<unsigned int VDimension>
class itk::PSMRBFCorrespondenceInterpolator< VDimension >

Definition at line 40 of file itkPSMRBFCorrespondenceInterpolator.h.

Member Typedef Documentation

template<unsigned int VDimension>

typedef PSMRBFCorrespondenceInterpolator itk::PSMRBFCorrespondenceInterpolator< VDimension >::Self

Standard class typedefs.

Definition at line 48 of file itkPSMRBFCorrespondenceInterpolator.h.

Member Function Documentation

template<unsigned int VDimension>

PSMRBFCorrespondenceInterpolator< VDimension >::PointType itk::PSMRBFCorrespondenceInterpolator< VDimension >::Evaluate ( const PointType &  pt ) const

virtual

Map a point from the surface domain defined by PointSetA to the domain defined by PointSetB.

Definition at line 157 of file itkPSMRBFCorrespondenceInterpolator.hxx.

{
  if (m_Initialized == false)
    {
      itkExceptionMacro("Function has not been initialized");
    }
  
 // N is the number of points
  const unsigned int N = m_PointSetA.size();

  // X is the input point
  vnl_vector<double> X(VDimension);
  for (unsigned int D = 0; D < VDimension; D++)
    {    X(D) = pt[D];    }

  // Convert PointSetA from points to vnl_vectors-- perhaps store this in the class?
  std::vector<vnl_vector<double> > points(N);
  for (unsigned int i = 0; i < N; i++)
    {
    points[i].set_size(VDimension);
    for (unsigned int D = 0; D < VDimension; D++)
      {
        points[i](D) = m_PointSetA[i][D];
      }

    }

  // Return value in which to accumulate
  PointType ret;
  
  // For each dimension ...
  for (unsigned int D = 0; D < VDimension; D++)
    {
    ret[D] = 0.0f; // initialize to zero
    
    // Compute RBF term in D
    for (unsigned int i = 0; i < N; i++)
      {      
        ret[D] += m_C(i,D) * (X - points[i]).magnitude();      
      }
    
    // Add the polynomial term
    for (unsigned int i = 0; i < VDimension; i++)
      {
        //  std::cout << "ret[D] += m_P(i+1,D) + X(i) ->" << ret[D] << " += " << m_P(i+1,D) << " + " << X(i) << std::endl;
        ret[D] += m_P(i+1,D) * X(i);
      }
    ret[D] += m_P(0,D);

    }

  return ret;
}

template<unsigned int VDimension>

void itk::PSMRBFCorrespondenceInterpolator< VDimension >::Initialize

(

)

Compute the interpolation function.

Definition at line 43 of file itkPSMRBFCorrespondenceInterpolator.hxx.

{ 
  // Make sure we have points on the input
  if (m_PointSetA.size() == 0 || m_PointSetB.size() == 0)
    {
      itkExceptionMacro("Surface point sets have not been specified.");
    }
  
  // Make sure the point set sizes match
  if (m_PointSetA.size() != m_PointSetB.size())
    {
      itkExceptionMacro("Surface point sets A and B are not the same size.");
    }

  // N is the number of points
  const unsigned int N = m_PointSetA.size();

  // TODO: Need algorithm references
  // TODO: Need to implement for at least 2D

  typedef vnl_vector<double> vectype;
  typedef vnl_matrix<double> matrixtype;
  
  // Transforms are from 1 to 2
  vectype b(N+4); // point data from file 2
  vectype x(N+4); // parameters to be solved for
  matrixtype A(N+4,N+4);
  matrixtype Phi(N,N);
  
  // Compute Phi values
  vectype Xi(3);
  vectype Xj(3);
  for (unsigned int i = 0; i < N; i++)
    {
    Xi[0] = (m_PointSetA[i])[0];
    Xi[1] = (m_PointSetA[i])[1];
    Xi[2] = (m_PointSetA[i])[2];
    
    for (unsigned int j = 0; j < N; j++)
      {
      Xj[0] = (m_PointSetA[j])[0];
      Xj[1] = (m_PointSetA[j])[1];
      Xj[2] = (m_PointSetA[j])[2];
      
      // Biharmonic Radial Basis Function
      Phi(i,j) = (Xi - Xj).magnitude();

      // Thin plate spline basis -- probably need some logic to switch on dimensionality
      // Phi(i,j) = dot_product(Xi - Xj,Xi - Xj) * log((Xi - Xj).magnitude() + 1.0e-6);
      }
    }
 
  // Construct the coefficient matrix A from PointSetA and Phi's
  for (unsigned int i = 0; i < 4; i++)
    {
    for (unsigned int j = 0; j < N +4; j++)
      {
      if      (j >= N) A(i,j) = 0.0;
      else if (i == 3) A(i,j) = 1.0;
      else             A(i,j) = (m_PointSetA[j])[i];
      }    
    }

  // Rest of the matrix A
  for (unsigned int i = 0; i < N; i++)
    {
    for (unsigned int j = 0; j < N +4; j++)
      {
      if      (j < N)    A(i+4,j) = Phi(i,j);
      else if (j == N+3) A(i+4,j) = 1.0;
      else               A(i+4,j) = (m_PointSetA[i])[2-(j-N)];
      }
    }

  // To keep things simple, solve in each dimension separately
  //  matrixtype P(4,3);
  //  matrixtype C(N,3);
  m_P.set_size(4,3);
  m_C.set_size(N,3);
  
  for (unsigned int D = 0; D < 3; D++)
    {
    // Construct the position data vector b from file 2 
    // [0 0 0 0 x1' x2' x3' ... xN' ]^T   
    b[0] = b[1] = b[2] = b[3] = 0.0;
    for(unsigned int i = 0; i < N; i++)
      {
      b[i+4] = (m_PointSetB[i])[D];
      }


    // Solve for parameters
    //    x = vnl_matrix_inverse<double>(A) * b;
    x = vnl_svd<double>(A).solve(b);

    // Store results
    for (unsigned int i = 0; i < N; i++)
      {
     m_C(i,D) = x[i];
      }
    
     m_P(0,D) = x[N+3];
     m_P(1,D) = x[N+2];
     m_P(2,D) = x[N+1];
     m_P(3,D) = x[N];

    } // end D

  m_Initialized = true;
}

template<unsigned int VDimension>

itk::PSMRBFCorrespondenceInterpolator< VDimension >::itkGetMacro	(	Initialized	,
		bool
	)

Returns true if the interpolation function has been computed.

template<unsigned int VDimension>

itk::PSMRBFCorrespondenceInterpolator< VDimension >::itkNewMacro

(

Self

)

Method to create through the object factory.

template<unsigned int VDimension>

itk::PSMRBFCorrespondenceInterpolator< VDimension >::itkStaticConstMacro	(	Dimension	,
		unsigned	int,
		VDimension
	)

Dimensionality of the points.

template<unsigned int VDimension>

itk::PSMRBFCorrespondenceInterpolator< VDimension >::itkTypeMacro	(	PSMRBFCorrespondenceInterpolator< VDimension >	,
		FunctionBase
	)

Run-time type information (and related methods).

template<unsigned int VDimension>

void itk::PSMRBFCorrespondenceInterpolator< VDimension >::SetPointSetA ( const std::vector< PointType > &  v )

inline

Set/Get the first point set "A". This is the point set that defines the domain from which points will be mapped by the function.

Definition at line 63 of file itkPSMRBFCorrespondenceInterpolator.h.

    {      
      m_PointSetA = v;
    }

template<unsigned int VDimension>

void itk::PSMRBFCorrespondenceInterpolator< VDimension >::SetPointSetB ( const std::vector< PointType > &  v )

inline

Set/Get the second point set "B". This is the point set surface that defines the domain into which points will be mapped by the function.

Definition at line 75 of file itkPSMRBFCorrespondenceInterpolator.h.

    {      
      m_PointSetB = v;
    }

---

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

• C:/Users/Brig/Documents/ShapeWorksStudio/src/ParticleShapeworks/include/itkPSMRBFCorrespondenceInterpolator.h
• C:/Users/Brig/Documents/ShapeWorksStudio/src/ParticleShapeworks/include/itkPSMRBFCorrespondenceInterpolator.hxx