https://github.com/cran/spatstat
Tip revision: adc4d1a73ce0bc27b15a9ab4b12d1a52f68a4a55 authored by Adrian Baddeley on 19 March 2010, 07:31:32 UTC
version 1.18-1
version 1.18-1
Tip revision: adc4d1a
poly2im.c
/*
poly2im.c
Conversion from (x,y) polygon to pixel image
poly2imI pixel value = 1{pixel centre is inside polygon}
poly2imA pixel value = area of intersection between pixel and polygon
$Revision: 1.1 $ $Date: 2009/08/20 17:00:09 $
*/
#undef DEBUG
#include <math.h>
#ifdef DEBUG
#include <stdio.h>
#endif
#define OUT(I,J) out[I + (J) * Ny]
void
poly2imI(xp, yp, np, nx, ny, out)
double *xp, *yp; /* polygon vertices, anticlockwise, CLOSED */
int *np;
int *nx, *ny; /* INTEGER raster points from (0,0) to (nx-1, ny-1) */
int *out; /* output matrix [ny, nx], initialised to zero */
{
int Np, Nx, Ny, Nxy;
int i, j, k;
double x0, y0, x1, y1, xleft, xright, yleft, yright;
double dx, dy, y, slope, intercept;
int jleft, jright, imax;
int sign;
Np = *np;
Nx = *nx;
Ny = *ny;
Nxy = Nx * Ny;
/* run through polygon edges */
for(k = 0; k < Np-1; k++) {
x0 = xp[k];
y0 = yp[k];
x1 = xp[k+1];
y1 = yp[k+1];
if(x0 < x1) {
xleft = x0;
xright = x1;
yleft = y0;
yright = y1;
sign = -1;
} else {
xleft = x1;
xright = x0;
yleft = y1;
yright = y0;
sign = +1;
}
/* determine relevant columns of pixels */
jleft = (int) ceil(xleft);
jright = (int) floor(xright);
if(jleft < Nx && jright >= 0 && jleft <= jright) {
if(jleft < 0) { jleft = 0; }
if(jright >= Nx) {jright = Nx - 1; }
/* equation of edge */
dx = xright - xleft;
dy = yright - yleft;
slope = dy/dx;
intercept = yleft - slope * xleft;
/* visit relevant columns */
for(j = jleft; j <= jright; j++) {
y = slope * ((double) j) + intercept;
imax = (int) floor(y);
if(imax >= Ny) imax = Ny-1;
if(imax >= 0) {
/* increment entries below edge in this column */
for(i = 0; i <= imax; i++) {
OUT(i,j) += sign;
}
}
}
}
}
}
#define BELOW -1
#define INSIDE 0
#define ABOVE 1
void
poly2imA(ncol, nrow, xpoly, ypoly, npoly, out, status)
int *ncol, *nrow; /* pixels are unit squares from (0,0) to (ncol,nrow) */
double *xpoly, *ypoly; /* vectors of coordinates of polygon vertices */
int *npoly;
double *out; /* double array [nrow, ncol] of pixel areas */
int *status;
{
double *xp, *yp;
int nx, ny, np;
int i, j, k;
double xcur, ycur, xnext, ynext, xleft, yleft, xright, yright;
int sgn, jmin, jmax, imin, imax;
double x0, y0, x1, y1, slope, yhi, ylo, area, xcut, xcutA, xcutB;
int klo, khi;
nx = *ncol;
ny = *nrow;
xp = xpoly;
yp = ypoly;
np = *npoly;
*status = 0;
/* initialise output array */
for(i = 0; i < ny; i++)
for(j = 0; j < nx; j++)
out[j + ny * i] = 0;
/* ............ loop over polygon edges ...................*/
for(k = 0; k < np - 1; k++) {
xcur = xp[k];
ycur = yp[k];
xnext = xp[k+1];
ynext = yp[k+1];
#ifdef DEBUG
fprintf(stderr,
"\nEdge %d from (%lf, %lf) to (%lf, %lf) .........\n",
k, xcur, ycur, xnext, ynext);
#endif
if(xcur != xnext) {
/* vertical edges are ignored */
if(xcur < xnext) {
#ifdef DEBUG
fprintf(stderr, "negative sign\n");
#endif
sgn = -1;
xleft = xcur;
yleft = ycur;
xright = xnext;
yright = ynext;
} else {
#ifdef DEBUG
fprintf(stderr, "positive sign\n");
#endif
sgn = 1;
xleft = xnext;
yleft = ynext;
xright = xcur;
yright = ycur;
}
/* we have now ensured xleft < xright */
slope = (yright - yleft)/(xright - xleft);
/* Find relevant columns of pixels */
jmin = floor(xleft);
jmin = (jmin < 0) ? 0 : jmin;
jmax = ceil(xright);
jmax = (jmax > nx - 1) ? nx - 1 : jmax;
/* Find relevant rows of pixels */
imin = floor((yleft < yright) ? yleft : yright);
imin = (imin < 0) ? 0 : imin;
imax = ceil((yleft < yright) ? yright : yleft);
imax = (imax > ny - 1) ? ny - 1 : imax;
#ifdef DEBUG
fprintf(stderr, "imin=%d, imax=%d, jmin=%d, jmax=%d\n",
imin, imax, jmin, jmax);
#endif
/* ........... loop over columns of pixels ..............*/
for(j = jmin; j <= jmax; j++) {
#ifdef DEBUG
fprintf(stderr, "\t j=%d:\n", j);
#endif
/*
Intersect trapezium with column of pixels
*/
if(xleft <= j+1 && xright >= j) {
if(xleft >= j) {
/* retain left corner */
#ifdef DEBUG
fprintf(stderr, "\tretain left corner\n");
#endif
x0 = xleft;
y0 = yleft;
} else {
/* trim left corner */
#ifdef DEBUG
fprintf(stderr, "\ttrim left corner\n");
#endif
x0 = (double) j;
y0 = yleft + slope * (x0 - xleft);
}
if(xright <= j+1) {
/* retain right corner */
#ifdef DEBUG
fprintf(stderr, "\tretain right corner\n");
#endif
x1 = xright;
y1 = yright;
} else {
/* trim right corner */
#ifdef DEBUG
fprintf(stderr, "\ttrim right corner\n");
#endif
x1 = (double) (j+1);
y1 = yright + slope * (x1 - xright);
}
/* save min and max y */
if(y0 < y1) {
#ifdef DEBUG
fprintf(stderr, "slope %lf > 0\n", slope);
#endif
ylo = y0;
yhi = y1;
} else {
#ifdef DEBUG
fprintf(stderr, "slope %lf <= 0\n", slope);
#endif
ylo = y1;
yhi = y0;
}
/* ............ loop over pixels within column ......... */
/* first part */
if(imin > 0) {
for(i = 0; i < imin; i++) {
#ifdef DEBUG
fprintf(stderr, "\ti=%d:\n", i);
#endif
/*
The trimmed pixel [x0, x1] * [i, i+1]
lies below the polygon edge.
*/
area = (x1 - x0);
#ifdef DEBUG
fprintf(stderr, "\tIncrementing area by %lf\n", sgn * area);
#endif
out[j + ny * i] += sgn * area;
}
}
/* second part */
for(i = imin; i <= imax; i++) {
#ifdef DEBUG
fprintf(stderr, "\ti=%d:\n", i);
#endif
/*
Compute area of intersection between trapezium
and trimmed pixel [x0, x1] x [i, i+1]
*/
klo = (ylo <= i) ? BELOW : (ylo >= (i+1))? ABOVE: INSIDE;
khi = (yhi <= i) ? BELOW : (yhi >= (i+1))? ABOVE: INSIDE;
if(klo == ABOVE) {
/* trapezium covers pixel */
#ifdef DEBUG
fprintf(stderr, "\t\ttrapezium covers pixel\n");
#endif
area = (x1-x0);
} else if(khi == BELOW) {
#ifdef DEBUG
fprintf(stderr, "\t\tpixel avoids trapezium\n");
#endif
/* pixel avoids trapezium */
area = 0.0;
} else if(klo == INSIDE && khi == INSIDE) {
/* polygon edge is inside pixel */
#ifdef DEBUG
fprintf(stderr, "\t\t polygon edge is inside pixel\n");
#endif
area = (x1-x0) * ((ylo + yhi)/2.0 - i);
} else if(klo == INSIDE && khi == ABOVE) {
/* polygon edge crosses upper edge of pixel */
#ifdef DEBUG
fprintf(stderr,
"\t\t polygon edge crosses upper edge of pixel\n");
#endif
xcut = x0 + ((i+1) - y0)/slope;
if(slope > 0)
area = (xcut - x0) * ((y0 + (i+1))/2 - i) + (x1 - xcut);
else
area = (x1 - xcut) * ((y1 + (i+1))/2 - i) + (xcut - x0);
} else if(klo == BELOW && khi == INSIDE) {
/* polygon edge crosses lower edge of pixel */
#ifdef DEBUG
fprintf(stderr, "\t\t polygon edge crosses lower edge of pixel\n");
#endif
xcut = x0 + (i - y0)/slope;
if(slope > 0)
area = (x1 - xcut) * ((y1 + i)/2 - i);
else
area = (xcut - x0) * ((y0 + i)/2 - i);
} else if(klo == BELOW && khi == ABOVE) {
/* polygon edge crosses upper and lower edges of pixel */
#ifdef DEBUG
fprintf(stderr,
"\t\t polygon edge crosses upper and lower edges of pixel\n");
#endif
xcutA = x0 + (i - y0)/slope;
xcutB = x0 + ((i+1) - y0)/slope;
if(slope > 0)
area = (xcutB - xcutA)/2 + (x1 - xcutB);
else
area = (xcutB - x0) + (xcutA - xcutB)/2;
} else {
/* control should not pass to here */
*status = 1;
return;
}
/* add contribution to area of pixel */
#ifdef DEBUG
fprintf(stderr, "\tIncrementing area by %lf\n", sgn * area);
#endif
out[j + ny * i] += sgn * area;
}
/* ............ end of loop over pixels within column ......... */
}
}
/* ........ end of loop over columns of pixels ...............*/
}
}
/* ......... end of loop over polygon edges ...................*/
}