Ice Sheet System Model  4.18
Code documentation
Functions
InterpFromGridToMeshx.cpp File Reference

"c" core code for interpolating values from a structured grid. More...

#include <cstring>
#include "./InterpFromGridToMeshx.h"
#include "../../shared/shared.h"
#include "../../shared/io/io.h"

Go to the source code of this file.

Functions

int InterpFromGridToMeshx (IssmSeqVec< IssmPDouble > **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, const char *interptype)
 
void * InterpFromGridToMeshxt (void *vpthread_handle)
 
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)
 

Detailed Description

"c" core code for interpolating values from a structured grid.

Definition in file InterpFromGridToMeshx.cpp.

Function Documentation

◆ InterpFromGridToMeshx()

int InterpFromGridToMeshx ( IssmSeqVec< IssmPDouble > **  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,
const char *  interptype 
)

Definition at line 17 of file InterpFromGridToMeshx.cpp.

17  {
18 
19  /*output: */
20  IssmSeqVec<IssmPDouble>* data_mesh=NULL;
21 
22  /*Intermediary*/
23  double* x=NULL;
24  double* y=NULL;
25  int i;
26 
27  /*Some checks on arguments: */
28  if ((M<2) || (N<2) || (nods<=0)){
29  _error_("nothing to be done according to the dimensions of input matrices and vectors.");
30  }
31  if (x_in[1]-x_in[0]<0){
32  _error_("x coordinate vector should be increasing.\n use Matlab's command x=flipud(x), also flip the data matrix data=fliplr(data)");
33  }
34  if (y_in[1]-y_in[0]<0){
35  _error_("y coordinate vector should be increasing.\n use Matlab's command y=flipud(y), also flip the data matrix data=flipud(data)");
36  }
37 
38  /*Allocate output vector: */
39  data_mesh=new IssmSeqVec<IssmPDouble>(nods);
40 
41  /*Find out what kind of coordinates (x_in,y_in) have been given is input*/
42  if(N==(x_rows-1) && M==(y_rows-1)){
43 
44  /*The coordinates given in input describe the contour of each pixel. Take the center of each pixel*/
45  x=xNew<double>(N);
46  y=xNew<double>(M);
47  for(i=0;i<N;i++) x[i]=(x_in[i]+x_in[i+1])/2.;
48  for(i=0;i<M;i++) y[i]=(y_in[i]+y_in[i+1])/2.;
49  x_rows=x_rows-1;
50  y_rows=y_rows-1;
51  }
52  else if (N==x_rows && M==y_rows){
53 
54  /*The coordinates given in input describe the center each pixel. Keep them*/
55  x=xNew<double>(N);
56  y=xNew<double>(M);
57  for(i=0;i<N;i++) x[i]=x_in[i];
58  for(i=0;i<M;i++) y[i]=y_in[i];
59  }
60  else{
61  _error_("x and y vectors length should be 1 or 0 more than data number of rows.");
62  }
63 
64  /*initialize thread parameters: */
65  InterpFromGridToMeshxThreadStruct gate;
66  gate.x_mesh = x_mesh;
67  gate.y_mesh = y_mesh;
68  gate.x_rows = x_rows;
69  gate.y_rows = y_rows;
70  gate.x = x;
71  gate.y = y;
72  gate.nods = nods;
73  gate.data_mesh = data_mesh;
74  gate.data = data;
75  gate.default_value = default_value;
76  gate.interp = interptype;
77  gate.M = M;
78  gate.N = N;
79 
80  /*launch the thread manager with InterpFromGridToMeshxt as a core: */
81  LaunchThread(InterpFromGridToMeshxt,(void*)&gate,_NUMTHREADS_);
82  //_printf_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"% \n");
83 
84  /*Assign output pointers:*/
85  *pdata_mesh=data_mesh;
86  return 1;
87 }

◆ InterpFromGridToMeshxt()

