Index: /issm/trunk/src/c/DataInterpx/DataInterpx.cpp
===================================================================
--- /issm/trunk/src/c/DataInterpx/DataInterpx.cpp	(revision 1167)
+++ /issm/trunk/src/c/DataInterpx/DataInterpx.cpp	(revision 1168)
@@ -11,27 +11,25 @@
 int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
 
-int DataInterpx( Vec* pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods) {
+int DataInterpx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods) {
 
-	int i;
-	int m,n;
-	int data_mesh_m; //local number of rows for vector data_mesh.
 
 	/*output: */
 	Vec data_mesh=NULL;
 	
+	/*Intermediary*/
+	double* x=NULL;
+	double* y=NULL;
+	int i,m,n;
+	int ind;
+	double dis;
 	double x_grid,y_grid;
 	double xi,eta;
 	double G1,G2,G3,G4,data_value;
-
+	double area_1,area_2,area_3;
+	double area;
 
 	/*Some checks on arguments: */
-	if ((M<=0) || (N<=0) || (nods<=0)){
+	if ((M<=2) || (N<=2) || (nods<=0)){
 		throw ErrorException(__FUNCT__,"nothing to be done according to the dimensions of input matrices and vectors.");
-	}
-	if(M!=(y_rows-1)){
-		throw ErrorException(__FUNCT__,"y vector length should be 1 more than data number of rows.");
-	}
-	if(N!=(x_rows-1)){
-		throw ErrorException(__FUNCT__,"x vector length should be 1 more than data number of columns.");
 	}
 
@@ -39,5 +37,28 @@
 	data_mesh=NewVec(nods);
 
-	/*Nearest interpolation method for now: */
+	/*Find out what kind of coordinates (x_in,y_in) have been given is input*/
+	if(N==(x_rows-1) && M==(y_rows-1)){
+
+		/*The coordinates given in input describe the contour of each pixel. Take the center of each pixel*/
+		x=(double*)xmalloc(N*sizeof(double));
+		y=(double*)xmalloc(M*sizeof(double));
+		for (i=0;i<N;i++) x[i]=(x_in[i]+x_in[i+1])/2;
+		for (i=0;i<M;i++) y[i]=(y_in[i]+y_in[i+1])/2;
+		x_rows=x_rows-1;
+		y_rows=y_rows-1;
+	}
+	else if (M==x_rows && N==x_rows){
+
+		/*The coordinates given in input describe the center each pixel. Keep them*/
+		x=(double*)xmalloc(N*sizeof(double));
+		y=(double*)xmalloc(M*sizeof(double));
+		for (i=0;i<N;i++) x[i]=x_in[i];
+		for (i=0;i<M;i++) y[i]=y_in[i];
+	}
+	else{
+		throw ErrorException(__FUNCT__,"x and y vectors length should be 1 or 0 more than data number of rows.");
+	}
+
+	/*Linear (triangle) interpolation: */
 	for ( i=MPI_Lowerrow(nods); i<MPI_Upperrow(nods); i++) {
 
@@ -48,15 +69,69 @@
 		if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
 			
-			/*Do a bilinear interpolation. First find parameter coordinates: */
-			xi=-1.0+2*(x_grid-x[n])/(x[n+1]-x[n]);
-			eta=-1.0+2*(y_grid-y[m])/(y[m+1]-y[m]);
-			/*Then find values of data at each grid: (xi,eta)=(0,0),(1,0),(1,1),(0,1) */
-			G1=*(data+m*N+n);
-			G2=*(data+m*N+n+1);
-			G3=*(data+(m+1)*N+n+1);
-			G4=*(data+(m+1)*N+n);
-			/*Find data_value, by bilinear interpolation: */
-			data_value=G1*(1+xi)*(1+eta)/4.0+ G2*(1-xi)*(1+eta)/4.0+G3*(1-xi)*(1-eta)/4.0+G4*(1+xi)*(1-eta)/4.0;
+			/*Get area*/
+			area=(x[n+1]-x[n])*(y[m+1]-y[m]);
 
+			/*is it the upper right triangle?*/
+			/*2'     3'
+			 *+-----+
+			 *1\    |
+			 *| \   |
+			 *|  \  |
+			 *|   \ |
+			 *|    \|
+			 *2----3+1' */
+			if ((x_grid-x[n])/(x[n+1]-x[n])<(y_grid-y[m])/(y[m+1]-y[m])){
+
+				/*Then find values of data at each summit*/
+				G1=*(data+m*N+n);
+				G2=*(data+(m+1)*N+n+1);
+				G3=*(data+(m+1)*N+n);
+
+				/*Get first area coordinate*/
+				area_1=((y[m+1]-y_grid)*(x[n+1]-x[n]))/area;
+				/*Get second area coordinate*/
+				area_2=((x_grid-x[n])*(y[m+1]-y[m]))/area;
+				/*Get third area coordinate = 1-area1-area2*/
+				area_3=1-area_1-area_2;
+
+				/*interpolate*/
+				data_value=area_1*G1+area_2*G2+area_3*G3;
+			}
+			else {
+
+				/*Then find values of data at each summit*/
+				G1=*(data+(m+1)*N+n+1);
+				G2=*(data+m*N+n);
+				G3=*(data+m*N+n+1);
+
+				/*Get first area coordinate*/
+				area_1=((y_grid-y[m])*(x[n+1]-x[n]))/area;
+				/*Get second area coordinate*/
+				area_2=((x[n+1]-x_grid)*(y[m+1]-y[m]))/area;
+				/*Get third area coordinate = 1-area1-area2*/
+				area_3=1-area_1-area_2;
+
+				/*interpolate*/
+				data_value=area_1*G1+area_2*G2+area_3*G3;
+			}
+
+			/*Treat NANs: take the closest non nan value avomg neighbors*/
+			if isnan(data_value){
+				ind=m*N+n;
+				dis=pow(x[n+1]-x[n],2);
+				if ( ((pow(y_grid-y[m],2)+pow(x_grid-x[n+1],2))<dis) && (!isnan(data[m*N+n+1]))){
+					ind=m*N+n+1;
+					dis=pow(y_grid-y[m],2)+pow(x_grid-x[n+1],2);
+				}
+				if ( ((pow(y_grid-y[m+1],2)+pow(x_grid-x[n],2))<dis) && (!isnan(data[(m+1)*N+n]))){
+					ind=(m+1)*N+n;
+					dis=pow(y_grid-y[m+1],2)+pow(x_grid-x[n],2);
+				}
+				if ( ((pow(y_grid-y[m+1],2)+pow(x_grid-x[n+1],2))<dis) && (!isnan(data[(m+1)*N+n+1]))){
+					ind=(m+1)*N+n+1;
+					dis=pow(y_grid-y[m+1],2)+pow(x_grid-x[n+1],2);
+				}
+				data_value=*(data+ind);
+			}
 		}
 		else{