[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