BSpline with nonuniform knot spacing. More...

#include <mitkItkNonUniformBSpline.h>

Public Types
typedef NonUniformBSpline	Self
typedef Object	Superclass
typedef SmartPointer< Self >	Pointer
typedef SmartPointer< const Self >	ConstPointer
typedef double	ScalarType
typedef itk::Point< ScalarType, TDimension >	PointType
typedef std::vector< PointType >	PointListType
typedef PointListType *	PointListPointer
typedef std::vector< double >	KnotListType
typedef std::vector< double >	CoordinateListType
typedef itk::Point< double, TDimension >	ControlPointType
typedef std::vector < ControlPointType >	ControlPointListType
typedef ControlPointListType *	ControlPointListPointer
typedef std::vector< double >	ChordLengthListType
Public Member Functions
virtual const char *	GetClassName () const
void	SetPoints (PointListType &newPoints)
const PointListType &	GetPoints () const
void	SetKnots (KnotListType &newKnots)
const KnotListType &	GetKnots () const
void	ComputeChordLengths ()
PointType	EvaluateSpline (const Array< double > &p) const
PointType	EvaluateSpline (double t) const
void	ComputeControlPoints ()
void	SetControlPoints (ControlPointListType &ctrlpts)
const ControlPointListType &	GetControlPoints () const
double	NonUniformBSplineFunctionRecursive (unsigned int order, unsigned int i, double t) const
virtual void	SetSplineOrder (unsigned int _arg)
virtual const unsigned int &	GetSplineOrder ()
Static Public Member Functions
static Pointer	New ()
Protected Member Functions
	NonUniformBSpline ()
virtual	~NonUniformBSpline ()
virtual void	PrintSelf (std::ostream &os, Indent indent) const
Protected Attributes
PointListType	m_Points
KnotListType	m_Knots
ControlPointListType	m_ControlPoints
ChordLengthListType	m_ChordLength
ChordLengthListType	m_CumulativeChordLength
unsigned int	m_SplineOrder
unsigned int	m_SpatialDimension

Detailed Description

template<unsigned int TDimension = 3>
class itk::NonUniformBSpline< TDimension >

BSpline with nonuniform knot spacing.

This class is a bspline with nonuniform knot spacing. The use can specify a set of points and a set of knots and the spline will attempt to find the control points which will cause the spline to interpolate the points.

CAUTION: THIS CLASS IS STILL UNDER DEVELOPMENT.

Definition at line 52 of file mitkItkNonUniformBSpline.h.

Member Typedef Documentation

template<unsigned int TDimension = 3>

typedef std::vector<double> itk::NonUniformBSpline< TDimension >::ChordLengthListType

Definition at line 72 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef SmartPointer< const Self > itk::NonUniformBSpline< TDimension >::ConstPointer

Definition at line 62 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef ControlPointListType* itk::NonUniformBSpline< TDimension >::ControlPointListPointer

Definition at line 71 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef std::vector< ControlPointType > itk::NonUniformBSpline< TDimension >::ControlPointListType

Definition at line 70 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef itk::Point<double, TDimension > itk::NonUniformBSpline< TDimension >::ControlPointType

Definition at line 69 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef std::vector<double> itk::NonUniformBSpline< TDimension >::CoordinateListType

Definition at line 68 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef std::vector< double > itk::NonUniformBSpline< TDimension >::KnotListType

Definition at line 67 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef SmartPointer< Self > itk::NonUniformBSpline< TDimension >::Pointer

Definition at line 61 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef PointListType* itk::NonUniformBSpline< TDimension >::PointListPointer

Definition at line 66 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef std::vector< PointType > itk::NonUniformBSpline< TDimension >::PointListType

Definition at line 65 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef itk::Point< ScalarType, TDimension > itk::NonUniformBSpline< TDimension >::PointType

Definition at line 64 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef double itk::NonUniformBSpline< TDimension >::ScalarType

Definition at line 63 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef NonUniformBSpline itk::NonUniformBSpline< TDimension >::Self

Typedefs

Definition at line 59 of file mitkItkNonUniformBSpline.h.

template<unsigned int TDimension = 3>

typedef Object itk::NonUniformBSpline< TDimension >::Superclass

Definition at line 60 of file mitkItkNonUniformBSpline.h.

Constructor & Destructor Documentation

template<unsigned int TDimension>

