kd_build.c

来自「NIST Handwriting OCR Testbed」· C语言 代码 · 共 310 行

/*
# proc: build_tree - create a tree and insert a specified number of floating
# proc:              point vectors of specified class and dimensionality
# proc:              into a tree, and return it.
# proc: tree_split - recursivly distribute the points in a tree into subtrees
# proc:              and attach them as children to the tree
# proc: tree_findsplitvalues - indentify a component to control the division
# proc:              of a tree, and call an assignment subroutine
# proc: tree_entiles - assign suitable ranges for the coefficient on
# proc:              which the tree's points are subdivided
# proc: tree_onevariance - compute the variance of a specified component of a
# proc:              tree's point vectors.
# proc: addtotree - add one point to a tree
# proc:
# proc: init_tree - initialize the elements of a tree
# proc:
# proc: tree_countleaves - count the leaves contained by a tree
# proc:
# proc: tree_whichinterval - count the leaves contained by a tree
# proc:
# proc: tree_order - reorder floating vectors pointed to by a tree into
# proc:		     a contiguous memory block, to avoid cache misses.
# proc: tree_free - release memory allocated and hel by a tree
# proc:
# proc: istreechildless - macro (used to be a function) that indicates whether
# proc:          tree is a node (which has children) or a leaf (which does not).
# proc:          a tree has children or points but not both.
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <values.h>
#include <kdtree.h>

static int nchild = 3;
static minpointspernode = 14;

TREE *build_tree(points, class, npoints, ni)
float *points; int *class, npoints, ni;
{
int i;
TREE *tree;
POINT *pts;

   if ((tree = (TREE *)malloc(sizeof(TREE))) == NULL)
      syserr("build_tree", "malloc", "space for a TREE");
   if ((pts = (POINT *)malloc(npoints * sizeof(POINT))) == NULL)
      syserr("build_tree", "malloc", "space for points TREE");

   init_tree(tree, (float *)points, npoints, ni);
   for ( i = 0 ; i < npoints ; i++, points += ni )
      pts[i].data = points,		/* just an address somewhere	*/
      pts[i].class = class[i];		/* integer class of the point	*/

   for ( i = 0 ; i < npoints ; i++ )
      addtotree(tree, &pts[i]);

   free(pts);
   tree_split(tree, minpointspernode, nchild);
   tree_order(tree);
   return tree;
}


void tree_split(tree, minpointspernode, nchild)
TREE *tree; int minpointspernode, nchild;
{
int i, oldc, allinsameinterval, c, d;

   if (tree->npoints <= minpointspernode)
      return;

   tree->nchild = nchild;
   if ((tree->children = (TREE **)malloc(nchild * sizeof(TREE *))) == NULL)
      syserr("tree_split", "malloc", "space for child pointers");

   for ( c = 0 ; c < nchild ; c++ )
   {
      if ((tree->children[c] = (TREE *)malloc(sizeof(TREE))) == NULL)
         syserr("tree_split", "malloc", "space for a child");
      init_tree(tree->children[c],
         tree->pointsdata, minpointspernode, tree->datadimension);
   }

   tree_findsplitvalues(tree);

   /* find which child this point should be passed down to, ie. look	*/
   /* for the interval in which the value lies				*/
   for ( i = 0, d = tree->splitindex ; i < tree->npoints ; i++ )
   {
      c = tree_whichinterval(tree, tree->points[i].data[d]);
      allinsameinterval = ( i == 0 || (allinsameinterval && oldc == c) );
      oldc = c;
   }

   if (allinsameinterval)
   {
      for ( c = 0 ; c < nchild ; c++ )
         free(tree->children[c]);
      free(tree->children);
      tree->children = NULL;
      tree->nchild = 0;
   }
   else
   {
      for ( i = 0, d = tree->splitindex ; i < tree->npoints ; i++ )
      {
         c = tree_whichinterval(tree, tree->points[i].data[d]);
         addtotree(tree->children[c], &tree->points[i]);
      }

      for ( c = 0 ; c < nchild ; c++ )
         tree_split(tree->children[c], minpointspernode, nchild);

      tree->allocpoints = 0; /* indicate that there is no allocated space, do not	*/
      free(tree->points);	  /* set npoints=0 since that is the number of points	*/
      tree->points = NULL;	  /* contained here PLUS in the children		*/
   }
}


void tree_findsplitvalues(tree)
TREE *tree;
{
int d;
float vari, maxvari = 0.0;

   for ( d = 0 ; d < tree->datadimension ; d++ )
      if (maxvari < (vari = tree_onevariance(tree, d)))
      {
         maxvari = vari;
         tree->splitindex = d;
      }

   tree_entiles(tree, tree->splitindex);
}


