[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