kdtree.cpp

一ＯＣＲ的相关资料。．希望对研究ＯＣＲ的朋友有所帮助．

 **	Parameters:
 **		Key1,Key2	search keys to be compared for equality
 **	Globals:
 **		N		number of parameters per key
 **	Operation:
 **		This routine returns TRUE if Key1 = Key2.
 **	Return:
 **		TRUE if Key1 = Key2, else FALSE.
 **	Exceptions:
 **		None
 **	History:
 **		3/11/89, DSJ, Created.
 */
  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) {
/*
 **	Parameters:
 **		Key	Access key for new node in KD tree
 **		Data	ptr to data to be stored in new node
 **		Index	index of Key to branch on
 **	Globals:
 **		KeyDesc	descriptions of key dimensions
 **	Operation:
 **		This routine allocates memory for a new K-D tree node
 **		and places the specified Key and Data into it.  The
 **		left and right subtree pointers for the node are
 **		initialized to empty subtrees.
 **	Return:
 **		pointer to new K-D tree node
 **	Exceptions:
 **		None
 **	History:
 **		3/11/89, DSJ, Created.
 */
  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) { 
/*
 **	Parameters:
 **		Node	ptr to node data structure to be freed
 **	Globals:
 **		None
 **	Operation:
 **		This routine frees up the memory allocated to Node.
 **	Return:
 **		None
 **	Exceptions:
 **		None
 **	History:
 **		3/13/89, DSJ, Created.
 */
  memfree ((char *) Node);
}                                /* FreeKDNode */


/*---------------------------------------------------------------------------*/
void Search(int Level, KDNODE *SubTree) { 
/*
 **	Parameters:
 **		Level		level in tree of sub-tree to be searched
 **		SubTree		sub-tree to be searched
 **	Globals:
 **		NumberOfNeighbors	# of neighbors found so far
 **		N			# of features in each key
 **		KeyDesc			description of key dimensions
 **		QueryPoint		point in D-space to find neighbors of
 **		MaxNeighbors		maximum # of neighbors to find
 **		Radius			current distance of furthest neighbor
 **		Furthest		index of furthest neighbor
 **		Neighbor		buffer of current neighbors
 **		Distance		buffer of neighbor distances
 **		SBMin			lower extent of small search region
 **		SBMax			upper extent of small search region
 **		LBMin			lower extent of large search region
 **		LBMax			upper extent of large search region
 **		QuickExit		quick exit from recursive search
 **	Operation:
 **		This routine searches SubTree for those entries which are
 **		possibly among the MaxNeighbors nearest neighbors of the
 **		QueryPoint and places their data in the Neighbor buffer and
 **		their distances from QueryPoint in the Distance buffer.
 **	Return: none
 **	Exceptions: none
 **	History:
 **		3/11/89, DSJ, Created.
 **		7/13/89, DSJ, Save node contents, not node, in neighbor buffer
 */
  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[]) {
/*
 **	Parameters:
 **		N		number of dimensions in K-D space
 **		Dim		descriptions of each dimension
 **		p1,p2		two different points in K-D space
 **	Globals:
 **		None
 **	Operation:
 **		This routine computes the euclidian distance
 **		between p1 and p2 in K-D space (an N dimensional space).
 **	Return:
 **		Distance between p1 and p2.
 **	Exceptions:
 **		None
 **	History:
 **		3/11/89, DSJ, Created.
 */
  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() { 
/*
 **	Parameters:
 **		None
 **	Globals:
 **		MaxNeighbors		maximum # of neighbors to find
 **		Radius			current distance of furthest neighbor
 **		Furthest		index of furthest neighbor
 **		Distance		buffer of neighbor distances
 **	Operation:
 **		This routine searches the Distance buffer for the maximum
 **		distance, places this distance in Radius, and places the
 **		index of this distance in Furthest.
 **	Return:
 **		None
 **	Exceptions:
 **		None
 **	History:
 **		3/11/89, DSJ, Created.
 */
  int i;

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


/*---------------------------------------------------------------------------*/
int QueryIntersectsSearch() { 
/*
 **	Parameters:
 **		None
 **	Globals:
 **		N			# of features in each key
 **		KeyDesc			descriptions of each dimension
 **		QueryPoint		point in D-space to find neighbors of
 **		Radius			current distance of furthest neighbor
 **		SBMin			lower extent of small search region
 **		SBMax			upper extent of small search region
 **	Operation:
 **		This routine returns TRUE if the query region intersects
 **		the current smallest search region.  The query region is
 **		the circle of radius Radius centered at QueryPoint.
 **		The smallest search region is the box (in N dimensions)
 **		whose edges in each dimension are specified by SBMin and SBMax.
 **		In the case of circular dimensions, we must also check the
 **		point which is one wrap-distance away from the query to
 **		see if it would intersect the search region.
 **	Return:
 **		TRUE if query region intersects search region, else FALSE
 **	Exceptions:
 **		None
 **	History:
 **		3/11/89, DSJ, Created.
 */
  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 this dimension is circular - check wraparound distance */
    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() { 
/*
 **	Parameters:
 **		None
 **	Globals:
 **		N			# of features in each key
 **		KeyDesc			descriptions of each dimension
 **		QueryPoint		point in D-space to find neighbors of
 **		Radius			current distance of furthest neighbor
 **		LBMin			lower extent of large search region
 **		LBMax			upper extent of large search region
 **	Operation:
 **		This routine returns TRUE if the current query region is
 **		totally contained in the current largest search region.
 **		The query region is the circle of
 **		radius Radius centered at QueryPoint.  The search region is
 **		the box (in N dimensions) whose edges in each
 **		dimension are specified by LBMin and LBMax.
 **	Return:
 **		TRUE if query region is inside search region, else FALSE
 **	Exceptions:
 **		None
 **	History:
 **		3/11/89, DSJ, Created.
 */
  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) { 
/*
 **	Parameters:
 **		SubTree		ptr to root of subtree to be walked
 **		Level		current level in the tree for this node
 **	Globals:
 **		WalkAction	action to be performed at every node
 **	Operation:
 **		This routine walks thru the specified SubTree and invokes
 **		WalkAction at each node.  WalkAction is invoked with three
 **		arguments as follows:
 **			WalkAction( NodeData, Order, Level )
 **		Data is the data contents of the node being visited,
 **		Order is either preorder,
 **		postorder, endorder, or leaf depending on whether this is
 **		the 1st, 2nd, or 3rd time a node has been visited, or
 **		whether the node is a leaf.  Level is the level of the node in
 **		the tree with the root being level 0.
 **	Return: none
 **	Exceptions: none
 **	History:
 **		3/13/89, DSJ, Created.
 **		7/13/89, DSJ, Pass node contents, not node, to WalkAction().
 */
  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) { 
/*
 **	Parameters:
 **		SubTree		ptr to root node of sub-tree to be freed
 **	Globals: none
 **	Operation:
 **		This routine recursively frees the memory allocated to
 **		to the specified subtree.
 **	Return: none
 **	Exceptions: none
 **	History: 7/13/89, DSJ, Created.
 */
  if (SubTree != NULL) {
    FreeSubTree (SubTree->Left);
    FreeSubTree (SubTree->Right);
    memfree(SubTree); 
  }
}                                /* FreeSubTree */