Changeset 5224


Ignore:
Timestamp:
08/13/10 09:06:05 (15 years ago)
Author:
Mathieu Morlighem
Message:

As usual

Location:
issm/trunk/src/c/objects/Bamg
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Bamg/QuadTree.cpp

    r5180 r5224  
    88
    99        /*MACROS {{{1*/
     10        /*
     11         *
     12         *    J    j
     13         *    ^    ^
     14         *    |    | +--------+--------+
     15         *    |    | |        |        |
     16         * 1X |    | |   2    |   3    |
     17         *    |    | |        |        |
     18         *    |    | +--------+--------+
     19         *    |    | |        |        |
     20         * 0X |    | |   0    |   1    |
     21         *    |    | |        |        |
     22         *    |    | +--------+--------+
     23         *    |    +-----------------------> i
     24         *    |         
     25         *    |----------------------------> I
     26         *              X0        X1 
     27         *
     28         * box 0 -> I=0 J=0 IJ=00  = 0
     29         * box 1 -> I=1 J=0 IJ=01  = 1
     30         * box 2 -> I=0 J=1 IJ=10  = 2
     31         * box 3 -> I=1 J=1 IJ=11  = 3
     32         */
    1033#define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b)))
    1134#define ABS(i) ((i)<0 ?-(i) :(i))
    1235#define MAX1(i,j) ((i)>(j) ?(i) :(j))
    1336#define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2)))
     37
    1438        //IJ(i,j,l) returns the box number of i and j with respect to l
    15         //if !j&l and !i$l -> 0 (box zero: lower left )
    16         //if !j&l and  i$l -> 1 (box one:  lower right)
    17         //if  j&l and !i$l -> 2 (box two:  upper left )
    18         //if  j&l and  i$l -> 3 (box three:upper right)
     39        //if !j&l and !i&l -> 0 (box zero: lower left )
     40        //if !j&l and  i&l -> 1 (box one:  lower right)
     41        //if  j&l and !i&l -> 2 (box two:  upper left )
     42        //if  j&l and  i&l -> 3 (box three:upper right)
    1943#define IJ(i,j,l)  ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 ))
     44
    2045        //I_IJ(k,l) returns l if first  bit of k is 1, else 0
    2146#define I_IJ(k,l)  ((k&1) ? l:0)
     
    121146        void  QuadTree::Add(BamgVertex &w){
    122147                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/Add)*/
    123 
    124148                QuadTreeBox** pb=NULL;
    125149                QuadTreeBox*  b=NULL;
     
    134158                pb = &root;
    135159
    136                 //Find the smallest box where w is located
     160                /*Find the smallest box where w is located*/
    137161                while((b=*pb) && (b->nbitems<0)){
    138162
     
    147171                }
    148172
    149                 //OK, we have found b, a Subbox holding vertices (might be full)
    150                 //check that the vertex is not already in the box
    151                 if  (b) {     
     173                /*OK, we have found b, a Subbox holding vertices (might be full)
     174                  check that the vertex is not already in the box*/
     175                if (b){     
    152176                        if (b->nbitems > 3 &&  b->v[3] == &w) return;
    153177                        if (b->nbitems > 2 &&  b->v[2] == &w) return;
     
    156180                }
    157181
    158                 //check that l is not 0 (this should not happen as MaxDeep = 30)
    159                 if (level==0){
    160                         ISSMERROR("level==0 cannot be true as it has been initialized as MaxISize = %i",MaxISize);
    161                 }
    162 
    163                 //Now, try to add the vertex, if the subbox is full (n=4), we have to divide it
    164                 //in 4 new subboxes
     182                /*check that l is not 0 (this should not happen as MaxDeep = 30)*/
     183                ISSMASSERT(level>0);
     184
     185                /*Now, try to add the vertex, if the subbox is full (nbitems=4), we have to divide it
     186                  in 4 new subboxes*/
    165187                while ((b= *pb) && (b->nbitems == 4)){ // the QuadTreeBox is full
    166188
    167                         //Copy the 4 vertices in the current QuadTreebox
     189                        /*Copy the 4 vertices in the current QuadTreebox*/
    168190                        BamgVertex* v4[4];
    169191                        v4[0]= b->v[0];
     
    172194                        v4[3]= b->v[3];
    173195
    174                         //set n as negative (box full -> holds 4 pointers toward subboxes and not 4 vertices)
     196                        /*set nbitems as negative
     197                         * (box full -> holds 4 pointers toward subboxes and not 4 vertices)*/
    175198                        b->nbitems = -b->nbitems;
    176199
    177                         //Initialize the 4 pointers toward the 4 subboxes
     200                        /*Initialize the 4 pointers toward the 4 subboxes*/
    178201                        b->b[0]=b->b[1]=b->b[2]=b->b[3]=NULL;
    179202
    180                         // div the size by 2
     203                        /*level = 0010000 -> 0001000*/
    181204                        level >>= 1;
    182205
    183                         //Put the four vertices in the new boxes
    184                         for (register int k=0;k<4;k++){
    185                                 register int ij;
    186                                 register QuadTreeBox* bb =  b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,level)];
    187 
    188                                 // alloc the QuadTreeBox
     206                        /*Put the four vertices in the new boxes*/
     207                        for (int k=0;k<4;k++){
     208
     209                                int          ij;
     210                                /*bb is the new "sub"box of b where v4[k] is located*/
     211                                QuadTreeBox *bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,level)];
     212
     213                                // alloc the QuadTreeBox
    189214                                if (!bb) bb=b->b[ij]=NewQuadTreeBox();
    190215
    191                                 //Copy the 4 vertices
     216                                /*Copy the current vertex*/
    192217                                bb->v[bb->nbitems++] = v4[k];
    193218                        }
    194219
    195                         //Get the subbox where w (i,j) is located
     220                        /*Get the subbox where w (i,j) is located*/
    196221                        pb = &b->b[IJ(i,j,level)];
    197222                }
    198223
    199                 //  alloc the QuadTreeBox
    200                 if (!(b = *pb)) b=*pb= NewQuadTreeBox();
    201 
    202                 //Add w
     224                /*alloc the QuadTreeBox if necessary*/
     225                if (!(b=*pb)) b=*pb= NewQuadTreeBox();
     226
     227                /*Add w*/
    203228                b->v[b->nbitems++]=&w;
    204229
    205230                //Increase NbVertices by one (we have one new vertex)
    206                 NbVertices++;   
     231                NbVertices++;
    207232        }
    208233        /*}}}1*/
     
    212237
    213238                /*Build QuadTree*/
    214                 QuadTreeBox* pb[MaxDeep];
    215                 int      pi[MaxDeep];
    216                 Icoor1   ii[MaxDeep];
    217                 Icoor1   jj[MaxDeep];
    218                 register int level=0; // levelevelevel
    219                 register long n0;
    220                 register QuadTreeBox* b;
    221                 long     h=MaxISize,h0;
    222                 long     hb=MaxISize;
    223                 Icoor1   i0=0,j0=0;
    224                 //iplus= i projected in [0,MaxISize-1] (example: if i<0, i=0)
    225                 Icoor1   iplus( i<MaxISize?(i<0?0:i):MaxISize-1);
    226                 //jplus= j projected in [0,MaxISize-1] (example: if j>=MaxISize, j=MaxISize-1)
    227                 Icoor1   jplus( j<MaxISize?(j<0?0:j):MaxISize-1);
    228                 //initial nearest vertex pointer
    229                 BamgVertex*  vn=NULL;
    230 
    231                 //Get initial Quadtree box (largest)
     239                QuadTreeBox          *pb[MaxDeep];
     240                int                   pi[MaxDeep];
     241                Icoor1                ii[MaxDeep];
     242                Icoor1                jj[MaxDeep];
     243                register int          level  = 0;
     244                register long         n0;
     245                register QuadTreeBox *b;
     246                long                  h      = MaxISize,h0;
     247                long                  hb     = MaxISize;
     248                Icoor1                i0     = 0,j0=0;
     249
     250                /*initial nearest vertex pointer*/
     251                BamgVertex*  nearest_v=NULL;
     252
     253                /*iplus= i projected in [0,MaxISize-1] (example: if i<0, i=0)*/
     254                Icoor1   iplus( i<MaxISize ? (i<0?0:i) : MaxISize-1);
     255                /*jplus= j projected in [0,MaxISize-1] (example: if j>=MaxISize, j=MaxISize-1)*/
     256                Icoor1   jplus( j<MaxISize ? (j<0?0:j) : MaxISize-1);
     257
     258                /*Get initial Quadtree box (largest)*/
    232259                b = root;
    233260
    234                 //if the tree is empty, return NULL pointer
    235                 if (!root->nbitems) return vn;
    236 
    237                 //else, find the non empty QuadTreeBox containing  the point (i,j)
     261                /*if the tree is empty, return NULL pointer*/
     262                if (!root->nbitems) return nearest_v;
     263
     264                /*else, find the non empty QuadTreeBox containing  the point (i,j)*/
    238265                while((n0=b->nbitems)<0){
    239266
     
    245272                        register QuadTreeBox* b0=b->b[k];
    246273
    247                         // break if NULL box or empty
     274                        /* break if NULL box or empty*/
    248275                        if (( b0 == NULL) || (b0->nbitems == 0)) break;
    249276
    250                         //Get next Qudtree box
     277                        /*Get next Qudtree box*/
    251278                        b=b0;   
    252279                        i0 += I_IJ(k,hb2); // i orign of QuadTreeBox (macro)
     
    255282                }
    256283
    257                 // if the current subbox is holding vertices, we are almost done
    258                 if ( n0>0 ){ 
     284                /* if the current subbox is holding vertices, we are almost done*/
     285                if (n0>0){ 
    259286                        //loop over the vertices of the box and find the closest vertex
    260                         for(register int k=0;k<n0;k++){
     287                        for(int k=0;k<n0;k++){
     288
     289                                /*get integer coordinates of current vertex*/
    261290                                I2 i2=b->v[k]->i;
     291
     292                                /*Compute norm with w*/
    262293                                h0=NORM(iplus,i2.x,jplus,i2.y);
     294
     295                                /*is it smaller than previous value*/
    263296                                if (h0<h){
    264297                                        h = h0;
    265                                         vn = b->v[k];
     298                                        nearest_v = b->v[k];
    266299                                }
    267300                        }
    268                         return vn;
     301                        return nearest_v;
    269302                }
    270303
     
    281314                h=hb;
    282315
    283                 //loop, until level=0
     316                /*loop, until level=0*/
    284317                do {
    285318                        //get current box
     
    299332                                        if (h0 <h){
    300333                                                h = h0;
    301                                                 vn = b->v[k];
     334                                                nearest_v = b->v[k];
    302335                                        }
    303336                                }
     
    338371                } while (level--);
    339372
    340                 //return vn, nearest vertex
    341                 return vn;
     373                //return nearest_v, nearest vertex
     374                return nearest_v;
    342375        }
    343376        /*}}}1*/
     
    453486        QuadTree::QuadTreeBox* QuadTree::NewQuadTreeBox(void){
    454487
    455                 /*if firstbox==lastbox or firstbox>lastbox (we have reach the end of the StorageQuadTreeBox)
     488                /*if currentbox==lastbox or currentbox>lastbox (we have reach the end of the StorageQuadTreeBox)
    456489                 * create a new StorageQuadTreeBox)*/
    457                 if(!(sb->firstbox<sb->lastbox)){
     490                if(!(sb->currentbox<sb->lastbox)){
    458491                        sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);
    459492                }
    460                 ISSMASSERT(sb && sb->firstbox->nbitems==0);
     493                ISSMASSERT(sb && sb->currentbox->nbitems==0);
    461494
    462495                /*Increase counter*/
    463496                NbQuadTreeBox++;
    464497
    465                 /*firstbox now points toward next quadtree box*/
    466                 return sb->firstbox++;
     498                /*currentbox now points toward next quadtree box*/
     499                return sb->currentbox++;
    467500        }/*}}}*/
    468501        /*FUNCTION QuadTree::SizeOf{{{1*/
     
    569602                }
    570603
    571                 firstbox = boxes;
     604                currentbox = boxes;
    572605                lastbox  = boxes + nbquadtreeboxes;
    573606
  • issm/trunk/src/c/objects/Bamg/QuadTree.h

    r5180 r5224  
    3939
    4040                                        /*Fields*/
    41                                         QuadTreeBox        *boxes,*firstbox,*lastbox;
     41                                        QuadTreeBox        *boxes,*currentbox,*lastbox;
    4242                                        long                nbquadtreeboxes;
    4343                                        StorageQuadTreeBox *nextsb;         // next StorageQuadTreeBox
Note: See TracChangeset for help on using the changeset viewer.