Files
SaraP 739088af9f Vroni 7.8 :
- aggiornamento versione.
2025-01-29 16:24:30 +01:00

528 lines
18 KiB
C++

/*****************************************************************************/
/* */
/* 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 <stdio.h>
#include <assert.h>
#include <math.h>
#include <float.h>
/* */
/* 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