[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 """
return quad(Speed, tmin, t, args=(cp,))[0]
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
More information about the SciPy-user
mailing list