Changeset 12225

Timestamp:

05/08/12 11:10:50 (13 years ago)

Author:

Mathieu Morlighem

Message:

Better Quadtree with new fields (box center and length)

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (3 diffs)
• Container/Observations.h (modified) (1 diff)
• matlab/io/FetchMatlabData.cpp (modified) (1 diff)
• matlab/io/OptionParse.cpp (modified) (2 diffs)
• objects/Kriging/Quadtree.cpp (modified) (7 diffs)
• objects/Kriging/Quadtree.h (modified) (3 diffs)

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

@@ -35 +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       xmin,xmax,ymin,ymax;
         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];
+        xmin=x[0]; ymin=y[0];
+        xmax=x[0]; ymax=y[0];
         for(i=1;i<n;i++){
                 xmin=min(xmin,x[i]); ymin=min(ymin,y[i]);
@@ -56 +52 @@
         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);
+        /*Initialize Quadtree*/
+        this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,30);
 
         /*Add observations one by one*/
         for(i=0;i<n;i++){
-                xi=int(coefficient*(x[i]-xmin));
-                yi=int(coefficient*(y[i]-ymin));
+                this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
                 observation = new Observation(x[i],y[i],xi,yi,observations_list[i]);
                 this->quadtree->Add(observation);
@@ -115 +105 @@
 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));
+        int xi,yi;
+
         for(int i=0;i<n;i++){
-                xi[i]=int(coefficient*(x[i]-xmin));
-                yi[i]=int(coefficient*(y[i]-ymin));
+                this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
+                this->quadtree->QuadtreeColoring(&A[i],xi,yi);
         }
 
