Revision d606122dc24b56ecf537d55eda38f4bf5ac4de1f authored by Adrian Baddeley on 25 October 2010, 10:40:51 UTC, committed by cran-robot on 25 October 2010, 10:40:51 UTC
1 parent 66bc933
areadiff.c
/*
areadiff.c
Area difference function
$Revision: 1.7 $ $Date: 2009/08/12 00:20:28 $
A(x,r) = area of disc b(0,r) not covered by discs b(x_i,r) for x_i in x
Area estimated by point-counting on a fine grid
For use in area-interaction model and related calculations
*/
#undef DEBUG
#ifdef DEBUG
#include <stdio.h>
#endif
/*
Original version areadiff()
1 point u
No trimming of discs
*/
void
areadiff(rad,x,y,nn,ngrid,answer)
/* inputs */
double *rad; /* radius */
double *x, *y; /* coordinate vectors for point pattern */
int *nn; /* length of vectors x and y */
int *ngrid; /* dimensions of point-counting grid */
/* output */
double *answer; /* computed area */
{
double dx, dy, xg, yg, r, r2, xdif, ydif;
int i, j, k, m, n, count, found;
r = *rad;
r2 = r * r;
n = *nn;
m = *ngrid;
dx = dy = 2 * r / (m-1);
count = 0;
/* run through grid points */
for(i = 0, xg = -r; i < m; i++, xg += dx)
for(j = 0, yg = -r; j < m; j++, yg += dy)
/* test for inside disc */
if(xg * xg + yg * yg < r2) {
#ifdef DEBUG
fprintf(stderr, "\n\n (xg,yg) = (%lf, %lf)\n", xg, yg);
#endif
/* run through data points seeking one close to (xy, yg) */
found = 0;
if(n > 0)
for(k = 0; k < n && found == 0; k++) {
#ifdef DEBUG
fprintf(stderr, "(x[%d],y[%d]) = (%lf,%lf)\n", k, k, x[k], y[k]);
#endif
xdif = x[k] - xg;
ydif = y[k] - yg;
if(xdif * xdif + ydif * ydif < r2) {
found = 1;
#ifdef DEBUG
fprintf(stderr, "(x[%d], y[%d]) = (%lf, %lf) covers!\n",
k, k, x[k], y[k]);
#endif
}
}
if(found == 0) {
++count;
#ifdef DEBUG
fprintf(stderr, "---------------Incrementing count\n");
#endif
}
}
#ifdef DEBUG
fprintf(stderr, "Count = %d\n", count);
#endif
/* calculate area */
*answer = ((double) count) * dx * dy;
}
/* similar function, handles multiple values of 'r' */
void
areadifs(rad,nrads,x,y,nxy,ngrid,answer)
/* inputs */
double *rad; /* vector of radii */
int *nrads; /* length of 'rads' */
double *x, *y; /* coordinate vectors for point pattern */
int *nxy; /* length of vectors x and y */
int *ngrid; /* dimensions of point-counting grid */
/* output */
double *answer; /* computed areas (vector of length 'nrads') */
{
double dx, dy, xg, yg, r, r2, xdif, ydif;
int i, j, k, l, m, n, nr, count, found;
n = *nxy;
nr = *nrads;
m = *ngrid;
/* run through radii */
for(l = 0; l < nr; l++) {
r = rad[l];
if(r == 0.0) {
answer[l] = 0.0;
} else {
r2 = r * r;
dx = dy = 2 * r / (m-1);
count = 0;
/* run through grid points */
for(i = 0, xg = -r; i < m; i++, xg += dx)
for(j = 0, yg = -r; j < m; j++, yg += dy)
/* test for inside disc */
if(xg * xg + yg * yg < r2) {
#ifdef DEBUG
fprintf(stderr, "\n\n (xg,yg) = (%lf, %lf)\n", xg, yg);
#endif
/* run through data points seeking one close to (xy, yg) */
found = 0;
if(n > 0)
for(k = 0; k < n && found == 0; k++) {
#ifdef DEBUG
fprintf(stderr, "(x[%d],y[%d]) = (%lf,%lf)\n", k, k, x[k], y[k]);
#endif
xdif = x[k] - xg;
ydif = y[k] - yg;
if(xdif * xdif + ydif * ydif < r2) {
found = 1;
#ifdef DEBUG
fprintf(stderr, "(x[%d], y[%d]) = (%lf, %lf) covers!\n",
k, k, x[k], y[k]);
#endif
}
}
if(found == 0) {
++count;
#ifdef DEBUG
fprintf(stderr, "---------------Incrementing count\n");
#endif
}
}
/* end of loop through data points */
#ifdef DEBUG
fprintf(stderr, "Count = %d\n", count);
#endif
/* calculate area for this value of r*/
answer[l] = ((double) count) * dx * dy;
}
/* end of if(r==0).. else {...} */
}
/* end of loop over r values */
}
/*
Modified version
multiple test points u
discs constrained inside a rectangle
*/
void
areaBdif(rad,nrads,x,y,nxy,ngrid,x0,y0,x1,y1,answer)
/* inputs */
double *rad; /* vector of radii */
int *nrads; /* length of 'rads' */
double *x, *y; /* coordinate vectors for point pattern */
int *nxy; /* length of vectors x and y */
int *ngrid; /* dimensions of point-counting grid */
double *x0,*y0,*x1,*y1; /* constraint rectangle */
/* output */
double *answer; /* computed areas (vector of length 'nrads') */
{
double dx, dy, xg, yg, r, r2, xdif, ydif, xmin, ymin, xmax, ymax;
int i, j, k, l, m, n, nr, count, found;
n = *nxy;
nr = *nrads;
m = *ngrid;
xmin = *x0;
ymin = *y0;
xmax = *x1;
ymax = *y1;
/* run through radii */
for(l = 0; l < nr; l++) {
r = rad[l];
if(r == 0.0) {
answer[l] = 0.0;
} else {
r2 = r * r;
dx = dy = 2 * r / (m-1);
count = 0;
/* run through grid points */
for(i = 0, xg = -r; i < m; i++, xg += dx) {
if(xg >= xmin && xg <= xmax) {
for(j = 0, yg = -r; j < m; j++, yg += dy) {
if(yg >= ymin && yg <= ymax) {
/* test for inside disc */
if(xg * xg + yg * yg < r2) {
#ifdef DEBUG
fprintf(stderr, "\n\n (xg,yg) = (%lf, %lf)\n", xg, yg);
#endif
/* run through data points seeking one close to (xy, yg) */
found = 0;
if(n > 0)
for(k = 0; k < n && found == 0; k++) {
#ifdef DEBUG
fprintf(stderr, "(x[%d],y[%d]) = (%lf,%lf)\n",
k, k, x[k], y[k]);
#endif
xdif = x[k] - xg;
ydif = y[k] - yg;
if(xdif * xdif + ydif * ydif < r2) {
found = 1;
#ifdef DEBUG
fprintf(stderr, "(x[%d], y[%d]) = (%lf, %lf) covers!\n",
k, k, x[k], y[k]);
#endif
}
}
if(found == 0) {
++count;
#ifdef DEBUG
fprintf(stderr, "---------------Incrementing count\n");
#endif
}
}
}
}
}
}
/* end of loop through data points */
#ifdef DEBUG
fprintf(stderr, "Count = %d\n", count);
#endif
/* calculate area for this value of r*/
answer[l] = ((double) count) * dx * dy;
}
/* end of if(r==0).. else {...} */
}
/* end of loop over r values */
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...