![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21623 for issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp
- Timestamp:
- 03/23/17 14:07:02 (8 years ago)
- File:
-
- 1 edited
-
issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/bamg/BamgQuadtree.cpp
r18064 r21623 7 7 8 8 namespace bamg { 9 10 /*MACROS {{{*/11 /*12 *13 * J j14 * ^ ^15 * | | +--------+--------+16 * | | | | |17 * 1X | | | 2 | 3 |18 * | | | | |19 * | | +--------+--------+20 * | | | | |21 * 0X | | | 0 | 1 |22 * | | | | |23 * | | +--------+--------+24 * | +-----------------------> i25 * |26 * |----------------------------> I27 * X0 X128 *29 * box 0 -> I=0 J=0 IJ=00 = 030 * box 1 -> I=1 J=0 IJ=01 = 131 * box 2 -> I=0 J=1 IJ=10 = 232 * box 3 -> I=1 J=1 IJ=11 = 333 */34 #define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b)))35 9 #define ABS(i) ((i)<0 ?-(i) :(i)) 36 #define MAX1(i,j) ((i)>(j) ?(i) :(j))37 #define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2)))38 39 //IJ(i,j,l) returns the box number of i and j with respect to l40 //if !j&l and !i&l -> 0 (box zero: lower left )41 //if !j&l and i&l -> 1 (box one: lower right)42 //if j&l and !i&l -> 2 (box two: upper left )43 //if j&l and i&l -> 3 (box three:upper right)44 #define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 ))45 46 //I_IJ(k,l) returns l if first bit of k is 1, else 047 #define I_IJ(k,l) ((k&1) ? l:0)48 //J_IJ(k,l) returns l if second bit of k is 1, else 049 #define J_IJ(k,l) ((k&2) ? l:0)50 /*}}}*/51 10 /*DOCUMENTATION What is a BamgQuadtree? {{{ 52 11 * A Quadtree is a very simple way to group vertices according … … 100 59 BamgQuadtree::BamgQuadtree(){/*{{{*/ 101 60 102 /*Number of boxes and vertices*/ 103 NbBamgQuadtreeBox=0; 104 NbVertices=0; 61 /*Initialize fields*/ 62 this->MaxDepth = 30; 63 this->MaxISize = 1073741824; //2^30 64 this->MaxICoord = 1073741823; //2^30 - 1 65 this->NbQuadtreeBox = 0; 66 this->NbVertices = 0; 105 67 106 68 /*Create container*/ 107 boxcontainer=new DataSet();69 this->boxcontainer=new DataSet(); 108 70 109 71 /*Create Root, pointer toward the main box*/ 110 root=NewBamgQuadtreeBox(); 111 112 } 113 /*}}}*/ 72 this->root=NewBamgQuadtreeBox(); 73 74 } /*}}}*/ 114 75 BamgQuadtree::BamgQuadtree(Mesh * t,long nbv){ /*{{{*/ 115 76 116 /*Number of boxes and vertices*/ 117 NbBamgQuadtreeBox=0; 118 NbVertices=0; 77 /*Initialize fields*/ 78 this->MaxDepth = 30; 79 this->MaxISize = 1073741824; //2^30 80 this->MaxICoord = 1073741823; //2^30 - 1 81 this->NbQuadtreeBox = 0; 82 this->NbVertices = 0; 119 83 120 84 /*Create container*/ 121 boxcontainer=new DataSet();85 this->boxcontainer=new DataSet(); 122 86 123 87 /*Create Root, pointer toward the main box*/ 124 root=NewBamgQuadtreeBox();88 this->root=NewBamgQuadtreeBox(); 125 89 126 90 /*Check Sizes*/ 127 _assert_( MaxISize>MaxICoor);91 _assert_(this->MaxISize>this->MaxICoord); 128 92 129 93 /*Add all vertices of the mesh*/ 130 if (nbv==-1) nbv=t->nbv; 131 for (int i=0;i<nbv;i++) Add(t->vertices[i]); 94 if(nbv==-1) nbv=t->nbv; 95 for(int i=0;i<nbv;i++){ 96 this->Add(t->vertices[i]); 97 } 132 98 133 99 } … … 140 106 141 107 /*Methods*/ 142 void BamgQuadtree::Add(BamgVertex &w){/*{{{*/108 void BamgQuadtree::Add(BamgVertex &w){/*{{{*/ 143 109 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/Add)*/ 144 110 BamgQuadtreeBox** pb=NULL; … … 164 130 165 131 //Get next subbox according to the bit value (level) 166 pb = &b->b [IJ(i,j,level)];132 pb = &b->box[BoxNumber(i,j,level)]; 167 133 } 168 134 … … 195 161 196 162 /*Initialize the 4 pointers toward the 4 subboxes*/ 197 b->b [0]=b->b[1]=b->b[2]=b->b[3]=NULL;163 b->box[0]=b->box[1]=b->box[2]=b->box[3]=NULL; 198 164 199 165 /*level = 0010000 -> 0001000*/ … … 205 171 int ij; 206 172 /*bb is the new "sub"box of b where v4[k] is located*/ 207 BamgQuadtreeBox *bb = b->b [ij=IJ(v4[k]->i.x,v4[k]->i.y,level)];173 BamgQuadtreeBox *bb = b->box[ij=BoxNumber(v4[k]->i.x,v4[k]->i.y,level)]; 208 174 209 175 // alloc the BamgQuadtreeBox 210 if (!bb) bb=b->b [ij]=NewBamgQuadtreeBox();176 if (!bb) bb=b->box[ij]=NewBamgQuadtreeBox(); 211 177 212 178 /*Copy the current vertex*/ … … 215 181 216 182 /*Get the subbox where w (i,j) is located*/ 217 pb = &b->b [IJ(i,j,level)];183 pb = &b->box[BoxNumber(i,j,level)]; 218 184 } 219 185 … … 228 194 } 229 195 /*}}}*/ 230 BamgVertex* BamgQuadtree::NearestVertex(Icoor1 i,Icoor1 j) {/*{{{*/ 231 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/NearestVertex)*/ 232 233 /*Intermediaries*/ 234 BamgQuadtreeBox *pb[MaxDepth]; 235 int pi[MaxDepth]; 236 Icoor1 ii[MaxDepth]; 237 Icoor1 jj[MaxDepth]; 238 int level; 239 long n0; 240 BamgQuadtreeBox *b; 241 long h0; 242 long h = MaxISize; 243 long hb= MaxISize; 244 Icoor1 i0=0,j0=0; 196 int BamgQuadtree::BoxNumber(int i,int j,int l){/*{{{*/ 197 /* 198 * 199 * J j 200 * ^ ^ 201 * | | +--------+--------+ 202 * | | | | | 203 * 1X | | | 2 | 3 | 204 * | | | | | 205 * | | +--------+--------+ 206 * | | | | | 207 * 0X | | | 0 | 1 | 208 * | | | | | 209 * | | +--------+--------+ 210 * | +-----------------------> i 211 * | 212 * |----------------------------> I 213 * X0 X1 214 * 215 * box 0 -> I=0 J=0 BoxNumber=00 = 0 216 * box 1 -> I=1 J=0 BoxNumber=01 = 1 217 * box 2 -> I=0 J=1 BoxNumber=10 = 2 218 * box 3 -> I=1 J=1 BoxNumber=11 = 3 219 */ 220 //BoxNumber(i,j,l) returns the box number of i and j with respect to l 221 //if !j&l and !i&l -> 0 (box zero: lower left ) 222 //if !j&l and i&l -> 1 (box one: lower right) 223 //if j&l and !i&l -> 2 (box two: upper left ) 224 //if j&l and i&l -> 3 (box three:upper right) 225 return ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )); 226 }/*}}}*/ 227 bool BamgQuadtree::Intersection(int a,int b,int x,int y){/*{{{*/ 228 /*is [x y] in [a b]*/ 229 return ((y) > (a)) && ((x) <(b)); 230 }/*}}}*/ 231 BamgVertex* BamgQuadtree::NearestVertex(int xi,int yi) {/*{{{*/ 245 232 246 233 /*initial output as NULL (no vertex found)*/ 247 234 BamgVertex* nearest_v=NULL; 248 235 249 /*Project w coordinates (i,j) onto [0,MaxISize-1] x [0,MaxISize-1] -> (iplus,jplus)*/250 Icoor1 iplus( i<MaxISize ? (i<0?0:i) : MaxISize-1);251 Icoor1 jplus( j<MaxISize ? (j<0?0:j) : MaxISize-1);252 253 /*Get initial Quadtree box (largest)*/254 b = root;255 256 236 /*if the tree is empty, return NULL pointer*/ 257 if (!root->nbitems) return nearest_v; 237 if(!this->root->nbitems) return nearest_v; 238 239 /*Project coordinates (xi,yi) onto [0,MaxICoord-1] x [0,MaxICoord-1]*/ 240 int xi2 = xi; 241 int yi2 = yi; 242 if(xi<0) xi2 = 0; 243 if(xi>MaxISize) xi2 = MaxICoord; 244 if(yi<0) yi2 = 0; 245 if(yi>MaxISize) yi2 = MaxICoord; 246 247 /*Get inital box (the largest)*/ 248 BamgQuadtreeBox* b = this->root; 249 250 /*Initialize level and sizes for largest box*/ 251 int levelbin = (1L<<this->MaxDepth);// = 2^30 252 int h = this->MaxISize; 253 int hb = this->MaxISize; 254 int i0 = 0; 255 int j0 = 0; 258 256 259 257 /*else, find the smallest non-empty BamgQuadtreeBox containing the point (i,j)*/ 260 while( (n0=b->nbitems)<0){261 262 Icoor1 hb2 = hb >> 1;//size of the current box263 int k = IJ(iplus,jplus,hb2);//box number (0,1,2 or 3)264 BamgQuadtreeBox *b0 = b->b[k];//pointer toward current box265 266 /* break if NULL box orempty (Keep previous box b)*/267 if (( b0 == NULL) || (b0->nbitems ==0)) break;258 while(b->nbitems<0){ 259 260 int hb2 = hb >> 1; //size of the current box 261 int k = BoxNumber(xi2,yi2,hb2); //box number (0,1,2 or 3) 262 BamgQuadtreeBox* b0 = b->box[k]; //pointer toward current box 263 264 /* break if box is empty (Keep previous box b)*/ 265 if((b0==NULL) || (b0->nbitems==0)) break; 268 266 269 267 /*Get next Quadtree box*/ 270 b=b0; 271 i0 += I_IJ(k,hb2); // i orign of BamgQuadtreeBox (macro) 272 j0 += J_IJ(k,hb2); // j orign of BamgQuadtreeBox 273 hb = hb2; // size of the box (in Int) 268 b = b0; 269 hb = hb2; 270 this->SubBoxCoords(&i0,&j0,k,hb2); 274 271 } 275 272 … … 280 277 * - i0: x coordinate of the lower left corner 281 278 * - j0: y coordinate of the lower left corner*/ 279 int n0 = b->nbitems; 282 280 283 281 /* if the current subbox is holding vertices, we are almost done*/ 284 if 285 / /loop over the vertices of the box and find the closest vertex282 if(n0>0){ 283 /*loop over the vertices of the box and find the closest vertex*/ 286 284 for(int k=0;k<n0;k++){ 287 288 /*get integer coordinates of current vertex*/ 289 I2 i2=b->v[k]->i; 290 291 /*Compute norm with w*/ 292 h0=NORM(iplus,i2.x,jplus,i2.y); 285 int xiv = b->v[k]->i.x; 286 int yiv = b->v[k]->i.y; 287 288 289 int h0 = Norm(xi2,xiv,yi2,yiv); 293 290 294 291 /*is it smaller than previous value*/ 295 if 292 if(h0<h){ 296 293 h = h0; 297 294 nearest_v = b->v[k]; … … 306 303 307 304 /*Initialize search variables*/ 305 BamgQuadtreeBox **pb = xNew<BamgQuadtreeBox*>(this->MaxDepth); 306 int* pi = xNew<int>(this->MaxDepth); 307 int* ii = xNew<int>(this->MaxDepth); 308 int* jj = xNew<int>(this->MaxDepth); 308 309 pb[0]=b; //pointer toward the box b 309 310 pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b … … 315 316 316 317 /*Main loop*/ 317 level=0; 318 do { 319 318 int level=0; 319 do{ 320 320 /*get current box*/ 321 b = pb[level];321 b = pb[level]; 322 322 323 323 /*Loop over the items in current box (if not empty!)*/ 324 while 324 while(pi[level]){ 325 325 326 326 /*We are looping now over the items of b. k is the current index (in [0 3])*/ … … 330 330 /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/ 331 331 if (b->nbitems>0){ 332 I2 i2=b->v[k]->i; 333 h0 = NORM(iplus,i2.x,jplus,i2.y); 332 int h0 = Norm(xi2,b->v[k]->i.x,yi2,b->v[k]->i.y); 334 333 if (h0<h){ 335 334 h=h0; … … 345 344 346 345 /*if the next box exists:*/ 347 if((b=b->b [k])){346 if((b=b->box[k])){ 348 347 349 348 /*Get size (hb) and coordinates of the current sub-box lowest left corner*/ 350 349 hb>>=1; 351 Icoor1 iii = ii[level]+I_IJ(k,hb); 352 Icoor1 jjj = jj[level]+J_IJ(k,hb); 353 354 /*if the current point (iplus,jplus) is in b (modulo h), this box is good: 350 int iii = ii[level]; 351 int jjj = jj[level]; 352 this->SubBoxCoords(&iii,&jjj,k,hb); 353 354 /*if the current point (xi2,yi2) is in b (modulo h), this box is good: 355 355 * it is holding vertices that are close to w */ 356 if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){356 if(this->Intersection(iii,iii+hb,xi2-h,xi2+h) && this->Intersection(jjj,jjj+hb,yi2-h,yi2+h)){ 357 357 level++; 358 pb[level] = b;359 pi[level] = b->nbitems>0 ?(int) b->nbitems :4 ;360 ii[level] = iii;361 jj[level] = jjj;358 pb[level] = b; 359 pi[level] = b->nbitems>0 ? b->nbitems:4 ; 360 ii[level] = iii; 361 jj[level] = jjj; 362 362 } 363 363 364 / /else go backwards364 /*else go backwards*/ 365 365 else{ 366 / /shifted righ by one bit: hb=001000000 -> 01000000366 /*shifted righ by one bit: hb=001000000 -> 01000000*/ 367 367 b=b0; 368 368 hb<<=1; … … 379 379 * in case there is a vertex closest to w that has not yet been tested*/ 380 380 hb <<= 1; 381 } while (level--); 382 383 /*return nearest_v, nearest vertex*/ 381 }while(level--); 382 383 /*return nearest vertex pointer*/ 384 xDelete<BamgQuadtreeBox*>(pb); 385 xDelete<int>(pi); 386 xDelete<int>(ii); 387 xDelete<int>(jj); 384 388 return nearest_v; 385 386 389 } 387 390 /*}}}*/ 388 BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){/*{{{*/ 389 390 /*Output*/ 391 BamgQuadtreeBox* newbox=NULL; 392 393 /*Create and initialize a new box*/ 394 newbox=new BamgQuadtreeBox; 395 newbox->nbitems=0; 396 newbox->b[0]=NULL; 397 newbox->b[1]=NULL; 398 newbox->b[2]=NULL; 399 newbox->b[3]=NULL; 400 401 /*Add root to the container*/ 402 boxcontainer->AddObject(newbox); 403 404 /*Increase counter*/ 405 NbBamgQuadtreeBox++; 406 407 /*currentbox now points toward next quadtree box*/ 408 return newbox; 391 int BamgQuadtree::Norm(int xi1,int xi2,int yi1,int yi2){/*{{{*/ 392 393 int deltax = xi2 - xi1; 394 int deltay = yi2 - yi1; 395 396 if(deltax<0) deltax = -deltax; 397 if(deltay<0) deltay = -deltay; 398 399 if(deltax> deltay){ 400 return deltax; 401 } 402 else{ 403 return deltay; 404 } 409 405 }/*}}}*/ 410 BamgVertex* BamgQuadtree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){/*{{{*/ 411 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/ToClose)*/ 412 413 const Icoor1 i=v.i.x; 414 const Icoor1 j=v.i.y; 415 const R2 X(v.r); 416 const Metric Mx(v.m); 417 418 BamgQuadtreeBox * pb[ MaxDepth ]; 419 int pi[ MaxDepth ]; 420 Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; 406 void BamgQuadtree::SubBoxCoords(int* pi,int*pj,int boxnumber,int length){/*{{{*/ 407 /* 408 * j (first bit) 409 * ^ 410 * | +--------+--------+ 411 * | | | | 412 * 1X | | 2 | 3 | 413 * | | | | 414 * | +--------+--------+ 415 * | | | | 416 * 0X | | 0 | 1 | 417 * | | | | 418 * | +--------+--------+ 419 * +-----------------------> i (second bit) 420 * X0 X1 421 */ 422 423 /*Add sub-box coordinate to i and j*/ 424 //_assert_(!(*pi & length)); 425 //_assert_(!(*pj & length)); 426 427 /*length if first bit of boxnumber is 1, elengthse 0*/ 428 *pi += ((boxnumber & 1) ? length:0); 429 /*length if second bit of boxnumber is 1, elengthse 0*/ 430 *pj += ((boxnumber & 2) ? length:0); 431 432 }/*}}}*/ 433 BamgVertex* BamgQuadtree::TooClose(BamgVertex* v,double threshold,int hx,int hy){/*{{{*/ 434 435 const int i=v->i.x; 436 const int j=v->i.y; 437 const double Xx = v->r.x; 438 const double Xy = v->r.y; 439 Metric* Mx = new Metric(v->m); 440 441 BamgQuadtreeBox *pb[MaxDepth]; 442 int pi[MaxDepth]; 443 int ii[MaxDepth], jj[MaxDepth]; 421 444 int l=0; // level 422 BamgQuadtreeBox * b; 423 Icoor1 hb = MaxISize; 424 Icoor1 i0=0,j0=0; 425 426 // BamgVertex *vn=0; 427 428 if (!root->nbitems) 429 return 0; // empty tree 430 431 // general case ----- 445 BamgQuadtreeBox* b; 446 int hb = MaxISize; 447 int i0=0,j0=0; 448 449 // BamgVertex *vn=0; 450 451 if(!root->nbitems) return 0; // empty tree 452 453 // general case 432 454 pb[0]=root; 433 pi[0]=root->nbitems>0 ?(int) root->nbitems : 4;455 pi[0]=root->nbitems>0 ?(int)root->nbitems:4; 434 456 ii[0]=i0; 435 457 jj[0]=j0; 436 do {458 do{ 437 459 b= pb[l]; 438 while 460 while(pi[l]--){ 439 461 int k = pi[l]; 440 462 441 if (b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty 442 I2 i2 = b->v[k]->i; 443 if ( ABS(i-i2.x) <hx && ABS(j-i2.y) <hy ) 444 { 445 R2 XY(X,b->v[k]->r); 446 if(LengthInterpole(Mx(XY), b->v[k]->m(XY)) < seuil){ 463 if(b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty 464 int i2x = b->v[k]->i.x; 465 int i2y = b->v[k]->i.y; 466 if (ABS(i-i2x)<hx && ABS(j-i2y) <hy ){ 467 double XYx = b->v[k]->r.x - Xx; 468 double XYy = b->v[k]->r.y - Xy; 469 if(LengthInterpole(Mx->Length(XYx,XYy),b->v[k]->m.Length(XYx,XYy)) < threshold){ 470 delete Mx; 447 471 return b->v[k]; 448 472 } 449 473 } 450 474 } 451 475 else{ // Pointer BamgQuadtreeBox 452 476 BamgQuadtreeBox *b0=b; 453 if ((b=b->b [k])){477 if ((b=b->box[k])){ 454 478 hb >>=1 ; // div by 2 455 long iii = ii[l]+I_IJ(k,hb); 456 long jjj = jj[l]+J_IJ(k,hb); 457 458 if (INTER_SEG(iii,iii+hb,i-hx,i+hx) && INTER_SEG(jjj,jjj+hb,j-hy,j+hy)){ 479 int iii = ii[l]; 480 int jjj = jj[l]; 481 this->SubBoxCoords(&iii,&jjj,k,hb); 482 483 if(this->Intersection(iii,iii+hb,i-hx,i+hx) && this->Intersection(jjj,jjj+hb,j-hy,j+hy)){ 459 484 pb[++l]= b; 460 485 pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; … … 474 499 } 475 500 hb <<= 1; // mul by 2 476 } while (l--); 477 501 }while(l--); 502 503 delete Mx; 478 504 return 0; 479 505 } 480 506 /*}}}*/ 507 508 BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){/*{{{*/ 509 510 /*Output*/ 511 BamgQuadtreeBox* newbox=NULL; 512 513 /*Create and initialize a new box*/ 514 newbox=new BamgQuadtreeBox; 515 newbox->nbitems=0; 516 newbox->box[0]=NULL; 517 newbox->box[1]=NULL; 518 newbox->box[2]=NULL; 519 newbox->box[3]=NULL; 520 521 /*Add root to the container*/ 522 boxcontainer->AddObject(newbox); 523 this->NbQuadtreeBox++; 524 525 /*currentbox now points toward next quadtree box*/ 526 return newbox; 527 }/*}}}*/ 481 528 }
Note:
See TracChangeset
for help on using the changeset viewer.
