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

yellowhat
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 https://stackoverflow.com/questions/22285994/how-to-generate-regular-points-on-cylindrical-surface, 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
        else:
            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(ovec.dot(cylvec))
       
        R = rotMatrixAxisAngle(oaxis, rot)
       
        return points.dot(R)

    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(orientation.dot(asse)) - 1) <= 1e-4, orientation.dot(asse)

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)
        # http://geometry.puzzles.narkive.com/2HaVJ3XF/geometry-equation-of-an-arbitrary-orientated-cylinder
        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


Loading...