Index: /issm/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12228) +++ /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12229) @@ -36,5 +36,5 @@ /*Intermediaries*/ - int i,maxdepth; + int i,maxdepth,level,counter; int xi,yi; double xmin,xmax,ymin,ymax; @@ -63,20 +63,20 @@ } - /*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;iquadtree->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,28 +93,44 @@ /*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;iGetObjectByOffset(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;iGetObjectByOffset(indices[i]); + observation->WriteXYObs(&x[i],&y[i],&obs[i]); + } } /*Assign output pointer*/ + xfree((void**)&indices); *px=x; *py=y; Index: /issm/trunk-jpl/src/c/Container/Observations.h =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.h (revision 12228) +++ /issm/trunk-jpl/src/c/Container/Observations.h (revision 12229) @@ -23,5 +23,5 @@ /*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); Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12228) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12229) @@ -23,4 +23,5 @@ /*Intermediaries*/ int i,j,n_obs; + double range; double numerator,denominator,ratio; double *x = NULL; @@ -33,4 +34,5 @@ double *gamma0 = NULL; double *ones = NULL; + char *output = NULL; Variogram *variogram = NULL; Observations *observations = NULL; @@ -38,4 +40,5 @@ /*Get Variogram from Options*/ ProcessVariogram(&variogram,options); + options->Get(&range,"searchrange",0.); /*Process observation dataset*/ @@ -44,9 +47,14 @@ /*Allocation output*/ predictions =(double*)xmalloc(n_interp*sizeof(double)); + for(i=0;iQuadtreeColoring(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,5 +64,5 @@ /*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,4 +113,5 @@ delete variogram; delete observations; + xfree((void**)&output); *ppredictions=predictions; return 1; Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h (revision 12228) +++ /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h (revision 12229) @@ -1,5 +1,5 @@ /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.h)*/ -#ifndef _QUADTREE_H -#define _QUADTREE_H +#ifndef _BAMGQUADTREE_H +#define _BAMGQUADTREE_H #include "./include.h" Index: /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp (revision 12228) +++ /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp (revision 12229) @@ -12,6 +12,6 @@ } /*}}}*/ -/*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,4 +19,5 @@ this->xi = xi_in; this->yi = yi_in; + this->index = index_in; this->value = value_in; @@ -36,4 +37,5 @@ printf("Observation\n"); + printf(" index : %i\n",this->index); printf(" x : %g\n",this->x); printf(" y : %g\n",this->y); Index: /issm/trunk-jpl/src/c/objects/Kriging/Observation.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Observation.h (revision 12228) +++ /issm/trunk-jpl/src/c/objects/Kriging/Observation.h (revision 12229) @@ -13,9 +13,10 @@ 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,3 +32,3 @@ void WriteXYObs(double* px,double* py,double* pobs); }; -#endif /* _EXPONENTIALVARIOGRAM_H */ +#endif /* _OBSERVATION_*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12228) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12229) @@ -393,4 +393,22 @@ *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,2 +423,91 @@ }/*}}}*/ +/*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;inbitems;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;inbitems;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; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12228) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12229) @@ -1,5 +1,5 @@ -#ifndef _QUADTREEK_H -#define _QUADTREEK_H +#ifndef _QUADTREE_H +#define _QUADTREE_H class Observation; @@ -32,4 +32,10 @@ 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,5 +45,4 @@ 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,4 +60,5 @@ 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