7

私は4つのポイントを持っていますが、これは1つの平面に非常に近いです-それは1,4-ジヒドロピリジンサイクルです。

C3とN1からC1-C2-C4-C5でできている平面までの距離を計算する必要があります。距離の計算はOKですが、平面のフィッティングは私にはかなり難しいです。

1,4-DHPサイクル：

1,4-DHPサイクル

1,4-DHPサイクル、別のビュー：

1,4-DHPサイクル、別のビュー

from array import *
from numpy import *
from scipy import *

# coordinates (XYZ) of C1, C2, C4 and C5
x = [0.274791784, -1.001679346, -1.851320839, 0.365840754]
y = [-1.155674199, -1.215133985, 0.053119249, 1.162878076]
z = [1.216239624, 0.764265677, 0.956099579, 1.198231236]

# plane equation Ax + By + Cz = D
# non-fitted plane
abcd = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

# creating distance variable
distance =  zeros(4, float)

# calculating distance from point to plane
for i in range(4):
    distance[i] = (x[i]*abcd[0]+y[i]*abcd[1]+z[i]*abcd[2]+abcd[3])/sqrt(abcd[0]**2 + abcd[1]**2 + abcd[2]**2)
    
print distance

# calculating squares
squares = distance**2

print squares

合計（平方）を最小化する方法は？最小二乗法を試しましたが、それも私にはありません。

4

6 に答える 6

19

それはほぼ正しいように聞こえますが、非線形最適化をSVDに置き換える必要があります。以下は、慣性モーメントテンソルMを作成し、次にSVDで平面の法線を取得します。これは最小二乗近似に近似し、はるかに高速で予測可能である必要があります。点群中心と法線を返します。

def planeFit(points):
    """
    p, n = planeFit(points)

    Given an array, points, of shape (d,...)
    representing points in d-dimensional space,
    fit an d-dimensional plane to the points.
    Return a point, p, on the plane (the point-cloud centroid),
    and the normal, n.
    """
    import numpy as np
    from numpy.linalg import svd
    points = np.reshape(points, (np.shape(points)[0], -1)) # Collapse trialing dimensions
    assert points.shape[0] <= points.shape[1], "There are only {} points in {} dimensions.".format(points.shape[1], points.shape[0])
    ctr = points.mean(axis=1)
    x = points - ctr[:,np.newaxis]
    M = np.dot(x, x.T) # Could also use np.cov(x) here.
    return ctr, svd(M)[0][:,-1]

例：（10、100）で、x方向に薄く、y方向に100倍大きい2Dクラウドを作成します。

>>> pts = np.diag((.1, 10)).dot(randn(2,1000)) + np.reshape((10, 100),(2,-1))

フィット平面は（10、100）に非常に近く、法線はx軸に非常に近くなっています。

>>> planeFit(pts)

    (array([ 10.00382471,  99.48404676]),
     array([  9.99999881e-01,   4.88824145e-04]))
于 2013-09-23T20:44:35.860 に答える
13

最小二乗法は平面に簡単にフィットする必要があります。平面の方程式は次のとおりです。ax+by+ c=z。したがって、すべてのデータを使用して次のようなマトリックスを設定します。

    x_0   y_0   1  
A = x_1   y_1   1  
          ... 
    x_n   y_n   1

と

    a  
x = b  
    c

と

    z_0   
B = z_1   
    ...   
    z_n

言い換えると、Ax =Bです。ここで、係数であるxを解きます。ただし、ポイントが3つを超えるため、システムが過剰決定されているため、左の疑似逆行列を使用する必要があります。したがって、答えは次のとおりです。

a 
b = (A^T A)^-1 A^T B
c

そして、ここに例を含むいくつかの簡単なPythonコードがあります：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET  = 5
EXTENTS = 5
NOISE = 5

# create random data
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
    zs.append(xs[i]*TARGET_X_SLOPE + \
              ys[i]*TARGET_y_SLOPE + \
              TARGET_OFFSET + np.random.normal(scale=NOISE))

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
    tmp_A.append([xs[i], ys[i], 1])
    tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)

print("solution: %f x + %f y + %f = z" % (fit[0], fit[1], fit[2]))
print("errors:")
print(errors)
print("residual: {}".format(residual))

# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

あなたのポイントの解決策：

0.143509 x + 0.057196 y + 1.129595 = z

平面フィット

于 2017-06-01T19:21:14.810 に答える
12

平面にフィットしているという事実は、ここではわずかに関連しています。あなたがやろうとしていることは、推測から始めて特定の機能を最小化することです。そのためにscipy.optimize。これがグローバルに最適なソリューションであるという保証はなく、ローカルでのみ最適であることに注意してください。異なる初期条件が異なる結果に収束する可能性があります。これは、求めている極小値の近くから開始する場合にうまく機能します。

