nuclear@48: /*
nuclear@48: This file is part of ``kdtree'', a library for working with kd-trees.
nuclear@48: Copyright (C) 2007-2011 John Tsiombikas <nuclear@member.fsf.org>
nuclear@48: 
nuclear@48: Redistribution and use in source and binary forms, with or without
nuclear@48: modification, are permitted provided that the following conditions are met:
nuclear@48: 
nuclear@48: 1. Redistributions of source code must retain the above copyright notice, this
nuclear@48:    list of conditions and the following disclaimer.
nuclear@48: 2. Redistributions in binary form must reproduce the above copyright notice,
nuclear@48:    this list of conditions and the following disclaimer in the documentation
nuclear@48:    and/or other materials provided with the distribution.
nuclear@48: 3. The name of the author may not be used to endorse or promote products
nuclear@48:    derived from this software without specific prior written permission.
nuclear@48: 
nuclear@48: THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
nuclear@48: WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
nuclear@48: MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
nuclear@48: EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
nuclear@48: EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
nuclear@48: OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
nuclear@48: INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
nuclear@48: CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
nuclear@48: IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
nuclear@48: OF SUCH DAMAGE.
nuclear@48: */
nuclear@48: /* single nearest neighbor search written by Tamas Nepusz <tamas@cs.rhul.ac.uk> */
nuclear@48: #include <stdio.h>
nuclear@48: #include <stdlib.h>
nuclear@48: #include <string.h>
nuclear@48: #include <math.h>
nuclear@48: #include "kdtree.h"
nuclear@48: 
nuclear@48: #define USE_LIST_NODE_ALLOCATOR
nuclear@48: #define NO_PTHREADS
nuclear@48: #define I_WANT_THREAD_BUGS
nuclear@48: 
nuclear@48: #if defined(WIN32) || defined(__WIN32__)
nuclear@48: #include <malloc.h>
nuclear@48: #endif
nuclear@48: 
nuclear@48: #ifdef USE_LIST_NODE_ALLOCATOR
nuclear@48: 
nuclear@48: #ifndef NO_PTHREADS
nuclear@48: #include <pthread.h>
nuclear@48: #else
nuclear@48: 
nuclear@48: #ifndef I_WANT_THREAD_BUGS
nuclear@48: #error "You are compiling with the fast list node allocator, with pthreads disabled! This WILL break if used from multiple threads."
nuclear@48: #endif	/* I want thread bugs */
nuclear@48: 
nuclear@48: #endif	/* pthread support */
nuclear@48: #endif	/* use list node allocator */
nuclear@48: 
nuclear@48: struct kdhyperrect {
nuclear@48: 	int dim;
nuclear@48: 	double *min, *max;              /* minimum/maximum coords */
nuclear@48: };
nuclear@48: 
nuclear@48: struct kdnode {
nuclear@48: 	double *pos;
nuclear@48: 	int dir;
nuclear@48: 	void *data;
nuclear@48: 
nuclear@48: 	struct kdnode *left, *right;	/* negative/positive side */
nuclear@48: };
nuclear@48: 
nuclear@48: struct res_node {
nuclear@48: 	struct kdnode *item;
nuclear@48: 	double dist_sq;
nuclear@48: 	struct res_node *next;
nuclear@48: };
nuclear@48: 
nuclear@48: struct kdtree {
nuclear@48: 	int dim;
nuclear@48: 	struct kdnode *root;
nuclear@48: 	struct kdhyperrect *rect;
nuclear@48: 	void (*destr)(void*);
nuclear@48: };
nuclear@48: 
nuclear@48: struct kdres {
nuclear@48: 	struct kdtree *tree;
nuclear@48: 	struct res_node *rlist, *riter;
nuclear@48: 	int size;
nuclear@48: };
nuclear@48: 
nuclear@48: #define SQ(x)			((x) * (x))
nuclear@48: 
nuclear@48: 
nuclear@48: static void clear_rec(struct kdnode *node, void (*destr)(void*));
nuclear@48: static int insert_rec(struct kdnode **node, const double *pos, void *data, int dir, int dim);
nuclear@48: static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq);
nuclear@48: static void clear_results(struct kdres *set);
nuclear@48: 
nuclear@48: static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max);
nuclear@48: static void hyperrect_free(struct kdhyperrect *rect);
nuclear@48: static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect);
nuclear@48: static void hyperrect_extend(struct kdhyperrect *rect, const double *pos);
nuclear@48: static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos);
nuclear@48: 
nuclear@48: #ifdef USE_LIST_NODE_ALLOCATOR
nuclear@48: static struct res_node *alloc_resnode(void);
nuclear@48: static void free_resnode(struct res_node*);
nuclear@48: #else
nuclear@48: #define alloc_resnode()		malloc(sizeof(struct res_node))
nuclear@48: #define free_resnode(n)		free(n)
nuclear@48: #endif
nuclear@48: 
nuclear@48: 
nuclear@48: 
nuclear@48: struct kdtree *kd_create(int k)
nuclear@48: {
nuclear@48: 	struct kdtree *tree;
nuclear@48: 
nuclear@48: 	if(!(tree = malloc(sizeof *tree))) {
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	tree->dim = k;
nuclear@48: 	tree->root = 0;
nuclear@48: 	tree->destr = 0;
nuclear@48: 	tree->rect = 0;
nuclear@48: 
nuclear@48: 	return tree;
nuclear@48: }
nuclear@48: 
nuclear@48: void kd_free(struct kdtree *tree)
nuclear@48: {
nuclear@48: 	if(tree) {
nuclear@48: 		kd_clear(tree);
nuclear@48: 		free(tree);
nuclear@48: 	}
nuclear@48: }
nuclear@48: 
nuclear@48: static void clear_rec(struct kdnode *node, void (*destr)(void*))
nuclear@48: {
nuclear@48: 	if(!node) return;
nuclear@48: 
nuclear@48: 	clear_rec(node->left, destr);
nuclear@48: 	clear_rec(node->right, destr);
nuclear@48: 	
nuclear@48: 	if(destr) {
nuclear@48: 		destr(node->data);
nuclear@48: 	}
nuclear@48: 	free(node->pos);
nuclear@48: 	free(node);
nuclear@48: }
nuclear@48: 
nuclear@48: void kd_clear(struct kdtree *tree)
nuclear@48: {
nuclear@48: 	clear_rec(tree->root, tree->destr);
nuclear@48: 	tree->root = 0;
nuclear@48: 
nuclear@48: 	if (tree->rect) {
nuclear@48: 		hyperrect_free(tree->rect);
nuclear@48: 		tree->rect = 0;
nuclear@48: 	}
nuclear@48: }
nuclear@48: 
nuclear@48: void kd_data_destructor(struct kdtree *tree, void (*destr)(void*))
nuclear@48: {
nuclear@48: 	tree->destr = destr;
nuclear@48: }
nuclear@48: 
nuclear@48: 
nuclear@48: static int insert_rec(struct kdnode **nptr, const double *pos, void *data, int dir, int dim)
nuclear@48: {
nuclear@48: 	int new_dir;
nuclear@48: 	struct kdnode *node;
nuclear@48: 
nuclear@48: 	if(!*nptr) {
nuclear@48: 		if(!(node = malloc(sizeof *node))) {
nuclear@48: 			return -1;
nuclear@48: 		}
nuclear@48: 		if(!(node->pos = malloc(dim * sizeof *node->pos))) {
nuclear@48: 			free(node);
nuclear@48: 			return -1;
nuclear@48: 		}
nuclear@48: 		memcpy(node->pos, pos, dim * sizeof *node->pos);
nuclear@48: 		node->data = data;
nuclear@48: 		node->dir = dir;
nuclear@48: 		node->left = node->right = 0;
nuclear@48: 		*nptr = node;
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	node = *nptr;
nuclear@48: 	new_dir = (node->dir + 1) % dim;
nuclear@48: 	if(pos[node->dir] < node->pos[node->dir]) {
nuclear@48: 		return insert_rec(&(*nptr)->left, pos, data, new_dir, dim);
nuclear@48: 	}
nuclear@48: 	return insert_rec(&(*nptr)->right, pos, data, new_dir, dim);
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_insert(struct kdtree *tree, const double *pos, void *data)
nuclear@48: {
nuclear@48: 	if (insert_rec(&tree->root, pos, data, 0, tree->dim)) {
nuclear@48: 		return -1;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	if (tree->rect == 0) {
nuclear@48: 		tree->rect = hyperrect_create(tree->dim, pos, pos);
nuclear@48: 	} else {
nuclear@48: 		hyperrect_extend(tree->rect, pos);
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	return 0;
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_insertf(struct kdtree *tree, const float *pos, void *data)
nuclear@48: {
nuclear@48: 	static double sbuf[16];
nuclear@48: 	double *bptr, *buf = 0;
nuclear@48: 	int res, dim = tree->dim;
nuclear@48: 
nuclear@48: 	if(dim > 16) {
nuclear@48: #ifndef NO_ALLOCA
nuclear@48: 		if(dim <= 256)
nuclear@48: 			bptr = buf = alloca(dim * sizeof *bptr);
nuclear@48: 		else
nuclear@48: #endif
nuclear@48: 			if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
nuclear@48: 				return -1;
nuclear@48: 			}
nuclear@48: 	} else {
nuclear@48: 		bptr = buf = sbuf;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	while(dim-- > 0) {
nuclear@48: 		*bptr++ = *pos++;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	res = kd_insert(tree, buf, data);
nuclear@48: #ifndef NO_ALLOCA
nuclear@48: 	if(tree->dim > 256)
nuclear@48: #else
nuclear@48: 	if(tree->dim > 16)
nuclear@48: #endif
nuclear@48: 		free(buf);
nuclear@48: 	return res;
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data)
nuclear@48: {
nuclear@48: 	double buf[3];
nuclear@48: 	buf[0] = x;
nuclear@48: 	buf[1] = y;
nuclear@48: 	buf[2] = z;
nuclear@48: 	return kd_insert(tree, buf, data);
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_insert3f(struct kdtree *tree, float x, float y, float z, void *data)
nuclear@48: {
nuclear@48: 	double buf[3];
nuclear@48: 	buf[0] = x;
nuclear@48: 	buf[1] = y;
nuclear@48: 	buf[2] = z;
nuclear@48: 	return kd_insert(tree, buf, data);
nuclear@48: }
nuclear@48: 
nuclear@48: static int find_nearest(struct kdnode *node, const double *pos, double range, struct res_node *list, int ordered, int dim)
nuclear@48: {
nuclear@48: 	double dist_sq, dx;
nuclear@48: 	int i, ret, added_res = 0;
nuclear@48: 
nuclear@48: 	if(!node) return 0;
nuclear@48: 
nuclear@48: 	dist_sq = 0;
nuclear@48: 	for(i=0; i<dim; i++) {
nuclear@48: 		dist_sq += SQ(node->pos[i] - pos[i]);
nuclear@48: 	}
nuclear@48: 	if(dist_sq <= SQ(range)) {
nuclear@48: 		if(rlist_insert(list, node, ordered ? dist_sq : -1.0) == -1) {
nuclear@48: 			return -1;
nuclear@48: 		}
nuclear@48: 		added_res = 1;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	dx = pos[node->dir] - node->pos[node->dir];
nuclear@48: 
nuclear@48: 	ret = find_nearest(dx <= 0.0 ? node->left : node->right, pos, range, list, ordered, dim);
nuclear@48: 	if(ret >= 0 && fabs(dx) < range) {
nuclear@48: 		added_res += ret;
nuclear@48: 		ret = find_nearest(dx <= 0.0 ? node->right : node->left, pos, range, list, ordered, dim);
nuclear@48: 	}
nuclear@48: 	if(ret == -1) {
nuclear@48: 		return -1;
nuclear@48: 	}
nuclear@48: 	added_res += ret;
nuclear@48: 
nuclear@48: 	return added_res;
nuclear@48: }
nuclear@48: 
nuclear@48: #if 0
nuclear@48: static int find_nearest_n(struct kdnode *node, const double *pos, double range, int num, struct rheap *heap, int dim)
nuclear@48: {
nuclear@48: 	double dist_sq, dx;
nuclear@48: 	int i, ret, added_res = 0;
nuclear@48: 
nuclear@48: 	if(!node) return 0;
nuclear@48: 	
nuclear@48: 	/* if the photon is close enough, add it to the result heap */
nuclear@48: 	dist_sq = 0;
nuclear@48: 	for(i=0; i<dim; i++) {
nuclear@48: 		dist_sq += SQ(node->pos[i] - pos[i]);
nuclear@48: 	}
nuclear@48: 	if(dist_sq <= range_sq) {
nuclear@48: 		if(heap->size >= num) {
nuclear@48: 			/* get furthest element */
nuclear@48: 			struct res_node *maxelem = rheap_get_max(heap);
nuclear@48: 
nuclear@48: 			/* and check if the new one is closer than that */
nuclear@48: 			if(maxelem->dist_sq > dist_sq) {
nuclear@48: 				rheap_remove_max(heap);
nuclear@48: 
nuclear@48: 				if(rheap_insert(heap, node, dist_sq) == -1) {
nuclear@48: 					return -1;
nuclear@48: 				}
nuclear@48: 				added_res = 1;
nuclear@48: 
nuclear@48: 				range_sq = dist_sq;
nuclear@48: 			}
nuclear@48: 		} else {
nuclear@48: 			if(rheap_insert(heap, node, dist_sq) == -1) {
nuclear@48: 				return =1;
nuclear@48: 			}
nuclear@48: 			added_res = 1;
nuclear@48: 		}
nuclear@48: 	}
nuclear@48: 
nuclear@48: 
nuclear@48: 	/* find signed distance from the splitting plane */
nuclear@48: 	dx = pos[node->dir] - node->pos[node->dir];
nuclear@48: 
nuclear@48: 	ret = find_nearest_n(dx <= 0.0 ? node->left : node->right, pos, range, num, heap, dim);
nuclear@48: 	if(ret >= 0 && fabs(dx) < range) {
nuclear@48: 		added_res += ret;
nuclear@48: 		ret = find_nearest_n(dx <= 0.0 ? node->right : node->left, pos, range, num, heap, dim);
nuclear@48: 	}
nuclear@48: 
nuclear@48: }
nuclear@48: #endif
nuclear@48: 
nuclear@48: static void kd_nearest_i(struct kdnode *node, const double *pos, struct kdnode **result, double *result_dist_sq, struct kdhyperrect* rect)
nuclear@48: {
nuclear@48: 	int dir = node->dir;
nuclear@48: 	int i;
nuclear@48: 	double dummy, dist_sq;
nuclear@48: 	struct kdnode *nearer_subtree, *farther_subtree;
nuclear@48: 	double *nearer_hyperrect_coord, *farther_hyperrect_coord;
nuclear@48: 
nuclear@48: 	/* Decide whether to go left or right in the tree */
nuclear@48: 	dummy = pos[dir] - node->pos[dir];
nuclear@48: 	if (dummy <= 0) {
nuclear@48: 		nearer_subtree = node->left;
nuclear@48: 		farther_subtree = node->right;
nuclear@48: 		nearer_hyperrect_coord = rect->max + dir;
nuclear@48: 		farther_hyperrect_coord = rect->min + dir;
nuclear@48: 	} else {
nuclear@48: 		nearer_subtree = node->right;
nuclear@48: 		farther_subtree = node->left;
nuclear@48: 		nearer_hyperrect_coord = rect->min + dir;
nuclear@48: 		farther_hyperrect_coord = rect->max + dir;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	if (nearer_subtree) {
nuclear@48: 		/* Slice the hyperrect to get the hyperrect of the nearer subtree */
nuclear@48: 		dummy = *nearer_hyperrect_coord;
nuclear@48: 		*nearer_hyperrect_coord = node->pos[dir];
nuclear@48: 		/* Recurse down into nearer subtree */
nuclear@48: 		kd_nearest_i(nearer_subtree, pos, result, result_dist_sq, rect);
nuclear@48: 		/* Undo the slice */
nuclear@48: 		*nearer_hyperrect_coord = dummy;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	/* Check the distance of the point at the current node, compare it
nuclear@48: 	 * with our best so far */
nuclear@48: 	dist_sq = 0;
nuclear@48: 	for(i=0; i < rect->dim; i++) {
nuclear@48: 		dist_sq += SQ(node->pos[i] - pos[i]);
nuclear@48: 	}
nuclear@48: 	if (dist_sq < *result_dist_sq) {
nuclear@48: 		*result = node;
nuclear@48: 		*result_dist_sq = dist_sq;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	if (farther_subtree) {
nuclear@48: 		/* Get the hyperrect of the farther subtree */
nuclear@48: 		dummy = *farther_hyperrect_coord;
nuclear@48: 		*farther_hyperrect_coord = node->pos[dir];
nuclear@48: 		/* Check if we have to recurse down by calculating the closest
nuclear@48: 		 * point of the hyperrect and see if it's closer than our
nuclear@48: 		 * minimum distance in result_dist_sq. */
nuclear@48: 		if (hyperrect_dist_sq(rect, pos) < *result_dist_sq) {
nuclear@48: 			/* Recurse down into farther subtree */
nuclear@48: 			kd_nearest_i(farther_subtree, pos, result, result_dist_sq, rect);
nuclear@48: 		}
nuclear@48: 		/* Undo the slice on the hyperrect */
nuclear@48: 		*farther_hyperrect_coord = dummy;
nuclear@48: 	}
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearest(struct kdtree *kd, const double *pos)
nuclear@48: {
nuclear@48: 	struct kdhyperrect *rect;
nuclear@48: 	struct kdnode *result;
nuclear@48: 	struct kdres *rset;
nuclear@48: 	double dist_sq;
nuclear@48: 	int i;
nuclear@48: 
nuclear@48: 	if (!kd) return 0;
nuclear@48: 	if (!kd->rect) return 0;
nuclear@48: 
nuclear@48: 	/* Allocate result set */
nuclear@48: 	if(!(rset = malloc(sizeof *rset))) {
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	if(!(rset->rlist = alloc_resnode())) {
nuclear@48: 		free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	rset->rlist->next = 0;
nuclear@48: 	rset->tree = kd;
nuclear@48: 
nuclear@48: 	/* Duplicate the bounding hyperrectangle, we will work on the copy */
nuclear@48: 	if (!(rect = hyperrect_duplicate(kd->rect))) {
nuclear@48: 		kd_res_free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	/* Our first guesstimate is the root node */
nuclear@48: 	result = kd->root;
nuclear@48: 	dist_sq = 0;
nuclear@48: 	for (i = 0; i < kd->dim; i++)
nuclear@48: 		dist_sq += SQ(result->pos[i] - pos[i]);
nuclear@48: 
nuclear@48: 	/* Search for the nearest neighbour recursively */
nuclear@48: 	kd_nearest_i(kd->root, pos, &result, &dist_sq, rect);
nuclear@48: 
nuclear@48: 	/* Free the copy of the hyperrect */
nuclear@48: 	hyperrect_free(rect);
nuclear@48: 
nuclear@48: 	/* Store the result */
nuclear@48: 	if (result) {
nuclear@48: 		if (rlist_insert(rset->rlist, result, -1.0) == -1) {
nuclear@48: 			kd_res_free(rset);
nuclear@48: 			return 0;
nuclear@48: 		}
nuclear@48: 		rset->size = 1;
nuclear@48: 		kd_res_rewind(rset);
nuclear@48: 		return rset;
nuclear@48: 	} else {
nuclear@48: 		kd_res_free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearestf(struct kdtree *tree, const float *pos)
nuclear@48: {
nuclear@48: 	static double sbuf[16];
nuclear@48: 	double *bptr, *buf = 0;
nuclear@48: 	int dim = tree->dim;
nuclear@48: 	struct kdres *res;
nuclear@48: 
nuclear@48: 	if(dim > 16) {
nuclear@48: #ifndef NO_ALLOCA
nuclear@48: 		if(dim <= 256)
nuclear@48: 			bptr = buf = alloca(dim * sizeof *bptr);
nuclear@48: 		else
nuclear@48: #endif
nuclear@48: 			if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
nuclear@48: 				return 0;
nuclear@48: 			}
nuclear@48: 	} else {
nuclear@48: 		bptr = buf = sbuf;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	while(dim-- > 0) {
nuclear@48: 		*bptr++ = *pos++;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	res = kd_nearest(tree, buf);
nuclear@48: #ifndef NO_ALLOCA
nuclear@48: 	if(tree->dim > 256)
nuclear@48: #else
nuclear@48: 	if(tree->dim > 16)
nuclear@48: #endif
nuclear@48: 		free(buf);
nuclear@48: 	return res;
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z)
nuclear@48: {
nuclear@48: 	double pos[3];
nuclear@48: 	pos[0] = x;
nuclear@48: 	pos[1] = y;
nuclear@48: 	pos[2] = z;
nuclear@48: 	return kd_nearest(tree, pos);
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z)
nuclear@48: {
nuclear@48: 	double pos[3];
nuclear@48: 	pos[0] = x;
nuclear@48: 	pos[1] = y;
nuclear@48: 	pos[2] = z;
nuclear@48: 	return kd_nearest(tree, pos);
nuclear@48: }
nuclear@48: 
nuclear@48: /* ---- nearest N search ---- */
nuclear@48: /*
nuclear@48: static kdres *kd_nearest_n(struct kdtree *kd, const double *pos, int num)
nuclear@48: {
nuclear@48: 	int ret;
nuclear@48: 	struct kdres *rset;
nuclear@48: 
nuclear@48: 	if(!(rset = malloc(sizeof *rset))) {
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	if(!(rset->rlist = alloc_resnode())) {
nuclear@48: 		free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	rset->rlist->next = 0;
nuclear@48: 	rset->tree = kd;
nuclear@48: 
nuclear@48: 	if((ret = find_nearest_n(kd->root, pos, range, num, rset->rlist, kd->dim)) == -1) {
nuclear@48: 		kd_res_free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	rset->size = ret;
nuclear@48: 	kd_res_rewind(rset);
nuclear@48: 	return rset;
nuclear@48: }*/
nuclear@48: 
nuclear@48: struct kdres *kd_nearest_range(struct kdtree *kd, const double *pos, double range)
nuclear@48: {
nuclear@48: 	int ret;
nuclear@48: 	struct kdres *rset;
nuclear@48: 
nuclear@48: 	if(!(rset = malloc(sizeof *rset))) {
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	if(!(rset->rlist = alloc_resnode())) {
nuclear@48: 		free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	rset->rlist->next = 0;
nuclear@48: 	rset->tree = kd;
nuclear@48: 
nuclear@48: 	if((ret = find_nearest(kd->root, pos, range, rset->rlist, 0, kd->dim)) == -1) {
nuclear@48: 		kd_res_free(rset);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	rset->size = ret;
nuclear@48: 	kd_res_rewind(rset);
nuclear@48: 	return rset;
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearest_rangef(struct kdtree *kd, const float *pos, float range)
nuclear@48: {
nuclear@48: 	static double sbuf[16];
nuclear@48: 	double *bptr, *buf = 0;
nuclear@48: 	int dim = kd->dim;
nuclear@48: 	struct kdres *res;
nuclear@48: 
nuclear@48: 	if(dim > 16) {
nuclear@48: #ifndef NO_ALLOCA
nuclear@48: 		if(dim <= 256)
nuclear@48: 			bptr = buf = alloca(dim * sizeof *bptr);
nuclear@48: 		else
nuclear@48: #endif
nuclear@48: 			if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
nuclear@48: 				return 0;
nuclear@48: 			}
nuclear@48: 	} else {
nuclear@48: 		bptr = buf = sbuf;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	while(dim-- > 0) {
nuclear@48: 		*bptr++ = *pos++;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	res = kd_nearest_range(kd, buf, range);
nuclear@48: #ifndef NO_ALLOCA
nuclear@48: 	if(kd->dim > 256)
nuclear@48: #else
nuclear@48: 	if(kd->dim > 16)
nuclear@48: #endif
nuclear@48: 		free(buf);
nuclear@48: 	return res;
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double z, double range)
nuclear@48: {
nuclear@48: 	double buf[3];
nuclear@48: 	buf[0] = x;
nuclear@48: 	buf[1] = y;
nuclear@48: 	buf[2] = z;
nuclear@48: 	return kd_nearest_range(tree, buf, range);
nuclear@48: }
nuclear@48: 
nuclear@48: struct kdres *kd_nearest_range3f(struct kdtree *tree, float x, float y, float z, float range)
nuclear@48: {
nuclear@48: 	double buf[3];
nuclear@48: 	buf[0] = x;
nuclear@48: 	buf[1] = y;
nuclear@48: 	buf[2] = z;
nuclear@48: 	return kd_nearest_range(tree, buf, range);
nuclear@48: }
nuclear@48: 
nuclear@48: void kd_res_free(struct kdres *rset)
nuclear@48: {
nuclear@48: 	clear_results(rset);
nuclear@48: 	free_resnode(rset->rlist);
nuclear@48: 	free(rset);
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_res_size(struct kdres *set)
nuclear@48: {
nuclear@48: 	return (set->size);
nuclear@48: }
nuclear@48: 
nuclear@48: void kd_res_rewind(struct kdres *rset)
nuclear@48: {
nuclear@48: 	rset->riter = rset->rlist->next;
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_res_end(struct kdres *rset)
nuclear@48: {
nuclear@48: 	return rset->riter == 0;
nuclear@48: }
nuclear@48: 
nuclear@48: int kd_res_next(struct kdres *rset)
nuclear@48: {
nuclear@48: 	rset->riter = rset->riter->next;
nuclear@48: 	return rset->riter != 0;
nuclear@48: }
nuclear@48: 
nuclear@48: void *kd_res_item(struct kdres *rset, double *pos)
nuclear@48: {
nuclear@48: 	if(rset->riter) {
nuclear@48: 		if(pos) {
nuclear@48: 			memcpy(pos, rset->riter->item->pos, rset->tree->dim * sizeof *pos);
nuclear@48: 		}
nuclear@48: 		return rset->riter->item->data;
nuclear@48: 	}
nuclear@48: 	return 0;
nuclear@48: }
nuclear@48: 
nuclear@48: void *kd_res_itemf(struct kdres *rset, float *pos)
nuclear@48: {
nuclear@48: 	if(rset->riter) {
nuclear@48: 		if(pos) {
nuclear@48: 			int i;
nuclear@48: 			for(i=0; i<rset->tree->dim; i++) {
nuclear@48: 				pos[i] = rset->riter->item->pos[i];
nuclear@48: 			}
nuclear@48: 		}
nuclear@48: 		return rset->riter->item->data;
nuclear@48: 	}
nuclear@48: 	return 0;
nuclear@48: }
nuclear@48: 
nuclear@48: void *kd_res_item3(struct kdres *rset, double *x, double *y, double *z)
nuclear@48: {
nuclear@48: 	if(rset->riter) {
nuclear@48: 		if(*x) *x = rset->riter->item->pos[0];
nuclear@48: 		if(*y) *y = rset->riter->item->pos[1];
nuclear@48: 		if(*z) *z = rset->riter->item->pos[2];
nuclear@48: 	}
nuclear@48: 	return 0;
nuclear@48: }
nuclear@48: 
nuclear@48: void *kd_res_item3f(struct kdres *rset, float *x, float *y, float *z)
nuclear@48: {
nuclear@48: 	if(rset->riter) {
nuclear@48: 		if(*x) *x = rset->riter->item->pos[0];
nuclear@48: 		if(*y) *y = rset->riter->item->pos[1];
nuclear@48: 		if(*z) *z = rset->riter->item->pos[2];
nuclear@48: 	}
nuclear@48: 	return 0;
nuclear@48: }
nuclear@48: 
nuclear@48: void *kd_res_item_data(struct kdres *set)
nuclear@48: {
nuclear@48: 	return kd_res_item(set, 0);
nuclear@48: }
nuclear@48: 
nuclear@48: /* ---- hyperrectangle helpers ---- */
nuclear@48: static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max)
nuclear@48: {
nuclear@48: 	size_t size = dim * sizeof(double);
nuclear@48: 	struct kdhyperrect* rect = 0;
nuclear@48: 
nuclear@48: 	if (!(rect = malloc(sizeof(struct kdhyperrect)))) {
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	rect->dim = dim;
nuclear@48: 	if (!(rect->min = malloc(size))) {
nuclear@48: 		free(rect);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	if (!(rect->max = malloc(size))) {
nuclear@48: 		free(rect->min);
nuclear@48: 		free(rect);
nuclear@48: 		return 0;
nuclear@48: 	}
nuclear@48: 	memcpy(rect->min, min, size);
nuclear@48: 	memcpy(rect->max, max, size);
nuclear@48: 
nuclear@48: 	return rect;
nuclear@48: }
nuclear@48: 
nuclear@48: static void hyperrect_free(struct kdhyperrect *rect)
nuclear@48: {
nuclear@48: 	free(rect->min);
nuclear@48: 	free(rect->max);
nuclear@48: 	free(rect);
nuclear@48: }
nuclear@48: 
nuclear@48: static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect)
nuclear@48: {
nuclear@48: 	return hyperrect_create(rect->dim, rect->min, rect->max);
nuclear@48: }
nuclear@48: 
nuclear@48: static void hyperrect_extend(struct kdhyperrect *rect, const double *pos)
nuclear@48: {
nuclear@48: 	int i;
nuclear@48: 
nuclear@48: 	for (i=0; i < rect->dim; i++) {
nuclear@48: 		if (pos[i] < rect->min[i]) {
nuclear@48: 			rect->min[i] = pos[i];
nuclear@48: 		}
nuclear@48: 		if (pos[i] > rect->max[i]) {
nuclear@48: 			rect->max[i] = pos[i];
nuclear@48: 		}
nuclear@48: 	}
nuclear@48: }
nuclear@48: 
nuclear@48: static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos)
nuclear@48: {
nuclear@48: 	int i;
nuclear@48: 	double result = 0;
nuclear@48: 
nuclear@48: 	for (i=0; i < rect->dim; i++) {
nuclear@48: 		if (pos[i] < rect->min[i]) {
nuclear@48: 			result += SQ(rect->min[i] - pos[i]);
nuclear@48: 		} else if (pos[i] > rect->max[i]) {
nuclear@48: 			result += SQ(rect->max[i] - pos[i]);
nuclear@48: 		}
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	return result;
nuclear@48: }
nuclear@48: 
nuclear@48: /* ---- static helpers ---- */
nuclear@48: 
nuclear@48: #ifdef USE_LIST_NODE_ALLOCATOR
nuclear@48: /* special list node allocators. */
nuclear@48: static struct res_node *free_nodes;
nuclear@48: 
nuclear@48: #ifndef NO_PTHREADS
nuclear@48: static pthread_mutex_t alloc_mutex = PTHREAD_MUTEX_INITIALIZER;
nuclear@48: #endif
nuclear@48: 
nuclear@48: static struct res_node *alloc_resnode(void)
nuclear@48: {
nuclear@48: 	struct res_node *node;
nuclear@48: 
nuclear@48: #ifndef NO_PTHREADS
nuclear@48: 	pthread_mutex_lock(&alloc_mutex);
nuclear@48: #endif
nuclear@48: 
nuclear@48: 	if(!free_nodes) {
nuclear@48: 		node = malloc(sizeof *node);
nuclear@48: 	} else {
nuclear@48: 		node = free_nodes;
nuclear@48: 		free_nodes = free_nodes->next;
nuclear@48: 		node->next = 0;
nuclear@48: 	}
nuclear@48: 
nuclear@48: #ifndef NO_PTHREADS
nuclear@48: 	pthread_mutex_unlock(&alloc_mutex);
nuclear@48: #endif
nuclear@48: 
nuclear@48: 	return node;
nuclear@48: }
nuclear@48: 
nuclear@48: static void free_resnode(struct res_node *node)
nuclear@48: {
nuclear@48: #ifndef NO_PTHREADS
nuclear@48: 	pthread_mutex_lock(&alloc_mutex);
nuclear@48: #endif
nuclear@48: 
nuclear@48: 	node->next = free_nodes;
nuclear@48: 	free_nodes = node;
nuclear@48: 
nuclear@48: #ifndef NO_PTHREADS
nuclear@48: 	pthread_mutex_unlock(&alloc_mutex);
nuclear@48: #endif
nuclear@48: }
nuclear@48: #endif	/* list node allocator or not */
nuclear@48: 
nuclear@48: 
nuclear@48: /* inserts the item. if dist_sq is >= 0, then do an ordered insert */
nuclear@48: /* TODO make the ordering code use heapsort */
nuclear@48: static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq)
nuclear@48: {
nuclear@48: 	struct res_node *rnode;
nuclear@48: 
nuclear@48: 	if(!(rnode = alloc_resnode())) {
nuclear@48: 		return -1;
nuclear@48: 	}
nuclear@48: 	rnode->item = item;
nuclear@48: 	rnode->dist_sq = dist_sq;
nuclear@48: 
nuclear@48: 	if(dist_sq >= 0.0) {
nuclear@48: 		while(list->next && list->next->dist_sq < dist_sq) {
nuclear@48: 			list = list->next;
nuclear@48: 		}
nuclear@48: 	}
nuclear@48: 	rnode->next = list->next;
nuclear@48: 	list->next = rnode;
nuclear@48: 	return 0;
nuclear@48: }
nuclear@48: 
nuclear@48: static void clear_results(struct kdres *rset)
nuclear@48: {
nuclear@48: 	struct res_node *tmp, *node = rset->rlist->next;
nuclear@48: 
nuclear@48: 	while(node) {
nuclear@48: 		tmp = node;
nuclear@48: 		node = node->next;
nuclear@48: 		free_resnode(tmp);
nuclear@48: 	}
nuclear@48: 
nuclear@48: 	rset->rlist->next = 0;
nuclear@48: }