Changeset 2909

Timestamp:

01/26/10 09:21:49 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added very useful comments

File:

• issm/trunk/src/c/Bamgx/objects/Geometry.cpp (modified) (6 diffs)

issm/trunk/src/c/Bamgx/objects/Geometry.cpp

@@ -492 +492 @@
                 Int4 i,j,k;
                 int jj;
-                Int4* hv=new Int4[nbv];
-                Int4* ev=new Int4[2*nbe];
+                Int4* head_v=new Int4[nbv];
+                Int4* next_p=new Int4[2*nbe];
                 float* eangle=new float[nbe];
                 double eps=1e-20;
@@ -560 +560 @@
                 }
 
-                //initialize hv as -1 for all vertices
-                for (i=0;i<nbv;i++) hv[i]=-1;
+                /* Here we use a powerful chaining algorithm
+                 *
+                 * 1. What is a chaining algorithm?
+                 *
+                 * If F is a function that goes from i in [0 n] to j in [0 m]
+                 * and we want to compute the reciprocal function F-1 of F
+                 * (what are the antecedents of a given j in Im(F) )
+                 * We use 2 lists:
+                 *    head_F[j] that holds the head of lists
+                 *    next_F[i] that holds the list of elements that have the same image
+                 *
+                 * Example:
+                 *    i1, i2, ..., ip in [0,n] are all antecedents of a given j in [0 m]
+                 *    head_F[j] = ip
+                 *    next_F[ip]= ip-1
+                 *    ....
+                 *    next_F[i2]= i1
+                 *    next_F[i1]= -1  //end of the list
+                 *
+                 * Algorithm:
+                 *    for(j=0;j<m;j++)  head_F[j] = -1 //initialization
+                 *    for(i=0;i<n;i++){
+                 *       j=F[i];
+                 *       next_F[i]= head_F[j];
+                 *       head_F[j]=i;
+                 *    }
+                 * 
+                 *    Then, we can go through all the elements that have for image j:
+                 *    for(i=head_F[j]; i!=-1; i=next_F[i])
+                 *    initialization of i by i=head_F[j]
+                 *    stop the loop when i=-1 (end of the chain)
+                 *    iterate using i=next_F[i] (next element that have for image j)
+                 * 
+                 * 2. How to use this algorithm here?
+                 * 
+                 * Here F is a function that associates two vertices v0 and v1 for a given edge E
+                 * We want to build the reciprocal function: what are the edges that contains
+                 * a vertex v?
+                 * To do so, we use the same chaining algorithm but there is a difficulty
+                 * coming from the fact that F we have a couple of vertices and not one 
+                 * vertices.
+                 * To overcome this difficulty, we use a global indices exactly like in 
+                 * C/C++ so that
+                 * a member of a 2-column-table can be described by one index p=i*2+j
+                 * i=(int)p/2 line number of p
+                 * j=p%2       column number of p
+                 *
+                 * Algorithm:
+                 *    for(i=0;i<nbv;i++)  head_v[i] = -1 //initialization
+                 *    for(i=0;i<nbe;i++){
+                 *       for(j=0;j<2;j++){
+                 *          p=2*i+j;
+                 *          v=edges(i,j);
+                 *          next_p[p]= head_v[v];
+                 *          head_v[v]=p;
+                 *       }
+                 *    }
+                 */
+
+                //initialize head_v as -1
+                for (i=0;i<nbv;i++) head_v[i]=-1;
                 k=0;
                 for (i=0;i<nbe;i++) {
@@ -573 +632 @@
                         //compute angle in [-Pi Pi]
                         eangle[i] = atan2(v10.y,v10.x);
-                        //build hv and ev
+                        //build chains head_v and next_p
                         for (j=0;j<2;j++){
                                 Int4 v=Number(edges[i].v[j]);
-                                ev[k]=hv[v];
-                                hv[v]=k++;
-                        }
-                }
-                // in our case, for 4 points and 4 edges (in parenthesis from now on):
-                // hv =  7 2 4 6 
-                // ev =  -1 -1 1 -1 3 -1 5 0 
-
-                //compute angle of edges
+                                next_p[k]=head_v[v];
+                                head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1;
+                        }
+                }
+
+                //sort head_v by order of increasing edges angle
                 for (i=0;i<nbv;i++) {
-
                         int exch=1, ord=0;      
+
+                        //exchange vertices position in head_v and next_p till tey are sorted
                         while (exch){
-                                Int4 *p=hv+i;                   // pointer to hv[vertex i]  (p=hv)
-                                Int4 *po=p;                     // copy of pointer p        (po=p)
-                                Int4  n=*p;                     // value hv[vertex i]       (n=7)
+                                Int4 *p=head_v+i;               // pointer toward head_v[vertex i]
+                                Int4 *po=p;                     // copy of pointer p
+                                Int4  n=*p;                     // next value of edge holding i
                                 register float angleold=-1000 ; // angle = - infinity
                                 ord=0; exch=0;
-                                while (n >=0){// loop as long as ...
+
+                                // loop over the edges that contain the vertex i
+                                while (n >=0){
                                         ord++;
-                                        register Int4  i1=n/2;       // i1 = floor (n/2)          (i1 = 7/2 = 3)
-                                        register Int4  j1=n%2;       // j1 = 1 if n is odd        (j1 = 1)
-                                        register Int4* pn=ev+n;      // pointer to ev[n]          (pn = &ev[7])
-                                        //compute angle from eangle:
-                                        // if n odd: angle = eangle - Pi [2*Pi]
-                                        // else    : angle = eangle
-                                        float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1];
-                                        n=*pn;                       //  n = ev[n]                (n = ev[7] = 0)
-                                        if (angleold > angle){       // exch to have : po -> pn -> p 
-                                                exch=1, *pn=*po, *po=*p, *p=n, po=pn;
+                                        register Int4  i1=n/2;       // i1 = floor (n/2)
+                                        register Int4  j1=n%2;       // j1 = 1 if n is odd
+                                        register Int4* pn=next_p+n;  // pointer to next_p[n]
+
+                                        //  n = next_p[n] = position in edge of next vertex i
+                                        n=*pn;                       
+
+                                        //compute angle between horizontal axis and v0->v1
+                                        float angle = j1 ? OppositeAngle(eangle[i1]):  eangle[i1]; 
+
+                                        //exchange if the current edge angle is smaller than the previous one
+                                        if (angleold > angle){
+                                                exch=1;
+                                                *pn=*po;  // next_p[n] = n + 1
+                                                *po=*p;   // 
+                                                *p=n;     // next_p[n+1] = n
+                                                po=pn;    // po now point toward pn (invert next and current)
                                         }
+
+                                        //else, continue going to the next edge position
                                         else{                        //  to have : po -> p -> pn
-                                                angleold=angle, po=p, p=pn;
+                                                angleold=angle; // update maximum angle
+                                                po=p;           // po now point toward p  (current position)
+                                                p=pn;           // p  now point toward pn (next position)
                                         }
                                 }
                         }
-                        //first iteration: i=0
-                        //hv 0 2 4 6 
-                        //ev 7 -1 1 -1 3 -1 5 -1 
-
-                        // angular test on current vertex to guess whether it is a corner
+                        printf("ord = %i\n",ord);
+
+                        // angular test on current vertex to guess whether it is a corner (ord = number of edges horlding i)
                         if(ord == 2) { 
-                                Int4  n1 = hv[i];
-                                Int4  n2 = ev[n1];
+                                Int4  n1 = head_v[i];
+                                Int4  n2 = next_p[n1];
                                 Int4  i1 = n1/2, i2 = n2/2; // edge number
                                 Int4  j1 = n1%2, j2 = n2%2; // vertex in the edge 
-                                /*//first iteration i=0
-                                //n1=hv[i]=0
-                                //n2=ev[i]=7
-                                //i1 = 0  i2 = 3
-                                //j1 = 0  j2 = 1
-                                printf("i=%i\n",i);
-                                printf("n1=hv[i]=%i\n",n1);
-                                printf("n2=ev[i]=%i\n",n2);
-                                printf("i1 = %i  i2 = %i\n",i1,i2);
-                                printf("j1 = %i  j2 = %i\n",j1,j2);*/
                                 float angle1=  j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
                                 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
@@ -651 +709 @@
 
                         // close the liste around the vertex 
-                        Int4 no=-1, ne = hv[i];
-                        while (ne >=0) ne = ev[no=ne];        
-                        if(no>=0) ev[no] = hv[i];
+                        Int4 no=-1, ne = head_v[i];
+                        while (ne >=0) ne = next_p[no=ne];        
+                        if(no>=0) next_p[no] = head_v[i];
                         // now the list around the vertex is circular
                 }
@@ -660 +718 @@
                 for (i=0;i<nbe;i++){
                         for (j=0;j<2;j++){
-                                Int4 n1 = ev[k++]; 
+                                Int4 n1 = next_p[k++]; 
                                 Int4 i1 = n1/2 ,j1=n1%2;
                                 if( edges[i1].v[j1] != edges[i].v[j]) {
@@ -765 +823 @@
 
                 /*clean up*/
-                delete []ev;
-                delete []hv;
+                delete []next_p;
+                delete []head_v;
                 delete []eangle;