Changeset 12229

Timestamp:

05/10/12 10:41:54 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added 2 options: boxlength and searchrange

Location:

issm/trunk-jpl/src/c

Files:

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

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

@@ -36 +36 @@
 
         /*Intermediaries*/
-        int          i,maxdepth;
+        int          i,maxdepth,level,counter;
         int          xi,yi;
         double       xmin,xmax,ymin,ymax;
@@ -63 +63 @@
         }
 
-
         /*Initialize Quadtree*/
-        printf("Generating quadtree with a maximum box size %g (depth=%i)...",minlength,maxdepth);
+        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;
+        counter = 0;
         for(i=0;i<n;i++){
                 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
                 this->quadtree->QuadtreeDepth2(&level,xi,yi);
-                if((int)level >= maxdepth){
+                if((int)level <= maxdepth){
+                        observation = new Observation(x[i],y[i],xi,yi,counter++,observations_list[i]);
+                        this->quadtree->Add(observation);
+                        this->AddObject(observation);
                 }
                 else{
-                        observation = new Observation(x[i],y[i],xi,yi,observations_list[i]);
-                        this->quadtree->Add(observation);
-                        this->AddObject(observation);
+                        //We need to average with the current observations (not done yet)
                 }
         }
@@ -93 +93 @@
 /*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,double range){
 
         /*Output and Intermediaries*/
-        int          nobs,i;
-        double      *x           = NULL;
-        double      *y           = NULL;
-        double      *obs         = NULL;
-        Observation *observation = NULL;
+        int          nobs,i,index;
+        double      *x            = NULL;
+        double      *y            = NULL;
+        double      *obs          = NULL;
+        Observation *observation  = NULL;
+        int         *indices      = NULL;
 
-        /*Get number of observations*/
-        nobs = this->Size();
+        /*Treat range*/
+        if(range==0){
+                range=this->quadtree->root->length;
+        }
 
-        /*Allocate vectors*/
-        x   = (double*)xmalloc(nobs*sizeof(double));
-        y   = (double*)xmalloc(nobs*sizeof(double));
-        obs = (double*)xmalloc(nobs*sizeof(double));
+        /*Find all observations that are in range*/
+        this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range);
 
-        /*Loop over all observations and fill in x, y and obs*/
-        for (i=0;i<nobs;i++){
-                observation=(Observation*)this->GetObjectByOffset(i);
-                observation->WriteXYObs(&x[i],&y[i],&obs[i]);
+        if(nobs==0){
+                /*No observation found, double range*/
+                printf("No observation found within range, doubling range\n");
+                xfree((void**)&indices);
+                this->ObservationList(&x,&y,&obs,&nobs,x_interp,y_interp,range*2);
+        }
+        else{
+                if(nobs>1000) printf("Taking more than 1000 observations\n");
+                /*Allocate vectors*/
+                x   = (double*)xmalloc(nobs*sizeof(double));
+                y   = (double*)xmalloc(nobs*sizeof(double));
+                obs = (double*)xmalloc(nobs*sizeof(double));
+
+                /*Loop over all observations and fill in x, y and obs*/
+                for (i=0;i<nobs;i++){
+                        observation=(Observation*)this->GetObjectByOffset(indices[i]);
+                        observation->WriteXYObs(&x[i],&y[i],&obs[i]);
+                }
         }
 
         /*Assign output pointer*/
+        xfree((void**)&indices);
         *px=x;
         *py=y;

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

@@ -23 +23 @@
 
                 /*Methods*/
-                void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp);
+                void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double range);
                 void QuadtreeColoring(double* A,double *x,double *y,int n);
 

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

@@ -23 +23 @@
         /*Intermediaries*/
         int           i,j,n_obs;
+        double        range;
         double        numerator,denominator,ratio;
         double       *x            = NULL;
@@ -33 +34 @@
         double       *gamma0       = NULL;
         double       *ones         = NULL;
+        char         *output       = NULL;
         Variogram    *variogram    = NULL;
         Observations *observations = NULL;
@@ -38 +40 @@
         /*Get Variogram from Options*/
         ProcessVariogram(&variogram,options);
+        options->Get(&range,"searchrange",0.);
 
         /*Process observation dataset*/
