Changeset 7439


Ignore:
Timestamp:
02/14/11 10:48:11 (14 years ago)
Author:
Mathieu Morlighem
Message:

Added bilinear interpolation

Location:
issm/trunk/src/c/modules/InterpFromGridToMeshx
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r6412 r7439  
    3434
    3535        /*Some checks on arguments: */
    36         if ((M<=2) || (N<=2) || (nods<=0)){
     36        if ((M<2) || (N<2) || (nods<=0)){
    3737                _error_("nothing to be done according to the dimensions of input matrices and vectors.");
    3838        }

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

    r3913 r7439  
    3636        double area_1,area_2,area_3;
    3737        double xi,eta;
     38        bool   triangle=true;
    3839
    3940
    … …  
    6162        PartitionRange(&i0,&i1,nods,num_threads,my_thread);
    6263
    63         /*Linear (triangle) interpolation: */
    6464        for ( i=i0; i<i1; i++) {
    6565
    … …  
    7070                if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
    7171                       
    72                         /*Get area*/
    73                         area=(x[n+1]-x[n])*(y[m+1]-y[m]);
    74 
    75                         /*is it the upper right triangle?*/
    76                         /*2'     3'
    77                          *+-----+
    78                          *1\    |
    79                          *| \   |
    80                          *|  \  |
    81                          *|   \ |
    82                          *|    \|
    83                          *2----3+1' */
    84                         if ((x_grid-x[n])/(x[n+1]-x[n])<(y_grid-y[m])/(y[m+1]-y[m])){
    85 
    86                                 /*Then find values of data at each summit*/
    87                                 G1=*(data+m*N+n);
    88                                 G2=*(data+(m+1)*N+n+1);
    89                                 G3=*(data+(m+1)*N+n);
    90 
    91                                 /*Get first area coordinate*/
    92                                 area_1=((y[m+1]-y_grid)*(x[n+1]-x[n]))/area;
    93                                 /*Get second area coordinate*/
    94                                 area_2=((x_grid-x[n])*(y[m+1]-y[m]))/area;
    95                                 /*Get third area coordinate = 1-area1-area2*/
    96                                 area_3=1-area_1-area_2;
    97 
    98                                 /*interpolate*/
    99                                 data_value=area_1*G1+area_2*G2+area_3*G3;
     72                        if(triangle){
     73                                /*Linear (triangle) interpolation: {{{1*/
     74                                /*Get area*/
     75                                area=(x[n+1]-x[n])*(y[m+1]-y[m]);
     76
     77                                /*is it the upper right triangle?*/
     78                                /*2'     3'
     79                                 *+-----+
     80                                 *1\    |
     81                                 *| \   |
     82                                 *|  \  |
     83                                 *|   \ |
     84                                 *|    \|
     85                                 *2----3+1' */
     86                                if ((x_grid-x[n])/(x[n+1]-x[n])<(y_grid-y[m])/(y[m+1]-y[m])){
     87
     88                                        /*Then find values of data at each summit*/
     89                                        G1=*(data+m*N+n);
     90                                        G2=*(data+(m+1)*N+n+1);
     91                                        G3=*(data+(m+1)*N+n);
     92
     93                                        /*Get first area coordinate*/
     94                                        area_1=((y[m+1]-y_grid)*(x[n+1]-x[n]))/area;
     95                                        /*Get second area coordinate*/
     96                                        area_2=((x_grid-x[n])*(y[m+1]-y[m]))/area;
     97                                        /*Get third area coordinate = 1-area1-area2*/
     98                                        area_3=1-area_1-area_2;
     99
     100                                        /*interpolate*/
     101                                        data_value=area_1*G1+area_2*G2+area_3*G3;
     102                                }
     103                                else {
     104
     105                                        /*Then find values of data at each summit*/
     106                                        G1=*(data+(m+1)*N+n+1);
     107                                        G2=*(data+m*N+n);
     108                                        G3=*(data+m*N+n+1);
     109
     110                                        /*Get first area coordinate*/
     111                                        area_1=((y_grid-y[m])*(x[n+1]-x[n]))/area;
     112                                        /*Get second area coordinate*/
     113                                        area_2=((x[n+1]-x_grid)*(y[m+1]-y[m]))/area;
     114                                        /*Get third area coordinate = 1-area1-area2*/
     115                                        area_3=1-area_1-area_2;
     116
     117                                        /*interpolate*/
     118                                        data_value=area_1*G1+area_2*G2+area_3*G3;
     119                                }
     120                                /*}}}*/
    100121                        }
    101                         else {
    102 
    103                                 /*Then find values of data at each summit*/
    104                                 G1=*(data+(m+1)*N+n+1);
    105                                 G2=*(data+m*N+n);
    106                                 G3=*(data+m*N+n+1);
    107 
    108                                 /*Get first area coordinate*/
    109                                 area_1=((y_grid-y[m])*(x[n+1]-x[n]))/area;
    110                                 /*Get second area coordinate*/
    111                                 area_2=((x[n+1]-x_grid)*(y[m+1]-y[m]))/area;
    112                                 /*Get third area coordinate = 1-area1-area2*/
    113                                 area_3=1-area_1-area_2;
    114 
    115                                 /*interpolate*/
    116                                 data_value=area_1*G1+area_2*G2+area_3*G3;
     122                        else{
     123                                /*Bilinear  interpolation: (http://en.wikipedia.org/wiki/Bilinear_interpolation) {{{1*/
     124
     125                                /*    Q12    R2        Q22
     126                                 * y2 x------x---------x
     127                                 *    |      |         |
     128                                 *    |      |         |
     129                                 *    |      +P        |
     130                                 *    |      |         |
     131                                 *    |Q11   R1        Q21
     132                                 * y1 x------x---------x
     133                                 *    x1               x2
     134                                 *
     135                                 */
     136
     137                                /*Coordinates*/
     138                                double x1=x[n];
     139                                double x2=x[n+1];
     140                                double y1=y[m];
     141                                double y2=y[m+1];
     142
     143                                /*Values on each summit*/
     144                                double Q11=*(data+m    *N+n);
     145                                double Q21=*(data+m    *N+n+1);
     146                                double Q12=*(data+(m+1)*N+n);
     147                                double Q22=*(data+(m+1)*N+n+1);
     148
     149                                /*Interpolate*/
     150                                data_value=
     151                                  +Q11*(x2-x_grid)*(y2-y_grid)/((x2-x1)*(y2-y1))
     152                                  +Q21*(x_grid-x1)*(y2-y_grid)/((x2-x1)*(y2-y1))
     153                                  +Q12*(x2-x_grid)*(y_grid-y1)/((x2-x1)*(y2-y1))
     154                                  +Q22*(x_grid-x1)*(y_grid-y1)/((x2-x1)*(y2-y1));
     155                                /*}}}*/
    117156                        }
    118157
Note: See TracChangeset for help on using the changeset viewer.