remap_bicubic.cc 10.3 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/
Oliver Heidmann's avatar
Oliver Heidmann committed
17

18
#include "process_int.h"
19
#include "cdo_wtime.h"
20
#include <mpim_grid.h>
21
#include "cdo_options.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
#include "remap.h"
23
#include "remap_store_link.h"
24
#include "progress.h"
25
#include "cimdOmp.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26

Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
// bicubic interpolation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28

29
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
bicubicSetWeights(double iw, double jw, double (&weights)[4][4])
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
{
32
33
34
35
36
37
  const auto iw1 = iw*iw*(iw-1.);
  const auto iw2 = iw*(iw-1.)*(iw-1.);
  const auto iw3 = iw*iw*(3.-2.*iw);
  const auto jw1 = jw*jw*(jw-1.);
  const auto jw2 = jw*(jw-1.)*(jw-1.);
  const auto jw3 = jw*jw*(3.-2.*jw);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
  weights[0][0] = (1.-jw3)  * (1.-iw3);
  weights[1][0] = (1.-jw3)  *     iw3;
  weights[2][0] =     jw3   *     iw3;
  weights[3][0] =     jw3   * (1.-iw3);
  weights[0][1] = (1.-jw3)  *     iw2;
  weights[1][1] = (1.-jw3)  *     iw1;
  weights[2][1] =     jw3   *     iw1;
  weights[3][1] =     jw3   *     iw2;
  weights[0][2] =     jw2   * (1.-iw3);
  weights[1][2] =     jw2   *     iw3;
  weights[2][2] =     jw1   *     iw3;
  weights[3][2] =     jw1   * (1.-iw3);
  weights[0][3] =     jw2   *     iw2;
  weights[1][3] =     jw2   *     iw1;
  weights[2][3] =     jw1   *     iw1;
  weights[3][3] =     jw1   *     iw2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
}

58
template <typename T>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
int num_src_points(const Varray<T> &mask, const size_t (&src_add)[4], double (&src_lats)[4]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60

61
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
renormalizeWeights(const double (&src_lats)[4], double (&weights)[4][4])
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
66
67
68
69
  double sum_weights = 0.0;  // sum of weights for normalization
  for (unsigned n = 0; n < 4; ++n) sum_weights += std::fabs(src_lats[n]);
  for (unsigned n = 0; n < 4; ++n) weights[n][0] = std::fabs(src_lats[n]) / sum_weights;
  for (unsigned n = 0; n < 4; ++n) weights[n][1] = 0.;
  for (unsigned n = 0; n < 4; ++n) weights[n][2] = 0.;
  for (unsigned n = 0; n < 4; ++n) weights[n][3] = 0.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
71
}

