Index: /issm/trunk/src/c/objects/Bamg/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/QuadTree.cpp (revision 5223) +++ /issm/trunk/src/c/objects/Bamg/QuadTree.cpp (revision 5224) @@ -8,14 +8,39 @@ /*MACROS {{{1*/ + /* + * + * J j + * ^ ^ + * | | +--------+--------+ + * | | | | | + * 1X | | | 2 | 3 | + * | | | | | + * | | +--------+--------+ + * | | | | | + * 0X | | | 0 | 1 | + * | | | | | + * | | +--------+--------+ + * | +-----------------------> i + * | + * |----------------------------> I + * X0 X1 + * + * box 0 -> I=0 J=0 IJ=00 = 0 + * box 1 -> I=1 J=0 IJ=01 = 1 + * box 2 -> I=0 J=1 IJ=10 = 2 + * box 3 -> I=1 J=1 IJ=11 = 3 + */ #define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b))) #define ABS(i) ((i)<0 ?-(i) :(i)) #define MAX1(i,j) ((i)>(j) ?(i) :(j)) #define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2))) + //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) + //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) @@ -121,5 +146,4 @@ void QuadTree::Add(BamgVertex &w){ /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/Add)*/ - QuadTreeBox** pb=NULL; QuadTreeBox* b=NULL; @@ -134,5 +158,5 @@ pb = &root; - //Find the smallest box where w is located + /*Find the smallest box where w is located*/ while((b=*pb) && (b->nbitems<0)){ @@ -147,7 +171,7 @@ } - //OK, we have found b, a Subbox holding vertices (might be full) - //check that the vertex is not already in the box - if (b) { + /*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->nbitems > 3 && b->v[3] == &w) return; if (b->nbitems > 2 && b->v[2] == &w) return; @@ -156,14 +180,12 @@ } - //check that l is not 0 (this should not happen as MaxDeep = 30) - if (level==0){ - ISSMERROR("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 + /*check that l is not 0 (this should not happen as MaxDeep = 30)*/ + ISSMASSERT(level>0); + + /*Now, try to add the vertex, if the subbox is full (nbitems=4), we have to divide it + in 4 new subboxes*/ while ((b= *pb) && (b->nbitems == 4)){ // the QuadTreeBox is full - //Copy the 4 vertices in the current QuadTreebox + /*Copy the 4 vertices in the current QuadTreebox*/ BamgVertex* v4[4]; v4[0]= b->v[0]; @@ -172,37 +194,40 @@ v4[3]= b->v[3]; - //set n as negative (box full -> holds 4 pointers toward subboxes and not 4 vertices) + /*set nbitems as negative + * (box full -> holds 4 pointers toward subboxes and not 4 vertices)*/ b->nbitems = -b->nbitems; - //Initialize the 4 pointers toward the 4 subboxes + /*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 = 0010000 -> 0001000*/ 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,level)]; - - // alloc the QuadTreeBox + /*Put the four vertices in the new boxes*/ + for (int k=0;k<4;k++){ + + int ij; + /*bb is the new "sub"box of b where v4[k] is located*/ + 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 + /*Copy the current vertex*/ bb->v[bb->nbitems++] = v4[k]; } - //Get the subbox where w (i,j) is located + /*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 + /*alloc the QuadTreeBox if necessary*/ + if (!(b=*pb)) b=*pb= NewQuadTreeBox(); + + /*Add w*/ b->v[b->nbitems++]=&w; //Increase NbVertices by one (we have one new vertex) - NbVertices++; + NbVertices++; } /*}}}1*/ @@ -212,28 +237,30 @@ /*Build QuadTree*/ - QuadTreeBox* pb[MaxDeep]; - int pi[MaxDeep]; - Icoor1 ii[MaxDeep]; - Icoor1 jj[MaxDeep]; - register int level=0; // levelevelevel - register long n0; - register QuadTreeBox* b; - long h=MaxISize,h0; - long hb=MaxISize; - Icoor1 i0=0,j0=0; - //iplus= i projected in [0,MaxISize-1] (example: if i<0, i=0) - Icoor1 iplus( i=MaxISize, j=MaxISize-1) - Icoor1 jplus( j=MaxISize, j=MaxISize-1)*/ + Icoor1 jplus( jnbitems) return vn; - - //else, find the non empty QuadTreeBox containing the point (i,j) + /*if the tree is empty, return NULL pointer*/ + if (!root->nbitems) return nearest_v; + + /*else, find the non empty QuadTreeBox containing the point (i,j)*/ while((n0=b->nbitems)<0){ @@ -245,8 +272,8 @@ register QuadTreeBox* b0=b->b[k]; - // break if NULL box or empty + /* break if NULL box or empty*/ if (( b0 == NULL) || (b0->nbitems == 0)) break; - //Get next Qudtree box + /*Get next Qudtree box*/ b=b0; i0 += I_IJ(k,hb2); // i orign of QuadTreeBox (macro) @@ -255,16 +282,22 @@ } - // if the current subbox is holding vertices, we are almost done - 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; + + /*Compute norm with w*/ h0=NORM(iplus,i2.x,jplus,i2.y); + + /*is it smaller than previous value*/ if (h0v[k]; + nearest_v = b->v[k]; } } - return vn; + return nearest_v; } @@ -281,5 +314,5 @@ h=hb; - //loop, until level=0 + /*loop, until level=0*/ do { //get current box @@ -299,5 +332,5 @@ if (h0 v[k]; + nearest_v = b->v[k]; } } @@ -338,6 +371,6 @@ } while (level--); - //return vn, nearest vertex - return vn; + //return nearest_v, nearest vertex + return nearest_v; } /*}}}1*/ @@ -453,16 +486,16 @@ QuadTree::QuadTreeBox* QuadTree::NewQuadTreeBox(void){ - /*if firstbox==lastbox or firstbox>lastbox (we have reach the end of the StorageQuadTreeBox) + /*if currentbox==lastbox or currentbox>lastbox (we have reach the end of the StorageQuadTreeBox) * create a new StorageQuadTreeBox)*/ - if(!(sb->firstboxlastbox)){ + if(!(sb->currentboxlastbox)){ sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); } - ISSMASSERT(sb && sb->firstbox->nbitems==0); + ISSMASSERT(sb && sb->currentbox->nbitems==0); /*Increase counter*/ NbQuadTreeBox++; - /*firstbox now points toward next quadtree box*/ - return sb->firstbox++; + /*currentbox now points toward next quadtree box*/ + return sb->currentbox++; }/*}}}*/ /*FUNCTION QuadTree::SizeOf{{{1*/ @@ -569,5 +602,5 @@ } - firstbox = boxes; + currentbox = boxes; lastbox = boxes + nbquadtreeboxes; Index: /issm/trunk/src/c/objects/Bamg/QuadTree.h =================================================================== --- /issm/trunk/src/c/objects/Bamg/QuadTree.h (revision 5223) +++ /issm/trunk/src/c/objects/Bamg/QuadTree.h (revision 5224) @@ -39,5 +39,5 @@ /*Fields*/ - QuadTreeBox *boxes,*firstbox,*lastbox; + QuadTreeBox *boxes,*currentbox,*lastbox; long nbquadtreeboxes; StorageQuadTreeBox *nextsb; // next StorageQuadTreeBox