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

Question

I have some geo data (the image below shows the path of a river as red dots) which I want to approximate using a multi segment cubic bezier curve. Through other questions on stackoverflow here and here I found the algorithm by Philip J. Schneider from "Graphics Gems". I successfully implemented it and can report that even with thousands of points it is very fast. Unfortunately that speed comes with some disadvantages, namely that the fitting is done quite sloppily. Consider the following graphic:

The red dots are my original data and the blue line is the multi segment bezier created by the algorithm by Schneider. As you can see, the input to the algorithm was a tolerance which is at least as high as the green line indicates. Nevertheless, the algorithm creates a bezier curve which has too many sharp turns. You see too of these unnecessary sharp turns in the image. It is easy to imagine a bezier curve with less sharp turns for the shown data while still maintaining the maximum tolerance condition (just push the bezier curve a bit into the direction of the magenta arrows). The problem seems to be that the algorithm picks data points from my original data as end points of the individual bezier curves (the magenta arrows point indicate some suspects). With the endpoints of the bezier curves restricted like that, it is clear that the algorithm will sometimes produce rather sharp curvatures.

What I am looking for is an algorithm which approximates my data with a multi segment bezier curve with two constraints:

• the multi segment bezier curve must never be more than a certain distance away from the data points (this is provided by the algorithm by Schneider)
• the multi segment bezier curve must never create curvatures that are too sharp. One way to check for this criteria would be to roll a circle with the minimum curvature radius along the multisegment bezier curve and check whether it touches all parts of the curve along its path. Though it seems there is a better method involving the cross product of the first and second derivative

The solutions I found which create better fits sadly either work only for single bezier curves (and omit the question of how to find good start and end points for each bezier curve in the multi segment bezier curve) or do not allow a minimum curvature contraint. I feel that the minimum curvature contraint is the tricky condition here.

Here another example (this is hand drawn and not 100% precise):

Lets suppose that figure one shows both, the curvature constraint (the circle must fit along the whole curve) as well as the maximum distance of any data point from the curve (which happens to be the radius of the circle in green). A successful approximation of the red path in figure two is shown in blue. That approximation honors the curvature condition (the circle can roll inside the whole curve and touches it everywhere) as well as the distance condition (shown in green). Figure three shows a different approximation to the path. While it honors the distance condition it is clear that the circle does not fit into the curvature any more. Figure four shows a path which is impossible to be approximated with the given constraints because it is too pointy. This example is supposed to illustrate that to properly approximate some pointy turns in the path, it is necessary that the algorithm chooses control points which are not part of the path. Figure three shows that if control points along the path were chosen, the curvature constraint cannot be fulfilled anymore. This example also shows that the algorithm must quit on some inputs as it is not possible to approximate it with the given constraints.

Does there exist a solution to this problem? The solution does not have to be fast. If it takes a day to process 1000 points, then that's fine. The solution does also not have to be optimal in the sense that it must result in a least squares fit.

In the end I will implement this in C and Python but I can read most other languages too.

Carlcarla

This warrants a very important question: do you want to use a Bezier curve, or do you actually want to fit a function to your data? Because you definitely don't want a Bezier curve here if you want to do good science. – Ursula

@Mike'Pomax'Kamermans functions relate each input to exactly one output. As you can see above, an x value can have multiple y values. How would a function help there? – Carlcarla

The same way you're using a Bezier "function". If you want an accurate depiction of the river, using poly-Beziers is kind of weird. You want curves through points, so at the very least we're looking at applying Catmull-Rom curves instead (fun fact: both are Hermite splines and are representations of each other, except Catmull-Rom curves go "through points" with a specific speed: pretty much exactly what you'd need for modeling river data) – Ursula

I dont want curves through points. As I noted in my question, the algorithm by Schneider produces results I do not want by the curve going through some points. I want to approximate the river, not trace it exactly. By using splines I am reusing the points I have as control points. I dont want that. I want a smooth approximation. – Carlcarla

fair enough, but that's where you hit the "you have requirements that are subjective". Your post doesn't actually formalise them, the requirements stay vague: "more than a certain distance", "curvatures that are too sharp" etc are not things we can help you with unless you get more precise. What, to you, is that "certain distance", what does "too sharp" mean, etc. You have a detailed question, but it has the wrong details. – Ursula

"certain distance" is just one of the inputs to the algorithm. It is the maximum distance of each input point from the curve in normal direction. And I explained what "too sharp" means with the example of rolling a sphere with a certain radius (representing the minimum curvature) along it. – Carlcarla

Perhaps you can take advantage of the fact that a Bezier segment always lies within the trapezoid formed by its 4 control points. – Cowslip

Answer 1

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 :)

