Ice Sheet System Model  4.18
Code documentation
Data Structures | Public Member Functions | Data Fields | Private Attributes
bamg::BamgQuadtree Class Reference

#include <BamgQuadtree.h>

Data Structures

class  BamgQuadtreeBox
 

Public Member Functions

 BamgQuadtree ()
 
 BamgQuadtree (Mesh *t, long nbv=-1)
 
 ~BamgQuadtree ()
 
void Add (BamgVertex &w)
 
int BoxNumber (int i, int j, int l)
 
bool Intersection (int a, int b, int x, int y)
 
BamgVertex * NearestVertex (int i, int j)
 
BamgQuadtreeBox * NewBamgQuadtreeBox (void)
 
int Norm (int xi1, int yi1, int xi2, int yi2)
 
void SubBoxCoords (int *pi, int *pj, int boxcoord, int length)
 
BamgVertex * TooClose (BamgVertex *, double, int, int)
 

Data Fields

BamgQuadtreeBox * root
 
long NbQuadtreeBox
 
long NbVertices
 

Private Attributes

std::vector< BamgQuadtreeBox * > boxcontainer
 

Detailed Description

Definition at line 10 of file BamgQuadtree.h.

Constructor & Destructor Documentation

◆ BamgQuadtree() [1/2]

bamg::BamgQuadtree::BamgQuadtree ( )

Definition at line 62 of file BamgQuadtree.cpp.

62  {/*{{{*/
63 
64  /*Initialize fields*/
65  this->NbQuadtreeBox = 0;
66  this->NbVertices = 0;
67 
68  /*Create Root, pointer toward the main box*/
69  this->root=NewBamgQuadtreeBox();
70 
71  } /*}}}*/

◆ BamgQuadtree() [2/2]

bamg::BamgQuadtree::BamgQuadtree ( Mesh *  t,
long  nbv = -1 
)

Definition at line 72 of file BamgQuadtree.cpp.

72  { /*{{{*/
73 
74  /*Initialize fields*/
75  this->NbQuadtreeBox = 0;
76  this->NbVertices = 0;
77 
78  /*Create Root, pointer toward the main box*/
79  this->root=NewBamgQuadtreeBox();
80 
81  /*Check Sizes*/
82  _assert_(MAXISIZE>MAXICOORD);
83 
84  /*Add all vertices of the mesh*/
85  if(nbv==-1) nbv=t->nbv;
86  for(int i=0;i<nbv;i++){
87  this->Add(t->vertices[i]);
88  }
89 
90  }

◆ ~BamgQuadtree()

bamg::BamgQuadtree::~BamgQuadtree ( )

Definition at line 92 of file BamgQuadtree.cpp.

92  {/*{{{*/
93 
94  vector<BamgQuadtreeBox*>::reverse_iterator object;
95  for(object=boxcontainer.rbegin() ; object <boxcontainer.rend(); object++ ){
96  delete (*object);
97  }
98  boxcontainer.clear();
99  root=NULL;
100  }

Member Function Documentation

◆ Add()

void bamg::BamgQuadtree::Add ( BamgVertex &  w)

Definition at line 104 of file BamgQuadtree.cpp.

