Changeset 12213

Timestamp:

05/04/12 16:53:41 (13 years ago)

Author:

Mathieu Morlighem

Message:

Fixed quadtree

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (4 diffs)
• Container/Observations.h (modified) (2 diffs)
• Container/Vertices.h (modified) (1 diff)
• modules/Krigingx/Krigingx.cpp (modified) (2 diffs)
• objects/Kriging/Observation.cpp (modified) (2 diffs)
• objects/Kriging/Quadtree.cpp (modified) (5 diffs)
• objects/Kriging/Quadtree.h (modified) (1 diff)
• shared/Elements/elements.h (modified) (1 diff)

issm/trunk-jpl/src/c/Container/Observations.cpp

@@ -28 +28 @@
 /*FUNCTION Observations::Observations(){{{*/
 Observations::Observations(){
+        this->quadtree = NULL;
         return;
 }
@@ -34 +35 @@
 Observations::Observations(double* observations_list,double* x,double* y,int n){
 
+        const int MaxICoor = 1073741823; // 2^30-1 =111...111 (29 times)
+
+        /*Intermediaries*/
+        int          i;
+        int          xi,yi;
+        double       offset;
+        Observation *observation = NULL;
+
+        /*Initialize Quadtree*/
+        this->quadtree = new Quadtree();
+
+        /*Get extrema*/
+        this->xmin=x[0]; this->ymin=y[0];
+        this->xmax=x[0]; this->ymax=y[0];
+        for(i=1;i<n;i++){
+                xmin=min(xmin,x[i]); ymin=min(ymin,y[i]);
+                xmax=max(xmax,x[i]); ymax=max(ymax,y[i]);
+        }
+        offset=0.05*(xmax-xmin); xmin-=offset; xmax+=offset;
+        offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
+
+        /*coeffIcoor is the coefficient used for integer coordinates:
+         *                (x-xmin)
+         * xi = (2^30 -1) ------------ 
+         *                   D
+         * coefficient = (2^30 -1)/D
+         */
+        this->coefficient = double(MaxICoor)/max(xmax-xmin,ymax-ymin); _assert_(coefficient>=0);
+
         /*Add observations one by one*/
-        for(int i=0;i<n;i++){
-                this->AddObject(new Observation(x[i],y[i],0,0,observations_list[i]));
+        for(i=0;i<n;i++){
+                xi=int(coefficient*(x[i]-xmin));
+                yi=int(coefficient*(y[i]-ymin));
+                observation = new Observation(x[i],y[i],xi,yi,observations_list[i]);
+                this->quadtree->Add(observation);
+                this->AddObject(observation);
         }
-
 }
 /*}}}*/
 /*FUNCTION Observations::~Observations(){{{*/
 Observations::~Observations(){
+        delete quadtree;
         return;
 }
@@ -49 +83 @@
 /*Methods*/
 /*FUNCTION Observations::ObservationList{{{*/
-void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp){/*{{{*/
+void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp){
 
         /*Output and Intermediaries*/
@@ -78 +112 @@
         *pnobs=nobs;
 }/*}}}*/
+/*FUNCTION Observations::QuadtreeColoring{{{*/
+void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){
+
+        /*Convert to integer coordinates*/
+        int *xi = (int*)xmalloc(n*sizeof(int));
+        int *yi = (int*)xmalloc(n*sizeof(int));
+        for(int i=0;i<n;i++){
+                xi[i]=int(coefficient*(x[i]-xmin));
+                yi[i]=int(coefficient*(y[i]-ymin));
+        }
+
+        /*Call quadtree method*/
+        this->quadtree->QuadtreeColoring(A,xi,yi,n);
+
+        /*clean-up*/
+        xfree((void**)&xi);
+        xfree((void**)&yi);
+
+}/*}}}*/

issm/trunk-jpl/src/c/Container/Observations.h

@@ -7 +7 @@
 
 class Obsevration;
+class Quadtree;
 
 class Observations: public DataSet{
+
+        private:
+                Quadtree* quadtree;
+                double    coefficient; //For integer coordinate conversion
+                double    xmin;
+                double    ymin;
+                double    xmax;
+                double    ymax;
 
         public:
@@ -19 +28 @@
                 /*Methods*/
                 void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp);
+                void QuadtreeColoring(double* A,double *x,double *y,int n);
 
 };

issm/trunk-jpl/src/c/Container/Vertices.h

@@ -32 +32 @@
 };
 
-
-
 #endif //ifndef _VERTICES_H_
-

issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

@@ -44 +44 @@
         /*Allocation output*/
         predictions =(double*)xmalloc(n_interp*sizeof(double));
+
+        for(i=0;i<n_interp;i++) predictions[i]=0;
+        observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
+        *ppredictions=predictions;
+        return 1;
 
         /*Loop over all interpolations*/
@@ -101 +106 @@
         delete observations;
         *ppredictions=predictions;
+        return 1;
 }
 

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

@@ -3 +3 @@
  */
 
+#include <stdlib.h>
 #include "../objects.h"
 
