円柱面上のポイントに最適な軸

Question

円筒形の表面上にあるポイントのうち、最も適合するものを見つけるには、Pythonを使用します。

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)）ランダムである複数の気筒を生成：

これはコードです。

おもしろいことに、svdは、orientationの最初のコンポーネントがゼロ（orientation[0] = 0）の場合、絶対に完璧な軸を返します。

この行にコメントすると、推定軸がオフになります。

アップデート1： でもシリンダー式のleastsqを使用して同じ動作を返します。svdを使用して、興味深いアプローチにもかかわらず

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

を、あなたも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)