Best-Fit axis of points on a cylindrical surface

Question

I would like to find the best-fit axis of points that are on a cylindrical surface, using python.

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 How to generate regular points on cylindrical surface, and estimate the axis.

This is the code:

def rotMatrixAxisAngle(axis, theta, theta2deg=False):

    # Load
    from math import radians, cos, sin
    from numpy import array

    # Convert to radians
    if theta2deg:
        theta = radians(theta)

    #
    a = cos(theta/2.0)
    b, c, d = -array(axis)*sin(theta/2.0)

    # Rotation matrix
    R = array([ [a*a+b*b-c*c-d*d,   2.0*(b*c-a*d),   2.0*(b*d+a*c)],
                [2.0*(b*c+a*d),   a*a+c*c-b*b-d*d,   2.0*(c*d-a*b)],
                [2.0*(b*d-a*c),     2.0*(c*d+a*b), a*a+d*d-b*b-c*c] ])

    return R

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

    # 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.

Update 1: 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

We don't have the `angoli` module, so we can't verify that `rotMatrixAxisAngle(oaxis, rot)` works as expected. — Warren Weckesser, Jul 14 '17 at 06:47

score 1 · Answer 1 · answered Jul 02 '17 at 12:55

Despite your interesting approach using svd, you could also do a more intuitive approach with scipy.optimize.leastsq.

This would need a function to calculate distance between the axis and your cloud of points in order to find the best-fitting axis.

The code could be something like shown below (distance_axis_points adapted from alg3dpy):

from numpy.linalg import norm
from numpy.random import rand
from scipy.optimize import leastsq

for i in range(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)

    def dist_axis_points(axis, points):
        axis_pt0 = points.mean(axis=0)
        axis = np.asarray(axis)
        x1 = axis_pt0[0]
        y1 = axis_pt0[1]
        z1 = axis_pt0[2]
        x2 = axis[0]
        y2 = axis[1]
        z2 = axis[2]
        x3 = points[:, 0]
        y3 = points[:, 1]
        z3 = points[:, 2]
        den = ((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)
        t = ((x1**2 + x2 * x3 - x1 * x3 - x1 * x2 +
              y1**2 + y2 * y3 - y1 * y3 - y1 * y2 +
              z1**2 + z2 * z3 - z1 * z3 - z1 * z2)/den)
        projected_pt = t[:, None]*(axis[None, :] - axis_pt0[None, :]) + axis_pt0[None, :]
        return np.sqrt(((points - projected_pt)**2).sum(axis=-1))

    popt, pconv = leastsq(dist_axis_points, x0=[1, 1, 1], args=(points,))
    popt /= norm(popt)

    assert abs(abs(orientation.dot(popt)) - 1) <= 1e-4, orientation.dot(popt)

Thanks for the reply. But I have the exactly same problem. – yellowhat Jul 02 '17 at 13:46 — yellowhat, Jul 02 '17 at 13:46

Best-Fit axis of points on a cylindrical surface

1 Answers1