classify/kdtree.cpp

Go to the documentation of this file.

/* =================
 Include Files and Type Defines
 ==================== */
#include "kdtree.h"
#include "const.h"
#include "emalloc.h"
#include "freelist.h"
#include <stdio.h>
#include <math.h>
#include <setjmp.h>

#define Magnitude(X)    ((X) < 0 ? -(X) : (X))
#define MIN(A,B)    ((A) < (B) ? (A) : (B))
#define NodeFound(N,K,D)  (( (N)->Key == (K) ) && ( (N)->Data == (D) ))

/* =================
        Global Data Definitions and Declarations
 ==================== */
#define MINSEARCH -MAX_FLOAT32
#define MAXSEARCH MAX_FLOAT32

static int NumberOfNeighbors;
static INT16 N;

static FLOAT32 *QueryPoint;
static int MaxNeighbors;
static FLOAT32 Radius;
static int Furthest;
static char **Neighbor;
static FLOAT32 *Distance;

static int MaxDimension = 0;
static FLOAT32 *SBMin;
static FLOAT32 *SBMax;
static FLOAT32 *LBMin;
static FLOAT32 *LBMax;

static PARAM_DESC *KeyDesc;

static jmp_buf QuickExit;

static void_proc WalkAction;

/* =================
              Public Code
 ==================== */
/* =============================== */
KDTREE *
MakeKDTree (INT16 KeySize, PARAM_DESC KeyDesc[]) {
  int i;
  void *NewMemory;
  KDTREE *KDTree;

  if (KeySize > MaxDimension) {
    NewMemory = Emalloc (KeySize * 4 * sizeof (FLOAT32));
    if (MaxDimension > 0) {
      memfree ((char *) SBMin);
      memfree ((char *) SBMax);
      memfree ((char *) LBMin);
      memfree ((char *) LBMax);
    }
    SBMin = (FLOAT32 *) NewMemory;
    SBMax = SBMin + KeySize;
    LBMin = SBMax + KeySize;
    LBMax = LBMin + KeySize;
  }

  KDTree =
    (KDTREE *) Emalloc (sizeof (KDTREE) +
    (KeySize - 1) * sizeof (PARAM_DESC));
  for (i = 0; i < KeySize; i++) {
    KDTree->KeyDesc[i].NonEssential = KeyDesc[i].NonEssential;
    KDTree->KeyDesc[i].Circular = KeyDesc[i].Circular;
    if (KeyDesc[i].Circular) {
      KDTree->KeyDesc[i].Min = KeyDesc[i].Min;
      KDTree->KeyDesc[i].Max = KeyDesc[i].Max;
      KDTree->KeyDesc[i].Range = KeyDesc[i].Max - KeyDesc[i].Min;
      KDTree->KeyDesc[i].HalfRange = KDTree->KeyDesc[i].Range / 2;
      KDTree->KeyDesc[i].MidRange = (KeyDesc[i].Max + KeyDesc[i].Min) / 2;
    }
    else {
      KDTree->KeyDesc[i].Min = MINSEARCH;
      KDTree->KeyDesc[i].Max = MAXSEARCH;
    }
  }
  KDTree->KeySize = KeySize;
  KDTree->Root.Left = NULL;
  KDTree->Root.Right = NULL;
  return (KDTree);
}                                /* MakeKDTree */


/* =============================== */
void KDStore(KDTREE *Tree, FLOAT32 *Key, void *Data) { 
  int Level;
  KDNODE *Node;
  KDNODE **PtrToNode;

  N = Tree->KeySize;
  KeyDesc = &(Tree->KeyDesc[0]);
  PtrToNode = &(Tree->Root.Left);
  Node = *PtrToNode;
  Level = 0;
  while (Node != NULL) {
    if (Key[Level] < Node->BranchPoint) {
      PtrToNode = &(Node->Left);
      if (Key[Level] > Node->LeftBranch)
        Node->LeftBranch = Key[Level];
    }
    else {
      PtrToNode = &(Node->Right);
      if (Key[Level] < Node->RightBranch)
        Node->RightBranch = Key[Level];
    }
    Level++;
    if (Level >= N)
      Level = 0;
    Node = *PtrToNode;
  }

  *PtrToNode = MakeKDNode (Key, (char *) Data, Level);
}                                /* KDStore */