104  {/*{{{*/
105  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.cpp/Add)*/
106  BamgQuadtreeBox** pb=NULL;
107  BamgQuadtreeBox* b=NULL;
108 
109  /*Get integer coodinate of current point w*/
110  long i=w.i.x, j=w.i.y;
111 
112  /*Initialize level*/
113  long level=MAXISIZE;
114 
115  /*Get inital box (the largest)*/
116  pb = &root;
117 
118  /*Find the smallest box where w is located*/
119  while((b=*pb) && (b->nbitems<0)){
120 
121  //shift b->nbitems by -1
122  b->nbitems--;
123 
124  //shifted righ by one bit: level=00000010 -> 00000001
125  level >>= 1;
126 
127  //Get next subbox according to the bit value (level)
128  pb = &b->box[BoxNumber(i,j,level)];
129  }
130 
131  /*OK, we have found b, a Subbox holding vertices (might be full)
132  check that the vertex is not already in the box*/
133  if (b){
134  if (b->nbitems > 3 && b->v[3] == &w) return;
135  if (b->nbitems > 2 && b->v[2] == &w) return;
136  if (b->nbitems > 1 && b->v[1] == &w) return;
137  if (b->nbitems > 0 && b->v[0] == &w) return;
138  }
139 
140  /*check that l is not 0 (this should not happen as MAXDEPTH = 30)*/
141  _assert_(level>0);
142 
143  /*Now, try to add the vertex, if the subbox is full (nbitems=4), we have to divide it
144  in 4 new subboxes*/
145  while ((b= *pb) && (b->nbitems == 4)){ // the BamgQuadtreeBox is full
146 
147  /*Copy the 4 vertices in the current BamgQuadtreebox*/
148  BamgVertex* v4[4];
149  v4[0]= b->v[0];
150  v4[1]= b->v[1];
151  v4[2]= b->v[2];
152  v4[3]= b->v[3];
153 
154  /*set nbitems as negative
155  * (box full -> holds 4 pointers toward subboxes and not 4 vertices)*/
156  b->nbitems = -b->nbitems;
157 
158  /*Initialize the 4 pointers toward the 4 subboxes*/
159  b->box[0]=b->box[1]=b->box[2]=b->box[3]=NULL;
160 
161  /*level = 0010000 -> 0001000*/
162  level >>= 1;
163 
164  /*Put the four vertices in the new boxes*/
165  for (int k=0;k<4;k++){
166 
167  int ij;
168  /*bb is the new "sub"box of b where v4[k] is located*/
169  BamgQuadtreeBox *bb = b->box[ij=BoxNumber(v4[k]->i.x,v4[k]->i.y,level)];
170 
171  // alloc the BamgQuadtreeBox
172  if (!bb) bb=b->box[ij]=NewBamgQuadtreeBox();
173 
174  /*Copy the current vertex*/
175  bb->v[bb->nbitems++] = v4[k];
176  }
177 
178  /*Get the subbox where w (i,j) is located*/
179  pb = &b->box[BoxNumber(i,j,level)];
180  }
181 
182  /*alloc the BamgQuadtreeBox if necessary*/
183  if (!(b=*pb)) b=*pb= NewBamgQuadtreeBox();
184 
185  /*Add w*/
186  b->v[b->nbitems++]=&w;
187 
188  //Increase NbVertices by one (we have one new vertex)
189  NbVertices++;
190  }

◆ BoxNumber()

int bamg::BamgQuadtree::BoxNumber ( int  i,
int  j,
int  l 
)

Definition at line 192 of file BamgQuadtree.cpp.

192  {/*{{{*/
193  /*
194  *
195  * J j
196  * ^ ^
197  * | | +--------+--------+
198  * | | | | |
199  * 1X | | | 2 | 3 |
200  * | | | | |
201  * | | +--------+--------+
202  * | | | | |
203  * 0X | | | 0 | 1 |
204  * | | | | |
205  * | | +--------+--------+
206  * | +-----------------------> i
207  * |
208  * |----------------------------> I
209  * X0 X1
210  *
211  * box 0 -> I=0 J=0 BoxNumber=00 = 0
212  * box 1 -> I=1 J=0 BoxNumber=01 = 1
213  * box 2 -> I=0 J=1 BoxNumber=10 = 2
214  * box 3 -> I=1 J=1 BoxNumber=11 = 3
215  */
216  //BoxNumber(i,j,l) returns the box number of i and j with respect to l
217  //if !j&l and !i&l -> 0 (box zero: lower left )
218  //if !j&l and i&l -> 1 (box one: lower right)
219  //if j&l and !i&l -> 2 (box two: upper left )
220  //if j&l and i&l -> 3 (box three:upper right)
221  return ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 ));
222  }/*}}}*/

◆ Intersection()

bool bamg::BamgQuadtree::Intersection ( int  a,
int  b,
int  x,
int  y 
)

Definition at line 223 of file BamgQuadtree.cpp.

223  {/*{{{*/
224  /*is [x y] in [a b]*/
225  return ((y) > (a)) && ((x) <(b));
226  }/*}}}*/

◆ NearestVertex()

BamgVertex * bamg::BamgQuadtree::NearestVertex ( int  i,
int  j 
)

Definition at line 227 of file BamgQuadtree.cpp.

