Accelerated code for distance from points to segments

Author: baddstats

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: spatstat.utils
-Version: 3.2-2.003
-Date: 2026-04-29
+Version: 3.2-2.004
+Date: 2026-05-02
 Title: Utility Functions for 'spatstat'
 Authors@R: c(person("Adrian", "Baddeley", 
                     role = c("aut", "cre"),

NEWS

@@ -1,5 +1,5 @@
 
-        CHANGES IN spatstat.utils VERSION 3.2-2.003
+        CHANGES IN spatstat.utils VERSION 3.2-2.004
 
 OVERVIEW
 

R/xysegment.R

@@ -1,7 +1,7 @@
 #
 #      xysegment.S
 #
-#     $Revision: 1.22 $    $Date: 2022/04/28 03:25:59 $
+#     $Revision: 1.23 $    $Date: 2026/05/02 06:04:37 $
 #
 # Low level utilities for analytic geometry for line segments
 #
@@ -25,8 +25,9 @@
 #
 # NNdist2segs
 #       distance to nearest line segment, from each point in a list.
-
-
+#
+# Copyright (c) Adrian Baddeley and Rob Foxall 1997-2026
+# GNU Public Licence (>= 2.0)
 
 distpl <- function(p, l) {
   xp <- p[1]

inst/doc/packagesizes.txt

@@ -33,4 +33,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
 "2025-09-20" "3.2-0" 42 195 0 3716 2462
 "2026-01-08" "3.2-1" 43 196 0 3754 2535
 "2026-03-10" "3.2-2" 44 197 0 3778 2535
-"2026-04-29" "3.2-2.003" 44 197 0 3784 2542
+"2026-05-02" "3.2-2.004" 44 197 0 3785 2583

inst/info/packagesizes.txt

@@ -33,4 +33,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
 "2025-09-20" "3.2-0" 42 195 0 3716 2462
 "2026-01-08" "3.2-1" 43 196 0 3754 2535
 "2026-03-10" "3.2-2" 44 197 0 3778 2535
-"2026-04-29" "3.2-2.003" 44 197 0 3784 2542
+"2026-05-02" "3.2-2.004" 44 197 0 3785 2583

src/distseg.h

@@ -13,18 +13,18 @@
 
   Author: Adrian Baddeley 2018-2021
 
-  Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner 2001-2021
+  Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner 2001-2026
   Licence: GNU Public Licence >= 2
 
-  $Revision: 1.3 $ $Date: 2022/10/19 02:58:46 $
+  $Revision: 1.5 $ $Date: 2026/05/02 07:03:14 $
 
 */
 
 void
 FNAME(
      /* input */
      double	*xp,
-     double     *yp,		        /* point/pixel coordinates */
+     double     *yp,		        /* point (or pixel) coordinates */
      int	*npoints,
      double	*x0,
      double     *y0,
@@ -33,7 +33,7 @@ FNAME(
      int	*nsegments,
      double     *epsilon,               /* tolerance for short segments */
      /* output */
-     double	*dist2		        /* squared distance from pixel 
+     double	*dist2		        /* squared distance from point
 					   to nearest line segment 
 					   INITIALISED TO LARGE VALUE */
 #ifdef WANT_INDEX
@@ -42,9 +42,15 @@ FNAME(
 ) {
   int	i,j, np, nseg, maxchunk;
   double dx,dy,leng,co,si;  /* parameters of segment */
+  double xpi, ypi, x0j, y0j, x1j, y1j; /* coordinates */
   double xdif0,ydif0,xdif1,ydif1,xpr,ypr; /* vectors */
   double dsq0,dsq1,dsq,dsqperp; /* squared distances */
   double eps;
+  double *xpp, *ypp;        /* pointers into data xp[], yp[] */
+  double *dist2p;           /* pointer into dist2[] */
+#ifdef WANT_INDEX
+  int    *indexp;           /* pointer into index[] */
+#endif
 
   np   = *npoints;
   nseg = *nsegments;
@@ -53,55 +59,90 @@ FNAME(
   OUTERCHUNKLOOP(j, nseg, maxchunk, 16384) {
     R_CheckUserInterrupt();
     INNERCHUNKLOOP(j, nseg, maxchunk, 16384) {
-      dx = x1[j] - x0[j];
-      dy = y1[j] - y0[j];
+      /* segment j coordinates */
+      x0j = x0[j];
+      y0j = y0[j];
+      x1j = x1[j];
+      y1j = y1[j];
+      /* segment j geometry */
+      dx = x1j - x0j;
+      dy = y1j - y0j;
       leng = hypot(dx, dy);
+      /* loop over points */
       if(leng > eps) {
 	/* normal case */
 	co = dx/leng;
 	si = dy/leng;
-	for(i = 0; i < np; i++) {
-	  /* vectors from pixel to segment endpoints */
-	  xdif0 =  xp[i] - x0[j];
-	  ydif0 =  yp[i] - y0[j];
-	  xdif1 =  xp[i] - x1[j];
-	  ydif1 =  yp[i] - y1[j];
+	for(
+#ifdef WANT_INDEX
+	    i = np, xpp = xp, ypp = yp, dist2p = dist2, indexp = index;
+	    i > 0;
+	    i--,    xpp++,    ypp++,    dist2p++,       indexp++
+#else	    
+	    i = np, xpp = xp, ypp = yp, dist2p = dist2;
+	    i > 0;
+	    i--,    xpp++,    ypp++,    dist2p++
+#endif	    
+	    ) {
+	  /* point i coordinates */
+	  xpi = *xpp;
+	  ypi = *ypp;
+	  /* vectors from point i to endpoints of segment j */
+	  xdif0 =  xpi - x0j;
+	  ydif0 =  ypi - y0j;
+	  xdif1 =  xpi - x1j;
+	  ydif1 =  ypi - y1j;
 	  /* squared distances to segment endpoints */
 	  dsq0 = xdif0 * xdif0 + ydif0 * ydif0;
 	  dsq1 = xdif1 * xdif1 + ydif1 * ydif1;
 	  dsq = (dsq0 < dsq1) ? dsq0 : dsq1;
 	  /* rotate pixel around 1st endpoint of segment
 	     so that line segment lies in x axis */
 	  xpr = xdif0 * co + ydif0 * si;
-	  ypr = -xdif0 * si + ydif0 * co;
 	  /* perpendicular distance applies only in perpendicular region */
 	  if(xpr >= 0.0 && xpr <= leng) {
+	    ypr = -xdif0 * si + ydif0 * co;
 	    dsqperp = ypr * ypr;
 	    if(dsqperp < dsq) dsq = dsqperp;
 	  }
-	  if(dist2[i] > dsq) {
-	    dist2[i] = dsq;
+	  /* update minimum distance for pixel i */
+	  if(*dist2p > dsq) {
+	    *dist2p = dsq;
 #ifdef WANT_INDEX	    
-	    index[i] = j;
+	    *indexp = j;
 #endif	    
 	  }
 	}
       } else {
 	/* short segment - use endpoints only */
-	for(i = 0; i < np; i++) {
-	  /* vectors from pixel to segment endpoints */
-	  xdif0 =  xp[i] - x0[j];
-	  ydif0 =  yp[i] - y0[j];
-	  xdif1 =  xp[i] - x1[j];
-	  ydif1 =  yp[i] - y1[j];
+	for(
+#ifdef WANT_INDEX
+	    i = np, xpp = xp, ypp = yp, dist2p = dist2, indexp = index;
+	    i > 0;
+	    i--,    xpp++,    ypp++,    dist2p++,       indexp++
+#else	    
+	    i = np, xpp = xp, ypp = yp, dist2p = dist2;
+	    i > 0;
+	    i--,    xpp++,    ypp++,    dist2p++
+#endif	    
+	    ) {
+	  /* point i coordinates */
+	  xpi = *xpp;
+	  ypi = *ypp;
+	  /* vectors from point i to endpoints of segment j */
+	  xdif0 =  xpi - x0j;
+	  ydif0 =  ypi - y0j;
+	  xdif1 =  xpi - x1j;
+	  ydif1 =  ypi - y1j;
 	  /* squared distances to segment endpoints */
 	  dsq0 = xdif0 * xdif0 + ydif0 * ydif0;
 	  dsq1 = xdif1 * xdif1 + ydif1 * ydif1;
 	  dsq = (dsq0 < dsq1) ? dsq0 : dsq1;
-	  if(dist2[i] > dsq) {
-	    dist2[i] = dsq;
+	  /* update minimum distance for pixel i */
+	  if(*dist2p > dsq) {
+	    *dist2p = dsq;
 #ifdef WANT_INDEX	    
-	    index[i] = j;
+	    *indexp = j;
 #endif	    
 	  }
 	}