Egalware/vroni
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
739088af9faac764e8d7a0a811cda5057b1cbc6b
vroni/grid.cc
T
SaraP 739088af9f Vroni 7.8 : 
- aggiornamento versione.
2025-01-29 16:24:30 +01:00

612 lines
18 KiB
C++
Raw Blame History

/*****************************************************************************/
/* */
/* 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 simple regular grid which is used for speeding up */
/* the search for nearest neighbors during the construction of point VDs. */
/* initially, a grid suitable for num_pnts * N_Percent many points is */
/* maintained. based on the average number of points per cell, either a grid */
/* or a quadtree-like structure is then constructed on the full set of */
/* num_pnts many points. the layout of both grids is adapted to the aspect */
/* ratio of the bounding box of the points. the number of cells of the */
/* second grid depends on the average number of points per cell in the first */
/* grid, and varies from num_pnts many cells to num_pnts * MAX_FACTOR */
/* many cells. for small datasets with 8 * num_pnts * N_Percent < N_Trivial */
/* only one grid with num_pnts * STD_FACTOR many cells is constructed. */
/* */
/* whether or not a tree or a grid is used depends on how non-uniformly the */
/* points are distributed. if tons of points would end up in only very few */
/* cells then the grid would be useless, and the nearest-neighbor search on */
/* K points would exhibit O(K^2) complexity. as stated above, a hopefully */
/* good guess of the distribution of the points is attempted by computing */
/* the average number of points per cell within the first grid, for a random */
/* sample of num_pnts * N_Percent points. let me clearly state that i do */
/* realize that this approach has no theoretical guarantee to avoid a */
/* quadratic complexity. (yes, it is possible to contrive datasets that fool */
/* my approach: put half of the points extremely close to the lower-left */
/* corner of the unit square, and the other half of the points extremely */
/* close to the upper-right corner of the unit square. then even the */
/* handling of the first grid for the first few percent of the points will */
/* exhibit a quadratic time complexity...) and, yes, i could complicate the */
/* code even further to catch such a situation. however, it seems that even */
/* highly non-uniform data is handled nicely by the grid, and extreme cases */
/* are handled by the tree. that is, i do not really feel a need to improve */
/* this approach... */
/* */
/* */
/* 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 "numerics.h"
#define PAGE_SIZE 32768
#define N_Percent 0.03
#define N_Trivial 1500
#define N_Occupy 10
#define MAX_FACTOR 4
#define STD_FACTOR 2
#define InsertIntoCell(k, ipnt) \
{\
assert(InGrid(k)); \
assert(InPntsList(ipnt)); \
\
if (num_pnt_list >= max_num_pnt_list) { \
max_num_pnt_list += PAGE_SIZE; \
pnt_list = (pnt_node*) ReallocateArray(pnt_list, max_num_pnt_list, \
sizeof(pnt_node), \
"grid:pnt_list"); \
} \
\
pnt_list[num_pnt_list].pnt = ipnt; \
pnt_list[num_pnt_list].next = the_grid[k]; \
the_grid[k] = num_pnt_list; \
++num_pnt_list; \
}
#define Convert(i, j, k) \
{ \
k = i * N_y + j; \
}
#define DetermineCell(x, y, i, j) \
{ \
assert(x > grid_min_x); \
i = REAL_TO_INT((x - grid_min_x) / grid_x); \
assert((i >= 0) && (i < max_num_raster_x)); \
if (x < raster_x[i]) { \
--i; \
} \
else if (x > raster_x[i+1]) { \
++i; \
} \
assert((i >= 0) && (i < N_x)); \
\
assert(y > grid_min_y); \
j = REAL_TO_INT((y - grid_min_y) / grid_y); \
if (y < raster_y[j]) { \
--j; \
} \
else if (y > raster_y[j+1]) { \
++j; \
} \
assert((j >= 0) && (j < N_y)); \
}
#define Invert(k, i, j) \
{ \
i = k / N_y; \
assert((0 <= i) && (i < N_x)); \
j = k - i * N_y; \
assert((0 <= j) && (j < N_y)); \
}
#define ScanCell(i, j, k, p, q, ind_pnt, pnt_i, dist, dist_min, i_min) \
{\
Convert(i, j, k); \
assert(InGrid(k)); \
ind_pnt = the_grid[k]; \
\
while (ind_pnt != NIL) { \
pnt_i = pnt_list[ind_pnt].pnt; \
q = GetPntCoords(pnt_i); \
dist = PntPntDistSq(p, q); \
if (dist < dist_min) {\
i_min = pnt_i; \
dist_min = dist; \
} \
ind_pnt = pnt_list[ind_pnt].next; \
} \
}
void vroniObject::FreeGrid(void)
{
FreeMemory((void**) &the_grid, "grid:grid");
raster_x.clear();
raster_y.clear();
FreeMemory((void**) &pnt_list, "grid:pnt_list");
FreeMemory((void**) &first_N_pnts, "grid:first_N_pnts");
N = 0;
N_x = 0;
N_y = 0;
num_pnt_list = 0;
max_num_pnt_list = 0;
max_num_grid = 0;
max_num_raster_x = 0;
max_num_raster_y = 0;
first_counter = 0;
grid_used = false;
return;
}
void vroniObject::InitGrid(void)
{
int k;
assert(N > 0);
assert(N_x * N_y <= N);
if (N > max_num_grid) {
max_num_grid = N;
the_grid = (int*) ReallocateArray(the_grid, max_num_grid, sizeof(int),
"grid:grid");
}
/* */
/* initially, all grid cells are empty. */
/* */
for (k = 0; k < N; ++k) the_grid[k] = NIL;
return;
}
vr_bool vroniObject::InGrid(int cell)
{
return ((0 <= cell) && (cell < max_num_grid));
}
void vroniObject::InitRaster(void)
{
int i, j;
if (N_x >= max_num_raster_x) {
max_num_raster_x = N_x + 1;
gentlyResizeSTLVector(raster_x, max_num_raster_x, "grid:raster_x");
}
if (N_y >= max_num_raster_y) {
max_num_raster_y = N_y + 1;
gentlyResizeSTLVector(raster_y, max_num_raster_y, "grid:raster_y");
}
for (i = 0; i < N_x; ++i) raster_x[i] = grid_min_x + i * grid_x;
raster_x[N_x] = grid_max_x;
for (j = 0; j < N_y; ++j) raster_y[j] = grid_min_y + j * grid_y;
raster_y[N_y] = grid_max_y;
return;
}
void vroniObject::InitPntList(int number)
{
if (number > max_num_pnt_list) {
max_num_pnt_list = number;
pnt_list = (pnt_node*) ReallocateArray(pnt_list, max_num_pnt_list,
sizeof(pnt_node),
"grid:pnt_list");
}
num_pnt_list = 0;
return;
}
#ifdef GRAPHICS
void vroniObject::DrawCell(int i, int j)
{
int i1, j1;
coord p, q;
i1 = i + 1;
j1 = j + 1;
p.x = raster_x[i];
q.x = raster_x[i1];
p.y = raster_y[j];
q.y = raster_y[j];
AddEdgeToBuffer(p.x, p.y, q.x, q.y, GridColor);
p.y = raster_y[j1];
q.y = raster_y[j1];
AddEdgeToBuffer(p.x, p.y, q.x, q.y, GridColor);
p.x = raster_x[i];
q.x = raster_x[i];
p.y = raster_y[j];
q.y = raster_y[j1];
AddEdgeToBuffer(p.x, p.y, q.x, q.y, GridColor);
p.x = raster_x[i1];
q.x = raster_x[i1];
AddEdgeToBuffer(p.x, p.y, q.x, q.y, GridColor);
return;
}
void vroniObject::DrawGrid(void)
{
int i, j;
for (i = 0; i < N_x; ++i) {
for (j = 0; j < N_y; ++j) {
DrawCell(i, j);
}
}
return;
}
#endif
void vroniObject::SetUpGrid(int number)
{
double n_pnt;
assert(number >= 0);
n_pnt = sqrt((double) number);
N_x = REAL_TO_INT(n_pnt * grid_a);
N_y = REAL_TO_INT(n_pnt * grid_b);
if (N_x < 1) N_x = 1;
if (N_y < 1) N_y = 1;
N = N_x * N_y;
grid_x = grid_delta_x / N_x;
grid_y = grid_delta_y / N_y;
InitGrid();
InitRaster();
InitPntList(num_pnts);
grid_used = true;
return;
}
void vroniObject::BuildGrid(void)
{
double eps;
int number;
eps = ZERO * 1000.0;
if (eq(eps, ZERO)) eps = 1.0e-10;
/* */
/* determine the grid size and location */
/* */
grid_delta_x = (bb_max.x - bb_min.x) * 0.01;
grid_delta_y = (bb_max.y - bb_min.y) * 0.01;
if (le(grid_delta_x, eps)) grid_delta_x = eps;
if (le(grid_delta_y, eps)) grid_delta_y = eps;
grid_min_x = bb_min.x - grid_delta_x;
grid_max_x = bb_max.x + grid_delta_x;
grid_min_y = bb_min.y - grid_delta_y;
grid_max_y = bb_max.y + grid_delta_y;
grid_delta_x = grid_max_x - grid_min_x;
grid_delta_y = grid_max_y - grid_min_y;
assert(gt(grid_delta_x, ZERO));
assert(gt(grid_delta_y, ZERO));
/* */
/* initialize the grid raster; we adapt the grid somewhat according to */
/* the aspect ratio of the bounding box of the points. roughly, the */
/* number of cells equals the number of points. */
/* */
grid_a = sqrt(grid_delta_x / grid_delta_y);
if (grid_a < 0.001) {
grid_a = 0.001;
grid_b = 1000.0;
}
else if (grid_a > 1000.0) {
grid_a = 1000.0;
grid_b = 0.001;
}
else {
grid_b = 1.0 / grid_a;
}
assert((N_Percent > 0.0) && (N_Percent <= 1.0));
assert(N_Trivial >= 0);
/* */
/* decide whether to use only one grid for all points, or a smaller grid */
/* for the first few percent of the points, and a larger grid (or a tree) */
/* for the full set of points. */
/* */
number = (int) (num_pnts * N_Percent);
if (number < N_Trivial) number *= 2;
if (number < N_Trivial) number *= 2;
if (number < N_Trivial) number *= 2;
if (number >= N_Trivial) {
N_Initial = number;
first_N_pnts = (int*) ReallocateArray(first_N_pnts, N_Initial + 2,
sizeof(int),
"grid:first_N_pnts");
}
else {
number = num_pnts;
N_Initial = -1;
}
SetUpGrid((int) (number * STD_FACTOR));
first_counter = 0;
return;
}
void vroniObject::DeletePntFromGrid(int index)
{
int i, j, k;
coord p;
assert(InPntsList(index));
p = GetPntCoords(index);
DetermineCell(p.x, p.y, i, j);
Convert(i, j, k);
assert(InGrid(k));
assert(pnt_list[the_grid[k]].pnt == index);
the_grid[k] = pnt_list[the_grid[k]].next;
return;
}
int vroniObject::InsertPntIntoGrid(int index, double *min_dist)
{
int i, j, k, n, m, i_min, pnt_i, ind_pnt;
double x, y, dist, dist_min;
int num;
double factor;
coord p, q;
assert(InPntsList(index));
p = GetPntCoords(index);
DetermineCell(p.x, p.y, i, j);
Convert(i, j, k);
assert(InGrid(k));
dist_min = DBL_MAX;
i_min = NIL;
ScanCell(i, j, k, p, q, ind_pnt, pnt_i, dist, dist_min, i_min);
*min_dist = dist_min;
if (gt(dist_min, ZERO2)) {
InsertIntoCell(k, index);
if (first_counter <= N_Initial) {
assert(first_counter >= 0);
if (first_counter < N_Initial) {
first_N_pnts[first_counter] = index;
}
else {
first_N_pnts[first_counter] = index;
num = 0;
for (n = 0; n < N; ++n) {
if (the_grid[n] != NIL) ++num;
}
assert(num > 0);
factor = N_Initial / ((double) num);
if (factor > N_Occupy) {
grid_used = false;
assert(N_Initial >= 0);
i_min = InitTree(num_pnts, first_N_pnts, N_Initial, min_dist);
FreeGrid();
}
else {
if (factor < STD_FACTOR) factor = STD_FACTOR;
else if (factor > MAX_FACTOR) factor = MAX_FACTOR;
SetUpGrid(REAL_TO_INT(num_pnts * factor));
for (n = 0; n < N_Initial; ++n) {
assert(InPntsList(first_N_pnts[n]));
m = first_N_pnts[n];
x = GetPntCoords(m).x;
y = GetPntCoords(m).y;
DetermineCell(x, y, i, j);
Convert(i, j, k);
assert(InGrid(k));
InsertIntoCell(k, m);
}
assert(index == first_N_pnts[N_Initial]);
DetermineCell(p.x, p.y, i, j);
Convert(i, j, k);
assert(InGrid(k));
ScanCell(i, j, k, p, q, ind_pnt, pnt_i, dist, dist_min, i_min);
*min_dist = dist_min;
if (gt(dist_min, ZERO2)) InsertIntoCell(k, index);
}
}
++first_counter;
}
}
return i_min;
}
int vroniObject::FindNearestNeighborGrid(int index, int i_min, double *min_dist)
{
int i, j, k, ind_pnt, pnt_i, i1, i2, j1, j2, m;
double dist, dist_min, delta;
vr_bool lft, rgt, top, bot;
coord p, q;
assert(InPntsList(index));
dist_min = *min_dist;
p = GetPntCoords(index);
/* */
/* find the cell that contains the point pnts[index] whose nearest */
/* neighbor is sought. */
/* */
DetermineCell(p.x, p.y, i, j);
/* */
/* search the cells around position i, j until a nearest neighbor is */
/* found. we terminate the search as soon as the nearest-neighbor circle */
/* does not intersect any neighboring cell that has not yet been searched.*/
/* */
i1 = i2 = i;
j1 = j2 = j;
lft = rgt = top = bot = true;
while (lft || rgt || top || bot) {
if (lft) {
if (i1 > 0) {
delta = p.x - raster_x[i1];
delta = delta * delta;
if (delta < dist_min) {
--i1;
for (j = j1; j <= j2; ++j) {
ScanCell(i1, j, k, p, q, ind_pnt, pnt_i, dist, dist_min,
i_min);
}
}
else {
lft = false;
}
}
else {
lft = false;
}
}
if (rgt) {
m = i2 + 1;
if (m < N_x) {
delta = p.x - raster_x[m];
delta = delta * delta;
if (delta < dist_min) {
i2 = m;
for (j = j1; j <= j2; ++j) {
ScanCell(i2, j, k, p, q, ind_pnt, pnt_i, dist, dist_min,
i_min);
}
}
else {
rgt = false;
}
}
else {
rgt = false;
}
}
if (bot) {
if (j1 > 0) {
delta = p.y - raster_y[j1];
delta = delta * delta;
if (delta < dist_min) {
--j1;
for (i = i1; i <= i2; ++i) {
ScanCell(i, j1, k, p, q, ind_pnt, pnt_i, dist, dist_min,
i_min);
}
}
else {
bot = false;
}
}
else {
bot = false;
}
}
if (top) {
m = j2 + 1;
if (m < N_y) {
delta = p.y - raster_y[m];
delta = delta * delta;
if (delta < dist_min) {
j2 = m;
for (i = i1; i <= i2; ++i) {
ScanCell(i, j2, k, p, q, ind_pnt, pnt_i, dist, dist_min,
i_min);
}
}
else {
top = false;
}
}
else {
top = false;
}
}
}
assert(i_min != NIL);
assert(dist_min < DBL_MAX);
*min_dist = dist_min;
return i_min;
}
Reference in New Issue View Git Blame Copy Permalink