@@ -44 +47 @@
         /*Allocation output*/
         predictions =(double*)xmalloc(n_interp*sizeof(double));
+        for(i=0;i<n_interp;i++) predictions[i]=0;
 
-        for(i=0;i<n_interp;i++) predictions[i]=0;
-        observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
-        *ppredictions=predictions;
-        return 1;
+        /*Get output*/
+        options->Get(&output,"output","quadtree");
+
+        if(strcmp(output,"quadtree")==0){
+                observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
+                *ppredictions=predictions;
+                return 1;
+        }
 
         /*Loop over all interpolations*/
@@ -56 +64 @@
 
                 /*Get list of observations for current point*/
-                observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx]);
+                observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],range);
 
                 /*Allocate intermediary matrix and vectors*/
@@ -105 +113 @@
         delete variogram;
         delete observations;
+        xfree((void**)&output);
         *ppredictions=predictions;
         return 1;

issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h ¶

@@ -1 +1 @@
 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, BamgQuadtree.h)*/
-#ifndef _QUADTREE_H
-#define _QUADTREE_H
+#ifndef _BAMGQUADTREE_H
+#define _BAMGQUADTREE_H
 
 #include "./include.h"

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

@@ -12 +12 @@
 }
 /*}}}*/
-/*FUNCTION Observation::Observation(double x,double y,int xi,int yi,double value){{{1*/
-Observation::Observation(double x_in,double y_in,int xi_in,int yi_in,double value_in){
+/*FUNCTION Observation::Observation(double x,double y,int xi,int yi,int index,double value){{{1*/
+Observation::Observation(double x_in,double y_in,int xi_in,int yi_in,int index_in,double value_in){
 
         this->x     = x_in;
@@ -19 +19 @@
         this->xi    = xi_in;
         this->yi    = yi_in;
+        this->index = index_in;
         this->value = value_in;
 
@@ -36 +37 @@
 
         printf("Observation\n");
+        printf("   index : %i\n",this->index);
         printf("   x     : %g\n",this->x);
         printf("   y     : %g\n",this->y);

issm/trunk-jpl/src/c/objects/Kriging/Observation.h ¶

@@ -13 +13 @@
                 double x,y;
                 int    xi,yi;
+                int    index;
                 double value;
 
                 /*Observation constructors, destructors*/
                 Observation();
-                Observation(double x_in,double y_in,int xi_in,int yi_in,double value_in);
+                Observation(double x_in,double y_in,int xi_in,int yi_in,int index_in,double value_in);
                 ~Observation();
 
@@ -31 +32 @@
                 void WriteXYObs(double* px,double* py,double* pobs);
 };
-#endif  /* _EXPONENTIALVARIOGRAM_H */
+#endif  /* _OBSERVATION_*/

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

@@ -393 +393 @@
         *A=level;
 }/*}}}*/
+/*FUNCTION Quadtree::RangeSearch{{{1*/
+void Quadtree::RangeSearch(int **pindices,int *pnobs,double x,double y,double range){
+
+        /*Intermediaries*/
+        int  nobs;
+        int *indices = NULL;
+
+        /*Allocate indices (maximum by default*/
+        indices = (int*)xmalloc(this->NbObs*sizeof(int));
+        nobs = 0;
+
+        this->root->RangeSearch(indices,&nobs,x,y,range);
+
+        /*Clean-up and return*/
+        *pnobs=nobs;
+        *pindices=indices;
+
+}/*}}}*/
 
 /*QuadtreeBox methos*/
@@ -405 +423 @@
 
 }/*}}}*/
