Index: /issm/trunk/src/c/Bamgx/QuadTree.h =================================================================== --- /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2979) +++ /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2980) @@ -7,9 +7,10 @@ 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: Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2979) +++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2980) @@ -23,10 +23,62 @@ //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,11 +128,25 @@ 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,29 +156,52 @@ } - //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,5 +215,6 @@ 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,18 +228,24 @@ 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,61 +255,91 @@ } - // 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;kv[k]->i; - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k]->i; + h0=NORM(iplus,i2.x,jplus,i2.y); + if (h0v[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 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; }