Download this example

Fit Spline Into Point Set

Note

Wish List Still needs additional work to finish proper creation of example.

Synopsis

Fit a spline to a point set.

Results

Note

Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>

Code

C++

#include "itkBSplineScatteredDataPointSetToImageFilter.h"
#include "itkPointSet.h"
#include "itkImage.h"
#include "itkVectorImage.h"
#include "itkImageFileWriter.h"

int
main()
{
  constexpr unsigned int ParametricDimension = 1;
  constexpr unsigned int DataDimension = 2;

  using DataType = itk::Vector<float, DataDimension>;

  using PointSetType = itk::PointSet<DataType, ParametricDimension>;

  auto pointSet = PointSetType::New();

  PointSetType::PointType param0, param1, param2;

  param0[0] = 0.0;
  DataType p0;
  p0[0] = 10.0;
  p0[1] = 10.0;

  pointSet->SetPoint(0, param0);
  pointSet->SetPointData(0, p0);

  param1[0] = 1.0;
  DataType p1;
  p1[0] = 80.0;
  p1[1] = 50.0;
  pointSet->SetPoint(1, param1);
  pointSet->SetPointData(1, p1);

  param2[0] = 2.0;
  DataType p2;
  p2[0] = 180.0;
  p2[1] = 180.0;
  pointSet->SetPoint(2, param2);
  pointSet->SetPointData(2, p2);

  using ImageType = itk::Image<DataType, ParametricDimension>;
  using SplineFilterType = itk::BSplineScatteredDataPointSetToImageFilter<PointSetType, ImageType>;
  auto splineFilter = SplineFilterType::New();

  int splineorder = 2; // complexity of the spline

  SplineFilterType::ArrayType ncontrol;
  ncontrol[0] = splineorder + 1;
  SplineFilterType::ArrayType closedim;
  closedim[0] = 0;

  ImageType::PointType parametricDomainOrigin;
  parametricDomainOrigin[0] = 0.0;

  ImageType::SpacingType parametricDomainSpacing;
  parametricDomainSpacing[0] = 0.0001; // this determines the sampling of the continuous B-spline object.

  ImageType::SizeType parametricDomainSize;
  parametricDomainSize[0] = 2.0 / parametricDomainSpacing[0] + 1;
  splineFilter->SetGenerateOutputImage(true); // the only reason to turn this off is if one only wants to use the
                                              // control point lattice for further processing
  splineFilter->SetInput(pointSet);
  splineFilter->SetSplineOrder(splineorder);
  splineFilter->SetNumberOfControlPoints(ncontrol);
  splineFilter->SetNumberOfLevels(3);
  splineFilter->SetCloseDimension(closedim);
  splineFilter->SetSize(parametricDomainSize);
  splineFilter->SetSpacing(parametricDomainSpacing);
  splineFilter->SetOrigin(parametricDomainOrigin);
  splineFilter->Update();

  // The output will consist of a 1-D image where each voxel contains the
  // (x,y,z) locations of the points
  using OutputImageType = itk::Image<unsigned char, 2>;
  auto                      outputImage = OutputImageType::New();
  OutputImageType::SizeType size;
  size.Fill(200);

  OutputImageType::IndexType start;
  start.Fill(0);

  OutputImageType::RegionType region(start, size);
  outputImage->SetRegions(region);
  outputImage->Allocate();
  outputImage->FillBuffer(0);

  for (unsigned int i = 0; i < splineFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]; ++i)
  {
    ImageType::IndexType splineIndex;
    splineIndex[0] = i;

    DataType outputPixel = splineFilter->GetOutput()->GetPixel(splineIndex);

    OutputImageType::IndexType index;
    index[0] = outputPixel[0];
    index[1] = outputPixel[1];

    outputImage->SetPixel(index, 255);
  }

  itk::WriteImage(outputImage, "spline.png");

  return EXIT_SUCCESS;
};

Classes demonstrated

template<typename TInputPointSet, typename TOutputImage>
class BSplineScatteredDataPointSetToImageFilter : public itk::PointSetToImageFilter<TInputPointSet, TOutputImage>

