knndist.h
/*
knndist.h
Code template for C functions supporting knndist and knnwhich
THE FOLLOWING CODE ASSUMES THAT y IS SORTED IN ASCENDING ORDER
This code is #included multiple times in knndistance.c
Variables used:
FNAME function name
DIST #defined if function returns distance to nearest neighbour
WHICH #defined if function returns id of nearest neighbour
Either or both DIST and WHICH may be defined.
Copyright (C) Adrian Baddeley, Jens Oehlschlagel and Rolf Turner 2000-2012
Licence: GPL >= 2
$Revision: 1.3 $ $Date: 2013/05/27 02:09:10 $
*/
void FNAME(n, kmax, x, y,
#ifdef DIST
nnd,
#endif
#ifdef WHICH
nnwhich,
#endif
huge)
/* inputs */
int *n, *kmax;
double *x, *y, *huge;
/* output matrices (npoints * kmax) in ROW MAJOR order */
#ifdef DIST
double *nnd;
#endif
#ifdef WHICH
int *nnwhich;
#endif
{
int npoints, maxchunk, nk, nk1, i, k, k1, left, right, unsorted;
double d2, d2minK, xi, yi, dx, dy, dy2, hu, hu2, tmp;
double *d2min;
#ifdef WHICH
int *which;
int itmp;
#endif
hu = *huge;
hu2 = hu * hu;
npoints = *n;
nk = *kmax;
nk1 = nk - 1;
/*
create space to store the nearest neighbour distances and indices
for the current point
*/
d2min = (double *) R_alloc((size_t) nk, sizeof(double));
#ifdef WHICH
which = (int *) R_alloc((size_t) nk, sizeof(int));
#endif
/* loop in chunks of 2^16 */
i = 0; maxchunk = 0;
while(i < npoints) {
R_CheckUserInterrupt();
maxchunk += 65536;
if(maxchunk > npoints) maxchunk = npoints;
for(; i < maxchunk; i++) {
#ifdef SPATSTAT_DEBUG
Rprintf("\ni=%d\n", i);
#endif
/* initialise nn distances and indices */
d2minK = hu2;
for(k = 0; k < nk; k++) {
d2min[k] = hu2;
#ifdef WHICH
which[k] = -1;
#endif
}
xi = x[i];
yi = y[i];
/* search backward */
for(left = i - 1; left >= 0; --left)
{
#ifdef SPATSTAT_DEBUG
Rprintf("L");
#endif
dy = yi - y[left];
dy2 = dy * dy;
if(dy2 > d2minK)
break;
dx = x[left] - xi;
d2 = dx * dx + dy2;
if (d2 < d2minK) {
/* overwrite last entry */
d2min[nk1] = d2;
#ifdef WHICH
which[nk1] = left;
#endif
/* bubble sort */
unsorted = YES;
for(k = nk1; unsorted && k > 0; k--) {
k1 = k - 1;
if(d2min[k] < d2min[k1]) {
/* swap entries */
tmp = d2min[k1];
d2min[k1] = d2min[k];
d2min[k] = tmp;
#ifdef WHICH
itmp = which[k1];
which[k1] = which[k];
which[k] = itmp;
#endif
} else {
unsorted = NO;
}
}
/* adjust maximum distance */
d2minK = d2min[nk1];
}
}
/* search forward */
for(right = i + 1; right < npoints; ++right)
{
#ifdef SPATSTAT_DEBUG
Rprintf("R");
#endif
dy = y[right] - yi;
dy2 = dy * dy;
if(dy2 > d2minK)
break;
dx = x[right] - xi;
d2 = dx * dx + dy2;
if (d2 < d2minK) {
/* overwrite last entry */
d2min[nk1] = d2;
#ifdef WHICH
which[nk1] = right;
#endif
/* bubble sort */
unsorted = YES;
for(k = nk1; unsorted && k > 0; k--) {
k1 = k - 1;
if(d2min[k] < d2min[k1]) {
/* swap entries */
tmp = d2min[k1];
d2min[k1] = d2min[k];
d2min[k] = tmp;
#ifdef WHICH
itmp = which[k1];
which[k1] = which[k];
which[k] = itmp;
#endif
} else {
unsorted = NO;
}
}
/* adjust maximum distance */
d2minK = d2min[nk1];
}
}
/* search finished for point i */
#ifdef SPATSTAT_DEBUG
Rprintf("\n");
#endif
/* copy nn distances for point i
to output matrix in ROW MAJOR order
*/
for(k = 0; k < nk; k++) {
#ifdef DIST
nnd[nk * i + k] = sqrt(d2min[k]);
#endif
#ifdef WHICH
nnwhich[nk * i + k] = which[k] + 1; /* R indexing */
#endif
}
/* end of i loop */
}
}
}