227  {/*{{{*/
228 
229  /*initial output as NULL (no vertex found)*/
230  BamgVertex* nearest_v=NULL;
231 
232  /*if the tree is empty, return NULL pointer*/
233  if(!this->root->nbitems) return nearest_v;
234 
235  /*Project coordinates (xi,yi) onto [0,MAXICOORD-1] x [0,MAXICOORD-1]*/
236  int xi2 = xi;
237  int yi2 = yi;
238  if(xi<0) xi2 = 0;
239  if(xi>MAXISIZE) xi2 = MAXICOORD;
240  if(yi<0) yi2 = 0;
241  if(yi>MAXISIZE) yi2 = MAXICOORD;
242 
243  /*Get inital box (the largest)*/
244  BamgQuadtreeBox* b = this->root;
245 
246  /*Initialize level and sizes for largest box*/
247  int levelbin = (1L<<MAXDEPTH);// = 2^30
248  int h = MAXISIZE;
249  int hb = MAXISIZE;
250  int i0 = 0;
251  int j0 = 0;
252 
253  /*else, find the smallest non-empty BamgQuadtreeBox containing the point (i,j)*/
254  while(b->nbitems<0){
255 
256  int hb2 = hb >> 1; //size of the current box
257  int k = BoxNumber(xi2,yi2,hb2); //box number (0,1,2 or 3)
258  BamgQuadtreeBox* b0 = b->box[k]; //pointer toward current box
259 
260  /* break if box is empty (Keep previous box b)*/
261  if((b0==NULL) || (b0->nbitems==0)) break;
262 
263  /*Get next Quadtree box*/
264  b = b0;
265  hb = hb2;
266  this->SubBoxCoords(&i0,&j0,k,hb2);
267  }
268 
269  /*The box b, is the smallest box containing the point (i,j) and
270  * has the following properties:
271  * - n0: number of items (>0 if vertices, else boxes)
272  * - hb: box size (int)
273  * - i0: x coordinate of the lower left corner
274  * - j0: y coordinate of the lower left corner*/
275  int n0 = b->nbitems;
276 
277  /* if the current subbox is holding vertices, we are almost done*/
278  if(n0>0){
279  /*loop over the vertices of the box and find the closest vertex*/
280  for(int k=0;k<n0;k++){
281  int xiv = b->v[k]->i.x;
282  int yiv = b->v[k]->i.y;
283 
284  int h0 = Norm(xi2,xiv,yi2,yiv);
285 
286  /*is it smaller than previous value*/
287  if(h0<h){
288  h = h0;
289  nearest_v = b->v[k];
290  }
291  }
292  /*return closest vertex*/
293  return nearest_v;
294  }
295 
296  /* general case: the current box is empty, we have to go backwards
297  and find the closest not-empty box and find the closest vertex*/
298 
299  /*Initialize search variables*/
300  BamgQuadtreeBox* pb[MAXDEPTH];
301  int pi[MAXDEPTH];
302  int ii[MAXDEPTH];
303  int jj[MAXDEPTH];
304  pb[0]=b; //pointer toward the box b
305  pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b
306  ii[0]=i0; //i coordinate of the box lowest left corner
307  jj[0]=j0; //j coordinate of the box lowest left corner
308 
309  /*initialize h: smallest box size, containing a vertex close to w*/
310  h=hb;
311 
312  /*Main loop*/
313  int level=0;
314  do{
315  /*get current box*/
316  b = pb[level];
317 
318  /*Loop over the items in current box (if not empty!)*/
319  while(pi[level]){
320 
321  /*We are looping now over the items of b. k is the current index (in [0 3])*/
322  pi[level]--;
323  int k=pi[level];
324 
325  /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/
326  if (b->nbitems>0){
327  int h0 = Norm(xi2,b->v[k]->i.x,yi2,b->v[k]->i.y);
328  if (h0<h){
329  h=h0;
330  nearest_v=b->v[k];
331  }
332  }
333  /*else: current box b is pointing toward 4 boxes
334  * test sub-box k and go deeper into the tree if it is non empty
335  * and contains the point w modulo a size h that is either the size of the smallest
336  * non empty box containing w, or the closest point to w (so far) */
337  else{
338  BamgQuadtreeBox* b0=b;
339 
340  /*if the next box exists:*/
341  if((b=b->box[k])){
342 
343  /*Get size (hb) and coordinates of the current sub-box lowest left corner*/
344  hb>>=1;
345  int iii = ii[level];
346  int jjj = jj[level];
347  this->SubBoxCoords(&iii,&jjj,k,hb);
348 
349  /*if the current point (xi2,yi2) is in b (modulo h), this box is good:
350  * it is holding vertices that are close to w */
351  if(this->Intersection(iii,iii+hb,xi2-h,xi2+h) && this->Intersection(jjj,jjj+hb,yi2-h,yi2+h)){
352  level++;
353  pb[level] = b;
354  pi[level] = b->nbitems>0 ? b->nbitems:4 ;
355  ii[level] = iii;
356  jj[level] = jjj;
357  }
358 
359  /*else go backwards*/
360  else{
361  /*shifted righ by one bit: hb=001000000 -> 01000000*/
362  b=b0;
363  hb<<=1;
364  }
365  }
366  else{
367  /*Current box is NULL, go to next subbox of b (k=k-1)*/
368  b=b0;
369  }
370  }
371  }
372 
373  /*We have found a vertex, now, let's try the other boxes of the previous level
374  * in case there is a vertex closest to w that has not yet been tested*/
375  hb <<= 1;
376  }while(level--);
377 
378  /*return nearest vertex pointer*/
379  return nearest_v;
380  }

