Changeset 7447

Timestamp:

02/15/11 11:03:15 (14 years ago)

Author:

Mathieu Morlighem

Message:

Added Bilinear interpolation by default. MUCH SMOOTHER

Location:

issm/trunk/src/c

Files:

• EnumDefinitions/EnumDefinitions.h (modified) (2 diffs)
• EnumDefinitions/EnumToString.cpp (modified) (1 diff)
• EnumDefinitions/StringToEnum.cpp (modified) (1 diff)
• Makefile.am (modified) (2 diffs)
• modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (modified) (5 diffs)
• modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h (modified) (1 diff)
• modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp (deleted)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -8 +8 @@
 
 enum definitions{
-
         /*Datasets {{{1*/
         ConstraintsEnum,
@@ -411 +410 @@
         SsetEnum,
         GroundingLineMigrationEnum,
-        YtsEnum
+        YtsEnum,
+        /*}}}*/
+        /*Interpolation {{{1*/
+        TriangleInterpEnum,
+        BilinearInterpEnum,
+        NearestInterpEnum
         /*}}}*/
 };

issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

@@ -366 +366 @@
                 case GroundingLineMigrationEnum : return "GroundingLineMigration";
                 case YtsEnum : return "Yts";
+                case TriangleInterpEnum : return "TriangleInterp";
+                case BilinearInterpEnum : return "BilinearInterp";
+                case NearestInterpEnum : return "NearestInterp";
                 default : return "unknown";
 

issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

@@ -364 +364 @@
         else if (strcmp(name,"GroundingLineMigration")==0) return GroundingLineMigrationEnum;
         else if (strcmp(name,"Yts")==0) return YtsEnum;
+        else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
+        else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
+        else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
         else _error_("Enum %s not found",name);
 

issm/trunk/src/c/Makefile.am

@@ -505 +505 @@
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
-                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
                                         ./modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp\
@@ -1084 +1083 @@
                                         ./modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp\
                                         ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\
-                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
                                         ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\

issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp

@@ -3 +3 @@
  */ 
 
+/*Include {{{*/
 #ifdef HAVE_CONFIG_H
         #include "config.h"
@@ -12 +13 @@
 #include "../../shared/shared.h"
 #include "../../include/include.h"
