remap_bicubic_scrip.c 11.1 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
#include "cdo.h"
#include "cdo_int.h"
#include "grid.h"
#include "remap.h"
5
#include "remap_store_link.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13


/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/*                                                                         */
/*      BICUBIC INTERPOLATION                                              */
/*                                                                         */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */

Uwe Schulzweida's avatar
Uwe Schulzweida committed
14
15
16
17
static
void set_bicubic_weights(double iw, double jw, double wgts[4][4])
{
  wgts[0][0] = (1.-jw*jw*(3.-2.*jw))  * (1.-iw*iw*(3.-2.*iw));
18
19
20
21
  wgts[1][0] = (1.-jw*jw*(3.-2.*jw))  *     iw*iw*(3.-2.*iw);
  wgts[2][0] =     jw*jw*(3.-2.*jw)   *     iw*iw*(3.-2.*iw);
  wgts[3][0] =     jw*jw*(3.-2.*jw)   * (1.-iw*iw*(3.-2.*iw));
  wgts[0][1] = (1.-jw*jw*(3.-2.*jw))  *     iw*(iw-1.)*(iw-1.);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
  wgts[1][1] = (1.-jw*jw*(3.-2.*jw))  *     iw*iw*(iw-1.);
23
24
25
26
  wgts[2][1] =     jw*jw*(3.-2.*jw)   *     iw*iw*(iw-1.);
  wgts[3][1] =     jw*jw*(3.-2.*jw)   *     iw*(iw-1.)*(iw-1.);
  wgts[0][2] =     jw*(jw-1.)*(jw-1.) * (1.-iw*iw*(3.-2.*iw));
  wgts[1][2] =     jw*(jw-1.)*(jw-1.) *     iw*iw*(3.-2.*iw);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
  wgts[2][2] =     jw*jw*(jw-1.)      *     iw*iw*(3.-2.*iw);
28
29
30
31
  wgts[3][2] =     jw*jw*(jw-1.)      * (1.-iw*iw*(3.-2.*iw));
  wgts[0][3] =     iw*(iw-1.)*(iw-1.) *     jw*(jw-1.)*(jw-1.);
  wgts[1][3] =     iw*iw*(iw-1.)      *     jw*(jw-1.)*(jw-1.);
  wgts[2][3] =     iw*iw*(iw-1.)      *     jw*jw*(jw-1.);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
33
34
  wgts[3][3] =     iw*(iw-1.)*(iw-1.) *     jw*jw*(jw-1.);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
int num_src_points(const int* restrict mask, const int src_add[4], double src_lats[4]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
void renormalize_weights(const double src_lats[4], double wgts[4][4])
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
41
42
{
  int n;
  double sum_wgts = 0.0; /* sum of weights for normalization */
  /* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
  for ( n = 0; n < 4; ++n ) sum_wgts  += fabs(src_lats[n]);
44
45
46
47
  for ( n = 0; n < 4; ++n ) wgts[n][0] = fabs(src_lats[n])/sum_wgts;
  for ( n = 0; n < 4; ++n ) wgts[n][1] = 0.;
  for ( n = 0; n < 4; ++n ) wgts[n][2] = 0.;
  for ( n = 0; n < 4; ++n ) wgts[n][3] = 0.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
}

static
void bicubic_warning(void)
{
  static int lwarn = TRUE;

  if ( cdoVerbose || lwarn )
    {
      lwarn = FALSE;
      // cdoWarning("Iteration for iw,jw exceed max iteration count of %d!", remap_max_iter);
      cdoWarning("Bicubic interpolation failed for some grid points - used a distance-weighted average instead!");
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
64
static
void bicubic_remap(double* restrict tgt_point, const double* restrict src_array, double wgts[4][4], const int src_add[4],
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
66
67
		   const double* restrict grad1, const double* restrict grad2, const double* restrict grad3)
{
  *tgt_point = 0.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
  for ( int n = 0; n < 4; ++n )
69
70
71
72
    *tgt_point += src_array[src_add[n]]*wgts[n][0] +
                      grad1[src_add[n]]*wgts[n][1] +
                      grad2[src_add[n]]*wgts[n][2] +
                      grad3[src_add[n]]*wgts[n][3];  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
74
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
78
79
80
81
/*
  -----------------------------------------------------------------------

  This routine computes the weights for a bicubic interpolation.

  -----------------------------------------------------------------------
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
void scrip_remap_bicubic_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
84
85
{
  /*   Local variables */
  int  search_result;
86
  long tgt_cell_add;        /*  destination addresss                   */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
88
89
90
91
  int src_add[4];      /*  address for the four source points     */
  double src_lats[4];  /*  latitudes  of four bilinear corners    */
  double src_lons[4];  /*  longitudes of four bilinear corners    */
  double wgts[4][4];   /*  bicubic weights for four corners       */
  double plat, plon;   /*  lat/lon coords of destination point    */
92
  extern int timer_remap_bic;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
94
95
96
  int remap_grid_type = src_grid->remap_grid_type;

  if ( cdoVerbose ) cdoPrint("Called %s()", __func__);

97
  if ( cdoTimer ) timer_start(timer_remap_bic);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98

99
  progressInit();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
102
103
104
105

  /* Compute mappings from source to target grid */

  if ( src_grid->rank != 2 )
    cdoAbort("Can not do bicubic interpolation when source grid rank != 2"); 

106
107
108
109
  long tgt_grid_size = tgt_grid->size;

  weightlinks4_t *weightlinks = (weightlinks4_t *) malloc(tgt_grid_size*sizeof(weightlinks4_t));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
111
  /* Loop over destination grid */

112
113
  double findex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
116
  shared(ompNumThreads, cdoVerbose, weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
117
  private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
#endif
119
  for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
121
    {
      int lprogress = 1;
122
123
      if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
#if defined(_OPENMP)
125
#include "pragma_omp_atomic_update.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
127
128
129
#endif
      findex++;
      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);

130
131
132
      weightlinks[tgt_cell_add].nlinks = 0;	

      if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133

134
135
      plat = tgt_grid->cell_center_lat[tgt_cell_add];
      plon = tgt_grid->cell_center_lon[tgt_cell_add];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

      /* Find nearest square of grid points on source grid  */
      if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
	search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons, 
					  plat, plon, src_grid->dims,
					  src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
      else
	search_result = grid_search(src_grid, src_add, src_lats, src_lons, 
				    plat, plon, src_grid->dims,
				    src_grid->cell_center_lat, src_grid->cell_center_lon,
				    src_grid->cell_bound_box, src_grid->bin_addr);

      /* Check to see if points are land points */
      if ( search_result > 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
	  for ( int n = 0; n < 4; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
153
154
155
156
157
158
159
	    if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
	}

      /* If point found, find local iw,jw coordinates for weights  */
      if ( search_result > 0 )
	{
	  double iw, jw;  /*  current guess for bilinear coordinate  */

160
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161

162
          if ( find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
164
	    {
	      /* Successfully found iw,jw - compute weights */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
	      set_bicubic_weights(iw, jw, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166

167
	      store_weightlinks4(4, src_add, wgts, tgt_cell_add, weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
169
170
	    }
          else
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171
	      bicubic_warning();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174
175
176
177
178

	      search_result = -1;
	    }
	}
	  
      /*
	Search for bicubic failed - use a distance-weighted average instead (this is typically near the pole)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
	Distance was stored in src_lats!
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
181
182
      */
      if ( search_result < 0 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
	  if ( num_src_points(src_grid->mask, src_add, src_lats) > 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
	      renormalize_weights(src_lats, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186

187
	      tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188

189
	      store_weightlinks4(4, src_add, wgts, tgt_cell_add, weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
191
	    }
        }
192
193
194
195
196
197
198
    }

  if ( cdoTimer ) timer_stop(timer_remap_bic);

  weightlinks2remaplinks4(tgt_grid_size, weightlinks, rv);

  if ( weightlinks ) free(weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199

Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
201
202
203
204
} /* scrip_remap_weights_bicubic */

/*
  -----------------------------------------------------------------------

Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
  This routine computes ans apply the weights for a bicubic interpolation.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
207
208
209
210
211
212

  -----------------------------------------------------------------------
*/
void scrip_remap_bicubic(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const double* restrict src_array, double* restrict tgt_array, double missval)
{
  /*   Local variables */
  int  search_result;
213
  long tgt_cell_add;        /*  destination addresss                 */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
215
216
217
218
  int src_add[4];      /*  address for the four source points   */
  double src_lats[4];  /*  latitudes  of four bilinear corners  */
  double src_lons[4];  /*  longitudes of four bilinear corners  */
  double wgts[4][4];   /*  bicubic weights for four corners     */
  double plat, plon;   /*  lat/lon coords of destination point  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
221
222
223
224
  int remap_grid_type = src_grid->remap_grid_type;

  if ( cdoVerbose ) cdoPrint("Called %s()", __func__);

  progressInit();

225
  long tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
227
228
229
230
231

  /* Compute mappings from source to target grid */

  if ( src_grid->rank != 2 )
    cdoAbort("Can not do bicubic interpolation when source grid rank != 2"); 

232
233
234
  double *grad1_lat    = (double*) malloc(src_grid->size*sizeof(double));
  double *grad1_lon    = (double*) malloc(src_grid->size*sizeof(double));
  double *grad1_latlon = (double*) malloc(src_grid->size*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
237

  remap_gradients(*src_grid, src_array, grad1_lat, grad1_lon, grad1_latlon);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
  /* Loop over destination grid */

240
241
  double findex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
#if defined(_OPENMP)
#pragma omp parallel for default(none) \
244
  shared(ompNumThreads, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, grad1_lat, grad1_lon, grad1_latlon, findex) \
245
  private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
#endif
247
  for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
    {
      int lprogress = 1;
250
251
      if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
#if defined(_OPENMP)
253
#include "pragma_omp_atomic_update.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
254
255
256
257
#endif
      findex++;
      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);

258
      tgt_array[tgt_cell_add] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259

260
      if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261

262
263
      plat = tgt_grid->cell_center_lat[tgt_cell_add];
      plon = tgt_grid->cell_center_lon[tgt_cell_add];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287

      /* Find nearest square of grid points on source grid  */
      if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
	search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons, 
					  plat, plon, src_grid->dims,
					  src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
      else
	search_result = grid_search(src_grid, src_add, src_lats, src_lons, 
				    plat, plon, src_grid->dims,
				    src_grid->cell_center_lat, src_grid->cell_center_lon,
				    src_grid->cell_bound_box, src_grid->bin_addr);

      /* Check to see if points are land points */
      if ( search_result > 0 )
	{
	  for ( int n = 0; n < 4; ++n )
	    if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
	}

      /* If point found, find local iw,jw coordinates for weights  */
      if ( search_result > 0 )
	{
	  double iw, jw;  /*  current guess for bilinear coordinate  */

288
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
291
292
293
294

          if ( find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw) )
	    {
	      /* Successfully found iw,jw - compute weights */
	      set_bicubic_weights(iw, jw, wgts);

295
	      sort_add_and_wgts4(4, src_add, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296

297
	      bicubic_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add, grad1_lat, grad1_lon, grad1_latlon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
	    }
          else
	    {
	      bicubic_warning();

	      search_result = -1;
	    }
	}
	  
      /*
	Search for bicubic failed - use a distance-weighted average instead (this is typically near the pole)
	Distance was stored in src_lats!
      */
      if ( search_result < 0 )
	{
	  if ( num_src_points(src_grid->mask, src_add, src_lats) > 0 )
	    {
	      renormalize_weights(src_lats, wgts);

317
	      tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318

319
	      sort_add_and_wgts4(4, src_add, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320

321
	      bicubic_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add, grad1_lat, grad1_lon, grad1_latlon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322
323
	    }
        }
324
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
326
327
328
329
330

  free(grad1_lat);
  free(grad1_lon);
  free(grad1_latlon);

} /* scrip_remap_bicubic */