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: 
[ n1 ] [u1*u1 u1*v1 v1*v1] [ a ] 
[ n2 ] = [u2*u2 u2*v2 v2*v2] [ b ] 
[... ] [ ... ... ... ] [ c ] 
[ nk ] [uk*uk uk*vk vk*vk] 

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

Replies to This Discussion

Permalink Reply by David Rutten on June 20, 2011 at 12:09pm

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

Permalink Reply by Nisha M on June 20, 2011 at 3:23pm

I have sent the file to your id

david@mcneel.com.

Permalink Reply by David Rutten on June 20, 2011 at 5:28pm

Got it. Looking into it now, but going to bed soon so don't hold out on an answer right away.

David Rutten

david@mcneel.com

Poprad, Slovakia

Permalink Reply by Nisha M on June 21, 2011 at 9:30pm

Hey David,

Did you get a chance to look into it?

Thanks,

Nisha

Permalink Reply by David Rutten on June 22, 2011 at 10:04am

Hi Nisha,

yes a bit. Nothing so far. I found a bunch of points away from the finger-print scan, I assume those can be safely ignored? My initial thought was to smooth the pointcloud in a smart way, then see which points are above and which below the smoothed mesh, this should then tell you which points belong to ridges and which to valleys.

But this is just a thought, I haven't tried to actually do this yet.

David Rutten

david@mcneel.com

Poprad, Slovakia

Permalink Reply by Daniel Hambleton on June 22, 2011 at 7:21am

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.

RSS