# [SciPy-user] Arc-length reparametrization. (newbie question ?)

Lionel Barret De Nazaris lionel@gamr7....
Wed Oct 22 04:18:15 CDT 2008

```Hi,

I answer to myself in the case that someone google this in the future.

Here follows my solution.
This is suboptimal and looks like a naive translation but it does works
and is short (i.e manageable)

# cp are the curve control points, it is an array

# Bernstein polynomial
# (see http://en.wikipedia.org/wiki/Bernstein_blending_function)
t  = poly1d([1, 0])    # t
s  = poly1d([-1, 1])   # 1-t

# Bezier curve as combination of Bernstein polynomial
# ( see http://en.wikipedia.org/wiki/Bezier_curve)
ps        = [s*s*s, 3*t*s*s, 3*s*t*t, t*t*t]  # position
ps_prime  = map(lambda p: p.deriv(), ps)      # first derivative
ps_second = map(lambda p: p.deriv(m = 2), ps) # second derivative

def Y(t, cp):
""" the position on the curve at time t """
ts = array( [p(t) for p in ps] ) # calculating the t coeff
i  = ts * cp                     # appliying the t coeff to the
control points
x = i[0,:].sum()                 # selection the xs and summing them
y = i[1,:].sum()                 # selection the ys and summing them
return x, y

def DY(t, cp):
""" the velocity vector on the curve at time t """
ts = array( [p(t) for p in ps_prime] )
i  = ts * cp
x = i[0,:].sum()
y = i[1,:].sum()
return x, y

def Speed(t, cp):
""" the speed (i.e velocity magnitude) on the curve at time t """
x, y = DY(t, cp)
return sqrt(x*x+y*y)

def ArcLength(t, cp, tmin=0):
""" the curve length  corresponding to time t """

def getCurveParameter(s, cp, max_tries, epsilon):
""" the time t corresponding to the curve length s """
L = ArcLength(1,cp)
t = 0 + s * (1-0)/L
for i in range(max_tries):
F = ArcLength(t, cp) - s
if abs(F) < epsilon:
return t
else:
DF = Speed(t, cp)
t -= F/DF
return t

def getCurveParameter2(s, cp):
""" the time t corresponding to the curve length s """
L = ArcLength(1, cp)
t = 0 + s * (1-0)/L
f = lambda t : ArcLength(t, cp) - s
return newton(f, t)

Lionel Barret De Nazaris wrote:
> Hello all,
>
> I've just inherited from a bunch of not so-good code about cubic
> splines, and as a relative scipy newbie, I was wondering what would be
> the right way to do it.
>
> I need to do Arc-length reparametrization.(
> http://www.math.hmc.edu/~gu/math142/mellon/Differential_Geometry/Geometry_of_curves/Parametric_Curves_and_arc.html
> <http://www.math.hmc.edu/%7Egu/math142/mellon/Differential_Geometry/Geometry_of_curves/Parametric_Curves_and_arc.html>
> )
>
> In the current code, I recognize the use of the newton's method to
> converge on a good estimate but this leave me non-plussed. The heavy use
> of integration to compute the curve length without keeping the
> intermediary results is really weird.
>
> So how would you do it, using the best of what scipy has to offer ?
>
> Note : I've looked as the splines in scipy , but this seems more focused
> on finding the spline that fits the samples, than on this kind of
> manipulation.
> Here we have the spline and its control points.
>
> this the code  (see http://pastebin.com/m7e276652 if indentation is not
> correct) :
>
> <code>
> def derivated_curvilign_abscissa(t, points):
>    return
> Vector(vectorial_derivated_cubic_interpolator(t,points)).magnitude()
>
> def curvilign_abscissa(t, points):
>    """
>     input :
>       points = spline control points [init_point, init_virtual_point,
> end_virtual_point, end_point]
>       t      = a float between 0 and 1
>     ouput :
>       the arc-length between f(t) and the start of the curve (aka
> curvilign_abscissa)
>     """
>    return Vector(points[0]).magnitude() + quad(lambda x
> :derivated_curvilign_abscissa( x, points), 0., t)[0]
>
> def create_curvy_grid(points, samples):
>    """
>    input :
>       points  = spline control points [init_point, init_virtual_point,
> end_virtual_point, end_point]
>       samples = list of of equally spaced t (float) between 0 and 1
>    output :
>        a list of float t  (float) where moving from tn to tn+1 means an
> advancing on the curve of an equal length s.
>    """
>    curve_length = curvilign_abscissa(1., points) -curvilign_abscissa(0.,
> points)
>    # ==
>    def function_to_solve(x, translation) :
>        return - translation + curvilign_abscissa(x, points)
> -curvilign_abscissa(0., points)
>    # ==
>    return [newton(lambda x : function_to_solve(x, curve_length *
> step_grid), 1 - PRECISION) for step_grid in grid]
>    # ==
> </code>
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
>

--
Best Regards,
lionel Barret de Nazaris,
Gamr7 Founder & CTO
=========================
Gamr7 : Cities for Games
http://www.gamr7.com

```