R/UNUImageProcessingFunctions.R

Defines functions NrrdSave NrrdFlip NrrdMerge NrrdProject NrrdCrc NrrdTestDataLength NrrdTestIntegrity NrrdQuantize NrrdHisto callunu NrrdResample NrrdMinMax Nrrd2op

Documented in Nrrd2op NrrdCrc NrrdFlip NrrdHisto NrrdMerge NrrdMinMax NrrdProject NrrdResample NrrdSave NrrdTestDataLength NrrdTestIntegrity

# Some Image Processing Functions making use of the unu command line tool
# unu = Utah Nrrd Utilities - ie they operate on nrrd files
# See http://teem.sourceforge.net/unrrdu/index.html


#' Binary operation on two nrrds, or on a nrrd and a constant
#'
#' @param infiles Input files
#' @param outfile Output files
#' @param fun Binary function or operator
#' @param gzip Whether to gzip the resultant file
#' @param CreateDirs Whether to create implied output directories
#' @param Verbose Print status messages
#' @param Force Whether to overwrite an existing file
#'
#' @return The number of input files
#' @export
#' @importFrom nat.utils RunCmdForNewerInput
Nrrd2op<-function(infiles,outfile,fun=c("max","min","+", "-", "x", "/"),
	gzip=FALSE,CreateDirs=TRUE,Verbose=TRUE,Force=FALSE){
	if(length(infiles)<1) return(NULL)
	# Do nothing if inputs are older than output unless Force=T
	if(!Force && !RunCmdForNewerInput(NULL,infiles,outfile)) return (NULL)
	if(CreateDirs && !file.exists(dirname(outfile))) dir.create(dirname(outfile))
  fun=match.arg(fun)

  # 1 input: f1 -> output
	if(length(infiles)==1) file.copy(infiles[1],outfile)
	else if(length(infiles)>1) {
		# 2 or more inputs: f1 x f2 -> output
		cmd=paste("unu 2op",fun,
			shQuote(infiles[1]),shQuote(infiles[2]),"-o",shQuote(outfile))
		system(cmd)
		if(Verbose) cat(".")
	}
	# 3 or more inputs: f3 ... x output -> output
	if(length(infiles)>2) {
		for (f in infiles[-(1:2)]){
			cmd=paste("unu 2op",fun,
				shQuote(f),shQuote(outfile),"-o",shQuote(outfile))
			system(cmd)
			if(Verbose) cat(".")
		}
	}
	if(length(infiles)>1 && gzip)
		system(paste("unu save --format nrrd --encoding gzip","-i",shQuote(outfile),"-o",shQuote(outfile)))
	return(length(infiles))
}

#' Return range of values in nrrd file (interface to unu minmax command)
#' @param filename nrrd file containing data
#' @param blind8 Assume 8 bit data has range 0-255 (or -127-127), default F.
#' @param ... passed to \code{\link{system}} function
#' @return c(min,max) or c(NA,NA) like R's \code{\link{range}} function
#' @author jefferis
#' @export
NrrdMinMax<-function(filename, blind8=FALSE, ...){
	args=shQuote(filename)
	if(blind8) args=c(args,"-blind8 true")
	minmax=.callunu("minmax",args,intern=TRUE,...)
	if(length(minmax)==0) return(c(NA,NA))
	as.numeric(sub("(min|max): ","",minmax[1:2]))
}

