https://stackoverflow.com/questions/1400213
-
05-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
3Dデータポイントのセットが与えられた場合、（x、y、z）空間で最小二乗平面を計算するアルゴリズムは何ですか？言い換えると、（1、2、3）、（4、5、6）、（7、8、9）などのような点がたくさんある場合、最適な平面fを計算するにはどうすればよいでしょうか（x、y）= ax + by + c？ 3Dポイントのセットからa、b、cを取得するアルゴリズムは何ですか？
解決
n個のデータポイント（x [i]、y [i]、z [i]）がある場合、エントリが次の3x3対称行列Aを計算します。
sum_i x[i]*x[i], sum_i x[i]*y[i], sum_i x[i]
sum_i x[i]*y[i], sum_i y[i]*y[i], sum_i y[i]
sum_i x[i], sum_i y[i], n
3要素ベクトルbも計算します：
{sum_i x[i]*z[i], sum_i y[i]*z[i], sum_i z[i]}
次に、与えられたAとbについてAx = bを解きます。解ベクトルの3つの成分は、最小二乗近似平面{a、b、c}の係数です。
これは「通常の最小二乗」であることに注意してください。 fitは、zがxおよびyの線形関数であると予想される場合にのみ適切です。より一般的に「最適な飛行機」を探している場合3空間では、「幾何学」について学びたいかもしれません。最小二乗。
また、例のポイントがそうであるように、ポイントがラインにある場合、これは失敗することに注意してください。
他のヒント
ここで方程式を入力する方法を教えてくれない限り、あなたがしなければならない最終的な計算を書き留めておきましょう：
まず、ポイントr_i \ n \ R、i = 1..Nを指定して、すべてのポイントの重心を計算します：
r_G = \frac{\sum_{i=1}^N r_i}{N}
次に、法線ベクトルnを計算します。これは、ベースベクトルr_Gとともに、3x3行列Aを次のように計算して平面を定義します
A = \sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T
この行列では、法線ベクトルnはAの最小固有値に対応するAの固有ベクトルによって与えられます。
固有ベクトル/固有値のペアについて調べるには、任意の線形代数ライブラリを使用します。
このソリューションは、エルミート行列Aのレイライト-リッツ定理に基づいています。
幾何学的な適合（点から平面までの直交距離）を最小にするためにこの方法を思いついた方法については、David Eberlyによる「データの最小二乗適合」を参照してください。
bool Geom_utils::Fit_plane_direct(const arma::mat& pts_in, Plane& plane_out)
{
bool success(false);
int K(pts_in.n_cols);
if(pts_in.n_rows == 3 && K > 2) // check for bad sizing and indeterminate case
{
plane_out._p_3 = (1.0/static_cast<double>(K))*arma::sum(pts_in,1);
arma::mat A(pts_in);
A.each_col() -= plane_out._p_3; //[x1-p, x2-p, ..., xk-p]
arma::mat33 M(A*A.t());
arma::vec3 D;
arma::mat33 V;
if(arma::eig_sym(D,V,M))
{
// diagonalization succeeded
plane_out._n_3 = V.col(0); // in ascending order by default
if(plane_out._n_3(2) < 0)
{
plane_out._n_3 = -plane_out._n_3; // upward pointing
}
success = true;
}
}
return success;
}
37ポイントで計時マイクロ秒平面を1000ポイントに合わせます（Windows 7、i7、32ビットプログラム）
最小二乗アプローチと同様に、次のように進みます。
コーディングを開始する前に
-
いくつかのパラメータ化で平面の方程式を書きます。たとえば、
0 = ax + by + z + dをeパラメータ（a、b、d）。 -
任意のポイント
\ vec {v}からの距離の式D（\ vec {v}; a、b、d）を見つけます。 -
合計
S = \ sigma_i = 0、n D ^ 2（\ vec {x} _i）を書き留め、単純な合計で表されるまで単純化します\ sigma v_x、\ sigma v_y ^ 2、\ sigma v_x * v_zなどのvのコンポーネント -
パラメータごとの最小化式を書き留めます
dS / dx_0 = 0、dS / dy_0 = 0... 3つのパラメーターと前のステップの合計。 -
パラメータのこの方程式のセットを解きます。
（または単純な場合は、フォームを検索してください）。シンボリック代数パッケージ（Mathematicaなど）を使用すると、作業がずっと楽になります。
コーディング
- 必要な合計を形成するコードを記述し、上記の最後のセットからパラメーターを見つけます。
代替案
実際に3つのポイントしか持っていない場合は、それらを通過する平面を見つけるだけでよいことに注意してください。
また、分析ソリューションが実行不可能な場合（飛行機の場合ではないが、一般的に可能）、ステップ1と2を実行し、ステップ3の合計のモンテカルロミニマイザー
CGAL::linear_least_squares_fitting_3
関数linear_least_squares_fitting_3は最適な3Dを計算します 次のような3Dオブジェクトのセットの線または平面（最小二乗の意味） ポイント、セグメント、三角形、球体、ボール、直方体または四面体として。
平面の方程式は、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:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print 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()
これは Total Least Squares の問題になります。これは SVD 分解。
OpenCVを使用したC ++コード：
float fitPlaneToSetOfPoints(const std::vector<cv::Point3f> &pts, cv::Point3f &p0, cv::Vec3f &nml) {
const int SCALAR_TYPE = CV_32F;
typedef float ScalarType;
// Calculate centroid
p0 = cv::Point3f(0,0,0);
for (int i = 0; i < pts.size(); ++i)
p0 = p0 + conv<cv::Vec3f>(pts[i]);
p0 *= 1.0/pts.size();
// Compose data matrix subtracting the centroid from each point
cv::Mat Q(pts.size(), 3, SCALAR_TYPE);
for (int i = 0; i < pts.size(); ++i) {
Q.at<ScalarType>(i,0) = pts[i].x - p0.x;
Q.at<ScalarType>(i,1) = pts[i].y - p0.y;
Q.at<ScalarType>(i,2) = pts[i].z - p0.z;
}
// Compute SVD decomposition and the Total Least Squares solution, which is the eigenvector corresponding to the least eigenvalue
cv::SVD svd(Q, cv::SVD::MODIFY_A|cv::SVD::FULL_UV);
nml = svd.vt.row(2);
// Calculate the actual RMS error
float err = 0;
for (int i = 0; i < pts.size(); ++i)
err += powf(nml.dot(pts[i] - p0), 2);
err = sqrtf(err / pts.size());
return err;
}
やらなければならないことは、連立方程式を解くことです。
それらがあなたのポイントである場合： （1、2、3）、（4、5、6）、（7、8、9）
これで方程式が得られます：
3=a*1 + b*2 + c
6=a*4 + b*5 + c
9=a*7 + b*8 + c
では、実際にあなたの質問は次のようになります。方程式系をどのように解くのですか？
したがって、これ SOの質問を読むことをお勧めします。
あなたの質問を誤解した場合はお知らせください。
編集：
おそらく何か他のものを意味していたので、私の答えを無視してください。
あなたがやりたいのは、2つの回帰変数を使った線形回帰だけです。件名のウィキペディアページは、あなたが知っておくべきすべてのこと、そしていくつかを教えてくれるはずです。
