Changeset 5224
- Timestamp:
- 08/13/10 09:06:05 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Bamg
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Bamg/QuadTree.cpp
r5180 r5224 8 8 9 9 /*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 */ 10 33 #define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b))) 11 34 #define ABS(i) ((i)<0 ?-(i) :(i)) 12 35 #define MAX1(i,j) ((i)>(j) ?(i) :(j)) 13 36 #define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2))) 37 14 38 //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) 19 43 #define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )) 44 20 45 //I_IJ(k,l) returns l if first bit of k is 1, else 0 21 46 #define I_IJ(k,l) ((k&1) ? l:0) … … 121 146 void QuadTree::Add(BamgVertex &w){ 122 147 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/Add)*/ 123 124 148 QuadTreeBox** pb=NULL; 125 149 QuadTreeBox* b=NULL; … … 134 158 pb = &root; 135 159 136 / /Find the smallest box where w is located160 /*Find the smallest box where w is located*/ 137 161 while((b=*pb) && (b->nbitems<0)){ 138 162 … … 147 171 } 148 172 149 / /OK, we have found b, a Subbox holding vertices (might be full)150 //check that the vertex is not already in the box151 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){ 152 176 if (b->nbitems > 3 && b->v[3] == &w) return; 153 177 if (b->nbitems > 2 && b->v[2] == &w) return; … … 156 180 } 157 181 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*/ 165 187 while ((b= *pb) && (b->nbitems == 4)){ // the QuadTreeBox is full 166 188 167 / /Copy the 4 vertices in the current QuadTreebox189 /*Copy the 4 vertices in the current QuadTreebox*/ 168 190 BamgVertex* v4[4]; 169 191 v4[0]= b->v[0]; … … 172 194 v4[3]= b->v[3]; 173 195 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)*/ 175 198 b->nbitems = -b->nbitems; 176 199 177 / /Initialize the 4 pointers toward the 4 subboxes200 /*Initialize the 4 pointers toward the 4 subboxes*/ 178 201 b->b[0]=b->b[1]=b->b[2]=b->b[3]=NULL; 179 202 180 / / div the size by 2203 /*level = 0010000 -> 0001000*/ 181 204 level >>= 1; 182 205 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 189 214 if (!bb) bb=b->b[ij]=NewQuadTreeBox(); 190 215 191 / /Copy the 4 vertices216 /*Copy the current vertex*/ 192 217 bb->v[bb->nbitems++] = v4[k]; 193 218 } 194 219 195 / /Get the subbox where w (i,j) is located220 /*Get the subbox where w (i,j) is located*/ 196 221 pb = &b->b[IJ(i,j,level)]; 197 222 } 198 223 199 / / alloc the QuadTreeBox200 if (!(b =*pb)) b=*pb= NewQuadTreeBox();201 202 / /Add w224 /*alloc the QuadTreeBox if necessary*/ 225 if (!(b=*pb)) b=*pb= NewQuadTreeBox(); 226 227 /*Add w*/ 203 228 b->v[b->nbitems++]=&w; 204 229 205 230 //Increase NbVertices by one (we have one new vertex) 206 NbVertices++; 231 NbVertices++; 207 232 } 208 233 /*}}}1*/ … … 212 237 213 238 /*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)*/ 232 259 b = root; 233 260 234 / /if the tree is empty, return NULL pointer235 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)*/ 238 265 while((n0=b->nbitems)<0){ 239 266 … … 245 272 register QuadTreeBox* b0=b->b[k]; 246 273 247 / / break if NULL box or empty274 /* break if NULL box or empty*/ 248 275 if (( b0 == NULL) || (b0->nbitems == 0)) break; 249 276 250 / /Get next Qudtree box277 /*Get next Qudtree box*/ 251 278 b=b0; 252 279 i0 += I_IJ(k,hb2); // i orign of QuadTreeBox (macro) … … 255 282 } 256 283 257 / / if the current subbox is holding vertices, we are almost done258 if ( n0>0){284 /* if the current subbox is holding vertices, we are almost done*/ 285 if (n0>0){ 259 286 //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*/ 261 290 I2 i2=b->v[k]->i; 291 292 /*Compute norm with w*/ 262 293 h0=NORM(iplus,i2.x,jplus,i2.y); 294 295 /*is it smaller than previous value*/ 263 296 if (h0<h){ 264 297 h = h0; 265 vn= b->v[k];298 nearest_v = b->v[k]; 266 299 } 267 300 } 268 return vn;301 return nearest_v; 269 302 } 270 303 … … 281 314 h=hb; 282 315 283 / /loop, until level=0316 /*loop, until level=0*/ 284 317 do { 285 318 //get current box … … 299 332 if (h0 <h){ 300 333 h = h0; 301 vn= b->v[k];334 nearest_v = b->v[k]; 302 335 } 303 336 } … … 338 371 } while (level--); 339 372 340 //return vn, nearest vertex341 return vn;373 //return nearest_v, nearest vertex 374 return nearest_v; 342 375 } 343 376 /*}}}1*/ … … 453 486 QuadTree::QuadTreeBox* QuadTree::NewQuadTreeBox(void){ 454 487 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) 456 489 * create a new StorageQuadTreeBox)*/ 457 if(!(sb-> firstbox<sb->lastbox)){490 if(!(sb->currentbox<sb->lastbox)){ 458 491 sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); 459 492 } 460 ISSMASSERT(sb && sb-> firstbox->nbitems==0);493 ISSMASSERT(sb && sb->currentbox->nbitems==0); 461 494 462 495 /*Increase counter*/ 463 496 NbQuadTreeBox++; 464 497 465 /* firstbox now points toward next quadtree box*/466 return sb-> firstbox++;498 /*currentbox now points toward next quadtree box*/ 499 return sb->currentbox++; 467 500 }/*}}}*/ 468 501 /*FUNCTION QuadTree::SizeOf{{{1*/ … … 569 602 } 570 603 571 firstbox = boxes;604 currentbox = boxes; 572 605 lastbox = boxes + nbquadtreeboxes; 573 606 -
issm/trunk/src/c/objects/Bamg/QuadTree.h
r5180 r5224 39 39 40 40 /*Fields*/ 41 QuadTreeBox *boxes,* firstbox,*lastbox;41 QuadTreeBox *boxes,*currentbox,*lastbox; 42 42 long nbquadtreeboxes; 43 43 StorageQuadTreeBox *nextsb; // next StorageQuadTreeBox
Note:
See TracChangeset
for help on using the changeset viewer.