#' Resample a nrrd with simple or complex filters
#'
#' @description wraps unu resample
#'
#'   When size is explicitly an integer (e.g. c(512L,512L,100L) then it is
#'   assumed to represent the size of the desired output in pixels.
#'
#'   unu resample defaults to cell centering in the asbence of information, but
#'   I prefer node.
#'
#'   * downsampling a **cell** image (e.g. x2) will result in \itemize{
#'
#'   \item an origin shift **iff** an origin was specified in the input file
#'
#'   \item a change in the space directions that is not an exact multiple of the
#'   downsampling factor } Downsampling a node image (e.g. x2) will result in:
#'   \itemize{
#'
#'   \item no change in origin
#'
#'   \item a doubling in the space directions field }
#' @param infile,outfile Input/output files
#' @param size Numeric vector of scale factors (NA=> don't touch) or integer
#'   vector of pixel dimensions.
#' @param voxdims Target voxel dimensions
#' @param centering Whether to assume node or cell centering if not specified in
#'   input file.
#' @param otherargs Passed to unu resample
#' @param gzip Whether or not to gzip encode the output file?
#' @param suffix A suffix added to the input file to construct an output
#'   filename
#' @param UseLock Whether to use file-based locking to enable simple
#'   parallelisation
#' @param ... Additional params passed to .callunu
#' @return The output file
#' @author jefferis
#' @export
#' @importFrom nat.utils makelock
#' @inheritParams Nrrd2op
NrrdResample<-function(infile,outfile,size,voxdims=NULL,
		centering=c("cell","node"),otherargs=NULL,gzip=TRUE, suffix=NULL,
		CreateDirs=TRUE,Verbose=TRUE,Force=FALSE,UseLock=FALSE,...){

	centering=match.arg(centering)
	if(!file.exists(infile)) stop("infile: ",infile," does not exist")
	if(!is.null(voxdims)){
		# we have a target voxel size instead of a standard size specification
		# first fetch existing voxel size
		size=nat::nrrd.voxdims(infile)/voxdims
	}
	if(is.integer(size)) size=paste("--size",paste(size,collapse=" "))
	else {
		size=paste("--size",paste("x",size,sep="",collapse=" "))
		# eg size=c(0.5,0.5,NA) => "--size x0.5 x0.5 ="
		size=sub("xNA","=",size)
	}

	if (missing(outfile)) {
		# we only want to add a default suffix if we are putting the output
		# file into the same directory
		if(is.null(suffix))
			suffix=gsub(" ","",sub("--size","",size))
		outfile=sub("\\.nrrd$",paste(suffix,"\\.nrrd",sep=""),infile)
	} else {
		# here we have been passed an output directory rather than a file
		if(file.exists(outfile) && file.info(outfile)$isdir)
			outfile=file.path(outfile,
				sub("\\.nrrd$",paste(suffix,"\\.nrrd",sep=""),basename(infile)))
	}

	# return outfile to signal output exists (we just didn't make it now)
	if(!Force && !RunCmdForNewerInput(NULL,infile,outfile)) return (outfile)
	# NB createdirs only makes sense if we have directly specified outfile
	if(CreateDirs && !file.exists(dirname(outfile))) dir.create(dirname(outfile),recursive = TRUE)
	lockfile=paste(outfile,sep=".","lock")
	# return FALSE to signal output doens't exist
	if(UseLock){
		if(makelock(lockfile))
			on.exit(unlink(lockfile))
		else
			return(FALSE)
	}
	otherargs=c(otherargs,"--center",centering)
	cmdargs=paste(size, paste(otherargs,collapse=" "),"-i",shQuote(infile))
	if(gzip)
		cmdargs=paste(cmdargs,"| unu save -f nrrd -e gzip")
	rval=.callunu("resample",paste(cmdargs,"-o",shQuote(outfile)),...)
	if(rval!=0) stop("error ",rval," in unu resample")
	return(outfile)
}

.callunu<-function(cmd,args,unu="unu",DryRun=FALSE,...){
	fullcmd=paste(unu,cmd,paste(args,collapse=" "))
	if(DryRun) {
		cat(fullcmd,"\n",sep="")
		return(0)
	}
	else system(fullcmd,...)
}

#' Title
#'
#' @details \code{bins} defaults to 1000 for float/double image data, to the
#'   number of levels between \code{min} and \code{max} for integer data. If
#'   \code{min} or \code{max} is not specified, they will be set to the minimum
#'   (maximum) value found in the input data.
#' @param maskfile Path to another image file defining which pixels in the input
#'   file to include.
#' @param bins Number of bins for the histogram (see details)
#' @param min,max Minimum and maximum values for the histogram (see details)
#' @param blind8 Whether to assume the full range (0,255) for 8 bit images
#'   without checking
#' @param ... Additional arguments passed to private \code{.callunu} function.
#'
#' @return Output file
#' @inheritParams NrrdResample
#' @importFrom nat read.nrrd.header
#' @export
NrrdHisto<-function(infile,outfile=sub("\\.([^.]+)$",".histo.\\1",infile),
	maskfile,bins,min,max,blind8=TRUE,...){
	h=read.nrrd.header(infile)
	nrrdType=nat:::.standardNrrdType(h$type)
	unuhistooptions=''
	if(nrrdType%in%c("uint8","int8") && blind8 && missing(min) && missing(max)){
		# this is an 8 bit image
		# and we're going to use unu histo's blind8 option (default) to
		# assume uchar [0,255], signed char is [-128,127]
		bins=256
	} else {
	    if (missing(min) || missing(max)) {
			# calculate the minimum and maximum
			r=NrrdMinMax(infile)
			if(any(is.na(r))) stop("Unable to find min and max from: ",infile)
			if(missing(min)) min=r[1]
			if(missing(max)) max=r[2]
    	}
		unuhistooptions=paste("-min",min,"-max",max)
	}
	if(missing(bins)){
		# check if this is a float data type
		if(nrrdType%in%c("float","double"))
			bins=1000
		else {
			bins=as.integer((max-min)+1)
			if(bins>2^16) bins=1000
		}
	}
	unuhistooptions=paste(unuhistooptions,"-b",bins)
	if(!missing(maskfile)) unuhistooptions=paste(unuhistooptions,"-w",shQuote(maskfile))
	.callunu("histo",paste(unuhistooptions,"-i",shQuote(infile),"-o",shQuote(outfile)),...)
	return(outfile)
}

