Grasshopper
3D FingerPrint - need to extract ridges and valleys?? Principal Curvatures Estimation
Hello,
I need to extract ridges and valleys of a 3D fingerprint. The output
should be a ply file which shows exactly where are the ridges and
valleys on 3d fingerprint using different colors.
**Input file** - ply file with just x,y,z locations. I got it from a
3d scanner.
Here is how first few lines of file looks like -
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
... ... ...
In case you need the data - please email me on
**nisha.m234@gmail.com**.
**Algorithm:**
I am trying to find principal curvatures for extracting the ridges and
valleys.
The steps I am following is:
1. Take a point x
2. Find its k nearest neighbors. I used k from 3 to 20.
3. average the k nearest neighbors => gives (_x, _y, _z)
4. compute covariance matrix
5. Now I take eigen values and eigen vectors of this covariance matrix
6. I get u, v and n here from eigen vectors.
u is a vector corresponding to largest eigen value
v corresponding to 2nd largest
n is 3rd smallest vector corresponding to smallest eigen value
7. Then for transforming the point(x,y,z) I compute matrix T
T =
[ui ] [u ] [x - _x]
[vi ] = [v ] x [y - _y]
[ni ] [n ] [z - _z]
8. for each i of the k nearest neighbors:<br>
[ n1 ] [u1*u1 u1*v1 v1*v1] [ a ]<br>
[ n2 ] = [u2*u2 u2*v2 v2*v2] [ b ] <br>
[... ] [ ... ... ... ] [ c ] <br>
[ nk ] [uk*uk uk*vk vk*vk]<br>
Solve this for a, b and c with least squares
9. this equations will give me a,b,c
10. now I compute eigen values of matrix
[a b
b a ]
11. This will give me 2 eigen values. one is Kmin and another Kmax.
**My Problem:**
The output is no where close to finding the correct Ridges and
Valleys. I am totally Stuck and frustrated. I am not sure where
exactly I am getting it wrong. I think the normal's are not computed
correctly. But I am not sure. I am very new to graphics programming
and so this maths, normals, shaders go way above my head. Any help
will be appreciated.
**PLEASE PLEASE HELP!!**
**Resources:**
I am using Visual Studio 2010 + Eigen Library + ANN Library.
**Other Options used**
I tried using MeshLab. I used ball pivoting triangles remeshing in
MeshLab and then applied the polkadot3d shader. If correctly
identifies the ridges and valleys. But I am not able to code it.
**My Function:**
//the function outputs to ply file
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.
I need to extract ridges and valleys of a 3D fingerprint. The output
should be a ply file which shows exactly where are the ridges and
valleys on 3d fingerprint using different colors.
**Input file** - ply file with just x,y,z locations. I got it from a
3d scanner.
Here is how first few lines of file looks like -
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
... ... ...
In case you need the data - please email me on
**nisha.m234@gmail.com**.
**Algorithm:**
I am trying to find principal curvatures for extracting the ridges and
valleys.
The steps I am following is:
1. Take a point x
2. Find its k nearest neighbors. I used k from 3 to 20.
3. average the k nearest neighbors => gives (_x, _y, _z)
4. compute covariance matrix
5. Now I take eigen values and eigen vectors of this covariance matrix
6. I get u, v and n here from eigen vectors.
u is a vector corresponding to largest eigen value
v corresponding to 2nd largest
n is 3rd smallest vector corresponding to smallest eigen value
7. Then for transforming the point(x,y,z) I compute matrix T
T =
[ui ] [u ] [x - _x]
[vi ] = [v ] x [y - _y]
[ni ] [n ] [z - _z]
8. for each i of the k nearest neighbors:<br>
[ n1 ] [u1*u1 u1*v1 v1*v1] [ a ]<br>
[ n2 ] = [u2*u2 u2*v2 v2*v2] [ b ] <br>
[... ] [ ... ... ... ] [ c ] <br>
[ nk ] [uk*uk uk*vk vk*vk]<br>
Solve this for a, b and c with least squares
9. this equations will give me a,b,c
10. now I compute eigen values of matrix
[a b
b a ]
11. This will give me 2 eigen values. one is Kmin and another Kmax.
**My Problem:**
The output is no where close to finding the correct Ridges and
Valleys. I am totally Stuck and frustrated. I am not sure where
exactly I am getting it wrong. I think the normal's are not computed
correctly. But I am not sure. I am very new to graphics programming
and so this maths, normals, shaders go way above my head. Any help
will be appreciated.
**PLEASE PLEASE HELP!!**
**Resources:**
I am using Visual Studio 2010 + Eigen Library + ANN Library.
**Other Options used**
I tried using MeshLab. I used ball pivoting triangles remeshing in
MeshLab and then applied the polkadot3d shader. If correctly
identifies the ridges and valleys. But I am not able to code it.
**My Function:**
//the function outputs to ply file
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();
}
Thanks,
NISHA
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();
}
Thanks,
NISHA
David Rutten
Hi Nisha,
can you post a Rhino file with the xyz point locations? I don't know much about Eigen vectors so I won't be able to help you there, but I might have some other ideas.
--
David Rutten
david@mcneel.com
Poprad, Slovakia
Jun 20, 2011
Daniel Hambleton
Meshlab is built on http://vcg.sourceforge.net/index.php/Main_Page (VCG Library). I don't know how many of their algorithms are accesible, but it might be worth a shot.
Jun 22, 2011