Index: /issm/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12224) +++ /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12225) @@ -35,18 +35,14 @@ Observations::Observations(double* observations_list,double* x,double* y,int n){ - const int MaxICoor = 1073741823; // 2^30-1 =111...111 (29 times) - /*Intermediaries*/ int i; int xi,yi; + double xmin,xmax,ymin,ymax; double offset; Observation *observation = NULL; - /*Initialize Quadtree*/ - this->quadtree = new Quadtree(); - /*Get extrema*/ - this->xmin=x[0]; this->ymin=y[0]; - this->xmax=x[0]; this->ymax=y[0]; + xmin=x[0]; ymin=y[0]; + xmax=x[0]; ymax=y[0]; for(i=1;icoefficient = double(MaxICoor)/max(xmax-xmin,ymax-ymin); _assert_(coefficient>=0); + /*Initialize Quadtree*/ + this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,30); /*Add observations one by one*/ for(i=0;iquadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]); observation = new Observation(x[i],y[i],xi,yi,observations_list[i]); this->quadtree->Add(observation); @@ -115,18 +105,10 @@ void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){ - /*Convert to integer coordinates*/ - int *xi = (int*)xmalloc(n*sizeof(int)); - int *yi = (int*)xmalloc(n*sizeof(int)); + int xi,yi; + for(int i=0;iquadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]); + this->quadtree->QuadtreeColoring(&A[i],xi,yi); } - /*Call quadtree method*/ - this->quadtree->QuadtreeColoring(A,xi,yi,n); - - /*clean-up*/ - xfree((void**)&xi); - xfree((void**)&yi); - }/*}}}*/ Index: /issm/trunk-jpl/src/c/Container/Observations.h =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.h (revision 12224) +++ /issm/trunk-jpl/src/c/Container/Observations.h (revision 12225) @@ -13,9 +13,4 @@ private: Quadtree* quadtree; - double coefficient; //For integer coordinate conversion - double xmin; - double ymin; - double xmax; - double ymax; public: Index: /issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp (revision 12224) +++ /issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp (revision 12225) @@ -404,16 +404,6 @@ else if (mxIsClass(dataref,"char") ){ - /*Check dataref is not pointing to NaN: */ - if ( mxIsNaN(*(mxGetPr(dataref))) && (mxGetNumberOfElements(dataref)==1) ){ - outmatrix_numel=0; - outmatrix_ndims=0; - outmatrix_size =NULL; - outmatrix=NULL; - } - else{ - - /*Convert matlab n-dim array to char* matrix: */ - MatlabNArrayToNArray(&outmatrix,&outmatrix_numel,&outmatrix_ndims,&outmatrix_size,dataref); - } + /*Convert matlab n-dim array to char* matrix: */ + MatlabNArrayToNArray(&outmatrix,&outmatrix_numel,&outmatrix_ndims,&outmatrix_size,dataref); } else{ Index: /issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp (revision 12224) +++ /issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp (revision 12225) @@ -41,5 +41,5 @@ /*check and parse the name */ ological=new OptionLogical; - ological->name =(char *) xmalloc((strlen(name)+1)*sizeof(char)); + ological->name =(char*)xmalloc((strlen(name)+1)*sizeof(char)); memcpy(ological->name,name,(strlen(name)+1)*sizeof(char)); @@ -59,6 +59,6 @@ /*check and parse the name */ - ochar=new OptionChar; - ochar->name =(char *) xmalloc((strlen(name)+1)*sizeof(char)); + ochar=new OptionChar(); + ochar->name =(char*)xmalloc((strlen(name)+1)*sizeof(char)); memcpy(ochar->name,name,(strlen(name)+1)*sizeof(char)); Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12224) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12225) @@ -80,20 +80,29 @@ /*Constructors/Destructors*/ - /*FUNCTION Quadtree::Quadtree(){{{1*/ - Quadtree::Quadtree(){ - - /*Initialize fields*/ - this->MaxDepth=30; - this->NbQuadtreeBox=0; - this->NbObs=0; - - /*Create container*/ - this->boxcontainer=new DataSet(); - - /*Create Root, pointer toward the main box*/ - this->root=NewQuadtreeBox(); - - } - /*}}}1*/ +/*FUNCTION Quadtree::Quadtree(){{{1*/ +Quadtree::Quadtree(){ + _error_("Constructor not supported"); + +} +/*}}}1*/ +/*FUNCTION Quadtree::Quadtree(double xmin,double xmax,double ymin,double ymax,int maxdepth){{{1*/ +Quadtree::Quadtree(double xmin,double xmax,double ymin,double ymax,int maxdepth){ + + /*Intermediaries*/ + double length; + + /*Initialize fields*/ + this->MaxDepth=maxdepth; + this->NbQuadtreeBox=0; + this->NbObs=0; + + /*Create container*/ + this->boxcontainer=new DataSet(); + + /*Create Root, pointer toward the main box*/ + length=max(xmax-xmin,ymax-ymin); + this->root=NewQuadtreeBox(xmin+length/2,ymin+length/2,length); +} +/*}}}1*/ /*FUNCTION Quadtree::~Quadtree(){{{1*/ Quadtree::~Quadtree(){ @@ -110,8 +119,9 @@ /*Intermediaries*/ - int xi,yi,ij,level; - QuadtreeBox **pb = NULL; - QuadtreeBox *b = NULL; - QuadtreeBox *bb = NULL; + int xi,yi,ij,level,levelbin; + QuadtreeBox **pbox = NULL; // pointer toward current box b + QuadtreeBox **pmaster = NULL; // pointer toward master of b + QuadtreeBox *box = NULL; // current box b + QuadtreeBox *slave = NULL; // suslaveox of b (if necessary) Observation *obs[4]; @@ -120,41 +130,42 @@ yi = observation->yi; - /*Initialize level*/ - level=(1L<MaxDepth);// = 2^30 + /*Initialize levels*/ + level = 0; + levelbin = (1L<MaxDepth);// = 2^30 /*Get inital box (the largest)*/ - pb=&root; + pmaster = &root; + pbox = &root; /*Find the smallest box where the observation is located*/ - while((b=*pb) && (b->nbitems<0)){ - - /*Go down one level = 00100 -> 00010*/ - level>>=1; - - /*Get next subbox according to the bit value (level)*/ - pb = &b->box[IJ(xi,yi,level)]; - } - _assert_(level>0); + while((box=*pbox) && (box->nbitems<0)){ + + /*Go down one level (levelbin = 00100 -> 00010)*/ + levelbin>>=1; level+=1; _assert_(levelMaxDepth); + + /*Get next box according to the bit value (levelbin)*/ + pmaster = pbox; + pbox = &box->box[IJ(xi,yi,levelbin)]; + } + _assert_(levelbin>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)){ + while((box=*pbox) && (box->nbitems==4)){ /*Copy the 4 observation in the current Quadtreebox*/ - 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; + obs[0] = box->obs[0]; + obs[1] = box->obs[1]; + obs[2] = box->obs[2]; + obs[3] = box->obs[3]; + + /*set nbitems as -1 (now holding boxes instead of observations)*/ + box->nbitems = -1; + box->box[0] = NULL; + box->box[1] = NULL; + box->box[2] = NULL; + box->box[3] = NULL; + + /*Go down one level (levelbin = 00010 -> 00001)*/ + levelbin>>=1; level+=1; _assert_(levelMaxDepth); /*Put the four observations in the new boxes*/ @@ -162,24 +173,26 @@ /*Get box for observation number k*/ - ij = IJ(obs[k]->xi,obs[k]->yi,level); - bb = b->box[ij]; - if(!bb){ - b->box[ij]=NewQuadtreeBox(); - bb=b->box[ij]; + ij = IJ(obs[k]->xi,obs[k]->yi,levelbin); + slave = box->box[ij]; + if(!slave){ + box->box[ij] = NewQuadtreeBox(box,ij); + slave = box->box[ij]; } - bb->obs[bb->nbitems++] = obs[k]; + slave->obs[slave->nbitems++] = obs[k]; } - /*Get the subbox where the current observation is located*/ - ij=IJ(xi,yi,level); - pb=&b->box[ij]; + /*Get the suslaveox where the current observation is located*/ + ij = IJ(xi,yi,levelbin); + pmaster = pbox; + pbox = &box->box[ij]; } /*alloc the QuadtreeBox if necessary and add current observation*/ - b=*pb; - if(!b){ - b=*pb=NewQuadtreeBox(); - } - b->obs[b->nbitems++]=observation; + box = *pbox; + if(!box){ + ij = IJ(xi,yi,levelbin); + box = *pbox = NewQuadtreeBox(*pmaster,ij); + } + box->obs[box->nbitems++]=observation; NbObs++; @@ -189,15 +202,77 @@ printf("Quadtree:\n"); - printf(" MaxDepth = %i\n",MaxDepth); - printf(" NbQuadtreeBox = %i\n",NbQuadtreeBox); - printf(" NbObs = %i\n",NbObs); - printf(" root = %p\n",root); - -}/*}}}*/ -/*FUNCTION Quadtree::NewQuadtreeBox {{{1*/ -Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){ + printf(" MaxDepth = %i\n",this->MaxDepth); + printf(" NbQuadtreeBox = %i\n",this->NbQuadtreeBox); + printf(" NbObs = %i\n",this->NbObs); + printf(" root = %p\n",this->root); + +}/*}}}*/ +/*FUNCTION Quadtree::DeepEcho{{{1*/ +void Quadtree::DeepEcho(void){ + + printf("Quadtree:\n"); + printf(" MaxDepth = %i\n",this->MaxDepth); + printf(" NbQuadtreeBox = %i\n",this->NbQuadtreeBox); + printf(" NbObs = %i\n",this->NbObs); + printf(" root = %p\n",this->root); + boxcontainer->Echo(); + +}/*}}}*/ +/*FUNCTION Quadtree::IntergerCoordinates{{{1*/ +void Quadtree::IntergerCoordinates(int *xi,int *yi,double x,double y){ + + /*Intermediaries*/ + double coefficient; + double xmin,ymin; + + /*Checks in debugging mode*/ + _assert_(xi && yi); + _assert_(this->root); + + /*coeffIcoor is the coefficient used for integer coordinates: + * (x-xmin) + * xi = (2^30 -1) --------- + * length + * coefficient = (2^30 -1)/length + */ + coefficient = double((1L<MaxDepth)-1)/(this->root->length); + xmin = this->root->xcenter - this->root->length/2; + ymin = this->root->ycenter - this->root->length/2; + + *xi=int(coefficient*(x - xmin)); + *yi=int(coefficient*(y - ymin)); +}/*}}}*/ +/*FUNCTION Quadtree::NewQuadtreeBox(double xcenter,double ycenter,double length){{{1*/ +Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(double xcenter,double ycenter,double length){ /*Output*/ QuadtreeBox* newbox=NULL; + + /*Create and initialize a new box*/ + newbox=new QuadtreeBox(); + newbox->nbitems=0; + newbox->xcenter=xcenter; + newbox->ycenter=ycenter; + newbox->length=length; + newbox->box[0]=NULL; + newbox->box[1]=NULL; + newbox->box[2]=NULL; + newbox->box[3]=NULL; + + /*Add to container*/ + this->boxcontainer->AddObject(newbox); + NbQuadtreeBox++; + + /*currentbox now points toward next quadtree box*/ + return newbox; +}/*}}}*/ +/*FUNCTION Quadtree::NewQuadtreeBox(QuadtreeBox* master,int index) {{{1*/ +Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(QuadtreeBox* master,int index){ + + /*Output*/ + QuadtreeBox* newbox=NULL; + + /*Checks in debugging mode*/ + _assert_(master); /*Create and initialize a new box*/ @@ -208,4 +283,25 @@ newbox->box[2]=NULL; newbox->box[3]=NULL; + switch(index){ + case 0: + newbox->xcenter=master->xcenter - master->length/4; + newbox->ycenter=master->ycenter - master->length/4; + break; + case 1: + newbox->xcenter=master->xcenter + master->length/4; + newbox->ycenter=master->ycenter - master->length/4; + break; + case 2: + newbox->xcenter=master->xcenter - master->length/4; + newbox->ycenter=master->ycenter + master->length/4; + break; + case 3: + newbox->xcenter=master->xcenter + master->length/4; + newbox->ycenter=master->ycenter + master->length/4; + break; + default: + _error_("Case %i not supported",index); + } + newbox->length=master->length/2; /*Add to container*/ @@ -217,28 +313,41 @@ }/*}}}*/ /*FUNCTION Quadtree::QuadtreeColoring{{{1*/ -void Quadtree::QuadtreeColoring(double* A,int *xi,int *yi,int n){ - - QuadtreeBox **pb = NULL; - QuadtreeBox *b = NULL; - int level; - - for(int i=0;iMaxDepth);// = 2^30 - - /*Get inital box (the largest)*/ - pb=&root; - - /*Find the smallest box where the observation is located*/ - while((b=*pb) && (b->nbitems<0)){ - - /*Color matrix onces more*/ - A[i]+=1; - - /*Go to one box deeper*/ - level>>=1; - pb = &b->box[IJ(xi[i],yi[i],level)]; - } - } -}/*}}}*/ +void Quadtree::QuadtreeColoring(double* A,int xi,int yi){ + + QuadtreeBox **pbox = NULL; + QuadtreeBox *box = NULL; + int level,levelbin; + + /*Initialize levels*/ + level = 0; + levelbin = (1L<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; _assert_(levelMaxDepth); + + pbox = &box->box[IJ(xi,yi,levelbin)]; + } + if(box && box->nbitems>0){ + /*This box is not empty, add one level*/ + level+=1; + } + + *A=level; +}/*}}}*/ + +/*QuadtreeBox methos*/ +/*FUNCTION QuadtreeBox::Echo{{{1*/ +void Quadtree::QuadtreeBox::Echo(void){ + + printf("QuadtreeBox:\n"); + printf(" nbitems = %i\n",this->nbitems); + printf(" xcenter = %g\n",this->xcenter); + printf(" ycenter = %g\n",this->ycenter); + printf(" length = %g\n",this->length); + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12224) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12225) @@ -16,5 +16,8 @@ class QuadtreeBox: public Object{ public: - int nbitems; // number of current vertices in the box + int nbitems; // number of current vertices in the box + double xcenter; // x position of the center (double) + double ycenter; // x position of the center (double) + double length; // width of the box union{ QuadtreeBox *box[4]; @@ -23,5 +26,5 @@ /*Object functions (Needed because the Quadtree uses a Container*/ - void Echo() {_error_("not implemented yet"); }; + void Echo(); void DeepEcho() {_error_("not implemented yet"); }; int Id() {_error_("not implemented yet"); }; @@ -35,15 +38,20 @@ public: - int MaxDepth; // maximum number of subdivision - QuadtreeBox *root; // main box - long NbQuadtreeBox; // total number of boxes - long NbObs; // number of points + int MaxDepth; // maximum number of subdivision + double coordconversion; // Coefficient to convert coordinates to integer + QuadtreeBox *root; // main box + long NbQuadtreeBox; // total number of boxes + long NbObs; // number of points Quadtree(); + Quadtree(double xmin,double xmax,double ymin,double ymax,int maxdepth_in); ~Quadtree(); + void Add(Observation *observation); + void DeepEcho(void); void Echo(void); - void Add(Observation *observation); - void QuadtreeColoring(double *A,int *xi,int *yi,int n); - QuadtreeBox *NewQuadtreeBox(void); + void IntergerCoordinates(int *xi,int *yi,double x,double y); + QuadtreeBox *NewQuadtreeBox(double xcenter,double ycenter,double length); + QuadtreeBox *NewQuadtreeBox(QuadtreeBox* master,int index); + void QuadtreeColoring(double *A,int xi,int yi); }; #endif //_QUADTREEK_H