Best-Fit axis of points on a cylindrical surface

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
Report Content as Inappropriate

Best-Fit axis of points on a cylindrical surface

This post has NOT been accepted by the mailing list yet.
I would like to find the best-fit axis of points that are on a cylindrical surface.

Seems that scipy.linalg.svd is the function to look for.

So to test out, I decide to generate some points, function makeCylinder, from this thread, and estimate the axis.

This is the code:

    def makeCylinder(radius, length, nlength, alpha, nalpha, center, orientation):
        # Load
        from numpy import array, allclose, linspace, tile, vstack
        from numpy import pi, cos, sin, arccos, cross
        from numpy.linalg import norm
        from angoli import rotMatrixAxisAngle
        # Create the length array
        I = linspace(0, length, nlength)
        # Create alpha array avoid duplication of endpoints
        if int(alpha) == 360:
            A = linspace(0, alpha, num=nalpha, endpoint=False)/180.0*pi
            A = linspace(0, alpha, num=nalpha)/180.0*pi
        # Calculate X and Y
        X = radius * cos(A)
        Y = radius * sin(A)
        # Tile/repeat indices so all unique pairs are present
        pz = tile(I, nalpha)
        px = X.repeat(nlength)
        py = Y.repeat(nlength)
        # Points
        points = vstack(( pz, px, py )).T
        ## Shift to center
        points += array(center) - points.mean(axis=0)
        # Orient tube to new vector
        ovec = orientation / norm(orientation)
        cylvec = array([1,0,0])
        if allclose(cylvec, ovec):
            return points
        # Get orthogonal axis and rotation
        oaxis = cross(ovec, cylvec)
        rot = arccos(
        R = rotMatrixAxisAngle(oaxis, rot)

    from numpy.linalg import norm
    from numpy.random import rand
    from scipy.linalg import svd
    for i in xrange(100):
        orientation = rand(3)
        orientation[0] = 0
        orientation /= norm(orientation)
        # Generate sample points
        points = makeCylinder(radius = 3.0,
                              length = 20.0, nlength = 20,
                              alpha = 360, nalpha = 30,
                              center = [0,0,0],
                              orientation = orientation)
        # Least Square
        uu, dd, vv = svd(points - points.mean(axis=0))
        asse = vv[0]
        assert abs( abs( - 1) <= 1e-4,

As you can see, I generate multiple cylinder whose axis is random (rand(3)).

The funny thing is that svd returns an axis that is absolutely perfect if the first component of orientation is zero (orientation[0] = 0).

If I comment this line the estimated axis is way off.

Even using leastsq on a cylinder equation returns the same behavior:

    def bestLSQ1(points):
        from numpy import array, sqrt
        from scipy.optimize import leastsq
        # Expand
        points = array(points)
        x = points[:,0]
        y = points[:,1]
        z = points[:,2]
        # Calculate the distance of each points from the center (xc, yc, zc)
        def calc_R(xc, yc, zc, u1, u2, u3):
            return sqrt( (x-xc)**2 + (y-yc)**2 + (z-zc)**2 - ( (x-xc)*u1 + (y-yc)*u2 + (z-zc)*u3 )**2 )
        # Calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc, zc)
        def dist(c):
            Ri = calc_R(*c)
            return Ri - Ri.mean()
        # Axes - Minimize residu
        xM, yM, zM = points.mean(axis=0)
        # Calculate the center
        center, ier = leastsq(dist, (xM, yM, zM, 0, 0, 1))
        xc, yc, zc, u1, u2, u3 = center
        asse = u1, u2, u3
        return asse