itk::NonUniformBSpline< TDimension >::NonUniformBSpline ( ) [protected]

Constructor

Definition at line 45 of file mitkItkNonUniformBSpline.txx.

{ 
  // Cubic bspline => 4th order
  m_SplineOrder = 3;
  m_SpatialDimension = TDimension;
}

template<unsigned int TDimension>

itk::NonUniformBSpline< TDimension >::~NonUniformBSpline ( ) [protected, virtual]

Virtual destructor

Destructor

Definition at line 55 of file mitkItkNonUniformBSpline.txx.

{ 
}

Member Function Documentation

template<unsigned int TDimension>

void itk::NonUniformBSpline< TDimension >::ComputeChordLengths ( )

Computes the chord lengths based on the points.

Definition at line 181 of file mitkItkNonUniformBSpline.txx.

{
  m_ChordLength.clear();
  m_CumulativeChordLength.clear();

  m_ChordLength.push_back(0);
  m_CumulativeChordLength.push_back(0);

  double total_chord_length = 0.0;
  ChordLengthListType temp;

  for (::size_t i = 0; i < m_Points.size()-1; i++)
    {
    PointType pt = m_Points[i];
    PointType pt2 = m_Points[i+1];

    double chord = pt.EuclideanDistanceTo(pt2);
    m_ChordLength.push_back(chord);
    total_chord_length = total_chord_length + chord;
    temp.push_back(total_chord_length);
    }

  for (ChordLengthListType::iterator aiter = temp.begin();
       aiter != temp.end();
       aiter++)
    {
    m_CumulativeChordLength.push_back(*aiter/total_chord_length);
    }

  //
  // Debug printouts
  //
#ifdef DEBUG_SPLINE
  std::cout << "Total chord length : " << total_chord_length << std::endl;

  std::cout << "Chord length : " << std::endl;
  for (ChordLengthListType::iterator aiter2 = m_ChordLength.begin();
       aiter2 != m_ChordLength.end();
       aiter2++)
    {
    std::cout << *aiter2 << std::endl;
    }

  std::cout << "Cumulative chord length : " << std::endl;
  for (ChordLengthListType::iterator aiter3 = m_CumulativeChordLength.begin();
       aiter3 != m_CumulativeChordLength.end();
       aiter3++)
    {
    std::cout << *aiter3 << std::endl;
    }
  std::cout << std::endl;
#endif
}

template<unsigned int TDimension>

void itk::NonUniformBSpline< TDimension >::ComputeControlPoints ( )

Compute the control points.

Definition at line 280 of file mitkItkNonUniformBSpline.txx.

References QuadProgPP::t().

{
  unsigned int dim = m_Points[0].GetPointDimension();

#ifdef DEBUG_SPLINE
  std::cout << "Points have dimension : " << dim  << std::endl;
#endif  

  //
  // +1 in cols for radius
  // 
  vnl_matrix<double> data_matrix(m_Points.size(), dim);

  //
  // Form data point matrix
  //
  int rr = 0;
  for (typename PointListType::iterator iter = m_Points.begin();
       iter != m_Points.end();
       iter++)
    {
    PointType pt = (*iter);
    for (unsigned int i = 0; i < dim; i++)
      {
      data_matrix(rr, i) = pt.GetVnlVector()[i];
      }
    rr++;
    }

#ifdef DEBUG_SPLINE
  std::cout << std::endl << "Data matrix" << std::endl;
  std::cout << data_matrix << std::endl;
#endif

  //
  // Form basis function matrix
  //
  //int num_basis_functions = 2 * m_SplineOrder - 1;
  //int num_basis_functions = m_Points.size();
  int num_rows = m_Points.size();

  //
  // Assumes multiplicity k (m_SplineOrder at the ends).
  //
  int num_cols = m_Knots.size() - m_SplineOrder;

  vnl_matrix<double> N_matrix(num_rows, num_cols);

  //N_matrix(0, 0) = 1.0;

  for (int r = 0; r < num_rows; r++)
    {
    for (int c = 0; c < num_cols; c++)
      {
      double t = m_CumulativeChordLength[r];
      N_matrix(r, c) = NonUniformBSplineFunctionRecursive(m_SplineOrder, c, t);
      }
    }

  N_matrix(num_rows-1, num_cols-1) = 1.0;

#ifdef DEBUG_SPLINE
  std::cout << "Basis function matrix : " << std::endl;
  std::cout << N_matrix << std::endl;
#endif

//FIXME: Use the LSQR linear solver here:
  vnl_matrix<double> B;

// = vnl_matrix_inverse<double>(N_matrix.transpose() * N_matrix) * N_matrix.transpose() * data_matrix;

//  vnl_linear_system ls( N_matrix.rows(), N_matrix.cols() );

//  vnl_lsqr solver( ls );


//#ifdef DEBUG_SPLINE
  std::cout << "Control point matrix : " << std::endl;
  std::cout << B << std::endl;
//#endif

  m_ControlPoints.clear();

  for ( unsigned int j = 0; j < B.rows(); j++ )
    {
    vnl_vector<double> v = B.get_row(j);
    itk::Vector<double> iv;
    iv.SetVnlVector(v);
    itk::Point<double, TDimension> pt;
    for ( unsigned int d = 0; d < dim; d++ )
      {
      pt[d] = v(d);
      }
    m_ControlPoints.push_back(pt);
    }
  return;
}