NrrdQuantize<-function(infile,outfile,min,max,bits=c("8","16","32"),
	gzip=FALSE,CreateDirs=TRUE,Verbose=TRUE,Force=FALSE,UseLock=FALSE){

	# Do nothing if inputs are older than output unless Force=T
	if(!Force && !RunCmdForNewerInput(NULL,infile,outfile)) return (FALSE)
	if(CreateDirs && !file.exists(dirname(outfile))) dir.create(dirname(outfile),recursive = TRUE)

	bits=match.arg(bits)

	lockfile=paste(outfile,sep=".","lock")
	if(UseLock && !makelock(lockfile)) return (FALSE)
	# unu quantize -b 8 -min 7397.386 -max 65535 -i AL-PNs_mask_MF_final_IS2_SAIA24-1_02.nrrd \
	# | unu save -f nrrd -e gz -o AL-PNs_mask_MF_final_IS2_SAIA24-1_02-quant.nrrd

	minmaxopts=paste(ifelse(missing(min),"",paste("-min",min)),
		ifelse(missing(max),"",paste("-max",max)))
	cmd=paste("unu quantize -b",bits,minmaxopts,"-i",shQuote(infile))
	if(gzip) cmd=paste(cmd,"| unu save -f nrrd -e gzip -o",shQuote(outfile))
	else cmd=paste(cmd,"-o",shQuote(outfile))
	system(cmd)
	if(Verbose) cat(".")
	if(UseLock) unlink(lockfile)
	return(TRUE)
}

#' Test integrity of a NRRD file using gzip/bzip2 checksum
#'
#' @param defaultReturnVal The value to return when file cannot be checked
#'
#' @return A logical indicating if the check passed
#' @export
#'
#' @inheritParams NrrdResample
#' @seealso \code{\link{NrrdTestDataLength}} which is much faster
NrrdTestIntegrity<-function(infile,defaultReturnVal=TRUE){
	# Tests integrity of a compressed nrrd file using crc check
	if(!file.exists(infile)) return(NA)
	h=nat::read.nrrd.header(infile)
	if(tolower(h$encoding)%in%c("gz","gzip")) testprog='gzip'
	else if(tolower(h$encoding)%in%c("bz2","bzip2")) testprog='bzip2'
	else {
		warning("Unable to test integrity of nrrd file with encoding",h$encoding)
		return(defaultReturnVal)
	}
	testval=system(paste("unu data",infile," | ",testprog,"-t"),ignore.stderr=TRUE)
	return(testval==0)
}