72
static void
73
bicubicWarning()
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
{
75
  static bool lwarn = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76

77
  if (Options::cdoVerbose || lwarn)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
    {
79
      lwarn = false;
80
      cdoWarning("Bicubic interpolation failed for some grid points - used a distance-weighted average instead!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
82
83
    }
}

84
static double
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
bicubicRemap(const double *restrict src_array, const double (&weights)[4][4], const size_t (&src_add)[4], const RemapGradients &gradients)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
{
87
88
89
  const auto &glat = gradients.grad_lat;
  const auto &glon = gradients.grad_lon;
  const auto &glatlon = gradients.grad_latlon;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

91
  double tgt_point = 0.;
92
  for (unsigned n = 0; n < 4; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
94
    tgt_point += src_array[src_add[n]] * weights[n][0] + glat[src_add[n]] * weights[n][1] + glon[src_add[n]] * weights[n][2]
                 + glatlon[src_add[n]] * weights[n][3];
95
96

  return tgt_point;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
98
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
101
102
103
104
105
/*
  -----------------------------------------------------------------------

  This routine computes the weights for a bicubic interpolation.

  -----------------------------------------------------------------------
*/
106
void
107
remapBicubicWeights(RemapSearch &rsearch, RemapVars &rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
{
109
110
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;
111

112
  if (Options::cdoVerbose) cdoPrint("Called %s()", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
113

Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
  if (src_grid->rank != 2) cdoAbort("Can't do bicubic interpolation when source grid rank != 2");

116
  auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117

118
  progress::init();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119

120
  // Compute mappings from source to target grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121

122
  auto tgt_grid_size = tgt_grid->size;
123

124
  std::vector<WeightLinks4> weightLinks(tgt_grid_size);
125
  weightLinks4Alloc(tgt_grid_size, weightLinks);
126

127
  auto findex = 0.0;
128

Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
  // Loop over destination grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130

131
132
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(findex, rsearch, weightLinks, tgt_grid_size, src_grid, tgt_grid, rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
#endif
134
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
    {
136
137
138
#ifdef _OPENMP
#pragma omp atomic
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
      findex++;
140
      if (cdo_omp_get_thread_num() == 0) progress::update(0, 1, findex / tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141

142
      weightLinks[tgt_cell_add].nlinks = 0;
143

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

146
      const auto llpoint = remapgrid_get_lonlat(tgt_grid, tgt_cell_add);
147

Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
149
150
151
      size_t 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 weights[4][4]; //  bicubic weights for four corners
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152

153
      // Find nearest square of grid points on source grid
154
      auto search_result = remapSearchSquare(rsearch, llpoint, src_add, src_lats, src_lons);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
155

156
      // Check to see if points are mask points
157
158
159
160
161
      if (search_result > 0)
        {
          for (unsigned n = 0; n < 4; ++n)
            if (!src_grid->mask[src_add[n]]) search_result = 0;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162

163
      // If point found, find local iw,jw coordinates for weights
164
165
      if (search_result > 0)
        {
166
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167

168
          double iw = 0.0, jw = 0.0;  // current guess for bilinear coordinate
169
          if (remapFindWeights(llpoint, src_lons, src_lats, &iw, &jw))
170
            {
171
              // Successfully found iw,jw - compute weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
              bicubicSetWeights(iw, jw, weights);
              store_weightlinks_bicubic(src_add, weights, tgt_cell_add, weightLinks);
174
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
          else
176
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
              bicubicWarning();
178
179
180
              search_result = -1;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
183
      // Search for bicubic failed - use a distance-weighted average instead
      // (this is typically near the pole) Distance was stored in src_lats!
184
185
      if (search_result < 0)
        {
186
          if (num_src_points(src_grid->mask, src_add, src_lats) > 0)
187
188
            {
              tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
190
              renormalizeWeights(src_lats, weights);
              store_weightlinks_bicubic(src_add, weights, tgt_cell_add, weightLinks);
191
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
        }
193
194
    }

195
  progress::update(0, 1, 1);
196

197
  weightLinks4ToRemapLinks(tgt_grid_size, weightLinks, rv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

199
  if (Options::cdoVerbose) cdoPrint("Square search nearest: %.2f seconds", cdo_get_wtime() - start);
200
}  // remapBicubicWeights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
202
203
204

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
  This routine computes and apply the weights for a bicubic interpolation.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
207
208

  -----------------------------------------------------------------------
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209

210
void
211
remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *restrict tgt_array, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
{
213
214
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;
215

216
  if (Options::cdoVerbose) cdoPrint("Called %s()", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217

Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
  if (src_grid->rank != 2) cdoAbort("Can't do bicubic interpolation when source grid rank != 2");

220
  auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221

222
  progress::init();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223

224
225
  auto tgt_grid_size = tgt_grid->size;
  auto src_grid_size = src_grid->size;
226

Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
  Varray<uint8_t> src_grid_mask(src_grid_size);
228
229
230
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(src_grid_size, src_array, src_grid_mask, missval)
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
  for (size_t i = 0; i < src_grid_size; ++i) src_grid_mask[i] = !DBL_IS_EQUAL(src_array[i], missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232

233
  // Compute mappings from source to target grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234

235
  RemapGradients gradients(src_grid->size);
236
  remapGradients(*src_grid, src_grid_mask, src_array, gradients);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237

238
  auto findex = 0.0;
239

Uwe Schulzweida's avatar
Uwe Schulzweida committed
240
  // Loop over destination grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241

242
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
#pragma omp parallel for default(none) \
244
    shared(findex, rsearch, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, src_grid_mask, gradients)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
#endif
246
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247
    {
248
249
250
#ifdef _OPENMP
#pragma omp atomic
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
      findex++;
252
      if (cdo_omp_get_thread_num() == 0) progress::update(0, 1, findex / tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253

254
      tgt_array[tgt_cell_add] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255

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

258
      const auto llpoint = remapgrid_get_lonlat(tgt_grid, tgt_cell_add);
259

Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
261
262
263
      size_t 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 weights[4][4]; //  bicubic weights for four corners
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264

265
      // Find nearest square of grid points on source grid
266
      auto search_result = remapSearchSquare(rsearch, llpoint, src_add, src_lats, src_lons);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267

268
      // Check to see if points are mask points
269
270
271
      if (search_result > 0)
        {
          for (unsigned n = 0; n < 4; ++n)
272
            if (!src_grid_mask[src_add[n]]) search_result = 0;
273
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274

275
      // If point found, find local iw,jw coordinates for weights
276
277
      if (search_result > 0)
        {
278
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279

280
          double iw = 0.0, jw = 0.0;  // current guess for bilinear coordinate
281
          if (remapFindWeights(llpoint, src_lons, src_lats, &iw, &jw))
282
            {
283
              // Successfully found iw,jw - compute weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
              bicubicSetWeights(iw, jw, weights);
285
              sort_weights_bicubic(src_add, weights);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
              tgt_array[tgt_cell_add] = bicubicRemap(src_array, weights, src_add, gradients);
287
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
          else
289
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
              bicubicWarning();
291
292
293
              search_result = -1;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294

295
296
      // Search for bicubic failed - use a distance-weighted average instead
      // (this is typically near the pole) Distance was stored in src_lats!
297
298
      if (search_result < 0)
        {
299
          if (num_src_points(src_grid_mask, src_add, src_lats) > 0)
300
301
            {
              tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
              renormalizeWeights(src_lats, weights);
303
              sort_weights_bicubic(src_add, weights);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
              tgt_array[tgt_cell_add] = bicubicRemap(src_array, weights, src_add, gradients);
305
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306
        }
307
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308

309
  progress::update(0, 1, 1);
310

311
  if (Options::cdoVerbose) cdoPrint("Square search nearest: %.2f seconds", cdo_get_wtime() - start);
312
}  // remapBicubic