The Cuik KD-Tree Library


cuik-kdtree.c
Go to the documentation of this file.
1 #include "cuik-kdtree.h"
2 
3 #include "definitions.h"
4 #include "vector.h"
5 #include "random.h"
6 
7 #include <assert.h>
8 #include <stdio.h>
9 #include <string.h>
10 
45 #define TRUE_MEDIAN 0
46 
55 #define KDTREE_PRINT_INDENT(f,indent) {unsigned int k; for(k=0;k<indent;k++) fprintf(f," ");}
56 
74 TKDtree *BuildKDTree(Trectangle *ambient,unsigned int *topology,
75  unsigned int nd,unsigned int m,double r,
76  unsigned int n,double **points,unsigned int *ids);
77 
91 void SearchKDtree(double dr2,double *p,double *d2,
92  unsigned int *id,TKDtree *kdt);
93 
108 void SearchInBallKDtree(double dr2,double *p,double r2,
109  unsigned int *n,unsigned int *m,
110  unsigned int **id,TKDtree *kdt);
111 
112 #if (PLOTKDT)
113 
127  void PlotKDtreeNode(FILE *f,unsigned int mode, TKDtree *kdt);
128 #endif
129 /************************************************************************/
130 /************************************************************************/
131 /************************************************************************/
132 TKDtree *BuildKDTree(Trectangle *ambient,unsigned int *topology,
133  unsigned int nd,unsigned int m,double r,
134  unsigned int n,double **points,unsigned int *ids)
135 {
136  TKDtree *kdt;
137 
138  NEW(kdt,1,TKDtree);
139 
140  kdt->nd=nd;
141  kdt->n=n;
142  kdt->m=m;
143  kdt->r=r;
144  CopyRectangle(&(kdt->ambient),ambient);
145  if (topology==NULL)
146  kdt->topology=NULL;
147  else
148  {
149  NEW(kdt->topology,kdt->nd,unsigned int);
150  memcpy(kdt->topology,topology,sizeof(unsigned int)*kdt->nd);
151  }
152 
153  if (n<m)
154  {
155  unsigned int i;
156 
157  kdt->isLeaf=TRUE;
158 
159  kdt->height=1;
160 
161  NEW(kdt->p,kdt->m+1,double *); /* +1 is useful when adding points */
162  NEW(kdt->id,kdt->m+1,unsigned int);
163 
164  if (n==0)
165  kdt->volume=0.0;
166  else
167  {
168  for(i=0;i<n;i++)
169  {
170  kdt->p[i]=points[i];
171  kdt->id[i]=ids[i];
172  if (i==0)
173  InitRectangleFromPoint(kdt->nd,points[i],&(kdt->boundingBox));
174  else
175  ExpandRectangle(points[i],&(kdt->boundingBox));
176  }
177  CopyRectangle(&(kdt->samplingBox),&(kdt->boundingBox)); /* Just to init samplingBox */
178  kdt->volume=EnlargeRectangleWithLimits(r,
179  &(kdt->ambient),
180  &(kdt->boundingBox),
181  &(kdt->samplingBox));
182  }
183  kdt->nextLeaf=NULL;
184 
185  #ifndef NDEBUG
186  /* Reset the internal node info */
187  kdt->splitDim=NO_UINT;
188  kdt->splitValue=INF;
189  kdt->left=NULL;
190  kdt->right=NULL;
191  kdt->leftLeaf=NULL;
192  kdt->rightLeaf=NULL;
193  #endif
194  }
195  else
196  {
197  signed int i,j,it;
198  double *ps;
199  unsigned int is;
200  TKDtree *kdtLeft,*kdtRight;
201  double u,l;
202 
203  kdt->isLeaf=FALSE;
204 
205  /* determine the split dimension and split point */
206  kdt->splitDim=GetRectangleSplitDim(&(kdt->ambient));
207  it=0;
208  do {
209  #if (TRUE_MEDIAN)
210  {
211  unsigned int med,right,left;
212  double val;
213 
214  /* We adapt the quickselect algorithm to find the median in O(n) */
215  med=(unsigned int)n/2; //position of the median
216  left=0;
217  right=n;
218  i=0; // do not start the loop if r==l==med
219  while (i!=med)
220  {
221  //we use the leftmost value as a pivot
222  //(can be any index betwee 'left' and 'right')
223  val=points[left][kdt->splitDim];
224 
225  /*
226  // this is only necessary if we pick a pivot value
227  // different from the left-most one
228  if (pivot!=left)
229  {
230  //Save the pivot value in the leftmost position
231  SWAP(points[left],points[i],ps);
232  SWAP(ids[left],ids[i],is);
233  }
234  */
235 
236  // for the rest of values, find the first value
237  // larger than the pivot value
238  i=left+1;
239  while ((i<right)&&(points[i][kdt->splitDim]<val))
240  i++;
241  // classify the values after 'i'
242  j=i+1;
243  while (j<right)
244  {
245  if (points[j][kdt->splitDim]<val)
246  {
247  SWAP(points[i],points[j],ps);
248  SWAP(ids[i],ids[j],is);
249  i++;
250  }
251  j++;
252  }
253  if (left!=i-1)
254  {
255  // insert the pivot just before the 1st value larger than it
256  SWAP(points[left],points[i-1],ps);
257  SWAP(ids[left],ids[i-1],is);
258  }
259  if (i<med)
260  left=i;
261  else
262  {
263  if (i>med)
264  right=i;
265  }
266  }
267  kdt->splitValue=points[i][kdt->splitDim];
268  }
269  #else
270  //approximate the median by the mean
271  kdt->splitValue=0.0;
272  for(i=0;i<n;i++)
273  kdt->splitValue+=points[i][kdt->splitDim];
274  kdt->splitValue/=(double)n;
275 
276 
277  //split the points in two sub-sets according to the split point
278  i=0;
279  j=n-1;
280  do {
281  while ((i<n)&&(points[i][kdt->splitDim]<kdt->splitValue))
282  i++;
283  while ((j>=0)&&(points[j][kdt->splitDim]>=kdt->splitValue))
284  j--;
285  if ((i<n)&&(j>=0)&&(i<j)&&
286  (points[i][kdt->splitDim]>=kdt->splitValue)&&
287  (points[j][kdt->splitDim]<kdt->splitValue))
288  {
289  SWAP(points[i],points[j],ps);
290  SWAP(ids[i],ids[j],is);
291  }
292  } while (i<j);
293  #endif
294 
295  /* at this point we have the split value and the vector of points is
296  sorted with 'i' elements with values below the split value and
297  n-i elements with values above it. */
298  if ((i==0)||(i==n))
299  {
300  /* if we picked a split dimension where all points have the same
301  value, just pick the next one. */
302  kdt->splitDim++;
303  kdt->splitDim%=kdt->nd;
304  /* but there are a limited number of dimensions :) */
305  it++;
306  assert(it<kdt->nd);
307  }
308  } while((i==0)||(i==n));
309 
310  /* Recursively build kd-trees */
311 
312  GetRectangleLimits(kdt->splitDim,&l,&u,&(kdt->ambient));
313 
314  SetRectangleUpperLimit(kdt->splitDim,kdt->splitValue,&(kdt->ambient));
315  kdtLeft=BuildKDTree(&(kdt->ambient),topology,nd,m,r,i,points,ids);
316 
317  SetRectangleLimits(kdt->splitDim,kdt->splitValue,u,&(kdt->ambient));
318  kdtRight=BuildKDTree(&(kdt->ambient),topology,nd,m,r,n-i,&(points[i]),&(ids[i]));
319 
320  SetRectangleLowerLimit(kdt->splitDim,l,&(kdt->ambient));
321 
322  /* Assemble the sub-trees */
323  kdt->left=(void *)kdtLeft;
324  kdt->right=(void *)kdtRight;
325 
326  kdt->leftLeaf=(kdtLeft->isLeaf?(void*)kdtLeft:kdtLeft->leftLeaf);
327  kdt->rightLeaf=(kdtRight->isLeaf?(void*)kdtRight:kdtRight->rightLeaf);
328 
329  /* Update the smapling area */
330  kdt->volume=kdtLeft->volume+kdtRight->volume;
331 
332  /* Update the height */
333  kdt->height=1+(kdtLeft->height>kdtRight->height?
334  kdtLeft->height:
335  kdtRight->height);
336 
337  /* Chain the sub-tree leaves */
338  if (kdtLeft->isLeaf)
339  {
340  if (kdtRight->isLeaf)
341  kdtLeft->nextLeaf=(void *)kdtRight;
342  else
343  kdtLeft->nextLeaf=kdtRight->leftLeaf;
344  }
345  else
346  {
347  if (kdtRight->isLeaf)
348  KDTPTR(kdtLeft->rightLeaf)->nextLeaf=(void *)kdtRight;
349  else
350  KDTPTR(kdtLeft->rightLeaf)->nextLeaf=kdtRight->leftLeaf;
351  }
352 
353  #ifndef NDEBUG
354  /* Reset the leaf info */
355  kdt->nextLeaf=NULL;
356  kdt->p=NULL;
357  kdt->id=NULL;
358  #endif
359  }
360 
361  return(kdt);
362 }
363 
364 void SearchKDtree(double dr2,double *p,double *d2,
365  unsigned int *id,TKDtree *kdt)
366 {
367  if (kdt->isLeaf)
368  {
369  if (SquaredDistanceToRectangle(*d2,p,kdt->topology,&(kdt->boundingBox))<*d2)
370  {
371  unsigned int i;
372  double dp;
373 
374  for(i=0;i<kdt->n;i++)
375  {
376  dp=VectorSquaredDistanceTopologyMin(*d2,kdt->nd,kdt->topology,p,kdt->p[i]);
377  if (dp<*d2)
378  {
379  *d2=dp;
380  *id=kdt->id[i];
381  }
382  }
383  }
384  }
385  else
386  {
387  TKDtree *t1,*t2;
388  double dt0,dt1,dt2,d,v;
389 
390  v=p[kdt->splitDim];
391  if (v<=kdt->splitValue)
392  {
393  /* The nearest point is most likely in the left branch */
394  t1=KDTPTR(kdt->left);
395  t2=KDTPTR(kdt->right);
396  }
397  else
398  {
399  /* The nearest point is most likely in the left branch */
400  t1=KDTPTR(kdt->right);
401  t2=KDTPTR(kdt->left);
402  }
403 
404  dt0=dr2-SquaredDistanceToRectangleDimension(kdt->splitDim,v,kdt->topology,
405  &(kdt->ambient));
406 
407  dt1=SquaredDistanceToRectangleDimension(kdt->splitDim,v,kdt->topology,
408  &(t1->ambient));
409  d=dt0+dt1;
410  if (d<*d2)
411  SearchKDtree(d,p,d2,id,t1);
412 
413  dt2=SquaredDistanceToRectangleDimension(kdt->splitDim,v,kdt->topology,
414  &(t2->ambient));
415  d=dt0+dt2;
416  if (d<*d2)
417  SearchKDtree(d,p,d2,id,t2);
418  }
419 }
420 
421 void SearchInBallKDtree(double dr2,double *p,double r2,
422  unsigned int *n,unsigned int *m,
423  unsigned int **id,TKDtree *kdt)
424 {
425  if (kdt->isLeaf)
426  {
427  if (SquaredDistanceToRectangle(r2,p,kdt->topology,&(kdt->boundingBox))<r2)
428  {
429  unsigned int i;
430 
431  for(i=0;i<kdt->n;i++)
432  {
433  if (VectorSquaredDistanceTopologyMin(r2,kdt->nd,kdt->topology,p,kdt->p[i])<r2)
434  {
435  if (*n==*m)
436  MEM_DUP(*id,*m,unsigned int);
437  (*id)[*n]=kdt->id[i];
438  (*n)++;
439  }
440  }
441  }
442  }
443  else
444  {
445  TKDtree *t;
446  double d,dt,v;
447 
448  v=p[kdt->splitDim];
449  dt=dr2-SquaredDistanceToRectangleDimension(kdt->splitDim,v,kdt->topology,
450  &(kdt->ambient));
451 
452  t=KDTPTR(kdt->left);
453  d=dt+SquaredDistanceToRectangleDimension(kdt->splitDim,v,kdt->topology,
454  &(t->ambient));
455  if (d<r2)
456  SearchInBallKDtree(d,p,r2,n,m,id,t);
457 
458  t=KDTPTR(kdt->right);
459  d=dt+SquaredDistanceToRectangleDimension(kdt->splitDim,v,kdt->topology,
460  &(t->ambient));
461  if (d<r2)
462  SearchInBallKDtree(d,p,r2,n,m,id,t);
463  }
464 }
465 
466 #if (PLOTKDT)
467  void PlotKDtreeNode(FILE *f,unsigned int mode, TKDtree *kdt)
468  {
469  if (kdt->isLeaf)
470  {
471  unsigned int i;
472  double lx,ux,ly,uy,lz,uz;
473 
474  /* Plot the points (in red) */
475  if (mode&1)
476  {
477  for(i=0;i<kdt->n;i++)
478  {
479  fprintf(f,"{VECT 1 2 1 \n 2 \n 1\n");
480  fprintf(f," %f %f %f\n",kdt->p[i][0]-0.01,kdt->p[i][1],kdt->p[i][2]);
481  fprintf(f," %f %f %f\n",kdt->p[i][0]+0.01,kdt->p[i][1],kdt->p[i][2]);
482  fprintf(f," 1 0 0 0}\n");
483  fprintf(f,"{VECT 1 2 1 \n 2 \n 1\n");
484  fprintf(f," %f %f %f\n",kdt->p[i][0],kdt->p[i][1]-0.01,kdt->p[i][2]);
485  fprintf(f," %f %f %f\n",kdt->p[i][0],kdt->p[i][1]+0.01,kdt->p[i][2]);
486  fprintf(f," 1 0 0 0}\n");
487  fprintf(f,"{VECT 1 2 1 \n 2 \n 1\n");
488  fprintf(f," %f %f %f\n",kdt->p[i][0],kdt->p[i][1],kdt->p[i][2]-0.01);
489  fprintf(f," %f %f %f\n",kdt->p[i][0],kdt->p[i][1],kdt->p[i][2]+0.01);
490  fprintf(f," 1 0 0 0}\n");
491  }
492  }
493 
494  if (mode&12)
495  {
496  Trectangle *r;
497 
498  if (mode&2)
499  r=&(kdt->samplingBox); /* Plot the sampling box */
500  else
501  {
502  if (mode&4)
503  r=&(kdt->boundingBox); /* Plot the bouding box */
504  else
505  r=&(kdt->ambient); /* Plot the ambient box */
506  }
507 
508  GetRectangleLimits(0,&lx,&ux,r);
509  GetRectangleLimits(1,&ly,&uy,r);
510  GetRectangleLimits(2,&lz,&uz,r);
511 
512  fprintf(f,"{OFF 8 6 1\n");
513  fprintf(f,"%f %f %f\n",lx,ly,lz);
514  fprintf(f,"%f %f %f\n",lx,uy,lz);
515  fprintf(f,"%f %f %f\n",ux,uy,lz);
516  fprintf(f,"%f %f %f\n",ux,ly,lz);
517 
518  fprintf(f,"%f %f %f\n",lx,ly,uz);
519  fprintf(f,"%f %f %f\n",lx,uy,uz);
520  fprintf(f,"%f %f %f\n",ux,uy,uz);
521  fprintf(f,"%f %f %f\n",ux,ly,uz);
522 
523  /*faces*/
524  fprintf(f,"4 3 2 1 0\n");
525  fprintf(f,"4 4 5 6 7\n");
526  fprintf(f,"4 0 1 5 4\n");
527  fprintf(f,"4 2 3 7 6\n");
528  fprintf(f,"4 3 0 4 7\n");
529  fprintf(f,"4 1 2 6 5}\n");
530  }
531  }
532  else
533  {
534  PlotKDtreeNode(f,mode,KDTPTR(kdt->left));
535  PlotKDtreeNode(f,mode,KDTPTR(kdt->right));
536  }
537  }
538 #endif
539 
540 /************************************************************************/
541 /************************************************************************/
542 /************************************************************************/
543 
544 TKDtree *InitKDTree(Trectangle *ambient,unsigned int *topology,unsigned int m,double r,
545  unsigned int n,double **points,unsigned int *ids)
546 {
547  unsigned int nd;
548 
549  assert(m>=2);
550 
551  /* Correct the ambient box using the topology */
552  nd=GetRectangleDim(ambient);
553  if (topology!=NULL)
554  {
555  unsigned int i;
556  double l,u;
557 
558  for(i=0;i<nd;i++)
559  {
560  if (topology[i]==TOPOLOGY_S)
561  {
562  GetRectangleLimits(i,&l,&u,ambient);
563 
564  assert(l<=u);
565 
566  /* If the given range is close to -Pi / Pi but
567  for rounding errors, we fix the values to avoid
568  numerical errors. */
569  if (fabs(l+M_PI)<ANGLE_ACCURACY)
570  l=-M_PI;
571  if (fabs(u-M_PI)<ANGLE_ACCURACY)
572  u=M_PI;
573 
574  SetRectangleLimits(i,l,u,ambient);
575  }
576  }
577  }
578 
579  return(BuildKDTree(ambient,topology,nd,m,r,n,points,ids));
580 }
581 
582 void AddPoint2KDtree(double *point,unsigned int id,TKDtree **kdt)
583 {
584  TKDtree *t;
585 
586  t=*kdt;
587 
588  if (t->isLeaf)
589  {
590  #ifndef NDEBUG
591  /* Check for the presense of repeated points in the tree */
592  unsigned int j;
593  double e;
594 
595  for(j=0;j<t->n;j++)
596  {
597  e=VectorSquaredDistanceTopologyMin(INF,t->nd,t->topology,t->p[j],point);
598  if (e<1e-12) /*threshold 1e-6*1e-6*/
599  {
600  fprintf(stderr,"Repeated point in kd-tree\n");
601  exit(0);
602  }
603  }
604  #endif
605 
606  t->p[t->n]=point;
607  t->id[t->n]=id;
608  t->n++;
609  if (t->n>t->m)
610  {
611  TKDtree *nt;
612 
613  nt=BuildKDTree(&(t->ambient),t->topology,t->nd,t->m,t->r,t->n,t->p,t->id);
614  DeleteKDtree(t);
615  *kdt=nt;
616  }
617  else
618  {
619  if (t->n==1)
620  {
621  /* adding the first point to an empty tree/leaf */
622  InitRectangleFromPoint(t->nd,point,&(t->boundingBox));
623  CopyRectangle(&(t->samplingBox),&(t->boundingBox)); /* just to allocate space */
624  }
625  else
626  ExpandRectangle(point,&(t->boundingBox));
627  t->volume=EnlargeRectangleWithLimits(t->r,
628  &(t->ambient),
629  &(t->boundingBox),
630  &(t->samplingBox));
631  }
632  }
633  else
634  {
635  TKDtree *kdtLeft,*kdtRight;
636 
637  if (point[t->splitDim]<=t->splitValue)
638  AddPoint2KDtree(point,id,(TKDtree **)&(t->left));
639  else
640  AddPoint2KDtree(point,id,(TKDtree **)&(t->right));
641 
642  t->n++;
643 
644  kdtLeft=KDTPTR(t->left);
645  kdtRight=KDTPTR(t->right);
646 
647  /* One of the sub-trees changed -> update the summary at this level */
648 
649  /* The chain of leaves must be updated even if we re-balance.
650  Before re-balancing, we need a correct chain of leaves. */
651 
652  t->leftLeaf=(kdtLeft->isLeaf?(void*)kdtLeft:kdtLeft->leftLeaf);
653  t->rightLeaf=(kdtRight->isLeaf?(void*)kdtRight:kdtRight->rightLeaf);
654 
655  /* Chain the sub-tree leaves */
656  if (kdtLeft->isLeaf)
657  {
658  if (kdtRight->isLeaf)
659  kdtLeft->nextLeaf=(void *)kdtRight;
660  else
661  kdtLeft->nextLeaf=kdtRight->leftLeaf;
662  }
663  else
664  {
665  if (kdtRight->isLeaf)
666  KDTPTR(kdtLeft->rightLeaf)->nextLeaf=(void *)kdtRight;
667  else
668  KDTPTR(kdtLeft->rightLeaf)->nextLeaf=kdtRight->leftLeaf;
669  }
670 
671  /* If necessary, rebalance the tree -> re-build the tree */
672  if ((KDTPTR(t->left)->height>KDTPTR(t->right)->height*2)||
673  (KDTPTR(t->right)->height>KDTPTR(t->left)->height*2))
674  {
675  double **p;
676  unsigned int *id;
677  TKDtree *l,*nt;
678  unsigned int i,j;
679 
680  /* Collect all the points in the leaves of the tree to rebuild */
681  NEW(p,t->n,double *);
682  NEW(id,t->n,unsigned int);
683  l=KDTPTR(t->leftLeaf);
684  i=0;
685  while(i<t->n)
686  {
687  for(j=0;j<l->n;j++)
688  {
689  p[i]=l->p[j];
690  id[i]=l->id[j];
691  i++;
692  }
693  l=KDTPTR(l->nextLeaf);
694  }
695 
696  /* Re-build the tree */
697  nt=BuildKDTree(&(t->ambient),t->topology,t->nd,t->m,t->r,t->n,p,id);
698  DeleteKDtree(t);
699  *kdt=nt;
700 
701  free(p);
702  free(id);
703  }
704  else
705  {
706  /* Some "summary" information is only valid if we do not re-balance */
707  t->volume=kdtLeft->volume+kdtRight->volume;
708  t->height=1+(kdtLeft->height>kdtRight->height?
709  kdtLeft->height:kdtRight->height);
710  }
711  }
712 }
713 
714 unsigned int NearestNeighbour(double *point,TKDtree *kdt)
715 {
716  unsigned int id;
717  double d;
718 
719  assert(kdt->n>0);
720 
721  id=NO_UINT;
722  d=INF;
723  SearchKDtree(0.0,point,&d,&id,kdt);
724 
725  assert(id!=NO_UINT);
726 
727  return(id);
728 }
729 
730 void NeighboursInBall(double *point,double r,
731  unsigned int *n,unsigned int **ids,TKDtree *kdt)
732 {
733  unsigned int m;
734 
735  assert(kdt->n>0);
736 
737  m=100;
738  NEW(*ids,m,unsigned int);
739  *n=0;
740 
741  SearchInBallKDtree(0.0,point,r*r,n,&m,ids,kdt);
742 
743  if (*n==0)
744  free(*ids);
745 }
746 
747 void SampleInKDtree(double *point,TKDtree *kdt)
748 {
749  if (kdt->isLeaf)
750  RandomPointInRectangle(point,&(kdt->samplingBox));
751  else
752  {
753  double r;
754 
755  r=randomDouble()*kdt->volume;
756 
757  if (r<KDTPTR(kdt->left)->volume)
758  SampleInKDtree(point,KDTPTR(kdt->left));
759  else
760  SampleInKDtree(point,KDTPTR(kdt->right));
761  }
762 }
763 
764 double KDtreeSamplingExpansion(TKDtree *kdt)
765 {
766  return(kdt->r);
767 }
768 
769 void PrintKDtree(FILE *f,unsigned int indent,TKDtree *kdt)
770 {
771  if (kdt->isLeaf)
772  {
773  unsigned int i,j;
774 
775  KDTREE_PRINT_INDENT(f,indent);PrintRectangle(f,&(kdt->boundingBox));fprintf(f,"\n");
776  for(i=0;i<kdt->n;i++)
777  {
778  KDTREE_PRINT_INDENT(f,indent);
779  fprintf(f,"%u: [",kdt->id[i]);
780  for(j=0;j<kdt->nd;j++)
781  fprintf(f,"%f ",kdt->p[i][j]);
782  fprintf(f,"]\n");
783  }
784  }
785  else
786  {
787  KDTREE_PRINT_INDENT(f,indent);fprintf(f,"Split dim: %u\n",kdt->splitDim);
788  KDTREE_PRINT_INDENT(f,indent);fprintf(f,"Split val: %f\n",kdt->splitValue);
789  KDTREE_PRINT_INDENT(f,indent);fprintf(f,"N. points: %u\n",kdt->n);
790  KDTREE_PRINT_INDENT(f,indent);fprintf(f,"Left tree:\n");
791  PrintKDtree(f,indent+1,KDTPTR(kdt->left));
792  KDTREE_PRINT_INDENT(f,indent);fprintf(f,"Right tree:\n");
793  PrintKDtree(f,indent+1,KDTPTR(kdt->right));
794  }
795 }
796 
797 #if (PLOTKDT)
798  void PlotKDtree(char *fname,boolean sphere,TKDtree *kdt)
799  {
800  FILE *f;
801 
802  assert(kdt->nd==3);
803 
804  f=fopen(fname,"w");
805  assert(f);
806 
807  fprintf(f,"(progn (normalization World none)\n(bbox-draw World no))\n");
808 
809  if (sphere)
810  fprintf(f,"(geometry ball {LIST {SPHERE 1 0 0 0}})\n");
811 
812  //points
813  fprintf(f,"(geometry kdtree-nodes {LIST \n");
814  PlotKDtreeNode(f,1,kdt);
815  fprintf(f,"})\n");
816 
817  //sampling volume
818  fprintf(f,"(geometry kdtree-sampling {LIST \n");
819  PlotKDtreeNode(f,2,kdt);
820  fprintf(f,"})\n");
821 
822  /*
823  //bouding box
824  fprintf(f,"(geometry kdtree-bounding {LIST \n");
825  PlotKDtreeNode(f,4,kdt);
826  fprintf(f,"})\n");
827  */
828 
829  fclose(f);
830  }
831 #endif
832 
833 void DeleteKDtree(TKDtree *kdt)
834 {
835  if (kdt->isLeaf)
836  {
837  free(kdt->p);
838  free(kdt->id);
839  if (kdt->n>0)
840  {
841  DeleteRectangle(&(kdt->boundingBox));
842  DeleteRectangle(&(kdt->samplingBox));
843  }
844  }
845  else
846  {
847  DeleteKDtree(KDTPTR(kdt->left));
848  DeleteKDtree(KDTPTR(kdt->right));
849  }
850  DeleteRectangle(&(kdt->ambient));
851  if (kdt->topology!=NULL)
852  free(kdt->topology);
853  free(kdt);
854 }
KDTPTR
#define KDTPTR(a)
Void pointer to KDtree pointer.
Definition: cuik-kdtree.h:38
KDtreeSamplingExpansion
double KDtreeSamplingExpansion(TKDtree *kdt)
Returns the sampling expansion used in the kd-tree.
Definition: cuik-kdtree.c:764
RandomPointInRectangle
void RandomPointInRectangle(double *c, Trectangle *b)
Returns the a random point along the selected dimensions.
Definition: rectangle.c:145
GetRectangleSplitDim
unsigned int GetRectangleSplitDim(Trectangle *b)
Computes the rectangle dimension for which it is better to split the rectangle.
Definition: rectangle.c:344
GetRectangleLimits
void GetRectangleLimits(unsigned int i, double *l, double *u, Trectangle *b)
Gets the limits of the rectangle along a given dimension.
Definition: rectangle.c:115
SampleInKDtree
void SampleInKDtree(double *point, TKDtree *kdt)
Generates a random sample.
Definition: cuik-kdtree.c:747
TKDtree::boundingBox
Trectangle boundingBox
Definition: cuik-kdtree.h:77
TKDtree::volume
double volume
Definition: cuik-kdtree.h:63
SquaredDistanceToRectangleDimension
double SquaredDistanceToRectangleDimension(unsigned int dim, double p, unsigned int *tp, Trectangle *b)
Squared distance from a value to a given rectangle dimension.
Definition: rectangle.c:153
TRUE
#define TRUE
TRUE.
Definition: definitions.h:161
TKDtree::splitValue
double splitValue
Definition: cuik-kdtree.h:67
TKDtree::samplingBox
Trectangle samplingBox
Definition: cuik-kdtree.h:78
InitRectangleFromPoint
void InitRectangleFromPoint(unsigned int dim, double *p, Trectangle *b)
Initializes a rectangle from a point.
Definition: rectangle.c:35
PlotKDtree
void PlotKDtree(char *fname, boolean sphere, TKDtree *kdt)
Plots a kd-tree.
SetRectangleLowerLimit
void SetRectangleLowerLimit(unsigned int i, double l, Trectangle *b)
Set the lower limit.
Definition: rectangle.c:123
TKDtree::height
unsigned int height
Definition: cuik-kdtree.h:62
CopyRectangle
void CopyRectangle(Trectangle *b_out, Trectangle *b_in)
Rectangle copy operator.
Definition: rectangle.c:91
TKDtree::nextLeaf
void * nextLeaf
Definition: cuik-kdtree.h:74
randomDouble
double randomDouble()
Returns a random double in the [0,1] interval.
Definition: random.c:33
TKDtree::r
double r
Definition: cuik-kdtree.h:59
SquaredDistanceToRectangle
double SquaredDistanceToRectangle(double t2, double *p, unsigned int *tp, Trectangle *b)
Squared distance from a point to a rectangle.
Definition: rectangle.c:230
ExpandRectangle
void ExpandRectangle(double *p, Trectangle *b)
Expands a rectangle so that it includes a given point.
Definition: rectangle.c:52
TKDtree::id
unsigned int * id
Definition: cuik-kdtree.h:76
SearchKDtree
void SearchKDtree(double dr2, double *p, double *d2, unsigned int *id, TKDtree *kdt)
Searches the nearest point on a KD-tree.
Definition: cuik-kdtree.c:364
KDTREE_PRINT_INDENT
#define KDTREE_PRINT_INDENT(f, indent)
Auxiliary macro to print the kd-tree.
Definition: cuik-kdtree.c:55
DeleteRectangle
void DeleteRectangle(Trectangle *b)
Destructor.
Definition: rectangle.c:375
DeleteKDtree
void DeleteKDtree(TKDtree *kdt)
Deletes a KD-tree.
Definition: cuik-kdtree.c:833
TKDtree::right
void * right
Definition: cuik-kdtree.h:69
SWAP
#define SWAP(a, b, c)
Swaps two values.
Definition: definitions.h:60
INF
#define INF
Infinite.
Definition: definitions.h:141
NEW
#define NEW(_var, _n, _type)
Allocates memory space.
Definition: definitions.h:87
TKDtree::m
unsigned int m
Definition: cuik-kdtree.h:58
definitions.h
Definition of constants and macros used in several parts of the library.
cuik-kdtree.h
Headers of the kd-tree implementation.
SetRectangleUpperLimit
void SetRectangleUpperLimit(unsigned int i, double u, Trectangle *b)
Set the upper limit.
Definition: rectangle.c:130
Trectangle
A rectangle.
Definition: rectangle.h:26
VectorSquaredDistanceTopologyMin
double VectorSquaredDistanceTopologyMin(double t2, unsigned int s, unsigned int *tp, double *v1, double *v2)
Computes the squared distance of two points.
Definition: vector.c:17
TKDtree::rightLeaf
void * rightLeaf
Definition: cuik-kdtree.h:71
TOPOLOGY_S
#define TOPOLOGY_S
One of the possible topologies.
Definition: definitions.h:188
TKDtree::p
double ** p
Definition: cuik-kdtree.h:75
BuildKDTree
TKDtree * BuildKDTree(Trectangle *ambient, unsigned int *topology, unsigned int nd, unsigned int m, double r, unsigned int n, double **points, unsigned int *ids)
Initializes a kd-tree from a set of points.
Definition: cuik-kdtree.c:132
vector.h
Definition of basic operations between points.
TKDtree::topology
unsigned int * topology
Definition: cuik-kdtree.h:61
TKDtree::isLeaf
boolean isLeaf
Definition: cuik-kdtree.h:53
TKDtree::left
void * left
Definition: cuik-kdtree.h:68
FALSE
#define FALSE
FALSE.
Definition: definitions.h:170
PrintRectangle
void PrintRectangle(FILE *f, Trectangle *b)
Prints a rectangle.
Definition: rectangle.c:365
NO_UINT
#define NO_UINT
Used to denote an identifier that has not been initialized.
Definition: definitions.h:133
TKDtree::nd
unsigned int nd
Definition: cuik-kdtree.h:56
TKDtree::leftLeaf
void * leftLeaf
Definition: cuik-kdtree.h:70
TKDtree::n
unsigned int n
Definition: cuik-kdtree.h:57
TKDtree::splitDim
unsigned int splitDim
Definition: cuik-kdtree.h:66
random.h
Definition of basic randomization functions.
PrintKDtree
void PrintKDtree(FILE *f, unsigned int indent, TKDtree *kdt)
Prints a kd-tree.
Definition: cuik-kdtree.c:769
NearestNeighbour
unsigned int NearestNeighbour(double *point, TKDtree *kdt)
Determines the nearest-neighbour for a point.
Definition: cuik-kdtree.c:714
NeighboursInBall
void NeighboursInBall(double *point, double r, unsigned int *n, unsigned int **ids, TKDtree *kdt)
Neighbours in a ball.
Definition: cuik-kdtree.c:730
M_PI
#define M_PI
Pi.
Definition: definitions.h:24
GetRectangleDim
unsigned int GetRectangleDim(Trectangle *b)
Returns the dimension of the rectangle.
Definition: rectangle.c:47
InitKDTree
TKDtree * InitKDTree(Trectangle *ambient, unsigned int *topology, unsigned int m, double r, unsigned int n, double **points, unsigned int *ids)
Initializes a kd-tree from a set of points.
Definition: cuik-kdtree.c:544
SetRectangleLimits
void SetRectangleLimits(unsigned int i, double l, double u, Trectangle *b)
Changes a rectangle along a given dimension.
Definition: rectangle.c:137
AddPoint2KDtree
void AddPoint2KDtree(double *point, unsigned int id, TKDtree **kdt)
Adds a point to a kd-tree.
Definition: cuik-kdtree.c:582
MEM_DUP
#define MEM_DUP(_var, _n, _type)
Duplicates a previously allocated memory space.
Definition: definitions.h:123
SearchInBallKDtree
void SearchInBallKDtree(double dr2, double *p, double r2, unsigned int *n, unsigned int *m, unsigned int **id, TKDtree *kdt)
Searches the near points on a KD-tree.
Definition: cuik-kdtree.c:421
TKDtree
A KD-tree.
Definition: cuik-kdtree.h:52
ANGLE_ACCURACY
#define ANGLE_ACCURACY
Accuracy used to compare angles.
Definition: definitions.h:44
TKDtree::ambient
Trectangle ambient
Definition: cuik-kdtree.h:60
EnlargeRectangleWithLimits
double EnlargeRectangleWithLimits(double r, Trectangle *limits, Trectangle *bIn, Trectangle *bOut)
Enlarges a box remaining in a given limits.
Definition: rectangle.c:68