#' Quick test that gzipped NRRD has correct data length
#'
#' @details Tests integrity of a nrrd file by checking that the data block is as
#'   long as it should be. For a gzip file, it checks that the last 4 bytes
#'   encode the length of the uncompressed data. This says nothing about the
#'   contents but does ensure that the file has not been truncated which is much
#'   the most common problem.
#' @return A logical indicating if the check passed
#' @export
#' @inheritParams NrrdTestIntegrity
#' @seealso \code{\link{NrrdTestIntegrity}} which is more comprehensive but much
#'   slower
NrrdTestDataLength<-function(infile,defaultReturnVal=TRUE){
	if(!file.exists(infile)) return(NA)
	# need to read this full header in order to get the size in bytes of
	# the uncompressed data
	fullh=nat::read.nrrd(infile,ReadData=F,AttachFullHeader=T)
	h=attr(fullh,'header')
	if(tolower(h$encoding)%in%c("gz","gzip")) enc='gzip'
	else if(tolower(h$encoding)=='raw') enc='raw'
	else {
		warning("Unable to test data length of nrrd file with encoding",h$encoding)
		return(defaultReturnVal)
	}
	# figure out how many bytes we are expecting
	dataLength=attr(fullh,"datablock")$n*attr(fullh,"datablock")$size
	datafile=nat:::nrrd.datafiles(h)
	if(length(datafile)>1){
		warning("Don't yet know how to handle detached headers with more than one data file")
		return(defaultReturnVal)
	}
	if(enc=='gzip'){
		# Check gzip file
		# last 4 bytes of gzip contain length of compressed stream
		con=file(datafile,open='rb')
		on.exit(close(con))
		seek(con,-4,origin='end')
		# how do we handle endianess?
		uncompressedBytes=readBin(con,what=integer(),size=4,n=1)
		if(uncompressedBytes>(2^32-1)) {
			warning("Don't know how to check integrity of files >= 2^32 bytes")
			return(defaultReturnVal)
		}
		return(dataLength==uncompressedBytes)
	} else {
		# TODO fix handling of detached nrrd files with raw encoding
		start=attr(fullh,'datablock')$datastartpos
		end=file.info(datafile)$size
		# TODO decide what we should say if we have MORE than enough data
		return(dataLength>=(end-start))
	}
}

#' Find the CRC (hash) of a gzip encoded NRRD file
#'
#' @param UseGzip Whether to use the gzip command line tool to extract CRC
#'   (default \code{FALSE}, see details)
#' @param FastHeader Parse only essential parts of NRRD header to save time.
#'
#' @details: Defaults to a quick method based on knowledge of gzip file format
#' from: \url{http://www.gzip.org/zlib/rfc-gzip.html} and assumption that there
#' is only one member in gzip data. Can also use gzip but this is much slower
#' since have to copy unu data to temporary file.
#' @inheritParams NrrdResample
NrrdCrc<-function(infile, UseGzip=FALSE, FastHeader=TRUE){

	if(!file.exists(infile)) return(NA)

	ext=tolower(sub(".*\\.([^.]+)$","\\1",basename(infile)))
	if(ext=='nhdr') FastHeader=FALSE

	h=NULL
	if(FastHeader){
		# quick and dirty reading of header
		con<-file(infile,'rb')
		on.exit(close(con))
		magic=readBin(con,what=raw(),5)
		if(!isTRUE(all.equal(magic,as.raw(c(0x4e, 0x52, 0x52, 0x44, 0x30))))){
			warning("This is not a nrrd")
			return(NA)
		}
		while( length(l<-readLines(con, 1))>0 && l!="" ){
			if(nchar(l)>=12 && substr(l, 1, 9)=="encoding:"){
				if(substr(l, 11, 12)!="gz"){
					warning("This is not a gzip encoded nrrd")
					return(NA)
				}
			}
		}
	} else {
		h=nat::read.nrrd.header(infile)
		if(tolower(h$encoding)%in%c("gz","gzip")) {
			# testprog='gzip'
		} else {
			warning("This is not a gzip encoded nrrd")
			return(NA)
		}
	}

	if(UseGzip){
		tmp=tempfile()
		on.exit(unlink(tmp),add=TRUE)
		system(paste("unu data ",shQuote(infile)," > ",shQuote(tmp)))
		x=system(paste("gzip -lv",shQuote(tmp)),intern=TRUE)
		crc=try(strsplit(x[2],"[ ]+")[[1]][[2]])
		if(inherits(crc,'try-error')) crc=NA
	} else {
		# TODO Fix handling of nhdr files
		if(!is.null(h) && !is.null(h$datafile)){
			if(length(h$datafile)>1)
				stop("Don't know how to handle more than 1 datafile")
			if(dirname(h$datafile)==".")
				h$datafile=file.path(dirname(infile),h$datafile)
			con=file(h$datafile,open='rb')
			on.exit(close(con),add=TRUE)
		} else if(!FastHeader){
			con=file(infile,open='rb')
			on.exit(close(con),add=TRUE)
		}
		seek(con,-8,origin='end')
		# TODO check endian issues (what happens if CRC was from opposite endian platform?)
		crc=readBin(con,integer(),size=4)
		crc=format(as.hexmode(crc),width=8)
	}
	crc
}


