3D Mínimos Quadrados plano
-
05-07-2019 - |
Pergunta
O que é o algoritmo para calcular um avião dos mínimos quadrados em (x, y, z) espaço, dado um conjunto de pontos de dados em 3D? Em outras palavras, se I um grupo de pontos como (1, 2, 3), (4, 5, 6), (7, 8, 9), etc, como é que um vá sobre calcular o melhor ajuste plano f (x, y) = ax + by + c? Qual é o algoritmo para obter a, b, c e fora de um conjunto de pontos 3D?
Solução
Se tiver n pontos de dados (x [i], y [i], z [i]), calcular a 3x3 simétrica matriz A cujas entradas são:
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
Também computar o 3 elemento vector b:
{sum_i x[i]*z[i], sum_i y[i]*z[i], sum_i z[i]}
Em seguida, resolver Ax = b para o dado A e b. Os três componentes do vector da solução são os coeficientes para o plano de ajuste de mínimos quadrados {a, b, c}.
Note que este é o "mínimos quadrados ordinários" em forma, o que é apropriado somente quando z é esperado para ser uma função linear de x e y. Se você está procurando uma forma mais geral para um "melhor avião ajuste" em 3-espaço, você pode querer aprender sobre "geométricas" mínimos quadrados.
Note-se também que este irá falhar se os pontos estão em uma linha, como os seus pontos de exemplo são.
Outras dicas
a menos que alguém me diz como digitar equações aqui, deixe-me anotar os cálculos finais você tem que fazer:
pontos primeiramente, dado r_i \ n \ R, i = 1..N, calcular o centro de massa de todos os pontos:
r_G = \frac{\sum_{i=1}^N r_i}{N}
, em seguida, calcular o vetor normal n, que, em conjunto com o vector base r_G define o plano através do cálculo da matriz 3x3 A, tal como
A = \sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T
com esta matriz, o vetor normal n é agora dada pelo vector próprio de um correspondente para o valor próprio mínima de A.
Para saber mais sobre os autovetores / pares de valores próprios, use qualquer biblioteca de álgebra linear de sua escolha.
Esta solução baseia-se na Rayleight-Ritz Teorema para a matriz Hermitiana A.
Veja 'Mínimos Quadrados ajuste dos dados' por David Eberly para como eu vim acima com este para minimizar o ajuste geométrica (distância ortogonal de pontos para o avião).
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;
}
Timed a 37 micro segundos montagem de um avião para 1000 pontos (Windows 7, i7, o programa de 32 bits)
Tal como acontece com qualquer dos mínimos quadrados aproximar, você proceder assim:
Antes de iniciar a codificação
-
Anote uma equação para um avião em algum parametrização, dizem
0 = ax + by + z + d
em ti parâmetros(a, b, d)
. -
Encontre um
D(\vec{v};a, b, d)
expressão para a distância de um\vec{v}
ponto arbitrário. -
Anote o
S = \sigma_i=0,n D^2(\vec{x}_i)
soma, e simplificar até que seja expressa em termos de somas simples dos componentes dov
como\sigma v_x
,\sigma v_y^2
,\sigma v_x*v_z
... -
Anote o per parâmetro expressões de minimização
dS/dx_0 = 0
,dS/dy_0 = 0
... o que lhe dá um conjunto de três equações em três parâmetros e as somas da etapa anterior. -
Resolva este conjunto de equações para os parâmetros.
(ou para casos simples, basta olhar para cima o formulário). Usando um pacote de álgebra simbólica (como Mathematica) poderia tornar a vida muito mais fácil.
A codificação
- Escrever código para formar as somas necessárias e encontrar os parâmetros do último conjunto acima.
Alternativas
Note que se você realmente tinha apenas três pontos, você seria melhor apenas encontrar o avião que passa por eles.
Além disso, se a solução analítica em inviável (não é o caso de um avião, mas possível em geral) que você pode fazer os passos 1 e 2, e usar um Monte Carlo minimizador na soma no passo 3.
CGAL::linear_least_squares_fitting_3
Função linear_least_squares_fitting_3 calcula o 3D que melhor se adapta A linha ou plano (no sentido dos mínimos quadrados) de um conjunto de objectos, tais 3D como pontos, segmentos, triângulos, esferas, bolas, cubóides ou tetraedros.
A equação para um plano é: ax + by + c = z. Assim configurar matrizes como este com todos os seus dados:
x_0 y_0 1
A = x_1 y_1 1
...
x_n y_n 1
E
a
x = b
c
E
z_0
B = z_1
...
z_n
Em outras palavras: Ax = B. Agora resolva para x que são os seus coeficientes. Mas desde que (suponho) você tem mais de 3 pontos, o sistema é over-determinada de modo que você precisa para usar o inverso pseudo esquerda. Portanto, a resposta é:
a
b = (A^T A)^-1 A^T B
c
E aqui está algumas simples código Python com um exemplo:
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()
Isso reduz a total Menos problema Squares , que pode ser resolvido usando SVD decomposição.
código C ++ usando OpenCV:
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;
}
Tudo que você tem a fazer é resolver o sistema de equações.
Se esses são seus pontos: (1, 2, 3), (4, 5, 6), (7, 8, 9)
Isso dá-lhe as equações:
3=a*1 + b*2 + c
6=a*4 + b*5 + c
9=a*7 + b*8 + c
Assim, a pergunta realmente deve ser:? Como faço para resolver um sistema de equações
Portanto, eu recomendo a leitura esta pergunta SO .
Se eu tenha entendido mal a sua pergunta nos avise.
Editar :
Ignorar a minha resposta como você provavelmente quis dizer outra coisa.
Parece que tudo o que você quer fazer é regressão linear com 2 regressores. A wikipedia página sobre o assunto deve dizer-lhe tudo o que você precisa saber e então alguns.