@@ -31 +32 @@
 /*FUNCTION Observation::Echo {{{1*/
 void Observation::Echo(void){
+
+        int  bit;
+
         printf("Observation\n");
         printf("   x     : %g\n",this->x);
         printf("   y     : %g\n",this->y);
-        printf("   xi    : %i\n",this->xi);
-        printf("   yi    : %i\n",this->yi);
+        printf("   xi    : "); printbinary(this->xi); printf("\n");
+        printf("   yi    : "); printbinary(this->yi); printf("\n");
         printf("   value : %g\n",this->value);
 }

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

@@ -110 +110 @@
         /*Intermediaries*/
         const int     MaxDeep  = 30;
-        int           k,ij;
-        long          i,j,level;
+        int          xi,yi,ij,level;
         QuadtreeBox **pb = NULL;
         QuadtreeBox  *b  = NULL;
+        QuadtreeBox  *bb = NULL;
+        Observation  *obs[4];
 
         /*Get integer coodinates*/
-        i = observation->xi;
-        j = observation->yi;
+        xi = observation->xi;
+        yi = observation->yi;
 
         /*Initialize level*/
@@ -128 +129 @@
         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)];
+                pb = &b->box[IJ(xi,yi,level)];
         }
         _assert_(level>0);
@@ -143 +141 @@
 
                 /*Copy the 4 observation in the current Quadtreebox*/
-                Observation* obs[4];
                 obs[0] = b->obs[0];
                 obs[1] = b->obs[1];
@@ -162 +159 @@
 
                 /*Put the four observations in the new boxes*/
-                for (k=0;k<4;k++){
+                for (int 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];
+                        ij = IJ(obs[k]->xi,obs[k]->yi,level);
+                        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);
+                ij=IJ(xi,yi,level);
                 pb=&b->box[ij];
         }
@@ -183 +178 @@
         /*alloc the QuadtreeBox if necessary and add current observation*/
         b=*pb;
-        if(!b) b=NewQuadtreeBox();
+        if(!b){
+                b=*pb=NewQuadtreeBox();
+        }
         b->obs[b->nbitems++]=observation;
         NbObs++;
 
 }/*}}}*/
-        /*FUNCTION Quadtree::NewQuadtreeBox {{{1*/
-        Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){
-
-                /*Output*/
-                QuadtreeBox* newbox=NULL;
-
-                /*Create and initialize a new box*/
-                newbox=new QuadtreeBox();
-                newbox->nbitems=0;
-                newbox->box[0]=NULL;
-                newbox->box[1]=NULL;
-                newbox->box[2]=NULL;
-                newbox->box[3]=NULL;
-
-                /*Add root to the container*/
-                boxcontainer->AddObject(newbox);
-
-                /*Increase counter*/
-                NbQuadtreeBox++;
-
-                /*currentbox now points toward next quadtree box*/
-                return newbox;
-        }/*}}}*/
+/*FUNCTION Quadtree::Echo{{{1*/
+void  Quadtree::Echo(void){
+
+        printf("Quadtree:\n");
+        printf("   NbQuadtreeBox = %i\n",NbQuadtreeBox);
+        printf("   NbObs         = %i\n",NbObs);
+        printf("   root          = %p\n",root);
+
+}/*}}}*/
+/*FUNCTION Quadtree::NewQuadtreeBox {{{1*/
+Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){
+
+        /*Output*/
+        QuadtreeBox* newbox=NULL;
+
+        /*Create and initialize a new box*/
+        newbox=new QuadtreeBox();
+        newbox->nbitems=0;
+        newbox->box[0]=NULL;
+        newbox->box[1]=NULL;
+        newbox->box[2]=NULL;
+        newbox->box[3]=NULL;
+
+        /*Add to container*/
+        this->boxcontainer->AddObject(newbox);
+        NbQuadtreeBox++;
+
+        /*currentbox now points toward next quadtree box*/
+        return newbox;
+}/*}}}*/
+/*FUNCTION Quadtree::QuadtreeColoring{{{1*/
+void Quadtree::QuadtreeColoring(double* A,int *xi,int *yi,int n){
+
+        const int     MaxDeep  = 30;
+        QuadtreeBox **pb = NULL;
+        QuadtreeBox  *b  = NULL;
+        int          level;
+
+        for(int i=0;i<n;i++){
+
+                /*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)){ 
+
+                        /*Color matrix onces more*/
+                        A[i]+=1;
+
+                        /*Go to one box deeper*/
+                        level>>=1;
+                        pb = &b->box[IJ(xi[i],yi[i],level)];
+                }
+        }
+}/*}}}*/

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

@@ -41 +41 @@
                 Quadtree();
                 ~Quadtree();
-                void Add(Observation* observation);
-                QuadtreeBox* NewQuadtreeBox(void);
+                void         Echo(void);
+                void         Add(Observation *observation);
+                void         QuadtreeColoring(double *A,int *xi,int *yi,int n);
+                QuadtreeBox *NewQuadtreeBox(void);
 };
 #endif //_QUADTREEK_H

issm/trunk-jpl/src/c/shared/Elements/elements.h

@@ -47 +47 @@
         printf("\n");
 }
+inline void printbinary(int n) {
+        unsigned int i;
+        i=1<<(sizeof(n)*8-1);
+
+        while (i>0) {
+                if (n&i)
+                 printf("1");
+                else
+                 printf("0");
+                i>>=1;
+        }
+}
 
 #endif //ifndef _SHARED_ELEMENTS_H_