Commit 78f68252 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Refactor sort_add_and_wgts_bic().

parent d91a8548
......@@ -29,23 +29,29 @@
static void
bicubicSetWeights(double iw, double jw, double (&wgts)[4][4])
{
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);
// clang-format off
wgts[0][0] = (1.-jw*jw*(3.-2.*jw)) * (1.-iw*iw*(3.-2.*iw));
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.);
wgts[1][1] = (1.-jw*jw*(3.-2.*jw)) * iw*iw*(iw-1.);
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);
wgts[2][2] = jw*jw*(jw-1.) * iw*iw*(3.-2.*iw);
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.);
wgts[3][3] = iw*(iw-1.)*(iw-1.) * jw*jw*(jw-1.);
wgts[0][0] = (1.-jw3) * (1.-iw3);
wgts[1][0] = (1.-jw3) * iw3;
wgts[2][0] = jw3 * iw3;
wgts[3][0] = jw3 * (1.-iw3);
wgts[0][1] = (1.-jw3) * iw2;
wgts[1][1] = (1.-jw3) * iw1;
wgts[2][1] = jw3 * iw1;
wgts[3][1] = jw3 * iw2;
wgts[0][2] = jw2 * (1.-iw3);
wgts[1][2] = jw2 * iw3;
wgts[2][2] = jw1 * iw3;
wgts[3][2] = jw1 * (1.-iw3);
wgts[0][3] = jw2 * iw2;
wgts[1][3] = jw2 * iw1;
wgts[2][3] = jw1 * iw1;
wgts[3][3] = jw1 * iw2;
// clang-format on
}
......@@ -276,7 +282,7 @@ remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *res
{
// Successfully found iw,jw - compute weights
bicubicSetWeights(iw, jw, wgts);
sort_add_and_wgts4(4, src_add, wgts);
sort_add_and_wgts_bic(src_add, wgts);
tgt_array[tgt_cell_add] = bicubicRemap(src_array, wgts, src_add, gradients);
}
else
......@@ -294,7 +300,7 @@ remapBicubic(RemapSearch &rsearch, const double *restrict src_array, double *res
{
tgt_grid->cell_frac[tgt_cell_add] = 1.;
renormalizeWeights(src_lats, wgts);
sort_add_and_wgts4(4, src_add, wgts);
sort_add_and_wgts_bic(src_add, wgts);
tgt_array[tgt_cell_add] = bicubicRemap(src_array, wgts, src_add, gradients);
}
}
......
......@@ -16,6 +16,7 @@
*/
#include <algorithm>
#include <array>
#include "dmemory.h"
#include "cdo_options.h"
......@@ -101,30 +102,28 @@ sort_add_and_wgts(size_t numWeights, size_t *src_add, double *wgts)
}
void
sort_add_and_wgts4(size_t numWeights, size_t *src_add, double wgts[4][4])
sort_add_and_wgts_bic(size_t *src_add, double (&wgts)[4][4])
{
constexpr size_t numWeights = 4;
size_t n;
for (n = 1; n < numWeights; ++n)
if (src_add[n] < src_add[n - 1]) break;
if (n == numWeights) return;
if (numWeights > 1)
{
std::vector<Addweight4> addweights(numWeights);
std::array<Addweight4, numWeights> addweights;
for (n = 0; n < numWeights; ++n)
{
addweights[n].add = src_add[n];
for (unsigned k = 0; k < 4; ++k) addweights[n].weight[k] = wgts[n][k];
}
for (n = 0; n < numWeights; ++n)
{
addweights[n].add = src_add[n];
for (unsigned k = 0; k < 4; ++k) addweights[n].weight[k] = wgts[n][k];
}
std::sort(addweights.begin(), addweights.end(), compareAdds4);
std::sort(addweights.begin(), addweights.end(), compareAdds4);
for (n = 0; n < numWeights; ++n)
{
src_add[n] = addweights[n].add;
for (unsigned k = 0; k < 4; ++k) wgts[n][k] = addweights[n].weight[k];
}
for (n = 0; n < numWeights; ++n)
{
src_add[n] = addweights[n].add;
for (unsigned k = 0; k < 4; ++k) wgts[n][k] = addweights[n].weight[k];
}
}
......
......@@ -56,6 +56,6 @@ void storeWeightLinks4(size_t *srch_add, double weights[4][4], size_t cell_add,
void weightLinksToRemapLinks(int lalloc, size_t gridSize, std::vector<WeightLinks> &weightLinks, RemapVars &rv);
void weightLinks4ToRemapLinks(size_t gridSize, std::vector<WeightLinks4> &weightLinks, RemapVars &rv);
void sort_add_and_wgts(size_t numWeights, size_t *src_add, double *wgts);
void sort_add_and_wgts4(size_t numWeights, size_t *src_add, double wgts[4][4]);
void sort_add_and_wgts_bic(size_t *src_add, double (&wgts)[4][4]);
#endif /* REMAP_STORE_LINK */
Markdown is supported
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