/* =============================== */
void
KDDelete (KDTREE * Tree, FLOAT32 Key[], void *Data) {
  int Level;
  KDNODE *Current;
  KDNODE *Father;
  KDNODE *Replacement;
  KDNODE *FatherReplacement;

  /* initialize search at root of tree */
  N = Tree->KeySize;
  KeyDesc = &(Tree->KeyDesc[0]);
  Father = &(Tree->Root);
  Current = Father->Left;
  Level = 0;

  /* search tree for node to be deleted */
  while ((Current != NULL) && (!NodeFound (Current, Key, Data))) {
    Father = Current;
    if (Key[Level] < Current->BranchPoint)
      Current = Current->Left;
    else
      Current = Current->Right;

    Level++;
    if (Level >= N)
      Level = 0;
  }

  if (Current != NULL) {         /* if node to be deleted was found */
    Replacement = Current;
    FatherReplacement = Father;

    /* search for replacement node (a leaf under node to be deleted */
    while (TRUE) {
      if (Replacement->Left != NULL) {
        FatherReplacement = Replacement;
        Replacement = Replacement->Left;
      }
      else if (Replacement->Right != NULL) {
        FatherReplacement = Replacement;
        Replacement = Replacement->Right;
      }
      else
        break;

      Level++;
      if (Level >= N)
        Level = 0;
    }

    /* compute level of replacement node's father */
    Level--;
    if (Level < 0)
      Level = N - 1;

    /* disconnect replacement node from it's father */
    if (FatherReplacement->Left == Replacement) {
      FatherReplacement->Left = NULL;
      FatherReplacement->LeftBranch = KeyDesc[Level].Min;
    }
    else {
      FatherReplacement->Right = NULL;
      FatherReplacement->RightBranch = KeyDesc[Level].Max;
    }

    /* replace deleted node with replacement (unless they are the same) */
    if (Replacement != Current) {
      Replacement->BranchPoint = Current->BranchPoint;
      Replacement->LeftBranch = Current->LeftBranch;
      Replacement->RightBranch = Current->RightBranch;
      Replacement->Left = Current->Left;
      Replacement->Right = Current->Right;

      if (Father->Left == Current)
        Father->Left = Replacement;
      else
        Father->Right = Replacement;
    }
    FreeKDNode(Current); 
  }
}                                /* KDDelete */


/* =============================== */
int
KDNearestNeighborSearch (KDTREE * Tree,
FLOAT32 Query[],
int QuerySize,
FLOAT32 MaxDistance,
void *NBuffer, FLOAT32 DBuffer[]) {
  int i;

  NumberOfNeighbors = 0;
  N = Tree->KeySize;
  KeyDesc = &(Tree->KeyDesc[0]);
  QueryPoint = Query;
  MaxNeighbors = QuerySize;
  Radius = MaxDistance;
  Furthest = 0;
  Neighbor = (char **) NBuffer;
  Distance = DBuffer;

  for (i = 0; i < N; i++) {
    SBMin[i] = KeyDesc[i].Min;
    SBMax[i] = KeyDesc[i].Max;
    LBMin[i] = KeyDesc[i].Min;
    LBMax[i] = KeyDesc[i].Max;
  }

  if (Tree->Root.Left != NULL) {
    if (setjmp (QuickExit) == 0)
      Search (0, Tree->Root.Left);
  }
  return (NumberOfNeighbors);
}                                /* KDNearestNeighborSearch */


/* =============================== */
void KDWalk(KDTREE *Tree, void_proc Action) { 
  WalkAction = Action;
  if (Tree->Root.Left != NULL)
    Walk (Tree->Root.Left, 0);
}                                /* KDWalk */


/* =============================== */
void FreeKDTree(KDTREE *Tree) { 
  FreeSubTree (Tree->Root.Left);
  memfree(Tree); 
}                                /* FreeKDTree */


/* =================
 Private Code
 ==================== */
