Commit 62c83442 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

remapnn/dis: replaced scrip search by kdtree

parent 12363d73
......@@ -2,9 +2,17 @@
#include <stdlib.h>
#include <math.h>
#include "dmemory.h"
#include "util.h"
#include "grid_search.h"
#include "kdtreelib/kdtree.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288 /* pi */
#endif
#define PI M_PI
#define PI2 (2.0*PI)
#define KDTREELIB
......@@ -19,17 +27,54 @@ static inline void LLtoXYZ_f(double lon, double lat, float *restrict xyz)
}
struct gridsearch {
unsigned n;
struct kdNode *kdt;
};
struct gridsearch *gridsearch_create_reg2d(unsigned nx, unsigned ny, const double *restrict lons, const double *restrict lats)
{
struct gridsearch *gs = (struct gridsearch *) calloc(1, sizeof(struct gridsearch));
gs->nx = nx;
gs->ny = ny;
double *reg2d_center_lon = (double *) malloc((nx+1)*sizeof(double));
double *reg2d_center_lat = (double *) malloc(ny*sizeof(double));
memcpy(reg2d_center_lon, lons, (nx+1)*sizeof(double));
memcpy(reg2d_center_lat, lats, ny*sizeof(double));
double *coslon = (double *) malloc(nx*sizeof(double));
double *sinlon = (double *) malloc(nx*sizeof(double));
double *coslat = (double *) malloc(ny*sizeof(double));
double *sinlat = (double *) malloc(ny*sizeof(double));
for ( unsigned n = 0; n < nx; ++n )
{
double rlon = lons[n];
if ( rlon > PI2 ) rlon -= PI2;
if ( rlon < 0 ) rlon += PI2;
coslon[n] = cos(rlon);
sinlon[n] = sin(rlon);
}
for ( unsigned n = 0; n < ny; ++n )
{
coslat[n] = cos(lats[n]);
sinlat[n] = sin(lats[n]);
}
gs->reg2d_center_lon = reg2d_center_lon;
gs->reg2d_center_lat = reg2d_center_lat;
gs->coslon = coslon;
gs->sinlon = sinlon;
gs->coslat = coslat;
gs->sinlat = sinlat;
return gs;
}
struct gridsearch *gridsearch_index_create(unsigned n, const double *restrict lons, const double *restrict lats, const unsigned *restrict index)
{
struct gridsearch *gs = (struct gridsearch *) malloc(sizeof(struct gridsearch));
struct gridsearch *gs = (struct gridsearch *) calloc(1, sizeof(struct gridsearch));
gs->n = n;
gs->kdt = NULL;
if ( n == 0 ) return gs;
......@@ -84,9 +129,9 @@ struct gridsearch *gridsearch_index_create(unsigned n, const double *restrict lo
struct gridsearch *gridsearch_create(unsigned n, const double *restrict lons, const double *restrict lats)
{
struct gridsearch *gs = (struct gridsearch *) malloc(sizeof(struct gridsearch));
struct gridsearch *gs = (struct gridsearch *) calloc(1, sizeof(struct gridsearch));
gs->n = n;
gs->kdt = NULL;
if ( n == 0 ) return gs;
......@@ -143,8 +188,16 @@ void gridsearch_delete(struct gridsearch *gs)
{
if ( gs )
{
kd_destroyTree(gs->kdt);
if ( gs->kdt ) kd_destroyTree(gs->kdt);
if ( gs->reg2d_center_lon ) free(gs->reg2d_center_lon);
if ( gs->reg2d_center_lat ) free(gs->reg2d_center_lat);
if ( gs->coslat ) free(gs->coslat);
if ( gs->coslon ) free(gs->coslon);
if ( gs->sinlat ) free(gs->sinlat);
if ( gs->sinlon ) free(gs->sinlon);
free(gs);
}
}
......@@ -164,7 +217,6 @@ void *gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double *
* than once around the sphere. The content
* of this variable is replaced with the
* distance to the NN squared. */
float range = range0;
#if defined(KD_SEARCH_SPH)
......@@ -200,7 +252,6 @@ struct pqueue *gridsearch_qnearest(struct gridsearch *gs, double lon, double lat
* than once around the sphere. The content
* of this variable is replaced with the
* distance to the NN squared. */
float range = range0;
LLtoXYZ_f(lon, lat, point);
......
......@@ -3,8 +3,17 @@
#include "kdtreelib/kdtree.h"
struct gridsearch;
struct gridsearch {
unsigned n;
unsigned nx, ny;
struct kdNode *kdt;
// reg2d search
double *reg2d_center_lon, *reg2d_center_lat;
double *coslat, *sinlat; // cosine, sine of grid lats (for distance)
double *coslon, *sinlon; // cosine, sine of grid lons (for distance)
};
struct gridsearch *gridsearch_create_reg2d(unsigned nx, unsigned ny, const double *restrict lons, const double *restrict lats);
struct gridsearch *gridsearch_create(unsigned n, const double *restrict lons, const double *restrict lats);
struct gridsearch *gridsearch_index_create(unsigned n, const double *restrict lons, const double *restrict lats, const unsigned *restrict index);
void gridsearch_delete(struct gridsearch *gs);
......
......@@ -274,6 +274,12 @@ kd_nearest(struct kdNode *node, float *p, float *max_dist_sq, int dim)
dist_sq = tmp_dist_sq;
*max_dist_sq = kd_min(dist_sq, *max_dist_sq);
}
// Uwe Schulzweida
else if (tmp_dist_sq <= dist_sq && tmp_nearest->index < nearest->index) {
nearest = tmp_nearest;
dist_sq = tmp_dist_sq;
*max_dist_sq = kd_min(dist_sq, *max_dist_sq);
}
}
return nearest;
}
......@@ -324,16 +330,13 @@ kd_qnearest(struct kdNode *node, float *p,
*
* return 1 if okay, zero in case of problems
*/
int
kd_doQnearest(struct kdNode *node, float *p,
// Uwe Schulzweida: extract kd_check_dist() from kd_doQnearest()
static int
kd_check_dist(struct kdNode *node, float *p,
float *max_dist_sq, unsigned int q, int dim, struct pqueue *res)
{
struct kdNode *nearer, *further;
struct resItem *point, *item;
float dist_sq, dx;
if (!node)
return 1;
float dist_sq;
dist_sq = kd_dist_sq(node->location, p, dim);
if (dist_sq < *max_dist_sq && kd_isleaf(node)) {
......@@ -344,6 +347,17 @@ kd_doQnearest(struct kdNode *node, float *p,
point->dist_sq = dist_sq;
pqinsert(res, point);
}
/*
else if (res->size > 1 && dist_sq <= *max_dist_sq && kd_isleaf(node)) {
pqpeek_max(res, &item);
if ( node->index < item->node->index )
{
item->node = node;
item->dist_sq = dist_sq;
}
}
*/
if (res->size > q) {
pqremove_max(res, &item);
free(item);
......@@ -361,6 +375,20 @@ kd_doQnearest(struct kdNode *node, float *p,
}
}
return 1;
}
int
kd_doQnearest(struct kdNode *node, float *p,
float *max_dist_sq, unsigned int q, int dim, struct pqueue *res)
{
struct kdNode *nearer, *further;
float dx;
if (!node) return 1;
if (!kd_check_dist(node, p, max_dist_sq, q, dim, res)) return 0;
if (p[node->split] < node->location[node->split]) {
nearer = node->left;
further = node->right;
......@@ -368,11 +396,9 @@ kd_doQnearest(struct kdNode *node, float *p,
nearer = node->right;
further = node->left;
}
if (!kd_doQnearest(nearer, p, max_dist_sq, q, dim, res))
return 0;
if (!kd_doQnearest(nearer, p, max_dist_sq, q, dim, res)) return 0;
if (!further)
return 1;
if (!further) return 1;
dx = kd_min(fabs(p[node->split] - further->min[node->split]),
fabs(p[node->split] - further->max[node->split]));
......@@ -381,34 +407,9 @@ kd_doQnearest(struct kdNode *node, float *p,
* some part of the further hyper-rectangle is in the search
* radius, search the further node
*/
if (!kd_doQnearest(further, p, max_dist_sq, q, dim, res))
return 0;
dist_sq = kd_dist_sq(node->location, p, dim);
if (!kd_doQnearest(further, p, max_dist_sq, q, dim, res)) return 0;
if (dist_sq < *max_dist_sq && kd_isleaf(node)) {
if ((point = kd_malloc(sizeof(struct resItem), "kd_doQnearest: "))
== NULL)
return 0;
point->node = node;
point->dist_sq = dist_sq;
pqinsert(res, point);
}
if (res->size > q) {
pqremove_max(res, &item);
free(item);
if (res->size > 1) {
/*
* Only inspect the queue if there are items left
*/
pqpeek_max(res, &item);
*max_dist_sq = item->dist_sq;
} else {
/*
* Nothing was found within the max search radius
*/
*max_dist_sq = 0;
}
}
if (!kd_check_dist(node, p, max_dist_sq, q, dim, res)) return 0;
}
return 1;
}
......
......@@ -102,7 +102,58 @@ _compPoints(const void *p1, const void *p2, int axis)
return 0;
}
}
/*
static int
_compPoints0(const void *p1, const void *p2)
{
struct kd_point *a = (struct kd_point *) p1;
struct kd_point *b = (struct kd_point *) p2;
if (a->point[0] > b->point[0]) return 1;
else if (a->point[0] < b->point[0]) return -1;
else
{
if ( a->index > b->index ) return 1;
else if ( a->index < b->index ) return -1;
return 0;
}
}
static int
_compPoints1(const void *p1, const void *p2)
{
struct kd_point *a = (struct kd_point *) p1;
struct kd_point *b = (struct kd_point *) p2;
if (a->point[1] > b->point[1]) return 1;
else if (a->point[1] < b->point[1]) return -1;
else
{
if ( a->index > b->index ) return 1;
else if ( a->index < b->index ) return -1;
return 0;
}
}
static int
_compPoints2(const void *p1, const void *p2)
{
struct kd_point *a = (struct kd_point *) p1;
struct kd_point *b = (struct kd_point *) p2;
if (a->point[2] > b->point[2]) return 1;
else if (a->point[2] < b->point[2]) return -1;
else
{
if ( a->index > b->index ) return 1;
else if ( a->index < b->index ) return -1;
return 0;
}
}
*/
void *
kd_doBuildTree(void *threadarg)
{
......@@ -131,6 +182,7 @@ kd_doBuildTree(void *threadarg)
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
sortaxis = depth % dim;
if (nPoints == 1) {
if ((node = kd_allocNode(points, 0, min, max, sortaxis, dim)) == NULL)
return NULL;
......@@ -142,6 +194,15 @@ kd_doBuildTree(void *threadarg)
* If this iteration is allowed to start more threads, we first
* use them to parallelize the sorting
*/
/*
if ( max_threads == 1 )
{
if ( sortaxis == 0 ) qsort(points, nPoints, sizeof(struct kd_point), _compPoints0);
else if ( sortaxis == 1 ) qsort(points, nPoints, sizeof(struct kd_point), _compPoints1);
else if ( sortaxis == 2 ) qsort(points, nPoints, sizeof(struct kd_point), _compPoints2);
}
else
*/
pmergesort(points, nPoints, sizeof(struct kd_point), _compPoints, sortaxis, max_threads);
pivot = nPoints / 2;
if ((node = kd_allocNode(points, pivot, min, max, sortaxis, dim)) == NULL)
......
This diff is collapsed.
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