![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 12210
- Timestamp:
- 05/04/12 12:26:38 (13 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 2 added
- 5 edited
-
Container/Observations.cpp (added)
-
Container/Observations.h (added)
-
Makefile.am (modified) (1 diff)
-
modules/Krigingx/Krigingx.cpp (modified) (5 diffs)
-
modules/Krigingx/Krigingx.h (modified) (1 diff)
-
objects/Kriging/Quadtree.cpp (modified) (3 diffs)
-
objects/Kriging/Quadtree.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Makefile.am
r12207 r12210 610 610 #}}} 611 611 #Kriging sources {{{1 612 kriging_sources = ./objects/Kriging/Variogram.h \ 613 ./objects/Kriging/GaussianVariogram.h\ 614 ./objects/Kriging/GaussianVariogram.cpp\ 615 ./objects/Kriging/ExponentialVariogram.h\ 616 ./objects/Kriging/ExponentialVariogram.cpp\ 617 ./objects/Kriging/SphericalVariogram.h\ 618 ./objects/Kriging/SphericalVariogram.cpp\ 619 ./objects/Kriging/Quadtree.h\ 620 ./objects/Kriging/Quadtree.cpp\ 621 ./objects/Kriging/Observation.h\ 622 ./objects/Kriging/Observation.cpp\ 623 ./modules/Krigingx/Krigingx.cpp\ 624 ./modules/Krigingx/Krigingx.h 612 kriging_sources = ./Container/Observations.h\ 613 ./Container/Observations.cpp\ 614 ./objects/Kriging/Variogram.h \ 615 ./objects/Kriging/GaussianVariogram.h\ 616 ./objects/Kriging/GaussianVariogram.cpp\ 617 ./objects/Kriging/ExponentialVariogram.h\ 618 ./objects/Kriging/ExponentialVariogram.cpp\ 619 ./objects/Kriging/SphericalVariogram.h\ 620 ./objects/Kriging/SphericalVariogram.cpp\ 621 ./objects/Kriging/Quadtree.h\ 622 ./objects/Kriging/Quadtree.cpp\ 623 ./objects/Kriging/Observation.h\ 624 ./objects/Kriging/Observation.cpp\ 625 ./modules/Krigingx/Krigingx.cpp\ 626 ./modules/Krigingx/Krigingx.h 625 627 626 628 #}}} -
issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
r12207 r12210 8 8 #include "../../toolkits/toolkits.h" 9 9 #include "../../objects/objects.h" 10 #include "../../Container/Observations.h" 10 11 #include "../modules.h" 11 12 … … 21 22 22 23 /*Intermediaries*/ 23 int i,j,n_obs;24 double numerator,denominator,ratio;25 double *x = NULL;26 double *y = NULL;27 double *obs = NULL;28 double *Gamma = NULL;29 double *GinvG0 = NULL;30 double *Ginv1 = NULL;31 double *GinvZ = NULL;32 double *gamma0 = NULL;33 double *ones = NULL;34 Variogram *variogram = NULL;35 DataSet*observations = NULL;24 int i,j,n_obs; 25 double numerator,denominator,ratio; 26 double *x = NULL; 27 double *y = NULL; 28 double *obs = NULL; 29 double *Gamma = NULL; 30 double *GinvG0 = NULL; 31 double *Ginv1 = NULL; 32 double *GinvZ = NULL; 33 double *gamma0 = NULL; 34 double *ones = NULL; 35 Variogram *variogram = NULL; 36 Observations *observations = NULL; 36 37 37 38 /*Get Variogram from Options*/ … … 39 40 40 41 /*Process observation dataset*/ 41 ProcessObservations(&observations,obs_list,obs_x,obs_y,obs_length);42 observations=new Observations(obs_list,obs_x,obs_y,obs_length); 42 43 43 44 /*Allocation output*/ … … 50 51 51 52 /*Get list of observations for current point*/ 52 ObservationList(&x,&y,&obs,&n_obs,observations,x_interp[idx],y_interp[idx]);53 observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx]); 53 54 54 55 /*Allocate intermediary matrix and vectors*/ … … 121 122 *pvariogram = variogram; 122 123 }/*}}}*/ 123 void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n){/*{{{*/124 125 int i;126 DataSet* observations = NULL;127 128 /*Initialize Observation Dataset*/129 observations = new DataSet();130 131 /*Add observations one by one*/132 for(i=0;i<n;i++){133 observations->AddObject(new Observation(x[i],y[i],0,0,observations_list[i]));134 }135 136 /*Assign output pointer*/137 *pobservations = observations;138 }/*}}}*/139 void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp){/*{{{*/140 141 /*Output and Intermediaries*/142 int nobs,i;143 double *x = NULL;144 double *y = NULL;145 double *obs = NULL;146 Observation *observation = NULL;147 148 /*Get number of observations*/149 nobs = observations->Size();150 151 /*Allocate vectors*/152 x = (double*)xmalloc(nobs*sizeof(double));153 y = (double*)xmalloc(nobs*sizeof(double));154 obs = (double*)xmalloc(nobs*sizeof(double));155 156 /*Loop over all observations and fill in x, y and obs*/157 for (i=0;i<observations->Size();i++){158 observation=(Observation*)observations->GetObjectByOffset(i);159 observation->WriteXYObs(&x[i],&y[i],&obs[i]);160 }161 162 /*Assign output pointer*/163 *px=x;164 *py=y;165 *pobs=obs;166 *pnobs=nobs;167 }/*}}}*/168 124 void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/ 169 125 #ifdef _HAVE_GSL_ -
issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h
r12207 r12210 11 11 int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options); 12 12 void ProcessVariogram(Variogram **pvariogram,Options* options); 13 void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n);14 void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp);15 13 void GslSolve(double** pX,double* A,double* B,int n); 16 14 -
issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp
r12207 r12210 47 47 * we just need to look at the bits from the left to the right (See ::Add) 48 48 }}}1*/ 49 /*MACROS {{{1*/ 50 /* 51 * 52 * J j 53 * ^ ^ 54 * | | +--------+--------+ 55 * | | | | | 56 * 1X | | | 2 | 3 | 57 * | | | | | 58 * | | +--------+--------+ 59 * | | | | | 60 * 0X | | | 0 | 1 | 61 * | | | | | 62 * | | +--------+--------+ 63 * | +-----------------------> i 64 * | 65 * |----------------------------> I 66 * X0 X1 67 * 68 * box 0 -> I=0 J=0 IJ=00 = 0 69 * box 1 -> I=1 J=0 IJ=01 = 1 70 * box 2 -> I=0 J=1 IJ=10 = 2 71 * box 3 -> I=1 J=1 IJ=11 = 3 72 */ 73 //IJ(i,j,l) returns the box number of i and j with respect to l 74 //if !j&l and !i&l -> 0 (box zero: lower left ) 75 //if !j&l and i&l -> 1 (box one: lower right) 76 //if j&l and !i&l -> 2 (box two: upper left ) 77 //if j&l and i&l -> 3 (box three:upper right) 78 #define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )) 79 /*}}}*/ 49 80 50 81 /*Constructors/Destructors*/ … … 54 85 /*Number of boxes and vertices*/ 55 86 NbQuadtreeBox=0; 56 Nb Points=0;87 NbObs=0; 57 88 58 89 /*Create container*/ … … 74 105 75 106 /*Methods*/ 107 /*FUNCTION Quadtree::Add{{{1*/ 108 void Quadtree::Add(Observation* observation){ 109 110 /*Intermediaries*/ 111 const int MaxDeep = 30; 112 int k,ij; 113 long i,j,level; 114 QuadtreeBox **pb = NULL; 115 QuadtreeBox *b = NULL; 116 117 /*Get integer coodinates*/ 118 i = observation->xi; 119 j = observation->yi; 120 121 /*Initialize level*/ 122 level=(1L<<MaxDeep);// = 2^30 123 124 /*Get inital box (the largest)*/ 125 pb=&root; 126 127 /*Find the smallest box where the observation is located*/ 128 while((b=*pb) && (b->nbitems<0)){ 129 130 /*shift b->nbitems by -1 as a counter*/ 131 b->nbitems--; 132 133 /*Go down one level = 00100 -> 00010*/ 134 level>>=1; 135 136 /*Get next subbox according to the bit value (level)*/ 137 pb = &b->box[IJ(i,j,level)]; 138 } 139 _assert_(level>0); 140 141 /*Now, try to add the vertex, if the box is full (nbitems=4), we have to divide it in 4 new boxes*/ 142 while((b=*pb) && (b->nbitems==4)){ 143 144 /*Copy the 4 observation in the current Quadtreebox*/ 145 Observation* obs[4]; 146 obs[0] = b->obs[0]; 147 obs[1] = b->obs[1]; 148 obs[2] = b->obs[2]; 149 obs[3] = b->obs[3]; 150 151 /*set nbitems as negative (now holding boxes instead of observations)*/ 152 b->nbitems = -b->nbitems; 153 154 /*Initialize the 4 pointers toward the 4 subboxes*/ 155 b->box[0]=NULL; 156 b->box[1]=NULL; 157 b->box[2]=NULL; 158 b->box[3]=NULL; 159 160 /*level = 00100 -> 00010*/ 161 level>>=1; 162 163 /*Put the four observations in the new boxes*/ 164 for (k=0;k<4;k++){ 165 166 /*Get box for observation number k*/ 167 ij=IJ(obs[k]->xi,obs[k]->yi,level); 168 QuadtreeBox *bb = b->box[ij]; 169 if(!bb){ 170 b->box[ij]=NewQuadtreeBox(); 171 bb=b->box[ij]; 172 } 173 174 /*Copy current observations*/ 175 bb->obs[bb->nbitems++] = obs[k]; 176 } 177 178 /*Get the subbox where the current observation is located*/ 179 ij=IJ(i,j,level); 180 pb=&b->box[ij]; 181 } 182 183 /*alloc the QuadtreeBox if necessary and add current observation*/ 184 b=*pb; 185 if(!b) b=NewQuadtreeBox(); 186 b->obs[b->nbitems++]=observation; 187 NbObs++; 188 189 }/*}}}*/ 76 190 /*FUNCTION Quadtree::NewQuadtreeBox {{{1*/ 77 191 Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){ -
issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h
r12207 r12210 2 2 #ifndef _QUADTREEK_H 3 3 #define _QUADTREEK_H 4 5 class Observation; 4 6 5 7 class Quadtree{ … … 17 19 union{ 18 20 QuadtreeBox *box[4]; 19 int index[4];21 Observation *obs[4]; 20 22 }; 21 23 … … 35 37 QuadtreeBox* root; // main box 36 38 long NbQuadtreeBox; // total number of boxes 37 long Nb Points; // number of points39 long NbObs; // number of points 38 40 39 41 Quadtree(); 40 42 ~Quadtree(); 43 void Add(Observation* observation); 41 44 QuadtreeBox* NewQuadtreeBox(void); 42 45 };
Note:
See TracChangeset
for help on using the changeset viewer.
