# Best-Fit axis of points on a cylindrical surface Classic List Threaded 1 message Reply | Threaded
Open this post in threaded view
|

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