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
exactPdist.c
/*
exactPdist.c
`Pseudoexact' distance transform of a discrete binary image
(the closest counterpart to `exactdist.c')
$Revision: 1.12 $ $Date: 2011/05/17 12:27:20 $
*/
#include <math.h>
#include "raster.h"
void dist_to_bdry();
void shape_raster();
void
ps_exact_dt(in, dist, row, col)
Raster *in; /* input: binary image */
Raster *dist; /* output: exact distance to nearest point */
Raster *row; /* output: row index of closest point */
Raster *col; /* output: column index of closest point */
/* rasters must have been dimensioned by shape_raster()
and must all have identical dimensions and margins */
{
int j,k;
double d, x, y;
int r, c;
double dnew;
double huge;
/* double bdiag; */
/* initialise */
#define UNDEFINED -1
#define Is_Defined(I) (I >= 0)
#define Is_Undefined(I) (I < 0)
Clear(*row,int,UNDEFINED)
Clear(*col,int,UNDEFINED)
huge = 2.0 * DistanceSquared(dist->xmin,dist->ymin,dist->xmax,dist->ymax);
Clear(*dist,double,huge)
/* if input pixel is TRUE, set distance to 0 and make pixel point to itself */
for(j = in->rmin; j <= in->rmax; j++)
for(k = in->cmin; k <= in->cmax; k++)
if(Entry(*in, j, k, int) != 0) {
Entry(*dist, j, k, double) = 0.0;
Entry(*row, j, k, int) = j;
Entry(*col, j, k, int) = k;
}
/* how to update the distance values */
#define GETVALUES(ROW,COL) \
x = Xpos(*in, COL); \
y = Ypos(*in, ROW); \
d = Entry(*dist,ROW,COL,double);
#define COMPARE(ROW,COL,RR,CC) \
r = Entry(*row,RR,CC,int); \
c = Entry(*col,RR,CC,int); \
if(Is_Defined(r) && Is_Defined(c) \
&& Entry(*dist,RR,CC,double) < d) { \
dnew = DistanceSquared(x, y, Xpos(*in,c), Ypos(*in,r)); \
if(dnew < d) { \
Entry(*row,ROW,COL,int) = r; \
Entry(*col,ROW,COL,int) = c; \
Entry(*dist,ROW,COL,double) = dnew; \
d = dnew; \
} \
}
/* bound on diagonal step distance squared */
/* bdiag = (in->xstep * in->xstep + in->ystep * in->ystep); */
/* forward pass */
for(j = in->rmin; j <= in->rmax; j++)
for(k = in->cmin; k <= in->cmax; k++) {
GETVALUES(j, k)
COMPARE(j,k, j-1,k-1)
COMPARE(j,k, j-1, k)
COMPARE(j,k, j-1,k+1)
COMPARE(j,k, j, k-1)
}
/* backward pass */
for(j = in->rmax; j >= in->rmin; j--)
for(k = in->cmax; k >= in->cmin; k--) {
GETVALUES(j, k)
COMPARE(j,k, j+1,k+1)
COMPARE(j,k, j+1, k)
COMPARE(j,k, j+1,k-1)
COMPARE(j,k, j, k+1)
}
/* take square roots of distances^2 */
for(j = in->rmax; j >= in->rmin; j--)
for(k = in->cmax; k >= in->cmin; k--)
Entry(*dist,j,k,double) = sqrt(Entry(*dist,j,k,double));
}
/* R interface */
void ps_exact_dt_R(xmin, ymin, xmax, ymax, nr, nc, mr, mc,
inp, distances, rows, cols, boundary)
double *xmin, *ymin, *xmax, *ymax; /* x, y dimensions */
int *nr, *nc; /* raster dimensions
EXCLUDING margins */
int *mr, *mc; /* margins */
int *inp; /* input: binary image */
double *distances; /* output: distance to nearest point */
int *rows; /* output: row of nearest point (start= 0) */
int *cols; /* output: column of nearest point (start = 0) */
double *boundary; /* output: distance to boundary of rectangle */
/* all images must have identical dimensions including a margin of 1 on each side */
{
Raster data, dist, row, col, bdist;
int mrow, mcol, nrow, ncol;
mrow = *mr;
mcol = *mc;
/* full dimensions */
nrow = *nr + 2 * mrow;
ncol = *nc + 2 * mcol;
shape_raster( &data, (void *) inp, *xmin,*ymin,*xmax,*ymax,
nrow, ncol, mrow, mcol);
shape_raster( &dist, (void *) distances, *xmin,*ymin,*xmax,*ymax,
nrow, ncol, mrow, mcol);
shape_raster( &row, (void *) rows, *xmin,*ymin,*xmax,*ymax,
nrow, ncol, mrow, mcol);
shape_raster( &col, (void *) cols, *xmin,*ymin,*xmax,*ymax,
nrow, ncol, mrow, mcol);
shape_raster( &bdist, (void *) boundary, *xmin,*ymin,*xmax,*ymax,
nrow, ncol, mrow, mcol);
ps_exact_dt(&data, &dist, &row, &col);
dist_to_bdry(&bdist);
}