Changeset 12210 for issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp

Timestamp:

05/04/12 12:26:38 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added Observations dataset

File:

• issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp

@@ -47 +47 @@
  * we just need to look at the bits from the left to the right (See ::Add)
  }}}1*/
+/*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
+ */
+//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)
+#define IJ(i,j,l)  ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 ))
+/*}}}*/
 
         /*Constructors/Destructors*/
@@ -54 +85 @@
                 /*Number of boxes and vertices*/
                 NbQuadtreeBox=0;
-                NbPoints=0;
+                NbObs=0;
 
                 /*Create container*/
@@ -74 +105 @@
 
         /*Methods*/
+/*FUNCTION Quadtree::Add{{{1*/
+void  Quadtree::Add(Observation* observation){
+
+        /*Intermediaries*/
+        const int     MaxDeep  = 30;
+        int           k,ij;
+        long          i,j,level;
+        QuadtreeBox **pb = NULL;
+        QuadtreeBox  *b  = NULL;
+
+        /*Get integer coodinates*/
+        i = observation->xi;
+        j = observation->yi;
+
+        /*Initialize level*/
+        level=(1L<<MaxDeep);// = 2^30
+
+        /*Get inital box (the largest)*/
+        pb=&root;
+
+        /*Find the smallest box where the observation is located*/
+        while((b=*pb) && (b->nbitems<0)){ 
+
+                /*shift b->nbitems by -1 as a counter*/
+                b->nbitems--;
+
+                /*Go down one level = 00100 -> 00010*/
+                level>>=1;
+
+                /*Get next subbox according to the bit value (level)*/
+                pb = &b->box[IJ(i,j,level)];
+        }
+        _assert_(level>0);
+
+        /*Now, try to add the vertex, if the box is full (nbitems=4), we have to divide it in 4 new boxes*/
+        while((b=*pb) && (b->nbitems==4)){
+
+                /*Copy the 4 observation in the current Quadtreebox*/
+                Observation* obs[4];
+                obs[0] = b->obs[0];
+                obs[1] = b->obs[1];
+                obs[2] = b->obs[2];
+                obs[3] = b->obs[3];
+
+                /*set nbitems as negative (now holding boxes instead of observations)*/
+                b->nbitems = -b->nbitems;
+
+                /*Initialize the 4 pointers toward the 4 subboxes*/
+                b->box[0]=NULL;
+                b->box[1]=NULL;
+                b->box[2]=NULL;
+                b->box[3]=NULL;
+
+                /*level = 00100 -> 00010*/
+                level>>=1;
+
+                /*Put the four observations in the new boxes*/
+                for (k=0;k<4;k++){
+
+                        /*Get box for observation number k*/
+                        ij=IJ(obs[k]->xi,obs[k]->yi,level);
+                        QuadtreeBox *bb = b->box[ij];
+                        if(!bb){
+                                b->box[ij]=NewQuadtreeBox();
+                                bb=b->box[ij];
+                        }
+
+                        /*Copy current observations*/
+                        bb->obs[bb->nbitems++] = obs[k];
+                }
+
+                /*Get the subbox where the current observation is located*/
+                ij=IJ(i,j,level);
+                pb=&b->box[ij];
+        }
+
+        /*alloc the QuadtreeBox if necessary and add current observation*/
+        b=*pb;
+        if(!b) b=NewQuadtreeBox();
+        b->obs[b->nbitems++]=observation;
+        NbObs++;
+
+}/*}}}*/
         /*FUNCTION Quadtree::NewQuadtreeBox {{{1*/
         Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){