int
Equal (FLOAT32 Key1[], FLOAT32 Key2[]) {
  int i;

  for (i = N; i > 0; i--, Key1++, Key2++)
    if (*Key1 != *Key2)
      return (FALSE);
  return (TRUE);
}                                /* Equal */


/* =============================== */
KDNODE *
MakeKDNode (FLOAT32 Key[], char *Data, int Index) {
  KDNODE *NewNode;

  NewNode = (KDNODE *) Emalloc (sizeof (KDNODE));

  NewNode->Key = Key;
  NewNode->Data = Data;
  NewNode->BranchPoint = Key[Index];
  NewNode->LeftBranch = KeyDesc[Index].Min;
  NewNode->RightBranch = KeyDesc[Index].Max;
  NewNode->Left = NULL;
  NewNode->Right = NULL;

  return (NewNode);
}                                /* MakeKDNode */


/* =============================== */
void FreeKDNode(KDNODE *Node) { 
  memfree ((char *) Node);
}                                /* FreeKDNode */


/* =============================== */
void Search(int Level, KDNODE *SubTree) { 
  FLOAT32 d;
  FLOAT32 OldSBoxEdge;
  FLOAT32 OldLBoxEdge;

  if (Level >= N)
    Level = 0;

  d = ComputeDistance (N, KeyDesc, QueryPoint, SubTree->Key);
  if (d < Radius) {
    if (NumberOfNeighbors < MaxNeighbors) {
      Neighbor[NumberOfNeighbors] = SubTree->Data;
      Distance[NumberOfNeighbors] = d;
      NumberOfNeighbors++;
      if (NumberOfNeighbors == MaxNeighbors)
        FindMaxDistance(); 
    }
    else {
      Neighbor[Furthest] = SubTree->Data;
      Distance[Furthest] = d;
      FindMaxDistance(); 
    }
  }
  if (QueryPoint[Level] < SubTree->BranchPoint) {
    OldSBoxEdge = SBMax[Level];
    SBMax[Level] = SubTree->LeftBranch;
    OldLBoxEdge = LBMax[Level];
    LBMax[Level] = SubTree->RightBranch;
    if (SubTree->Left != NULL)
      Search (Level + 1, SubTree->Left);
    SBMax[Level] = OldSBoxEdge;
    LBMax[Level] = OldLBoxEdge;
    OldSBoxEdge = SBMin[Level];
    SBMin[Level] = SubTree->RightBranch;
    OldLBoxEdge = LBMin[Level];
    LBMin[Level] = SubTree->LeftBranch;
    if ((SubTree->Right != NULL) && QueryIntersectsSearch ())
      Search (Level + 1, SubTree->Right);
    SBMin[Level] = OldSBoxEdge;
    LBMin[Level] = OldLBoxEdge;
  }
  else {
    OldSBoxEdge = SBMin[Level];
    SBMin[Level] = SubTree->RightBranch;
    OldLBoxEdge = LBMin[Level];
    LBMin[Level] = SubTree->LeftBranch;
    if (SubTree->Right != NULL)
      Search (Level + 1, SubTree->Right);
    SBMin[Level] = OldSBoxEdge;
    LBMin[Level] = OldLBoxEdge;
    OldSBoxEdge = SBMax[Level];
    SBMax[Level] = SubTree->LeftBranch;
    OldLBoxEdge = LBMax[Level];
    LBMax[Level] = SubTree->RightBranch;
    if ((SubTree->Left != NULL) && QueryIntersectsSearch ())
      Search (Level + 1, SubTree->Left);
    SBMax[Level] = OldSBoxEdge;
    LBMax[Level] = OldLBoxEdge;
  }
  if (QueryInSearch ())
    longjmp (QuickExit, 1);
}                                /* Search */