#' Make projections along a nrrd axis, typically also rescaling
#'
#' @details gamma: Just as in xv, the gamma value here is actually the
#'   reciprocal of the exponent actually used to transform the values. Note also
#'   that for \code{cropmin,cropmax} the special value M can be used to indicate
#'   the maximum index for that axis (i.e. n-1 when there are n samples).
#'
#'   scale: For each axis, information about how many samples in output:
#'   \itemize{
#'
#'   \item "=": leave this axis completely untouched: no resampling whatsoever
#'   \item "x<float>": multiply the number of input samples by <float>, and
#'   round to the nearest integer, to get the number of output samples. Use "x1"
#'   to resample the axis but leave the number of samples unchanged
#'
#'   \item "/<float>": divide number of samples by <float>
#'
#'   \item "+=<uint>", "-=<uint>": add <uint> to or subtract <uint> from number
#'   input samples to get number output samples
#'
#'   \item "<uint>": exact number of output samples
#'
#'   \item "a": resample this axis to whatever number of samples preserves the
#'   aspect ratio of other resampled axes. Currently needs to be used on all but
#'   one of the resampled axes, if at all. (1 or more sampling specifications) }
#'
#'   Kernel: The kernel to use for resampling. Kernels logically live in the
#'   input index space for upsampling, and in the output index space for
#'   downsampling. Possibilities include: \itemize{ \item "box": nearest
#'   neighbor interpolation on upsampling, and uniform averaging on downsampling
#'   \item "cheap": nearest neighbor interpolation for upsampling, and
#'   non-blurring sub-sampling (pick subset of input samples) on downsampling
#'
#'   \item "tent": linear interpolation
#'
#'   \item "cubic:B,C": Mitchell/Netravali BC-family of cubics: "cubic:1,0":
#'   B-spline; maximal blurring "cubic:0,0.5": Catmull-Rom; good interpolating
#'   kernel
#'
#'   \item "c4h": 6-sample-support, C^4 continuous, accurate
#'
#'   \item "c4hai": discrete pre-filter to make c4h interpolate
#'
#'   \item "bspl3", "bspl5", "bspl7": cubic (same as cubic:1,0), quintic, and
#'   7th order B-spline
#'
#'   \item "bspl3ai", "bspl5ai", "bspl7ai": discrete pre-filters to make bspl3,
#'   bspl5, bspl7 interpolate
#'
#'   \item "hann:R": Hann (cosine bell) windowed sinc, radius R
#'
#'   \item "black:R": Blackman windowed sinc, radius R
#'
#'   \item "gauss:S,C": Gaussian blurring, with standard deviation S and cut-off
#'   at C standard deviations
#'
#'   \item "dgauss:S,C": Lindeberg's discrete Gaussian. default: "cubic:0,0.5" }
#'
#'   Crop max, min: \itemize{
#'
#'   \item <int> gives 0-based index
#'
#'   \item M, M+<int>, M-<int> give index relative to the last sample on the
#'   axis (M == #samples-1).
#'
#'   \item m+<int> give index relative to minimum axis (only relevant for max)
#'
#'   }
#' @param scale How to scale each axis (see details)
#' @param kernel Type of smoothing kernel to use (see details)
#' @param gamma Exponent for the mapping of input to output values (see
#'   details)
#' @param outfile Optional path to output file (constructed automatically when
#'   missing)
#' @param axis Number indicating 0-indexed axis or character "x", "y" or "z"
#' @param measure Character vector indicating summary function to apply to
#'   values in each column. Choose from \code{c("max", "min", "mean", "median",
#'   "mode", "variance", "skew", "intc", "slope", "error", "sd", "product",
#'   "sum", "L1", "L2", "Linf")}
#' @param cropmin,cropmax triples defining the lower and higher corners of
#'   bounding box
#' @return Logical indicating success
#' @export
#' @seealso \code{\link{NrrdSave}}, \code{\link{NrrdMerge}}
#' @inheritParams NrrdResample
#' @examples
#' \dontrun{
#' # NB cropping from 50th slice to 50 before final slice
#' NrrdProject(infile,axis='z',cropmin='0 0 50',cropmax='0 0 M-50')
#' }
NrrdProject<-function(infile,outfile,axis,
	measure=c("max", "min", "mean", "median", "mode", "variance", "skew",
	"intc", "slope", "error", "sd", "product", "sum", "L1", "L2", "Linf"),
	scale="x0.3333 x0.333",kernel='cheap',gamma=NA,cropmin=NULL,cropmax=NULL,
	suffix=NULL,
	CreateDirs=TRUE,Verbose=TRUE,Force=FALSE,UseLock=FALSE){
	measure=match.arg(measure)
	if (missing(outfile)) {
		if(is.null(suffix)) suffix=paste("-",axis,measure,sep="")
		outfile=sub("\\.nrrd$",paste(suffix,".png",sep=""),infile)
	}
	if(!file.exists(infile)) stop("infile: ",infile," does not exist")
	if(is.character(axis)) {
		axis=match(tolower(axis),c("x",'y','z'))-1
		if(is.na(axis)) stop("Unable to recognise nrrd axis")
	}
	# return TRUE to signal output exists (we just didn't make it now)
	if(!Force && !RunCmdForNewerInput(NULL,infile,outfile)) return (TRUE)
	if(CreateDirs && !file.exists(dirname(outfile))) dir.create(dirname(outfile),recursive = TRUE)
	lockfile=paste(outfile,sep=".","lock")
	# return FALSE to signal output doens't exist
	if(UseLock && !makelock(lockfile)) return (FALSE)
	if(is.numeric(scale)) scale=paste(scale,collapse=" ")
	if(!is.null(cropmax)||!is.null(cropmin)){
		cmd=paste("unu crop",
			ifelse(is.null(cropmax),'',paste("-max",cropmax)),
			ifelse(is.null(cropmin),'',paste("-min",cropmin)),
			"-i",shQuote(infile),
			"| unu resample -s",scale," = -k ",kernel,
			"| unu project -a",axis,"-m ",measure," | unu quantize -b 8 ")
	} else cmd=paste("unu resample -s",scale," = -k ",kernel," -i",shQuote(infile),
		"| unu project -a",axis,"-m ",measure," | unu quantize -b 8 ")

	if(!is.na(gamma)){
		cmd=paste(cmd,"| unu gamma --gamma",gamma)
	}
	cmd=paste(cmd,"| unu save -f png -o",shQuote(outfile))
	rval = system(cmd)
	if(Verbose) cat(".")
	if(UseLock) unlink(lockfile)
	if(rval!=0) stop("unu error ",rval," in NrrdProject")
	return(TRUE)
}