◆ NewBamgQuadtreeBox()

BamgQuadtree::BamgQuadtreeBox * bamg::BamgQuadtree::NewBamgQuadtreeBox ( void  )

Definition at line 499 of file BamgQuadtree.cpp.

499  {/*{{{*/
500 
501  /*Output*/
502  BamgQuadtreeBox* newbox=NULL;
503 
504  /*Create and initialize a new box*/
505  newbox=new BamgQuadtreeBox;
506  newbox->nbitems=0;
507  newbox->box[0]=NULL;
508  newbox->box[1]=NULL;
509  newbox->box[2]=NULL;
510  newbox->box[3]=NULL;
511 
512  /*Add root to the container*/
513  boxcontainer.push_back(newbox);
514  this->NbQuadtreeBox++;
515 
516  /*currentbox now points toward next quadtree box*/
517  return newbox;
518  }/*}}}*/

◆ Norm()

int bamg::BamgQuadtree::Norm ( int  xi1,
int  yi1,
int  xi2,
int  yi2 
)

Definition at line 382 of file BamgQuadtree.cpp.

382  {/*{{{*/
383 
384  int deltax = xi2 - xi1;
385  int deltay = yi2 - yi1;
386 
387  if(deltax<0) deltax = -deltax;
388  if(deltay<0) deltay = -deltay;
389 
390  if(deltax> deltay){
391  return deltax;
392  }
393  else{
394  return deltay;
395  }
396  }/*}}}*/

◆ SubBoxCoords()

void bamg::BamgQuadtree::SubBoxCoords ( int *  pi,
int *  pj,
int  boxcoord,
int  length 
)

Definition at line 397 of file BamgQuadtree.cpp.

397  {/*{{{*/
398  /*
399  * j (first bit)
400  * ^
401  * | +--------+--------+
402  * | | | |
403  * 1X | | 2 | 3 |
404  * | | | |
405  * | +--------+--------+
406  * | | | |
407  * 0X | | 0 | 1 |
408  * | | | |
409  * | +--------+--------+
410  * +-----------------------> i (second bit)
411  * X0 X1
412  */
413 
414  /*Add sub-box coordinate to i and j*/
415  //_assert_(!(*pi & length));
416  //_assert_(!(*pj & length));
417 
418  /*length if first bit of boxnumber is 1, elengthse 0*/
419  *pi += ((boxnumber & 1) ? length:0);
420  /*length if second bit of boxnumber is 1, elengthse 0*/
421  *pj += ((boxnumber & 2) ? length:0);
422 
423  }/*}}}*/

◆ TooClose()

BamgVertex * bamg::BamgQuadtree::TooClose ( BamgVertex *  v,
double  threshold,
int  hx,
int  hy 
)

Definition at line 424 of file BamgQuadtree.cpp.

