diff --git a/DESCRIPTION b/DESCRIPTION
index 7bd38877..43a88d23 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -40,7 +40,8 @@ Imports:
     checkmate,
     stringr,
     pbapply,
-    zip
+    zip,
+    matrixStats
 Suggests:
     spelling,
     Rvcg (>= 0.17),
diff --git a/R/coordinates.R b/R/coordinates.R
index b5f74d2b..8dbbe51e 100644
--- a/R/coordinates.R
+++ b/R/coordinates.R
@@ -82,19 +82,21 @@ 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 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 
-#' 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
@@ -104,22 +106,22 @@ 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 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,
-                            Clamp=FALSE, CheckRanges=!Clamp, ...){
+                            Clamp=FALSE, CheckRanges=!Clamp, ...) {
   if(is.object(imdims)){
     if(!inherits(imdims, "im3d"))
       imdims=as.im3d(imdims)
@@ -137,20 +139,41 @@ 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)
-  
-  pixcoords=t(round(t(coords)/voxdims))+1
-  
+  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")
+  }
+  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, 
+                               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)
+  } 
+  # 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]))
     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)
+    ranges=t(matrixStats::colRanges(pixcoords))
     if(any(ranges[2,]>imdims) || any(ranges[1,]<1))
       stop("pixcoords out of range")
   }
@@ -168,13 +191,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.2'))
+    return(natcpp::c_sub2ind(dims, indices))
+  
   k=cumprod(c(1,dims[-length(dims)]))
   ndx=1
   for(i in 1:length(dims)){
diff --git a/R/im3d.R b/R/im3d.R
index 44b0b674..9ca61ec7 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.2')) {
+    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/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
 }
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
diff --git a/man/coord2ind.Rd b/man/coord2ind.Rd
index 4139a27d..6a3e641f 100644
--- a/man/coord2ind.Rd
+++ b/man/coord2ind.Rd
@@ -20,24 +20,26 @@ 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.}
 
-\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{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.}
 }
@@ -46,15 +48,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
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
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)
+}
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)