#' Merge multiple NRRD images into a multichannel PNG
#'
#' @param outdir Output directory
#' @param axis 0-indexed axis to join (merge) along (defaults to 0)
#' @param ... Currently ignored
#'
#' @inheritParams Nrrd2op
#' @inheritParams NrrdResample
#' @return Logical value indicating whether output exists
#' @export
NrrdMerge<-function(infiles,outdir=NULL,outfile=NULL,axis=0,
	CreateDirs=TRUE,Verbose=TRUE,Force=FALSE,UseLock=FALSE,...){
	if(is.null(outfile) && is.null(outdir)) {
		outfile=sub("\\.[^.]+$","-merge.png",infiles[1])
	} else if (is.null(outfile)){
		outfile=file.path(outdir,basename(infiles[1]))
		outfile=sub("\\.[^.]+$",".png",outfile)
	} else {
		# make sure this is a png
		outfile=sub("\\.[^.]+$",".png",outfile)
	}

	if(!all(file.exists(infiles)))
	  stop("one or more infiles: ",infiles," do not exist")
	# return TRUE to signal output exists (we just didn't make it now)
	if(!Force && !RunCmdForNewerInput(NULL,infiles,outfile)) return (TRUE)
	if(CreateDirs && !file.exists(dirname(outfile)))
	  dir.create(dirname(outfile),recursive = TRUE)
	lockfile=paste(outfile,sep=".","lock")
	# return FALSE to signal output doens't exist
	if(UseLock && !makelock(lockfile)) return (FALSE)

	cmd=paste("unu join -a", axis, "-incr -i",
	          paste(shQuote(infiles),collapse=" "),
	          "-o",shQuote(outfile))
	cat(cmd, "\n")
	rval = system(cmd)
	if(Verbose) cat(".")
	if(UseLock) unlink(lockfile)
	if(rval!=0) stop("unu error ",rval," in NrrdProject")
	return(TRUE)
}

