1 /***************************************************************************/
2 /* */
3 /* ROUTINES FOR FINDING NEAREST NEIGHBORS */
4 /* */
5 /* (Based on Jon Bentley's paper "K-d trees for semidynamic point sets") */
6 /* */
7 /* TSP CODE */
8 /* */
9 /* */
10 /* Written by: Applegate, Bixby, Chvatal, and Cook */
11 /* Date: February 24, 1995 (cofeb24) */
12 /* Changes: August 6 (bico) - added wcoord to fixed radius search */
13 /* */
14 /* */
15 /* EXPORTED FUNCTIONS: */
16 /* int CCkdtree_k_nearest (CCkdtree *kt, int ncount, int k, */
17 /* CCdatagroup *dat, double *wcoord, int wantlist, int *ocount, */
18 /* int **olist) */
19 /* RETURNS the k-nearest neighbor graph. */
20 /* -kt can be NULL, otherwise it should point to a CCkdtree buildt */
21 /* by a call to kdbuild () */
22 /* -ncount is the number of points. */
23 /* -k is the number of nearest neighbors wanted. */
24 /* -wcoord is an array of node weights (like Held-Karp), it can */
25 /* be NULL. The weights should be nonnegative. */
26 /* -wantlist is 1 if you want the function to return the edges. */
27 /* -ocount returns the number of edges (if wantlist is 1) and */
28 /* olist returns the edgelist is end1 end2 format. */
29 /* int CCkdtree_quadrant_k_nearest (CCkdtree *kt, int ncount, int k, */
30 /* CCdatagroup *dat, double *wcoord, */
31 /* int wantlist, int *ocount, int **olist) */
32 /* RETURNS the quadrant k-nearest neighbor graph. */
33 /* -see CCkdtree_k_nearest. */
34 /* int CCkdtree_node_k_nearest (CCkdtree *kt, int ncount, int n, int k, */
35 /* CCdatagroup *dat, double *wcoord, int *list) */
36 /* RETURNS the k nearest points to point n. */
37 /* -The k points are return in list (and list must be allocated by */
38 /* calling routine. */
39 /* -kt is a pointer to a CCkdtree previously built by */
40 /* CCkdtree_build. */
41 /* int CCkdtree_node_quadrant_k_nearest (CCkdtree *kt, int ncount, */
42 /* int n, int k, CCdatagroup *dat, double *wcoord, int *list) */
43 /* RETURNS the quadrant k nearest point to point n. */
44 /* -see CCkdtree_node_k_nearest. */
45 /* int CCkdtree_node_nearest (ktree *kt, int n, CCdatagroup *dat, */
46 /* double *wcoord) */
47 /* RETURNS the nearest point to point n. */
48 /* -kt CANNOT be NULL. */
49 /* -The point is returned as the function value. kt is a pointer */
50 /* to a CCkdtree (previously buildt by a call to CCkdtree_build) */
51 /* int CCkdtree_fixed_radius_nearest (CCkdtree *kt, CCdatagroup *dat, */
52 /* double *wcoord, int n, double rad, */
53 /* int (*doit_fn) (int, int, void *), void *pass_param) */
54 /* ACTION: Calls the function doit_fn (n, a, void *), where a ranges */
55 /* over all points within distance rad of the point n. The */
56 /* void * field can be used to bundle a group of parmeters */
57 /* into pass_param that will be passed to doit_fn. */
58 /* -kt CANNOT be NULL. */
59 /* -doit_fn can also call CCkdtree_fixed_radius_nearest (no globals */
60 /* are set by the function calls) */
61 /* -pass_param can be NULL or used to point to a structure with */
62 /* with parameters for doit_fn. */
63 /* int CCkdtree_nearest_neighbor_tour (CCkdtree *kt, int ncount, */
64 /* int start, CCdatagroup *dat, int *outcycle, double *val) */
65 /* -kt can be NULL. */
66 /* -Node weights are not used. */
67 /* -start is the starting node for the tour. */
68 /* -if outcycle is not NULL, then it should point to a array of */
69 /* length at least ncount (allocated by the calling routine). The */
70 /* cycle will be returned in the array in node node node format. */
71 /* -the length of the tour is return in val. */
72 /* int CCkdtree_nearest_neighbor_2match (CCkdtree *kt, int ncount, */
73 /* int start, CCdatagroup *dat, int *outmatch, double *val) */
74 /* -Like CCkdtree_nearest_neighbor_tour. If outmatch is not NULL */
75 /* then it should point to an array of length at leaast 2*ncount. */
76 /* */
77 /* NOTES: */
78 /* If memory is tight, use CCkdtree_node_k_nearest to get the edges */
79 /* one node at a time. (CCkdtree_k_nearest () builds a hash table to */
80 /* avoid duplicate edges, and it will use 8 * nedges bytes.) */
81 /* CCkdtree_node_nearest returns the nearest point as the function */
82 /* value; CCkdtree_fixed_radius_nearest returns 1 if doit_fn returns a */
83 /* a nonzero value, otherwise it returns 0; all other routines return 0 */
84 /* if successful and 1 otherwise. */
85 /* */
86 /***************************************************************************/
87
88
89 #include "machdefs.h"
90 #include "util.h"
91 #include "kdtree.h"
92
93 #define BIGDOUBLE 100000000.0
94 #define SWAP(x,y,t) ((t)=(x),(x)=(y),(y)=(t))
95 #define NEAR_HEAP_CUTOFF 100 /* When to switch from list to heap */
96 /* On an RS6000, the heap started winning */
97 /* at 100 (by 200 it was 3 times faster) */
98
99 typedef struct shortedge {
100 int end;
101 double length;
102 } shortedge;
103
104 typedef struct intptr {
105 int this;
106 struct intptr *next;
107 } intptr;
108
109 #define Fedgelen(n1, n2) \
110 (datw != (double *) NULL ? \
111 CCutil_dat_edgelen ((n1), (n2), dat) \
112 + datw[(n1)] + datw[(n2)] : \
113 CCutil_dat_edgelen ((n1), (n2), dat))
114
115 #ifdef CC_PROTOTYPE_ANSI
116
117 static void
118 node_k_nearest_work (CCkdtree *thetree, CCdatagroup *dat, double *datw,
119 CCkdnode *p, CCdheap *near_heap, int *heap_names, int *heap_count,
120 int target, int num, shortedge *nearlist, double *worst_on_list,
121 CCkdbnds *box),
122 node_nearest_work (CCkdtree *thetree, CCdatagroup *dat, double *datw,
123 CCkdnode *p, int target, double *ndist, int *nnode);
124 static int
125 run_kdtree_k_nearest (CCkdtree *kt, int ncount, int k, CCdatagroup *dat,
126 double *wcoord, int wantlist, int *ocount, int **olist, int doquad),
127 put_in_table (int i, int j, int *added, intptr **table),
128 q_run_it (CCkdtree *thetree, CCdatagroup *dat, double *datw, int *llist,
129 int *lcount, int *list, int target, int num, CCkdbnds *box),
130 run_kdtree_node_k_nearest (CCkdtree *thetree, CCdatagroup *dat,
131 double *datw,
132 int *list, int target, int num, CCkdbnds *box),
133 ball_in_bounds (CCdatagroup *dat, CCkdbnds *bnds, int n, double dist),
134 fixed_radius_nearest_work (CCkdtree *thetree, CCkdnode *p,
135 int (*doit_fn) (int, int, void *),
136 int target, double dist, CCdatagroup *dat, double *wcoord,
137 double xtarget, double ytarget, void *pass_param);
138
139 #else
140
141 static void
142 node_k_nearest_work (),
143 node_nearest_work ();
144 static int
145 run_kdtree_k_nearest (),
146 put_in_table (),
147 q_run_it (),
148 run_kdtree_node_k_nearest (),
149 ball_in_bounds (),
150 fixed_radius_nearest_work ();
151
152 #endif
153
CC_PTR_ALLOC_ROUTINE(intptr,intptralloc,intptrchunklist,intptrfreelist)154 CC_PTR_ALLOC_ROUTINE (intptr, intptralloc, intptrchunklist, intptrfreelist)
155 CC_PTR_FREE_ROUTINE (intptr, intptrfree, intptrfreelist)
156 CC_PTR_FREE_WORLD_ROUTINE (intptr, intptrfree_world, intptrchunklist,
157 intptrfreelist)
158 CC_PTR_LEAKS_ROUTINE (intptr, intptr_check_leaks, intptrchunklist,
159 intptrfreelist, this, int)
160
161 #ifdef CC_PROTOTYPE_ANSI
162 int CCkdtree_k_nearest (CCkdtree *kt, int ncount, int k, CCdatagroup *dat,
163 double *wcoord, int wantlist, int *ocount, int **olist)
164 #else
165 int CCkdtree_k_nearest (kt, ncount, k, dat, wcoord, wantlist, ocount, olist)
166 CCkdtree *kt;
167 int ncount, k;
168 CCdatagroup *dat;
169 double *wcoord;
170 int wantlist, *ocount, **olist;
171 #endif
172 {
173 return run_kdtree_k_nearest (kt, ncount, k, dat, wcoord,
174 wantlist, ocount, olist, 0);
175 }
176
177 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_quadrant_k_nearest(CCkdtree * kt,int ncount,int k,CCdatagroup * dat,double * wcoord,int wantlist,int * ocount,int ** olist)178 int CCkdtree_quadrant_k_nearest (CCkdtree *kt, int ncount, int k,
179 CCdatagroup *dat, double *wcoord, int wantlist, int *ocount,
180 int **olist)
181 #else
182 int CCkdtree_quadrant_k_nearest (kt, ncount, k, dat, wcoord, wantlist, ocount,
183 olist)
184 CCkdtree *kt;
185 int ncount, k;
186 CCdatagroup *dat;
187 double *wcoord;
188 int wantlist, *ocount, **olist;
189 #endif
190 {
191 return run_kdtree_k_nearest (kt, ncount, k, dat, wcoord,
192 wantlist, ocount, olist, 1);
193 }
194
195
196 #ifdef CC_PROTOTYPE_ANSI
run_kdtree_k_nearest(CCkdtree * kt,int ncount,int k,CCdatagroup * dat,double * wcoord,int wantlist,int * ocount,int ** olist,int doquad)197 static int run_kdtree_k_nearest (CCkdtree *kt, int ncount, int k,
198 CCdatagroup *dat, double *wcoord, int wantlist, int *ocount,
199 int **olist, int doquad)
200 #else
201 static int run_kdtree_k_nearest (kt, ncount, k, dat, wcoord, wantlist, ocount,
202 olist, doquad)
203 CCkdtree *kt;
204 int ncount, k;
205 CCdatagroup *dat;
206 double *wcoord;
207 int wantlist, *ocount, **olist;
208 int doquad;
209 #endif
210 {
211 int i, n;
212 intptr *ip, *ipnext;
213 int total, onlist;
214 CCkdtree localkt, *mykt;
215 int added, ntotal = 0;
216 int rval = 0;
217 int *list = (int *) NULL;
218 int goal = (doquad ? (4 * k) : k);
219 int newtree = 0;
220 intptr **table = (intptr **) NULL;
221
222 if (wcoord != (double *) NULL) {
223 for (i = 0; i < ncount; i++) {
224 if (wcoord[i] < -0.00000001) {
225 fprintf (stderr, "Cannot CCkdtree with negative node weights\n");
226 return 1;
227 }
228 }
229 }
230
231 if (wantlist) {
232 *ocount = 0;
233 *olist = (int *) NULL;
234 }
235
236 if (kt == (CCkdtree *) NULL) {
237 if (CCkdtree_build (&localkt, ncount, dat, wcoord)) {
238 fprintf (stderr, "Unable to build CCkdtree\n");
239 return 1;
240 }
241 mykt = &localkt;
242 newtree = 1;
243 } else {
244 mykt = kt;
245 }
246
247
248 table = CC_SAFE_MALLOC (ncount, intptr *);
249 if (!table) {
250 rval = 1;
251 goto CLEANUP;
252 }
253 for (i = 0; i < ncount; i++)
254 table[i] = (intptr *) NULL;
255 list = CC_SAFE_MALLOC (goal, int);
256 if (!list) {
257 rval = 1;
258 goto CLEANUP;
259 }
260
261 for (n = 0; n < ncount; n++) {
262 if (doquad) {
263 if (CCkdtree_node_quadrant_k_nearest (mykt, ncount, n, k, dat,
264 wcoord, list)) {
265 rval = 1;
266 goto CLEANUP;
267 }
268 } else {
269 if (CCkdtree_node_k_nearest (mykt, ncount, n, k, dat, wcoord,
270 list)) {
271 rval = 1;
272 goto CLEANUP;
273 }
274 }
275 for (i = 0; i < goal; i++) {
276 if (list[i] != -1) {
277 if (put_in_table (n, list[i], &added, table)) {
278 rval = 1;
279 goto CLEANUP;
280 } else {
281 ntotal += added;
282 }
283 }
284 }
285 /*
286 if (n == 0) {
287 printf ("Neighbors of Node %d (%d, %d) :\n", n,
288 (int) dat->x[n], (int) dat->y[n]);
289 for (i = 0; i < goal; i++) {
290 if (list[i] != -1) {
291 printf ("%d %d (%d, %d)\n", list[i],
292 CCutil_dat_edgelen (n, list[i], dat),
293 (int) dat->x[list[i]], (int) dat->y[list[i]]);
294 }
295 }
296 }
297 */
298 if (n % 1000 == 999) {
299 printf (".");
300 fflush (stdout);
301 }
302 }
303 printf (" %d edges\n", ntotal);
304 fflush (stdout);
305
306 if (wantlist) {
307 int j = 0;
308 *olist = CC_SAFE_MALLOC (2 * ntotal, int);
309 if (!(*olist)) {
310 rval = 1;
311 goto CLEANUP;
312 }
313 *ocount = ntotal;
314 for (i = 0; i < ncount; i++) {
315 for (ip = table[i]; ip; ip = ipnext) {
316 ipnext = ip->next;
317 (*olist)[j++] = i;
318 (*olist)[j++] = ip->this;
319 intptrfree (ip);
320 }
321 table[i] = (intptr *) NULL;
322 }
323 } else {
324 for (i = 0; i < ncount; i++) {
325 for (ip = table[i]; ip; ip = ipnext) {
326 ipnext = ip->next;
327 intptrfree (ip);
328 }
329 table[i] = (intptr *) NULL;
330 }
331 }
332 if (intptr_check_leaks (&total, &onlist)) {
333 fprintf (stderr, "WARNING: %d outstanding intptrs in kdnear\n",
334 total - onlist);
335 }
336
337 CLEANUP:
338
339 intptrfree_world ();
340 if (table)
341 CC_FREE (table, intptr *);
342 if (list)
343 CC_FREE (list, int);
344 if (newtree)
345 CCkdtree_free (&localkt);
346
347 return rval;
348 }
349
350 #ifdef CC_PROTOTYPE_ANSI
put_in_table(int i,int j,int * added,intptr ** table)351 static int put_in_table (int i, int j, int *added, intptr **table)
352 #else
353 static int put_in_table (i, j, added, table)
354 int i, j;
355 int *added;
356 intptr **table;
357 #endif
358 {
359 intptr *ip;
360
361 if (j < i) {
362 int temp;
363 SWAP(i, j, temp);
364 }
365
366 for (ip = table[i]; ip; ip = ip->next)
367 if (ip->this == j) {
368 *added = 0;
369 return 0;
370 }
371 ip = intptralloc ();
372 if (!ip) {
373 *added = 0;
374 return 1;
375 }
376 ip->this = j;
377 ip->next = table[i];
378 table[i] = ip;
379 *added = 1;
380 return 0;
381 }
382
383 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_node_quadrant_k_nearest(CCkdtree * kt,int ncount,int n,int k,CCdatagroup * dat,double * wcoord,int * list)384 int CCkdtree_node_quadrant_k_nearest (CCkdtree *kt, int ncount, int n, int k,
385 CCdatagroup *dat, double *wcoord, int *list)
386 #else
387 int CCkdtree_node_quadrant_k_nearest (kt, ncount, n, k, dat, wcoord, list)
388 CCkdtree *kt;
389 int ncount, n, k;
390 CCdatagroup *dat;
391 double *wcoord;
392 int *list;
393 #endif
394 {
395 CCkdbnds localbnds;
396 int i, lcount = 0;
397 int *llist = (int *) NULL;
398 int rval = 0;
399 CCkdtree localkt;
400 CCkdtree *thetree;
401 int newtree = 0;
402
403 if (kt == (CCkdtree *) NULL) {
404 if (CCkdtree_build (&localkt, ncount, dat, wcoord)) {
405 fprintf (stderr, "Unable to build CCkdtree\n");
406 return 1;
407 }
408 thetree = &localkt;
409 newtree = 1;
410 } else {
411 thetree = kt;
412 }
413
414 llist = CC_SAFE_MALLOC (k, int);
415 if (!llist) {
416 rval = 1;
417 goto CLEANUP;
418 }
419
420 localbnds.x[0] = dat->x[n];
421 localbnds.x[1] = BIGDOUBLE;
422 localbnds.y[0] = dat->y[n];
423 localbnds.y[1] = BIGDOUBLE;
424 if (q_run_it (thetree, dat, wcoord, llist, &lcount, list, n, k,
425 &localbnds)) {
426 fprintf (stderr, "run_kdtree_node_k_nearest failed\n");
427 rval = 1;
428 goto CLEANUP;
429 }
430
431 localbnds.x[0] = dat->x[n];
432 localbnds.x[1] = BIGDOUBLE;
433 localbnds.y[0] = -BIGDOUBLE;
434 localbnds.y[1] = dat->y[n];
435 if (q_run_it (thetree, dat, wcoord, llist, &lcount, list, n, k,
436 &localbnds)) {
437 fprintf (stderr, "run_kdtree_node_k_nearest failed\n");
438 rval = 1;
439 goto CLEANUP;
440 }
441
442 localbnds.x[0] = -BIGDOUBLE;
443 localbnds.x[1] = dat->x[n];
444 localbnds.y[0] = -BIGDOUBLE;
445 localbnds.y[1] = dat->y[n];
446 if (q_run_it (thetree, dat, wcoord, llist, &lcount, list, n, k,
447 &localbnds)) {
448 fprintf (stderr, "run_kdtree_node_k_nearest failed\n");
449 rval = 1;
450 goto CLEANUP;
451 }
452
453 localbnds.x[0] = -BIGDOUBLE;
454 localbnds.x[1] = dat->x[n];
455 localbnds.y[0] = dat->y[n];
456 localbnds.y[1] = BIGDOUBLE;
457 if (q_run_it (thetree, dat, wcoord, llist, &lcount, list, n, k,
458 &localbnds)) {
459 fprintf (stderr, "run_kdtree_node_k_nearest failed\n");
460 rval = 1;
461 goto CLEANUP;
462 }
463
464 for (i = lcount; i < (4 * k); i++)
465 list[i] = -1;
466
467 CLEANUP:
468
469 CC_FREE (llist, int);
470 if (newtree)
471 CCkdtree_free (&localkt);
472
473 return rval;
474 }
475
476 #ifdef CC_PROTOTYPE_ANSI
q_run_it(CCkdtree * thetree,CCdatagroup * dat,double * datw,int * llist,int * lcount,int * list,int target,int num,CCkdbnds * box)477 static int q_run_it (CCkdtree *thetree, CCdatagroup *dat, double *datw,
478 int *llist, int *lcount, int *list, int target, int num, CCkdbnds *box)
479 #else
480 static int q_run_it (thetree, dat, datw, llist, lcount, list, target, num, box)
481 CCkdtree *thetree;
482 CCdatagroup *dat;
483 double *datw;
484 int *llist;
485 int *lcount;
486 int *list;
487 int target;
488 int num;
489 CCkdbnds *box;
490 #endif
491 {
492 int i, j;
493
494 if (run_kdtree_node_k_nearest (thetree, dat, datw, llist, target, num,
495 box))
496 return 1;
497 for (i = 0; i < num; i++) {
498 if (llist[i] != -1) {
499 for (j = 0; j < *lcount; j++)
500 if (list[j] == llist[i])
501 break;
502 if (j == *lcount)
503 list[(*lcount)++] = llist[i];
504 }
505 }
506 return 0;
507 }
508
509 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_node_k_nearest(CCkdtree * kt,int ncount,int n,int k,CCdatagroup * dat,double * wcoord,int * list)510 int CCkdtree_node_k_nearest (CCkdtree *kt, int ncount, int n, int k,
511 CCdatagroup *dat, double *wcoord, int *list)
512 #else
513 int CCkdtree_node_k_nearest (kt, ncount, n, k, dat, wcoord, list)
514 CCkdtree *kt;
515 int ncount, n, k;
516 CCdatagroup *dat;
517 double *wcoord;
518 int *list;
519 #endif
520 {
521 CCkdtree localkt;
522 CCkdtree *thetree;
523 int newtree = 0;
524 int rval = 0;
525
526 if (kt == (CCkdtree *) NULL) {
527 if (CCkdtree_build (&localkt, ncount, dat, wcoord)) {
528 fprintf (stderr, "Unable to build CCkdtree\n");
529 return 1;
530 }
531 thetree = &localkt;
532 newtree = 1;
533 } else {
534 thetree = kt;
535 }
536
537 rval = run_kdtree_node_k_nearest (thetree, dat, wcoord, list, n, k,
538 (CCkdbnds *) NULL);
539 if (newtree)
540 CCkdtree_free (&localkt);
541 return rval;
542 }
543
544 #ifdef CC_PROTOTYPE_ANSI
run_kdtree_node_k_nearest(CCkdtree * thetree,CCdatagroup * dat,double * datw,int * list,int target,int num,CCkdbnds * box)545 static int run_kdtree_node_k_nearest (CCkdtree *thetree, CCdatagroup *dat,
546 double *datw, int *list, int target, int num, CCkdbnds *box)
547 #else
548 static int run_kdtree_node_k_nearest (thetree, dat, datw, list, target, num,
549 box)
550 CCkdtree *thetree;
551 CCdatagroup *dat;
552 double *datw;
553 int *list;
554 int target;
555 int num;
556 CCkdbnds *box;
557 #endif
558 {
559 int i;
560 CCkdnode *p, *lastp;
561 double diff;
562 CCdheap near_heap;
563 int *heap_names = (int *) NULL;
564 int heap_count = 0;
565 shortedge *nearlist = (shortedge *) NULL;
566 double worst_on_list = BIGDOUBLE;
567
568 if (num >= NEAR_HEAP_CUTOFF) {
569 if (CCutil_dheap_init (&near_heap, num))
570 return 1;
571 heap_names = CC_SAFE_MALLOC (num, int);
572 if (!heap_names) {
573 CCutil_dheap_free (&near_heap);
574 return 1;
575 }
576 heap_count = 0;
577 } else {
578 nearlist = CC_SAFE_MALLOC (num + 1, shortedge);
579 if (!nearlist) {
580 CCutil_dheap_free (&near_heap);
581 CC_FREE (heap_names, int);
582 return 1;
583 }
584 for (i = 0; i < num; i++)
585 nearlist[i].length = BIGDOUBLE;
586 nearlist[num].length = -BIGDOUBLE;
587 }
588
589 /*
590 To do top down search just use:
591
592 node_k_nearest_work (thetree->root);
593 */
594
595 p = thetree->bucketptr[target];
596 node_k_nearest_work (thetree, dat, datw, p, &near_heap, heap_names,
597 &heap_count, target, num, nearlist, &worst_on_list,
598 box);
599 while (1) {
600 lastp = p;
601 p = p->father;
602 if (p == (CCkdnode *) NULL)
603 break;
604 switch (p->cutdim) {
605 case 0:
606 diff = p->cutval - dat->x[target];
607 if (lastp == p->loson) { /* So target is on low side */
608 if (worst_on_list > (double) ((int) diff))
609 if (box == (CCkdbnds *) NULL || p->cutval <= box->x[1])
610 node_k_nearest_work (thetree, dat, datw, p->hison,
611 &near_heap, heap_names, &heap_count, target,
612 num, nearlist, &worst_on_list, box);
613 } else {
614 if (worst_on_list > (double) ((int) -diff))
615 if (box == (CCkdbnds *) NULL || p->cutval >= box->x[0])
616 node_k_nearest_work (thetree, dat, datw, p->loson,
617 &near_heap, heap_names, &heap_count, target,
618 num, nearlist, &worst_on_list, box);
619 }
620 break;
621 case 1:
622 diff = p->cutval - dat->y[target];
623 if (lastp == p->loson) {
624 if (worst_on_list > (double) ((int) diff))
625 if (box == (CCkdbnds *) NULL || p->cutval <= box->y[1])
626 node_k_nearest_work (thetree, dat, datw, p->hison,
627 &near_heap, heap_names, &heap_count, target,
628 num, nearlist, &worst_on_list, box);
629 } else {
630 if (worst_on_list > (double) ((int) -diff))
631 if (box == (CCkdbnds *) NULL || p->cutval >= box->y[0])
632 node_k_nearest_work (thetree, dat, datw, p->loson,
633 &near_heap, heap_names, &heap_count, target,
634 num, nearlist, &worst_on_list, box);
635 }
636 break;
637 case 2:
638 if (lastp == p->loson) {
639 if (worst_on_list > p->cutval + datw[target])
640 node_k_nearest_work (thetree, dat, datw, p->hison,
641 &near_heap, heap_names, &heap_count, target,
642 num, nearlist, &worst_on_list, box);
643 } else {
644 node_k_nearest_work (thetree, dat, datw, p->loson, &near_heap,
645 heap_names, &heap_count, target, num, nearlist,
646 &worst_on_list, box);
647 }
648 break;
649 }
650 if (datw == (double *) NULL && p->bnds &&
651 ball_in_bounds (dat, p->bnds, target, worst_on_list))
652 /* Doing extra check for box with quad-nearest appears to slow */
653 /* things down. */
654 break;
655 }
656
657 if (num >= NEAR_HEAP_CUTOFF) {
658 if (heap_count < num) {
659 if (box == (CCkdbnds *) NULL) {
660 fprintf (stderr, "WARNING: There do not exist %d neighbors\n",
661 num);
662 }
663 for (i = 0; i < heap_count; i++) {
664 list[i] = heap_names[i];
665 }
666 for (; i < num; i++)
667 list[i] = -1;
668 } else {
669 for (i = 0; i < num; i++)
670 list[i] = heap_names[i];
671 }
672 } else {
673 int ntot = 0;
674 for (i = 0; i < num; i++) {
675 if (nearlist[i].length < BIGDOUBLE)
676 list[ntot++] = nearlist[i].end;
677 }
678 if (ntot < num) {
679 if (box == (CCkdbnds *) NULL) {
680 fprintf (stderr, "WARNING: There do not exist %d neighbors\n",
681 num);
682 }
683 for (i = ntot; i < num; i++)
684 list[i] = -1;
685 }
686 }
687
688 if (num >= NEAR_HEAP_CUTOFF) {
689 CC_FREE (heap_names, int);
690 CCutil_dheap_free (&near_heap);
691 } else {
692 CC_FREE (nearlist, shortedge);
693 }
694 return 0;
695 }
696
697 #ifdef CC_PROTOTYPE_ANSI
node_k_nearest_work(CCkdtree * thetree,CCdatagroup * dat,double * datw,CCkdnode * p,CCdheap * near_heap,int * heap_names,int * heap_count,int target,int num,shortedge * nearlist,double * worst_on_list,CCkdbnds * box)698 static void node_k_nearest_work (CCkdtree *thetree, CCdatagroup *dat,
699 double *datw, CCkdnode *p, CCdheap *near_heap, int *heap_names,
700 int *heap_count, int target, int num, shortedge *nearlist,
701 double *worst_on_list, CCkdbnds *box)
702 #else
703 static void node_k_nearest_work (thetree, dat, datw, p, near_heap, heap_names,
704 heap_count, target, num, nearlist, worst_on_list, box)
705 CCkdtree *thetree;
706 CCdatagroup *dat;
707 double *datw;
708 CCkdnode *p;
709 CCdheap *near_heap;
710 int *heap_names;
711 int *heap_count;
712 int target, num;
713 shortedge *nearlist;
714 double *worst_on_list;
715 CCkdbnds *box;
716 #endif
717 {
718 int i, h, k;
719 double val, thisx, thisdist;
720
721 if (p->bucket) {
722 if (num >= NEAR_HEAP_CUTOFF) {
723 for (i = p->lopt; i <= p->hipt; i++) {
724 if (thetree->perm[i] != target) {
725 if (box == (CCkdbnds *) NULL ||
726 (dat->x[thetree->perm[i]] >= box->x[0] &&
727 dat->x[thetree->perm[i]] <= box->x[1] &&
728 dat->y[thetree->perm[i]] >= box->y[0] &&
729 dat->y[thetree->perm[i]] <= box->y[1])) {
730 thisdist = Fedgelen (thetree->perm[i], target);
731 if (*heap_count < num) {
732 near_heap->key[*heap_count] = -thisdist;
733 heap_names[*heap_count] = thetree->perm[i];
734 CCutil_dheap_insert (near_heap, *heap_count);
735 (*heap_count)++;
736 } else if (*worst_on_list > thisdist) {
737 h = CCutil_dheap_deletemin (near_heap);
738 heap_names[h] = thetree->perm[i];
739 near_heap->key[h] = -thisdist;
740 CCutil_dheap_insert (near_heap, h);
741 h = CCutil_dheap_findmin (near_heap);
742 *worst_on_list = -near_heap->key[h];
743 }
744 }
745 }
746 }
747 } else {
748 for (i = p->lopt; i <= p->hipt; i++) {
749 if (thetree->perm[i] != target) {
750 if (box == (CCkdbnds *) NULL ||
751 (dat->x[thetree->perm[i]] >= box->x[0] &&
752 dat->x[thetree->perm[i]] <= box->x[1] &&
753 dat->y[thetree->perm[i]] >= box->y[0] &&
754 dat->y[thetree->perm[i]] <= box->y[1])) {
755 thisdist = Fedgelen (thetree->perm[i], target);
756 if (*worst_on_list > thisdist) {
757 for (k = 0; nearlist[k+1].length > thisdist; k++) {
758 nearlist[k].end = nearlist[k + 1].end;
759 nearlist[k].length = nearlist[k + 1].length;
760 }
761 nearlist[k].length = thisdist;
762 nearlist[k].end = thetree->perm[i];
763 *worst_on_list = nearlist[0].length;
764 }
765 }
766 }
767 }
768 }
769 } else {
770 val = p->cutval;
771 switch (p->cutdim) {
772 case 0:
773 thisx = dat->x[target];
774 if (thisx < val) {
775 node_k_nearest_work (thetree, dat, datw, p->loson, near_heap,
776 heap_names, heap_count, target, num, nearlist,
777 worst_on_list, box);
778 /* Truncation for floating point coords */
779 if (*worst_on_list > (double) ((int) (val - thisx)))
780 if (box == (CCkdbnds *) NULL || val >= box->x[0])
781 node_k_nearest_work (thetree, dat, datw, p->hison,
782 near_heap, heap_names, heap_count, target,
783 num, nearlist, worst_on_list, box);
784 } else {
785 node_k_nearest_work (thetree, dat, datw, p->hison, near_heap,
786 heap_names, heap_count, target, num, nearlist,
787 worst_on_list, box);
788 if (*worst_on_list > (double) ((int) (thisx - val)))
789 if (box == (CCkdbnds *) NULL || val <= box->x[1])
790 node_k_nearest_work (thetree, dat, datw, p->loson,
791 near_heap, heap_names, heap_count, target,
792 num, nearlist, worst_on_list, box);
793 }
794 break;
795 case 1:
796 thisx = dat->y[target];
797 if (thisx < val) {
798 node_k_nearest_work (thetree, dat, datw, p->loson, near_heap,
799 heap_names, heap_count, target, num, nearlist,
800 worst_on_list, box);
801 if (*worst_on_list > (double) ((int) (val - thisx)))
802 if (box == (CCkdbnds *) NULL || val >= box->y[0])
803 node_k_nearest_work (thetree, dat, datw, p->hison,
804 near_heap, heap_names, heap_count, target,
805 num, nearlist, worst_on_list, box);
806 } else {
807 node_k_nearest_work (thetree, dat, datw, p->hison, near_heap,
808 heap_names, heap_count, target, num, nearlist,
809 worst_on_list, box);
810 if (*worst_on_list > (double) ((int) (thisx - val)))
811 if (box == (CCkdbnds *) NULL || val <= box->y[1])
812 node_k_nearest_work (thetree, dat, datw, p->loson,
813 near_heap, heap_names, heap_count, target,
814 num, nearlist, worst_on_list, box);
815 }
816 break;
817 case 2:
818 thisx = datw[target];
819 node_k_nearest_work (thetree, dat, datw, p->loson, near_heap,
820 heap_names, heap_count, target, num, nearlist,
821 worst_on_list, box);
822 if (*worst_on_list > val + thisx)
823 node_k_nearest_work (thetree, dat, datw, p->hison, near_heap,
824 heap_names, heap_count, target, num, nearlist,
825 worst_on_list, box);
826 break;
827 }
828 }
829 }
830
831
832 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_node_nearest(CCkdtree * kt,int n,CCdatagroup * dat,double * wcoord)833 int CCkdtree_node_nearest (CCkdtree *kt, int n, CCdatagroup *dat,
834 double *wcoord)
835 #else
836 int CCkdtree_node_nearest (kt, n, dat, wcoord)
837 CCkdtree *kt;
838 int n;
839 CCdatagroup *dat;
840 double *wcoord;
841 #endif
842 {
843 CCkdnode *p, *lastp;
844 double diff;
845 double ndist = BIGDOUBLE;
846 int nnode;
847 CCkdtree *thetree = (CCkdtree *) NULL;
848
849 if (kt == (CCkdtree *) NULL) {
850 fprintf (stderr, "ERROR: kt cannot be NULL in CCkdtree_node_nearest)\n");
851 return n;
852 }
853
854 thetree = kt;
855
856 ndist = BIGDOUBLE;
857 nnode = n;
858
859 /*
860 To do top down search just use:
861
862 node_nearest_work (kt->root);
863 return nnode;
864 */
865
866 p = kt->bucketptr[n];
867 node_nearest_work (thetree, dat, wcoord, p, n, &ndist, &nnode);
868 while (1) {
869 lastp = p;
870 p = p->father;
871 if (p == (CCkdnode *) NULL)
872 break;
873 switch (p->cutdim) {
874 case 0:
875 diff = p->cutval - dat->x[n];
876 if (lastp == p->loson) {
877 if (ndist > (double) ((int) diff))
878 node_nearest_work (thetree, dat, wcoord, p->hison, n,
879 &ndist, &nnode);
880 } else {
881 if (ndist > (double) ((int) (-diff)))
882 node_nearest_work (thetree, dat, wcoord, p->loson, n,
883 &ndist, &nnode);
884 }
885 break;
886 case 1:
887 diff = p->cutval - dat->y[n];
888 if (lastp == p->loson) {
889 if (ndist > (double) ((int) diff))
890 node_nearest_work (thetree, dat, wcoord, p->hison, n,
891 &ndist, &nnode);
892 } else {
893 if (ndist > (double) ((int) (-diff)))
894 node_nearest_work (thetree, dat, wcoord, p->loson, n,
895 &ndist, &nnode);
896 }
897 break;
898 case 2:
899 if (lastp == p->loson) {
900 if (ndist > p->cutval + wcoord[n])
901 node_nearest_work (thetree, dat, wcoord, p->hison, n,
902 &ndist, &nnode);
903 } else {
904 node_nearest_work (thetree, dat, wcoord, p->loson, n,
905 &ndist, &nnode);
906 }
907 break;
908 }
909 if (wcoord == (double *) NULL && p->bnds &&
910 ball_in_bounds (dat, p->bnds, n, ndist))
911 break;
912 }
913 return nnode;
914 }
915
916 #ifdef CC_PROTOTYPE_ANSI
ball_in_bounds(CCdatagroup * dat,CCkdbnds * bnds,int n,double dist)917 static int ball_in_bounds (CCdatagroup *dat, CCkdbnds *bnds, int n, double dist)
918 #else
919 static int ball_in_bounds (dat, bnds, n, dist)
920 CCdatagroup *dat;
921 CCkdbnds *bnds;
922 int n;
923 double dist;
924 #endif
925 {
926 if ((double) ((int) (dat->x[n] - bnds->x[0])) < dist ||
927 (double) ((int) (bnds->x[1] - dat->x[n])) < dist ||
928 (double) ((int) (dat->y[n] - bnds->y[0])) < dist ||
929 (double) ((int) (bnds->y[1] - dat->y[n])) < dist)
930 return 0;
931 return 1;
932 }
933
934 #ifdef CC_PROTOTYPE_ANSI
node_nearest_work(CCkdtree * thetree,CCdatagroup * dat,double * datw,CCkdnode * p,int target,double * ndist,int * nnode)935 static void node_nearest_work (CCkdtree *thetree, CCdatagroup *dat,
936 double *datw, CCkdnode *p, int target, double *ndist, int *nnode)
937 #else
938 static void node_nearest_work (thetree, dat, datw, p, target, ndist, nnode)
939 CCkdtree *thetree;
940 CCdatagroup *dat;
941 double *datw;
942 CCkdnode *p;
943 int target;
944 double *ndist;
945 int *nnode;
946 #endif
947 {
948 int i;
949 double val, thisx, thisdist;
950
951 if (!p->empty) {
952 if (p->bucket) {
953 for (i = p->lopt; i <= p->hipt; i++) {
954 if (thetree->perm[i] != target) {
955 thisdist = Fedgelen (thetree->perm[i], target);
956 if (*ndist > thisdist) {
957 *ndist = thisdist;
958 *nnode = thetree->perm[i];
959 }
960 }
961 }
962 } else {
963 val = p->cutval;
964 switch (p->cutdim) {
965 case 0:
966 thisx = dat->x[target];
967 if (thisx < val) {
968 node_nearest_work (thetree, dat, datw, p->loson, target,
969 ndist, nnode);
970 if (*ndist > (double) ((int) (val - thisx)))
971 node_nearest_work (thetree, dat, datw, p->hison,
972 target, ndist, nnode);
973 } else {
974 node_nearest_work (thetree, dat, datw, p->hison, target,
975 ndist, nnode);
976 if (*ndist > (double) ((int) (thisx - val)))
977 node_nearest_work (thetree, dat, datw, p->loson,
978 target, ndist, nnode);
979 }
980 break;
981 case 1:
982 thisx = dat->y[target];
983 if (thisx < val) {
984 node_nearest_work (thetree, dat, datw, p->loson, target,
985 ndist, nnode);
986 if (*ndist > (double) ((int) (val - thisx)))
987 node_nearest_work (thetree, dat, datw, p->hison,
988 target, ndist, nnode);
989 } else {
990 node_nearest_work (thetree, dat, datw, p->hison, target,
991 ndist, nnode);
992 if (*ndist > (double) ((int) (thisx - val)))
993 node_nearest_work (thetree, dat, datw, p->loson,
994 target, ndist, nnode);
995 }
996 break;
997 case 2:
998 thisx = datw[target];
999 node_nearest_work (thetree, dat, datw, p->loson, target, ndist,
1000 nnode);
1001 if (*ndist > val + thisx)
1002 node_nearest_work (thetree, dat, datw, p->hison, target,
1003 ndist, nnode);
1004 break;
1005 }
1006 }
1007 }
1008 }
1009
1010 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_fixed_radius_nearest(CCkdtree * kt,CCdatagroup * dat,double * datw,int n,double rad,int (* doit_fn)(int,int,void *),void * pass_param)1011 int CCkdtree_fixed_radius_nearest (CCkdtree *kt, CCdatagroup *dat,
1012 double *datw, int n, double rad, int (*doit_fn) (int, int, void *),
1013 void *pass_param)
1014 #else
1015 int CCkdtree_fixed_radius_nearest (kt, dat, datw, n, rad, doit_fn, pass_param)
1016 CCkdtree *kt;
1017 CCdatagroup *dat;
1018 double *datw;
1019 int n;
1020 double rad;
1021 int (*doit_fn)();
1022 void *pass_param;
1023 #endif
1024 {
1025 CCkdnode *p, *lastp;
1026 double dist, diff, xtarget, ytarget;
1027 int target;
1028
1029 if (kt == (CCkdtree *) NULL) {
1030 fprintf (stderr, "ERROR: fixed_radius_nearest needs a CCkdtree\n");
1031 return 0;
1032 }
1033
1034 dist = (double) rad;
1035 target = n;
1036 xtarget = dat->x[target];
1037 ytarget = dat->y[target];
1038
1039 p = kt->bucketptr[target];
1040 if (fixed_radius_nearest_work (kt, p, doit_fn, target, dist, dat, datw,
1041 xtarget, ytarget, pass_param))
1042 return 1;
1043
1044 if (datw) {
1045 double wdist = dist - datw[target];
1046 while (1) {
1047 lastp = p;
1048 p = p->father;
1049 if (p == (CCkdnode *) NULL)
1050 return 0;
1051 if (p->cutdim == 0)
1052 diff = p->cutval - xtarget;
1053 else if (p->cutdim == 1)
1054 diff = p->cutval - ytarget;
1055 else
1056 diff = p->cutval;
1057 if (lastp == p->loson) {
1058 if (wdist > (double) ((int) diff)) {
1059 if (fixed_radius_nearest_work (kt, p->hison, doit_fn,
1060 target, dist, dat, datw, xtarget, ytarget,
1061 pass_param))
1062 return 1;
1063 }
1064 } else {
1065 if (wdist > (double) ((int) -diff)) {
1066 if (fixed_radius_nearest_work (kt, p->loson, doit_fn,
1067 target, dist, dat, datw, xtarget, ytarget,
1068 pass_param))
1069 return 1;
1070 }
1071 }
1072 if (p->bnds && /* ball_in_bounds */
1073 !((double) ((int) (xtarget - p->bnds->x[0])) < wdist ||
1074 (double) ((int) (p->bnds->x[1] - xtarget)) < wdist ||
1075 (double) ((int) (ytarget - p->bnds->y[0])) < wdist ||
1076 (double) ((int) (p->bnds->y[1] - ytarget)) < wdist))
1077 return 0;
1078 }
1079 } else {
1080 while (1) {
1081 lastp = p;
1082 p = p->father;
1083 if (p == (CCkdnode *) NULL)
1084 return 0;
1085 if (p->cutdim == 0)
1086 diff = p->cutval - xtarget;
1087 else
1088 diff = p->cutval - ytarget;
1089
1090 if (lastp == p->loson) {
1091 if (dist > (double) ((int) diff)) {
1092 if (fixed_radius_nearest_work (kt, p->hison, doit_fn,
1093 target, dist, dat, datw, xtarget, ytarget,
1094 pass_param))
1095 return 1;
1096 }
1097 } else {
1098 if (dist > (double) ((int) -diff) || p->cutdim == 2) {
1099 if (fixed_radius_nearest_work (kt, p->loson, doit_fn,
1100 target, dist, dat, datw, xtarget, ytarget,
1101 pass_param))
1102 return 1;
1103 }
1104 }
1105 if (p->bnds && /* ball_in_bounds */
1106 !((double) ((int) (xtarget - p->bnds->x[0])) < dist ||
1107 (double) ((int) (p->bnds->x[1] - xtarget)) < dist ||
1108 (double) ((int) (ytarget - p->bnds->y[0])) < dist ||
1109 (double) ((int) (p->bnds->y[1] - ytarget)) < dist))
1110 return 0;
1111 }
1112 }
1113 }
1114
1115 #ifdef CC_PROTOTYPE_ANSI
fixed_radius_nearest_work(CCkdtree * thetree,CCkdnode * p,int (* doit_fn)(int,int,void *),int target,double dist,CCdatagroup * dat,double * datw,double xtarget,double ytarget,void * pass_param)1116 static int fixed_radius_nearest_work (CCkdtree *thetree, CCkdnode *p,
1117 int (*doit_fn) (int, int, void *), int target, double dist,
1118 CCdatagroup *dat, double *datw, double xtarget, double ytarget,
1119 void *pass_param)
1120 #else
1121 static int fixed_radius_nearest_work (thetree, p, doit_fn, target, dist, dat,
1122 datw, xtarget, ytarget, pass_param)
1123 CCkdtree *thetree;
1124 CCkdnode *p;
1125 int (*doit_fn)();
1126 int target;
1127 double dist;
1128 CCdatagroup *dat;
1129 double *datw;
1130 double xtarget, ytarget;
1131 void *pass_param;
1132 #endif
1133 {
1134 int i, c;
1135 double val, thisx, thisdist;
1136
1137 if (p->empty)
1138 return 0;
1139
1140 if (p->bucket) {
1141 for (i = p->lopt; i <= p->hipt; i++) {
1142 c = thetree->perm[i];
1143 if (c != target) {
1144 thisdist = Fedgelen (c, target);
1145 if (thisdist < dist) {
1146 if (doit_fn (target, c, pass_param)) {
1147 return 1;
1148 }
1149 }
1150 }
1151 }
1152 return 0;
1153 } else {
1154 if (datw) {
1155 double wdist = dist - datw[target];
1156
1157 val = p->cutval;
1158 switch (p->cutdim) {
1159 case 0:
1160 thisx = xtarget;
1161 break;
1162 case 1:
1163 thisx = ytarget;
1164 break;
1165 case 2:
1166 if (fixed_radius_nearest_work (thetree, p->loson, doit_fn,
1167 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1168 return 1;
1169 }
1170 if (p->cutval <= wdist) {
1171 if (fixed_radius_nearest_work (thetree, p->hison, doit_fn,
1172 target, dist, dat, datw, xtarget, ytarget,
1173 pass_param)) {
1174 return 1;
1175 }
1176 }
1177 return 0;
1178 default:
1179 return 0;
1180 }
1181 if (thisx < val) {
1182 if (fixed_radius_nearest_work (thetree, p->loson, doit_fn,
1183 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1184 return 1;
1185 }
1186 if (wdist > (double) ((int) (val - thisx))) {
1187 if (fixed_radius_nearest_work (thetree, p->hison, doit_fn,
1188 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1189 return 1;
1190 }
1191 }
1192 } else {
1193 if (fixed_radius_nearest_work (thetree, p->hison, doit_fn,
1194 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1195 return 1;
1196 }
1197 if (wdist > (double) ((int) (thisx - val))) {
1198 if (fixed_radius_nearest_work (thetree, p->loson, doit_fn,
1199 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1200 return 1;
1201 }
1202 }
1203 }
1204 } else {
1205 val = p->cutval;
1206 switch (p->cutdim) {
1207 case 0:
1208 thisx = xtarget;
1209 break;
1210 case 1:
1211 thisx = ytarget;
1212 break;
1213 case 2:
1214 default:
1215 fprintf (stderr, "ERROR: split on w without node weights\n");
1216 return 0;
1217 }
1218 if (thisx < val) {
1219 if (fixed_radius_nearest_work (thetree, p->loson, doit_fn,
1220 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1221 return 1;
1222 }
1223 if (dist > (double) ((int) (val - thisx))) {
1224 if (fixed_radius_nearest_work (thetree, p->hison, doit_fn,
1225 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1226 return 1;
1227 }
1228 }
1229 } else {
1230 if (fixed_radius_nearest_work (thetree, p->hison, doit_fn,
1231 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1232 return 1;
1233 }
1234 if (dist > (double) ((int) (thisx - val))) {
1235 if (fixed_radius_nearest_work (thetree, p->loson, doit_fn,
1236 target, dist, dat, datw, xtarget, ytarget, pass_param)) {
1237 return 1;
1238 }
1239 }
1240 }
1241 }
1242 }
1243 return 0;
1244 }
1245
1246 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_nearest_neighbor_tour(CCkdtree * kt,int ncount,int start,CCdatagroup * dat,int * outcycle,double * val)1247 int CCkdtree_nearest_neighbor_tour (CCkdtree *kt, int ncount, int start,
1248 CCdatagroup *dat, int *outcycle, double *val)
1249 #else
1250 int CCkdtree_nearest_neighbor_tour (kt, ncount, start, dat, outcycle, val)
1251 CCkdtree *kt;
1252 int ncount, start;
1253 CCdatagroup *dat;
1254 int *outcycle;
1255 double *val;
1256 #endif
1257 {
1258 double len;
1259 int i, current, next;
1260 CCkdtree localkt, *mykt;
1261 int newtree = 0;
1262
1263 if (ncount < 3) {
1264 fprintf (stderr, "Cannot find tour in an %d node graph\n", ncount);
1265 return 1;
1266 }
1267
1268 if (kt == (CCkdtree *) NULL) {
1269 if (CCkdtree_build (&localkt, ncount, dat, (double *) NULL)) {
1270 fprintf (stderr, "Unable to build CCkdtree\n");
1271 return 1;
1272 }
1273 mykt = &localkt;
1274 newtree = 1;
1275 } else {
1276 mykt = kt;
1277 }
1278
1279 /*
1280 printf ("Grow nearest neighbor tour from node %d\n", start);
1281 fflush (stdout);
1282 */
1283
1284 len = 0.0;
1285 current = start;
1286 if (outcycle != (int *) NULL)
1287 outcycle[0] = start;
1288 for (i = 1; i < ncount; i++) {
1289 CCkdtree_delete (mykt, current);
1290 next = CCkdtree_node_nearest (mykt, current, dat, (double *) NULL);
1291 if (outcycle != (int *) NULL)
1292 outcycle [i] = next;
1293 len += (double) CCutil_dat_edgelen (current, next, dat);
1294 current = next;
1295 }
1296 len += (double) CCutil_dat_edgelen (current, start, dat);
1297 *val = len;
1298 if (newtree)
1299 CCkdtree_free (&localkt);
1300 else
1301 CCkdtree_undelete_all (kt, ncount);
1302 return 0;
1303 }
1304
1305 #ifdef CC_PROTOTYPE_ANSI
CCkdtree_nearest_neighbor_2match(CCkdtree * kt,int ncount,int start,CCdatagroup * dat,int * outmatch,double * val)1306 int CCkdtree_nearest_neighbor_2match (CCkdtree *kt, int ncount, int start,
1307 CCdatagroup *dat, int *outmatch, double *val)
1308 #else
1309 int CCkdtree_nearest_neighbor_2match (kt, ncount, start, dat, outmatch, val)
1310 CCkdtree *kt;
1311 int ncount, start;
1312 CCdatagroup *dat;
1313 int *outmatch;
1314 double *val;
1315 #endif
1316 {
1317 double len;
1318 int i, j, cur, next;
1319 CCkdtree localkt, *mykt;
1320 double szeit;
1321 int count = 0, cyccount = 0;
1322 char *mark = (char *) NULL;
1323 int newtree = 0;
1324 int rval = 0;
1325
1326 if (ncount < 3) {
1327 fprintf (stderr, "Cannot find 2-matching in an %d node graph\n",
1328 ncount);
1329 return 1;
1330 }
1331
1332 if (kt == (CCkdtree *) NULL) {
1333 if (CCkdtree_build (&localkt, ncount, dat, (double *) NULL)) {
1334 fprintf (stderr, "Unable to build CCkdtree\n");
1335 return 1;
1336 }
1337 mykt = &localkt;
1338 newtree = 1;
1339 } else {
1340 mykt = kt;
1341 }
1342
1343 mark = CC_SAFE_MALLOC (ncount, char);
1344 if (!mark) {
1345 rval = 1;
1346 goto CLEANUP;
1347 }
1348 for (i = 0 ; i < ncount; i++)
1349 mark[i] = 0;
1350
1351 printf ("Grow nearest neighbor 2-matching from node %d\n", start);
1352 fflush (stdout);
1353 szeit = CCutil_zeit ();
1354 len = 0.0;
1355
1356 while (count < ncount) {
1357 for (j = start; j < ncount && mark[j]; j++);
1358 if (j == ncount) {
1359 for (j = 0; j < start && mark[j]; j++);
1360 if (j == start) {
1361 fprintf (stderr, "ERROR in near-2match\n");
1362 rval = 1;
1363 goto CLEANUP;
1364 }
1365 }
1366 start = j;
1367 mark[start] = 1;
1368 CCkdtree_delete (mykt, start);
1369 next = CCkdtree_node_nearest (mykt, start, dat, (double *) NULL);
1370 mark[next] = 1;
1371 len += (double) CCutil_dat_edgelen (start, next, dat);
1372 if (outmatch != (int *) NULL) {
1373 outmatch[2 * count] = start;
1374 outmatch[(2 * count) + 1] = next;
1375 }
1376 count++;
1377 CCkdtree_delete (mykt, next);
1378 cur = CCkdtree_node_nearest (mykt, next, dat, (double *) NULL);
1379 len += (double) CCutil_dat_edgelen (next, cur, dat);
1380 if (outmatch != (int *) NULL) {
1381 outmatch[2 * count] = next;
1382 outmatch[(2 * count) + 1] = cur;
1383 }
1384 count++;
1385 CCkdtree_undelete (mykt, start);
1386 while (cur != start && count < ncount - 3) {
1387 mark[cur] = 1;
1388 CCkdtree_delete (mykt, cur);
1389 next = CCkdtree_node_nearest (mykt, cur, dat, (double *) NULL);
1390 len += (double) (*CCutil_dat_edgelen) (cur, next, dat);
1391 if (outmatch != (int *) NULL) {
1392 outmatch[2 * count] = cur;
1393 outmatch[(2 * count) + 1] = next;
1394 }
1395 count++;
1396 cur = next;
1397 }
1398 CCkdtree_delete (mykt, start);
1399
1400 if (cur != start) { /* Not enough nodes for another circuit */
1401 while (count < ncount - 1) {
1402 mark[cur] = 1;
1403 CCkdtree_delete (mykt, cur);
1404 next = CCkdtree_node_nearest (mykt, cur, dat, (double *) NULL);
1405 len += (double) CCutil_dat_edgelen (cur, next, dat);
1406 if (outmatch != (int *) NULL) {
1407 outmatch[2 * count] = cur;
1408 outmatch[(2 * count) + 1] = next;
1409 }
1410 count++;
1411 cur = next;
1412 }
1413 len += (double) CCutil_dat_edgelen (cur, start, dat);
1414 if (outmatch != (int *) NULL) {
1415 outmatch[2 * count] = cur;
1416 outmatch[(2 * count) + 1] = start;
1417 }
1418 count++;
1419 }
1420 cyccount++;
1421 }
1422
1423 *val = len;
1424 printf ("%d cycles in 2-matching\n", cyccount);
1425 printf ("Running time for Nearest Neighbor 2-match: %.2f\n",
1426 CCutil_zeit () - szeit);
1427 fflush (stdout);
1428
1429 CLEANUP:
1430
1431 if (newtree)
1432 CCkdtree_free (&localkt);
1433 else
1434 CCkdtree_undelete_all (kt, ncount);
1435 if (mark)
1436 CC_FREE (mark, char);
1437 return rval;
1438 }
1439