Commit 2ffd85b0 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapBicubic: optimized for changing masks.

parent 1b236ea5
......@@ -3,6 +3,10 @@
* Using CDI library version 1.9.8
* Version 1.9.8 release
2019-07-18 Uwe Schulzweida
* remapBicubic: optimized for changing masks
2019-07-17 Uwe Schulzweida
* remapBilinear: optimized for changing masks
......
......@@ -929,7 +929,7 @@ Remap(void *argument)
if (!remap_genweights && (mapType == RemapMethod::CONSERV_SCRIP || mapType == RemapMethod::CONSERV)) remap_genweights = true;
bool useMask = !(!remap_genweights && mapType == RemapMethod::BILINEAR);
bool useMask = !(!remap_genweights && (mapType == RemapMethod::BILINEAR || mapType == RemapMethod::BICUBIC));
remap_set_int(REMAP_GENWEIGHTS, (int) remap_genweights);
......
......@@ -165,6 +165,7 @@ void remapConserv(RemapSearch &rsearch, const double *src_array, double *tgt_arr
void remapStat(int remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv, const double *array1, const double *array2,
double missval);
void remapGradients(RemapGrid &grid, const std::vector<bool> &mask, const double *array, RemapGradients &gradients);
void remapGradients(RemapGrid &grid, const double *array, RemapGradients &gradients);
void sort_add(size_t num_links, size_t num_wts, size_t *add1, size_t *add2, double *weights);
......
......@@ -49,6 +49,7 @@ bicubicSetWeights(double iw, double jw, double wgts[4][4])
// clang-format on
}
unsigned num_src_points(const std::vector<bool> &mask, const size_t src_add[4], double src_lats[4]);
unsigned num_src_points(const int *restrict mask, const size_t src_add[4], double src_lats[4]);
static void
......@@ -146,7 +147,7 @@ remapBicubicWeights(RemapSearch &rsearch, RemapVars &rv)
// Find nearest square of grid points on source grid
int search_result = remapSearchSquare(rsearch, plon, plat, src_add, src_lats, src_lons);
// Check to see if points are land points
// Check to see if points are mask points
if (search_result > 0)
{
for (unsigned n = 0; n < 4; ++n)
......@@ -217,11 +218,19 @@ remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *res
Progress::init();
size_t tgt_grid_size = tgt_grid->size;
size_t src_grid_size = src_grid->size;
std::vector<bool> src_grid_mask(src_grid_size);
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(src_grid_size, src_array, src_grid_mask, missval)
#endif
for (size_t i = 0; i < src_grid_size; ++i)
src_grid_mask[i] = !DBL_IS_EQUAL(src_array[i], missval);
// Compute mappings from source to target grid
RemapGradients gradients(src_grid->size);
remapGradients(*src_grid, src_array, gradients);
remapGradients(*src_grid, src_grid_mask, src_array, gradients);
double findex = 0;
......@@ -229,7 +238,7 @@ remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *res
#ifdef _OPENMP
#pragma omp parallel for default(none) \
shared(findex, rsearch, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, gradients)
shared(findex, rsearch, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, src_grid_mask, gradients)
#endif
for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
{
......@@ -254,11 +263,11 @@ remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *res
// Find nearest square of grid points on source grid
int search_result = remapSearchSquare(rsearch, plon, plat, src_add, src_lats, src_lons);
// Check to see if points are land points
// Check to see if points are mask points
if (search_result > 0)
{
for (unsigned n = 0; n < 4; ++n)
if (!src_grid->mask[src_add[n]]) search_result = 0;
if (!src_grid_mask[src_add[n]]) search_result = 0;
}
// If point found, find local iw,jw coordinates for weights
......@@ -287,7 +296,7 @@ remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *res
*/
if (search_result < 0)
{
if (num_src_points(&src_grid->mask[0], src_add, src_lats) > 0)
if (num_src_points(src_grid_mask, src_add, src_lats) > 0)
{
tgt_grid->cell_frac[tgt_cell_add] = 1.;
renormalizeWeights(src_lats, wgts);
......
......@@ -104,7 +104,7 @@ bilinearSetWeights(double iw, double jw, double wgts[4])
// clang-format on
}
static unsigned
unsigned
num_src_points(const std::vector<bool> &mask, const size_t src_add[4], double src_lats[4])
{
unsigned icount = 0;
......
......@@ -846,7 +846,7 @@ remapStat(int remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &r
/*****************************************************************************/
void
remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gradients)
remapGradients(RemapGrid &grid, const std::vector<bool> &mask, const double *restrict array, RemapGradients &gradients)
{
if (grid.rank != 2) cdoAbort("Internal problem (remapGradients), grid rank = %d!", grid.rank);
......@@ -859,7 +859,7 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
size_t ny = grid.dims[1];
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(grid_size, grad_lat, grad_lon, grad_latlon, grid, nx, ny, array)
#pragma omp parallel for default(none) shared(grid_size, grad_lat, grad_lon, grad_latlon, grid, nx, ny, array, mask)
#endif
for (size_t n = 0; n < grid_size; ++n)
{
......@@ -867,7 +867,7 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
grad_lon[n] = 0.0;
grad_latlon[n] = 0.0;
if (grid.mask[n])
if (mask[n])
{
double delew = 0.5;
double delns = 0.5;
......@@ -904,12 +904,12 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
size_t isw = (jm1 - 1) * nx + im1 - 1;
// Compute i-gradient
if (!grid.mask[ie])
if (!mask[ie])
{
ie = n;
delew = 1.0;
}
if (!grid.mask[iw])
if (!mask[iw])
{
iw = n;
delew = 1.0;
......@@ -918,12 +918,12 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
grad_lat[n] = delew * (array[ie] - array[iw]);
// Compute j-gradient
if (!grid.mask[in])
if (!mask[in])
{
in = n;
delns = 1.0;
}
if (!grid.mask[is])
if (!mask[is])
{
is = n;
delns = 1.0;
......@@ -935,7 +935,7 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
delew = 0.5;
delns = (jp1 == j || jm1 == j) ? 1.0 : 0.5;
if (!grid.mask[ine])
if (!mask[ine])
{
if (in != n)
{
......@@ -958,7 +958,7 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
}
}
if (!grid.mask[inw])
if (!mask[inw])
{
if (in != n)
{
......@@ -983,7 +983,7 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
double grad_lat_zero = delew * (array[ine] - array[inw]);
if (!grid.mask[ise])
if (!mask[ise])
{
if (is != n)
{
......@@ -1006,7 +1006,7 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
}
}
if (!grid.mask[isw])
if (!mask[isw])
{
if (is != n)
{
......@@ -1036,6 +1036,20 @@ remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gr
}
} /* remapGradients */
void
remapGradients(RemapGrid &grid, const double *restrict array, RemapGradients &gradients)
{
size_t grid_size = grid.size;
std::vector<bool> mask(grid_size);
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(grid_size, grid, mask)
#endif
for (size_t i = 0; i < grid_size; ++i)
mask[i] = grid.mask[i] > 0;
remapGradients(grid, mask, array, gradients);
}
/*****************************************************************************/
void
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment