Changeset 2980

Timestamp:

02/08/10 11:00:45 (15 years ago)

Author:

Mathieu Morlighem

Message:

minor comments

Location:

issm/trunk/src/c/Bamgx

Files:

• QuadTree.h (modified) (1 diff)
• objects/QuadTree.cpp (modified) (6 diffs)

issm/trunk/src/c/Bamgx/QuadTree.h

@@ -7 +7 @@
         const int MaxDeep = 30;
         typedef  long  IntQuad;
-        const IntQuad MaxISize = ( 1L << MaxDeep); //long int 1, bitwise operation: 8L = 00001000 << 2L -> 00100000 shifted left by 2
+        //long int 1, bitwise operation: 8L = 00001000 << 2L -> 00100000 shifted left by 2
+        const IntQuad MaxISize = ( 1L << MaxDeep); 
         class Triangles;
         class Vertex;
         class QuadTree {
-                public:
+                private:
                         class QuadTreeBox { 
                                 public:

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

@@ -23 +23 @@
         //NORM(i1,j1,i2,j2) returns max(|i1-j1|,|i2-j2|)
 #define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2)))
-        //IJ(i,j,l) returns 3 if (1,1,1), 2 if (0,1,1), 1 if (1,0,1), else 0
-#define IJ(i,j,l) ( ( j & l) ? (( i & l) ? 3 : 2 ) :( ( i & l)? 1 : 0 ))
-        //I_IJ(k,l) retruns l if first  bit of k is 1, else 0
-#define I_IJ(k,l)  (( k&1) ? l : 0)
-        //J_IJ(k,l) retruns l if second bit of k is 1, else 0
-#define J_IJ(k,l)  (( k&2) ? l : 0)
+        //IJ(i,j,l) returns the box number of i and j with respect to l
+        //if !j&l and !i$l -> 0 (box zero: lower left )
+        //if !j&l and  i$l -> 1 (box one:  lower right)
+        //if  j&l and !i$l -> 2 (box two:  upper left )
+        //if  j&l and  i$l -> 3 (box three:upper right)
+#define IJ(i,j,l)  ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 ))
+        //I_IJ(k,l) returns l if first  bit of k is 1, else 0
+#define I_IJ(k,l)  ((k&1) ? l:0)
+        //J_IJ(k,l) returns l if second bit of k is 1, else 0
+#define J_IJ(k,l)  ((k&2) ? l:0)
+
+        /*What is a QuadTree?
+         * A Quadtree is a very simple way to group the vertices according
+         * to their location. A square that holds all the points of the mesh
+         * (or the geometry) is divided into 4 boxes. As soom as one box
+         * hold more than 4 vertices, it is divided into 4 new boxes, etc...
+         * There cannot be more than MAXDEEP (=30) subdivision.
+         * This process is like a Dichotomy in dimension 2
+         *
+         *  + - -  -    - -    -    - - + -   - + - + - + - -     - - +
+         *  |                           |       |   | X |             |
+         *                                      + - + - +
+         *  |                           |       |   |   |             |
+         *                              + -   - + - + - +             +
+         *  |                           |       |       |             |
+         *                         
+         *  |                           |       |       |             |
+         *  + - -  -    - -    -    - - + -   - + -   - + - -     - - +
+         *  |                           |               |             |
+         *                         
+         *  |                           |               |             |
+         *                         
+         *  |                           |               |             |
+         *  |                           |               |             |
+         *  + - -  -    - -    -    - - + -   -   -   - + - -     - - +
+         *  |                           |                             |
+         *                         
+         *  |                           |                             |
+         *                         
+         *  |                           |                             |
+         *                         
+         *  |                           |                             |
+         *  |                           |                             |
+         *  |                           |                             |
+         *  |                           |                             |
+         *  |                           |                             |
+         *  + - -  -    - -    -    - - + -   -   -   -   - -     - - +
+         *
+         * The coordinate system used in a quadtree are integers to avoid
+         * round-off errors. The vertex in the lower left box has the coordinates
+         * (0 0) 
+         * The upper right vertex has the follwing coordinates:
+         * 2^30 -1           2^30 -1        in decimal
+         * 0 1 1 1 .... 1    0 1 1 1 .... 1 in binary
+         *  \--   29  --/     \--   29  --/
+         * Using the binaries is therefor very easy to locate a vertex in a box:
+         * we just need to look at the bits from the left to the right (See ::Add)
+         */
 
         /*Constructors/Destructors*/
@@ -76 +128 @@
                 QuadTreeBox** pb;
                 QuadTreeBox*  b;
-                register long i=w.i.x, j=w.i.y,l=MaxISize;
+                register long i=w.i.x, j=w.i.y;
+                register long level=MaxISize;
+
+                //Get inital box (the largest)
                 pb = &root;
+
+                //Find the smallest box where w is located
                 while((b=*pb) && (b->n<0)){ 
-                        b->n--;                //n=n-1 in b->n
-                        l >>= 1;               //shifted righ by one bit: l=00000010 -> 00000001
-                        pb = &b->b[IJ(i,j,l)]; //pointer toward b
-                }
+
+                        //shift b->n by -1
+                        b->n--;
+
+                        //shifted righ by one bit: level=00000010 -> 00000001
+                        level >>= 1;
+
+                        //Get next subbox according to the bit value (level)
+                        pb = &b->b[IJ(i,j,level)];
+                }
+
+                //OK, we have found b, a Subbox holding vertices (might be full)
+                //check that the vertex is not already in the box
                 if  (b) {      
                         if (b->n > 3 &&  b->v[3] == &w) return;
@@ -90 +156 @@
                 }
 
-                //check that l is still non zero
-                if (l==0){
-                        throw ErrorException(__FUNCT__,exprintf("l==0 cannot be true as it has been initialized as MaxISize = %i",MaxISize));
-                }
+                //check that l is not 0 (this should not happen as MaxDeep = 30)
+                if (level==0){
+                        throw ErrorException(__FUNCT__,exprintf("level==0 cannot be true as it has been initialized as MaxISize = %i",MaxISize));
+                }
+
+                //Now, try to add the vertex, if the subbox is full (n=4), we have to divide it
+                //in 4 new subboxes
                 while ((b= *pb) && (b->n == 4)){ // the QuadTreeBox is full
-                        Vertex *v4[4]; // copy of the QuadTreeBox vertices
-
+
+                        //Copy the 4 vertices in the current QuadTreebox
+                        Vertex* v4[4];
                         v4[0]= b->v[0];
                         v4[1]= b->v[1];
                         v4[2]= b->v[2];
                         v4[3]= b->v[3];
-                        b->n = -b->n; // mark is pointer QuadTreeBox
-                        b->b[0]=b->b[1]=b->b[2]=b->b[3]=0; // set empty QuadTreeBox ptr
-                        l >>= 1;    // div the size by 2
-                        for (register int k=0;k<4;k++){ // for the 4 vertices find the sub QuadTreeBox ij
+
+                        //set n as negative (box full -> holds 4 pointers toward subboxes and not 4 vertices)
+                        b->n = -b->n;
+
+                        //Initialize the 4 pointers toward the 4 subboxes
+                        b->b[0]=b->b[1]=b->b[2]=b->b[3]=NULL;
+
+                        // div the size by 2
+                        level >>= 1;
+
+                        //Put the four vertices in the new boxes
+                        for (register int k=0;k<4;k++){
                                 register int ij;
-                                register QuadTreeBox * bb =  b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,l)];
-                                if (!bb) 
-                                 bb=b->b[ij]=NewQuadTreeBox(); // alloc the QuadTreeBox 
+                                register QuadTreeBox* bb =  b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,level)];
+
+                                // alloc the QuadTreeBox 
+                                if (!bb) bb=b->b[ij]=NewQuadTreeBox(); 
+
+                                //Copy the 4 vertices
                                 bb->v[bb->n++] = v4[k];
                         }