void* InterpFromGridToMeshxt ( void *  vpthread_handle)

Definition at line 90 of file InterpFromGridToMeshx.cpp.

90  {
91 
92  /*gate variables :*/
93  InterpFromGridToMeshxThreadStruct *gate = NULL;
94  pthread_handle *handle = NULL;
95  int my_thread;
96  int num_threads;
97  int i0,i1;
98 
99  /*intermediary: */
100  int i,m,n;
101  double x_grid;
102  double y_grid;
103  double data_value;
104  double x1,x2,y1,y2;
105  double Q11,Q12,Q21,Q22;
106 
107  /*recover handle and gate: */
108  handle=(pthread_handle*)vpthread_handle;
109  gate=(InterpFromGridToMeshxThreadStruct*)handle->gate;
110  my_thread=handle->id;
111  num_threads=handle->num;
112 
113  /*recover parameters :*/
114  double *x_mesh = gate->x_mesh;
115  double *y_mesh = gate->y_mesh;
116  int x_rows = gate->x_rows;
117  int y_rows = gate->y_rows;
118  double *x = gate->x;
119  double *y = gate->y;
120  int nods = gate->nods;
121  IssmSeqVec<IssmPDouble>*data_mesh = gate->data_mesh;
122  double *data = gate->data;
123  double default_value = gate->default_value;
124  const char* interptype = gate->interp;
125  int M = gate->M;
126  int N = gate->N;
127 
128  bool debug = M*N>1? true:false;
129 
130  /*partition loop across threads: */
131  PartitionRange(&i0,&i1,nods,num_threads,my_thread);
132  for (i=i0;i<i1;i++) {
133 
134  //if(debug && my_thread==0) _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(i-i0)/double(i1-i0)*100<<"% ");
135  x_grid=*(x_mesh+i);
136  y_grid=*(y_mesh+i);
137 
138  /*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)*/
139  if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
140 
141  /* Q12 Q22
142  * y2 x---------+-----x
143  * | | |
144  * | |P |
145  * |---------+-----|
146  * | | |
147  * | | |
148  * y1 x---------+-----x Q21
149  * x1 x2
150  *
151  */
152  x1=x[n]; x2=x[n+1];
153  y1=y[m]; y2=y[m+1];
154  Q11=data[m*N+n];
155  Q12=data[(m+1)*N+n];
156  Q21=data[m*N+n+1];
157  Q22=data[(m+1)*N+n+1];
158 
159  if(strcmp(interptype,"triangle")==0){
160  data_value=triangleinterp(x1,x2,y1,y2,Q11,Q12,Q21,Q22,x_grid,y_grid);
161  }
162  else if(strcmp(interptype,"bilinear")==0){
163  data_value=bilinearinterp(x1,x2,y1,y2,Q11,Q12,Q21,Q22,x_grid,y_grid);
164  }
165  else if(strcmp(interptype,"nearest")==0){
166  data_value=nearestinterp(x1,x2,y1,y2, Q11,Q12,Q21,Q22,x_grid,y_grid);
167  }
168  else{
169  _printf_("Interpolation " << interptype << " not supported yet (supported intepolations are: triangle, bilinear and nearest)\n");
170  return NULL; /*WARNING: no error because it would blow up the multithreading!*/
171  }
172  if(xIsNan<IssmPDouble>(data_value)) data_value=default_value;
173  }
174  else{
175  data_value=default_value;
176  }
177 
178  data_mesh->SetValue(i,data_value,INS_VAL);
179  }
180 
181  return NULL;
182 }/*}}}*/

◆ findindices()

bool findindices ( int *  pn,
int *  pm,
double *  x,
int  x_rows,
double *  y,
int  y_rows,
double  xgrid,
double  ygrid 
)

Definition at line 185 of file InterpFromGridToMeshx.cpp.

