Get points in the plane of an Euler spiral given by curvature

Question

I am currently working on a visualization of a road network. I am using OpenDrive as my standard for road description. I now have a problem with visualizing the curved parts of roads. These are given as either line, arc or spiral. Line meaning a simple line (curvature = 0), arc is a segment of a circle (curvature is constant).

I now hae a problem with the spiral type. It is given with a start point in the plane (x and y), it's length and a linear curvature function given by 2 parameters: curvatureStart and curvatureEnd. The first is the curvature at the start of the spiral, the second at the end.

From this information I found out, that it is an Euler spiral (also called a clothoid).

Now I am stuck at the point of getting plottable points from these informations about this curve. I want to get a fixed number of points along the spiral. So let's say I want 3 points would mean one at the start of the curve, one at the end and one exactly in the middle along the way.

I hope you can understand my problem. Sadly I haven't found a good solution yet. It also does not have to be totally accurate, an approximation would suffice, as I can only plot points to a certain precision.

Answer 1

There are lots of tricks to doing this that aren't published. Don't implement the stuff in papers (except for Levien's papers), use existing code. IIRC, the auto industry has a "standard" clothoid code library. I can't remember what it's called, though.

Anyway, it's not terribly hard to integrate a two-parameter clothoid. Below is some quick example code I've jotted down. Exactly how you fit the clothoid between two points varies a bit, see the second function for one example.

#k1, k2 are curvature values at endpoints
#t is between 0 and 1
def quadrature(t, k1, k2, steps=32):
  dt = t / dt;
  t2 = 0

  x = 0
  y = 0

  for i in range(steps):
    k = k1*(1.0-t2) + k2*t2;

    #integral of k
    th = (-((t2-2)*k1-k2*t2)*t2)/2

    dx1 = cos(th)
    dy1 = sin(th)

    #this bit exploits the fact that the second derivative of an
    #arc length curve is it's normal, i.e. the perpenticular vector whose magnitude is equal
    #to the curvature.

    dx2 = sin(th)*k
    dy2 = -cos(th)*k

    #taylor
    x += dx1*dt + dx2*0.5*dt*dt
    y += dy1*dt + dy2*0.5*dt*dt

    t2 += dt;

  return (x, y);

#p1, p2 are vectors

def fit_clothoid(t, k1, k2, p1, p2):
   #find start and end points
   start = quadrature(0.0, k1, k2)
   end = quadrature(1.0, k1, k2)
   vec = (end - start)

   #calculate rotation angle to fit clothoid

   th = atan2(p2.y-p1.y, p2.x-p1.x) - atan2(vec.y, vec.x)

   #calculate point at t
   p = quadrature(t, k1, k2) - start

   #scale. . .
   p *= length(p2 - p1) / length(end - start)

   #rotate. . .
   p = rotate(p, th)

   #and finally translate
   return p + p1
```