Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2908)
+++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp	(revision 2909)
@@ -492,6 +492,6 @@
 		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,6 +560,65 @@
 		}
 
-		//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,63 +632,62 @@
 			//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,7 +709,7 @@
 
 			// 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,5 +718,5 @@
 		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,6 +823,6 @@
 
 		/*clean up*/
-		delete []ev;
-		delete []hv;
+		delete []next_p;
+		delete []head_v;
 		delete []eangle;