python - python optimize.leastsq: 円を 3D セットの点に当てはめる

Question

3D データ セットに円フィッティング コードを使用しようとしています。必要に応じてz座標を追加するだけで、3Dポイント用に修正しました。私の変更は、あるポイントのセットではうまく機能し、別のポイントではうまく機能しません。コードにエラーがある場合は、コードを確認してください。

import trig_items
import numpy as np
from trig_items import *
from numpy import *
from matplotlib import pyplot as p
from scipy import optimize

# Coordinates of the 3D points
##x = r_[36, 36, 19, 18, 33, 26]
##y = r_[14, 10, 28, 31, 18, 26]
##z = r_[0, 1, 2, 3, 4, 5]

x = r_[ 2144.18908574,  2144.26880854,  2144.05552972,  2143.90303742,  2143.62520676,
  2143.43628579,  2143.14005775,  2142.79919654,  2142.51436023,  2142.11240866,
  2141.68564346,  2141.29333828,  2140.92596405,  2140.3475612,   2139.90848046,
  2139.24661021,  2138.67384709,  2138.03313547,  2137.40301734,  2137.40908256,
  2137.06611224,  2136.50943781,  2136.0553113,   2135.50313189,  2135.07049922,
  2134.62098139,  2134.10459535,  2133.50838433,  2130.6600465,   2130.03537342,
  2130.04047644,  2128.83522468,  2127.79827542,  2126.43513385,  2125.36700593,
  2124.00350543,  2122.68564431,  2121.20709478,  2119.79047011,  2118.38417647,
  2116.90063343,  2115.52685778,  2113.82246629,  2112.21159431,  2110.63180117,
  2109.00713198,  2108.94434529,  2106.82777156,  2100.62343757,  2098.5090226,
  2096.28787738,  2093.91550703,  2091.66075061,  2089.15316429,  2086.69753869,
  2084.3002414,   2081.87590579,  2079.19141866,  2076.5394574,   2073.89128676,
  2071.18786213]
y = r_[ 725.74913818,  724.43874065,  723.15226506,  720.45950581,  717.77827954,
  715.07048092,  712.39633862,  709.73267688,  707.06039438,  704.43405908,
  701.80074596,  699.15371526,  696.5309022,   693.96109921,  691.35585501,
  688.83496327,  686.32148661,  683.80286662,  681.30705568,  681.30530975,
  679.66483676,  678.01922321,  676.32721779,  674.6667554,   672.9658024,
  671.23686095,  669.52021535,  667.84999077,  659.19757984,  657.46179949,
  657.45700508,  654.46901086,  651.38177517,  648.41739432,  645.32356976,
  642.39034578,  639.42628453,  636.51107198,  633.57732055,  630.63825133,
  627.75308356,  624.80162215,  622.01980232,  619.18814892,  616.37688894,
  613.57400131,  613.61535723,  610.4724493,   600.98277781,  597.84782844,
  594.75983001,  591.77946964,  588.74874068,  585.84525834,  582.92311166,
  579.99564481,  577.06666417,  574.30782762,  571.54115037,  568.79760614,
  566.08551098]
z = r_[ 339.77146775,  339.60021095,  339.47645894,  339.47130963,  339.37216218,
  339.4126132,   339.67942046,  339.40917728,  339.39500353,  339.15041461,
  339.38959195,  339.3358209,   339.47764895,  339.17854867,  339.14624071,
  339.16403926,  339.02308811,  339.27011082,  338.97684183,  338.95087698,
  338.97321177,  339.02175448,  339.02543922,  338.88725411,  339.06942374,
  339.0557553,   339.04414618,  338.89234303,  338.95572249,  339.00880416,
  339.00413073,  338.91080374,  338.98214758,  339.01135789,  338.96393537,
  338.73446188,  338.62784913,  338.72443217,  338.74880562,  338.69090173,
  338.50765186,  338.49056867,  338.57353355,  338.6196255,   338.43754399,
  338.27218569,  338.10587265,  338.43880881,  338.28962141,  338.14338705,
  338.25784154,  338.49792568,  338.15572139,  338.52967693,  338.4594245,
  338.1511823,   338.03711207,  338.19144663,  338.22022045,  338.29032321,
  337.8623197 ]

# coordinates of the barycenter
xm = mean(x)
ym = mean(y)
zm = mean(z)

### Basic usage of optimize.leastsq

def calc_R(xc, yc, zc):
    """ calculate the distance of each 3D points from the center (xc, yc, zc) """
    return sqrt((x - xc) ** 2 + (y - yc) ** 2 + (z - zc) ** 2)

def func(c):
    """ calculate the algebraic distance between the 3D points and the mean circle centered at c=(xc, yc, zc) """
    Ri = calc_R(*c)
    return Ri - Ri.mean()

