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>;
PointSetType::Pointer 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>;
SplineFilterType::Pointer 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 Bspline 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 1D image where each voxel contains the
// (x,y,z) locations of the points
using OutputImageType = itk::Image<unsigned char, 2>;
OutputImageType::Pointer 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);
}
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
writer>SetFileName("spline.png");
writer>SetInput(outputImage);
writer>Update();
return EXIT_SUCCESS;
};
Classes demonstrated¶

template<typename
TInputPointSet
, typenameTOutputImage
>
classBSplineScatteredDataPointSetToImageFilter
: public itk::PointSetToImageFilter<TInputPointSet, TOutputImage> Image filter which provides a Bspline output approximation.
Given an nD image with scattered data, this filter finds a fast approximation to that irregularly spaced data using uniform Bsplines. The traditional method of inverting the observation matrix to find a leastsquares 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 multithreaded. This class generalizes from Lee’s original paper to encompass nD data in mD parametric space and any feasible Bspline 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 Bspline nD 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 Bspline object outputs.
This filter is general in that it accepts nD scattered data in mD parametric dimensions. Input to this filter is an mD 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 Bspline 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 3D displacement field which resides in 4D 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 Bspline parametric domain where each pixel houses the sampled Bspline object value. For a curve fit to 3D points, the output is a 1D 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 2D deformation on 2D points, the output is a 2D 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: “ND C^k BSpline Scattered Data Approximation” by Nicholas J. Tustison, James C. Gee
https://www.insightjournal.org/browse/publication/57 Author
Nicholas J. Tustison
 REFERENCE
S. Lee, G. Wolberg, and S. Y. Shin, “Scattered Data Interpolation
with Multilevel BSplines”, IEEE Transactions on Visualization and Computer Graphics, 3(3):228244, 1997.
 REFERENCE
N.J. Tustison and J.C. Gee, “Generalized nD C^k Scattered Data Approximation
with Confidence Values”, Proceedings of the MIAR conference, August 2006.
 ITK Sphinx Examples: