diff --git a/DESCRIPTION b/DESCRIPTION index 5b7bbc7e..808b9909 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -63,6 +63,8 @@ Collate: 'dist3D_Segment_to_Segment.R' 'neuron.R' 'dotprops.R' + 'elastix.R' + 'elastix_io.R' 'geometry.R' 'graph-nodes.R' 'hxsurf.R' diff --git a/NAMESPACE b/NAMESPACE index 73fe8d86..a1fff03a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -183,6 +183,7 @@ S3method(xformimage,reglist) S3method(xformpoints,character) S3method(xformpoints,cmtkreg) S3method(xformpoints,default) +S3method(xformpoints,elastixreg) S3method(xformpoints,reglist) S3method(xformpoints,tpsreg) S3method(xyzmatrix,default) diff --git a/R/elastix.R b/R/elastix.R new file mode 100644 index 00000000..49ca239e --- /dev/null +++ b/R/elastix.R @@ -0,0 +1,233 @@ +elastixreg <- function(x, ...) { + if(!inherits(x,'elastixreg')) + class(x)=c("elastixreg",class(x)) + x +} + +#' @export +#' @rdname xformpoints +xformpoints.elastixreg<-function(reg, points, transformtype=c('warp','affine'), + direction=NULL, ...) { + if(is.list(reg)){ + stop("transformation using the in memory representation is not yet supported") + } + + # By default, or if swap=FALSE, we will use elastix's inverse direction + if(is.null(direction)) { + if(!is.null(swap<-attr(reg,'swap'))) { + direction=ifelse(swap, 'forward', 'inverse') + } else { + direction="inverse" + } + } else { + direction=match.arg(direction,c("inverse",'forward'), several.ok=TRUE) + } + + if (length(reg) > 1) { + # need to recycle manually + if (length(direction) == 1) + direction = rep(direction, length(reg)) + for (i in seq_along(reg)) { + points = xformpoints.elastixreg(reg[[i]], direction = direction[[i]], points = + points, ...) + } + return(points) + } + # check for NAs + nas=is.na(points[,1]) + if(sum(nas)) { + origpoints=points + points=points[!nas, , drop=FALSE] + } + rawout=transformix.points(points, reg) + pointst=xyzmatrix(rawout) + # if(transformtype=='warp'){ + # naPoints=is.na(pointst[,1]) + # if(any(naPoints)){ + # if(FallBackToAffine) { + # stop("FallBackToAffine not currently supported for elastix transforms!") + # affpoints = xformpoints(reg, points[naPoints,,drop=FALSE],transformtype='affine') + # pointst[naPoints, ] = affpoints + # } else { + # pointst[naPoints, ] = NA_real_ + # } + # } + # } + + # if(sum(nas)){ + # origpoints[!nas, ]=pointst + # origpoints + # } else { + # dimnames(pointst)=dimnames(points) + # pointst + # } + dimnames(pointst)=dimnames(points) + pointst +} + +# function to handle the details of calling the transformix command line tool +# nb reg must be a specific transform file (not a directory) +transformix.points <- function(points, reg, stderr=FALSE) { + pointsfile=tempfile(fileext=".txt") + on.exit(unlink(pointsfile)) + outdir=tempfile() + dir.create(outdir) + on.exit(unlink(outdir, recursive = TRUE), add = TRUE) + + writeLines(text = c("point", nvertices(points)), con = pointsfile) + write.table(points, file=pointsfile, row.names=FALSE, col.names=FALSE, append = TRUE) + + rval <- elastix.call('transformix', out=outdir, tp=reg, def=pointsfile) + + if(rval!=0) + stop("Error running transformix!", + "See log file in output directory for details:\n", outdir) + + outfile=dir(outdir, pattern = '\\.txt$', full.names = TRUE) + if(isFALSE(length(outfile)==1)) + stop("transformix error: could not file exactly one output file in ", + "the directory:", outdir) + transformixOut <- read.table(outfile, row.names=NULL) + rawout <- transformixOut[c(2,7:9,15:17,23:25,31:33, 39:41)] + colnames(rawout) <- c("id", "Ii", "Ij", "Ik", "Ix", "Iy","Iz", + "i", "j", "k", "x", "y", "z", "Dx", "Dy", "Dz") + rawout +} + +elastix.call <- function(tool=c("transformix", "elastix"), ..., + PROCESSED.ARGS=NULL, + FINAL.ARGS=NULL, stderr=FALSE, stdout=FALSE){ + tool=match.arg(tool) + cmd=file.path(elastix.bindir(check=TRUE), tool) + elastixargs=character() + + if(!is.null(PROCESSED.ARGS)){ + elastixargs=c(elastixargs, PROCESSED.ARGS) + } + + if(!missing(...)){ + xargs=pairlist(...) + for(n in names(xargs)){ + arg=xargs[[n]] + elastixarg=elastix.arg.name(n) + if(is.character(arg)){ + if(length(arg)!=1) stop("character arguments must have length 1") + elastixargs=c(elastixargs, elastixarg, arg) + } else if(is.logical(arg)){ + if(isTRUE(arg)) elastixargs=c(elastixargs, elastixarg) + } else if(is.numeric(arg)){ + arg=paste(arg, collapse=',') + elastixargs=c(elastixargs, elastixarg, arg) + } else if(is.null(arg)){ + # just ignore null arguemnts + } else { + stop("unrecognised argument type") + } + } + } + + if(!is.null(FINAL.ARGS)){ + elastixargs=c(elastixargs, FINAL.ARGS) + } + call=list(command=cmd, args=elastixargs) + cmtk.system2(call, stderr=stderr, stdout=stdout) +} + +elastix.arg.name<-function(x) { + # elastix sometimes uses -- but mostly uses - + prefix <- if(x %in% c("version", "extended.version", "help")) "--" else "-" + paste(prefix, gsub("\\.", '-', x), sep='') +} + +elastix.bindir<-function(firstdir=getOption('nat.elastix.bindir'), + extradirs=c('~/bin', '/usr/local/bin','/opt/local/bin'), + set=FALSE, check=FALSE, elastixtool='transformix'){ + # TODO check pure Windows vs cygwin + if(isTRUE(.Platform$OS.type=="windows" && tools::file_ext(elastixtool)!="exe")) + cmtktool=paste0(elastixtool,".exe") + bindir=NULL + if(!is.null(firstdir)) { + bindir=firstdir + if(check && !file.exists(file.path(bindir,elastixtool))) + stop("cmtk is _not_ installed at:", bindir, + "\nPlease check value of options('nat.elastix.bindir')") + } + # note the use of while loop + break to avoid heavy if nesting + while(is.null(bindir)){ + if(nzchar(cmtktoolpath<-Sys.which(elastixtool))){ + bindir=dirname(cmtktoolpath) + break + } + # otherwise check some plausible locations + for(d in extradirs){ + if(file.exists(file.path(d, elastixtool))) { + bindir=d + break + } + } + # we're out of luck but we still need to break out of the while loop + break + } + if(!is.null(bindir)){ + if(is.na(bindir)) bindir=NULL + else bindir=path.expand(bindir) + } + if(check && is.null(bindir)) + stop("Cannot find Elastix. Please install from ", + "http://elastix.isi.uu.nl and make sure that it is in your path!") + + if(set) { + options(nat.elastix.bindir=bindir) + } + bindir +} + +plot3d.elastixreg <- function(x, type=c("p","l","b"), ..., plotengine = getOption('nat.plotengine')) { + plotengine <- check_plotengine(plotengine) + if (plotengine == 'plotly') { + psh <- openplotlyscene()$plotlyscenehandle + params=list(...) + opacity <- if("alpha" %in% names(params)) params$alpha else 1 + } + + reg=NULL + if(is.list(x)) { + if(is.null(x$TransformParameters)) + stop("This looks like an in memory elastix registration but I can't find the TransformParameters field!") + } else if(is.character(x)) { + x <- read.elastixreg(x) + } else stop("Don't know what to do with this input!") + + type=match.arg(type) + + # make a fake im3d object using the grid information + gridim=im3d(dims = x$GridSize, voxdims = x$GridSpacing, origin = x$GridOrigin) + # now make an Nx3 matrix + grid=imexpand.grid(gridim) + if(!isTRUE(all.equal(dim(grid), dim(x$TransformParameters)))) + stop("Mismatch between transform parameters and expected grid size!") + gridt=grid+x$TransformParameters + + if (plotengine == 'rgl'){ + if(type %in% c("b", "p")) { + plot3d(gridt, xlab='X', ylab = "Y", zlab = "Z", ...) + } else { + plot3d(gridt, xlab='X', ylab = "Y", zlab = "Z", type='n', ...) + } + interleaved=matrix(t(cbind(grid,gridt)),ncol=3,byrow=T) + if(type %in% c('l','b')) { + segments3d(interleaved, ...) + } + + } else{ + plotdata <- as.data.frame(gridt) + names(plotdata) <- c('X','Y','Z') + psh <- psh %>% + plotly::add_trace(data = plotdata, x = ~X, y = ~Y , z = ~Z, + hoverinfo = "none",type = 'scatter3d', mode = 'markers', + opacity = opacity, marker=list(color = 'black', size = 3)) %>% + plotly::layout(showlegend = FALSE, scene=list(camera=.plotly3d$camera)) + assign("plotlyscenehandle", psh, envir=.plotly3d) + psh + } +} diff --git a/R/elastix_io.R b/R/elastix_io.R new file mode 100644 index 00000000..b48db9d2 --- /dev/null +++ b/R/elastix_io.R @@ -0,0 +1,87 @@ +read.elastixreg <- function(x, ...) { + l=read.elastix(x) + if(isTRUE(l$Transform=='BSplineTransform')) { + l$TransformParameters=matrix(l$TransformParameters, ncol=3, byrow=FALSE) + } + # reginmem is supposed to help handle the situation when we have read the + # registration into memory + class(l)=union(c('elastixreg', 'reginmem', 'reg'), class(l)) + l +} + +# internal function to read the elastix files into an R list +read.elastix<-function(con, CheckLabel=TRUE){ + + l=list() + + if(is.character(con)) { + filename=con + con=file(filename,'rt') + on.exit(close(con)) + t=readLines(filename,1) + if( !any(grep("(Transform",t[1],fixed=TRUE,useBytes=TRUE)) ) + stop(paste("This doesn't appear to be an elastix registration:",filename)) + } + + checkLabel=function(label) { + if( any(names(l)==label) ){ + newlabel=make.unique(c(names(l),label))[length(l)+1] + warning(paste("Duplicate item",label,"renamed",newlabel)) + label=newlabel + } + label + } + + removeBrackets <- function(x) { + if(nchar(x)<2) stop("Invalid line!") + if(substr(x,1,1)!="(") stop("No opening bracket!") + if(substr(x,nchar(x),nchar(x))!=")") stop("No closing bracket!") + substr(x, 2, nchar(x)-1L) + } + + # Should this check to see if the connection still exists? + # in case we want to bail out sooner + while ( isTRUE(isOpen(con)) ){ + thisLine<-readLines(con,1) + # no lines returned - ie end of file + if(length(thisLine)==0) break + + # trim and split it up by white space + thisLine=trimws(thisLine) + + # skip if this is a blank line + if(nchar(thisLine)==0) next + # skip if this is a comment + if(isTRUE(substr(thisLine, 1, 2) == "//")) next + + thisLine=removeBrackets(thisLine) + items=strsplit(thisLine," ",fixed=TRUE)[[1]] + + if(length(items)==0) next + # get the label and items + label=items[1]; items=items[-1] + + firstItemFirstChar=substr(items[1],1,1) + if (firstItemFirstChar=="\""){ + # dequote quoted string + # can do this by using a textConnection + tc=textConnection(thisLine) + items=scan(tc,what="",quiet=TRUE)[-1] + close(tc) + # attr(items,"quoted")=TRUE + } else { + items=as.numeric(items) + } + + # set the list element! + if(CheckLabel) + label=checkLabel(label) + l[[label]]=items + } + + if(isTRUE(try(file.exists(filename)))){ + attr(l,"file.info")=file.info(filename) + } + return(l) +} + diff --git a/man/xformpoints.Rd b/man/xformpoints.Rd index 0979b5fd..81bb977d 100644 --- a/man/xformpoints.Rd +++ b/man/xformpoints.Rd @@ -1,6 +1,7 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/xformpoints.R -\name{xformpoints} +% Please edit documentation in R/elastix.R, R/xformpoints.R +\name{xformpoints.elastixreg} +\alias{xformpoints.elastixreg} \alias{xformpoints} \alias{xformpoints.character} \alias{xformpoints.cmtkreg} @@ -8,6 +9,14 @@ \alias{xformpoints.default} \title{Transform 3D points using a registration, affine matrix or function} \usage{ +\method{xformpoints}{elastixreg}( + reg, + points, + transformtype = c("warp", "affine"), + direction = NULL, + ... +) + xformpoints(reg, points, ...) \method{xformpoints}{character}(reg, points, ...) @@ -33,8 +42,6 @@ path(s) to registrations on disk (see details).} \item{points}{Nx3 matrix of points} -\item{...}{Additional arguments passed to methods} - \item{transformtype}{Which transformation to use when the CMTK file contains both warp (default) and affine} @@ -42,6 +49,8 @@ both warp (default) and affine} space (called \strong{inverse} by CMTK) or from reference to sample space (called \strong{forward} by CMTK). Default (when \code{NULL} is inverse).} +\item{...}{Additional arguments passed to methods} + \item{FallBackToAffine}{Whether to use the affine transformation for points that fail to transform under a warping transformation.} }