center_estimate = xm, ym, zm
center, ier = optimize.leastsq(func, center_estimate)
##print center

xc, yc, zc = center
Ri       = calc_R(xc, yc, zc)
R        = Ri.mean()
residu   = sum((Ri - R)**2)
print 'R =', R

したがって、x, y, z(コードでコメントされている) の最初のセットでは、うまく機能します。出力はR = 39.0097846735. 2 番目のポイント セット (コメントなし) を使用してコードを実行すると、結果の半径はR = 108576.859834になり、ほぼ直線になります。ラストを描きました。

青い点は指定されたデータ セットで、赤い点は結果の半径の円弧ですR = 108576.859834。与えられたデータセットの半径が結果よりもはるかに小さいことは明らかです。

その他のポイントまとめはこちら。

最小二乗法が正しく機能しないことは明らかです。

この問題の解決を手伝ってください。

アップデート

これが私の解決策です：

### fit 3D arc into a set of 3D points             ###
### output is the centre and the radius of the arc ###
def fitArc3d(arr, eps = 0.0001):
    # Coordinates of the 3D points
    x = numpy.array([arr[k][0] for k in range(len(arr))])
    y = numpy.array([arr[k][4] for k in range(len(arr))])
    z = numpy.array([arr[k][5] for k in range(len(arr))])
    # coordinates of the barycenter
    xm = mean(x)
    ym = mean(y)
    zm = mean(z)
    ### gradient descent minimisation method ###
    pnts = [[x[k], y[k], z[k]] for k in range(len(x))]
    meanP = Point(xm, ym, zm) # mean point
    Ri = [Point(*meanP).distance(Point(*pnts[k])) for k in range(len(pnts))] # radii to the points
    Rm = math.fsum(Ri) / len(Ri) # mean radius
    dR = Rm + 10 # difference between mean radii
    alpha = 0.1
    c = meanP
    cArr = []
    while dR  > eps:
        cArr.append(c)
        Jx = math.fsum([2 * (x[k] - c[0]) * (Ri[k] - Rm) / Ri[k] for k in range(len(Ri))])
        Jy = math.fsum([2 * (y[k] - c[1]) * (Ri[k] - Rm) / Ri[k] for k in range(len(Ri))])
        Jz = math.fsum([2 * (z[k] - c[2]) * (Ri[k] - Rm) / Ri[k] for k in range(len(Ri))])
        gradJ = [Jx, Jy, Jz] # find gradient
        c = [c[k] + alpha * gradJ[k] for k in range(len(c)) if len(c) == len(gradJ)] # find new centre point
        Ri = [Point(*c).distance(Point(*pnts[k])) for k in range(len(pnts))] # calculate new radii
        RmOld = Rm
        Rm = math.fsum(Ri) / len(Ri) # calculate new mean radius
        dR = abs(Rm - RmOld) # new difference between mean radii

    return Point(*c), Rm

これはあまり最適なコードではありません (微調整する時間がありません) が、動作します。

Answer 1

問題はデータと対応するアルゴリズムだと思います。最小二乗法は、単純な勾配法がほぼ方向最小になるように、局所的な放物線最小値を生成する場合にうまく機能します。残念ながら、これは必ずしもデータに当てはまるとは限りません。xcこれは、との概算をいくつかyc固定し、残差の平方和をzcとの関数としてプロットすることで確認できます。R. ブーメラン型の最小値が得られます。開始パラメータによっては、実際の最小値から離れた分岐の 1 つに終わる可能性があります。谷に入ると、これは非常に平坦になり、最大反復数を超えるか、アルゴリズムの許容範囲内で受け入れられる何かが得られます。いつものように、最初のパラメータが良いほど良いと思います。残念ながら、円の弧が小さいため、改善するのは困難です。私は Python の専門家ではありませんが、Pythonleastsqのおかげでヤコビアン法と勾配法を試すことができると思います。許容範囲も試してみてください。要するに、コードは基本的に問題ないように見えますが、データは病的であり、コードをその種のデータに適応させる必要があります。Karimäki からの 2D の非反復ソリューションがあります。適応できるかもしれません。 このメソッドを 3D に変換します。これも見ることができます。きっとあなたはより多くの文献を見つけるでしょう。

シンプレックス アルゴリズムを使用してデータを確認しました。最低なのは、私が言ったように、行儀が悪いことです。ここで、残差関数のカットをいくつか見てください。xy 平面でのみ、合理的な動作が得られます。zr- および xr- 平面の特性により、検索プロセスが非常に困難になります。