185  {
186 
187  bool foundx=false,foundy=false;
188  int m=-1,n=-1;
189  int i;
190 
191  for (i=0;i<x_rows-1;i++){
192  if ((x[i]<=xgrid) && (xgrid<x[i+1])){
193  n=i;
194  foundx=true;
195  break;
196  }
197  }
198  if(xgrid==x[x_rows-1]){
199  n=x_rows-2;
200  foundx=true;
201  }
202 
203  for (i=0;i<y_rows-1;i++){
204  if ((y[i]<=ygrid) && (ygrid<y[i+1])){
205  m=i;
206  foundy=true;
207  break;
208  }
209  }
210  if(ygrid==y[y_rows-1]){
211  m=y_rows-2;
212  foundy=true;
213  }
214 
215  /*Assign output pointers:*/
216  *pm=m; *pn=n;
217  return (foundx && foundy);
218 }/*}}}*/

◆ triangleinterp()

double triangleinterp ( double  x1,
double  x2,
double  y1,
double  y2,
double  Q11,
double  Q12,
double  Q21,
double  Q22,
double  x,
double  y 
)

Definition at line 220 of file InterpFromGridToMeshx.cpp.

220  {
221  /*split the rectangle in 2 triangle and
222  * use Lagrange P1 interpolation
223  *
224  * +3----+2,3' Q12----Q22
225  * | /| | /|
226  * | / | | / |
227  * | / | | / |
228  * | / | | / |
229  * |/ | |/ |
230  * 1-----2' Q11---Q21 */
231 
232  /*Intermediaries*/
233  double area,area_1,area_2,area_3;
234 
235  /*Checks*/
236  _assert_(x2>x1 && y2>y1);
237  _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
238 
239  /*area of the rectangle*/
240  area=(x2-x1)*(y2-y1);
241 
242  /*is it the upper left triangle?*/
243  if ((x-x1)/(x2-x1)<(y-y1)/(y2-y1)){
244 
245  area_1=((y2-y)*(x2-x1))/area;
246  area_2=((x-x1)*(y2-y1))/area;
247  area_3=1-area_1-area_2;
248 
249  return area_1*Q11+area_2*Q22+area_3*Q12;
250  }
251  else {
252 
253  area_1=((y-y1)*(x2-x1))/area;
254  area_2=((x2-x)*(y2-y1))/area;
255  area_3=1-area_1-area_2;
256 
257  return area_1*Q22+area_2*Q11+area_3*Q21;
258  }
259 }/*}}}*/

◆ bilinearinterp()

double bilinearinterp ( double  x1,
double  x2,
double  y1,
double  y2,
double  Q11,
double  Q12,
double  Q21,
double  Q22,
double  x,
double  y 
)

Definition at line 261 of file InterpFromGridToMeshx.cpp.

261  {
262  /*Bilinear interpolation: (http://en.wikipedia.org/wiki/Bilinear_interpolation) */
263 
264  /* Q12 R2 Q22
265  * y2 x------x---------x
266  * | | |
267  * | | |
268  * | +P |
269  * | | |
270  * |Q11 R1 Q21
271  * y1 x------x---------x
272  * x1 x2
273  *
274  */
275 
276  /*Checks*/
277  _assert_(x2>x1 && y2>y1);
278  _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
279 
280  return
281  +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1))
282  +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1))
283  +Q12*(x2-x)*(y-y1)/((x2-x1)*(y2-y1))
284  +Q22*(x-x1)*(y-y1)/((x2-x1)*(y2-y1));
285 }

◆ nearestinterp()

double nearestinterp ( double  x1,
double  x2,
double  y1,
double  y2,
double  Q11,
double  Q12,
double  Q21,
double  Q22,
double  x,
double  y 
)

Definition at line 288 of file InterpFromGridToMeshx.cpp.

