https://github.com/cran/spatstat
Tip revision: 32c7daeb36b6e48fd0356bdcec9580ae124fee5e authored by Adrian Baddeley on 29 December 2015, 22:08:27 UTC
version 1.44-1
version 1.44-1
Tip revision: 32c7dae
areapair.c
/*
areapair.c
$Revision: 1.6 $ $Date: 2013/09/18 04:11:42 $
Specialised code for the second order conditional intensity
of the area-interaction process
*/
#include <R.h>
#include <math.h>
#include "yesno.h"
/* computes area of b(A, r) \int b(B, r) \setminus \bigcup_i b(X[i], r) */
void delta2area(xa, ya, xb, yb,
nother, xother, yother,
radius, epsilon, pixcount)
double *xa, *ya, *xb, *yb;
int *nother;
double *xother, *yother;
double *radius, *epsilon;
int *pixcount;
{
int Ni, Nj, Nk, i, j, k, count, covered;
double xA, yA, xB, yB, r, eps, r2;
double xmin, xmax, ymin, ymax, xi, yj;
double dxA, dyA;
double dxB, dyB;
double dx, dy;
Nk = *nother;
xA = *xa;
yA = *ya;
xB = *xb;
yB = *yb;
r = *radius;
eps = *epsilon;
r2 = r * r;
/* find intersection of squares centred on A and B */
if(xA < xB) {
xmin = xB - r;
xmax = xA + r;
} else {
xmin = xA - r;
xmax = xB + r;
}
if(xmin > xmax) return;
if(yA < yB) {
ymin = yB - r;
ymax = yA + r;
} else {
ymin = yA - r;
ymax = yB + r;
}
if(ymin > ymax) return;
/* set up grid */
Ni = (int) ceil((xmax - xmin)/eps) + 1;
Nj = (int) ceil((ymax - ymin)/eps) + 1;
count = 0;
for(i = 0, xi = xmin; i < Ni; i++, xi += eps) {
dxA = xi - xA;
for(j = 0, yj = ymin; j < Nj; j++, yj += eps) {
dyA = yj - yA;
if(dxA * dxA + dyA * dyA <= r2) {
/* grid point belongs to b(A, r) */
dxB = xi - xB;
dyB = yj - yB;
if(dxB * dxB + dyB * dyB <= r2) {
/* grid point belongs to b(A,r) \cap b(B,r) */
covered = NO;
/* test whether it is covered by another b(X[k], r) */
for(k = 0; k < Nk; k++) {
dx = xi - xother[k];
dy = yj - yother[k];
if(dx * dx + dy * dy <= r2) {
covered = YES;
break;
}
}
if(!covered) {
++count;
}
}
}
}
}
*pixcount = count;
}