float tree_onevariance(tree, d)
TREE *tree; int d;
{
int i;
float vari = 0.0, mean = 0.0;
float *point;

   for ( i = 0 ; i < tree->npoints; i++ )
   {
      point = tree->points[i].data;
      mean += point[d];
      vari += point[d] * point[d];
   }

   mean /= (float)tree->npoints;
   vari /= (float)tree->npoints;
   vari -= mean * mean;
   return vari;
}

int comv(a, b)
float *a, *b;
{
float x = *a, y = *b;
   return (x < y) ? -1 : ((x > y) ? 1 : 0);
}

void tree_entiles(tree, d)
TREE *tree; int d;
{
int i, bin;
float *vals;

   bin = tree->npoints / tree->nchild;
   malloc_flt(&vals, tree->npoints, "tree_entiles");

   for ( i = 0 ; i < tree->npoints ; i++ )
      vals[i] = tree->points[i].data[d];

   qsort((void *)vals, (size_t)tree->npoints, sizeof(float), comv);

   tree->children[0]->vmin = -MAXFLOAT;
   for ( i = 1 ; i < tree->nchild ; i++ )
      tree->children[i]->vmin = tree->children[i-1]->vmax = vals[i*bin];
   tree->children[i-1]->vmax = MAXFLOAT;

   free(vals);
}



void addtotree(tree, point)
TREE *tree; POINT *point;
{
   if (tree->npoints >= tree->allocpoints)
   {
      tree->allocpoints += 30;
      tree->points = (tree->points == NULL) ?
          (POINT *)malloc(tree->allocpoints * sizeof(POINT)) :
          (POINT *)realloc(tree->points, tree->allocpoints * sizeof(POINT));

      if (tree->points == NULL)
         syserr("addtotree", "realloc", "space for 30 more POINTs");
   }
   tree->points[tree->npoints] = *point;
   tree->npoints++;
}

void init_tree(tree, points, npoints, nidim)
TREE *tree; int npoints, nidim; float *points;
{
    tree->datadimension = nidim;
    tree->splitindex = -1;
    tree->vmin = MAXFLOAT; tree->vmax = -MAXFLOAT; /* sic! wrong way round */

    tree->pointsdata = points;
    tree->points = NULL;
    tree->npoints = 0;
    tree->allocpoints = 0;

    tree->nchild = 0;
    tree->children = NULL;
}

int tree_countleaves(tree)
TREE *tree;
{
int i, count = 0;

   if (istreechildless(tree))
      return 1;

   for ( i = 0 ; i < tree->nchild ; i++ )
      count += tree_countleaves(tree->children[i]);

   return count;
}

/* actually get under the hood and copy the pointed to floating vectors */
/* (not just the pointers to them) to the local leaf, this allows 	*/
/* contiguous access by the machine, keeping things in the cache and	*/
/* avoiding paging, which is useful for CPU's to get speed. icky but...	*/
void tree_order(tree)
TREE *tree;
{
int i;

   if (istreechildless(tree))
   {
      if (tree->npoints > 0)
      {
         float *dat, *dst;
         POINT *pts;

         if ((pts = (POINT *)malloc(tree->npoints * sizeof(POINT))) == NULL)
            syserr("tree_order", "malloc", "space for POINTs");
         malloc_flt(&dat, tree->datadimension * tree->npoints, "tree_order");

         for ( i = 0, dst = dat ; i < tree->npoints ; i++, dst += tree->datadimension )
         {
            pts[i].class = tree->points[i].class;
            memcpy(dst, tree->points[i].data, tree->datadimension * sizeof(float));
            pts[i].data  = dst;
         }

         free(tree->points);
         tree->pointsdata = dat;
         tree->points = pts;
         tree->allocpoints = tree->npoints;
      }
   }
   else
      for ( i = 0 ; i < tree->nchild ; i++ )
         tree_order(tree->children[i]);
}


/* find which interval the number value lies in, the index found	*/
/* indicates the child tree which should contain such a point 		*/
int tree_whichinterval(t, val)
TREE *t; float val;
{
int i;
   for ( i = 0 ; i < t->nchild ; i++ )
      if ( val >= t->children[i]->vmin && val <= t->children[i]->vmax )
         return i;

   fatalerr("tree_whichinterval", "value", "not in any interval");
}


void treefree(tree)
TREE *tree;
{
int i;

   if (istreechildless(tree))
   {
      if (tree->pointsdata != NULL) free(tree->pointsdata);
      if (tree->points != NULL) free(tree->points);
   }
   else
   {
      for ( i = 0 ; i < tree->nchild ; i++ )
         treefree(tree->children[i]);
      free(tree->children);
   }
   free(tree);
}