Index: /issm/trunk-jpl/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 12212)
+++ /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 12213)
@@ -28,4 +28,5 @@
 /*FUNCTION Observations::Observations(){{{*/
 Observations::Observations(){
+	this->quadtree = NULL;
 	return;
 }
@@ -34,13 +35,46 @@
 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       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];
+	for(i=1;i<n;i++){
+		xmin=min(xmin,x[i]); ymin=min(ymin,y[i]);
+		xmax=max(xmax,x[i]); ymax=max(ymax,y[i]);
+	}
+	offset=0.05*(xmax-xmin); xmin-=offset; xmax+=offset;
+	offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
+
+	/*coeffIcoor is the coefficient used for integer coordinates:
+	 *                (x-xmin)
+	 * xi = (2^30 -1) ------------ 
+	 *                   D
+	 * coefficient = (2^30 -1)/D
+	 */
+	this->coefficient = double(MaxICoor)/max(xmax-xmin,ymax-ymin); _assert_(coefficient>=0);
+
 	/*Add observations one by one*/
-	for(int i=0;i<n;i++){
-		this->AddObject(new Observation(x[i],y[i],0,0,observations_list[i]));
+	for(i=0;i<n;i++){
+		xi=int(coefficient*(x[i]-xmin));
+		yi=int(coefficient*(y[i]-ymin));
+		observation = new Observation(x[i],y[i],xi,yi,observations_list[i]);
+		this->quadtree->Add(observation);
+		this->AddObject(observation);
 	}
-
 }
 /*}}}*/
 /*FUNCTION Observations::~Observations(){{{*/
 Observations::~Observations(){
+	delete quadtree;
 	return;
 }
@@ -49,5 +83,5 @@
 /*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){
 
 	/*Output and Intermediaries*/
@@ -78,2 +112,21 @@
 	*pnobs=nobs;
 }/*}}}*/
+/*FUNCTION Observations::QuadtreeColoring{{{*/
+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));
+	for(int i=0;i<n;i++){
+		xi[i]=int(coefficient*(x[i]-xmin));
+		yi[i]=int(coefficient*(y[i]-ymin));
+	}
+
+	/*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 12212)
+++ /issm/trunk-jpl/src/c/Container/Observations.h	(revision 12213)
@@ -7,6 +7,15 @@
 
 class Obsevration;
+class Quadtree;
 
 class Observations: public DataSet{
+
+	private:
+		Quadtree* quadtree;
+		double    coefficient; //For integer coordinate conversion
+		double    xmin;
+		double    ymin;
+		double    xmax;
+		double    ymax;
 
 	public:
@@ -19,4 +28,5 @@
 		/*Methods*/
 		void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp);
+		void QuadtreeColoring(double* A,double *x,double *y,int n);
 
 };
Index: /issm/trunk-jpl/src/c/Container/Vertices.h
===================================================================
--- /issm/trunk-jpl/src/c/Container/Vertices.h	(revision 12212)
+++ /issm/trunk-jpl/src/c/Container/Vertices.h	(revision 12213)
@@ -32,6 +32,3 @@
 };
 
-
-
 #endif //ifndef _VERTICES_H_
-
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12212)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12213)
@@ -44,4 +44,9 @@
 	/*Allocation output*/
 	predictions =(double*)xmalloc(n_interp*sizeof(double));
+
+	for(i=0;i<n_interp;i++) predictions[i]=0;
+	observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
+	*ppredictions=predictions;
+	return 1;
 
 	/*Loop over all interpolations*/
@@ -101,4 +106,5 @@
 	delete observations;
 	*ppredictions=predictions;
+	return 1;
 }
 
Index: /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp	(revision 12212)
+++ /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp	(revision 12213)
@@ -3,4 +3,5 @@
  */
 
+#include <stdlib.h>
 #include "../objects.h"
 
@@ -31,9 +32,12 @@
 /*FUNCTION Observation::Echo {{{1*/
 void Observation::Echo(void){
+
+	int  bit;
+
 	printf("Observation\n");
 	printf("   x     : %g\n",this->x);
 	printf("   y     : %g\n",this->y);
-	printf("   xi    : %i\n",this->xi);
-	printf("   yi    : %i\n",this->yi);
+	printf("   xi    : "); printbinary(this->xi); printf("\n");
+	printf("   yi    : "); printbinary(this->yi); printf("\n");
 	printf("   value : %g\n",this->value);
 }
Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp	(revision 12212)
+++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp	(revision 12213)
@@ -110,12 +110,13 @@
 	/*Intermediaries*/
 	const int     MaxDeep  = 30;