Image filter which provides a B-spline output approximation.

Given an n-D image with scattered data, this filter finds a fast approximation to that irregularly spaced data using uniform B-splines. The traditional method of inverting the observation matrix to find a least-squares fit is made obsolete. Therefore, memory issues are not a concern and inverting large matrices is not applicable. In addition, this allows fitting to be multi-threaded. This class generalizes from Lee’s original paper to encompass n-D data in m-D parametric space and any feasible B-spline order as well as the option of specifying a confidence value for each point.

In addition to specifying the input point set, one must specify the number of control points. The specified number of control points must be greater than m_SplineOrder. If one wishes to use the multilevel component of this algorithm, one must also specify the number of levels in the hierarchy. If this is desired, the number of control points becomes the number of control points for the coarsest level. The algorithm then increases the number of control points at each level so that the B-spline n-D grid is refined to twice the previous level.

There are two parts to fitting scattered data: the parameterization assignment problem and the fitting problem given a parameterization. This filter only addresses the second problem in that the user must provide a parametric value for each scattered datum. Different parametric assignment schemes result in different B-spline object outputs.

This filter is general in that it accepts n-D scattered data in m-D parametric dimensions. Input to this filter is an m-D point set with a Vector data type of n dimensions. This means that the parametric values are stored in the points container of the point set whereas the scattered data are stored in the points data container of the point set.

Typical B-spline objects include curves, which have a parametric dimension of 1 and a data dimension of 2 or 3 (depending on the space in which the curve resides) and deformation fields which commonly have parametric and data dimensions of 2 or 3 (again depending on the space of the field). As an example, a curve through a set of 2D points has data dimension 2 and parametric dimension 1. The univariate curve could be represented as: <x(u),y(u)> Another example is a 3D deformation of 3D points, which has parametric dimension 3 and data dimension 3 and can be represented as: <dx(u,v,w), dy(u,v,w), dz(u,v,w)>. However, as mentioned before, the code is general such that, if the user wanted, she could model a time varying 3-D displacement field which resides in 4-D space as <dx(u, v, w, t), dy(u, v, w, t), dz(u, v, w, t)>.

The output is an image defining the sampled B-spline parametric domain where each pixel houses the sampled B-spline object value. For a curve fit to 3-D points, the output is a 1-D image where each voxel contains a vector with the approximated (x,y,z) location. The continuous, finite, rectilinear domain (as well as the sampling rate) is specified via the combination of the SetSpacing() and SetSize() functions. For a 2-D deformation on 2-D points, the output is a 2-D image where each voxel contains the approximated (dx, dy) vector.

The parameterization must be specified using SetPoint, where the actual coordinates of the point are set via SetPointData. For example, to compute a spline through the (ordered) 2D points (5,6) and (7,8), you should use:

using DataType = itk::Vector< float, 2 >;
PointSetType::PointType param0;
param0[0] = 0.0;
DataType p0;
p0[0] =  10.0; p0[1]= 10.0;
pointSet->SetPoint(0, param0);
pointSet->SetPointData( 0, p0 );

PointSetType::PointType param1;
param1[0] = 1.0;
DataType p1;
p1[0] =  80.0; p1[1]= 50.0;
pointSet->SetPoint(1, param1);
pointSet->SetPointData( 1, p1 );

This code was contributed in the Insight Journal paper: “N-D C^k B-Spline Scattered Data Approximation” by Nicholas J. Tustison, James C. Gee

https://www.insight-journal.org/browse/publication/57

Author

Nicholas J. Tustison

REFERENCE

S. Lee, G. Wolberg, and S. Y. Shin, “Scattered Data Interpolation

with Multilevel B-Splines”, IEEE Transactions on Visualization and Computer Graphics, 3(3):228-244, 1997.

REFERENCE

N.J. Tustison and J.C. Gee, “Generalized n-D C^k Scattered Data Approximation

with Confidence Values”, Proceedings of the MIAR conference, August 2006.

ITK Sphinx Examples:

• All ITK Sphinx Examples

See itk::BSplineScatteredDataPointSetToImageFilter for additional documentation.