-
-int InterpFromGridToMeshx( 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,double default_value) {
-
+/*}}}*/
+
+/*InterpFromGridToMeshx{{{*/
+int InterpFromGridToMeshx( 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,double default_value, int interpenum){
 
         /*output: */
@@ -28 +30 @@
         InterpFromGridToMeshxThreadStruct gate;
         int num=1;
-
         #ifdef _MULTITHREADING_
         num=_NUMTHREADS_;
@@ -82 +83 @@
         gate.M=M;
         gate.N=N;
+        gate.interp=interpenum;
 
         /*launch the thread manager with InterpFromGridToMeshxt as a core: */
@@ -89 +91 @@
         *pdata_mesh=data_mesh;
 }
+/*}}}*/
+/*InterpFromGridToMeshxt {{{1*/
+void* InterpFromGridToMeshxt(void* vpthread_handle){
+
+        /*gate variables :*/
+        InterpFromGridToMeshxThreadStruct *gate    = NULL;
+        pthread_handle                    *handle  = NULL;
+        int my_thread;
+        int num_threads;
+        int i0,i1;
+
+        /*intermediary: */
+        int    i,m,n;
+        double x_grid;
+        double y_grid;
+        double data_value;
+        double x1,x2,y1,y2;
+        double Q11,Q12,Q21,Q22;
+
+        /*recover handle and gate: */
+        handle=(pthread_handle*)vpthread_handle;
+        gate=(InterpFromGridToMeshxThreadStruct*)handle->gate;
+        my_thread=handle->id;
+        num_threads=handle->num;
+
+        /*recover parameters :*/
+        double *x_mesh        = gate->x_mesh;
+        double *y_mesh        = gate->y_mesh;
+        int     x_rows        = gate->x_rows;
+        int     y_rows        = gate->y_rows;
+        double *x             = gate->x;
+        double *y             = gate->y;
+        int     nods          = gate->nods;
+        Vec     data_mesh     = gate->data_mesh;
+        double *data          = gate->data;
+        double  default_value = gate->default_value;
+        int     interpenum    = gate->interp;
+        int     M             = gate->M;
+        int     N             = gate->N;
+
+        /*partition loop across threads: */
+        PartitionRange(&i0,&i1,nods,num_threads,my_thread);
+        for (i=i0;i<i1;i++) {
+
+                x_grid=*(x_mesh+i);
+                y_grid=*(y_mesh+i);
+
+                /*Find indices m and n into y and x, for which  y(m)<=y_grids<=y(m+1) and x(n)<=x_grid<=x(n+1)*/
+                if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
+
+                        /*    Q12             Q22
+                         * y2 x---------+-----x
+                         *    |         |     |
+                         *    |         |P    |
+                         *    |---------+-----|
+                         *    |         |     |
+                         *    |         |     |
+                         * y1 x---------+-----x Q21
+                         *    x1                 x2       
+                         *
+                         */
+                        x1=x[n]; x2=x[n+1];
+                        y1=y[m]; y2=y[m+1];
+                        Q11=data[m*N+n];
+                        Q12=data[(m+1)*N+n];
+                        Q21=data[m*N+n+1];
+                        Q22=data[(m+1)*N+n+1];
+
+                        switch(interpenum){
+                                case TriangleInterpEnum:
+                                        data_value=triangleinterp(x1,x2,y1,y2,Q11,Q12,Q21,Q22,x_grid,y_grid);
+                                        break;
+                                case BilinearInterpEnum:
+                                        data_value=bilinearinterp(x1,x2,y1,y2,Q11,Q12,Q21,Q22,x_grid,y_grid);
+                                        break;
+                                case NearestInterpEnum:
+                                        data_value=nearestinterp(x1,x2,y1,y2, Q11,Q12,Q21,Q22,x_grid,y_grid);
+                                        break;
+                                default:
+                                        printf("Interpolation %s not supported yet\n",EnumToString(interpenum));
+                                        return NULL; /*WARNING: no error because it would blow up the multithreading!*/
+                        }
+                        if(isnan(data_value)) data_value=default_value;
+                }
+                else{
+                        data_value=default_value;
+                }
+
+                VecSetValue(data_mesh,i,data_value,INSERT_VALUES);
+        }
+}/*}}}*/
+
+/*findindices {{{1*/
+bool findindices(int* pn,int* pm,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid){
+
+        bool foundx=false,foundy=false;
+        int m=-1,n=-1;
+        int i;
+
+        for (i=0;i<x_rows-1;i++){
+                if ((x[i]<=xgrid) && (xgrid<x[i+1])){
+                        n=i;
+                        foundx=true;
+                        break;
+                }
+        }
+        if(xgrid==x[x_rows-1]){
+                n=x_rows-2;
+                foundx=true;
+        }
+
+        for (i=0;i<y_rows-1;i++){
+                if ((y[i]<=ygrid) && (ygrid<y[i+1])){
+                        m=i;
+                        foundy=true;
+                        break;
+                }
+        }
+        if(ygrid==y[y_rows-1]){
+                m=y_rows-2;
+                foundy=true;
+        }
+
+        /*Assign output pointers:*/
+        *pm=m; *pn=n;
+        return (foundx && foundy);
+}/*}}}*/
+/*triangleinterp{{{1*/
+double triangleinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
+        /*split the rectangle in 2 triangle and
+         * use Lagrange P1 interpolation
+         *       
+         *   +3----+2,3' Q12----Q22
+         *   |    /|     |    /|
+         *   |   / |     |   / |
+         *   |  /  |     |  /  |
+         *   | /   |     | /   |
+         *   |/    |     |/    |
+         *   1-----2'    Q11---Q21        */
+
+        /*Intermediaries*/
+        double area,area_1,area_2,area_3;
+
+        /*Checks*/
+        _assert_(x2>x1 && y2>y1);
+        _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
+
+        /*area of the rectangle*/
+        area=(x2-x1)*(y2-y1);
+
+        /*is it the upper left triangle?*/
+        if ((x-x1)/(x2-x1)<(y-y1)/(y2-y1)){
+
+                area_1=((y2-y)*(x2-x1))/area;
+                area_2=((x-x1)*(y2-y1))/area;
+                area_3=1-area_1-area_2;
+
+                return area_1*Q11+area_2*Q22+area_3*Q12;
+        }
+        else {
+
+                area_1=((y-y1)*(x2-x1))/area;
+                area_2=((x2-x)*(y2-y1))/area;
+                area_3=1-area_1-area_2;
+
+                return area_1*Q22+area_2*Q11+area_3*Q21;
+        }
+}/*}}}*/
+/*bilinearinterp{{{1*/
+double bilinearinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
+        /*Bilinear  interpolation: (http://en.wikipedia.org/wiki/Bilinear_interpolation) */
+
+        /*    Q12    R2        Q22
+         * y2 x------x---------x
+         *    |      |         |
+         *    |      |         |
+         *    |      +P        |
+         *    |      |         |
+         *    |Q11   R1        Q21
+         * y1 x------x---------x
+         *    x1               x2
+         *
+         */
+
+        /*Checks*/
+        _assert_(x2>x1 && y2>y1);
+        _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
+
+        return 
+          +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1))
+          +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1))
+          +Q12*(x2-x)*(y-y1)/((x2-x1)*(y2-y1))
+          +Q22*(x-x1)*(y-y1)/((x2-x1)*(y2-y1));
+}
+/*}}}*/
+/*nearestinterp{{{1*/
+double nearestinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
+        /*Nearest neighbor interpolation*/
+
+        /*    Q12             Q22
+         * y2 x--------x---------x
+         *    |        |         |
+         *    |        |  xP     |
+         * ym |--------+---------|
+         *    |        |         |
+         *    |        |         |
+         * y1 x--------x---------x Q21
+         *    x1       xm        x2 
+         *
+         */
+        /*Checks*/
+        _assert_(x2>x1 && y2>y1);
+        _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
+        
+        double xm=(x2-x1)/2;
+        double ym=(y2-y1)/2;
+
+        if (x<xm && y<ym) return Q11;
+        if (x<xm && y>ym) return Q12;
+        if (x>xm && y<ym) return Q21;
+        else return Q22;
+}
+/*}}}*/

issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h

@@ -7 +7 @@
 
 #include "../../toolkits/toolkits.h"
-
-
+#include "../../EnumDefinitions/EnumDefinitions.h"
 
 /*threading: */
 typedef struct{
-
+        double* x;
+        int     x_rows;
+        double* y;
+        int     y_rows;
+        double* data;
+        double  default_value;
+        int     interp;
+        int     M;
+        int     N;
+        int     nods;
         double* x_mesh;
         double* y_mesh;
-        int     x_rows;
-        int     y_rows;
-        double* x;
-        double* y;
-        int     nods;
-        double  default_value;
         Vec     data_mesh;
-        double* data;
-        int     M;
-        int     N;
-
 } InterpFromGridToMeshxThreadStruct;
 
-int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
-
-int InterpFromGridToMeshx( 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, double default_value);
-
-void* InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
-
+int    InterpFromGridToMeshx( 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, double default_value, int interpenum=BilinearInterpEnum);
+void*  InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
+bool   findindices(int* pn,int* pm,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
+double triangleinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y);
+double bilinearinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y);
+double nearestinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y);
 
 #endif /* _INTERPFROMGRIDTOMESHX_H */
-