-	int           k,ij;
-	long          i,j,level;
+	int          xi,yi,ij,level;
 	QuadtreeBox **pb = NULL;
 	QuadtreeBox  *b  = NULL;
+	QuadtreeBox  *bb = NULL;
+	Observation  *obs[4];
 
 	/*Get integer coodinates*/
-	i = observation->xi;
-	j = observation->yi;
+	xi = observation->xi;
+	yi = observation->yi;
 
 	/*Initialize level*/
@@ -128,12 +129,9 @@
 	while((b=*pb) && (b->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)];
+		pb = &b->box[IJ(xi,yi,level)];
 	}
 	_assert_(level>0);
@@ -143,5 +141,4 @@
 
 		/*Copy the 4 observation in the current Quadtreebox*/
-		Observation* obs[4];
 		obs[0] = b->obs[0];
 		obs[1] = b->obs[1];
@@ -162,20 +159,18 @@
 
 		/*Put the four observations in the new boxes*/
-		for (k=0;k<4;k++){
+		for (int 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];
+			ij = IJ(obs[k]->xi,obs[k]->yi,level);
+			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);
+		ij=IJ(xi,yi,level);
 		pb=&b->box[ij];
 	}
@@ -183,30 +178,67 @@
 	/*alloc the QuadtreeBox if necessary and add current observation*/
 	b=*pb;
-	if(!b) b=NewQuadtreeBox();
+	if(!b){
+		b=*pb=NewQuadtreeBox();
+	}
 	b->obs[b->nbitems++]=observation;
 	NbObs++;
 
 }/*}}}*/
-	/*FUNCTION Quadtree::NewQuadtreeBox {{{1*/
-	Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){
-
-		/*Output*/
-		QuadtreeBox* newbox=NULL;
-
-		/*Create and initialize a new box*/
-		newbox=new QuadtreeBox();
-		newbox->nbitems=0;
-		newbox->box[0]=NULL;
-		newbox->box[1]=NULL;
-		newbox->box[2]=NULL;
-		newbox->box[3]=NULL;
-
-		/*Add root to the container*/
-		boxcontainer->AddObject(newbox);
-
-		/*Increase counter*/
-		NbQuadtreeBox++;
-
-		/*currentbox now points toward next quadtree box*/
-		return newbox;
-	}/*}}}*/
+/*FUNCTION Quadtree::Echo{{{1*/
+void  Quadtree::Echo(void){
+
+	printf("Quadtree:\n");
+	printf("   NbQuadtreeBox = %i\n",NbQuadtreeBox);
+	printf("   NbObs         = %i\n",NbObs);
+	printf("   root          = %p\n",root);
+
+}/*}}}*/
+/*FUNCTION Quadtree::NewQuadtreeBox {{{1*/
+Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){
+
+	/*Output*/
+	QuadtreeBox* newbox=NULL;
+
+	/*Create and initialize a new box*/
+	newbox=new QuadtreeBox();
+	newbox->nbitems=0;
+	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::QuadtreeColoring{{{1*/
+void Quadtree::QuadtreeColoring(double* A,int *xi,int *yi,int n){
+
+	const int     MaxDeep  = 30;
+	QuadtreeBox **pb = NULL;
+	QuadtreeBox  *b  = NULL;
+	int          level;
+
+	for(int i=0;i<n;i++){
+
+		/*Initialize level*/
+		level=(1L<<MaxDeep);// = 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)];
+		}
+	}
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h	(revision 12212)
+++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h	(revision 12213)
@@ -41,6 +41,8 @@
 		Quadtree();
 		~Quadtree();
-		void Add(Observation* observation);
-		QuadtreeBox* NewQuadtreeBox(void);
+		void         Echo(void);
+		void         Add(Observation *observation);
+		void         QuadtreeColoring(double *A,int *xi,int *yi,int n);
+		QuadtreeBox *NewQuadtreeBox(void);
 };
 #endif //_QUADTREEK_H
Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/elements.h	(revision 12212)
+++ /issm/trunk-jpl/src/c/shared/Elements/elements.h	(revision 12213)
@@ -47,4 +47,16 @@
 	printf("\n");
 }
+inline void printbinary(int n) {
+	unsigned int i;
+	i=1<<(sizeof(n)*8-1);
+
+	while (i>0) {
+		if (n&i)
+		 printf("1");
+		else
+		 printf("0");
+		i>>=1;
+	}
+}
 
 #endif //ifndef _SHARED_ELEMENTS_H_