Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
5 changed files
with
336 additions
and
4 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 | ||
} | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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) | ||
} | ||
|
Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.
Oops, something went wrong.