Changeset 5300

Timestamp:

08/17/10 08:46:16 (15 years ago)

Author:

Mathieu Morlighem

Message:

As usual

File:

• issm/trunk/src/c/objects/Bamg/QuadTree.cpp (modified) (8 diffs)

issm/trunk/src/c/objects/Bamg/QuadTree.cpp

@@ -235 +235 @@
                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertex)*/
 
-                /*Build QuadTree*/
-                QuadTreeBox          *pb[MaxDeep];
-                int                   pi[MaxDeep];
-                Icoor1                ii[MaxDeep];
-                Icoor1                jj[MaxDeep];
-                register int          level;
-                register long         n0;
-                register QuadTreeBox *b;
-                long                  h      = MaxISize,h0;
-                long                  hb     = MaxISize;
-                Icoor1                i0     = 0,j0=0;
-
-                /*initial nearest vertex pointer*/
+                /*Intermediaries*/
+                QuadTreeBox *pb[MaxDeep];
+                int          pi[MaxDeep];
+                Icoor1       ii[MaxDeep];
+                Icoor1       jj[MaxDeep];
+                int          level;
+                long         n0;
+                QuadTreeBox *b;
+                long         h0;
+                long         h = MaxISize;
+                long         hb= MaxISize;
+                Icoor1       i0=0,j0=0;
+
+                /*initial output as NULL (no vertex found)*/
                 BamgVertex*  nearest_v=NULL;
 
-                /*iplus= i projected in [0,MaxISize-1] (example: if i<0, i=0)*/
-                Icoor1   iplus( i<MaxISize ? (i<0?0:i) : MaxISize-1);
-                /*jplus= j projected in [0,MaxISize-1] (example: if j>=MaxISize, j=MaxISize-1)*/
-                Icoor1   jplus( j<MaxISize ? (j<0?0:j) : MaxISize-1);
+                /*Project w coordinates (i,j) onto [0,MaxISize-1] x [0,MaxISize-1] -> (iplus,jplus)*/
+                Icoor1 iplus( i<MaxISize ? (i<0?0:i) : MaxISize-1);
+                Icoor1 jplus( j<MaxISize ? (j<0?0:j) : MaxISize-1);
 
                 /*Get initial Quadtree box (largest)*/
@@ -261 +261 @@
                 if (!root->nbitems) return nearest_v; 
 
-                /*else, find the non empty QuadTreeBox containing  the point (i,j)*/
+                /*else, find the smallest non-empty QuadTreeBox containing  the point (i,j)*/
                 while((n0=b->nbitems)<0){
 
-                        //shifted righ by one bit: hb2=01000000 -> 00100000
-                        register Icoor1 hb2 = hb >> 1;
-                        //Get QuadTreeBox number of size hb2 containing i;j (Macro)
-                        register int      k = IJ(iplus,jplus,hb2);
-                        //Get the corresponding box b0
-                        register QuadTreeBox* b0=b->b[k];
-
-                        /* break if NULL box or empty*/
+                        Icoor1       hb2 = hb >> 1;             //size of the current box
+                        int          k   = IJ(iplus,jplus,hb2); //box number (0,1,2 or 3)
+                        QuadTreeBox *b0  = b->b[k];             //pointer toward current box
+
+                        /* break if NULL box or empty (Keep previous box b)*/
                         if (( b0 == NULL) || (b0->nbitems == 0)) break;
 
-                        /*Get next Qudtree box*/
+                        /*Get next Quadtree box*/
                         b=b0;   
                         i0 += I_IJ(k,hb2); // i orign of QuadTreeBox (macro)
@@ -280 +277 @@
                         hb = hb2;          // size of the box (in Int)
                 }
+
+                /*The box b, is the smallest box containing the point (i,j) and
+                 * has the following properties:
+                 * - n0: number of items (>0 if vertices, else boxes)
+                 * - hb: box size (int)
+                 * - i0: x coordinate of the lower left corner
+                 * - j0: y coordinate of the lower left corner*/
 
                 /* if the current subbox is holding vertices, we are almost done*/
@@ -298 +302 @@
                                 }
                         }
+                        /*return closest vertex*/
                         return nearest_v;
                 }
 
-                /* general case: the current box is empty, we have to get backwards
+                /* general case: the current box is empty, we have to go backwards
                         and find the closest not-empty box and find the closest vertex*/
 
-                //initialize pb pi ii and jj
-                pb[0]=b;                  //pointer toward the box b
-                pi[0]=b->nbitems>0? (int)b->nbitems:4;//number of boxes in b
-                ii[0]=i0;                 // i coordinate of the box
-                jj[0]=j0;                 // j coordinate of the box
-
-                //initialize h as hb
+                /*Initialize search variables*/
+                pb[0]=b;                             //pointer toward the box b
+                pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b
+                ii[0]=i0;                            //i coordinate of the box lowest left corner
+                jj[0]=j0;                            //j coordinate of the box lowest left corner
+
+                /*initialize h: smallest box size, containing a vertex close to w*/
                 h=hb;
 
-                /*loop, until level=0 (main quadtree box)*/
+                /*Main loop*/
                 level=0;
                 do {
@@ -321 +326 @@
 
                         /*Loop over the items in current box (if not empty!)*/
-                        while (pi[level]--){
-
-                                /*k = number of items in the box*/
+                        while (pi[level]){
+
+                                /*We are looping now over the items of b. k is the current index (in [0 3])*/
+                                pi[level]--;
                                 int k=pi[level];
 
-                                /*if the current subbox is holding vertices*/
+                                /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/
                                 if (b->nbitems>0){
                                         I2 i2=b->v[k]->i;
@@ -335 +341 @@
                                         }
                                 }
-                                /*else: current box b is pointing toward 4 boxes (There must be a vertex at some point)*/
+                                /*else: current box b is pointing toward 4 boxes
+                                 * test sub-box k and go deeper into the tree if it is non empty
+                                 * and contains the point w modulo a size h that is either the size of the smallest
+                                 * non empty box containing w, or the closest point to w (so far) */
                                 else{
                                         QuadTreeBox* b0=b;
 
-                                        //if the next box exists:
+                                        /*if the next box exists:*/
                                         if (b=b->b[k]){
 
-                                                /*Get size (hb) and coordinates of the current sub-box*/
+                                                /*Get size (hb) and coordinates of the current sub-box lowest left corner*/
                                                 hb>>=1;
                                                 Icoor1 iii = ii[level]+I_IJ(k,hb);
                                                 Icoor1 jjj = jj[level]+J_IJ(k,hb);
 
-                                                /*if the current point (iplus,jplus) is in b,go to next box*/
+                                                /*if the current point (iplus,jplus) is in b (modulo h), this box is good:
+                                                 * it is holding vertices that are close to w */
                                                 if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){
-                                                        pb[++level]=  b;
+                                                        level++;
+                                                        pb[level]= b;
                                                         pi[level]= b->nbitems>0 ?(int)  b->nbitems : 4  ;
                                                         ii[level]= iii;
@@ -362 +373 @@
                                                 }
                                         }
-                                        /*Go backwards*/
                                         else{
+                                                /*Current box is NULL, go to next subbox of b (k=k-1)*/
                                                 b=b0;
                                         }
@@ -369 +380 @@
                         }
 
-                        //else go backwards
-                        //shifted righ by one bit: hb=001000000 -> 01000000
+                        /*We have found a vertex, now, let's try the other boxes of the previous level
+                         * in case there is a vertex closest to w that has not yet been tested*/
                         hb <<= 1;
                 } while (level--);