-                        pb = &b->b[IJ(i,j,l)];
-                }
-                if (!(b = *pb)) b=*pb= NewQuadTreeBox(); //  alloc the QuadTreeBox 
-                b->v[b->n++]=&w; // we add the vertex 
+
+                        //Get the subbox where w (i,j) is located
+                        pb = &b->b[IJ(i,j,level)];
+                }
+
+                //  alloc the QuadTreeBox 
+                if (!(b = *pb)) b=*pb= NewQuadTreeBox();
+
+                //Add w
+                b->v[b->n++]=&w;
+
+                //Increase NbVertices by one (we have one new vertex)
                 NbVertices++;    
         }
@@ -126 +215 @@
                 Icoor1   ii[MaxDeep];
                 Icoor1   jj[MaxDeep];
-                register int l=0; // level
+                register int level=0; // levelevelevel
+                register Int4 n0;
                 register QuadTreeBox* b;
                 IntQuad  h=MaxISize,h0;
@@ -138 +228 @@
                 Vertex*  vn=NULL;
 
-                //initialization for optimization
+                //Get initial Quadtree box (largest)
                 b = root;
-                register Int4 n0;
 
                 //if the tree is empty, return NULL pointer
                 if (!root->n) return vn; 
 
-                //else: find the index n of the non empty
-                //QuadTreeBox containing  the point (i,j)
-                while( (n0 = b->n) < 0){
-                        register Icoor1 hb2 = hb >> 1 ;
-                        register  int k = IJ(iplus,jplus,hb2);// QuadTreeBox number of size hb2 containing i;j (Macro)
-                        register QuadTreeBox* b0= b->b[k];
-                        if (( b0 == 0) || (b0->n == 0)) break;// null box or empty   => break       
+                //else, find the non empty QuadTreeBox containing  the point (i,j)
+                while((n0=b->n)<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
+                        if (( b0 == NULL) || (b0->n == 0)) break;
+
+                        //Get next Qudtree box
                         NbQuadTreeBoxSearch++;
                         b=b0;   
@@ -159 +255 @@
                 }
 
-                // if n0>0 ???
-                if ( n0 > 0){  
+                // if the current subbox is holding vertices, we are almost done
+                if ( n0>0 ){  
+                        //loop over the vertices of the box and find the closest vertex
                         for(register int k=0;k<n0;k++){
-                                I2 i2 =  b->v[k]->i;
-                                h0 = NORM(iplus,i2.x,jplus,i2.y);
-                                if (h0 <h) {
+                                I2 i2=b->v[k]->i;
+                                h0=NORM(iplus,i2.x,jplus,i2.y);
+                                if (h0<h){
                                         h = h0;
-                                        vn = b->v[k];}
-                                        NbVerticesSearch++;
+                                        vn = b->v[k];
+                                }
+                                NbVerticesSearch++;
                         }
                         return vn;
                 }
 
-                // general case
-                pb[0]= b;
-                pi[0]=b->n>0 ?(int)  b->n : 4  ;
-                ii[0]=i0;
-                jj[0]=j0;
+                /* general case: the current box is empty, we have to get 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->n>0? (int)b->n:4;//number of vertices in the box
+                ii[0]=i0;                 // i coordinate of the box
+                jj[0]=j0;                 // j coordinate of the box
+
+                //initialize h as hb
                 h=hb;
-                do {    
-                        b= pb[l];
-                        while (pi[l]--){              
-                                register int k = pi[l];
-
+
+                //loop, until level=0
+                do {
+                        //get current box
+                        b= pb[level];
+
+                        //Loop over the vertices in current box (if not empty!)
+                        while (pi[level]--){
+
+                                //k = number of vertices in the box if there are vertices
+                                //k = 4 if the current box is pointing toward 4 other boxes
+                                register int k=pi[level];
+
+                                //if the current subbox is holding vertices,
                                 if (b->n>0){ // Vertex QuadTreeBox not empty
                                         NbVerticesSearch++;
                                         I2 i2 =  b->v[k]->i;
                                         h0 = NORM(iplus,i2.x,jplus,i2.y);
-                                        if (h0 <h) 
-                                          {
+                                        if (h0 <h){
                                                 h = h0;
                                                 vn = b->v[k];
-                                          }
+                                        }
                                 }
-                                else{ // Pointer QuadTreeBox
-                                        register QuadTreeBox *b0=b;
+
+                                else{
+                                        register QuadTreeBox* b0=b;
                                         NbQuadTreeBoxSearch++;
+
+                                        //if the next box exists:
                                         if ((b=b->b[k])){
-                                                hb >>=1 ; // div by 2
-                                                register Icoor1 iii = ii[l]+I_IJ(k,hb);
-                                                register Icoor1 jjj = jj[l]+J_IJ(k,hb);
-
+                                                //shifted righ by one bit: hb2=01000000 -> 00100000
+                                                hb>>=1;
+                                                register Icoor1 iii = ii[level]+I_IJ(k,hb);
+                                                register Icoor1 jjj = jj[level]+J_IJ(k,hb);
+
+                                                //if the current point is in b,go to next box
                                                 if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){
-                                                        pb[++l]=  b;
-                                                        pi[l]= b->n>0 ?(int)  b->n : 4  ;
-                                                        ii[l]= iii;
-                                                        jj[l]= jjj;
+                                                        pb[++level]=  b;
+                                                        pi[level]= b->n>0 ?(int)  b->n : 4  ;
+                                                        ii[level]= iii;
+                                                        jj[level]= jjj;
                                                 }
+
+                                                //else go backwards
                                                 else{
-                                                        b=b0, hb <<=1 ;
+                                                        //shifted righ by one bit: hb=001000000 -> 01000000
+                                                        b=b0;
+                                                        hb<<=1;
                                                 }
                                         }
+                                        //Go backwards
                                         else b=b0;
                                 }
                         }
-                        hb <<= 1; // mul by 2 
-                } while (l--);
+
+                        //else go backwards
+                        //shifted righ by one bit: hb=001000000 -> 01000000
+                        hb <<= 1;
+                } while (level--);
+
+                //return vn, nearest vertex
                 return vn;
         }