[SciPy-user] Arc-length reparametrization. (newbie question ?)
Lionel Barret De Nazaris
lionel@gamr7....
Mon Oct 20 04:37:06 CDT 2008
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>
More information about the SciPy-user
mailing list