-        /*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

@@ -13 +13 @@
         private:
                 Quadtree* quadtree;
-                double    coefficient; //For integer coordinate conversion
-                double    xmin;
-                double    ymin;
-                double    xmax;
-                double    ymax;
 
         public:

issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp

@@ -404 +404 @@
         else if (mxIsClass(dataref,"char") ){
 
-                /*Check dataref is not pointing to NaN: */
-                if ( mxIsNaN(*(mxGetPr(dataref))) && (mxGetNumberOfElements(dataref)==1) ){
-                        outmatrix_numel=0;
-                        outmatrix_ndims=0;
-                        outmatrix_size =NULL;
-                        outmatrix=NULL;
-                }
-                else{
-
-                        /*Convert matlab n-dim array to char* matrix: */
-                        MatlabNArrayToNArray(&outmatrix,&outmatrix_numel,&outmatrix_ndims,&outmatrix_size,dataref);
-                }
+                /*Convert matlab n-dim array to char* matrix: */
+                MatlabNArrayToNArray(&outmatrix,&outmatrix_numel,&outmatrix_ndims,&outmatrix_size,dataref);
         }
         else{

issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp

@@ -41 +41 @@
         /*check and parse the name  */
         ological=new OptionLogical;
-        ological->name =(char *) xmalloc((strlen(name)+1)*sizeof(char));
+        ological->name =(char*)xmalloc((strlen(name)+1)*sizeof(char));
         memcpy(ological->name,name,(strlen(name)+1)*sizeof(char));
 
@@ -59 +59 @@
 
         /*check and parse the name  */
-        ochar=new OptionChar;
-        ochar->name =(char *) xmalloc((strlen(name)+1)*sizeof(char));
+        ochar=new OptionChar();
+        ochar->name =(char*)xmalloc((strlen(name)+1)*sizeof(char));
         memcpy(ochar->name,name,(strlen(name)+1)*sizeof(char));
 

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

@@ -80 +80 @@
 
         /*Constructors/Destructors*/
-        /*FUNCTION Quadtree::Quadtree(){{{1*/
-        Quadtree::Quadtree(){
-
-                /*Initialize fields*/
-                this->MaxDepth=30;
-                this->NbQuadtreeBox=0;
-                this->NbObs=0;
-
-                /*Create container*/
-                this->boxcontainer=new DataSet();
-
-                /*Create Root, pointer toward the main box*/
-                this->root=NewQuadtreeBox();
-
-                }
-        /*}}}1*/
+/*FUNCTION Quadtree::Quadtree(){{{1*/
+Quadtree::Quadtree(){
+        _error_("Constructor not supported");
+
+}
+/*}}}1*/
+/*FUNCTION Quadtree::Quadtree(double xmin,double xmax,double ymin,double ymax,int maxdepth){{{1*/
+Quadtree::Quadtree(double xmin,double xmax,double ymin,double ymax,int maxdepth){
+
+        /*Intermediaries*/
+        double length;
+
+        /*Initialize fields*/
+        this->MaxDepth=maxdepth;
+        this->NbQuadtreeBox=0;
+        this->NbObs=0;
+
+        /*Create container*/
+        this->boxcontainer=new DataSet();
+
+        /*Create Root, pointer toward the main box*/
+        length=max(xmax-xmin,ymax-ymin);
+        this->root=NewQuadtreeBox(xmin+length/2,ymin+length/2,length);
+}
+/*}}}1*/
         /*FUNCTION Quadtree::~Quadtree(){{{1*/
         Quadtree::~Quadtree(){
@@ -110 +119 @@
 
         /*Intermediaries*/
-        int          xi,yi,ij,level;
-        QuadtreeBox **pb = NULL;
-        QuadtreeBox  *b  = NULL;
-        QuadtreeBox  *bb = NULL;
+        int          xi,yi,ij,level,levelbin;
+        QuadtreeBox **pbox    = NULL; // pointer toward current box b
+        QuadtreeBox **pmaster = NULL; // pointer toward master of b
+        QuadtreeBox  *box     = NULL; // current box b
+        QuadtreeBox  *slave   = NULL; // suslaveox of b (if necessary)
         Observation  *obs[4];
 
@@ -120 +130 @@
         yi = observation->yi;
 
-        /*Initialize level*/
-        level=(1L<<this->MaxDepth);// = 2^30
+        /*Initialize levels*/
+        level    = 0;
+        levelbin = (1L<<this->MaxDepth);// = 2^30
 
         /*Get inital box (the largest)*/
-        pb=&root;
+        pmaster = &root;
+        pbox    = &root;
 
         /*Find the smallest box where the observation is located*/
-        while((b=*pb) && (b->nbitems<0)){ 
-
-                /*Go down one level = 00100 -> 00010*/
-                level>>=1;
-
-                /*Get next subbox according to the bit value (level)*/
-                pb = &b->box[IJ(xi,yi,level)];
-        }
-        _assert_(level>0);
+        while((box=*pbox) && (box->nbitems<0)){ 
+
+                /*Go down one level (levelbin = 00100 -> 00010)*/
+                levelbin>>=1; level+=1; _assert_(level<this->MaxDepth);
+
+                /*Get next box according to the bit value (levelbin)*/
+                pmaster = pbox;
+                pbox    = &box->box[IJ(xi,yi,levelbin)];
+        }
+        _assert_(levelbin>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)){
+        while((box=*pbox) && (box->nbitems==4)){
 
                 /*Copy the 4 observation in the current Quadtreebox*/
-                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;
+                obs[0] = box->obs[0];
+                obs[1] = box->obs[1];
+                obs[2] = box->obs[2];
+                obs[3] = box->obs[3];
+
+                /*set nbitems as -1 (now holding boxes instead of observations)*/
+                box->nbitems = -1;
+                box->box[0]  = NULL;
+                box->box[1]  = NULL;
+                box->box[2]  = NULL;
+                box->box[3]  = NULL;
+
+                /*Go down one level (levelbin = 00010 -> 00001)*/
+                levelbin>>=1; level+=1; _assert_(level<this->MaxDepth);
 
                 /*Put the four observations in the new boxes*/
@@ -162 +173 @@
 
                         /*Get box for observation number k*/
-                        ij = IJ(obs[k]->xi,obs[k]->yi,level);
-                        bb = b->box[ij];
-                        if(!bb){
-                                b->box[ij]=NewQuadtreeBox();
-                                bb=b->box[ij];
+                        ij    = IJ(obs[k]->xi,obs[k]->yi,levelbin);
+                        slave = box->box[ij];
+                        if(!slave){
+                                box->box[ij] = NewQuadtreeBox(box,ij);
+                                slave        = box->box[ij];
                         }
-                        bb->obs[bb->nbitems++] = obs[k];
+                        slave->obs[slave->nbitems++] = obs[k];
                 }
 
-                /*Get the subbox where the current observation is located*/
-                ij=IJ(xi,yi,level);
-                pb=&b->box[ij];
+                /*Get the suslaveox where the current observation is located*/
+                ij      = IJ(xi,yi,levelbin);
+                pmaster = pbox;
+                pbox    = &box->box[ij];
         }
 
         /*alloc the QuadtreeBox if necessary and add current observation*/
-        b=*pb;
-        if(!b){
-                b=*pb=NewQuadtreeBox();
-        }
-        b->obs[b->nbitems++]=observation;
+        box = *pbox;
+        if(!box){
+                ij  = IJ(xi,yi,levelbin);
+                box = *pbox = NewQuadtreeBox(*pmaster,ij);
+        }
+        box->obs[box->nbitems++]=observation;
         NbObs++;
 
@@ -189 +202 @@
 
         printf("Quadtree:\n");
-        printf("   MaxDepth      = %i\n",MaxDepth);
-        printf("   NbQuadtreeBox = %i\n",NbQuadtreeBox);
-        printf("   NbObs         = %i\n",NbObs);
-        printf("   root          = %p\n",root);
-
-}/*}}}*/
-/*FUNCTION Quadtree::NewQuadtreeBox {{{1*/
-Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){
+        printf("   MaxDepth      = %i\n",this->MaxDepth);
+        printf("   NbQuadtreeBox = %i\n",this->NbQuadtreeBox);
+        printf("   NbObs         = %i\n",this->NbObs);
+        printf("   root          = %p\n",this->root);
+
+}/*}}}*/
+/*FUNCTION Quadtree::DeepEcho{{{1*/
+void  Quadtree::DeepEcho(void){
+
+        printf("Quadtree:\n");
+        printf("   MaxDepth      = %i\n",this->MaxDepth);
+        printf("   NbQuadtreeBox = %i\n",this->NbQuadtreeBox);
+        printf("   NbObs         = %i\n",this->NbObs);
+        printf("   root          = %p\n",this->root);
+        boxcontainer->Echo();
+
+}/*}}}*/
+/*FUNCTION Quadtree::IntergerCoordinates{{{1*/
+void  Quadtree::IntergerCoordinates(int *xi,int *yi,double x,double y){
+
+        /*Intermediaries*/
+        double coefficient;
+        double xmin,ymin;
+
+        /*Checks in debugging mode*/
+        _assert_(xi && yi);
+        _assert_(this->root);
+
+        /*coeffIcoor is the coefficient used for integer coordinates:
+         *                (x-xmin)
+         * xi = (2^30 -1) --------- 
+         *                 length
+         * coefficient = (2^30 -1)/length
+         */
+        coefficient = double((1L<<this->MaxDepth)-1)/(this->root->length);
+        xmin        = this->root->xcenter - this->root->length/2;
+        ymin        = this->root->ycenter - this->root->length/2;
+
+        *xi=int(coefficient*(x - xmin));
+        *yi=int(coefficient*(y - ymin));
+}/*}}}*/
+/*FUNCTION Quadtree::NewQuadtreeBox(double xcenter,double ycenter,double length){{{1*/
+Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(double xcenter,double ycenter,double length){
 
         /*Output*/
         QuadtreeBox* newbox=NULL;
+
+        /*Create and initialize a new box*/
+        newbox=new QuadtreeBox();
+        newbox->nbitems=0;
+        newbox->xcenter=xcenter;
+        newbox->ycenter=ycenter;
+        newbox->length=length;
+        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::NewQuadtreeBox(QuadtreeBox* master,int index) {{{1*/
+Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(QuadtreeBox* master,int index){
+
+        /*Output*/
+        QuadtreeBox* newbox=NULL;
+
+        /*Checks in debugging mode*/
+        _assert_(master);
 
         /*Create and initialize a new box*/
@@ -208 +283 @@
         newbox->box[2]=NULL;
         newbox->box[3]=NULL;
+        switch(index){
+                case 0:
+                        newbox->xcenter=master->xcenter - master->length/4;
+                        newbox->ycenter=master->ycenter - master->length/4;
+                        break;
+                case 1:
+                        newbox->xcenter=master->xcenter + master->length/4;
+                        newbox->ycenter=master->ycenter - master->length/4;
+                        break;
+                case 2:
+                        newbox->xcenter=master->xcenter - master->length/4;
+                        newbox->ycenter=master->ycenter + master->length/4;
+                        break;
+                case 3:
+                        newbox->xcenter=master->xcenter + master->length/4;
+                        newbox->ycenter=master->ycenter + master->length/4;
+                        break;
+                default:
+                        _error_("Case %i not supported",index);
+        }
+        newbox->length=master->length/2;
 
         /*Add to container*/
@@ -217 +313 @@
 }/*}}}*/
 /*FUNCTION Quadtree::QuadtreeColoring{{{1*/
-void Quadtree::QuadtreeColoring(double* A,int *xi,int *yi,int n){
-
-        QuadtreeBox **pb = NULL;
-        QuadtreeBox  *b  = NULL;
-        int          level;
-
-        for(int i=0;i<n;i++){
-
-                /*Initialize level*/
-                level=(1L<<this->MaxDepth);// = 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)];
-                }
-        }
-}/*}}}*/
+void Quadtree::QuadtreeColoring(double* A,int xi,int yi){
+
+        QuadtreeBox **pbox = NULL;
+        QuadtreeBox  *box  = NULL;
+        int           level,levelbin;
+
+        /*Initialize levels*/
+        level    = 0;
+        levelbin = (1L<<this->MaxDepth);// = 2^30
+
+        /*Get inital box (the largest)*/
+        pbox=&root;
+
+        /*Find the smallest box where this point is located*/
+        while((box=*pbox) && (box->nbitems<0)){ 
+
+                levelbin>>=1; level+=1; _assert_(level<this->MaxDepth);
+
+                pbox = &box->box[IJ(xi,yi,levelbin)];
+        }
+        if(box && box->nbitems>0){
+                /*This box is not empty, add one level*/
+                level+=1;
+        }
+
+        *A=level;
+}/*}}}*/
+
+/*QuadtreeBox methos*/
+/*FUNCTION QuadtreeBox::Echo{{{1*/
+void  Quadtree::QuadtreeBox::Echo(void){
+
+        printf("QuadtreeBox:\n");
+        printf("   nbitems = %i\n",this->nbitems);
+        printf("   xcenter = %g\n",this->xcenter);
+        printf("   ycenter = %g\n",this->ycenter);
+        printf("   length  = %g\n",this->length);
+
+}/*}}}*/

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

@@ -16 +16 @@
                 class QuadtreeBox: public Object{ 
                         public:
-                                int nbitems; // number of current vertices in the box
+                                int    nbitems; // number of current vertices in the box
+                                double xcenter; // x position of the center (double)
+                                double ycenter; // x position of the center (double)
+                                double length;  // width of the box
                                 union{
                                         QuadtreeBox *box[4];
@@ -23 +26 @@
 
                                 /*Object functions (Needed because the Quadtree uses a Container*/
-                                void    Echo()      {_error_("not implemented yet"); };
+                                void    Echo();
                                 void    DeepEcho()  {_error_("not implemented yet"); };
                                 int     Id()        {_error_("not implemented yet"); };
@@ -35 +38 @@
 
         public:
-                int          MaxDepth;      // maximum number of subdivision
-                QuadtreeBox *root;          // main box
-                long         NbQuadtreeBox; // total number of boxes
-                long         NbObs;         // number of points
+                int          MaxDepth;          // maximum number of subdivision
+                double       coordconversion;   // Coefficient to convert coordinates to integer
+                QuadtreeBox *root;              // main box
+                long         NbQuadtreeBox;     // total number of boxes
+                long         NbObs;             // number of points
 
                 Quadtree();
+                Quadtree(double xmin,double xmax,double ymin,double ymax,int maxdepth_in);
                 ~Quadtree();
+                void         Add(Observation *observation);
+                void         DeepEcho(void);
                 void         Echo(void);
-                void         Add(Observation *observation);
-                void         QuadtreeColoring(double *A,int *xi,int *yi,int n);
-                QuadtreeBox *NewQuadtreeBox(void);
+                void         IntergerCoordinates(int *xi,int *yi,double x,double y);
+                QuadtreeBox *NewQuadtreeBox(double xcenter,double ycenter,double length);
+                QuadtreeBox *NewQuadtreeBox(QuadtreeBox* master,int index);
+                void         QuadtreeColoring(double *A,int xi,int yi);
 };
 #endif //_QUADTREEK_H