/* =============================== */
FLOAT32
ComputeDistance (register int N,
register PARAM_DESC Dim[],
register FLOAT32 p1[], register FLOAT32 p2[]) {
  register FLOAT32 TotalDistance;
  register FLOAT32 DimensionDistance;
  FLOAT32 WrapDistance;

  TotalDistance = 0;
  for (; N > 0; N--, p1++, p2++, Dim++) {
    if (Dim->NonEssential)
      continue;

    DimensionDistance = *p1 - *p2;

    /* if this dimension is circular - check wraparound distance */
    if (Dim->Circular) {
      DimensionDistance = Magnitude (DimensionDistance);
      WrapDistance = Dim->Max - Dim->Min - DimensionDistance;
      DimensionDistance = MIN (DimensionDistance, WrapDistance);
    }

    TotalDistance += DimensionDistance * DimensionDistance;
  }
  return ((FLOAT32) sqrt ((FLOAT64) TotalDistance));
}                                /* ComputeDistance */


/* =============================== */
void FindMaxDistance() { 
  int i;

  Radius = Distance[Furthest];
  for (i = 0; i < MaxNeighbors; i++) {
    if (Distance[i] > Radius) {
      Radius = Distance[i];
      Furthest = i;
    }
  }
}                                /* FindMaxDistance */


/* =============================== */
int QueryIntersectsSearch() { 
  register int i;
  register FLOAT32 *Query;
  register FLOAT32 *Lower;
  register FLOAT32 *Upper;
  register FLOAT64 TotalDistance;
  register FLOAT32 DimensionDistance;
  register FLOAT64 RadiusSquared;
  register PARAM_DESC *Dim;
  register FLOAT32 WrapDistance;

  RadiusSquared = Radius * Radius;
  Query = QueryPoint;
  Lower = SBMin;
  Upper = SBMax;
  TotalDistance = 0.0;
  Dim = KeyDesc;
  for (i = N; i > 0; i--, Dim++, Query++, Lower++, Upper++) {
    if (Dim->NonEssential)
      continue;

    if (*Query < *Lower)
      DimensionDistance = *Lower - *Query;
    else if (*Query > *Upper)
      DimensionDistance = *Query - *Upper;
    else
      DimensionDistance = 0;

    if (Dim->Circular) {
      if (*Query < *Lower)
        WrapDistance = *Query + Dim->Max - Dim->Min - *Upper;
      else if (*Query > *Upper)
        WrapDistance = *Lower - (*Query - (Dim->Max - Dim->Min));
      else
        WrapDistance = MAX_FLOAT32;

      DimensionDistance = MIN (DimensionDistance, WrapDistance);
    }

    TotalDistance += DimensionDistance * DimensionDistance;
    if (TotalDistance >= RadiusSquared)
      return (FALSE);
  }
  return (TRUE);
}                                /* QueryIntersectsSearch */


/* =============================== */
int QueryInSearch() { 
  register int i;
  register FLOAT32 *Query;
  register FLOAT32 *Lower;
  register FLOAT32 *Upper;
  register PARAM_DESC *Dim;

  Query = QueryPoint;
  Lower = LBMin;
  Upper = LBMax;
  Dim = KeyDesc;

  for (i = N - 1; i >= 0; i--, Dim++, Query++, Lower++, Upper++) {
    if (Dim->NonEssential)
      continue;

    if ((*Query < *Lower + Radius) || (*Query > *Upper - Radius))
      return (FALSE);
  }
  return (TRUE);
}                                /* QueryInSearch */


/* =============================== */
void Walk(KDNODE *SubTree, INT32 Level) { 
  if ((SubTree->Left == NULL) && (SubTree->Right == NULL))
    (*WalkAction) (SubTree->Data, leaf, Level);
  else {
    (*WalkAction) (SubTree->Data, preorder, Level);
    if (SubTree->Left != NULL)
      Walk (SubTree->Left, Level + 1);
    (*WalkAction) (SubTree->Data, postorder, Level);
    if (SubTree->Right != NULL)
      Walk (SubTree->Right, Level + 1);
    (*WalkAction) (SubTree->Data, endorder, Level);
  }
}                                /* Walk */


/* =============================== */
void FreeSubTree(KDNODE *SubTree) { 
  if (SubTree != NULL) {
    FreeSubTree (SubTree->Left);
    FreeSubTree (SubTree->Right);
    memfree(SubTree); 
  }
}                                /* FreeSubTree */

---