Skip to content

Commit

Permalink
Extended distance transform code
Browse files Browse the repository at this point in the history
  • Loading branch information
baddstats committed Aug 29, 2023
1 parent dabf607 commit 1a77225
Show file tree
Hide file tree
Showing 10 changed files with 191 additions and 93 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: spatstat.geom
Version: 3.2-4.002
Date: 2023-08-15
Version: 3.2-4.003
Date: 2023-08-29
Title: Geometrical Functionality of the 'spatstat' Family
Authors@R: c(person("Adrian", "Baddeley",
role = c("aut", "cre", "cph"),
Expand Down
7 changes: 6 additions & 1 deletion NEWS
Original file line number Diff line number Diff line change
@@ -1,13 +1,18 @@
CHANGES IN spatstat.geom VERSION 3.2-4.002
CHANGES IN spatstat.geom VERSION 3.2-4.003

OVERVIEW

o Extension of distance transform algorithm.

o Suppress annoying warnings.

o Further bug fix in quadratcount

SIGNIFICANT USER-VISIBLE CHANGES

o distmap.owin
New argument 'connect'.

o unnormdensity
Suppress annoying warning messages from density.default.
This affects many functions in the spatstat family of packages.
Expand Down
9 changes: 7 additions & 2 deletions R/distmap.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
#
# distmap.R
#
# $Revision: 1.32 $ $Date: 2022/05/21 09:52:11 $
# $Revision: 1.34 $ $Date: 2023/08/28 06:38:54 $
#
#
# Distance transforms
Expand Down Expand Up @@ -45,7 +45,8 @@ distmap.ppp <- function(X, ..., clip=FALSE, metric=NULL) {
return(V)
}

distmap.owin <- function(X, ..., discretise=FALSE, invert=FALSE, metric=NULL) {
distmap.owin <- function(X, ..., discretise=FALSE, invert=FALSE,
connect=8, metric=NULL) {
verifyclass(X, "owin")
uni <- unitname(X)
if(!is.null(metric)) {
Expand Down Expand Up @@ -79,6 +80,9 @@ distmap.owin <- function(X, ..., discretise=FALSE, invert=FALSE, metric=NULL) {
Dist[complement.owin(X, bigbox)] <- 0
}
} else {
check.1.integer(connect)
if(!(connect %in% c(8, 24)))
stop("Argument 'connect' must equal 8 or 24", call.=FALSE)
X <- as.mask(X, ...)
if(invert)
X <- complement.owin(X)
Expand All @@ -92,6 +96,7 @@ distmap.owin <- function(X, ..., discretise=FALSE, invert=FALSE, metric=NULL) {
mat <- rbind(FALSE, mat, FALSE)
## call C routine
res <- .C(SG_distmapbin,
connect=as.integer(connect),
xmin=as.double(xc[1L]),
ymin=as.double(yr[1L]),
xmax=as.double(xc[nc]),
Expand Down
2 changes: 1 addition & 1 deletion inst/doc/packagesizes.txt
Original file line number Diff line number Diff line change
Expand Up @@ -21,4 +21,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
"2023-07-03" "3.2-2" 450 1202 0 35864 15621
"2023-07-20" "3.2-3" 450 1203 0 35915 15747
"2023-07-20" "3.2-4" 450 1203 0 35915 15747
"2023-08-15" "3.2-4.002" 450 1203 0 35919 15747
"2023-08-29" "3.2-4.003" 450 1203 0 35924 15822
7 changes: 6 additions & 1 deletion man/distmap.owin.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,8 @@
in the given window.
}
\usage{
\method{distmap}{owin}(X, \dots, discretise=FALSE, invert=FALSE, metric=NULL)
\method{distmap}{owin}(X, \dots, discretise=FALSE, invert=FALSE,
connect=8, metric=NULL)
}
\arguments{
\item{X}{
Expand All @@ -24,6 +25,10 @@
If \code{TRUE}, compute the distance transform of the
complement of the window.
}
\item{connect}{
Neighbourhood connectivity for the discrete distance transform
algorithm. Either 8 or 24.
}
\item{metric}{
Optional. A distance metric
(object of class \code{"metric"}, see \code{\link{metric.object}})
Expand Down
105 changes: 22 additions & 83 deletions src/distmapbin.c
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,11 @@
distmapbin.c
Distance transform of a discrete binary image
(8-connected path metric)
(8-connected or 24-connected path metric)
$Revision: 1.10 $ $Date: 2022/10/22 09:29:51 $
$Revision: 1.12 $ $Date: 2023/08/28 06:27:24 $
Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner 2001-2018
Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner 2001-2023
Licence: GNU Public Licence >= 2
*/
Expand All @@ -22,89 +21,24 @@ void shape_raster(Raster *ras, void *data,
double xmin, double ymin, double xmax, double ymax,
int nrow, int ncol, int mrow, int mcol);

void
distmap_bin(
Raster *in, /* input: binary image */
Raster *dist /* output: distance to nearest point */
/* rasters must have been dimensioned by shape_raster()
and must all have identical dimensions and margins */
) {
int j,k;
double d, dnew;
double xstep, ystep, diagstep, huge;
int rmin, rmax, cmin, cmax;

/* distances between neighbouring pixels */
xstep = in->xstep;
ystep = in->ystep;
diagstep = sqrt(xstep * xstep + ystep * ystep);
if(xstep < 0) xstep = -xstep;
if(ystep < 0) ystep = -ystep;

/* effectively infinite distance */
huge = 2.0 * Distance(dist->xmin,dist->ymin,dist->xmax,dist->ymax);

/* image boundaries */
rmin = in->rmin;
rmax = in->rmax;
cmin = in->cmin;
cmax = in->cmax;

#define DISTANCE(ROW, COL) Entry(*dist, ROW, COL, double)
#define MASKTRUE(ROW, COL) (Entry(*in, ROW, COL, int) != 0)
#define MASKFALSE(ROW, COL) (Entry(*in, ROW, COL, int) == 0)
#define UPDATE(D, ROW, COL, STEP) \
dnew = STEP + DISTANCE(ROW, COL); \
if(D > dnew) D = dnew

/* initialise edges to boundary condition */
for(j = rmin-1; j <= rmax+1; j++) {
DISTANCE(j, cmin-1) = (MASKTRUE(j, cmin-1)) ? 0.0 : huge;
DISTANCE(j, cmax+1) = (MASKTRUE(j, cmax+1)) ? 0.0 : huge;
}
for(k = cmin-1; k <= cmax+1; k++) {
DISTANCE(rmin-1, k) = (MASKTRUE(rmin-1, k)) ? 0.0 : huge;
DISTANCE(rmax+1, k) = (MASKTRUE(rmax+1, k)) ? 0.0 : huge;
}

/* forward pass */
/* define core algorithms, using template 'distmapbin.h' */

for(j = rmin; j <= rmax; j++) {
R_CheckUserInterrupt();
for(k = cmin; k <= cmax; k++) {
if(MASKTRUE(j, k))
d = DISTANCE(j, k) = 0.0;
else {
d = huge;
UPDATE(d, j-1, k-1, diagstep);
UPDATE(d, j-1, k, ystep);
UPDATE(d, j-1, k+1, diagstep);
UPDATE(d, j, k-1, xstep);
DISTANCE(j,k) = d;
}
}
}
#define FNAME distmap_bin
#define CONNECT 8
#include "distmapbin.h"
#undef FNAME
#undef CONNECT

/* backward pass */

for(j = rmax; j >= rmin; j--) {
R_CheckUserInterrupt();
for(k = cmax; k >= cmin; k--) {
if(MASKFALSE(j,k)) {
d = DISTANCE(j,k);
UPDATE(d, j+1, k+1, diagstep);
UPDATE(d, j+1, k, ystep);
UPDATE(d, j+1, k-1, diagstep);
UPDATE(d, j, k+1, xstep);
DISTANCE(j,k) = d;
}
}
}
}
#define FNAME dist24map_bin
#define CONNECT 24
#include "distmapbin.h"
#undef FNAME
#undef CONNECT

/* R interface */

void distmapbin(
int *connect, /* connectivity: 8 or 24 */
double *xmin,
double *ymin,
double *xmax,
Expand All @@ -120,15 +54,20 @@ void distmapbin(
Raster data, dist, bdist;

void distmap_bin(Raster *in, Raster *dist);
void dist24map_bin(Raster *in, Raster *dist);

shape_raster( &data, (void *) inp, *xmin,*ymin,*xmax,*ymax,
*nr+2, *nc+2, 1, 1);
shape_raster( &dist, (void *) distances,*xmin,*ymin,*xmax,*ymax,
*nr+2,*nc+2,1,1);
shape_raster( &bdist, (void *) boundary, *xmin,*ymin,*xmax,*ymax,
*nr+2,*nc+2,1,1);

distmap_bin(&data, &dist);

if(*connect != 24) {
distmap_bin(&data, &dist);
} else {
dist24map_bin(&data, &dist);
}

dist_to_bdry(&bdist);
}
136 changes: 136 additions & 0 deletions src/distmapbin.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
/*
distmapbin.h
Distance transform of a discrete binary image
Template for core algorithm
This file is #included in 'distmapbin.c' several times
with different values of the macros:
FNAME Name of function
CONNECT Connectivity (8 or 24)
$Revision: 1.2 $ $Date: 2023/08/28 07:37:38 $
Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner 2001-2023
Licence: GNU Public Licence >= 2
*/

#define DISTANCE(ROW, COL) Entry(*dist, ROW, COL, double)
#define MASKTRUE(ROW, COL) (Entry(*in, ROW, COL, int) != 0)
#define MASKFALSE(ROW, COL) (Entry(*in, ROW, COL, int) == 0)
#define UPDATE(D, ROW, COL, STEP) \
dnew = STEP + DISTANCE(ROW, COL); \
if(D > dnew) D = dnew

void
FNAME(
Raster *in, /* input: binary image */
Raster *dist /* output: distance to nearest point */
/* rasters must have been dimensioned by shape_raster()
and must all have identical dimensions and margins */
) {
int j,k;
double d, dnew;
double xstep, ystep, diagstep, huge;
int rmin, rmax, cmin, cmax;

#if (CONNECT == 24)
double LstepXXy, LstepYYx;
#endif

/* distances between neighbouring pixels */
xstep = in->xstep;
ystep = in->ystep;
if(xstep < 0) xstep = -xstep;
if(ystep < 0) ystep = -ystep;
diagstep = sqrt(xstep * xstep + ystep * ystep);
#if (CONNECT == 24)
LstepXXy = sqrt(4.0 * xstep * xstep + ystep * ystep);
LstepYYx = sqrt(xstep * xstep + 4.0 * ystep * ystep);
#endif
/* effectively infinite distance */
huge = 2.0 * Distance(dist->xmin,dist->ymin,dist->xmax,dist->ymax);

/* image boundaries */
rmin = in->rmin;
rmax = in->rmax;
cmin = in->cmin;
cmax = in->cmax;

/* initialise edges to boundary condition */
for(j = rmin-1; j <= rmax+1; j++) {
DISTANCE(j, cmin-1) = (MASKTRUE(j, cmin-1)) ? 0.0 : huge;
DISTANCE(j, cmax+1) = (MASKTRUE(j, cmax+1)) ? 0.0 : huge;
}
for(k = cmin-1; k <= cmax+1; k++) {
DISTANCE(rmin-1, k) = (MASKTRUE(rmin-1, k)) ? 0.0 : huge;
DISTANCE(rmax+1, k) = (MASKTRUE(rmax+1, k)) ? 0.0 : huge;
}

/* forward pass */

for(j = rmin; j <= rmax; j++) {
R_CheckUserInterrupt();
for(k = cmin; k <= cmax; k++) {
if(MASKTRUE(j, k))
d = DISTANCE(j, k) = 0.0;
else {
d = huge;
UPDATE(d, j-1, k-1, diagstep);
UPDATE(d, j-1, k, ystep);
UPDATE(d, j-1, k+1, diagstep);
UPDATE(d, j, k-1, xstep);
#if (CONNECT == 24)
if(j > rmin) {
UPDATE(d, j-2, k-1, LstepYYx);
UPDATE(d, j-2, k+1, LstepYYx);
}
if(k > cmin) {
UPDATE(d, j-1, k-2, LstepXXy);
}
if(k < cmax) {
UPDATE(d, j-1, k+2, LstepXXy);
}
#endif
DISTANCE(j,k) = d;
}
}
}

/* backward pass */

for(j = rmax; j >= rmin; j--) {
R_CheckUserInterrupt();
for(k = cmax; k >= cmin; k--) {
if(MASKFALSE(j,k)) {
d = DISTANCE(j,k);
UPDATE(d, j+1, k+1, diagstep);
UPDATE(d, j+1, k, ystep);
UPDATE(d, j+1, k-1, diagstep);
UPDATE(d, j, k+1, xstep);
#if (CONNECT == 24)
if(j < rmax) {
UPDATE(d, j+2, k-1, LstepYYx);
UPDATE(d, j+2, k+1, LstepYYx);
}
if(k > cmin) {
UPDATE(d, j+1, k-2, LstepXXy);
}
if(k < cmax) {
UPDATE(d, j+1, k+2, LstepXXy);
}
#endif
DISTANCE(j,k) = d;
}
}
}
}

#undef DISTANCE
#undef MASKTRUE
#undef MASKFALSE
#undef UPDATE

2 changes: 1 addition & 1 deletion src/init.c
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ static const R_CMethodDef CEntries[] = {
{"dinfty_R", (DL_FUNC) &dinfty_R, 3},
{"discareapoly", (DL_FUNC) &discareapoly, 12},
{"discs2grid", (DL_FUNC) &discs2grid, 11},
{"distmapbin", (DL_FUNC) &distmapbin, 9},
{"distmapbin", (DL_FUNC) &distmapbin, 10},
{"dwpure", (DL_FUNC) &dwpure, 6},
{"exact_dt_R", (DL_FUNC) &exact_dt_R, 14},
{"fardist2grid", (DL_FUNC) &fardist2grid, 10},
Expand Down
2 changes: 1 addition & 1 deletion src/proto.h
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ void nnXwMD(int *, int *, double *, int *, double *, double *, int *, double *);
void nnXxMD(int *, int *, double *, int *, int *, double *, int *, double *, int *, double *);
void knnXwMD(int *, int *, double *, int *, double *, int *, double *, int *, double *);
void knnXxMD(int *, int *, double *, int *, int *, double *, int *, int *, double *, int *, double *);
void distmapbin(double *, double *, double *, double *, int *, int *, int *, double *, double *);
void distmapbin(int *, double *, double *, double *, double *, int *, int *, int *, double *, double *);
void tabsumweight(int *, double *, double *, int *, double *, double *);
void exact_dt_R(double *, double *, int *, double *, double *, double *, double *, int *, int *, int *, int *, double *, int *, double *);
void ps_exact_dt_R(double *, double *, double *, double *, int *, int *, int *, int *, int *, double *, int *, int *, double *);
Expand Down
Loading

0 comments on commit 1a77225

Please sign in to comment.