3D FingerPrint - need to extract ridges and valleys?? Visual Studio 2010 and C++

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
...     ...     开发者_JAVA技巧 ...       

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();
}

Its part of my Thesis project. I need to do it using the 3D point cloud data. I don't have scanner with me. Its of 3rd party company and they just provide me with the 3D points. I have to work on this 3D points only.

Thanks,

NISHA

@Tom - Thanks. The polkadot3d shader of MeshLab was not accurate but it gave me rough idea of where ridges and valleys are. I think the ANN library is not giving me the correct neigbours to start with which leads to the wrong normals. But I am not sure how to fix this. As this is part of thesis, I and my professor came up with this algorithm to extract ridges and valleys. As per my research and other papers I read this method does work for extracting rigeys and valleys. I am just not getting it correctly in code :( I am quite sure that the method you have suggested will also work, but I may have to stick with my current algorithm and if it does not work at all should be able to say why it is not working! But, currently the problem seems to be with my code than the method I am using or I am missing some step here.

Have you considered turning the problem into an image analysis problem? If you the data is scanned by pressing the finger to a flat surface, then you can turn the scanning into an image where the gray level is determined by the z-depth of each point (distance between surface and point). You can then just paint the triangles by Gourad shading. You can probably do something similar by creating a 3D convex hull and then measure z-distances from it.

Once you have an image you can easily find ridges and valleys.