From 2457a652d1858eeeb28a794089055a0b2ce5b29c Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sun, 12 Oct 2025 21:52:13 +0100 Subject: [PATCH 01/16] doc: whitespace changes --- R/coordinates.R | 12 ++++++------ man/coord2ind.Rd | 8 ++++---- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index b5f74d2b..99cafa40 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -82,19 +82,19 @@ ind2coord.im3d<-function(inds, voxdims=NULL, origin=NULL, ...){ #' Find 1D or 3D voxel indices into a 3D image given spatial coordinates -#' +#' #' @details \code{coord2ind} is designed to cope with any user-defined class for -#' which an as.im3d method exists. Presently the only example in the nat.* -#' ecosystem is \code{nat.templatebrains::as.im3d.templatebrain}. The -#' existence of an \code{as.im3d} method implies that +#' which an as.im3d method exists. Presently the only example in the nat.* +#' ecosystem is \code{nat.templatebrains::as.im3d.templatebrain}. The +#' existence of an \code{as.im3d} method implies that #' \code{voxdims},\code{origin}, and \code{dim} functions can be called. This #' is the necessary information required to convert i,j,k logical indices into #' x,y,z spatial indices. #' @param coords spatial coordinates of image voxels. #' @param ... extra arguments passed to methods. #' @export -#' @examples -#' coord2ind(cbind(1,2,3), imdims = c(1024,512,218), +#' @examples +#' coord2ind(cbind(1,2,3), imdims = c(1024,512,218), #' voxdims = c(0.622088, 0.622088, 0.622088), origin = c(0,0,0)) #' \dontrun{ #' ## repeat but using a templatebrain object to specify the coordinate system diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index 4139a27d..1900343f 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -46,15 +46,15 @@ Find 1D or 3D voxel indices into a 3D image given spatial coordinates } \details{ \code{coord2ind} is designed to cope with any user-defined class for - which an as.im3d method exists. Presently the only example in the nat.* - ecosystem is \code{nat.templatebrains::as.im3d.templatebrain}. The - existence of an \code{as.im3d} method implies that + which an as.im3d method exists. Presently the only example in the nat.* + ecosystem is \code{nat.templatebrains::as.im3d.templatebrain}. The + existence of an \code{as.im3d} method implies that \code{voxdims},\code{origin}, and \code{dim} functions can be called. This is the necessary information required to convert i,j,k logical indices into x,y,z spatial indices. } \examples{ -coord2ind(cbind(1,2,3), imdims = c(1024,512,218), +coord2ind(cbind(1,2,3), imdims = c(1024,512,218), voxdims = c(0.622088, 0.622088, 0.622088), origin = c(0,0,0)) \dontrun{ ## repeat but using a templatebrain object to specify the coordinate system From 9e83c3ae8cb53a9021c5ceb30935129d877b5150 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sun, 12 Oct 2025 23:14:41 +0100 Subject: [PATCH 02/16] clarify coords doesn't have to be a matrix * and doesn't have to be coerced to one either --- R/coordinates.R | 15 ++++++++++----- man/coord2ind.Rd | 4 +++- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index 99cafa40..e8e98f1a 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -90,7 +90,9 @@ ind2coord.im3d<-function(inds, voxdims=NULL, origin=NULL, ...){ #' \code{voxdims},\code{origin}, and \code{dim} functions can be called. This #' is the necessary information required to convert i,j,k logical indices into #' x,y,z spatial indices. -#' @param coords spatial coordinates of image voxels. +#' @param coords spatial coordinates of image voxels. Must be an Nx3 matrix or +#' data.frame or an object for which \code{ncol} and column subscripting +#' `[,1]` work. #' @param ... extra arguments passed to methods. #' @export #' @examples @@ -137,10 +139,13 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, if(length(imdims) != 3) stop('coord2ind only handles 3D data') - if(!is.matrix(coords)) - coords=matrix(coords,byrow=TRUE,ncol=length(coords)) - if(!missing(origin)) - coords=t(t(coords)-origin) + if(is.null(dim(coords))) { + if(length(coords)==3) coords=xyzmatrix(coords) + else stop("coordinates should be an N x 3 matrix") + } else { + if(ncol(coords)!=3) + stop("coordinates should be an N x 3 matrix-like object") + } pixcoords=t(round(t(coords)/voxdims))+1 diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index 1900343f..5ac0c3d7 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -20,7 +20,9 @@ coord2ind(coords, ...) ) } \arguments{ -\item{coords}{spatial coordinates of image voxels.} +\item{coords}{spatial coordinates of image voxels. Must be an Nx3 matrix or +data.frame or an object for which \code{ncol} and column subscripting +`[,1]` work.} \item{...}{extra arguments passed to methods.} From 4f9d69889f723d8d1f5e947d9be6f72348cee076 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Sun, 12 Oct 2025 23:24:43 +0100 Subject: [PATCH 03/16] coord2ind: alternative method with reduced memory use MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * for testing, provide arg to switch between new and original algorithm > bench::mark(coord2ind(m, testImage), coord2ind(m, testImage, use.scale = T), min_time = 3) # A tibble: 2 × 13 expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc 1 coord2ind(m… 63.7ms 73.4ms 1.16106306906514e1 252MB 2.25577967704085e1 35 68 2 coord2ind(m… 56.3ms 62.8ms 1.39956214910579e1 206MB 2.13266613197073e1 42 64 --- R/coordinates.R | 13 +++++++++++-- man/coord2ind.Rd | 3 +++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index e8e98f1a..fe645b68 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -112,6 +112,7 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @param origin the origin of the 3D image. #' @param linear.indices Whether or not to convert the voxel indices into a #' linear 1D form (the default) or to keep as 3D indices. +#' @param use.scale For testing an alternative algorithm #' @param aperm permutation order for axes. #' @param Clamp Whether or not to map out of range coordinates to the nearest #' in range index (default \code{FALSE}) @@ -120,7 +121,7 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @export #' @rdname coord2ind coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, - linear.indices=TRUE, aperm=NULL, + linear.indices=TRUE, aperm=NULL, use.scale=FALSE, Clamp=FALSE, CheckRanges=!Clamp, ...){ if(is.object(imdims)){ if(!inherits(imdims, "im3d")) @@ -147,7 +148,15 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, stop("coordinates should be an N x 3 matrix-like object") } - pixcoords=t(round(t(coords)/voxdims))+1 + if(use.scale) { + if(missing(origin) || is.null(origin) || all(origin==0)) + origin=FALSE + pixcoords=round(scale(coords, center = origin, scale = voxdims))+1L + } else { + if(!missing(origin)) + coords=t(t(coords)-origin) + pixcoords=t(round(t(coords)/voxdims))+1 + } # make sure no points are out of range if(Clamp){ diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index 5ac0c3d7..4fdeec74 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -14,6 +14,7 @@ coord2ind(coords, ...) origin = NULL, linear.indices = TRUE, aperm = NULL, + use.scale = FALSE, Clamp = FALSE, CheckRanges = !Clamp, ... @@ -38,6 +39,8 @@ linear 1D form (the default) or to keep as 3D indices.} \item{aperm}{permutation order for axes.} +\item{use.scale}{For testing an alternative algorithm} + \item{Clamp}{Whether or not to map out of range coordinates to the nearest in range index (default \code{FALSE})} From 879a185acb9b04dd61ff19de3e5ad6fda0001bbb Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Mon, 13 Oct 2025 09:56:53 +0100 Subject: [PATCH 04/16] use matrixStats::colRanges to save more memory MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit > bench::mark(v1=coord2ind(mdf, testImage, version=1), v2=coord2ind(mdf, testImage, version=2), v3=coord2ind(mdf, testImage, version=3), min_time = 1) # A tibble: 3 × 13 expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time 1 v1 82.9ms 90.2ms 7.93571545940208e0 298MB 1.78553597836547e1 8 18 1.01s 2 v2 68.6ms 80ms 9.53099788536415e0 252MB 1.90619957707283e1 10 20 1.05s 3 v3 48.6ms 60.1ms 1.44604097192445e1 172MB 1.9280546292326 e1 15 20 1.04s # ℹ 4 more variables: result , memory , time , gc --- DESCRIPTION | 3 ++- R/coordinates.R | 12 ++++++++---- man/coord2ind.Rd | 4 ++-- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7bd38877..7e1f33e7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -56,7 +56,8 @@ Suggests: readobj, brotli, qs, - jsonlite + jsonlite, + matrixStats Enhances: natcpp License: GPL-3 diff --git a/R/coordinates.R b/R/coordinates.R index fe645b68..ae168b08 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -112,7 +112,7 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @param origin the origin of the 3D image. #' @param linear.indices Whether or not to convert the voxel indices into a #' linear 1D form (the default) or to keep as 3D indices. -#' @param use.scale For testing an alternative algorithm +#' @param version For testing alternative algorithms #' @param aperm permutation order for axes. #' @param Clamp Whether or not to map out of range coordinates to the nearest #' in range index (default \code{FALSE}) @@ -121,7 +121,7 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @export #' @rdname coord2ind coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, - linear.indices=TRUE, aperm=NULL, use.scale=FALSE, + linear.indices=TRUE, aperm=NULL, version=1L, Clamp=FALSE, CheckRanges=!Clamp, ...){ if(is.object(imdims)){ if(!inherits(imdims, "im3d")) @@ -148,7 +148,7 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, stop("coordinates should be an N x 3 matrix-like object") } - if(use.scale) { + if(version>=2) { if(missing(origin) || is.null(origin) || all(origin==0)) origin=FALSE pixcoords=round(scale(coords, center = origin, scale = voxdims))+1L @@ -164,7 +164,11 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, pixcoords[,2]=pmin(imdims[2],pmax(1,pixcoords[,2])) pixcoords[,3]=pmin(imdims[3],pmax(1,pixcoords[,3])) } else if(CheckRanges){ - ranges=apply(pixcoords,2,range) + if(version>=3 && requireNamespace('matrixStats', quietly = T)) + ranges=t(matrixStats::colRanges(pixcoords)) + else { + ranges=apply(pixcoords,2,range) + } if(any(ranges[2,]>imdims) || any(ranges[1,]<1)) stop("pixcoords out of range") } diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index 4fdeec74..c89b7167 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -14,7 +14,7 @@ coord2ind(coords, ...) origin = NULL, linear.indices = TRUE, aperm = NULL, - use.scale = FALSE, + version = 1L, Clamp = FALSE, CheckRanges = !Clamp, ... @@ -39,7 +39,7 @@ linear 1D form (the default) or to keep as 3D indices.} \item{aperm}{permutation order for axes.} -\item{use.scale}{For testing an alternative algorithm} +\item{version}{For testing alternative algorithms} \item{Clamp}{Whether or not to map out of range coordinates to the nearest in range index (default \code{FALSE})} From 232ca9baa39d02cc6ede3e2b5b3153409ec37b6a Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Mon, 13 Oct 2025 14:28:53 +0100 Subject: [PATCH 05/16] More use of matrixStats to reduce memory footprint * empirically t_tx_OP_y seems to be better than scale (benefit is more limited with data.frame input mind you) --- R/coordinates.R | 32 +++++++++++++++++++++----------- man/coord2ind.Rd | 13 +++++++------ 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index ae168b08..5a72e5fc 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -106,23 +106,25 @@ ind2coord.im3d<-function(inds, voxdims=NULL, origin=NULL, ...){ coord2ind <- function(coords, ...) UseMethod("coord2ind") -#' @param imdims array dimensions of 3D image \emph{OR} an object for which a +#' @param imdims array dimensions of 3D image \emph{OR} an object for which a #' \code{\link{as.im3d}} object has been defined (see Details). #' @param voxdims vector of 3 voxels dimensions (width, height, depth). #' @param origin the origin of the 3D image. -#' @param linear.indices Whether or not to convert the voxel indices into a +#' @param linear.indices Whether or not to convert the voxel indices into a #' linear 1D form (the default) or to keep as 3D indices. -#' @param version For testing alternative algorithms +#' @param version For testing alternative algorithms. Currently 1-4 where 1 is +#' the original version and 4 seems to be the best so far. #' @param aperm permutation order for axes. -#' @param Clamp Whether or not to map out of range coordinates to the nearest -#' in range index (default \code{FALSE}) +#' @param Clamp Whether or not to map out of range coordinates to the nearest in +#' range index (default \code{FALSE}) #' @param CheckRanges whether to check if coordinates are out of range. #' @seealso \code{\link{ind2coord}}, \code{\link{sub2ind}}, \code{\link{ijkpos}} #' @export #' @rdname coord2ind coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, - linear.indices=TRUE, aperm=NULL, version=1L, - Clamp=FALSE, CheckRanges=!Clamp, ...){ + linear.indices=TRUE, aperm=NULL, version=4L, + Clamp=FALSE, CheckRanges=!Clamp, ...) { + checkmate::assertIntegerish(version, lower = 1L, upper = 4L) if(is.object(imdims)){ if(!inherits(imdims, "im3d")) imdims=as.im3d(imdims) @@ -148,10 +150,18 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, stop("coordinates should be an N x 3 matrix-like object") } - if(version>=2) { - if(missing(origin) || is.null(origin) || all(origin==0)) - origin=FALSE - pixcoords=round(scale(coords, center = origin, scale = voxdims))+1L + if(version>=4) { + if(missing(origin) || is.null(origin)) origin=c(0,0,0) + origin=origin-voxdims + coords=matrixStats::t_tx_OP_y(as.matrix(coords), origin, OP = '-') + coords=matrixStats::t_tx_OP_y(coords, voxdims, OP = '/') + pixcoords=round(coords) + } else if(version>=2) { + if(missing(origin) || is.null(origin)) origin=c(0,0,0) + origin=origin-voxdims + # if(missing(origin) || is.null(origin) || all(origin==0)) + # origin=FALSE + pixcoords=round(scale(coords, center = origin, scale = voxdims)) } else { if(!missing(origin)) coords=t(t(coords)-origin) diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index c89b7167..1d792d19 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -14,7 +14,7 @@ coord2ind(coords, ...) origin = NULL, linear.indices = TRUE, aperm = NULL, - version = 1L, + version = 4L, Clamp = FALSE, CheckRanges = !Clamp, ... @@ -27,22 +27,23 @@ data.frame or an object for which \code{ncol} and column subscripting \item{...}{extra arguments passed to methods.} -\item{imdims}{array dimensions of 3D image \emph{OR} an object for which a +\item{imdims}{array dimensions of 3D image \emph{OR} an object for which a \code{\link{as.im3d}} object has been defined (see Details).} \item{voxdims}{vector of 3 voxels dimensions (width, height, depth).} \item{origin}{the origin of the 3D image.} -\item{linear.indices}{Whether or not to convert the voxel indices into a +\item{linear.indices}{Whether or not to convert the voxel indices into a linear 1D form (the default) or to keep as 3D indices.} \item{aperm}{permutation order for axes.} -\item{version}{For testing alternative algorithms} +\item{version}{For testing alternative algorithms. Currently 1-4 where 1 is +the original version and 4 seems to be the best so far.} -\item{Clamp}{Whether or not to map out of range coordinates to the nearest -in range index (default \code{FALSE})} +\item{Clamp}{Whether or not to map out of range coordinates to the nearest in +range index (default \code{FALSE})} \item{CheckRanges}{whether to check if coordinates are out of range.} } From 45aed2a1a8809138262b96c4cf3a83efb17ced82 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 10:37:50 +0100 Subject: [PATCH 06/16] Another 4 versions of testing coord2ind using natcpp * the final version is ~ 40x more memory efficient and 10x faster than the original implementation --- R/coordinates.R | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index 5a72e5fc..0d50bd5b 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -122,9 +122,9 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @export #' @rdname coord2ind coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, - linear.indices=TRUE, aperm=NULL, version=4L, + linear.indices=TRUE, aperm=NULL, version=7L, Clamp=FALSE, CheckRanges=!Clamp, ...) { - checkmate::assertIntegerish(version, lower = 1L, upper = 4L) + checkmate::assertIntegerish(version, lower = 1L, upper = 7L) if(is.object(imdims)){ if(!inherits(imdims, "im3d")) imdims=as.im3d(imdims) @@ -149,8 +149,26 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, if(ncol(coords)!=3) stop("coordinates should be an N x 3 matrix-like object") } - - if(version>=4) { + if(version>=7 && use_natcpp(version='0.1.1') && linear.indices && is.null(aperm)) { + if(missing(origin) || is.null(origin)) origin=c(0,0,0) + res=natcpp::c_coords21dindex(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) + return(res) + } else if(version>=6 && use_natcpp(version='0.1.1')) { + if(missing(origin) || is.null(origin)) origin=c(0,0,0) + pixcoords=natcpp::c_ijkpos(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) + if (!is.null(aperm)) + imdims=imdims[aperm] + res=if(isTRUE(linear.indices)) natcpp::c_sub2ind(imdims, pixcoords) else pixcoords + return(res) + } else if(version>=5 && use_natcpp(version='0.1.1')) { + coords=as.matrix(coords) + if(missing(origin) || is.null(origin)) origin=c(0,0,0) + pixcoords=natcpp::c_ijkpos(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) + if (!is.null(aperm)) + imdims=imdims[aperm] + res=if(isTRUE(linear.indices)) natcpp::c_sub2ind(imdims, pixcoords) else pixcoords + return(res) + } else if(version>=4) { if(missing(origin) || is.null(origin)) origin=c(0,0,0) origin=origin-voxdims coords=matrixStats::t_tx_OP_y(as.matrix(coords), origin, OP = '-') From d80cd618b35c99842732f23057270db947965578 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 10:39:19 +0100 Subject: [PATCH 07/16] give use_natcpp a (min) version argument --- R/utils.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/utils.R b/R/utils.R index 246642fb..84231457 100644 --- a/R/utils.R +++ b/R/utils.R @@ -19,7 +19,8 @@ nat_progress <- function (x, max = 100, message = NULL) { # to check if we should use natcpp # always=TRUE => use if installed even if option says otherwise -use_natcpp <- function(always=FALSE) { +# version = minimum required version +use_natcpp <- function(always=FALSE, version='0') { opcheck <- isTRUE(always) || !isFALSE(getOption('nat.use_natcpp')) - opcheck && requireNamespace('natcpp', quietly = TRUE) + opcheck && requireNamespace('natcpp', quietly = TRUE) && packageVersion('natcpp')>=version } From 91e08e1d1c9b79bda3a2a4703124057a007106b5 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 11:28:15 +0100 Subject: [PATCH 08/16] roxygenise --- man/coord2ind.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index 1d792d19..61ca590d 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -14,7 +14,7 @@ coord2ind(coords, ...) origin = NULL, linear.indices = TRUE, aperm = NULL, - version = 4L, + version = 7L, Clamp = FALSE, CheckRanges = !Clamp, ... From 09275887a9ae35d6fa333b1c272b7222b0f4c1b9 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 11:28:37 +0100 Subject: [PATCH 09/16] implement faster version of natcpp --- R/im3d.R | 18 +++++++++++++++--- man/im3d-coords.Rd | 5 ++++- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/R/im3d.R b/R/im3d.R index 44b0b674..5899d009 100644 --- a/R/im3d.R +++ b/R/im3d.R @@ -1166,6 +1166,8 @@ xyzpos<-function(d, ijk) #' @param xyz Nx3 matrix of physical coordinates #' @param roundToNearestPixel Whether to round calculated pixel coordinates to #' nearest integer value (i.e. nearest pixel). default: \code{TRUE} +#' @param clamp Whether to clamp any pixel coordinates to the dimensions of the +#' image #' @return Nx3 matrix of physical or pixel coordinates #' @rdname im3d-coords #' @aliases ijkpos @@ -1175,7 +1177,7 @@ xyzpos<-function(d, ijk) #' d=im3d(,dim=c(20,30,40),origin=c(10,20,30),voxdims=c(1,2,3)) #' # check round trip for origin #' stopifnot(all.equal(ijkpos(d,xyzpos(d,c(1,1,1))), c(1,1,1))) -ijkpos<-function(d, xyz, roundToNearestPixel=TRUE) +ijkpos<-function(d, xyz, roundToNearestPixel=TRUE, clamp=TRUE) { # return the ijk position for a physical location (x,y,z) # This will be the pixel centre based on the bounding box @@ -1183,9 +1185,19 @@ ijkpos<-function(d, xyz, roundToNearestPixel=TRUE) # transpose if we have received a matrix (with 3 cols x,y,z) so that # multiplication below doesn not need to be changed + + od=origin(d) + checkmate::assert_numeric(od, len = 3) + vd=voxdims(d) + checkmate::assert_numeric(vd, len = 3) + dims=dim(d) + checkmate::assert_integer(dims, len = 3) + if(!is.null(dim(xyz)) && roundToNearestPixel && use_natcpp(version = '0.1.1.9000')) { + res=natcpp::c_ijkpos(xyz, dims = dims, origin = od, voxdims = vd, clamp = clamp) + return(res) + } if(is.matrix(xyz)) xyz=t(xyz) - - ijk=(xyz-origin(d))/voxdims(d)+1 + ijk=(xyz-od)/voxdims(d)+1 if(roundToNearestPixel) { ijk=round(ijk) if(any(ijk<1) || any(ijk>dim(d))) warning("pixel coordinates outside image data") diff --git a/man/im3d-coords.Rd b/man/im3d-coords.Rd index 82a7ec10..00f7e106 100644 --- a/man/im3d-coords.Rd +++ b/man/im3d-coords.Rd @@ -8,7 +8,7 @@ \usage{ xyzpos(d, ijk) -ijkpos(d, xyz, roundToNearestPixel = TRUE) +ijkpos(d, xyz, roundToNearestPixel = TRUE, clamp = TRUE) } \arguments{ \item{d}{An \code{im3d} object defining a physical space} @@ -19,6 +19,9 @@ ijkpos(d, xyz, roundToNearestPixel = TRUE) \item{roundToNearestPixel}{Whether to round calculated pixel coordinates to nearest integer value (i.e. nearest pixel). default: \code{TRUE}} + +\item{clamp}{Whether to clamp any pixel coordinates to the dimensions of the +image} } \value{ Nx3 matrix of physical or pixel coordinates From 915c8b932f691a66a01a797f66f1fb8db4b26c36 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 11:39:24 +0100 Subject: [PATCH 10/16] dev: use test_path for paths * just a convenience for running tests interactively --- tests/testthat/test-coordinates.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-coordinates.R b/tests/testthat/test-coordinates.R index 7299b4ef..e1cbef02 100644 --- a/tests/testthat/test-coordinates.R +++ b/tests/testthat/test-coordinates.R @@ -9,7 +9,7 @@ test_that("sub2ind returns correct indices", { }) test_that("ind2coord returns correct coordinates when given an im3d", { - testImage <- read.im3d("testdata/nrrd/LHMask.nrrd") + testImage <- read.im3d(test_path("testdata/nrrd/LHMask.nrrd")) coord <- ind2coord(testImage) coord <- coord[1:10] coord.expected <- c(40.6, 42, 36.4, 37.8, 39.2, 40.6, 42, 36.4, 37.8, 39.2) @@ -17,7 +17,7 @@ test_that("ind2coord returns correct coordinates when given an im3d", { }) test_that("coord2ind returns correct coordinates", { - testImage <- read.im3d("testdata/nrrd/LHMask.nrrd", ReadData = F) + testImage <- read.im3d(test_path("testdata/nrrd/LHMask.nrrd"), ReadData = F) ind <- coord2ind(matrix(c(10, 20, 30, 11, 20, 30), nrow=2, byrow=TRUE), testImage) ind.expected <- c(53208, 53209) expect_equal(ind, ind.expected) From a36fb8db9078eee9edcaa8b346e240ff935fc4e6 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 11:42:13 +0100 Subject: [PATCH 11/16] remove redundant implementations --- R/coordinates.R | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index 0d50bd5b..b9f6fafb 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -149,19 +149,11 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, if(ncol(coords)!=3) stop("coordinates should be an N x 3 matrix-like object") } - if(version>=7 && use_natcpp(version='0.1.1') && linear.indices && is.null(aperm)) { + if(version>=7 && use_natcpp(version='0.1.1.9000') && linear.indices && is.null(aperm)) { if(missing(origin) || is.null(origin)) origin=c(0,0,0) res=natcpp::c_coords21dindex(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) return(res) - } else if(version>=6 && use_natcpp(version='0.1.1')) { - if(missing(origin) || is.null(origin)) origin=c(0,0,0) - pixcoords=natcpp::c_ijkpos(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) - if (!is.null(aperm)) - imdims=imdims[aperm] - res=if(isTRUE(linear.indices)) natcpp::c_sub2ind(imdims, pixcoords) else pixcoords - return(res) - } else if(version>=5 && use_natcpp(version='0.1.1')) { - coords=as.matrix(coords) + } else if(version>=6 && use_natcpp(version='0.1.1.9000')) { if(missing(origin) || is.null(origin)) origin=c(0,0,0) pixcoords=natcpp::c_ijkpos(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) if (!is.null(aperm)) @@ -174,12 +166,6 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, coords=matrixStats::t_tx_OP_y(as.matrix(coords), origin, OP = '-') coords=matrixStats::t_tx_OP_y(coords, voxdims, OP = '/') pixcoords=round(coords) - } else if(version>=2) { - if(missing(origin) || is.null(origin)) origin=c(0,0,0) - origin=origin-voxdims - # if(missing(origin) || is.null(origin) || all(origin==0)) - # origin=FALSE - pixcoords=round(scale(coords, center = origin, scale = voxdims)) } else { if(!missing(origin)) coords=t(t(coords)-origin) From 01b44bf8bb8d692c62c41411a5d19eefe318c90f Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 11:53:47 +0100 Subject: [PATCH 12/16] try to work around CMAKE version issue CMake Error at CMakeLists.txt:38 (CMAKE_MINIMUM_REQUIRED): Compatibility with CMake < 3.5 has been removed from CMake. Update the VERSION argument value. Or, use the ... syntax to tell CMake that the project requires at least but has been updated to work with policies introduced by or earlier. Or, add -DCMAKE_POLICY_VERSION_MINIMUM=3.5 to try configuring anyway. --- build-cmtk.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build-cmtk.sh b/build-cmtk.sh index 2e1e158c..3f31618c 100755 --- a/build-cmtk.sh +++ b/build-cmtk.sh @@ -5,7 +5,7 @@ if [ ! -d "$HOME/usr/local/bin" ]; then mkdir -p $HOME/src && cd $HOME/src git clone --depth 10 --branch natdev https://github.com/jefferis/cmtk cd cmtk && git checkout natdev - cd core && mkdir build && cd build && cmake .. && make DESTDIR=$HOME/ all install + cd core && mkdir build && cd build && cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 .. && make DESTDIR=$HOME/ all install else echo 'Using cached $HOME/usr/local directory.'; fi From 4af8e5dc83498122d0f31b83e2fccd830c4dbae4 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 15:38:43 +0100 Subject: [PATCH 13/16] use natcpp for sub2ind * with example and some additional checks --- R/coordinates.R | 15 ++++++++++++--- man/sub2ind.Rd | 8 ++++++++ 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index b9f6fafb..e2c10475 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -200,13 +200,22 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, #' @param dims vector of dimensions of object to index into. #' @param indices vector of n-dimensional indices. #' @export -sub2ind<-function(dims,indices){ +#' @details +#' There is a *much* more efficient version implemented in natcpp > 0.1.1 +#' @examples +#' dims <- 3:5 +#' ijk <- matrix(c(1L, 1L, 1L, 3L, 4L, 5L), ncol = 3, byrow = TRUE) +#' sub2ind(dims, ijk) +sub2ind<-function(dims, indices){ # convert vector containing 1 coordinate into matrix if(!is.matrix(indices)) indices=matrix(indices,byrow=TRUE,ncol=length(indices)) - if(length(dims)!=ncol(indices)){ + if(length(dims)!=ncol(indices)) stop("indices must have the same number of columns as dimensions in dims") - } + dims <- as.integer(checkmate::assert_integerish(dims)) + if(use_natcpp(version = '0.1.1.9000')) + return(nat::c_sub2ind(dims, indices)) + k=cumprod(c(1,dims[-length(dims)])) ndx=1 for(i in 1:length(dims)){ diff --git a/man/sub2ind.Rd b/man/sub2ind.Rd index ea6a3ac0..d5397aa1 100644 --- a/man/sub2ind.Rd +++ b/man/sub2ind.Rd @@ -14,3 +14,11 @@ sub2ind(dims, indices) \description{ Emulates the MATLAB function \code{sub2ind}. } +\details{ +There is a *much* more efficient version implemented in natcpp > 0.1.1 +} +\examples{ +dims <- 3:5 +ijk <- matrix(c(1L, 1L, 1L, 3L, 4L, 5L), ncol = 3, byrow = TRUE) +sub2ind(dims, ijk) +} From a8cd4d49222fdabcca4338aca9ee1844ed183ce1 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 15:43:09 +0100 Subject: [PATCH 14/16] wire up the new natcpp code by default --- R/coordinates.R | 44 ++++++++++++++++++++------------------------ 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index e2c10475..c1f3e451 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -112,8 +112,6 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @param origin the origin of the 3D image. #' @param linear.indices Whether or not to convert the voxel indices into a #' linear 1D form (the default) or to keep as 3D indices. -#' @param version For testing alternative algorithms. Currently 1-4 where 1 is -#' the original version and 4 seems to be the best so far. #' @param aperm permutation order for axes. #' @param Clamp Whether or not to map out of range coordinates to the nearest in #' range index (default \code{FALSE}) @@ -122,7 +120,7 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") #' @export #' @rdname coord2ind coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, - linear.indices=TRUE, aperm=NULL, version=7L, + linear.indices=TRUE, aperm=NULL, Clamp=FALSE, CheckRanges=!Clamp, ...) { checkmate::assertIntegerish(version, lower = 1L, upper = 7L) if(is.object(imdims)){ @@ -149,29 +147,27 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, if(ncol(coords)!=3) stop("coordinates should be an N x 3 matrix-like object") } - if(version>=7 && use_natcpp(version='0.1.1.9000') && linear.indices && is.null(aperm)) { + if(use_natcpp(version='0.1.1.9000')) { if(missing(origin) || is.null(origin)) origin=c(0,0,0) - res=natcpp::c_coords21dindex(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) - return(res) - } else if(version>=6 && use_natcpp(version='0.1.1.9000')) { - if(missing(origin) || is.null(origin)) origin=c(0,0,0) - pixcoords=natcpp::c_ijkpos(coords, dims = imdims, origin = origin, voxdims = voxdims, clamp = Clamp) - if (!is.null(aperm)) - imdims=imdims[aperm] - res=if(isTRUE(linear.indices)) natcpp::c_sub2ind(imdims, pixcoords) else pixcoords + if(linear.indices && is.null(aperm)) + res=natcpp::c_coords21dindex(coords, dims = imdims, origin = origin, + voxdims = voxdims, clamp = Clamp) + else { + pixcoords=natcpp::c_ijkpos(coords, dims = imdims, origin = origin, + voxdims = voxdims, clamp = Clamp) + if (!is.null(aperm)) + imdims=imdims[aperm] + res=if(linear.indices) natcpp::c_sub2ind(imdims, pixcoords) else pixcoords + } return(res) - } else if(version>=4) { - if(missing(origin) || is.null(origin)) origin=c(0,0,0) - origin=origin-voxdims - coords=matrixStats::t_tx_OP_y(as.matrix(coords), origin, OP = '-') - coords=matrixStats::t_tx_OP_y(coords, voxdims, OP = '/') - pixcoords=round(coords) - } else { - if(!missing(origin)) - coords=t(t(coords)-origin) - pixcoords=t(round(t(coords)/voxdims))+1 - } - + } + # base R / favoured package implementation + if(missing(origin) || is.null(origin)) origin=c(0,0,0) + origin=origin-voxdims + coords=matrixStats::t_tx_OP_y(as.matrix(coords), origin, OP = '-') + coords=matrixStats::t_tx_OP_y(coords, voxdims, OP = '/') + pixcoords=round(coords) + # make sure no points are out of range if(Clamp){ pixcoords[,1]=pmin(imdims[1],pmax(1,pixcoords[,1])) From beae42625543ec87fdf0c8b60af26579aa77f779 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 15:55:25 +0100 Subject: [PATCH 15/16] eliminate thinko version checks * and add matrixStats to imports --- DESCRIPTION | 6 +++--- R/coordinates.R | 7 +------ man/coord2ind.Rd | 4 ---- 3 files changed, 4 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7e1f33e7..43a88d23 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -40,7 +40,8 @@ Imports: checkmate, stringr, pbapply, - zip + zip, + matrixStats Suggests: spelling, Rvcg (>= 0.17), @@ -56,8 +57,7 @@ Suggests: readobj, brotli, qs, - jsonlite, - matrixStats + jsonlite Enhances: natcpp License: GPL-3 diff --git a/R/coordinates.R b/R/coordinates.R index c1f3e451..f26c9a16 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -122,7 +122,6 @@ coord2ind <- function(coords, ...) UseMethod("coord2ind") coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, linear.indices=TRUE, aperm=NULL, Clamp=FALSE, CheckRanges=!Clamp, ...) { - checkmate::assertIntegerish(version, lower = 1L, upper = 7L) if(is.object(imdims)){ if(!inherits(imdims, "im3d")) imdims=as.im3d(imdims) @@ -174,11 +173,7 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, pixcoords[,2]=pmin(imdims[2],pmax(1,pixcoords[,2])) pixcoords[,3]=pmin(imdims[3],pmax(1,pixcoords[,3])) } else if(CheckRanges){ - if(version>=3 && requireNamespace('matrixStats', quietly = T)) - ranges=t(matrixStats::colRanges(pixcoords)) - else { - ranges=apply(pixcoords,2,range) - } + ranges=t(matrixStats::colRanges(pixcoords)) if(any(ranges[2,]>imdims) || any(ranges[1,]<1)) stop("pixcoords out of range") } diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd index 61ca590d..6a3e641f 100644 --- a/man/coord2ind.Rd +++ b/man/coord2ind.Rd @@ -14,7 +14,6 @@ coord2ind(coords, ...) origin = NULL, linear.indices = TRUE, aperm = NULL, - version = 7L, Clamp = FALSE, CheckRanges = !Clamp, ... @@ -39,9 +38,6 @@ linear 1D form (the default) or to keep as 3D indices.} \item{aperm}{permutation order for axes.} -\item{version}{For testing alternative algorithms. Currently 1-4 where 1 is -the original version and 4 seems to be the best so far.} - \item{Clamp}{Whether or not to map out of range coordinates to the nearest in range index (default \code{FALSE})} From 376cb13376e13808bb6f36dc721b832b47395ff1 Mon Sep 17 00:00:00 2001 From: Gregory Jefferis Date: Tue, 14 Oct 2025 22:48:56 +0100 Subject: [PATCH 16/16] fix namespace bug and update natcpp version * v0.2 is now merged to master --- R/coordinates.R | 6 +++--- R/im3d.R | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index f26c9a16..8dbbe51e 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -146,7 +146,7 @@ coord2ind.default<-function(coords, imdims, voxdims=NULL, origin=NULL, if(ncol(coords)!=3) stop("coordinates should be an N x 3 matrix-like object") } - if(use_natcpp(version='0.1.1.9000')) { + if(use_natcpp(version='0.2')) { if(missing(origin) || is.null(origin)) origin=c(0,0,0) if(linear.indices && is.null(aperm)) res=natcpp::c_coords21dindex(coords, dims = imdims, origin = origin, @@ -204,8 +204,8 @@ sub2ind<-function(dims, indices){ if(length(dims)!=ncol(indices)) stop("indices must have the same number of columns as dimensions in dims") dims <- as.integer(checkmate::assert_integerish(dims)) - if(use_natcpp(version = '0.1.1.9000')) - return(nat::c_sub2ind(dims, indices)) + if(use_natcpp(version = '0.2')) + return(natcpp::c_sub2ind(dims, indices)) k=cumprod(c(1,dims[-length(dims)])) ndx=1 diff --git a/R/im3d.R b/R/im3d.R index 5899d009..9ca61ec7 100644 --- a/R/im3d.R +++ b/R/im3d.R @@ -1192,7 +1192,7 @@ ijkpos<-function(d, xyz, roundToNearestPixel=TRUE, clamp=TRUE) checkmate::assert_numeric(vd, len = 3) dims=dim(d) checkmate::assert_integer(dims, len = 3) - if(!is.null(dim(xyz)) && roundToNearestPixel && use_natcpp(version = '0.1.1.9000')) { + if(!is.null(dim(xyz)) && roundToNearestPixel && use_natcpp(version = '0.2')) { res=natcpp::c_ijkpos(xyz, dims = dims, origin = od, voxdims = vd, clamp = clamp) return(res) }