template<unsigned int TDimension>

NonUniformBSpline< TDimension >::PointType itk::NonUniformBSpline< TDimension >::EvaluateSpline ( double t ) const

Definition at line 391 of file mitkItkNonUniformBSpline.txx.

References QuadProgPP::sum().

{
  int i = 0;

  vnl_vector<double> result(TDimension);
  result.fill(0);

  for (typename ControlPointListType::const_iterator cpiter = m_ControlPoints.begin();
       cpiter != m_ControlPoints.end();
       cpiter++)
    {
    ControlPointType pt = *cpiter;
    vnl_vector<double> v = pt.GetVnlVector();

    const double N =  this->NonUniformBSplineFunctionRecursive(m_SplineOrder, i, t);

    for( unsigned j = 0; j < TDimension; j++ )
      {
      result[j] += N * v[j];
      }

    i++;
    }

  double array[TDimension];
  for ( unsigned int d = 0; d < TDimension; d++ )
    {
    array[d] = result[d];
    }

  ControlPointType sum(array);
#ifdef DEBUG_SPLINE
  std::cout << "Result : " << result << std::endl;
  std::cout << "Sum : " << sum << std::endl;
#endif
  
  return sum;
}

template<unsigned int TDimension = 3>

PointType itk::NonUniformBSpline< TDimension >::EvaluateSpline ( const Array< double > & p ) const

Methods for evaluating the spline. The parameterization is always between 0 and 1.

template<unsigned int TDimension = 3>

virtual const char* itk::NonUniformBSpline< TDimension >::GetClassName ( ) const [virtual]

Method for creation through the object factory.

template<unsigned int TDimension>

const NonUniformBSpline< TDimension >::ControlPointListType & itk::NonUniformBSpline< TDimension >::GetControlPoints ( ) const

Get the control points for the spline

Definition at line 254 of file mitkItkNonUniformBSpline.txx.

{
  return this->m_ControlPoints;
}

template<unsigned int TDimension>

const NonUniformBSpline< TDimension >::KnotListType & itk::NonUniformBSpline< TDimension >::GetKnots ( ) const

Get the knot vector.

Definition at line 263 of file mitkItkNonUniformBSpline.txx.

{
  return this->m_Knots;
}

template<unsigned int TDimension>

const NonUniformBSpline< TDimension >::PointListType & itk::NonUniformBSpline< TDimension >::GetPoints ( ) const

Get the points the spline is trying to interpolate.

Definition at line 272 of file mitkItkNonUniformBSpline.txx.

{
  return this->m_Points;
}

template<unsigned int TDimension = 3>

virtual const unsigned int& itk::NonUniformBSpline< TDimension >::GetSplineOrder ( ) [virtual]

template<unsigned int TDimension = 3>

static Pointer itk::NonUniformBSpline< TDimension >::New ( ) [static]

Method for creation through the object factory.

Referenced by mitk::VirtualTrackingTool::VirtualTrackingTool().

template<unsigned int TDimension>

double itk::NonUniformBSpline< TDimension >::NonUniformBSplineFunctionRecursive	(	unsigned int	order,
		unsigned int	i,
		double	t
	)		const

Evaluate the basis functions directly. order - order of the basis function, i.e. 3 = cubic. i - basis function number, zero based. t - parameter of the spline.

