Posted on 2011-07-05 18:17
魔のkyo 阅读(3010)
评论(1) 编辑 收藏 引用 所属分类:
Math
首先一个维随机变量的均值记为E[X]
两个随机变量的协方差定义为cov(X,Y) = E[(X-E[X])(Y-E[Y])]
两个随机向量的协方差矩阵由随机向量分量间的协方差定义而成。
协方差矩阵元素Aij = cov(Xi,Yj)
协方差矩阵的阶=随机向量的维数，而与随机向量的个数无关。
协方差矩阵是对称矩阵。
其主对角线元素为点的各个分量的方差，可以反映点集在各个轴向的离散程度。
其他元素反映了不同分量的相关性，当两个分量距各自均值发生同向偏离时，协方差（即相应的协方差矩阵元素）趋向正值，当两个分量距各自均值发生异向偏离时，协方差趋向负值。
通过计算协方差矩阵的特征值和特征向量，可以反映出点集离散的主要方向。若特征向量为最大特征向量，则特征值即为沿着该轴的顶点数据具有的最大方差值；若特征向量为最小特征向量，则特征值即为沿着该轴的顶点数据具有的最小方差值。
特征值和特征向量的计算可以使用Jacobi方法，可以参考《实时碰撞检测算法技术》Page63~65，《Real-Time Collision Dectection》Page95~97
实用意义：
估计一个模型（一个点集）的最小包围球。
估计一个模型（一个点集）的最小OBB。
估计一个2D图形（二维点集）的最小面积包围矩形。
求2D点集的拟合直线。
求3D点集的拟合平面。
//3维向量协方差矩阵
void CovarianceMatrix(__out Matrix3f &covarianceMatrix, __in const Vector3f pt[], __in int numPts)
{
float invNum = 1.0f / (float)numPts;
Vector3f sum(0.0f);
for(int i=0;i<numPts;i++)
{
sum += pt[i];
}
Vector3f u = sum * invNum;
for(int i=0;i<3;i++)
{
for(int j=i;j<3;j++)
{
covarianceMatrix[i][j] = 0.0f;
for(int k=0;k<numPts;k++)
{
covarianceMatrix[i][j] += (pt[k][i] - u[i]) * (pt[k][j] - u[j]);
}
covarianceMatrix[i][j] *= invNum;
}
}
for(int i=0;i<3;i++)
{
for(int j=0;j<i;j++)
{
covarianceMatrix[i][j] = covarianceMatrix[j][i];
}
}
}
// 2-by-2 Symmetric Schur decomposition. Given an n-by-n symmetric matrix
// and indices p, q such that 1 <= p < q <= n, computes a sine-cosine pair
// (s, c) that will serve to form a Jacobi rotation matrix.
//
// See Golub, Van Loan, Matrix Computations, 3rd ed, p428
void SymSchur2(__in const Matrix3f &a, __in int p, __in int q, __out float* c, __out float* s)
{
if (abs(a[p][q]) > 0.0001f)
{
float r = (a[q][q] - a[p][p]) / (2.0f * a[p][q]);
float t;
if (r >= 0.0f)
t = 1.0f / (r + sqrt(1.0f + r*r));
else
t = -1.0f / (-r + sqrt(1.0f + r*r));
*c = 1.0f / sqrt(1.0f + t*t);
*s = t * *c;
}
else
{
*c = 1.0f;
*s = 0.0f;
}
}
// Computes the eigenvectors and eigenvalues of the symmetric matrix A using
// the classic Jacobi method of iteratively updating A as A = J∧T * A * J,
// where J = J(p, q, theta) is the Jacobi rotation matrix.
//
// On exit, v will contain the eigenvectors, and the diagonal elements
// of a are the corresponding eigenvalues.
//
// See Golub, Van Loan, Matrix Computations, 3rd ed, p428
void Jacobi(__inout Matrix3f &a, __out Matrix3f &v)
{
v.SetIdentity();
float prevoff;
// Repeat for some maximum number of iterations
const int MAX_ITERATIONS = 50;
for (int k = 0; k < MAX_ITERATIONS; k++)
{
// Find largest off-diagonal absolute element a[p][q]
int p = 0, q = 1;
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
if (i != j && abs(a[i][j]) > abs(a[p][q]))
{
p = i;
q = j;
}
}
}
// Compute the Jacobi rotation matrix J(p, q, theta)
// (This code can be optimized for the three different cases of rotation)
float c,s;
SymSchur2(a, p, q, &c, &s);
Matrix3f J;
J.SetIdentity();
J[p][p] = c; J[p][q] = s;
J[q][p] = -s; J[q][q] = c;
// Cumulate rotations into what will contain the eigenvectors
v = v * J;
// Make ’a’ more diagonal, until just eigenvalues remain on diagonal
a = (J.Transpose() * a) * J;
// Compute "norm" of off-diagonal elements
float off = 0.0f;
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
if (i != j)
{
off += a[i][j] * a[i][j];
}
}
}
/* off = sqrt(off); not needed for norm comparison */
// Stop when norm no longer decreasing
if (k > 2 && off >= prevoff)
return;
prevoff = off;
}
}
需要注意的是，
Jacobi(__inout Matrix3f
&a, __out Matrix3f &v);
参数a要求是对称矩阵。
函数返回时a是一个与传入矩阵相似的对角矩阵，v中的特征向量是按列存储的。
三个特征向量分别是
Vector3f(v[0][0],v[1][0],v[2][0])
Vector3f(v[0][1],v[1][1],v[2][1])
Vector3f(v[0][2],v[1][2],v[2][2])