そのため、最初にシンプレックス アルゴリズムはいくつかのほぼ安定した解を見つけます。以下のグラフでは、それらを平らなステップとして見ることができます (青の x、紫の y、黄色の z、緑の R)。最後に、アルゴリズムはほぼ平坦だが非常に伸びた谷を下らなければならず、その結果、z と R の最終的な変換が行われます。それでも、許容範囲が不十分な場合、解のように見える多くの領域が予想されます。10^-5 の標準許容誤差では、アルゴリズムは約 350 回の反復後に停止しました。この解、つまり [1899.32, 741.874, 298.696, 248.956] を得るには、10^-10 に設定する必要がありました。

アップデート

前述のように、ソリューションは作業精度と要求される精度によって異なります。したがって、これらの値は組み込みの最小二乗法と比較して異なるため、手作りの勾配法はおそらくよりうまく機能します。それにもかかわらず、これは 2 つのステップを適合させる私のバージョンです。まず、平面をデータに当てはめます。次のステップでは、この平面内に円を当てはめます。どちらのステップも最小二乗法を使用します。今回は、各ステップが重要な形状の最小値を回避するため、機能します。(当然、円弧セグメントが小さくなり、データが実質的に直線上にある場合、平面の適合に問題が発生します。ただし、これはすべてのアルゴリズムで発生します)。

from math import *
from matplotlib import pyplot as plt
from scipy import optimize
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import pprint as pp

dataTupel=zip(xs,ys,zs) #your data from above

# Fitting a plane first
# let the affine plane be defined by two vectors, 
# the zero point P0 and the plane normal n0
# a point p is member of the plane if (p-p0).n0 = 0 

def distanceToPlane(p0,n0,p):
    return np.dot(np.array(n0),np.array(p)-np.array(p0))    

def residualsPlane(parameters,dataPoint):
    px,py,pz,theta,phi = parameters
    nx,ny,nz =sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta)
    distances = [distanceToPlane([px,py,pz],[nx,ny,nz],[x,y,z]) for x,y,z in dataPoint]
    return distances

estimate = [1900, 700, 335,0,0] # px,py,pz and zeta, phi
#you may automize this by using the center of mass data
# note that the normal vector is given in polar coordinates
bestFitValues, ier = optimize.leastsq(residualsPlane, estimate, args=(dataTupel))
xF,yF,zF,tF,pF = bestFitValues

point  = [xF,yF,zF]
normal = [sin(tF)*cos(pF),sin(tF)*sin(pF),cos(tF)]

# Fitting a circle inside the plane
#creating two inplane vectors
sArr=np.cross(np.array([1,0,0]),np.array(normal))#assuming that normal not parallel x!
sArr=sArr/np.linalg.norm(sArr)
rArr=np.cross(sArr,np.array(normal))
rArr=rArr/np.linalg.norm(rArr)#should be normalized already, but anyhow


def residualsCircle(parameters,dataPoint):
    r,s,Ri = parameters
    planePointArr = s*sArr + r*rArr + np.array(point)
    distance = [ np.linalg.norm( planePointArr-np.array([x,y,z])) for x,y,z in dataPoint]
    res = [(Ri-dist) for dist in distance]
    return res

estimateCircle = [0, 0, 335] # px,py,pz and zeta, phi
bestCircleFitValues, ier = optimize.leastsq(residualsCircle, estimateCircle,args=(dataTupel))

rF,sF,RiF = bestCircleFitValues
print bestCircleFitValues

# Synthetic Data
centerPointArr=sF*sArr + rF*rArr + np.array(point)
synthetic=[list(centerPointArr+ RiF*cos(phi)*rArr+RiF*sin(phi)*sArr) for phi in np.linspace(0, 2*pi,50)]
[cxTupel,cyTupel,czTupel]=[ x for x in zip(*synthetic)]

### Plotting
d = -np.dot(np.array(point),np.array(normal))# dot product
# create x,y mesh
xx, yy = np.meshgrid(np.linspace(2000,2200,10), np.linspace(540,740,10))
# calculate corresponding z
# Note: does not work if normal vector is without z-component
z = (-normal[0]*xx - normal[1]*yy - d)/normal[2]

# plot the surface, data, and synthetic circle
fig = plt.figure()
ax = fig.add_subplot(211, projection='3d')
ax.scatter(xs, ys, zs, c='b', marker='o')
ax.plot_wireframe(xx,yy,z)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
bx = fig.add_subplot(212, projection='3d')
bx.scatter(xs, ys, zs, c='b', marker='o')
bx.scatter(cxTupel,cyTupel,czTupel, c='r', marker='o')
bx.set_xlabel('X Label')
bx.set_ylabel('Y Label')
bx.set_zlabel('Z Label')
plt.show()

これは 245 の半径を与えます。これは、他のアプローチが与えたもの (249) に近いです。したがって、誤差範囲内で同じ結果が得られます。

プロットされた結果は妥当に見えます。お役に立てれば。