Changeset 1168

Timestamp:

06/29/09 17:01:17 (16 years ago)

Author:

Mathieu Morlighem

Message:

Brand new DataInterp with a nice linear interpolation

File:

• issm/trunk/src/c/DataInterpx/DataInterpx.cpp (modified) (3 diffs)

issm/trunk/src/c/DataInterpx/DataInterpx.cpp

@@ -11 +11 @@
 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 +37 @@
         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 +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{