Index: /issm/trunk/src/c/objects/Bamg/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/QuadTree.cpp (revision 5299) +++ /issm/trunk/src/c/objects/Bamg/QuadTree.cpp (revision 5300) @@ -235,23 +235,23 @@ /*Original code from Frederic Hecht (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, j=MaxISize-1)*/ - Icoor1 jplus( j (iplus,jplus)*/ + Icoor1 iplus( inbitems) 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,4 +277,11 @@ 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,20 +302,21 @@ } } + /*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,10 +326,11 @@ /*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,19 +341,24 @@ } } - /*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,6 +373,6 @@ } } - /*Go backwards*/ else{ + /*Current box is NULL, go to next subbox of b (k=k-1)*/ b=b0; } @@ -369,6 +380,6 @@ } - //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--);