https://github.com/cran/gstat
Tip revision: a4ba0967b3c7f05b70d7e3f625228786cc4290e5 authored by Edzer Pebesma on 06 April 2023, 09:32:40 UTC
version 2.1-1
version 2.1-1
Tip revision: a4ba096
select.c
/*
* select.c: neighborhood selection for local prediction
*/
#include <stdio.h>
#include <stdlib.h> /* qsort() */
#include <string.h> /* memcpy() */
#include <math.h> /* sqrt(), fabs() */
#include "defs.h"
void Rprintf(const char *,...);
#include "userio.h"
#include "data.h"
#include "utils.h"
#include "debug.h"
#include "vario.h"
#include "defaults.h"
#include "glvars.h"
#include "nsearch.h"
#include "select.h"
static int octant_select(DATA *d, DPOINT *where);
static int which_octant(DPOINT *where, DPOINT *p, int mode);
int CDECL dist_cmp(const DPOINT **ap, const DPOINT **bp);
static void print_selection(DATA *d, DPOINT *where);
/* beware-of-side-effect macro: don't call with ++/--'s */
#define DPSWAP(a,b) { if (a != b) { tmp = a; a = b; b = tmp; }}
static int select_qtree(DATA *d, DPOINT *where)
{
if (d->n_list <= 0 || d->id < 0 || d->sel_max == 0)
return (d->n_sel = 0);
/*
* global neighbourhood?
*/
if (IS_GLOBAL(d) || where == NULL) {
d->sel = d->list;
d->n_sel = d->n_sel_max = d->n_list;
if (DEBUG_SEL) {
print_selection(d, where);
}
return d->n_sel;
}
/* set up or resize selection list, d->sel */
if (d->sel == NULL || d->sel == d->list) {
/* first time or after forced global selection: */
d->sel = (DPOINT **) emalloc(d->n_list * sizeof(DPOINT *));
d->n_sel_max = d->n_list;
} else if (d->n_list > d->n_sel_max) { /* no, something has changed... */
d->n_sel_max += MAX(MAX_DATA, d->n_list - d->n_sel_max);
d->sel = (DPOINT **) erealloc(d->sel, d->n_sel_max * sizeof(DPOINT *));
}
if (d->id > 0) { /* 2nd, 3d,..., n_vars-th variable */
if (gl_coincide == DEF_coincide)
gl_coincide = decide_on_coincide(); /* establish first ... */
if (gl_coincide) {
int i;
DATA **data = get_gstat_data();
d->n_sel = data[0]->n_sel;
for (i = 0; i < d->n_sel; i++) /* copy previous selection: */
d->sel[i] = d->list[GET_INDEX(data[0]->sel[i])];
if (DEBUG_SEL)
print_selection(d, where);
return d->n_sel; /* and we're done!! */
}
}
/*
* So far, no selection has been done.
* Let's see if today's job is easily done:
*/
memcpy(d->sel, d->list, d->n_list * sizeof(DPOINT *));
if (d->sel_rad >= DBL_MAX && d->sel_max >= d->n_list && d->oct_max == 0) {
d->n_sel = d->n_list;
if (DEBUG_SEL)
print_selection(d, where);
return d->n_sel;
}
/*
* now we're sure that (a) sel_rad is set, or (b) sel_max < n_list,
* or (c) oct_max is set, so let's do the smart thing:
*/
qtree_select(where, d);
return -1; /* more work to do */
}
int select_at(DATA *d, DPOINT *where) {
/*
* fill the array d->sel appropriatly given x,y,z,d->sel_min,d->sel_max
* and d->sel_rad: possibly by corresponding semivariance value
* first select all points within a distance d->sel_rad from where
* then select at max the d->sel_max nearest and return no selection
* if there are less then d->sel_min
* if "FORCE", then select ALWAYS the d->sel_min nearest points.
*
* corrected variogram distance feb. 16th 1993
* changed search to normalizing to (0,0,0) first, aug 1993
*/
if (d->what_is_u == U_UNKNOWN)
d->what_is_u = U_ISDIST; /* we're going to fill this right now */
else if (d->what_is_u != U_ISDIST)
ErrMsg(ER_IMPOSVAL, "select_at() needs distances");
if (select_qtree(d,where) != -1)
return d->n_sel;
/*
* so, now we're at the stage where one of the following conditions holds:
* (a) we selected all points within d->sel_rad
* (b) we selected (at least) the nearest d->sel_max
* (c) we selected (forced) at least d->sel_min, possibly beyond d->sel_rad
* Now, should we select further, sorting on distance?
*
*/
if (d->vdist) { /* use variogram distance as sort criterium */
int i;
for (i = 0; i < d->n_sel; i++)
d->sel[i]->u.dist2 = get_semivariance(get_vgm(LTI(d->id,d->id)),
where->x - d->sel[i]->x,
where->y - d->sel[i]->y,
where->z - d->sel[i]->z);
}
if (d->oct_max) { /* do octant selection */
d->oct_filled = octant_select(d, where);
/* sorts, adjusts n_sel and destroys distance order, so only */
if (get_method() == SPREAD) /* then we've got to re-order them */
qsort(d->sel, (size_t) d->n_sel, sizeof(DPOINT *), (int
CDECL (*)(const void *,const void *)) dist_cmp);
}
if (d->vdist) {
qsort(d->sel, (size_t) d->n_sel, sizeof(DPOINT *), (int
CDECL (*)(const void *, const void *)) dist_cmp);
/* pick d->sel_[max|min] nearest: */
if (d->sel_min && d->n_sel == d->sel_min
&& d->sel[d->n_sel]->u.dist2 > d->sel_rad) /* we forced: */
d->n_sel = d->sel_min;
if (d->n_sel > d->sel_max)
d->n_sel = d->sel_max;
}
if (DEBUG_SEL)
print_selection(d, where);
return d->n_sel;
}
static int octant_select(DATA *d, DPOINT *where) {
/*
* selects the d->oct_max nearest points per octant/quadrant/secant,
* using euclidian or variogram distance
*/
int i, j, noct = 1, n, start = 0, end = 0, n_notempty = 0;
DPOINT **sel, *tmp;
if (d->mode & Z_BIT_SET)
noct = 8;
else if (d->mode & Y_BIT_SET)
noct = 4;
else
noct = 2;
sel = d->sel;
for (i = 0; i < noct; i++) { /* for each octant: */
/* start == end: */
for (j = start; j < d->n_sel; j++) { /* for the remaining of sel: */
if (which_octant(where, sel[j], d->mode) == i) {
DPSWAP(sel[end], sel[j]);
end++; /* j >= end */
}
}
n = end - start; /* # of pts in octant i */
if (n > 0)
n_notempty++; /* Yahoo, another non-empty octant! */
if (n > d->oct_max) {
/* to get the closest n: sort sel from start to end: */
qsort(sel + start, (size_t) n, sizeof(DPOINT *),
(int CDECL (*)(const void *, const void *)) dist_cmp);
/* swap the remaining ones to the end of sel and forget about'm: */
for (j = start + d->oct_max; j < end; j++) {
d->n_sel--;
DPSWAP(sel[j], sel[d->n_sel]);
}
/* accept the first d->oct_max: */
start += d->oct_max;
/* proceed with the next octant: */
end = start;
} else /* accept all n: */
start = end;
}
if (end != d->n_sel) {
Rprintf("end: %d, n_sel: %d\n", end, d->n_sel);
ErrMsg(ER_IMPOSVAL, "octant_select(): remaining points");
}
return n_notempty; /* # non-empty octants */
}
static int which_octant(DPOINT *where, DPOINT *p, int mode) {
/*
* it's pretty hard to get this right for perfectly alligned 3D data:
* in case of omax=1 and only one point at exactly equal distances
* in each of the 6 directions, not all 6 are taken.
* For 2D (omax=1, 4 points) this works, however.
*/
double dx, dy, dz;
int x = 0, y = 0, z = 0;
dx = p->x - where->x; /* < 0 : p in west half */
dy = p->y - where->y; /* < 0 : p in south half */
dz = p->z - where->z; /* < 0 : p in lower half */
if (mode & Z_BIT_SET)
z = (dz < 0);
if (mode & Y_BIT_SET) {
x = (dy < 0 ? dx > 0 : dx >= 0);
y = (dx < 0 ? dy >= 0 : dy > 0);
} else
x = (where->x < p->x);
return (x | (y << 1) | (z << 2));
}
int CDECL dist_cmp(const DPOINT **pa, const DPOINT **pb) {
/* ANSI qsort() conformant dist_cmp */
if ( (*pa)->u.dist2 < (*pb)->u.dist2 )
return -1;
if ( (*pa)->u.dist2 > (*pb)->u.dist2 )
return 1;
return 0;
}
static void print_selection(DATA *d, DPOINT *where) {
/* Add this statement to filter out
* empty selections
if (!d->n_sel)
return;
*/
if (where) {
printlog("selection at ");
logprint_point(where, d);
} else
printlog("(NULL selection location)");
print_data_selection(d);
}