Index: /issm/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12210) +++ /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12210) @@ -0,0 +1,79 @@ +/* + * \file Observations.c + * \brief: implementation of the Observations class, derived from DataSet class + */ + +/*Headers: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include +#include + +#include "./DataSet.h" +#include "./Observations.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +using namespace std; +/*}}}*/ + +/*Object constructors and destructor*/ +/*FUNCTION Observations::Observations(){{{*/ +Observations::Observations(){ + return; +} +/*}}}*/ +/*FUNCTION Observations::Observations(double* observations_list,double* x,double* y,int n){{{*/ +Observations::Observations(double* observations_list,double* x,double* y,int n){ + + /*Add observations one by one*/ + for(int i=0;iAddObject(new Observation(x[i],y[i],0,0,observations_list[i])); + } + +} +/*}}}*/ +/*FUNCTION Observations::~Observations(){{{*/ +Observations::~Observations(){ + return; +} +/*}}}*/ + +/*Methods*/ +/*FUNCTION Observations::ObservationList{{{*/ +void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp){/*{{{*/ + + /*Output and Intermediaries*/ + int nobs,i; + double *x = NULL; + double *y = NULL; + double *obs = NULL; + Observation *observation = NULL; + + /*Get number of observations*/ + nobs = this->Size(); + + /*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;iGetObjectByOffset(i); + observation->WriteXYObs(&x[i],&y[i],&obs[i]); + } + + /*Assign output pointer*/ + *px=x; + *py=y; + *pobs=obs; + *pnobs=nobs; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/Container/Observations.h =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.h (revision 12210) +++ /issm/trunk-jpl/src/c/Container/Observations.h (revision 12210) @@ -0,0 +1,24 @@ +/*!\file: Observations.h + * \brief prototypes for Observations.h + */ + +#ifndef _CONTAINER_OBSERVATIONS_H_ +#define _CONTAINER_OBSERVATIONS_H_ + +class Obsevration; + +class Observations: public DataSet{ + + public: + + /*constructors, destructors*/ + Observations(); + Observations(double* observations_list,double* x,double* y,int n); + ~Observations(); + + /*Methods*/ + void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp); + +}; +#endif //ifndef _OBSERVATIONS_H_ + Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12209) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12210) @@ -610,17 +610,19 @@ #}}} #Kriging sources {{{1 -kriging_sources = ./objects/Kriging/Variogram.h \ - ./objects/Kriging/GaussianVariogram.h\ - ./objects/Kriging/GaussianVariogram.cpp\ - ./objects/Kriging/ExponentialVariogram.h\ - ./objects/Kriging/ExponentialVariogram.cpp\ - ./objects/Kriging/SphericalVariogram.h\ - ./objects/Kriging/SphericalVariogram.cpp\ - ./objects/Kriging/Quadtree.h\ - ./objects/Kriging/Quadtree.cpp\ - ./objects/Kriging/Observation.h\ - ./objects/Kriging/Observation.cpp\ - ./modules/Krigingx/Krigingx.cpp\ - ./modules/Krigingx/Krigingx.h +kriging_sources = ./Container/Observations.h\ + ./Container/Observations.cpp\ + ./objects/Kriging/Variogram.h \ + ./objects/Kriging/GaussianVariogram.h\ + ./objects/Kriging/GaussianVariogram.cpp\ + ./objects/Kriging/ExponentialVariogram.h\ + ./objects/Kriging/ExponentialVariogram.cpp\ + ./objects/Kriging/SphericalVariogram.h\ + ./objects/Kriging/SphericalVariogram.cpp\ + ./objects/Kriging/Quadtree.h\ + ./objects/Kriging/Quadtree.cpp\ + ./objects/Kriging/Observation.h\ + ./objects/Kriging/Observation.cpp\ + ./modules/Krigingx/Krigingx.cpp\ + ./modules/Krigingx/Krigingx.h #}}} Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12209) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12210) @@ -8,4 +8,5 @@ #include "../../toolkits/toolkits.h" #include "../../objects/objects.h" +#include "../../Container/Observations.h" #include "../modules.h" @@ -21,17 +22,17 @@ /*Intermediaries*/ - int i,j,n_obs; - double numerator,denominator,ratio; - double *x = NULL; - double *y = NULL; - double *obs = NULL; - double *Gamma = NULL; - double *GinvG0 = NULL; - double *Ginv1 = NULL; - double *GinvZ = NULL; - double *gamma0 = NULL; - double *ones = NULL; - Variogram *variogram = NULL; - DataSet *observations = NULL; + int i,j,n_obs; + double numerator,denominator,ratio; + double *x = NULL; + double *y = NULL; + double *obs = NULL; + double *Gamma = NULL; + double *GinvG0 = NULL; + double *Ginv1 = NULL; + double *GinvZ = NULL; + double *gamma0 = NULL; + double *ones = NULL; + Variogram *variogram = NULL; + Observations *observations = NULL; /*Get Variogram from Options*/ @@ -39,5 +40,5 @@ /*Process observation dataset*/ - ProcessObservations(&observations,obs_list,obs_x,obs_y,obs_length); + observations=new Observations(obs_list,obs_x,obs_y,obs_length); /*Allocation output*/ @@ -50,5 +51,5 @@ /*Get list of observations for current point*/ - ObservationList(&x,&y,&obs,&n_obs,observations,x_interp[idx],y_interp[idx]); + observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx]); /*Allocate intermediary matrix and vectors*/ @@ -121,49 +122,4 @@ *pvariogram = variogram; }/*}}}*/ -void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n){/*{{{*/ - - int i; - DataSet* observations = NULL; - - /*Initialize Observation Dataset*/ - observations = new DataSet(); - - /*Add observations one by one*/ - for(i=0;iAddObject(new Observation(x[i],y[i],0,0,observations_list[i])); - } - - /*Assign output pointer*/ - *pobservations = observations; -}/*}}}*/ -void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp){/*{{{*/ - - /*Output and Intermediaries*/ - int nobs,i; - double *x = NULL; - double *y = NULL; - double *obs = NULL; - Observation *observation = NULL; - - /*Get number of observations*/ - nobs = observations->Size(); - - /*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;iSize();i++){ - observation=(Observation*)observations->GetObjectByOffset(i); - observation->WriteXYObs(&x[i],&y[i],&obs[i]); - } - - /*Assign output pointer*/ - *px=x; - *py=y; - *pobs=obs; - *pnobs=nobs; -}/*}}}*/ void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/ #ifdef _HAVE_GSL_ Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12209) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12210) @@ -11,6 +11,4 @@ int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options); void ProcessVariogram(Variogram **pvariogram,Options* options); -void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n); -void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp); void GslSolve(double** pX,double* A,double* B,int n); Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12209) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12210) @@ -47,4 +47,35 @@ * we just need to look at the bits from the left to the right (See ::Add) }}}1*/ +/*MACROS {{{1*/ +/* + * + * J j + * ^ ^ + * | | +--------+--------+ + * | | | | | + * 1X | | | 2 | 3 | + * | | | | | + * | | +--------+--------+ + * | | | | | + * 0X | | | 0 | 1 | + * | | | | | + * | | +--------+--------+ + * | +-----------------------> i + * | + * |----------------------------> I + * X0 X1 + * + * box 0 -> I=0 J=0 IJ=00 = 0 + * box 1 -> I=1 J=0 IJ=01 = 1 + * box 2 -> I=0 J=1 IJ=10 = 2 + * box 3 -> I=1 J=1 IJ=11 = 3 + */ +//IJ(i,j,l) returns the box number of i and j with respect to l +//if !j&l and !i&l -> 0 (box zero: lower left ) +//if !j&l and i&l -> 1 (box one: lower right) +//if j&l and !i&l -> 2 (box two: upper left ) +//if j&l and i&l -> 3 (box three:upper right) +#define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )) +/*}}}*/ /*Constructors/Destructors*/ @@ -54,5 +85,5 @@ /*Number of boxes and vertices*/ NbQuadtreeBox=0; - NbPoints=0; + NbObs=0; /*Create container*/ @@ -74,4 +105,87 @@ /*Methods*/ +/*FUNCTION Quadtree::Add{{{1*/ +void Quadtree::Add(Observation* observation){ + + /*Intermediaries*/ + const int MaxDeep = 30; + int k,ij; + long i,j,level; + QuadtreeBox **pb = NULL; + QuadtreeBox *b = NULL; + + /*Get integer coodinates*/ + i = observation->xi; + j = observation->yi; + + /*Initialize level*/ + level=(1L<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)]; + } + _assert_(level>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)){ + + /*Copy the 4 observation in the current Quadtreebox*/ + Observation* obs[4]; + 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; + + /*Put the four observations in the new boxes*/ + for (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]; + 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); + pb=&b->box[ij]; + } + + /*alloc the QuadtreeBox if necessary and add current observation*/ + b=*pb; + if(!b) b=NewQuadtreeBox(); + b->obs[b->nbitems++]=observation; + NbObs++; + +}/*}}}*/ /*FUNCTION Quadtree::NewQuadtreeBox {{{1*/ Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){ Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12209) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12210) @@ -2,4 +2,6 @@ #ifndef _QUADTREEK_H #define _QUADTREEK_H + +class Observation; class Quadtree{ @@ -17,5 +19,5 @@ union{ QuadtreeBox *box[4]; - int index[4]; + Observation *obs[4]; }; @@ -35,8 +37,9 @@ QuadtreeBox* root; // main box long NbQuadtreeBox; // total number of boxes - long NbPoints; // number of points + long NbObs; // number of points Quadtree(); ~Quadtree(); + void Add(Observation* observation); QuadtreeBox* NewQuadtreeBox(void); };