/*****************************************************************************/ /* */ /* Copyright (C) 2003-2025 M. Held */ /* */ /* This code is not in the public domain. All rights reserved! Please make */ /* sure to read the full copyright statement contained in "README.txt" or in */ /* the "main" file of this code, such as "main.cc". */ /* */ /*****************************************************************************/ /* */ /* Written by: Martin Held */ /* */ /* E-Mail: held@cs.sbg.ac.at */ /* Fax Mail: (+4 662) 8044-611 */ /* Voice Mail: (+4 662) 8044-6304 */ /* Snail Mail: Martin Held */ /* FB Informatik */ /* Universitaet Salzburg */ /* A-5020 Salzburg, Austria */ /* */ /*****************************************************************************/ /* */ /* this file maintains a quadtree-like data structure which is used for */ /* speeding up the search for nearest neighbors during the construction of */ /* VDs of highly non-uniformly distributed points. the tree is generated */ /* "on the fly" with no additional preprocessing. when a new point is to be */ /* inserted into the tree, the cell/node that contains it is located. if its */ /* squared distance to the nearest neighbor is less than ZERO2 = ZERO^2 */ /* then the point is not inserted. otherwise, the cell that contains it is */ /* split into four sub-cells of (typically) non-equal size by using this */ /* point as the splitter. cells that cannot immediately searched during the */ /* tree traversal are maintained in a heap that is ordered according to the */ /* distance of the cell from the point to be inserted. */ /* */ /* please note that a random insertion of the points is assumed. (this is */ /* the default compile-time option for VRONI!) inserting the points in some */ /* prearranged form, such as sorted by x-coordinates, is likely to cause the */ /* depth of the tree to grow substantially, which may once again result in */ /* an overall quadratic time complexity. in order to enforce fairly balanced */ /* splits during the first few levels, good splitters are determined by an */ /* explicit search among all points that are to be inserted at once. (recall */ /* that this tree is only constructed once the handling of the first few */ /* percent of the points by means of a grid indicated a highly non-uniform */ /* distribution of the points.) */ /* */ /* */ /* get standard libraries */ /* */ #include #include #include #include /* */ /* get my header files */ /* */ #include "fpkernel.h" #include "vronivector.h" #include "vroni_object.h" #include "defs.h" #include "numerics.h" #define Initial_Depth 2 #define GetSplitter(node_idx) \ (\ assert(InTree(node_idx)), \ tree[node_idx].splitter \ ) #define SetChild_1(node_idx, ind) \ {\ assert(InTree(node_idx)); \ tree[node_idx].ch1 = ind; \ } #define GetChild_1(node_idx) \ (\ assert(InTree(node_idx)), \ tree[node_idx].ch1 \ ) #define SetChild_2(node_idx, ind) \ {\ assert(InTree(node_idx)); \ tree[node_idx].ch2 = ind; \ } #define GetChild_2(node_idx) \ (\ assert(InTree(node_idx)), \ tree[node_idx].ch2 \ ) #define SetChild_3(node_idx, ind) \ {\ assert(InTree(node_idx)); \ tree[node_idx].ch3 = ind; \ } #define GetChild_3(node_idx) \ (\ assert(InTree(node_idx)), \ tree[node_idx].ch3 \ ) #define SetChild_4(node_idx, ind) \ {\ assert(InTree(node_idx)); \ tree[node_idx].ch4 = ind; \ } #define GetChild_4(node_idx) \ (\ assert(InTree(node_idx)), \ tree[node_idx].ch4 \ ) #define UpdateDistance(split, diff, min_ind, min_dist) \ {\ if (diff < min_dist) { \ min_dist = diff; \ min_ind = split; \ } \ } #define MakeTreeNode(node_idx, ind, addr) \ { \ if (num_tree >= max_num_tree) { \ max_num_tree += HEAP_BLOCK_SIZE; \ tree = (tree_node*) ReallocateArray(tree, max_num_tree, \ sizeof(tree_node), "tree:tree"); \ } \ node_idx = num_tree; \ addr = &(tree[node_idx]); \ addr->splitter = ind; \ addr->ch1 = addr->ch2 = addr->ch3 = addr->ch4 = NIL; \ \ ++num_tree; \ } #define RememberCell(node_idx, dist, distance) \ {\ if (node_idx != NIL) {\ if ((dist) <= distance) InsertIntoHeap(dist, node_idx); \ }\ } #define RecoverNextCell(node_idx, diff) \ (\ DeleteFromHeap(&diff, &node_idx) \ ) void vroniObject::RecursiveInitTree(double x_min, double y_min, double x_max, double y_max, int *first_N_points, int *m, int level) { double q_x, q_y, p_x, p_y, min_dist, min_diff, diff_x, diff_y, diff; int i_diff, n, ind; tree_node *addr; /* */ /* looking for a good splitter: we determine the point that is closest to */ /* the center of this cell. */ /* */ /* */ q_x = (x_min + x_max) / 2.0; q_y = (y_min + y_max) / 2.0; min_diff = DBL_MAX; i_diff = NIL; for (n = 0; n < *m; ++n) { p_x = GetPntCoords(first_N_points[n]).x; p_y = GetPntCoords(first_N_points[n]).y; if ((x_min <= p_x) && (p_x <= x_max) && (y_min <= p_y) && (p_y <= y_max)) { diff_x = q_x - p_x; diff_y = q_y - p_y; diff = diff_x * diff_x + diff_y * diff_y; if (diff < min_diff) { i_diff = n; min_diff = diff; } } } if (i_diff == NIL) return; assert(InPntsList(first_N_points[i_diff])); if (root_of_tree == NIL) { /* */ /* use this splitter as the first point to be inserted into the tree. */ /* */ ind = first_N_points[i_diff]; MakeTreeNode(root_of_tree, ind, addr); assert(root_of_tree == 0); } else { (void) InsertPntIntoTree(first_N_points[i_diff], &min_dist, true); } p_x = GetPntCoords(first_N_points[i_diff]).x; p_y = GetPntCoords(first_N_points[i_diff]).y; --(*m); first_N_points[i_diff] = first_N_points[*m]; if (level < Initial_Depth) { /* */ /* look for good splitters in the four sub-cells. */ /* */ RecursiveInitTree(x_min, y_min, p_x, p_y, first_N_points, m, level + 1); RecursiveInitTree(x_min, p_y, p_x, y_max, first_N_points, m, level + 1); RecursiveInitTree(p_x, y_min, x_max, p_y, first_N_points, m, level + 1); RecursiveInitTree(p_x, p_y, x_max, y_max, first_N_points, m, level + 1); } return; } int vroniObject::InitTree(int number, int *first_N_points, int N_Initial, double *min_dist) { double eps; int m, n, i_min; if (number>= max_num_tree) { max_num_tree = number; tree = (tree_node*) ReallocateArray(tree, max_num_tree, sizeof(tree_node), "tree:tree"); } max_num_tree = 0; num_tree = 0; root_of_tree = NIL; eps = ZERO_MAX; /* */ /* determine the location of the tree */ /* */ tree_delta_x = (bb_max.x - bb_min.x) * 0.01; tree_delta_y = (bb_max.y - bb_min.y) * 0.01; if (le(tree_delta_x, eps)) tree_delta_x = eps; if (le(tree_delta_y, eps)) tree_delta_y = eps; tree_min_x = bb_min.x - tree_delta_x; tree_max_x = bb_max.x + tree_delta_x; tree_min_y = bb_min.y - tree_delta_y; tree_max_y = bb_max.y + tree_delta_y; tree_delta_x = tree_max_x - tree_min_x; tree_delta_y = tree_max_y - tree_min_y; assert(gt(tree_delta_x, ZERO)); assert(gt(tree_delta_y, ZERO)); /* */ /* insert all points that have been processed so far into the tree. we */ /* start finding a few good splitters for some top levels of the tree. */ /* */ assert(N_Initial >= 1); m = N_Initial; RecursiveInitTree(tree_min_x, tree_min_y, tree_max_x, tree_max_y, first_N_points, &m, 0); assert(root_of_tree != NIL); /* */ /* insert all remaining points */ /* */ for (n = 0; n < m; ++n) { i_min = InsertPntIntoTree(first_N_points[n], min_dist, true); } assert(InPntsList(first_N_points[N_Initial])); i_min = InsertPntIntoTree(first_N_points[N_Initial], min_dist, false); return i_min; } vr_bool vroniObject::InTree(int idx) { return ((idx >= 0) && (idx < num_tree)); } void vroniObject::FreeTree(void) { FreeHeap(); FreeMemory((void**) &tree, "tree:tree"); max_num_tree = 0; num_tree = 0; root_of_tree = NIL; return; } /* */ /* insert pnts[ind] into the tree, and determine its nearest neighbor */ /* (unless the flag no_distance is true). */ /* */ int vroniObject::InsertPntIntoTree(int ind, double *dist_min, vr_bool no_distance) { coord p, q; int root, min_ind, split; int child, next; double min_dist, diff, diff_x, diff_y, diff2_x, diff2_y; tree_node *addr; vr_bool remember_insert; int root_insert, child_insert; p = GetPntCoords(ind); min_dist = DBL_MAX; min_ind = NIL; next = root_of_tree; remember_insert = true; root_insert = NIL; child_insert = NIL; assert(p.x > tree_min_x); assert(p.y > tree_min_y); assert(p.x < tree_max_x); assert(p.y < tree_max_y); assert(next == 0); ResetHeap(); do { /* */ /* search entire tree for a nearest neighbor */ /* */ do { /* */ /* search one subtree for a nearest neighbor */ /* */ root = next; split = GetSplitter(root); q = GetPntCoords(split); diff_x = p.x - q.x; diff_y = p.y - q.y; diff2_x = diff_x * diff_x; diff2_y = diff_y * diff_y; diff = diff2_x + diff2_y; UpdateDistance(split, diff, min_ind, min_dist); if (diff_x < 0.0) { /* */ /* point is left of the vertical line through this splitter */ /* */ if (diff_y < 0.0) { /* */ /* point is below the horizontal line through this splitter */ /* */ next = GetChild_2(root); RememberCell(next, diff2_y, min_dist); next = GetChild_3(root); RememberCell(next, diff2_x, min_dist); next = GetChild_1(root); child = 1; } else { /* */ /* point is above the horizontal line through this splitter */ /* */ next = GetChild_1(root); RememberCell(next, diff2_y, min_dist); next = GetChild_4(root); RememberCell(next, diff2_x, min_dist); next = GetChild_2(root); child = 2; } } else { /* */ /* point is right of the vertical line through this splitter */ /* */ if (diff_y < 0.0) { /* */ /* point is below the horizontal line through this splitter */ /* */ next = GetChild_4(root); RememberCell(next, diff2_y, min_dist); next = GetChild_1(root); RememberCell(next, diff2_x, min_dist); next = GetChild_3(root); child = 3; } else { /* */ /* point is above the horizontal line through this splitter */ /* */ next = GetChild_3(root); RememberCell(next, diff2_y, min_dist); next = GetChild_2(root); RememberCell(next, diff2_x, min_dist); next = GetChild_4(root); child = 4; } } } while (next != NIL); if (remember_insert) { remember_insert = false; root_insert = root; child_insert = child; if (no_distance) ResetHeap(); } while ((next == NIL) && RecoverNextCell(root, diff)) { if (diff <= min_dist) next = root; } } while (next != NIL); if (gt(min_dist, ZERO2)) { assert(root_insert != NIL); assert(child_insert != NIL); MakeTreeNode(next, ind, addr); switch(child_insert) { case 1: SetChild_1(root_insert, next); break; case 2: SetChild_2(root_insert, next); break; case 3: SetChild_3(root_insert, next); break; case 4: SetChild_4(root_insert, next); break; default: assert((1 <= child_insert) && (child_insert <= 4)); } } assert(min_ind != NIL); *dist_min = min_dist; return min_ind; } #ifdef GRAPHICS void vroniObject::RecursiveDrawTree(double x_min, double y_min, double_arg x_max, double_arg y_max, int root) { coord p1, p2, q; int split, next; while (root != NIL) { split = GetSplitter(root); q = GetPntCoords(split); p1.x = x_min; p1.y = q.y; p2.x = x_max; p2.y = q.y; AddEdgeToBuffer(p1.x, p1.y, p2.x, p2.y, TreeColor); p1.x = q.x; p1.y = y_min; p2.x = q.x; p2.y = y_max; AddEdgeToBuffer(p1.x, p1.y, p2.x, p2.y, TreeColor); next = GetChild_1(root); if (next != NIL) { RecursiveDrawTree(x_min, y_min, q.x, q.y, next); } next = GetChild_2(root); if (next != NIL) { RecursiveDrawTree(x_min, q.y, q.x, y_max, next); } next = GetChild_3(root); if (next != NIL) { RecursiveDrawTree(q.x, y_min, x_max, q.y, next); } root = GetChild_4(root); if (root != NIL) { x_min = q.x; y_min = q.y; } } return; } void vroniObject::DrawTree(void) { coord p1, p2, p3, p4; p1.x = tree_min_x; p1.y = tree_min_y; p2.x = tree_max_x; p2.y = tree_min_y; p3.x = tree_max_x; p3.y = tree_max_y; p4.x = tree_min_x; p4.y = tree_max_y; AddEdgeToBuffer(p1.x, p1.y, p2.x, p2.y, TreeColor); AddEdgeToBuffer(p2.x, p2.y, p3.x, p3.y, TreeColor); AddEdgeToBuffer(p3.x, p3.y, p4.x, p4.y, TreeColor); AddEdgeToBuffer(p4.x, p4.y, p1.x, p1.y, TreeColor); RecursiveDrawTree(tree_min_x, tree_min_y, tree_max_x, tree_max_y, root_of_tree); return; } #endif