424  {/*{{{*/
425 
426  const int i=v->i.x;
427  const int j=v->i.y;
428  const double Xx = v->r.x;
429  const double Xy = v->r.y;
430  Metric* Mx = new Metric(v->m);
431 
432  BamgQuadtreeBox *pb[MAXDEPTH];
433  int pi[MAXDEPTH];
434  int ii[MAXDEPTH], jj[MAXDEPTH];
435  int l=0; // level
436  BamgQuadtreeBox* b;
437  int hb = MAXISIZE;
438  int i0=0,j0=0;
439 
440  // BamgVertex *vn=0;
441 
442  if(!root->nbitems) return 0; // empty tree
443 
444  // general case
445  pb[0]=root;
446  pi[0]=root->nbitems>0 ?(int)root->nbitems:4;
447  ii[0]=i0;
448  jj[0]=j0;
449  do{
450  b= pb[l];
451  while(pi[l]--){
452  int k = pi[l];
453 
454  if(b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty
455  int i2x = b->v[k]->i.x;
456  int i2y = b->v[k]->i.y;
457  if (ABS(i-i2x)<hx && ABS(j-i2y) <hy ){
458  double XYx = b->v[k]->r.x - Xx;
459  double XYy = b->v[k]->r.y - Xy;
460  if(LengthInterpole(Mx->Length(XYx,XYy),b->v[k]->m.Length(XYx,XYy)) < threshold){
461  delete Mx;
462  return b->v[k];
463  }
464  }
465  }
466  else{ // Pointer BamgQuadtreeBox
467  BamgQuadtreeBox *b0=b;
468  if ((b=b->box[k])){
469  hb >>=1 ; // div by 2
470  int iii = ii[l];
471  int jjj = jj[l];
472  this->SubBoxCoords(&iii,&jjj,k,hb);
473 
474  if(this->Intersection(iii,iii+hb,i-hx,i+hx) && this->Intersection(jjj,jjj+hb,j-hy,j+hy)){
475  pb[++l]= b;
476  pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ;
477  ii[l]= iii;
478  jj[l]= jjj;
479 
480  }
481  else{
482  b=b0;
483  hb <<=1 ;
484  }
485  }
486  else{
487  b=b0;
488  }
489  }
490  }
491  hb <<= 1; // mul by 2
492  }while(l--);
493 
494  delete Mx;
495  return 0;
496  }

Field Documentation

◆ boxcontainer

std::vector<BamgQuadtreeBox*> bamg::BamgQuadtree::boxcontainer
private

Definition at line 27 of file BamgQuadtree.h.

◆ root

BamgQuadtreeBox* bamg::BamgQuadtree::root

Definition at line 32 of file BamgQuadtree.h.

◆ NbQuadtreeBox

long bamg::BamgQuadtree::NbQuadtreeBox

Definition at line 33 of file BamgQuadtree.h.

◆ NbVertices

long bamg::BamgQuadtree::NbVertices

Definition at line 34 of file BamgQuadtree.h.

The documentation for this class was generated from the following files:
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
bamg::BamgQuadtree::BamgQuadtreeBox::nbitems
int nbitems
Definition: BamgQuadtree.h:21
MAXISIZE
#define MAXISIZE
Definition: BamgQuadtree.cpp:10
bamg::BamgQuadtree::NbVertices
long NbVertices
Definition: BamgQuadtree.h:34
bamg::BamgQuadtree::BoxNumber
int BoxNumber(int i, int j, int l)
Definition: BamgQuadtree.cpp:192
bamg::BamgQuadtree::root
BamgQuadtreeBox * root
Definition: BamgQuadtree.h:32
bamg::BamgQuadtree::SubBoxCoords
void SubBoxCoords(int *pi, int *pj, int boxcoord, int length)
Definition: BamgQuadtree.cpp:397
ABS
#define ABS(i)
Definition: BamgQuadtree.cpp:8
bamg::BamgQuadtree::boxcontainer
std::vector< BamgQuadtreeBox * > boxcontainer
Definition: BamgQuadtree.h:27
bamg::LengthInterpole
double LengthInterpole(const Metric &Ma, const Metric &Mb, R2 AB)
Definition: Metric.cpp:158
bamg::BamgQuadtree::Add
void Add(BamgVertex &w)
Definition: BamgQuadtree.cpp:104
MAXICOORD
#define MAXICOORD
Definition: BamgQuadtree.cpp:11
bamg::BamgQuadtree::Intersection
bool Intersection(int a, int b, int x, int y)
Definition: BamgQuadtree.cpp:223
bamg::BamgQuadtree::Norm
int Norm(int xi1, int yi1, int xi2, int yi2)
Definition: BamgQuadtree.cpp:382
bamg::BamgQuadtree::NewBamgQuadtreeBox
BamgQuadtreeBox * NewBamgQuadtreeBox(void)
Definition: BamgQuadtree.cpp:499
bamg::BamgQuadtree::NbQuadtreeBox
long NbQuadtreeBox
Definition: BamgQuadtree.h:33
MAXDEPTH
#define MAXDEPTH
Definition: BamgQuadtree.cpp:9
Generated on Thu Jul 2 2020 08:09:31 for Ice Sheet System Model by   doxygen 1.8.19