+/*FUNCTION QuadtreeBox::IsWithinRange{{{1*/
+int Quadtree::QuadtreeBox::IsWithinRange(double x,double y,double range){
+
+        /*Return 0 if the 2 boxes do not overlap*/
+        if(this->xcenter+this->length/2 < x-range) return 0;
+        if(this->xcenter-this->length/2 > x+range) return 0;
+        if(this->ycenter+this->length/2 < y-range) return 0;
+        if(this->ycenter-this->length/2 > y+range) return 0;
+
+        /*Return 2 if the this box is included in the range*/
+        if(
+                                this->xcenter+this->length/2 <= x+range &&
+                                this->ycenter+this->length/2 <= y+range &&
+                                this->xcenter-this->length/2 >= x-range &&
+                                this->ycenter-this->length/2 >= y-range
+          ) return 2;
+
+        /*This is a simple overlap*/
+        return 1;
+
+}/*}}}*/
+/*FUNCTION QuadtreeBox::RangeSearch{{{1*/
+void Quadtree::QuadtreeBox::RangeSearch(int* indices,int *pnobs,double x,double y,double range){
+
+        /*Intermediaries*/
+        int i,nobs;
+
+        /*Recover current number of observations*/
+        nobs = *pnobs;
+
+        switch(this->IsWithinRange(x,y,range)){
+                case 0:
+                        /*If this box is not within range, return*/
+                        break;
+                case 2:
+                        /*This box is included in range*/
+                        this->WriteObservations(indices,&nobs);
+                        break;
+                case 1:
+                        /*This box is partly included*/
+                        if(this->nbitems>0){
+                                /*If this box has only observations, add indices that are within range*/
+                                for(i=0;i<this->nbitems;i++){
+                                        if(fabs(this->obs[i]->x-x) <= range && fabs(this->obs[i]->y-y) <= range){
+                                                indices[nobs++]=this->obs[i]->index;
+                                        }
+                                }
+                        }
+                        else{
+                                /*This box points toward boxes*/
+                                if(this->box[0]) this->box[0]->RangeSearch(indices,&nobs,x,y,range);
+                                if(this->box[1]) this->box[1]->RangeSearch(indices,&nobs,x,y,range);
+                                if(this->box[2]) this->box[2]->RangeSearch(indices,&nobs,x,y,range);
+                                if(this->box[3]) this->box[3]->RangeSearch(indices,&nobs,x,y,range);
+                        }
+                        break;
+                default:
+                        _error_("Case %i not supported",this->IsWithinRange(x,y,range));
+        }
+
+        /*Assign output pointers: */
+        *pnobs=nobs;
+}/*}}}*/
+/*FUNCTION QuadtreeBox::WriteObservations{{{1*/
+void Quadtree::QuadtreeBox::WriteObservations(int* indices,int *pnobs){
+
+        /*Intermediaries*/
+        int i,nobs;
+
+        /*Recover current number of observations*/
+        nobs = *pnobs;
+
+        if(this->nbitems>0){
+                /*If this box has only observations, add all indices*/
+                for(i=0;i<this->nbitems;i++){
+                        indices[nobs++]=this->obs[i]->index;
+                }
+        }
+        else{
+                /*This box points toward boxes, */
+                if(this->box[0]) this->box[0]->WriteObservations(indices,&nobs);
+                if(this->box[1]) this->box[1]->WriteObservations(indices,&nobs);
+                if(this->box[2]) this->box[2]->WriteObservations(indices,&nobs);
+                if(this->box[3]) this->box[3]->WriteObservations(indices,&nobs);
+        }
+
+        /*Assign output pointers: */
+        *pnobs=nobs;
+}/*}}}*/

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

@@ -1 +1 @@
 
-#ifndef _QUADTREEK_H
-#define _QUADTREEK_H
+#ifndef _QUADTREE_H
+#define _QUADTREE_H
 
 class Observation;
@@ -32 +32 @@
                                 int     ObjectEnum(){_error_("not implemented yet"); };
                                 Object *copy()      {_error_("not implemented yet"); };
+
+                                /*Methods*/
+                                int IsWithinRange(double x,double y,double range);
+                                void RangeSearch(int* indices,int *pnobs,double x,double y,double range);
+                                void WriteObservations(int* indices,int *pnobs);
+
                 };
 
@@ -39 +45 @@
         public:
                 int          MaxDepth;          // maximum number of subdivision
-                double       coordconversion;   // Coefficient to convert coordinates to integer
                 QuadtreeBox *root;              // main box
                 long         NbQuadtreeBox;     // total number of boxes
@@ -55 +60 @@
                 void         QuadtreeDepth(int *A,int xi,int yi);
                 void         QuadtreeDepth2(int *A,int xi,int yi);
+                void         RangeSearch(int **pindices,int *pnobs,double x,double y,double range);
 };
-#endif //_QUADTREEK_H
+#endif //_QUADTREE_H