Approximating data with a multi segment cubic bezier curve and a distance as well as a curvature contraint

I found the solution that fulfills my criterea. The solution is to first find a B-Spline that approximates the points in the least square sense and then convert that spline into a multi segment bezier curve. B-Splines do have the advantage that in contrast to bezier curves they will not pass through the control points as well as providing a way to specify a desired "smoothness" of the approximation curve. The needed functionality to generate such a spline is implemented in the FITPACK library to which scipy offers a python binding. Lets suppose I read my data into the lists x and y, then I can do:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
tck,u = interpolate.splprep([x,y],s=3)
unew = np.arange(0,1.01,0.01)
out = interpolate.splev(unew,tck)
plt.figure()
plt.plot(x,y,out[0],out[1])
plt.show()

The result then looks like this:

If I want the curve more smooth, then I can increase the s parameter to splprep. If I want the approximation closer to the data I can decrease the s parameter for less smoothness. By going through multiple s parameters programatically I can find a good parameter that fits the given requirements.

The question though is how to convert that result into a bezier curve. The answer in this email by Zachary Pincus. I will replicate his solution here to give a complete answer to my question:

def b_spline_to_bezier_series(tck, per = False):
  """Convert a parametric b-spline into a sequence of Bezier curves of the same degree.

  Inputs:
    tck : (t,c,k) tuple of b-spline knots, coefficients, and degree returned by splprep.
    per : if tck was created as a periodic spline, per *must* be true, else per *must* be false.

  Output:
    A list of Bezier curves of degree k that is equivalent to the input spline. 
    Each Bezier curve is an array of shape (k+1,d) where d is the dimension of the
    space; thus the curve includes the starting point, the k-1 internal control 
    points, and the endpoint, where each point is of d dimensions.
  """
  from fitpack import insert
  from numpy import asarray, unique, split, sum
  t,c,k = tck
  t = asarray(t)
  try:
    c[0][0]
  except:
    # I can't figure out a simple way to convert nonparametric splines to 
    # parametric splines. Oh well.
    raise TypeError("Only parametric b-splines are supported.")
  new_tck = tck
  if per:
    # ignore the leading and trailing k knots that exist to enforce periodicity 
    knots_to_consider = unique(t[k:-k])
  else:
    # the first and last k+1 knots are identical in the non-periodic case, so
    # no need to consider them when increasing the knot multiplicities below
    knots_to_consider = unique(t[k+1:-k-1])
  # For each unique knot, bring it's multiplicity up to the next multiple of k+1
  # This removes all continuity constraints between each of the original knots, 
  # creating a set of independent Bezier curves.
  desired_multiplicity = k+1
  for x in knots_to_consider:
    current_multiplicity = sum(t == x)
    remainder = current_multiplicity%desired_multiplicity
    if remainder != 0:
      # add enough knots to bring the current multiplicity up to the desired multiplicity
      number_to_insert = desired_multiplicity - remainder
      new_tck = insert(x, new_tck, number_to_insert, per)
  tt,cc,kk = new_tck
  # strip off the last k+1 knots, as they are redundant after knot insertion
  bezier_points = numpy.transpose(cc)[:-desired_multiplicity]
  if per:
    # again, ignore the leading and trailing k knots
    bezier_points = bezier_points[k:-k]
  # group the points into the desired bezier curves
  return split(bezier_points, len(bezier_points) / desired_multiplicity, axis = 0)

So B-Splines, FITPACK, numpy and scipy saved my day :)

1. polygonize data

   find the order of points so you just find the closest points to each other and try them to connect 'by lines'. Avoid to loop back to origin point

2. compute derivation along path

   it is the change of direction of the 'lines' where you hit local min or max there is your control point ... Do this to reduce your input data (leave just control points).

3. curve

   now use these points as control points. I strongly recommend interpolation polynomial for both x and y separately for example something like this:

   x=a0+a1*t+a2*t*t+a3*t*t*t
   y=b0+b1*t+b2*t*t+b3*t*t*t
   

   where a0..a3 are computed like this:

   d1=0.5*(p2.x-p0.x);
   d2=0.5*(p3.x-p1.x);
   a0=p1.x;
   a1=d1;
   a2=(3.0*(p2.x-p1.x))-(2.0*d1)-d2;
   a3=d1+d2+(2.0*(-p2.x+p1.x));
   

   • b0 .. b3 are computed in same way but use y coordinates of course
   • p0..p3 are control points for cubic interpolation curve
   • t =<0.0,1.0> is curve parameter from p1 to p2

   this ensures that position and first derivation is continuous (c1) and also you can use BEZIER but it will not be as good match as this.

[edit1] too sharp edges is a BIG problem

To solve it you can remove points from your dataset before obtaining the control points. I can think of two ways to do it right now ... choose what is better for you

1. remove points from dataset with too high first derivation

   dx/dl or dy/dl where x,y are coordinates and l is curve length (along its path). The exact computation of curvature radius from curve derivation is tricky

2. remove points from dataset that leads to too small curvature radius

   compute intersection of neighboring line segments (black lines) midpoint. Perpendicular axises like on image (red lines) the distance of it and the join point (blue line) is your curvature radius. When the curvature radius is smaller then your limit remove that point ...

   

   now if you really need only BEZIER cubics then you can convert my interpolation cubic to BEZIER cubic like this:

//  ---------------------------------------------------------------------------
//  x=cx[0]+(t*cx[1])+(tt*cx[2])+(ttt*cx[3]); // cubic x=f(t), t = <0,1>
//  ---------------------------------------------------------------------------
//  cubic matrix                           bz4 = it4
//  ---------------------------------------------------------------------------
//  cx[0]=                            (    x0) =                    (    X1)
//  cx[1]=                   (3.0*x1)-(3.0*x0) =           (0.5*X2)         -(0.5*X0)
//  cx[2]=          (3.0*x2)-(6.0*x1)+(3.0*x0) = -(0.5*X3)+(2.0*X2)-(2.5*X1)+(    X0)
//  cx[3]= (    x3)-(3.0*x2)+(3.0*x1)-(    x0) =  (0.5*X3)-(1.5*X2)+(1.5*X1)-(0.5*X0)
//  ---------------------------------------------------------------------------
    const double m=1.0/6.0;
    double x0,y0,x1,y1,x2,y2,x3,y3;
    x0 = X1;           y0 = Y1;
    x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
    x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
    x3 = X2;           y3 = Y2;

In case you need the reverse conversion see:

• Bezier curve with control points within the curve