Carlcarla

Your original post indicate that the desired algorithm needs to be able to take curvature constraints. However, I did not see any such constraints in using FITPACK in your solutions. – Greensboro

@Greensboro the smoothness parameter seems to be enough to control curvature for my use case. The larger s, the less sharp will the curves be. – Carlcarla

OK. Your original post talks about using a circle with minimum curvature radius r to check if the spline has sharp turns. It makes me think that you would like to have quantitative control for spline's curvature (namely, you want the spline's maximum curvature to be smaller than a given value). If qualitative control for curvature is good for you, indeed that parameter "s" you mentioned suits the purpose. – Greensboro

Answer 2

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

Felicitasfelicitate

Hi, thanks for your detailed answer! I dont think it answers my question though. To make my question more clear I added another graphical example. More precisely I dont think that your solution can honor a curvature constraint because it chooses points of the path as control points for the curve. My example shows that for some very pointy paths it is necessary to choose control points other than from the input path. – Carlcarla

@Carlcarla added [edit1]. btw that interpolation cubic is convertibile to Bezier cubic but in interpolation form it is much more manageable ... – Felicitasfelicitate

thanks! this is an incredibly helpful answer, I'll check back after I tried out your suggestions – Carlcarla

I just encountered this post and notice the statement "position and first derivation is continuous (c2)". This is not correct. If the curve is only continuous in first derivative, it is only C1. C2 continuity requires the second derivatives to be continuous. – Greensboro

@Greensboro in english it is coded from C0? I saw both encodings from C0 and from C1 in literature (not in English) so are you certain it should be from C0 (therefore I allways state the variables that are continuous)? in that case I will edit ... – Felicitasfelicitate

@Spektre: I am certain that it is C0 for positional continuity, C1 for first derivative continuity and C2 for 2nd derivative continuity for all the English literatures I have seen. However, I have to admit that I don't know if there are any non-English literatures that adopt a different convention. – Greensboro

Answer 3

The question was posted long ago, but here is a simple solution based on splprep, finding the minimal value of s allowing to fulfill a minimum curvature radius criteria.

route is the set of input points, the first dimension being the number of points.

import numpy as np
from scipy.interpolate import splprep, splev

#The minimum curvature radius we want to enforce
minCurvatureConstraint = 2000
#Relative tolerance on the radius
relTol = 1.e-6

#Initial values for bisection search, should bound the solution
s_0 = 0
minCurvature_0 = 0
s_1 = 100000000 #Should be high enough to produce curvature radius larger than constraint
s_1 *= 2
minCurvature_1 = np.float('inf')

while np.abs(minCurvature_0 - minCurvature_1)>minCurvatureConstraint*relTol:
    s = 0.5 * (s_0 + s_1)
    tck, u  = splprep(np.transpose(route), s=s)
    smoothed_route = splev(u, tck)
    
    #Compute radius of curvature
    derivative1 = splev(u, tck, der=1)
    derivative2 = splev(u, tck, der=2) 
    xprim = derivative1[0]
    xprimprim = derivative2[0]
    yprim = derivative1[1]
    yprimprim = derivative2[1]
    curvature = 1.0 / np.abs((xprim*yprimprim - yprim* xprimprim) / np.power(xprim*xprim + yprim*yprim, 3 / 2))
    minCurvature = np.min(curvature)
    print("s is %g => Minimum curvature radius is %g"%(s,np.min(curvature)))

    #Perform bisection
    if minCurvature > minCurvatureConstraint:
        s_1 = s
        minCurvature_1 = minCurvature
    else:
        s_0 = s
        minCurvature_0 = minCurvature

It may require some refinements such as iterations to find a suitable s_1, but works.

Quarto

Thanks for this answer, in my opinion it is the only one that really answers the question. But there is a small problem to the application of bisection search: although curvature generally decreases for higher s, this is not always the case (based on me trying this approach on random points) – Feudality