288  {
289  /*Nearest neighbor interpolation*/
290 
291  /* Q12 Q22
292  * y2 x--------x---------x
293  * | | |
294  * | | xP |
295  * ym |--------+---------|
296  * | | |
297  * | | |
298  * y1 x--------x---------x Q21
299  * x1 xm x2
300  *
301  */
302  /*Checks*/
303  _assert_(x2>x1 && y2>y1);
304  _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
305 
306  double xm=(x2-x1)/2;
307  double ym=(y2-y1)/2;
308 
309  if (x<xm && y<ym) return Q11;
310  if (x<xm && y>ym) return Q12;
311  if (x>xm && y<ym) return Q21;
312  else return Q22;
313 }
InterpFromGridToMeshxThreadStruct::M
int M
Definition: InterpFromGridToMeshx.h:21
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
InterpFromGridToMeshxThreadStruct::interp
const char * interp
Definition: InterpFromGridToMeshx.h:20
pthread_handle::id
int id
Definition: issm_threads.h:12
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
LaunchThread
void LaunchThread(void *function(void *), void *gate, int num_threads)
Definition: LaunchThread.cpp:25
pthread_handle::num
int num
Definition: issm_threads.h:13
bilinearinterp
double bilinearinterp(double x1, double x2, double y1, double y2, double Q11, double Q12, double Q21, double Q22, double x, double y)
Definition: InterpFromGridToMeshx.cpp:261
InterpFromGridToMeshxThreadStruct::N
int N
Definition: InterpFromGridToMeshx.h:22
IssmSeqVec< IssmPDouble >
InterpFromGridToMeshxThreadStruct::x_mesh
double * x_mesh
Definition: InterpFromGridToMeshx.h:24
InterpFromGridToMeshxThreadStruct::y_mesh
double * y_mesh
Definition: InterpFromGridToMeshx.h:25
InterpFromGridToMeshxThreadStruct::data_mesh
IssmSeqVec< IssmPDouble > * data_mesh
Definition: InterpFromGridToMeshx.h:26
nearestinterp
double nearestinterp(double x1, double x2, double y1, double y2, double Q11, double Q12, double Q21, double Q22, double x, double y)
Definition: InterpFromGridToMeshx.cpp:288
triangleinterp
double triangleinterp(double x1, double x2, double y1, double y2, double Q11, double Q12, double Q21, double Q22, double x, double y)
Definition: InterpFromGridToMeshx.cpp:220
InterpFromGridToMeshxThreadStruct::y
double * y
Definition: InterpFromGridToMeshx.h:16
pthread_handle
Definition: issm_threads.h:10
PartitionRange
void PartitionRange(int *pi0, int *pi1, int num_el, int num_threads, int my_thread)
Definition: PartitionRange.cpp:13
InterpFromGridToMeshxt
void * InterpFromGridToMeshxt(void *vpthread_handle)
Definition: InterpFromGridToMeshx.cpp:90
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
InterpFromGridToMeshxThreadStruct::default_value
double default_value
Definition: InterpFromGridToMeshx.h:19
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
InterpFromGridToMeshxThreadStruct::x
double * x
Definition: InterpFromGridToMeshx.h:14
InterpFromGridToMeshxThreadStruct::data
double * data
Definition: InterpFromGridToMeshx.h:18
pthread_handle::gate
void * gate
Definition: issm_threads.h:11
InterpFromGridToMeshxThreadStruct::y_rows
int y_rows
Definition: InterpFromGridToMeshx.h:17
InterpFromGridToMeshxThreadStruct
Definition: InterpFromGridToMeshx.h:13
InterpFromGridToMeshxThreadStruct::x_rows
int x_rows
Definition: InterpFromGridToMeshx.h:15
InterpFromGridToMeshxThreadStruct::nods
int nods
Definition: InterpFromGridToMeshx.h:23
findindices
bool findindices(int *pn, int *pm, double *x, int x_rows, double *y, int y_rows, double xgrid, double ygrid)
Definition: InterpFromGridToMeshx.cpp:185
Generated on Thu Jul 2 2020 08:09:23 for Ice Sheet System Model by   doxygen 1.8.19