Commit a1bb211a authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added function scrip_remap_distwgt()

parent 33c688d2
......@@ -582,10 +582,6 @@ src/interpol.c -text
src/interpol.h -text
src/job.c -text
src/juldate.c -text
src/kdtree.c -text
src/kdtree.h -text
src/kdtree_fast.c -text
src/kdtree_fast.h -text
src/kdtreelib/COPYING.DJ -text
src/kdtreelib/kdtree.h -text
src/kdtreelib/kdtree_cartesian.c -text
......
......@@ -350,12 +350,17 @@ void setmisstonn(field_t *field1, field_t *field2, int maxfill)
start = clock();
double search_radius = M_PI;
double range = SQR(2*search_radius);
unsigned index;
void *gs_result;
#pragma omp parallel for private(gs_result) shared(mindex, array1, array2, xvals, yvals) private(index)
#pragma omp parallel for private(gs_result, index) shared(mindex, array1, array2, xvals, yvals, range)
for ( unsigned i = 0; i < nmiss; ++i )
{
gs_result = gridsearch_nearest(gs, xvals[mindex[i]], yvals[mindex[i]]);
double prange = range;
gs_result = gridsearch_nearest(gs, xvals[mindex[i]], yvals[mindex[i]], &prange);
if ( gs_result )
index = gridsearch_item(gs_result);
else
......
......@@ -70,10 +70,6 @@ libcdo_la_SOURCES = \
interpol.h \
job.c \
juldate.c \
kdtree.c \
kdtree.h \
kdtree_fast.c \
kdtree_fast.h \
grid_search.c \
grid_search.h \
kvlist.c \
......
......@@ -135,7 +135,6 @@ am_libcdo_la_OBJECTS = libcdo_la-cdo_pthread.lo libcdo_la-cdo_vlist.lo \
libcdo_la-griddes_h5.lo libcdo_la-griddes_nc.lo \
libcdo_la-hetaeta.lo libcdo_la-institution.lo \
libcdo_la-interpol.lo libcdo_la-job.lo libcdo_la-juldate.lo \
libcdo_la-kdtree.lo libcdo_la-kdtree_fast.lo \
libcdo_la-grid_search.lo libcdo_la-kvlist.lo libcdo_la-list.lo \
libcdo_la-merge_sort2.lo libcdo_la-modules.lo \
libcdo_la-namelist.lo libcdo_la-normal.lo \
......@@ -562,14 +561,13 @@ libcdo_la_SOURCES = cdo_int.h compare.h cdo_pthread.c cdo_vlist.c \
gradsdeslib.c gradsdeslib.h grid.c grid.h grid_area.c \
grid_gme.c grid_lcc.c grid_rot.c gridreference.c griddes.c \
griddes.h griddes_h5.c griddes_nc.c hetaeta.c hetaeta.h \
institution.c interpol.c interpol.h job.c juldate.c kdtree.c \
kdtree.h kdtree_fast.c kdtree_fast.h grid_search.c \
grid_search.h kvlist.c kvlist.h list.c list.h merge_sort2.c \
merge_sort2.h modules.c modules.h namelist.c namelist.h \
normal.c nth_element.c nth_element.h operator_help.h par_io.c \
par_io.h percentiles.c percentiles.h pipe.c pipe.h \
pragma_omp_atomic_update.h printinfo.h process.c process.h \
pstream.c pstream.h pstream_write.h pstream_int.h \
institution.c interpol.c interpol.h job.c juldate.c \
grid_search.c grid_search.h kvlist.c kvlist.h list.c list.h \
merge_sort2.c merge_sort2.h modules.c modules.h namelist.c \
namelist.h normal.c nth_element.c nth_element.h \
operator_help.h par_io.c par_io.h percentiles.c percentiles.h \
pipe.c pipe.h pragma_omp_atomic_update.h printinfo.h process.c \
process.h pstream.c pstream.h pstream_write.h pstream_int.h \
pthread_debug.c pthread_debug.h readline.c realtime.c remap.h \
remaplib.c remapsort.c remap_scrip_io.c remap_search_reg2d.c \
remap_search_latbins.c remap_store_link.c remap_store_link.h \
......@@ -1048,8 +1046,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-interpol.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-job.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-juldate.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-kdtree.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-kdtree_fast.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-kvlist.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-list.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-merge_sort2.Plo@am__quote@
......@@ -1417,20 +1413,6 @@ libcdo_la-juldate.lo: juldate.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libcdo_la-juldate.lo `test -f 'juldate.c' || echo '$(srcdir)/'`juldate.c
libcdo_la-kdtree.lo: kdtree.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libcdo_la-kdtree.lo -MD -MP -MF $(DEPDIR)/libcdo_la-kdtree.Tpo -c -o libcdo_la-kdtree.lo `test -f 'kdtree.c' || echo '$(srcdir)/'`kdtree.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libcdo_la-kdtree.Tpo $(DEPDIR)/libcdo_la-kdtree.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='kdtree.c' object='libcdo_la-kdtree.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libcdo_la-kdtree.lo `test -f 'kdtree.c' || echo '$(srcdir)/'`kdtree.c
libcdo_la-kdtree_fast.lo: kdtree_fast.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libcdo_la-kdtree_fast.lo -MD -MP -MF $(DEPDIR)/libcdo_la-kdtree_fast.Tpo -c -o libcdo_la-kdtree_fast.lo `test -f 'kdtree_fast.c' || echo '$(srcdir)/'`kdtree_fast.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libcdo_la-kdtree_fast.Tpo $(DEPDIR)/libcdo_la-kdtree_fast.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='kdtree_fast.c' object='libcdo_la-kdtree_fast.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libcdo_la-kdtree_fast.lo `test -f 'kdtree_fast.c' || echo '$(srcdir)/'`kdtree_fast.c
libcdo_la-grid_search.lo: grid_search.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libcdo_la-grid_search.lo -MD -MP -MF $(DEPDIR)/libcdo_la-grid_search.Tpo -c -o libcdo_la-grid_search.lo `test -f 'grid_search.c' || echo '$(srcdir)/'`grid_search.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libcdo_la-grid_search.Tpo $(DEPDIR)/libcdo_la-grid_search.Plo
......
......@@ -919,7 +919,9 @@ void *Remap(void *argument)
get_map_type(operfunc, &map_type, &submap_type, &num_neighbors, &remap_order);
}
if ( remap_genweights == FALSE && map_type != MAP_TYPE_BILINEAR && map_type != MAP_TYPE_BICUBIC && map_type != MAP_TYPE_CONSERV_YAC )
if ( remap_genweights == FALSE &&
map_type != MAP_TYPE_BILINEAR && map_type != MAP_TYPE_BICUBIC &&
map_type != MAP_TYPE_DISTWGT && map_type != MAP_TYPE_CONSERV_YAC )
remap_genweights = TRUE;
remap_set_int(REMAP_GENWEIGHTS, remap_genweights);
......@@ -1182,6 +1184,7 @@ void *Remap(void *argument)
{
if ( map_type == MAP_TYPE_BILINEAR ) scrip_remap_bilinear(&remaps[r].src_grid, &remaps[r].tgt_grid, array1, array2, missval);
else if ( map_type == MAP_TYPE_BICUBIC ) scrip_remap_bicubic(&remaps[r].src_grid, &remaps[r].tgt_grid, array1, array2, missval);
else if ( map_type == MAP_TYPE_DISTWGT ) scrip_remap_distwgt(num_neighbors, &remaps[r].src_grid, &remaps[r].tgt_grid, array1, array2, missval);
else if ( map_type == MAP_TYPE_CONSERV_YAC ) remap_conserv(&remaps[r].src_grid, &remaps[r].tgt_grid, array1, array2, missval);
}
......
......@@ -21,7 +21,6 @@ static inline void LLtoXYZ_f(double lon, double lat, float *restrict xyz)
struct gridsearch {
unsigned n;
struct kdNode *kdt;
};
......@@ -142,25 +141,32 @@ struct gridsearch *gridsearch_create(unsigned n, const double *restrict lons, co
void gridsearch_delete(struct gridsearch *gs)
{
kd_destroyTree(gs->kdt);
if ( gs )
{
kd_destroyTree(gs->kdt);
free(gs);
free(gs);
}
}
#define SQR(a) ((a)*(a))
void *gridsearch_nearest(struct gridsearch *gs, double lon, double lat)
void *gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double *prange)
{
if ( gs->kdt == NULL ) return NULL;
float point[3];
float range0 = SQR(2 * M_PI); /* This has to be bigger than the presumed
* maximum distance to the NN but smaller
* than once around the sphere. The content
* of this variable is replaced with the
* distance to the NN squared. */
//range0 = SQR(1./180. * M_PI);;
float range0;
if ( prange )
range0 = *prange;
else
range0 = SQR(2 * M_PI); /* This has to be bigger than the presumed
* maximum distance to the NN but smaller
* 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)
point[0] = lon;
point[1] = lat;
......@@ -171,9 +177,11 @@ void *gridsearch_nearest(struct gridsearch *gs, double lon, double lat)
kdNode *node = kd_nearest(gs->kdt, point, &range, 3);
// printf("range %g %g %g %p\n", lon, lat, range, node);
#endif
if ( !(range < range0) ) node = NULL;
#endif
if ( prange ) *prange = range;
return (void *) node;
}
......
......@@ -6,7 +6,7 @@ struct gridsearch;
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);
void *gridsearch_nearest(struct gridsearch *gs, double lon, double lat);
void *gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double *range);
unsigned gridsearch_item(void *gs_result);
#endif
/*
This file is part of ``kdtree'', a library for working with kd-trees.
Copyright (C) 2007-2011 John Tsiombikas <nuclear@member.fsf.org>
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. The name of the author may not be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
OF SUCH DAMAGE.
*/
/* single nearest neighbor search written by Tamas Nepusz <tamas@cs.rhul.ac.uk> */
/* U. Schulzweida, 20150627: removed float interface, set MAX_DIM=3 */
/* U. Schulzweida, 20150628: changed void *data to unsiged index */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "kdtree.h"
#define MAX_DIM 3
struct kdhyperrect {
int dim;
double min[MAX_DIM], max[MAX_DIM]; /* minimum/maximum coords */
};
struct kdnode {
double pos[MAX_DIM];
int dir;
unsigned index;
struct kdnode *left, *right; /* negative/positive side */
};
struct res_node {
struct kdnode *item;
double dist_sq;
struct res_node *next;
};
struct kdtree {
int dim;
struct kdnode *root;
struct kdhyperrect *rect;
};
struct kdres {
struct kdtree *tree;
struct res_node *rlist, *riter;
int size;
};
#define SQ(x) ((x) * (x))
static int insert_rec(struct kdnode **node, const double *pos, unsigned index, int dir, int dim);
static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq);
static void clear_results(struct kdres *set);
static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max);
static void hyperrect_free(struct kdhyperrect *rect);
static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect);
static void hyperrect_extend(struct kdhyperrect *rect, const double *pos);
static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos);
#define alloc_resnode() malloc(sizeof(struct res_node))
#define free_resnode(n) free(n)
struct kdtree *kd_create(int k)
{
struct kdtree *tree;
if(!(tree = malloc(sizeof *tree))) return 0;
tree->dim = k;
tree->root = 0;
tree->rect = 0;
return tree;
}
void kd_free(struct kdtree *tree)
{
if(tree)
{
kd_clear(tree);
free(tree);
}
}
static void clear_rec(struct kdnode *node)
{
if(!node) return;
clear_rec(node->left);
clear_rec(node->right);
free(node);
}
void kd_clear(struct kdtree *tree)
{
clear_rec(tree->root);
tree->root = 0;
if (tree->rect)
{
hyperrect_free(tree->rect);
tree->rect = 0;
}
}
static int insert_rec(struct kdnode **nptr, const double *pos, unsigned index, int dir, int dim)
{
int new_dir;
struct kdnode *node;
if (!*nptr)
{
if(!(node = malloc(sizeof *node))) return -1;
memcpy(node->pos, pos, dim * sizeof *node->pos);
node->index = index;
node->dir = dir;
node->left = node->right = 0;
*nptr = node;
return 0;
}
node = *nptr;
new_dir = (node->dir + 1) % dim;
if (pos[node->dir] < node->pos[node->dir])
return insert_rec(&(*nptr)->left, pos, index, new_dir, dim);
return insert_rec(&(*nptr)->right, pos, index, new_dir, dim);
}
int kd_insert(struct kdtree *tree, const double *pos, unsigned index)
{
if (insert_rec(&tree->root, pos, index, 0, tree->dim)) return -1;
if (tree->rect == 0)
tree->rect = hyperrect_create(tree->dim, pos, pos);
else
hyperrect_extend(tree->rect, pos);
return 0;
}
int kd_insert3(struct kdtree *tree, double x, double y, double z, unsigned index)
{
double buf[3];
buf[0] = x;
buf[1] = y;
buf[2] = z;
return kd_insert(tree, buf, index);
}
static
int find_nearest(struct kdnode *node, const double *pos, double range, struct res_node *list, int ordered, int dim)
{
double dx;
int i, ret, added_res = 0;
if ( !node ) return 0;
double dist_sq = 0;
for ( i=0; i<dim; i++ )
{
dist_sq += SQ(node->pos[i] - pos[i]);
}
if ( dist_sq <= SQ(range) )
{
if (rlist_insert(list, node, ordered ? dist_sq : -1.0) == -1) return -1;
added_res = 1;
}
dx = pos[node->dir] - node->pos[node->dir];
ret = find_nearest(dx <= 0.0 ? node->left : node->right, pos, range, list, ordered, dim);
if ( ret >= 0 && fabs(dx) < range )
{
added_res += ret;
ret = find_nearest(dx <= 0.0 ? node->right : node->left, pos, range, list, ordered, dim);
}
if ( ret == -1 ) return -1;
added_res += ret;
return added_res;
}
#if 0
static int find_nearest_n(struct kdnode *node, const double *pos, double range, int num, struct rheap *heap, int dim)
{
double dist_sq, dx;
int i, ret, added_res = 0;
if(!node) return 0;
/* if the photon is close enough, add it to the result heap */
dist_sq = 0;
for(i=0; i<dim; i++) {
dist_sq += SQ(node->pos[i] - pos[i]);
}
if(dist_sq <= range_sq) {
if(heap->size >= num) {
/* get furthest element */
struct res_node *maxelem = rheap_get_max(heap);
/* and check if the new one is closer than that */
if(maxelem->dist_sq > dist_sq) {
rheap_remove_max(heap);
if(rheap_insert(heap, node, dist_sq) == -1) {
return -1;
}
added_res = 1;
range_sq = dist_sq;
}
} else {
if(rheap_insert(heap, node, dist_sq) == -1) {
return =1;
}
added_res = 1;
}
}
/* find signed distance from the splitting plane */
dx = pos[node->dir] - node->pos[node->dir];
ret = find_nearest_n(dx <= 0.0 ? node->left : node->right, pos, range, num, heap, dim);
if(ret >= 0 && fabs(dx) < range) {
added_res += ret;
ret = find_nearest_n(dx <= 0.0 ? node->right : node->left, pos, range, num, heap, dim);
}
}
#endif
static void kd_nearest_i(struct kdnode *node, const double *pos, struct kdnode **result, double *result_dist_sq, struct kdhyperrect* rect)
{
int dir = node->dir;
int i;
double dummy, dist_sq;
struct kdnode *nearer_subtree, *farther_subtree;
double *nearer_hyperrect_coord, *farther_hyperrect_coord;
/* Decide whether to go left or right in the tree */
dummy = pos[dir] - node->pos[dir];
if (dummy <= 0) {
nearer_subtree = node->left;
farther_subtree = node->right;
nearer_hyperrect_coord = rect->max + dir;
farther_hyperrect_coord = rect->min + dir;
} else {
nearer_subtree = node->right;
farther_subtree = node->left;
nearer_hyperrect_coord = rect->min + dir;
farther_hyperrect_coord = rect->max + dir;
}
if (nearer_subtree) {
/* Slice the hyperrect to get the hyperrect of the nearer subtree */
dummy = *nearer_hyperrect_coord;
*nearer_hyperrect_coord = node->pos[dir];
/* Recurse down into nearer subtree */
kd_nearest_i(nearer_subtree, pos, result, result_dist_sq, rect);
/* Undo the slice */
*nearer_hyperrect_coord = dummy;
}
/* Check the distance of the point at the current node, compare it
* with our best so far */
dist_sq = 0;
for(i=0; i < rect->dim; i++) {
dist_sq += SQ(node->pos[i] - pos[i]);
}
if (dist_sq < *result_dist_sq) {
*result = node;
*result_dist_sq = dist_sq;
}
if (farther_subtree) {
/* Get the hyperrect of the farther subtree */
dummy = *farther_hyperrect_coord;
*farther_hyperrect_coord = node->pos[dir];
/* Check if we have to recurse down by calculating the closest
* point of the hyperrect and see if it's closer than our
* minimum distance in result_dist_sq. */
if (hyperrect_dist_sq(rect, pos) < *result_dist_sq) {
/* Recurse down into farther subtree */
kd_nearest_i(farther_subtree, pos, result, result_dist_sq, rect);
}
/* Undo the slice on the hyperrect */
*farther_hyperrect_coord = dummy;
}
}
struct kdres *kd_nearest(struct kdtree *kd, const double *pos)
{
struct kdhyperrect *rect;
struct kdnode *result;
struct kdres *rset;
double dist_sq;
int i;
if (!kd) return 0;
if (!kd->rect) return 0;
/* Allocate result set */
if(!(rset = malloc(sizeof *rset))) {
return 0;
}
if(!(rset->rlist = alloc_resnode())) {
free(rset);
return 0;
}
rset->rlist->next = 0;
rset->tree = kd;
/* Duplicate the bounding hyperrectangle, we will work on the copy */
if (!(rect = hyperrect_duplicate(kd->rect))) {
kd_res_free(rset);
return 0;
}
/* Our first guesstimate is the root node */
result = kd->root;
dist_sq = 0;
for (i = 0; i < kd->dim; i++)
dist_sq += SQ(result->pos[i] - pos[i]);
/* Search for the nearest neighbour recursively */
kd_nearest_i(kd->root, pos, &result, &dist_sq, rect);
/* Free the copy of the hyperrect */
hyperrect_free(rect);
/* Store the result */
if (result) {
if (rlist_insert(rset->rlist, result, -1.0) == -1) {
kd_res_free(rset);
return 0;
}
rset->size = 1;
kd_res_rewind(rset);
return rset;
} else {
kd_res_free(rset);
return 0;
}
}
struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z)
{
double pos[3];
pos[0] = x;
pos[1] = y;
pos[2] = z;
return kd_nearest(tree, pos);
}
/* ---- nearest N search ---- */
/*
static kdres *kd_nearest_n(struct kdtree *kd, const double *pos, int num)
{
int ret;
struct kdres *rset;
if(!(rset = malloc(sizeof *rset))) {
return 0;
}
if(!(rset->rlist = alloc_resnode())) {
free(rset);
return 0;
}
rset->rlist->next = 0;
rset->tree = kd;
if((ret = find_nearest_n(kd->root, pos, range, num, rset->rlist, kd->dim)) == -1) {
kd_res_free(rset);
return 0;
}
rset->size = ret;
kd_res_rewind(rset);
return rset;
}*/
struct kdres *kd_nearest_range(struct kdtree *kd, const double *pos, double range)
{
int ret;
struct kdres *rset;
if(!(rset = malloc(sizeof *rset))) {
return 0;
}
if(!(rset->rlist = alloc_resnode())) {
free(rset);
return 0;
}
rset->rlist->next = 0;
rset->tree = kd;
if((ret = find_nearest(kd->root, pos, range, rset->rlist, 0, kd->dim)) == -1) {
kd_res_free(rset);
return 0;
}
rset->size = ret;
kd_res_rewind(rset);