#' Flip a NRRD image along one axis
#'
#' @param axes Which axes to flip (0-indexed integer)
#' @param suffix Suffix to add to input file to construct output file.
#' @param endian Whether output image should be big (default) or little endian
#'   for multibyte data types (essentially all other than 8 bit)
#' @param OverWrite Whether to overwrite (yes or no) or update output image.
#'
#' @return \code{TRUE} if output exists
#' @export
#'
#' @inheritParams NrrdResample
#' @inheritParams NrrdProject
NrrdFlip<-function(infile,outfile,axes,suffix=NULL,endian=c("big","little"),
	CreateDirs=TRUE,Verbose=TRUE,UseLock=FALSE, OverWrite=c("no","update","yes")){
	# TODO would be nice if we could
	# a) have an absolute flip mode that checks the nrrd content field
	# b) similarly checks whether output image has been flipped accordingly

	if(is.logical(OverWrite)) OverWrite=ifelse(OverWrite,"yes","no")
	else OverWrite=match.arg(OverWrite)

	endian=match.arg(endian)
	if (missing(outfile)) {
		if(is.null(suffix)) suffix=paste("-flip",paste(axes,collapse=""),sep="")
		outfile=sub("\\.nrrd$",paste(suffix,".nrrd",sep=""),infile)
	}

	if(!file.exists(infile)) stop("infile: ",infile," does not exist")

	# return TRUE to signal output exists (whether or not we made it)
	if(file.exists(outfile)){
		if(OverWrite=="no"){
			if(Verbose) cat("Output",outfile,"already exists; use OverWrite=\"yes\" or \"update\" to overwrite or update\n")
			return(TRUE)
		} else if(OverWrite=="update"){
			# check modification times
			if(!RunCmdForNewerInput(NULL,infile,outfile)) return (TRUE)
		} else if(Verbose) cat("Overwriting",outfile,"because OverWrite=\"yes\"\n")
	}

	if(CreateDirs && !file.exists(dirname(outfile))) dir.create(dirname(outfile),recursive = TRUE)
	lockfile=paste(outfile,sep=".","lock")
	# return FALSE to signal output doesn't (yet) exist
	if(UseLock && !makelock(lockfile)) return (FALSE)
	on.exit(unlink(lockfile))

	# First axis
	cmd=paste("unu flip -a",axes[1],"-i",shQuote(infile))
	# any additional axes
	for(axis in axes[-1]){
		cmd=paste(cmd,"| unu flip -a",axis)
	}
	# save
	cmd=paste(cmd," | unu save -f nrrd -e gz -en",endian,"-o",shQuote(outfile))

	rval = system(cmd)
	if(Verbose) cat(".")
	if(rval!=0) stop("unu error ",rval," in NrrdProject")
	return(TRUE)
}

#' Save a nrrd file in a different encoding / format
#'
#' In outfile is missing, will overwrite infile
#' @param infile,outfile Paths to input and output nrrds
#' @param format (one of 6 supported by unu save, default nrrd)
#' @param encoding (one of 5 supported by unu save, default gz when nrrd)
#' @param CreateDirs Whether to make missing direcories implied by output path
#' @param UseLock Whether to make a lockfile (useful for parallel processing)
#' @param DryRun Return command instead of running it
#' @return path to ouput file
#' @export
#' @seealso \code{\link{NrrdCrc}}
NrrdSave<-function(infile,outfile,format=c("nrrd","pnm","text","vtk","png","eps"),
  encoding=ifelse(format=='nrrd','gzip','raw'),
  CreateDirs=TRUE,UseLock=FALSE,DryRun=FALSE){
  format=match.arg(format)
  encodings=c("raw","ascii","hex","gzip","bzip2")
  encoding=encodings[pmatch(encoding,encodings)]
  if(is.na(encoding)) stop("Invalid encoding")

  if(missing(outfile)) {
    outfile=infile
    samefile=TRUE
  } else samefile=FALSE

  cmd=paste('unu save','--format',format,'-e',encoding,'-i',shQuote(infile),'-o',shQuote(outfile))
  if(DryRun) return(cmd)

  if(CreateDirs && !file.exists(dirname(outfile))) dir.create(dirname(outfile),recursive = TRUE)
  lockfile=paste(outfile,sep=".","lock")
  # return FALSE to signal output doesn't (yet) exist
  if(UseLock && !nat.utils::makelock(lockfile)) return (FALSE)
  on.exit(unlink(lockfile))

  rval=system(cmd)==0
  attr(rval,'outfile')=outfile
  rval
}
jefferis/runu documentation built on July 26, 2017, 4:05 a.m.