Definition at line 141 of file mitkItkNonUniformBSpline.txx.

References QuadProgPP::t().

{
  if (order == 1)
    {
    if (m_Knots[i] <= t && t < m_Knots[i+1])
      {
      return 1;
      }
    else
      {
      return 0;
      }
    }

  //
  // Be careful, we must use the passed in parameter for the order since this
  // function is recursive.
  //
  double numer1 = (t - m_Knots[i]) * NonUniformBSplineFunctionRecursive(order-1, i, t);
  double denom1 = (m_Knots[i+order-1] - m_Knots[i]);
  double val1 = numer1 / denom1;
  if (denom1 == 0 && numer1 == 0)
    val1 = 0;
  else if (denom1 == 0)
    std::cout << "Error : " << denom1 << ", " << numer1 << std::endl;

  double numer2 = (m_Knots[i+order] - t) * NonUniformBSplineFunctionRecursive(order-1, i+1, t);
  double denom2 = (m_Knots[i + order] - m_Knots[i+1]);
  double val2 = numer2 / denom2;
  if (denom2 == 0 && numer2 == 0)
    val2 = 0;
  else if (denom2 == 0)
    std::cout << "Error : " << denom2 << ", " << numer2 << std::endl;

  return val1 + val2;
}

template<unsigned int TDimension>

void itk::NonUniformBSpline< TDimension >::PrintSelf	(	std::ostream &	os,
		Indent	indent
	)		const `[protected, virtual]`

Method to print the object.

Print the object

Definition at line 63 of file mitkItkNonUniformBSpline.txx.

{ 
  Superclass::PrintSelf( os, indent ); 
  os << indent << "NonUniformBSpline(" << this << ")" << std::endl; 

  os << indent << "Chord lengths : " << std::endl;
  for (ChordLengthListType::const_iterator iter = m_CumulativeChordLength.begin();
       iter != m_CumulativeChordLength.end();
       iter++)
    {
    os << indent << indent << *iter << std::endl;
    }
  os << indent << "Knots : " << std::endl;
  for (KnotListType::const_iterator kiter = m_Knots.begin();
       kiter != m_Knots.end();
       kiter++)
    {
    os << indent << indent << *kiter << std::endl;
    }
  os << indent << "Control Points : " << std::endl;
  for (typename ControlPointListType::const_iterator cpiter = m_ControlPoints.begin();
       cpiter != m_ControlPoints.end();
       cpiter++)
    {
    os << indent << indent << *cpiter << std::endl;
    }
}

template<unsigned int TDimension>

void itk::NonUniformBSpline< TDimension >::SetControlPoints ( ControlPointListType & ctrlpts )

Set the control points for the spline.

Definition at line 238 of file mitkItkNonUniformBSpline.txx.

{
  m_ControlPoints.clear();
  for (typename ControlPointListType::iterator iter = ctrlpts.begin();
       iter != ctrlpts.end();
       iter++)
    {
    m_ControlPoints.push_back(*iter);   
    } 
  this->Modified();
}

template<unsigned int TDimension>

void itk::NonUniformBSpline< TDimension >::SetKnots ( KnotListType & knots )

Set the knot vector. Knots may be nonuniformly spaced. Knots will be rescaled to be between 0 and 1.

Set the list of points composing the tube

Definition at line 116 of file mitkItkNonUniformBSpline.txx.

{
  m_Knots.clear();
  
  int len = knots.size();
  double max_knot = knots[len - 1];
         
  typename KnotListType::iterator it;
  typename KnotListType::iterator end;

  it = knots.begin();
  end = knots.end();

  while(it != end)
    {
    m_Knots.push_back(*it/max_knot);
    it++;
    }
      
  this->Modified();
}

template<unsigned int TDimension>

void itk::NonUniformBSpline< TDimension >::SetPoints ( PointListType & points )

Set points which the spline will attempt to interpolate.

Set the list of points composing the tube

Definition at line 96 of file mitkItkNonUniformBSpline.txx.

{
  m_Points.clear();
         
  typename PointListType::iterator it,end;
  it = points.begin();
  end = points.end();
  while(it != end)
    {
    m_Points.push_back(*it);
    it++;
    }
      
  this->Modified();
}

template<unsigned int TDimension = 3>

virtual void itk::NonUniformBSpline< TDimension >::SetSplineOrder ( unsigned int _arg ) [virtual]

Set the order of the spline.