bspline-regression

This repository provides the code necessary to fit explicit uniform B-splines of any degree to unstructured 2D/3D/ND point data. In the above example, the data points (red) are approximated by a uniform cubic B-spline (blue) with 14 control points (black).

The primary purpose of this repository is to provide a general B-spline solver which is easy to use, has minimal dependencies, but is still performant. The secondary purpose is to provide a starting point for people learning about B-splines, sparse non-linear least squares optimisation, and the damped Newton and Levenberg-Marquardt algorithms in particular.

Author: Richard Stebbing

License: MIT (refer to LICENSE)

Dependencies

This repository is tested to work under Python 2.7 and Python 3.4.

The required dependencies are:

• Numpy
• Scipy
• Sympy
• matplotlib

Getting Started

To generate the above example problem and initialisation:

python generate_example.py 512 1 1e-1 3 14 Example_1.json --seed=0 --frequency=3

The arguments passed to generate_example.py are:

Value	Argument	Description
512	num_data_points	The number of data points to generate.
1	w	The weight applied to each squared residual.
1e-1	lambda_	The regularisation weight.
3	degree	The degree of the uniform B-spline.
14	num_control_points	The number of control points.
Example_1.json	output_path	The output path.
--seed	0	The random number generator seed (optional).
--frequency	3	The frequency of the sin (optional).

Note:

• w can be a pair (2D) or triple (3D) so that w[i] is the weight applied to the ith dimension of each squared residual. This is useful for simulating time-series data where the uncertainty of the measurement (y-axis) is much greater than that of the reported time (x-axis).
• Increasing lambda_ increases the "force" pulling adjacent B-spline control points together.

The output Example_1.json is a dictionary of the form:

Key	Value Type	Description
degree	int	The degree of the uniform B-spline.
num_control_points	int	The number of control points.
dim	int	The dimension of the problem.
is_closed	bool	True if the B-spline is closed, False otherwise.
Y	array_like, shape = (num_data_points, dim)	The array of data points.
w	array_like, shape = (num_data_points, dim)	The weight of each squared residual on each dimension.
lambda_	float	The regularisation weight.
X	array_like, shape = (num_control_points, dim)	The array of B-spline control points.
u	array_like, shape = (num_data_points, 1)	The array of correspondences (preimages).
where u[i] gives the coordinate of the point on the B-spline that is closest to Y[i].

In summary, the first four keys define the structure of the B-spline, the next three keys specify the data points and problem configuration, and the final two keys define the state. In generate_example.py, X is initialised to a straight line and u is set approximately.

Visualising the initialisation is straightforward:

python visualise.py Example_1.json --empty

where `--empty` generates the plot *without* axis labels or a title, and the correspondences are shown in orange. (All figures in this document were generated with the additional arguments `--width=8` and `--height=4` but this is omitted for brevity.)

To solve for X and u:

python fit_uniform_bspline.py Example_1.json Example_1_Output.json

and visualise the output:

python visualise.py Example_1_Output.json --empty

Alternatively, to save and visualise all optimisation steps:

python fit_uniform_bspline.py Example_1.json Example_1_Output --output-all
python visualise.py Example_1_Output Example_1_Output_Visualisation

Additional arguments to visualise.py and fit_uniform_bspline.py can be found with --help. Further details regarding the solver implementation are provided in the docstring for UniformBSplineLeastSquaresOptimiser.minimise in fit_uniform_bspline.py.

Additional Examples

• Fitting a uniform quadratic B-spline with 5 control points and penalising errors in the x- direction more heavily:

python generate_example.py 256 "1e2,1" 1e-1 2 5 Example_2.json --seed=0 --frequency=0.75 --sigma=0.1 --alpha=-0.1
python fit_uniform_bspline.py Example_2.json Example_2_Output.json
python visualise.py Example_2_Output.json --empty

Solving with Levenberg-Marquardt instead of damped Newton:

python fit_uniform_bspline.py Example_2.json Example_2_Output_LM.json lm
python visualise.py Example_2_Output_LM.json --empty

(For comparison, damped Newton converges at 4.18 whereas Levenberg-Marquardt converges at 6.96. The initial energy is 1.97e3.)

• Fitting a uniform quintic B-spline with 9 control points to 65536 data points:

python generate_example.py 65536 1 1e-1 5 9 Example_3.json --seed=0 --frequency=2 --alpha=0
python fit_uniform_bspline.py Example_3.json Example_3_Output.json
python visualise.py Example_3_Output.json -d u -d X --empty

• Fitting a uniform quartic B-spline with 10 control points in 3D:

python generate_example.py 128 1 1e-1 4 10 Example_4.json --seed=0 --frequency=3 --dim=3
python fit_uniform_bspline.py Example_4.json Example_4_Output.json
python visualise.py Example_4_Output.json --empty