Changeset 12226

Timestamp:

05/08/12 13:35:24 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added boxlength option

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (3 diffs)
• Container/Observations.h (modified) (2 diffs)
• modules/Krigingx/Krigingx.cpp (modified) (1 diff)
• objects/Kriging/Quadtree.cpp (modified) (2 diffs)
• objects/Kriging/Quadtree.h (modified) (1 diff)

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

@@ -32 +32 @@
 }
 /*}}}*/
-/*FUNCTION Observations::Observations(double* observations_list,double* x,double* y,int n){{{*/
-Observations::Observations(double* observations_list,double* x,double* y,int n){
+/*FUNCTION Observations::Observations(double* observations_list,double* x,double* y,int n,Options* options){{{*/
+Observations::Observations(double* observations_list,double* x,double* y,int n,Options* options){
 
         /*Intermediaries*/
-        int          i;
+        int          i,maxdepth;
         int          xi,yi;
         double       xmin,xmax,ymin,ymax;
-        double       offset;
+        double       offset,minlength;
         Observation *observation = NULL;
 
@@ -52 +52 @@
         offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
 
+        /*Get Minimum box size*/
+        if(options->GetOption("boxlength")){
+                options->Get(&minlength,"boxlength");
+                if(minlength<=0)_error_("boxlength should be a positive number");
+                maxdepth=int(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2));
+        }
+        else{
+                maxdepth = 30;
+                minlength=max(xmax-xmin,ymax-ymin)/double((1L<<maxdepth)-1);
+        }
+
+
         /*Initialize Quadtree*/
-        this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,30);
+        printf("Generating quadtree with a maximum box size %g (depth=%i)...",minlength,maxdepth);
+        this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,maxdepth);
 
         /*Add observations one by one*/
+        int level;
         for(i=0;i<n;i++){
                 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);
-                this->AddObject(observation);
+                this->quadtree->QuadtreeDepth2(&level,xi,yi);
+                if((int)level >= maxdepth){
+                }
+                else{
+                        observation = new Observation(x[i],y[i],xi,yi,observations_list[i]);
+                        this->quadtree->Add(observation);
+                        this->AddObject(observation);
+                }
         }
+        printf("done\n");
 }
 /*}}}*/
@@ -105 +125 @@
 void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){
 
-        int xi,yi;
+        int xi,yi,level;
 
         for(int i=0;i<n;i++){
                 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
-                this->quadtree->QuadtreeColoring(&A[i],xi,yi);
+                this->quadtree->QuadtreeDepth(&level,xi,yi);
+                A[i]=(double)level;
         }
 

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

@@ -8 +8 @@
 class Obsevration;
 class Quadtree;
+class Options;
 
 class Observations: public DataSet{
@@ -18 +19 @@
                 /*constructors, destructors*/
                 Observations();
-                Observations(double* observations_list,double* x,double* y,int n);
+                Observations(double* observations_list,double* x,double* y,int n,Options* options);
                 ~Observations();
 

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

@@ -40 +40 @@
 
         /*Process observation dataset*/
-        observations=new Observations(obs_list,obs_x,obs_y,obs_length);
+        observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
 
         /*Allocation output*/

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

@@ -312 +312 @@
         return newbox;
 }/*}}}*/
-/*FUNCTION Quadtree::QuadtreeColoring{{{1*/
-void Quadtree::QuadtreeColoring(double* A,int xi,int yi){
+/*FUNCTION Quadtree::QuadtreeDepth{{{1*/
+void Quadtree::QuadtreeDepth(int* A,int xi,int yi){
 
         QuadtreeBox **pbox = NULL;
@@ -336 +336 @@
                 /*This box is not empty, add one level*/
                 level+=1;
+        }
+
+        *A=level;
+}/*}}}*/
+/*FUNCTION Quadtree::QuadtreeDepth2{{{1*/
+void Quadtree::QuadtreeDepth2(int* 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; 
+
+                pbox = &box->box[IJ(xi,yi,levelbin)];
+        }
+        if(box && box->nbitems>0){
+                /*This box is not empty, add one level*/
+                level+=1;
+        }
+
+        /*If we were to add the vertex, get level*/
+        if(box && box->nbitems==4){
+                int ij;
+                bool flag=true;
+                while(flag){
+
+                        levelbin>>=1; level+=1;
+                        if(level>this->MaxDepth){
+                                level+=1;
+                                break;
+                        }
+
+                        /*loop over the four observations*/
+                        ij=IJ(box->obs[0]->xi,box->obs[0]->yi,levelbin);
+                        for (int k=1;k<4;k++){
+                                if(IJ(box->obs[k]->xi,box->obs[k]->yi,levelbin) != ij){
+                                        flag = false;
+                                }
+                        }
+                        if(IJ(xi,yi,levelbin)!=ij){
+                                flag = false;
+                        }
+                }
         }
 

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

@@ -53 +53 @@
                 QuadtreeBox *NewQuadtreeBox(double xcenter,double ycenter,double length);
                 QuadtreeBox *NewQuadtreeBox(QuadtreeBox* master,int index);
-                void         QuadtreeColoring(double *A,int xi,int yi);
+                void         QuadtreeDepth(int *A,int xi,int yi);
+                void         QuadtreeDepth2(int *A,int xi,int yi);
 };
 #endif //_QUADTREEK_H