3D 최소 제곱 평면
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를 해결하십시오. 용액 벡터의 세 가지 구성 요소는 최소 제곱 피트 평면 {A, B, C}에 대한 계수입니다.
z가 x와 y의 선형 함수 일 때만 적절한 "일반 최소 제곱"적합도입니다. 3 공간에서 "가장 적합한 비행기"를 더 일반적으로 찾고 있다면 "기하학적"최소 제곱에 대해 배우고 싶을 수도 있습니다.
또한 예제 포인트와 같이 포인트가 한 줄에 있으면 실패합니다.
다른 팁
누군가가 여기에 방정식을 입력하는 방법을 알려주지 않는 한, 당신이해야 할 최종 계산을 기록하겠습니다.
첫째, 지점이 주어진 포인트 R_i n r, i = 1..n, 모든 지점의 질량 중심을 계산합니다.
r_G = \frac{\sum_{i=1}^N r_i}{N}
그런 다음 기본 벡터 R_G와 함께 3x3 행렬 A를 다음과 같이 계산하여 평면을 정의하는 일반 벡터 N을 계산합니다.
A = \sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T
이 매트릭스를 사용하여, 정상 벡터 N은 이제 A의 최소 고유 값에 해당하는 A의 고유 벡터에 의해 주어진다.
고유 벡터/고유 값 쌍에 대해 알아 보려면 선택한 선형 대수 라이브러리를 사용하십시오.
이 솔루션은 Hermitian Matrix A의 Rayleight-Ritz 정리를 기반으로합니다.
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당신의 매개 변수에서(a, b, d).표현을 찾으십시오
D(\vec{v};a, b, d)임의의 지점에서의 거리\vec{v}.합계를 기록하십시오
S = \sigma_i=0,n D^2(\vec{x}_i), 그리고 구성 요소의 단순한 합으로 표현 될 때까지 단순화합니다.v처럼\sigma v_x,\sigma v_y^2,\sigma v_x*v_z...매개 변수 최소화 표현식을 기록하십시오
dS/dx_0 = 0,dS/dy_0 = 0... 3 개의 매개 변수와 이전 단계의 합계로 세 가지 방정식 세트를 제공합니다.매개 변수에 대한이 방정식 세트를 해결하십시오.
(또는 간단한 경우에는 양식을 찾아보십시오). 상징적 대수 패키지 (Mathematica와 같은)를 사용하면 인생이 훨씬 쉬워 질 수 있습니다.
코딩
- 코드를 작성하여 필요한 합을 형성하고 위의 마지막 세트에서 매개 변수를 찾으십시오.
대안
실제로 3 점 밖에되지 않으면 비행기를 통과하는 비행기를 찾는 것이 좋습니다.
또한 불가능한 분석 솔루션 (평면의 경우는 아니지만 일반적으로 가능)에서 1 단계와 2 단계를 수행하고 Monte Carlo Minimizer 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
예를 들어 간단한 파이썬 코드가 있습니다.
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()
이것은로 감소합니다 총 최소 제곱 문제를 해결할 수 있습니다 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
그래서 당신의 질문은 실제로 : 방정식 시스템을 어떻게 해결합니까?
그러므로 나는 읽는 것이 좋습니다 이것 그래서 질문.
귀하의 질문을 오해했다면 알려주십시오.
편집하다:
아마 당신이 다른 것을 의미했을 때 내 대답을 무시하십시오.
당신이 원하는 것은 2 개의 회귀기가있는 선형 회귀뿐입니다. 그만큼 Wikipedia 페이지 주제는 당신이 알아야 할 모든 것을 말해야합니다.
