Extended 'is.close'

isclose.h Accelerated in 3D case, and avoided GCC 323. Added code for periodic distance isclose.h isclose.c isclose.R is.close.Rd New arguments 'periodic' and 'sorted'. tests/nndist.R Added tests for validity of is.close DESCRIPTION NEWS updated as version 1.47-0.002
spatstat · Oct 15, 2016 · b08a942 · b08a942
1 parent 31a16a6
commit b08a942
Show file tree

Hide file tree

Showing 9 changed files with 377 additions and 84 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: spatstat
-Version: 1.47-0.001
+Version: 1.47-0.002
 Nickname: Recurring Dream
-Date: 2016-10-14
+Date: 2016-10-15
 Title: Spatial Point Pattern Analysis, Model-Fitting, Simulation, Tests
 Author: Adrian Baddeley <Adrian.Baddeley@curtin.edu.au>,
 	Rolf Turner <r.turner@auckland.ac.nz> 

diff --git a/NAMESPACE b/NAMESPACE
@@ -113,8 +113,11 @@ useDynLib(spatstat,
 "inxyp",
 "irevcumsum",
 "isX3close",
+"isX3pclose",
 "isXclose",
+"isXpclose",
 "isXYclose",
+"isXYpclose",
 "KborderD",
 "KborderI",
 "knnd3D",

diff --git a/NEWS b/NEWS
@@ -1,5 +1,5 @@
 
-        CHANGES IN spatstat VERSION 1.47-0.001
+        CHANGES IN spatstat VERSION 1.47-0.002
 
 OVERVIEW
 

diff --git a/R/isclose.R b/R/isclose.R
@@ -3,92 +3,173 @@
 #'
 #'    Determine whether each point has a close neighbour
 #'
-#'    $Revision: 1.4 $  $Date: 2016/10/14 07:36:22 $
+#'    $Revision: 1.7 $  $Date: 2016/10/15 09:25:25 $
 
-is.close <- function(X, r, Y=NULL) {
+is.close <- function(X, r, Y=NULL, ...) {
   UseMethod("is.close")
 }
 
-is.close.default <- function(X, r, Y=NULL) {
-  nd <- if(is.null(Y)) nndist(X) else nncross(X, Y, what="dist")
-  return(nd <= r)
+is.close.default <- function(X, r, Y=NULL, ..., periodic=FALSE) {
+  trap.extra.arguments(...)
+  if(!periodic) {
+    nd <- if(is.null(Y)) nndist(X) else nncross(X, Y, what="dist")
+    return(nd <= r)
+  }
+  if(is.null(Y)) {
+    pd <- pairdist(X, periodic=TRUE)
+    diag(pd) <- Inf
+  } else {
+    pd <- crossdist(X, Y, periodic=TRUE)
+  }
+  return(apply(pd <= r, 1, any))
 }
 
-is.close.ppp <- function(X, r, Y=NULL) {
+is.close.ppp <- function(X, r, Y=NULL, ..., periodic=FALSE, sorted=FALSE) {
+  trap.extra.arguments(...)
   nX <- npoints(X)
   if(nX <= 1) return(logical(nX))
   #' sort by increasing x coordinate
   cX <- coords(X)
-  oo <- order(cX$x)
-  cX <- cX[oo,,drop=FALSE]
+  if(!sorted) {
+    oo <- order(cX$x)
+    cX <- cX[oo,,drop=FALSE]
+  }
   if(is.null(Y)) {
-    zz <- .C("isXclose",
-             n = as.integer(nX),
-             x = as.double(cX$x),
-             y = as.double(cX$y),
-             r = as.double(r),
-             t = as.integer(integer(nX)))
-    tt <- as.logical(zz$t)
+    if(!periodic) {
+      zz <- .C("isXclose",
+               n = as.integer(nX),
+               x = as.double(cX$x),
+               y = as.double(cX$y),
+               r = as.double(r),
+               t = as.integer(integer(nX)))
+    } else {
+      b <- sidelengths(Frame(X))
+      zz <- .C("isXpclose",
+               n = as.integer(nX),
+               x = as.double(cX$x),
+               y = as.double(cX$y),
+               r = as.double(r),
+               b = as.double(b),
+               t = as.integer(integer(nX)))
+    }
   } else {
     stopifnot(is.ppp(Y))
-    #' sort Y by increasing x coordinate
-    cY <- coords(Y)
-    ooY <- order(cY$x)
-    cY <- cY[ooY, , drop=FALSE]
     nY <- npoints(Y)
-    zz <- .C("isXYclose",
-             n1 = as.integer(nX),
-             x1 = as.double(cX$x),
-             y1 = as.double(cX$y),
-             n2 = as.integer(nY),
-             x2 = as.double(cY$x),
-             y2 = as.double(cY$y),
-             r = as.double(r),
-             t = as.integer(integer(nX)))
-    tt <- as.logical(zz$t)
+    if(nY == 0) return(logical(nX))
+    cY <- coords(Y)
+    #' sort Y by increasing x coordinate
+    if(!sorted) {
+      ooY <- order(cY$x)
+      cY <- cY[ooY, , drop=FALSE]
+    }
+    if(!periodic) {
+      zz <- .C("isXYclose",
+               n1 = as.integer(nX),
+               x1 = as.double(cX$x),
+               y1 = as.double(cX$y),
+               n2 = as.integer(nY),
+               x2 = as.double(cY$x),
+               y2 = as.double(cY$y),
+               r = as.double(r),
+               t = as.integer(integer(nX)))
+    } else {
+      bX <- sidelengths(Frame(X))
+      bY <- sidelengths(Frame(Y))
+      if(any(bX != bY))
+        warning("Windows are not equal: periodic distance may be erroneous")
+      zz <- .C("isXYpclose",
+               n1 = as.integer(nX),
+               x1 = as.double(cX$x),
+               y1 = as.double(cX$y),
+               n2 = as.integer(nY),
+               x2 = as.double(cY$x),
+               y2 = as.double(cY$y),
+               r = as.double(r),
+               b = as.double(bX),
+               t = as.integer(integer(nX)))
+    }
   }
+  tt <- as.logical(zz$t)
+  if(sorted) return(tt)
   #' reinstate original order
   ans <- logical(nX)
   ans[oo] <- tt
   return(ans)
 }
 
-is.close.pp3 <- function(X, r, Y=NULL) {
+is.close.pp3 <- function(X, r, Y=NULL, ..., periodic=FALSE, sorted=FALSE) {
+  trap.extra.arguments(...)
   nX <- npoints(X)
   if(nX <= 1) return(logical(nX))
-  #' sort by increasing x coordinate
   cX <- coords(X)
-  oo <- order(cX$x)
-  cX <- cX[oo,,drop=FALSE]
+  if(!sorted) {
+    #' sort by increasing x coordinate
+    oo <- order(cX$x)
+    cX <- cX[oo,,drop=FALSE]
+  }
   if(is.null(Y)) {
-    zz <- .C("isX3close",
-             n = as.integer(nX),
-             x = as.double(cX$x),
-             y = as.double(cX$y),
-             z = as.double(cX$z),
-             r = as.double(r),
-             t = as.integer(integer(nX)))
-    tt <- as.logical(zz$t)
+    if(!periodic) {
+      zz <- .C("isX3close",
+               n = as.integer(nX),
+               x = as.double(cX$x),
+               y = as.double(cX$y),
+               z = as.double(cX$z),
+               r = as.double(r),
+               t = as.integer(integer(nX)))
+    } else {
+      b <- sidelengths(as.box3(X))
+      zz <- .C("isX3pclose",
+               n = as.integer(nX),
+               x = as.double(cX$x),
+               y = as.double(cX$y),
+               z = as.double(cX$z),
+               r = as.double(r),
+               b = as.double(b), 
+               t = as.integer(integer(nX)))
+    }
   } else {
     stopifnot(is.pp3(Y))
-    #' sort Y by increasing x coordinate
-    cY <- coords(Y)
-    ooY <- order(cY$x)
-    cY <- cY[ooY, , drop=FALSE]
     nY <- npoints(Y)
-    zz <- .C("isXYclose",
-             n1 = as.integer(nX),
-             x1 = as.double(cX$x),
-             y1 = as.double(cX$y),
-             z1 = as.double(cX$z),
-             n2 = as.integer(nY),
-             x2 = as.double(cY$x),
-             y2 = as.double(cY$y),
-             z2 = as.double(cY$z),
-             r = as.double(r),
-             t = as.integer(integer(nX)))
-    tt <- as.logical(zz$t)
+    if(nY == 0) return(logical(nX))
+    cY <- coords(Y)
+    if(!sorted) {
+      #' sort Y by increasing x coordinate
+      ooY <- order(cY$x)
+      cY <- cY[ooY, , drop=FALSE]
+    }
+    if(!periodic) {
+      zz <- .C("isXYclose",
+               n1 = as.integer(nX),
+               x1 = as.double(cX$x),
+               y1 = as.double(cX$y),
+               z1 = as.double(cX$z),
+               n2 = as.integer(nY),
+               x2 = as.double(cY$x),
+               y2 = as.double(cY$y),
+               z2 = as.double(cY$z),
+               r = as.double(r),
+               t = as.integer(integer(nX)))
+    } else {
+      bX <- sidelengths(as.box3(X))
+      bY <- sidelengths(as.box3(Y))
+      if(any(bX != bY))
+        warning("Domains are not equal: periodic distance may be erroneous")
+      zz <- .C("isXYclose",
+               n1 = as.integer(nX),
+               x1 = as.double(cX$x),
+               y1 = as.double(cX$y),
+               z1 = as.double(cX$z),
+               n2 = as.integer(nY),
+               x2 = as.double(cY$x),
+               y2 = as.double(cY$y),
+               z2 = as.double(cY$z),
+               r = as.double(r),
+               b = as.double(bX),
+               t = as.integer(integer(nX)))
+    }
   }
+  tt <- as.logical(zz$t)
+  if(sorted) return(tt)
   #' reinstate original order
   ans <- logical(nX)
   ans[oo] <- tt

diff --git a/man/is.close.Rd b/man/is.close.Rd
@@ -11,13 +11,13 @@
   point has a close neighbour in the same pattern.
 }
 \usage{
-  is.close(X, r, Y=NULL)
+  is.close(X, r, Y=NULL, \dots)
 
-  \method{is.close}{default}(X,r, Y=NULL)
+  \method{is.close}{default}(X,r, Y=NULL, \dots, periodic=FALSE)
 
-  \method{is.close}{ppp}(X,r, Y=NULL)
+  \method{is.close}{ppp}(X,r, Y=NULL, \dots, periodic=FALSE, sorted=FALSE)
 
-  \method{is.close}{pp3}(X,r, Y=NULL)
+  \method{is.close}{pp3}(X,r, Y=NULL, \dots, periodic=FALSE, sorted=FALSE)
 }
 \arguments{
   \item{X,Y}{
@@ -26,6 +26,17 @@
   \item{r}{
     Threshold distance: a number greater than zero.
   }
+  \item{periodic}{
+    Logical value indicating whether to measure distances in the
+    periodic sense, so that opposite sides of the (rectangular) window
+    are treated as identical.
+  }
+  \item{sorted}{
+    Logical value, indicating whether the points of \code{X}
+    (and \code{Y}, if given) are already sorted into increasing order of the
+    \eqn{x} coordinates.
+  }
+  \item{\dots}{Other arguments are ignored.}
 }
 \details{
   This is simply a faster version of \code{(nndist(X) <= r)}

diff --git a/src/isclose.c b/src/isclose.c
@@ -2,7 +2,7 @@
 
   isclose.c
 
-  $Revision: 1.1 $ $Date: 2016/10/14 01:22:03 $
+  $Revision: 1.2 $ $Date: 2016/10/15 06:50:24 $
 
   Determine whether a point has a neighbour closer than 'r'
   
@@ -11,17 +11,38 @@
 
 #include <R.h>
 
+#undef TORUS
+
+#undef ZCOORD
+
 #define CLOSEFUN isXclose
 #define CROSSFUN isXYclose
-#undef ZCOORD
 #include "isclose.h"
 #undef CLOSEFUN
 #undef CROSSFUN
 
+#define ZCOORD
+
 #define CLOSEFUN isX3close
 #define CROSSFUN isXY3close
-#define ZCOORD
 #include "isclose.h"
 #undef CLOSEFUN
 #undef CROSSFUN
+
+#define TORUS
+
 #undef ZCOORD
+
+#define CLOSEFUN isXpclose
+#define CROSSFUN isXYpclose
+#include "isclose.h"
+#undef CLOSEFUN
+#undef CROSSFUN
+
+#define ZCOORD
+
+#define CLOSEFUN isX3pclose
+#define CROSSFUN isXY3pclose
+#include "isclose.h"
+#undef CLOSEFUN
+#undef CROSSFUN