1 /* -copyright-
2 #-#
3 #-# xsnow: let it snow on your desktop
4 #-# Copyright (C) 1984,1988,1990,1993-1995,2000-2001 Rick Jansen
5 #-# 2019,2020,2021 Willem Vermin
6 #-#
7 #-# This program is free software: you can redistribute it and/or modify
8 #-# it under the terms of the GNU General Public License as published by
9 #-# the Free Software Foundation, either version 3 of the License, or
10 #-# (at your option) any later version.
11 #-#
12 #-# This program is distributed in the hope that it will be useful,
13 #-# but WITHOUT ANY WARRANTY; without even the implied warranty of
14 #-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 #-# GNU General Public License for more details.
16 #-#
17 #-# You should have received a copy of the GNU General Public License
18 #-# along with this program. If not, see <http://www.gnu.org/licenses/>.
19 #-#
20 */
21 /*
22 This file is part of ``kdtree'', a library for working with kd-trees.
23 Copyright (C) 2007-2011 John Tsiombikas <nuclear@member.fsf.org>
24
25 Redistribution and use in source and binary forms, with or without
26 modification, are permitted provided that the following conditions are met:
27
28 1. Redistributions of source code must retain the above copyright notice, this
29 list of conditions and the following disclaimer.
30 2. Redistributions in binary form must reproduce the above copyright notice,
31 this list of conditions and the following disclaimer in the documentation
32 and/or other materials provided with the distribution.
33 3. The name of the author may not be used to endorse or promote products
34 derived from this software without specific prior written permission.
35
36 THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
37 WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
38 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
39 EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
40 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
41 OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
42 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
43 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
44 IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
45 OF SUCH DAMAGE.
46 */
47 /* single nearest neighbor search written by Tamas Nepusz <tamas@cs.rhul.ac.uk> */
48
49 #define USE_LIST_NODE_ALLOCATOR
50 #define NO_ALLOCA
51
52 #ifdef HAVE_CONFIG_H
53 #include <config.h>
54 #endif
55
56 #include <stdio.h>
57 #include <stdlib.h>
58 #ifndef NO_ALLOCA
59 #ifdef HAVE_ALLOCA_H
60 #include <alloca.h>
61 #endif
62 #endif
63 #include <string.h>
64 #include <math.h>
65 #include <assert.h>
66 #include "kdtree.h"
67
68 #if defined(WIN32) || defined(__WIN32__)
69 #include <malloc.h>
70 #endif
71
72 #ifdef USE_LIST_NODE_ALLOCATOR
73
74 #ifndef NO_PTHREADS
75 #include <pthread.h>
76 #else
77
78 #ifndef I_WANT_THREAD_BUGS
79 #error "You are compiling with the fast list node allocator, with pthreads disabled! This WILL break if used from multiple threads."
80 #endif /* I want thread bugs */
81
82 #endif /* pthread support */
83 #endif /* use list node allocator */
84
85 struct kdhyperrect {
86 int dim;
87 double *min, *max; /* minimum/maximum coords */
88 };
89
90 struct kdnode {
91 double *pos;
92 int dir;
93 void *data;
94
95 struct kdnode *left, *right; /* negative/positive side */
96 };
97
98 struct res_node {
99 struct kdnode *item;
100 double dist_sq;
101 struct res_node *next;
102 };
103
104 struct kdtree {
105 int dim;
106 struct kdnode *root;
107 struct kdhyperrect *rect;
108 void (*destr)(void*);
109 };
110
111 struct kdres {
112 struct kdtree *tree;
113 struct res_node *rlist, *riter;
114 int size;
115 };
116
117 #define SQ(x) ((x) * (x))
118
119
120 static void clear_rec(struct kdnode *node, void (*destr)(void*));
121 static int insert_rec(struct kdnode **node, const double *pos, void *data, int dir, int dim);
122 static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq);
123 static void clear_results(struct kdres *set);
124
125 static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max);
126 static void hyperrect_free(struct kdhyperrect *rect);
127 static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect);
128 static void hyperrect_extend(struct kdhyperrect *rect, const double *pos);
129 static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos);
130
131 #ifdef USE_LIST_NODE_ALLOCATOR
132 static struct res_node *alloc_resnode(void);
133 static void free_resnode(struct res_node*);
134 #else
135 #define alloc_resnode() (res_node*) malloc(sizeof(struct res_node))
136 #define free_resnode(n) free(n)
137 #endif
138
139
kd_create(int k)140 struct kdtree *kd_create(int k)
141 {
142 struct kdtree *tree;
143
144 if(!(tree = (struct kdtree *)malloc(sizeof *tree))) {
145 return NULL;
146 }
147
148 tree->dim = k;
149 tree->root = NULL;
150 tree->destr = NULL;
151 tree->rect = NULL;
152
153 return tree;
154 }
155
kd_free(struct kdtree * tree)156 void kd_free(struct kdtree *tree)
157 {
158 if(tree) {
159 kd_clear(tree);
160 free(tree);
161 }
162 }
163
clear_rec(struct kdnode * node,void (* destr)(void *))164 static void clear_rec(struct kdnode *node, void (*destr)(void*))
165 {
166 if(!node) return;
167
168 clear_rec(node->left, destr);
169 clear_rec(node->right, destr);
170
171 if(destr) {
172 destr(node->data);
173 }
174 free(node->pos);
175 free(node);
176 }
177
kd_clear(struct kdtree * tree)178 void kd_clear(struct kdtree *tree)
179 {
180 clear_rec(tree->root, tree->destr);
181 tree->root = NULL;
182
183 if (tree->rect) {
184 hyperrect_free(tree->rect);
185 tree->rect = NULL;
186 }
187 }
188
kd_data_destructor(struct kdtree * tree,void (* destr)(void *))189 void kd_data_destructor(struct kdtree *tree, void (*destr)(void*))
190 {
191 tree->destr = destr;
192 }
193
194
insert_rec(struct kdnode ** nptr,const double * pos,void * data,int dir,int dim)195 static int insert_rec(struct kdnode **nptr, const double *pos, void *data, int dir, int dim)
196 {
197 int new_dir;
198 struct kdnode *node;
199
200 if(!*nptr) {
201 if(!(node = (struct kdnode *)malloc(sizeof *node))) {
202 return -1;
203 }
204 assert(dim>0);
205 if(!(node->pos = (double *)malloc(dim * sizeof *node->pos))) {
206 free(node);
207 return -1;
208 }
209 memcpy(node->pos, pos, dim * sizeof *node->pos);
210 node->data = data;
211 node->dir = dir;
212 node->left = node->right = NULL;
213 *nptr = node;
214 return 0;
215 }
216
217 node = *nptr;
218 new_dir = (node->dir + 1) % dim;
219 if(pos[node->dir] < node->pos[node->dir]) {
220 return insert_rec(&(*nptr)->left, pos, data, new_dir, dim);
221 }
222 return insert_rec(&(*nptr)->right, pos, data, new_dir, dim);
223 }
224
kd_insert(struct kdtree * tree,const double * pos,void * data)225 int kd_insert(struct kdtree *tree, const double *pos, void *data)
226 {
227 if (insert_rec(&tree->root, pos, data, 0, tree->dim)) {
228 return -1;
229 }
230
231 if (tree->rect == NULL) {
232 tree->rect = hyperrect_create(tree->dim, pos, pos);
233 } else {
234 hyperrect_extend(tree->rect, pos);
235 }
236
237 return 0;
238 }
239
kd_insertf(struct kdtree * tree,const float * pos,void * data)240 int kd_insertf(struct kdtree *tree, const float *pos, void *data)
241 {
242 static double sbuf[16];
243 double *bptr, *buf = NULL;
244 int res, dim = tree->dim;
245 assert(dim>0);
246
247 if(dim > 16) {
248 #ifndef NO_ALLOCA
249 if(dim <= 256)
250 bptr = buf = (double *)alloca(dim * sizeof *bptr);
251 else
252 #endif
253 if(!(bptr = buf = (double *)malloc(dim * sizeof *bptr))) {
254 return -1;
255 }
256 } else {
257 bptr = buf = sbuf;
258 }
259
260 while(dim-- > 0) {
261 *bptr++ = *pos++;
262 }
263
264 res = kd_insert(tree, buf, data);
265 #ifndef NO_ALLOCA
266 if(tree->dim > 256)
267 #else
268 if(tree->dim > 16)
269 #endif
270 free(buf);
271 return res;
272 }
273
kd_insert3(struct kdtree * tree,double x,double y,double z,void * data)274 int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data)
275 {
276 double buf[3];
277 buf[0] = x;
278 buf[1] = y;
279 buf[2] = z;
280 return kd_insert(tree, buf, data);
281 }
282
kd_insert3f(struct kdtree * tree,float x,float y,float z,void * data)283 int kd_insert3f(struct kdtree *tree, float x, float y, float z, void *data)
284 {
285 double buf[3];
286 buf[0] = x;
287 buf[1] = y;
288 buf[2] = z;
289 return kd_insert(tree, buf, data);
290 }
291
find_nearest(struct kdnode * node,const double * pos,double range,struct res_node * list,int ordered,int dim)292 static int find_nearest(struct kdnode *node, const double *pos, double range, struct res_node *list, int ordered, int dim)
293 {
294 double dist_sq, dx;
295 int i, ret, added_res = 0;
296
297 if(!node) return 0;
298
299 dist_sq = 0;
300 for(i=0; i<dim; i++) {
301 dist_sq += SQ(node->pos[i] - pos[i]);
302 }
303 if(dist_sq <= SQ(range)) {
304 if(rlist_insert(list, node, ordered ? dist_sq : -1.0) == -1) {
305 return -1;
306 }
307 added_res = 1;
308 }
309
310 dx = pos[node->dir] - node->pos[node->dir];
311
312 ret = find_nearest(dx <= 0.0 ? node->left : node->right, pos, range, list, ordered, dim);
313 if(ret >= 0 && fabs(dx) < range) {
314 added_res += ret;
315 ret = find_nearest(dx <= 0.0 ? node->right : node->left, pos, range, list, ordered, dim);
316 }
317 if(ret == -1) {
318 return -1;
319 }
320 added_res += ret;
321
322 return added_res;
323 }
324
325 #if 0
326 static int find_nearest_n(struct kdnode *node, const double *pos, double range, int num, struct rheap *heap, int dim)
327 {
328 double dist_sq, dx;
329 int i, ret, added_res = 0;
330
331 if(!node) return 0;
332
333 /* if the photon is close enough, add it to the result heap */
334 dist_sq = 0;
335 for(i=0; i<dim; i++) {
336 dist_sq += SQ(node->pos[i] - pos[i]);
337 }
338 if(dist_sq <= range_sq) {
339 if(heap->size >= num) {
340 /* get furthest element */
341 struct res_node *maxelem = rheap_get_max(heap);
342
343 /* and check if the new one is closer than that */
344 if(maxelem->dist_sq > dist_sq) {
345 rheap_remove_max(heap);
346
347 if(rheap_insert(heap, node, dist_sq) == -1) {
348 return -1;
349 }
350 added_res = 1;
351
352 range_sq = dist_sq;
353 }
354 } else {
355 if(rheap_insert(heap, node, dist_sq) == -1) {
356 return =1;
357 }
358 added_res = 1;
359 }
360 }
361
362
363 /* find signed distance from the splitting plane */
364 dx = pos[node->dir] - node->pos[node->dir];
365
366 ret = find_nearest_n(dx <= 0.0 ? node->left : node->right, pos, range, num, heap, dim);
367 if(ret >= 0 && fabs(dx) < range) {
368 added_res += ret;
369 ret = find_nearest_n(dx <= 0.0 ? node->right : node->left, pos, range, num, heap, dim);
370 }
371
372 }
373 #endif
374
kd_nearest_i(struct kdnode * node,const double * pos,struct kdnode ** result,double * result_dist_sq,struct kdhyperrect * rect)375 static void kd_nearest_i(struct kdnode *node, const double *pos, struct kdnode **result, double *result_dist_sq, struct kdhyperrect* rect)
376 {
377 int dir = node->dir;
378 int i;
379 double dummy, dist_sq;
380 struct kdnode *nearer_subtree, *farther_subtree;
381 double *nearer_hyperrect_coord, *farther_hyperrect_coord;
382
383 /* Decide whether to go left or right in the tree */
384 dummy = pos[dir] - node->pos[dir];
385 if (dummy <= 0) {
386 nearer_subtree = node->left;
387 farther_subtree = node->right;
388 nearer_hyperrect_coord = rect->max + dir;
389 farther_hyperrect_coord = rect->min + dir;
390 } else {
391 nearer_subtree = node->right;
392 farther_subtree = node->left;
393 nearer_hyperrect_coord = rect->min + dir;
394 farther_hyperrect_coord = rect->max + dir;
395 }
396
397 if (nearer_subtree) {
398 /* Slice the hyperrect to get the hyperrect of the nearer subtree */
399 dummy = *nearer_hyperrect_coord;
400 *nearer_hyperrect_coord = node->pos[dir];
401 /* Recurse down into nearer subtree */
402 kd_nearest_i(nearer_subtree, pos, result, result_dist_sq, rect);
403 /* Undo the slice */
404 *nearer_hyperrect_coord = dummy;
405 }
406
407 /* Check the distance of the point at the current node, compare it
408 * with our best so far */
409 dist_sq = 0;
410 for(i=0; i < rect->dim; i++) {
411 dist_sq += SQ(node->pos[i] - pos[i]);
412 }
413 if (dist_sq < *result_dist_sq) {
414 *result = node;
415 *result_dist_sq = dist_sq;
416 }
417
418 if (farther_subtree) {
419 /* Get the hyperrect of the farther subtree */
420 dummy = *farther_hyperrect_coord;
421 *farther_hyperrect_coord = node->pos[dir];
422 /* Check if we have to recurse down by calculating the closest
423 * point of the hyperrect and see if it's closer than our
424 * minimum distance in result_dist_sq. */
425 if (hyperrect_dist_sq(rect, pos) < *result_dist_sq) {
426 /* Recurse down into farther subtree */
427 kd_nearest_i(farther_subtree, pos, result, result_dist_sq, rect);
428 }
429 /* Undo the slice on the hyperrect */
430 *farther_hyperrect_coord = dummy;
431 }
432 }
433
kd_nearest(struct kdtree * kd,const double * pos)434 struct kdres *kd_nearest(struct kdtree *kd, const double *pos)
435 {
436 struct kdhyperrect *rect;
437 struct kdnode *result;
438 struct kdres *rset;
439 double dist_sq;
440 int i;
441
442 if (!kd) return NULL;
443 if (!kd->rect) return NULL;
444
445 /* Allocate result set */
446 if(!(rset = (struct kdres *)malloc(sizeof *rset))) {
447 return NULL;
448 }
449 if(!(rset->rlist = alloc_resnode())) {
450 free(rset);
451 return NULL;
452 }
453 rset->rlist->next = NULL;
454 rset->tree = kd;
455
456 /* Duplicate the bounding hyperrectangle, we will work on the copy */
457 if (!(rect = hyperrect_duplicate(kd->rect))) {
458 kd_res_free(rset);
459 return NULL;
460 }
461
462 /* Our first guesstimate is the root node */
463 result = kd->root;
464 dist_sq = 0;
465 for (i = 0; i < kd->dim; i++)
466 dist_sq += SQ(result->pos[i] - pos[i]);
467
468 /* Search for the nearest neighbour recursively */
469 kd_nearest_i(kd->root, pos, &result, &dist_sq, rect);
470
471 /* Free the copy of the hyperrect */
472 hyperrect_free(rect);
473
474 /* Store the result */
475 if (result) {
476 if (rlist_insert(rset->rlist, result, -1.0) == -1) {
477 kd_res_free(rset);
478 return NULL;
479 }
480 rset->size = 1;
481 kd_res_rewind(rset);
482 return rset;
483 } else {
484 kd_res_free(rset);
485 return NULL;
486 }
487 }
488
kd_nearestf(struct kdtree * tree,const float * pos)489 struct kdres *kd_nearestf(struct kdtree *tree, const float *pos)
490 {
491 static double sbuf[16];
492 double *bptr, *buf = NULL;
493 int dim = tree->dim;
494 struct kdres *res;
495 assert(dim>0);
496
497 if(dim > 16) {
498 #ifndef NO_ALLOCA
499 if(dim <= 256)
500 bptr = buf = (double *)alloca(dim * sizeof *bptr);
501 else
502 #endif
503 if(!(bptr = buf = (double *)malloc(dim * sizeof *bptr))) {
504 return NULL;
505 }
506 } else {
507 bptr = buf = sbuf;
508 }
509
510 while(dim-- > 0) {
511 *bptr++ = *pos++;
512 }
513
514 res = kd_nearest(tree, buf);
515 #ifndef NO_ALLOCA
516 if(tree->dim > 256)
517 #else
518 if(tree->dim > 16)
519 #endif
520 free(buf);
521 return res;
522 }
523
kd_nearest3(struct kdtree * tree,double x,double y,double z)524 struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z)
525 {
526 double pos[3];
527 pos[0] = x;
528 pos[1] = y;
529 pos[2] = z;
530 return kd_nearest(tree, pos);
531 }
532
kd_nearest3f(struct kdtree * tree,float x,float y,float z)533 struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z)
534 {
535 double pos[3];
536 pos[0] = x;
537 pos[1] = y;
538 pos[2] = z;
539 return kd_nearest(tree, pos);
540 }
541
542 /* ---- nearest N search ---- */
543 /*
544 static kdres *kd_nearest_n(struct kdtree *kd, const double *pos, int num)
545 {
546 int ret;
547 struct kdres *rset;
548
549 if(!(rset = malloc(sizeof *rset))) {
550 return 0;
551 }
552 if(!(rset->rlist = alloc_resnode())) {
553 free(rset);
554 return 0;
555 }
556 rset->rlist->next = 0;
557 rset->tree = kd;
558
559 if((ret = find_nearest_n(kd->root, pos, range, num, rset->rlist, kd->dim)) == -1) {
560 kd_res_free(rset);
561 return 0;
562 }
563 rset->size = ret;
564 kd_res_rewind(rset);
565 return rset;
566 }*/
567
kd_nearest_range(struct kdtree * kd,const double * pos,double range)568 struct kdres *kd_nearest_range(struct kdtree *kd, const double *pos, double range)
569 {
570 int ret;
571 struct kdres *rset;
572
573 if(!(rset = (struct kdres *)malloc(sizeof *rset))) {
574 return NULL;
575 }
576 if(!(rset->rlist = alloc_resnode())) {
577 free(rset);
578 return NULL;
579 }
580 rset->rlist->next = NULL;
581 rset->tree = kd;
582
583 if((ret = find_nearest(kd->root, pos, range, rset->rlist, 0, kd->dim)) == -1) {
584 kd_res_free(rset);
585 return NULL;
586 }
587 rset->size = ret;
588 kd_res_rewind(rset);
589 return rset;
590 }
591
kd_nearest_rangef(struct kdtree * kd,const float * pos,float range)592 struct kdres *kd_nearest_rangef(struct kdtree *kd, const float *pos, float range)
593 {
594 static double sbuf[16];
595 double *bptr, *buf = NULL;
596 int dim = kd->dim;
597 struct kdres *res;
598 assert(dim>0);
599
600 if(dim > 16) {
601 #ifndef NO_ALLOCA
602 if(dim <= 256)
603 bptr = buf = (double *)alloca(dim * sizeof *bptr);
604 else
605 #endif
606 if(!(bptr = buf = (double *)malloc(dim * sizeof *bptr))) {
607 return NULL;
608 }
609 } else {
610 bptr = buf = sbuf;
611 }
612
613 while(dim-- > 0) {
614 *bptr++ = *pos++;
615 }
616
617 res = kd_nearest_range(kd, buf, range);
618 #ifndef NO_ALLOCA
619 if(kd->dim > 256)
620 #else
621 if(kd->dim > 16)
622 #endif
623 free(buf);
624 return res;
625 }
626
kd_nearest_range3(struct kdtree * tree,double x,double y,double z,double range)627 struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double z, double range)
628 {
629 double buf[3];
630 buf[0] = x;
631 buf[1] = y;
632 buf[2] = z;
633 return kd_nearest_range(tree, buf, range);
634 }
635
kd_nearest_range3f(struct kdtree * tree,float x,float y,float z,float range)636 struct kdres *kd_nearest_range3f(struct kdtree *tree, float x, float y, float z, float range)
637 {
638 double buf[3];
639 buf[0] = x;
640 buf[1] = y;
641 buf[2] = z;
642 return kd_nearest_range(tree, buf, range);
643 }
644
kd_res_free(struct kdres * rset)645 void kd_res_free(struct kdres *rset)
646 {
647 clear_results(rset);
648 free_resnode(rset->rlist);
649 free(rset);
650 }
651
kd_res_size(struct kdres * set)652 int kd_res_size(struct kdres *set)
653 {
654 return (set->size);
655 }
656
kd_res_rewind(struct kdres * rset)657 void kd_res_rewind(struct kdres *rset)
658 {
659 rset->riter = rset->rlist->next;
660 }
661
kd_res_end(struct kdres * rset)662 int kd_res_end(struct kdres *rset)
663 {
664 return rset->riter == NULL;
665 }
666
kd_res_next(struct kdres * rset)667 int kd_res_next(struct kdres *rset)
668 {
669 rset->riter = rset->riter->next;
670 return rset->riter != NULL;
671 }
672
kd_res_item(struct kdres * rset,double * pos)673 void *kd_res_item(struct kdres *rset, double *pos)
674 {
675 if(rset->riter) {
676 if(pos) {
677 memcpy(pos, rset->riter->item->pos, rset->tree->dim * sizeof *pos);
678 }
679 return rset->riter->item->data;
680 }
681 return NULL;
682 }
683
kd_res_itemf(struct kdres * rset,float * pos)684 void *kd_res_itemf(struct kdres *rset, float *pos)
685 {
686 if(rset->riter) {
687 if(pos) {
688 int i;
689 for(i=0; i<rset->tree->dim; i++) {
690 pos[i] = rset->riter->item->pos[i];
691 }
692 }
693 return rset->riter->item->data;
694 }
695 return NULL;
696 }
697
kd_res_item3(struct kdres * rset,double * x,double * y,double * z)698 void *kd_res_item3(struct kdres *rset, double *x, double *y, double *z)
699 {
700 if(rset->riter) {
701 if(x) *x = rset->riter->item->pos[0];
702 if(y) *y = rset->riter->item->pos[1];
703 if(z) *z = rset->riter->item->pos[2];
704 return rset->riter->item->data;
705 }
706 return NULL;
707 }
708
kd_res_item3f(struct kdres * rset,float * x,float * y,float * z)709 void *kd_res_item3f(struct kdres *rset, float *x, float *y, float *z)
710 {
711 if(rset->riter) {
712 if(x) *x = rset->riter->item->pos[0];
713 if(y) *y = rset->riter->item->pos[1];
714 if(z) *z = rset->riter->item->pos[2];
715 return rset->riter->item->data;
716 }
717 return NULL;
718 }
719
kd_res_item_data(struct kdres * set)720 void *kd_res_item_data(struct kdres *set)
721 {
722 return kd_res_item(set, NULL);
723 }
724
725 /* ---- hyperrectangle helpers ---- */
hyperrect_create(int dim,const double * min,const double * max)726 static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max)
727 {
728 assert(dim>0);
729 size_t size = dim * sizeof(double);
730 struct kdhyperrect* rect = NULL;
731
732 if (!(rect = (struct kdhyperrect *)malloc(sizeof(struct kdhyperrect)))) {
733 return NULL;
734 }
735
736 rect->dim = dim;
737 if (!(rect->min = (double *)malloc(size))) {
738 free(rect);
739 return NULL;
740 }
741 if (!(rect->max = (double *)malloc(size))) {
742 free(rect->min);
743 free(rect);
744 return NULL;
745 }
746 memcpy(rect->min, min, size);
747 memcpy(rect->max, max, size);
748
749 return rect;
750 }
751
hyperrect_free(struct kdhyperrect * rect)752 static void hyperrect_free(struct kdhyperrect *rect)
753 {
754 free(rect->min);
755 free(rect->max);
756 free(rect);
757 }
758
hyperrect_duplicate(const struct kdhyperrect * rect)759 static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect)
760 {
761 return hyperrect_create(rect->dim, rect->min, rect->max);
762 }
763
hyperrect_extend(struct kdhyperrect * rect,const double * pos)764 static void hyperrect_extend(struct kdhyperrect *rect, const double *pos)
765 {
766 int i;
767
768 for (i=0; i < rect->dim; i++) {
769 if (pos[i] < rect->min[i]) {
770 rect->min[i] = pos[i];
771 }
772 if (pos[i] > rect->max[i]) {
773 rect->max[i] = pos[i];
774 }
775 }
776 }
777
hyperrect_dist_sq(struct kdhyperrect * rect,const double * pos)778 static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos)
779 {
780 int i;
781 double result = 0;
782
783 for (i=0; i < rect->dim; i++) {
784 if (pos[i] < rect->min[i]) {
785 result += SQ(rect->min[i] - pos[i]);
786 } else if (pos[i] > rect->max[i]) {
787 result += SQ(rect->max[i] - pos[i]);
788 }
789 }
790
791 return result;
792 }
793
794 /* ---- static helpers ---- */
795
796 #ifdef USE_LIST_NODE_ALLOCATOR
797 /* special list node allocators. */
798 static struct res_node *free_nodes;
799
800 #ifndef NO_PTHREADS
801 static pthread_mutex_t alloc_mutex = PTHREAD_MUTEX_INITIALIZER;
802 #endif
803
alloc_resnode(void)804 static struct res_node *alloc_resnode(void)
805 {
806 struct res_node *node;
807
808 #ifndef NO_PTHREADS
809 pthread_mutex_lock(&alloc_mutex);
810 #endif
811
812 if(!free_nodes) {
813 node = (struct res_node *)malloc(sizeof *node);
814 } else {
815 node = free_nodes;
816 free_nodes = free_nodes->next;
817 node->next = NULL;
818 }
819
820 #ifndef NO_PTHREADS
821 pthread_mutex_unlock(&alloc_mutex);
822 #endif
823
824 return node;
825 }
826
free_resnode(struct res_node * node)827 static void free_resnode(struct res_node *node)
828 {
829 #ifndef NO_PTHREADS
830 pthread_mutex_lock(&alloc_mutex);
831 #endif
832
833 node->next = free_nodes;
834 free_nodes = node;
835
836 #ifndef NO_PTHREADS
837 pthread_mutex_unlock(&alloc_mutex);
838 #endif
839 }
840 #endif /* list node allocator or not */
841
842
843 /* inserts the item. if dist_sq is >= 0, then do an ordered insert */
844 /* TODO make the ordering code use heapsort */
rlist_insert(struct res_node * list,struct kdnode * item,double dist_sq)845 static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq)
846 {
847 struct res_node *rnode;
848
849 if(!(rnode = alloc_resnode())) {
850 return -1;
851 }
852 rnode->item = item;
853 rnode->dist_sq = dist_sq;
854
855 if(dist_sq >= 0.0) {
856 while(list->next && list->next->dist_sq < dist_sq) {
857 list = list->next;
858 }
859 }
860 rnode->next = list->next;
861 list->next = rnode;
862 return 0;
863 }
864
clear_results(struct kdres * rset)865 static void clear_results(struct kdres *rset)
866 {
867 struct res_node *tmp, *node = rset->rlist->next;
868
869 while(node) {
870 tmp = node;
871 node = node->next;
872 free_resnode(tmp);
873 }
874
875 rset->rlist->next = NULL;
876 }
877