圆柱面上点的最佳拟合轴

Question

我想使用 python 找到最适合的axis 圆柱表面上的点。

似乎scipy.linalg.svd 是要查找的函数。

所以为了测试，我决定从这个线程How to generate regular points on cylindrical surface生成一些点，函数makeCylinder，并估计轴。

这是代码：

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)

如您所见，我生成了多个轴是随机的圆柱体 (rand(3))。

有趣的是，如果orientation 的第一个分量为零 (orientation[0] = 0)，svd 返回一个绝对完美的轴。

如果我评论这条线，估计的轴就会偏离。

更新 1： 即使在圆柱方程上使用 leastsq 也会返回相同的行为：

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

Answer 1

尽管您使用svd 的方法很有趣，但您也可以使用scipy.optimize.leastsq 进行更直观的方法。

这需要一个函数来计算轴和点云之间的距离，以便找到最合适的轴。

代码可能如下所示（distance_axis_points 改编自 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)