Learn differential geometry in sympy (it's a CAS)

Differential geometry is hard, especially when you should put it into practice. "Minimum Curvature Variation Curves, Networks and Surfaces for Fair Free-Form Shape Design" is exactly what I'd have wanted to do for a while, but I just didn't grok the contents until now.

The paper presents an algorithm for optimizing a polynomial spline into a fair spline, then extends it to surfaces. I hope to demonstrate this with an image:


On the left you've got a subdivision surface. On the right there's a quarter of cylinder. The curvature across the cylinder doesn't change, so it makes neater reflection than the surface on the left.

Most paraboloids look like shit, at least more shit than spheres or cylinders. For that reason I'd like to have modelling primitives that naturally form them. The funny thing is. I will likely need a CAS to implement it well.

Without CAS, I would have to differentiate the spline curvature integral with more than 12 variables by hand. It'd be several hours of work that would have to be redone if there's an error. Of course, first you'd sort of have to know what needs to be done. So it could take easily weeks to initially implement the algorithm without a CAS of any kind.

As a bonus letting the computer do the rote work helps you concentrate on what is actually going on.


Sympy is a small-featured computer algebra system. It's easy to install in unix if you've got python installed. Type pip install sympy and you have it. You can also use the online shell

I'm using Sympy by treating the equations as immutable python values. I like that way of thinking because it's simple. It's also feasible to think of embedding it into an application.

An Example: curvature integral for parabel

The paper I had trouble learning did lot of things with integrals and curvature functions. So lets calculate a curvature integral for a simple parabel.

I use f(x) = x², but I guess you can use any other polynomial too. First we got to parametrize it and give it a vector representation. In sympy it may look like this:

from sympy import *
t = symbols('t')
function = (t, t**2)
print "C(t) =", function

Now, say you want to compute what this is at point '0', you can do:

[axis.subs(t, 0.0) for axis in function]

You could also just do 0.0**2. Next we got to know how to get a curvature function for parabel. I'm not entirely sure about how to do it yet, but I found something and put it into curvature_2d. Now I can do:

kf = curvature_2d(function)
print "k(t) =", kf

At this point I can calculate curvatures, kf.subs(t, 0.0) gives me the curvature at point 0.0.

Minimum Energy Curve (MEC) would attempt to minimize the integral of a curvature on the spline. The value for current curve would be computed like this:

mec = Integral(kf, (t, -1.0, 1.0)).evalf()

Minimum Variation Curvature Curve would attempt to minimize the variation in curvature. That would be computed like this:

mvcc = Integral(diff(kf, t), (t, -1.0, +1.0)).evalf()
print "mvcc =", mvcc

The evalf in the end means to switch into numerical integration: We want to know the value even though the integral cannot be evaluated.

As a bonus, lets calculate length of the curve:

x_d = diff(function[0], t)
y_d = diff(function[1], t)
length = Integral(sqrt(x_d**2 + y_d**2), (t, -1.0, +1.0)).evalf()
print "length =", length

Finally, it's common to need the equations extracted so you can run a larger computation. Sympy has customizable pretty printing for the equations. It's been also documented how to manipulate the expression tree directly.

Another example: cubic hermite splines

Hermite curves are parametric polynomial curves that are specified by values and first derivatives at the end points of the spline. Here's how to calculate the hermite basis functions from the definition:

from sympy import *
a0, a1, a2, a3, t = symbols('a0 a1 a2 a3 t')
polynomial = a0 + a1*t + a2*t**2 + a3*t**3
polynomial_d = diff(polynomial, t)
p0, v0, p1, v1 = symbols('p0 v0 p1 v1')
spline = [
    Eq(p0, polynomial.subs(t, 0)),
    Eq(v0, polynomial_d.subs(t, 0)),
    Eq(p1, polynomial.subs(t, 1)),
    Eq(v1, polynomial_d.subs(t, 1))]
solution = solve(spline, (a0, a1, a2, a3))
cf = polynomial.subs(solution)
print 'c(t) =', cf

Here we got a system of equations that describe the control points. From there we solve for the elements of a polynomial and substitute them in the original clause. Now we already have the result, but it is desirable to have the basis functions extracted. It happens by taking the coefficients apart and subtracting them from the clause:

for label, var in [('h0', p0), ('h1', v0), ('h2', v1), ('h3', p1)]:
    cf = apart(cf, var)
    coeff = cf.coeff(var)
    cf = simplify(cf - coeff * var)
    print label, '=', coeff
assert cf == 0

To get quintic hermite curves, add two more degrees, the second derivative and the associated control points for accelerations. Well.. That's what was proposed in the paper, but I didn't manage to get the basis functions described for quintic curves.


You may like to see the full code samples in gist. There are some more samples there. I need to learn bit more differential geometry to implement the interesting paper, but I will document the progress with numpy sheets.

I've understood many people believe that technology detracts from learning mathematics because the computer is doing the rote work. I've observed that technology can help you understand things you couldn't understand otherwise.

Similar posts