Skip to content
GitLab
Commit 2b5531f7 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added quad_cross_products().

parent 21afd5ba
Hide whitespace changes
Inline Side-by-side
src/remap.h
View file @ 2b5531f7
......@@ -226,6 +226,8 @@ int grid_search_reg2d(remapgrid_t *src_grid, size_t *restrict src_add, double *r
double *restrict src_lons, double plat, double plon, const size_t *restrict src_grid_dims,
const double *restrict src_center_lat, const double *restrict src_center_lon);
unsigned quad_cross_products(double plon, double plat, double lons[4], double lats[4]);
int grid_search(remapgrid_t *src_grid, size_t *restrict src_add, double *restrict src_lats,
double *restrict src_lons, double plat, double plon, const size_t *restrict src_grid_dims,
const double *restrict src_center_lat, const double *restrict src_center_lon,
......
src/remap_bilinear_scrip.cc
View file @ 2b5531f7
......@@ -282,6 +282,98 @@ void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid,
-----------------------------------------------------------------------
*/
//#define TEST_KDTREE
#ifdef TEST_KDTREE
#include "grid_search.h"
int grid_search_test(struct gridsearch *gs, size_t *restrict src_add, double *restrict src_lats,
double *restrict src_lons, double plat, double plon, const size_t *restrict src_grid_dims,
double *restrict src_center_lat, double *restrict src_center_lon)
{
/*
Output variables:
int src_add[4] ! address of each corner point enclosing P
double src_lats[4] ! latitudes of the four corner points
double src_lons[4] ! longitudes of the four corner points
Input variables:
double plat ! latitude of the search point
double plon ! longitude of the search point
int src_grid_dims[2] ! size of each src grid dimension
*/
bool is_cyclic = true;
int search_result = 0;
for ( unsigned n = 0; n < 4; ++n ) src_add[n] = 0;
/* Now perform a more detailed search */
size_t nx = src_grid_dims[0];
size_t ny = src_grid_dims[1];
double search_radius = gs->search_radius;
const double range0 = SQR(search_radius);
double range = range0;
size_t add = gridsearch_nearest(gs, plon, plat, &range);
// printf("plon, plat, add, range %g %g %g %g %zu %g\n", plon*RAD2DEG, plat*RAD2DEG,
// src_center_lon[add]*RAD2DEG, src_center_lat[add]*RAD2DEG,add, range);
if ( add != GS_NOT_FOUND )
{
size_t idx[4];
for ( unsigned k = 0; k < 4; ++k )
{
/* Determine neighbor addresses */
size_t j = add/nx;
size_t i = add - j*nx;
if ( k == 1 || k == 3 ) i = (i > 0) ? i - 1 : (is_cyclic) ? nx-1 : 0;
if ( k == 2 || k == 3 ) j = (j > 0) ? j - 1 : 0;
size_t ip1 = (i < (nx-1)) ? i + 1 : (is_cyclic) ? 0 : i;
size_t jp1 = (j < (ny-1)) ? j + 1 : j;
idx[0] = j *nx + i;
idx[1] = j *nx + ip1; // east
idx[2] = jp1*nx + ip1; // north-east
idx[3] = jp1*nx + i; // north
for ( unsigned j = 0; j < 4; ++j ) src_lons[j] = src_center_lon[idx[j]];
for ( unsigned j = 0; j < 4; ++j ) src_lats[j] = src_center_lat[idx[j]];
unsigned n = quad_cross_products(plon, plat, src_lons, src_lats);
/* If cross products all same sign, we found the location */
if ( n >= 4 )
{
for ( unsigned j = 0; j < 4; ++j ) src_add[j] = idx[j];
search_result = 1;
return search_result;
}
/* Otherwise move on to next cell */
}
/*
If no cell found, point is likely either in a box that straddles either pole or is outside
the grid. Fall back to a distance-weighted average of the four closest points.
Go ahead and compute weights here, but store in src_lats and return -add to prevent the
parent routine from computing bilinear weights.
*/
// if ( !src_grid->lextrapolate ) return search_result;
/*
printf("Could not find location for %g %g\n", plat*RAD2DEG, plon*RAD2DEG);
printf("Using nearest-neighbor for this point\n");
*/
search_result = add;
}
return search_result;
} /* grid_search_test */
#endif
void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const double *restrict src_array, double *restrict tgt_array, double missval)
{
extern int timer_remap_bil;
......@@ -293,6 +385,12 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
progressInit();
#ifdef TEST_KDTREE
struct gridsearch *gs = NULL;
if (remap_grid_type != REMAP_GRID_TYPE_REG2D )
gs = gridsearch_create_nn(src_grid->size, src_grid->cell_center_lon, src_grid->cell_center_lat);
#endif
size_t tgt_grid_size = tgt_grid->size;
/* Compute mappings from source to target grid */
......@@ -335,10 +433,15 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
plat, plon, src_grid->dims,
src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
else
#ifdef TEST_KDTREE
search_result = grid_search_test(gs, src_add, src_lats, src_lons, plat, plon, src_grid->dims,
src_grid->cell_center_lat, src_grid->cell_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);
#endif
// Check to see if points are mask points
if ( search_result > 0 )
......@@ -389,5 +492,9 @@ void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const do
}
}
#ifdef TEST_KDTREE
if ( gs ) gridsearch_delete(gs);
#endif
if ( cdoTimer ) timer_stop(timer_remap_bil);
} // scrip_remap_bilinear
src/remap_search_latbins.cc
View file @ 2b5531f7
......@@ -222,6 +222,73 @@ int grid_search_nn(size_t min_add, size_t max_add, size_t *restrict nbr_add, dou
}
unsigned quad_cross_products(double plon, double plat, double lons[4], double lats[4])
{
unsigned n;
int scross[4], scross_last = 0;
// Vectors for cross-product check
double vec1_lat, vec1_lon;
double vec2_lat, vec2_lon;
/* For consistency, we must make sure all lons are in same 2pi interval */
vec1_lon = lons[0] - plon;
if ( vec1_lon > PI ) lons[0] -= PI2;
else if ( vec1_lon < -PI ) lons[0] += PI2;
for ( n = 1; n < 4; ++n )
{
vec1_lon = lons[n] - lons[0];
if ( vec1_lon > PI ) lons[n] -= PI2;
else if ( vec1_lon < -PI ) lons[n] += PI2;
}
/* corner_loop */
for ( n = 0; n < 4; ++n )
{
unsigned next_n = (n+1)%4;
/*
Here we take the cross product of the vector making up each box side
with the vector formed by the vertex and search point.
If all the cross products are positive, the point is contained in the box.
*/
vec1_lat = lats[next_n] - lats[n];
vec1_lon = lons[next_n] - lons[n];
vec2_lat = plat - lats[n];
vec2_lon = plon - lons[n];
/* Check for 0,2pi crossings */
if ( vec1_lon > THREE*PIH ) vec1_lon -= PI2;
else if ( vec1_lon < -THREE*PIH ) vec1_lon += PI2;
if ( vec2_lon > THREE*PIH ) vec2_lon -= PI2;
else if ( vec2_lon < -THREE*PIH ) vec2_lon += PI2;
double cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat;
/* If cross product is less than ZERO, this cell doesn't work */
/* 2008-10-16 Uwe Schulzweida: bug fix for cross_product eq zero */
scross[n] = cross_product < 0 ? -1 : cross_product > 0 ? 1 : 0;
if ( n == 0 ) scross_last = scross[n];
if ( (scross[n] < 0 && scross_last > 0) || (scross[n] > 0 && scross_last < 0) ) break;
scross_last = scross[n];
}
if ( n >= 4 )
{
n = 0;
if ( scross[0]>=0 && scross[1]>=0 && scross[2]>=0 && scross[3]>=0 ) n = 4;
else if ( scross[0]<=0 && scross[1]<=0 && scross[2]<=0 && scross[3]<=0 ) n = 4;
}
return n;
}
int grid_search(remapgrid_t *src_grid, size_t *restrict src_add, double *restrict src_lats,
double *restrict src_lons, double plat, double plon, const size_t *restrict src_grid_dims,
const double *restrict src_center_lat, const double *restrict src_center_lon,
......@@ -249,11 +316,7 @@ int grid_search(remapgrid_t *src_grid, size_t *restrict src_add, double *restric
int src_bin_add[][2] ! latitude bins for restricting
*/
/* Local variables */
size_t n2, next_n, srch_add, srch_add4; /* dummy indices */
/* Vectors for cross-product check */
double vec1_lat, vec1_lon;
double vec2_lat, vec2_lon;
int scross[4], scross_last = 0;
size_t n2, srch_add, srch_add4; /* dummy indices */
int search_result = 0;
restr_t *bin_lats = src_grid->bin_lats;
......@@ -285,6 +348,8 @@ int grid_search(remapgrid_t *src_grid, size_t *restrict src_add, double *restric
size_t nx = src_grid_dims[0];
size_t ny = src_grid_dims[1];
size_t idx[4];
/* srch_loop */
for ( srch_add = min_add; srch_add <= max_add; ++srch_add )
{
......@@ -295,8 +360,6 @@ int grid_search(remapgrid_t *src_grid, size_t *restrict src_add, double *restric
rlat >= src_grid_bound_box[srch_add4 ] &&
rlat <= src_grid_bound_box[srch_add4+1])
{
unsigned n;
/* We are within bounding box so get really serious */
/* Determine neighbor addresses */
......@@ -306,86 +369,22 @@ int grid_search(remapgrid_t *src_grid, size_t *restrict src_add, double *restric
size_t ip1 = (i < (nx-1)) ? i + 1 : (src_grid->is_cyclic) ? 0 : i;
size_t jp1 = (j < (ny-1)) ? j + 1 : j;
size_t n_add = jp1*nx + i;
size_t e_add = j *nx + ip1;
size_t ne_add = jp1*nx + ip1;
src_lons[0] = src_center_lon[srch_add];
src_lons[1] = src_center_lon[e_add];
src_lons[2] = src_center_lon[ne_add];
src_lons[3] = src_center_lon[n_add];
src_lats[0] = src_center_lat[srch_add];
src_lats[1] = src_center_lat[e_add];
src_lats[2] = src_center_lat[ne_add];
src_lats[3] = src_center_lat[n_add];
/* For consistency, we must make sure all lons are in same 2pi interval */
vec1_lon = src_lons[0] - plon;
if ( vec1_lon > PI ) src_lons[0] -= PI2;
else if ( vec1_lon < -PI ) src_lons[0] += PI2;
for ( n = 1; n < 4; ++n )
{
vec1_lon = src_lons[n] - src_lons[0];
if ( vec1_lon > PI ) src_lons[n] -= PI2;
else if ( vec1_lon < -PI ) src_lons[n] += PI2;
}
/* corner_loop */
for ( n = 0; n < 4; ++n )
{
next_n = (n+1)%4;
/*
Here we take the cross product of the vector making up each box side
with the vector formed by the vertex and search point.
If all the cross products are positive, the point is contained in the box.
*/
vec1_lat = src_lats[next_n] - src_lats[n];
vec1_lon = src_lons[next_n] - src_lons[n];
vec2_lat = plat - src_lats[n];
vec2_lon = plon - src_lons[n];
/* Check for 0,2pi crossings */
if ( vec1_lon > THREE*PIH ) vec1_lon -= PI2;
else if ( vec1_lon < -THREE*PIH ) vec1_lon += PI2;
if ( vec2_lon > THREE*PIH ) vec2_lon -= PI2;
else if ( vec2_lon < -THREE*PIH ) vec2_lon += PI2;
double cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat;
idx[0] = j *nx + i;
idx[1] = j *nx + ip1; // east
idx[2] = jp1*nx + ip1; // north-east
idx[3] = jp1*nx + i; // north
/* If cross product is less than ZERO, this cell doesn't work */
/* 2008-10-16 Uwe Schulzweida: bug fix for cross_product eq zero */
scross[n] = cross_product < 0 ? -1 : cross_product > 0 ? 1 : 0;
for ( unsigned j = 0; j < 4; ++j ) src_lons[j] = src_center_lon[idx[j]];
for ( unsigned j = 0; j < 4; ++j ) src_lats[j] = src_center_lat[idx[j]];
if ( n == 0 ) scross_last = scross[n];
if ( (scross[n] < 0 && scross_last > 0) || (scross[n] > 0 && scross_last < 0) ) break;
scross_last = scross[n];
} /* corner_loop */
if ( n >= 4 )
{
n = 0;
if ( scross[0]>=0 && scross[1]>=0 && scross[2]>=0 && scross[3]>=0 ) n = 4;
else if ( scross[0]<=0 && scross[1]<=0 && scross[2]<=0 && scross[3]<=0 ) n = 4;
}
unsigned n = quad_cross_products(plon, plat, src_lons, src_lats);
/* If cross products all same sign, we found the location */
if ( n >= 4 )
{
src_add[0] = srch_add;
src_add[1] = e_add;
src_add[2] = ne_add;
src_add[3] = n_add;
for ( unsigned j = 0; j < 4; ++j ) src_add[j] = idx[j];
search_result = 1;
return search_result;
}
......
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

Imprint | Privacy Policy