Changeset 12273

Timestamp:

05/22/12 10:54:55 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added observation averaging

Location:

issm/trunk-jpl/src/c

Files:

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

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

@@ -78 +78 @@
                 }
                 else{
-                        //We need to average with the current observations (not done yet)
+                        /*We need to average with the current observations*/
+                        this->quadtree->AddAndAverage(x[i],y[i],observations_list[i]);
                 }
         }

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

@@ -15 +15 @@
 Observation::Observation(double x_in,double y_in,int xi_in,int yi_in,int index_in,double value_in){
 
-        this->x     = x_in;
-        this->y     = y_in;
-        this->xi    = xi_in;
-        this->yi    = yi_in;
-        this->index = index_in;
-        this->value = value_in;
+        this->x      = x_in;
+        this->y      = y_in;
+        this->xi     = xi_in;
+        this->yi     = yi_in;
+        this->index  = index_in;
+        this->value  = value_in;
+        this->weight = 1.;
 
 }
@@ -42 +43 @@
         printf("   xi    : "); printbinary(this->xi); printf("\n");
         printf("   yi    : "); printbinary(this->yi); printf("\n");
+        printf("   weight: %g\n",this->weight);
         printf("   value : %g\n",this->value);
 }

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

@@ -14 +14 @@
                 int    xi,yi;
                 int    index;
+                double weight;
                 double value;
 

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

@@ -198 +198 @@
 
 }/*}}}*/
+/*FUNCTION Quadtree::AddAndAverage{{{1*/
+void Quadtree::AddAndAverage(double x,double y,double value){
+
+        QuadtreeBox **pbox = NULL;
+        QuadtreeBox  *box  = NULL;
+        int           xi,yi;
+        int           level,levelbin;
+        int           index;
+        double        length,length2;
+
+        /*Get integer coodinates*/
+        this->IntergerCoordinates(&xi,&yi,x,y);
+
+        /*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)];
+        }
+
+        /*Add obervation in this box (should be full)*/
+        if(box && box->nbitems==4){
+                index  = 0;
+                length = pow(box->obs[0]->x - x,2.) + pow(box->obs[0]->y - y,2.);
+                for(int i=1;i<4;i++){
+                        length2 = pow(box->obs[i]->x - x,2.) + pow(box->obs[i]->y - y,2.);
+                        if(length2<length){
+                                index  = i;
+                                length = length2;
+                        }
+                }
+
+                /*We found the closest observation, now average observation (do not change xi and yi to avoid round off errors*/
+                box->obs[index]->x = (box->obs[index]->weight*box->obs[index]->x + x)/(box->obs[index]->weight+1);
+                box->obs[index]->y = (box->obs[index]->weight*box->obs[index]->y + y)/(box->obs[index]->weight+1);
+                box->obs[index]->xi= int((box->obs[index]->weight*double(box->obs[index]->xi) + double(xi))/(box->obs[index]->weight+1));
+                box->obs[index]->yi= int((box->obs[index]->weight*double(box->obs[index]->yi) + double(yi))/(box->obs[index]->weight+1));
+                box->obs[index]->value   = (box->obs[index]->weight*box->obs[index]->value + value)/(box->obs[index]->weight+1);
+                box->obs[index]->weight += 1;
+        }
+        else{
+                _error_("Box is not full");
+        }
+}/*}}}*/
 /*FUNCTION Quadtree::Echo{{{1*/
 void  Quadtree::Echo(void){

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

@@ -53 +53 @@
                 ~Quadtree();
                 void         Add(Observation *observation);
+                void         AddAndAverage(double x,double y,double value);
                 void         DeepEcho(void);
                 void         Echo(void);