我需要提取 3D 指纹的脊和谷。输出应该是一个 ply 文件，它使用不同的颜色准确显示 3d 指纹上的脊和谷在哪里。
输入文件- 只有 x、y、z 位置的层文件。我是从 3d 扫描仪中得到的。这是文件的前几行的样子 -
ply
format ascii 1.0
comment VCGLIB generated
element vertex 6183
property float x
property float y
property float z
end_header
-32.3271 -43.9859 11.5124
-32.0631 -43.983 11.4945
12.9266 -44.4913 28.2031
13.1701 -44.4918 28.2568
13.4138 -44.4892 28.2531
13.6581 -44.4834 28.1941
13.9012 -44.4851 28.2684
... ... ...
如果您需要数据 - 请发送电子邮件至nisha.m234@gmail.com给我。
算法： 我正在尝试找到用于提取山脊和山谷的主曲率。
我遵循的步骤是：
- 取点 x
- 找到它的 k 个最近邻。我使用了从 3 到 20 的 k。
- 平均 k 个最近的邻居 => 给出 (_x, _y, _z)
- 计算协方差矩阵
- 现在我取这个协方差矩阵的特征值和特征向量
- 我从特征向量中得到 u、v 和 n。u 是对应于最大特征值的向量 v 对应于第二大 n 是对应于最小特征值的第三小向量
- 然后为了转换 point(x,y,z) 我计算矩阵 T
T = | ui | |u | |x - _x| | vi | = |v | x |y - _y| | ni | |n | |z - _z| - 对于 k 个最近邻居中的每个 i：
用最小二乘法求解 a、b 和 c| n1 | |u1*u1 u1*v1 v1*v1| | a |
| n2 | = |u2*u2 u2*v2 v2*v2| | b |
|... | | ... ... ... | | c |
| nk | |uk*uk uk*vk vk*vk| - 这个方程会给我 a,b,c
- 现在我计算矩阵的特征值
| a b | | b a | - 这会给我2个特征值。一个是 Kmin，另一个是 Kmax。
我的问题： 输出远未找到正确的山脊和山谷。我完全陷入困境和沮丧。我不确定我到底在哪里弄错了。我认为法线的计算不正确。但我不确定。我对图形编程非常陌生，所以这些数学、法线、着色器超出了我的想象。任何帮助将不胜感激。 请帮忙！！
资源： 我正在使用 Visual Studio 2010 + Eigen 库 + ANN 库。
使用的其他选项 我尝试使用 MeshLab。我在 MeshLab 中使用了球旋转三角形重新网格化，然后应用了 polkadot3d 着色器。如果正确识别山脊和山谷。但我无法对其进行编码。
My Function: //函数输出到ply文件
void getEigen()
{
int nPts; // actual number of data points
ANNpointArray dataPts; // data points
ANNpoint queryPt; // query point
ANNidxArray nnIdx;// near neighbor indices
ANNdistArray dists; // near neighbor distances
ANNkd_tree* kdTree; // search structure
//for k = 25 and esp = 2, seems to got few ridges
queryPt = annAllocPt(dim); // allocate query point
dataPts = annAllocPts(maxPts, dim); // allocate data points
nnIdx = new ANNidx[k]; // allocate near neigh indices
dists = new ANNdist[k]; // allocate near neighbor dists
nPts = 0; // read data points
ifstream dataStream;
dataStream.open(inputFile, ios::in);// open data file
dataIn = &dataStream;
ifstream queryStream;
queryStream.open("input/query.pts", ios::in);// open data file
queryIn = &queryStream;
while (nPts < maxPts && readPt(*dataIn, dataPts[nPts])) nPts++;
kdTree = new ANNkd_tree( // build search structure
dataPts, // the data points
nPts, // number of points
dim); // dimension of space
while (readPt(*queryIn, queryPt)) // read query points
{
kdTree->annkSearch( // search
queryPt, // query point
k, // number of near neighbors
nnIdx, // nearest neighbors (returned)
dists, // distance (returned)
eps); // error bound
double x = queryPt[0];
double y = queryPt[1];
double z = queryPt[2];
double _x = 0.0;
double _y = 0.0;
double _z = 0.0;
#pragma region Compute covariance matrix
for (int i = 0; i < k; i++)
{
_x += dataPts[nnIdx[i]][0];
_y += dataPts[nnIdx[i]][1];
_z += dataPts[nnIdx[i]][2];
}
_x = _x/k; _y = _y/k; _z = _z/k;
double A[3][3] = {0,0,0,0,0,0,0,0,0};
for (int i = 0; i < k; i++)
{
double X = dataPts[nnIdx[i]][0];
double Y = dataPts[nnIdx[i]][1];
double Z = dataPts[nnIdx[i]][2];
A[0][0] += (X-_x) * (X-_x);
A[0][1] += (X-_x) * (Y-_y);
A[0][2] += (X-_x) * (Z-_z);
A[1][0] += (Y-_y) * (X-_x);
A[1][1] += (Y-_y) * (Y-_y);
A[1][2] += (Y-_y) * (Z-_z);
A[2][0] += (Z-_z) * (X-_x);
A[2][1] += (Z-_z) * (Y-_y);
A[2][2] += (Z-_z) * (Z-_z);
}
MatrixXd C(3,3);
C <<A[0][0]/k, A[0][1]/k, A[0][2]/k,
A[1][0]/k, A[1][1]/k, A[1][2]/k,
A[2][0]/k, A[2][1]/k, A[2][2]/k;
#pragma endregion
EigenSolver<MatrixXd> es(C);
MatrixXd Eval = es.eigenvalues().real().asDiagonal();
MatrixXd Evec = es.eigenvectors().real();
MatrixXd u,v,n;
double a = Eval.row(0).col(0).value();
double b = Eval.row(1).col(1).value();
double c = Eval.row(2).col(2).value();
#pragma region SET U V N
if(a>b && a>c)
{
u = Evec.row(0);
if(b>c)
{ v = Eval.row(1); n = Eval.row(2);}
else
{ v = Eval.row(2); n = Eval.row(1);}
}
else
if(b>a && b>c)
{
u = Evec.row(1);
if(a>c)
{ v = Eval.row(0); n = Eval.row(2);}
else
{ v = Eval.row(2); n = Eval.row(0);}
}
else
{
u = Eval.row(2);
if(a>b)
{ v = Eval.row(0); n = Eval.row(1);}
else
{ v = Eval.row(1); n = Eval.row(0);}
}
#pragma endregion
MatrixXd O(3,3);
O <<u,
v,
n;
MatrixXd UV(k,3);
VectorXd N(k,1);
for( int i=0; i<k; i++)
{
double x = dataPts[nnIdx[i]][0];;
double y = dataPts[nnIdx[i]][1];;
double z = dataPts[nnIdx[i]][2];;
MatrixXd X(3,1);
X << x-_x,
y-_y,
z-_z;
MatrixXd T = O * X;
double ui = T.row(0).col(0).value();
double vi = T.row(1).col(0).value();
double ni = T.row(2).col(0).value();
UV.row(i) << ui * ui, ui * vi, vi * vi;
N.row(i) << ni;
}
Vector3d S = UV.colPivHouseholderQr().solve(N);
MatrixXd II(2,2);
II << S.row(0).value(), S.row(1).value(),
S.row(1).value(), S.row(2).value();
EigenSolver<MatrixXd> es2(II);
MatrixXd Eval2 = es2.eigenvalues().real().asDiagonal();
MatrixXd Evec2 = es2.eigenvectors().real();
double kmin, kmax;
if(Eval2.row(0).col(0).value() < Eval2.row(1).col(1).value())
{
kmin = Eval2.row(0).col(0).value();
kmax = Eval2.row(1).col(1).value();
}
else
{
kmax = Eval2.row(0).col(0).value();
kmin = Eval2.row(1).col(1).value();
}
double thresh = 0.0020078;
if (kmin < thresh && kmax > thresh )
cout << x << " " << y << " " << z << " "
<< 255 << " " << 0 << " " << 0
<< endl;
else
cout << x << " " << y << " " << z << " "
<< 255 << " " << 255 << " " << 255
<< endl;
}
delete [] nnIdx;
delete [] dists;
delete kdTree;
annClose();
}
它是我论文项目的一部分。我需要使用 3D 点云数据来完成。我没有随身携带扫描仪。它是第 3 方公司的，他们只是为我提供 3D 点。我只需要处理这个 3D 点。
谢谢，
尼沙
@Tom - 谢谢。MeshLab 的 polkadot3d 着色器并不准确，但它让我大致了解了山脊和山谷的位置。我认为 ANN 库没有给我正确的邻居，这会导致错误的法线。但我不知道如何解决这个问题。由于这是论文的一部分，我和我的教授提出了这个算法来提取山脊和山谷。根据我的研究和其他论文，我读到这种方法确实适用于提取 rigeys 和 valleys。我只是在代码中没有正确理解它:(我很确定你建议的方法也可以，但我可能必须坚持我目前的算法，如果它根本不起作用应该能够说出原因不工作！但是，目前问题似乎出在我的代码上，而不是我使用的方法，或者我在这里遗漏了一些步骤。