numpyのブロードキャストを利用して、自由にコードをクリーンアップしました。

import numpy as np

# coordinates (XYZ) of C1, C2, C4 and C5
XYZ = np.array([
        [0.274791784, -1.001679346, -1.851320839, 0.365840754],
        [-1.155674199, -1.215133985, 0.053119249, 1.162878076],
        [1.216239624, 0.764265677, 0.956099579, 1.198231236]])

# Inital guess of the plane
p0 = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

def f_min(X,p):
    plane_xyz = p[0:3]
    distance = (plane_xyz*X.T).sum(axis=1) + p[3]
    return distance / np.linalg.norm(plane_xyz)

def residuals(params, signal, X):
    return f_min(X, params)

from scipy.optimize import leastsq
sol = leastsq(residuals, p0, args=(None, XYZ))[0]

print("Solution: ", sol)
print("Old Error: ", (f_min(XYZ, p0)**2).sum())
print("New Error: ", (f_min(XYZ, sol)**2).sum())

これは与える：

Solution:  [  14.74286241    5.84070802 -101.4155017   114.6745077 ]
Old Error:  0.441513295404
New Error:  0.0453564286112
于 2012-09-06T13:48:23.887 に答える
3

外れ値を処理しながら（大きなデータセットがある場合）、svd以外の別の方法で解決策にすばやく到達する方法は、ransacです。

def fit_plane(voxels, iterations=50, inlier_thresh=10):  # voxels : x,y,z
    inliers, planes = [], []
    xy1 = np.concatenate([voxels[:, :-1], np.ones((voxels.shape[0], 1))], axis=1)
    z = voxels[:, -1].reshape(-1, 1)
    for _ in range(iterations):
        random_pts = voxels[np.random.choice(voxels.shape[0], voxels.shape[1] * 10, replace=False), :]
        plane_transformation, residual = fit_pts_to_plane(random_pts)
        inliers.append(((z - np.matmul(xy1, plane_transformation)) <= inlier_thresh).sum())
        planes.append(plane_transformation)
    return planes[np.array(inliers).argmax()]


def fit_pts_to_plane(voxels):  # x y z  (m x 3)
    # https://math.stackexchange.com/questions/99299/best-fitting-plane-given-a-set-of-points
    xy1 = np.concatenate([voxels[:, :-1], np.ones((voxels.shape[0], 1))], axis=1)
    z = voxels[:, -1].reshape(-1, 1)
    fit = np.matmul(np.matmul(np.linalg.inv(np.matmul(xy1.T, xy1)), xy1.T), z)
    errors = z - np.matmul(xy1, fit)
    residual = np.linalg.norm(errors)
    return fit, residual
于 2018-05-09T13:54:52.403 に答える
2

これにより、3D平面係数と近似のRMSEが返されます。

平面は同次座標表現で提供されます。つまり、点の同次座標との内積により、2つの間の距離が生成されます。

def fit_plane(points):
    assert points.shape[1] == 3
    centroid = points.mean(axis=0)
    x = points - centroid[None, :]
    U, S, Vt = np.linalg.svd(x.T @ x)
    normal = U[:, -1]
    origin_distance = normal @ centroid
    rmse = np.sqrt(S[-1] / len(points))
    return np.hstack([normal, -origin_distance]), rmse

マイナーノート：SVDは、外積行列の代わりにポイントに直接適用することもできますが、NumPyのSVD実装では遅くなることがわかりました。

U, S, Vt = np.linalg.svd(x.T, full_matrices=False)
rmse = S[-1] / np.sqrt(len(points))
于 2021-03-11T07:24:16.663 に答える
1

これが1つの方法です。ポイントがP[1].. P [n]の場合、これらの平均Mを計算し、それぞれから減算して、ポイントp [1] ..p[n]を取得します。次に、C = Sum {p [i] * p [i]'}（点の「共分散」行列）を計算します。次にCを対角線化します。つまり、C = U * E *U'となるように直交Uと対角線Eを見つけます。ポイントが実際に平面上にある場合、固有値の1つ（つまり、Eの対角要素）は非常に小さくなります（完全な算術では0になります）。いずれにせよ、これらのj番目の列が最小の場合は、Uのj番目の列を（A、B、C）とし、D =-M'*Nを計算します。これらのパラメータは、「最良の」平面を定義します。これは、P[]から平面までの距離の2乗の合計が最小になるようなものです。

于 2012-09-06T14:03:39.457 に答える