R/dotprops.R

Defines functions prune.default prune.neuronlist prune.dotprops prune.neuron prune subset.dotprops plot.dotprops plot3d.dotprops all.equal.dotprops dotprops2swc dotprops.default dotprops.neuron dotprops.neuronlist dotprops.list dotprops.im3d dotprops.dotprops dotprops.character dotprops scale.dotprops Ops.dotprops as.dotprops is.dotprops

Documented in all.equal.dotprops as.dotprops dotprops dotprops dotprops.character dotprops.default dotprops.dotprops dotprops.im3d dotprops.neuron dotprops.neuronlist is.dotprops Ops.dotprops plot3d.dotprops plot.dotprops prune prune.default prune.dotprops prune.neuron prune.neuronlist scale.dotprops subset.dotprops

#' dotprops: Neurons as point clouds with tangent vectors (but no connectivity)
#' @param x Object to be tested/converted
#' @name dotprops
#' @export
is.dotprops<-function(x) inherits(x,"dotprops")

#' @rdname dotprops
#' @export
as.dotprops<-function(x, ...){
  if(is.null(x)) return (NULL)
  if(!is.dotprops(x)) class(x)=c("dotprops",class(x))
  if(is.null(colnames(x$points))) colnames(x$points) <-c("X","Y","Z") 
  x
}

#' Arithmetic for nat dotprops and surface objects
#' 
#' @param e1 A dotprops or surface (hxsurf, mesh3d) object
#' @param e2 A scalar or 3-vector that will be applied to the dotprops object
#' @return A new dotprops / surface object
#' @export
#' @rdname dotprops-arithmetic
#' @seealso \code{\link{scale.dotprops}}, \code{\link{Ops.neuron}}
#' @examples 
#' kcs20.shift=kcs20+c(2,3,4)
#' \donttest{
#' plot3d(kcs20, col='grey')
#' plot3d(kcs20.shift, col='red')
#' }
Ops.dotprops <- function(e1, e2=NULL) {
  r=e1
  e1=xyzmatrix(e1)
  # I don't exactly know why it is necessary to change this directly, but if not
  # NextMethod dispatches on original class of e1 even when I specify object=
  .Class=class(e1)
  lx=length(e2)
  if(lx==3) e1=t(e1) else if(lx>1) stop("expects a numeric vector of length 0, 1 or 3")
  res <- NextMethod(generic=.Generic)
  if(lx==3) res=t(res)
  xyzmatrix(r)=res
  r
}

#' @export
#' @rdname dotprops-arithmetic
Ops.mesh3d <- Ops.dotprops

#' @export
#' @rdname dotprops-arithmetic
Ops.hxsurf <- Ops.dotprops

#' @rdname scale.neuron
#' @include neuron.R
#' @export
#' @aliases scale.dotprops
#' @description note that \code{scale.dotprops} recalculates the tangent vectors
#'   after scaling the 3D coords. See \code{\link{dotprops}} for details.
scale.dotprops<-function(x, center=TRUE, scale=TRUE){
  xyzmatrix(x)<-scale(xyzmatrix(x),scale=scale,center=center)
  dotprops(x)
}

#' @description \code{dotprops} makes dotprops representation from raw 3D points
#'   (extracting vertices from S3 objects that have them)
#' @details \code{k} will default to 20 nearest neighbours when unset (i.e. when
#'   it has default value of NA) unless \code{x} is a dotprops object (when the
#'   original value of \code{k} is reused).
#' @param ... Additional arguments passed to methods
#' @export
#' @rdname dotprops
dotprops<-function(x, ...) UseMethod('dotprops')

#' @export
#' @rdname dotprops
dotprops.character <- function(x, ...) {
  fileName <- x
  x <- read.im3d(x)
  l <- dotprops(x, ...)
  attr(l,'file') <- fileName
  fi <- file.info(fileName)
  attr(l, 'mtime') <- fi$mtime
  attr(l, 'size') <- fi$size
  l
}

#' @description \code{dotprops.dotprops} will default to the original vale of 
#'   \code{k} and copy over all attributes that are not set by
#'   \code{dotprops.default}.
#' @method dotprops dotprops
#' @export
#' @export
#' @rdname dotprops
dotprops.dotprops<-function(x, k=attr(x,'k'), ...) {
  y=dotprops(xyzmatrix(x), k=k, ...)
  # copy over attributes, taking care not to overwrite any
  attin=attributes(x)
  attout=attributes(y)
  attributes(y)<-c(attout, attin[setdiff(names(attin), names(attout))])
  y
}

#' @export
#' @rdname dotprops
dotprops.im3d <- function(x, ...) {
  l <- ind2coord(x)
  l <- dotprops(l, ...)
  l
}

#' @export
dotprops.list<-function(x, ...) {
  # FIXME - change to an abstract base class for objects with 3D vertices
  # rather than the completely generic list
  dotprops(xyzmatrix(x), ...)
}

#' @description \code{dotprops.neuronlist} will run for every object in the 
#'   neuronlist using \code{\link{nlapply}}. \code{...} arguments will be passed to 
#'   \code{nlapply} in addition to the named argument \code{OmitFailures}.
#' @export
#' @method dotprops neuronlist
#' @rdname dotprops
#' @inheritParams nlapply
#' @seealso \code{\link{nlapply}}
dotprops.neuronlist<-function(x, ..., OmitFailures=NA) {
  nlapply(x, dotprops, ..., OmitFailures=OmitFailures)
}

#' @export
#' @method dotprops neuron
#' @param resample When finite, a new length to which all segmented edges will
#'   be resampled. See \code{\link{resample.neuron}}.
#' @rdname dotprops
dotprops.neuron<-function(x, Labels=NULL, resample=NA, ...) {
  if(is.finite(resample)) x=resample(x, stepsize = resample)
  if(is.null(Labels) || isTRUE(Labels)) Labels=x$d$Label
  else if(is.logical(labels) && labels==FALSE) Labels=NULL
  dotprops(xyzmatrix(x), Labels=Labels, ...)
}

#' @method dotprops default
#' @export
#' @rdname dotprops
#' @param k Number of nearest neighbours to use for tangent vector calculation 
#'   (set to k=20 when passed NULL)
#' @param Labels Vector of labels for each point e.g. identifying axon vs 
#'   dendrite. The default value \code{NULL} will produce class-specific default
#'   behaviour for different classes of input object, \code{TRUE} always uses 
#'   labels when an incoming object has them and \code{FALSE} never uses labels.
#' @param na.rm Whether to remove \code{NA} points (default FALSE)
#' @importFrom nabor knn
#' @references The dotprops format is essentially identical to that developed 
#'   in:
#'   
#'   Masse N.Y., Cachero S., Ostrovsky A., and Jefferis G.S.X.E. (2012). A 
#'   mutual information approach to automate identification of neuronal clusters
#'   in \emph{Drosophila} brain images. Frontiers in Neuroinformatics 6 (00021).
#'   \href{http://dx.doi.org/10.3389/fninf.2012.00021}{doi: 10.3389/fninf.2012.00021}
dotprops.default<-function(x, k=NULL, Labels=NULL, na.rm=FALSE, ...){
  # store labels from SWC format data if this is a neuron
  x=xyzmatrix(x)
  if(is.null(k)) k=20

  if(length(Labels) && length(Labels)!=nrow(x))
  stop("Length of Labels does not match number of points!")

  if(na.rm){
    narows=rowSums(is.na(x))>0
    if(any(narows)) {
      x=x[!narows,]
      # don't forget to remove labels for NA points as well
      if(length(Labels)) Labels=Labels[!narows]
    }
  }
  npoints=nrow(x)
  if(npoints<k) stop("Too few points to calculate properties")
  if(ncol(x)!=3) stop("points must be a N x 3 matrix")
  
  alpha=rep(0,npoints)
  vect=matrix(0,ncol=3,nrow=npoints)
  
  nns=knn(x, k=k)
  # transpose points to 3xN because 
  # R arithemtic of matric / vector operates column-wise
  pointst=t(x)
  for(i in 1:npoints){
    indNN=nns$nn.idx[i,]
    pt=pointst[,indNN]
    cpt=pt-rowMeans(pt)
    
    inertia=cpt%*%t(cpt)
    v1d1<-eigen(inertia,symmetric=TRUE)
    
    alpha[i]=(v1d1$values[1]-v1d1$values[2])/sum(v1d1$values)
    vect[i,]=v1d1$vectors[,1]
  }
  rlist=list(points=x,alpha=alpha,vect=vect)
  
  rlist$labels=Labels
  
  attr(rlist,'k')=k
  return(as.dotprops(rlist))
}

# internal function to convert a dotprops object to SWC
# representation.
dotprops2swc<-function(x, label=0L, veclength=1, radius=0) {
  
  # compute segments based on points and tangent vectors
  halfvect=x$vect/2*veclength
  starts=x$points-halfvect
  stops=x$points+halfvect
  interleaved=matrix(t(cbind(starts,stops)),ncol=3,byrow=T)
  
  n=nrow(x$points)
  # make sequence of parent ids
  # should look like -1, 1, -1, 2, -1, 3 ...
  parents=as.vector(t(cbind(-1L, 
                            seq.int(from=1L, length.out=n, by=2L)
  )))
  df=data.frame(PointNo=seq.int(2*n), Label=label)
  df[,c("X","Y","Z")]=interleaved
  df$R=radius
  df$Parent=parents
  df
}

#' all.equal method tailored to dotprops objects
#' 
#' @details This method is required because the direction vectors are computed
#'   using an eigenvector decomposition where the sign of the eigenvector is
#'   essentially random and subject to small numerical instabilities. Therefore
#'   it does not usually make sense to check the value of vect exactly.
#' @method all.equal dotprops
#' @param target,current dotprops objects to compare
#' @param check.attributes Whether to check attributes (false by default)
#' @param absoluteVectors Whether to check only the absolute value of
#'   eigenvectors for equality (default TRUE, see details)
#' @param ... Additional arguments passed to base \code{all.equal}.
#' @export
#' @examples
#' # equal using default 
#' kc1=kcs20[[1]]
#' kc1.recalc=dotprops(kc1)
#' # not equal due to differences in attributes and vectors
#' all.equal.default(kc1.recalc, kc1)
#' # still not equal because of tangent vector flipping
#' all.equal.default(kc1.recalc, kc1, check.attributes=FALSE)
#' # equal using appropriate method
#' stopifnot(isTRUE(all.equal(kc1.recalc, kc1)))
#' # NB identical when recalculated on same setup from same data
#' stopifnot(isTRUE(all.equal.default(kc1.recalc, dotprops(kc1))))
all.equal.dotprops<-function(target, current, check.attributes=FALSE, 
                             absoluteVectors=TRUE, ...){
  if(!is.dotprops(current)) stop("current is not a dotprop object")
  if(absoluteVectors){
    target$vect=abs(target$vect)
    current$vect=abs(current$vect)
  }
  NextMethod(check.attributes=check.attributes, ...)
}

#' 3D plots of dotprops objects using rgl package
#' 
#' @details Tangent vectors are plotted by \code{segments3d} and centered on the
#'   relevant point. Points are plotted by \code{points3d}.
#'   
#'   \code{color} will be recycled by \code{points3d} and \code{segments3d}. 
#'   However in the special case that \code{color} has length equal to the 
#'   number of points in \code{x}, then it will be duplicated before being 
#'   passed to \code{segments3d} so that the result is that each vector is 
#'   coloured uniformly according to \code{color} (since segments3d expects 2
#'   colours for each line segment, blending them if they are different).
#' @param x A dotprops object
#' @param scalevecs Factor by which to scale unit vectors (numeric, default: 
#'   1.0)
#' @param alpharange Restrict plotting to points with \code{alpha} values in 
#'   this range to plot (default: null => all points). See 
#'   \code{\link{dotprops}} for definition of \code{alpha}.
#' @param color Character or numeric vector specifying colours for 
#'   points/vectors. See details.
#' @param PlotPoints,PlotVectors Whether to plot points and/or tangent vectors 
#'   (logical, default: tangent vectors only)
#' @param UseAlpha Whether to scale tangent vector length by the value of 
#'   \code{alpha}
#' @param ... Additional arguments passed to \code{points3d} and/or 
#'   \code{segments3d}
#' @return invisible list of results of rgl plotting commands
#' @method plot3d dotprops
#' @export
#' @seealso \code{\link{dotprops}, \link[rgl]{plot3d}, \link[rgl]{points3d}, 
#'   \link[rgl]{segments3d}}
#' @examples
#' \donttest{
#' open3d()
#' plot3d(kcs20[[1]])
#' clear3d()
#' plot3d(kcs20[[1]],col='red')
#' clear3d()
#' plot3d(kcs20[[1]],col='red',lwd=2)
#' plot3d(kcs20[[2]],col='green',lwd=2)
#' }
plot3d.dotprops<-function(x, scalevecs=1.0, alpharange=NULL, color='black', 
                          PlotPoints=FALSE, PlotVectors=TRUE, UseAlpha=FALSE, ...){
  # rgl's generic plot3d will dispatch on this
  if (!is.null(alpharange))
    x=subset(x,x$alpha<=alpharange[2] & x$alpha>=alpharange[1])
  rlist=list()
  if(PlotPoints){
    rlist$points=points3d(x$points, color=color, ...)
  }
    
  if(PlotVectors){
    if(length(color)>1 && length(color)==nrow(x$points)){
      color=rep(color,rep.int(2,length(color)))
    }
    halfvect=x$vect/2*scalevecs
    if(UseAlpha) halfvect=halfvect*x$alpha
    starts=x$points-halfvect
    stops=x$points+halfvect
    interleaved=matrix(t(cbind(starts,stops)),ncol=3,byrow=T)
    rlist$segments=segments3d(interleaved, color=color, ...)
  }
  invisible(rlist)
}

#' @rdname plot.neuron
#' @description \code{plot.dotprops} plots a 2D projection of a
#'   \code{\link{dotprops}} format object
#' @details \code{plot.dotprops} is limited in that 1) it cannot plot somata
#'   directly (this is handled by \code{\link{plot.neuronlist}}) and 2) it can
#'   only plot a frontal (XY) view.
#' @inheritParams plot.neuron
#' @inheritParams plot3d.dotprops
#' @export
#' @examples
#'
#' plot(kcs20[[1]], col='red')
#' # NB soma ignored
#' plot(kcs20[[1]], col='red', soma=TRUE)
#' plot(kcs20[1], col='red', soma=TRUE)
plot.dotprops<-function(x, scalevecs=1.0, alpharange=NULL, col='black', 
                        PlotPoints=FALSE, PlotVectors=TRUE, UseAlpha=FALSE,
                        asp=1, add=FALSE, axes=TRUE, tck=NA,
                        boundingbox=NULL, xlim=NULL, ylim=NULL, 
                        soma=FALSE, ...){
  if (!is.null(alpharange))
    x=subset(x,x$alpha<=alpharange[2] & x$alpha>=alpharange[1])
  
  
  if(!add) {
    if(is.null(boundingbox))
      boundingbox=boundingbox(x, na.rm=TRUE)
    if(is.null(xlim)) xlim=boundingbox[,1]
    if(is.null(ylim)) ylim=rev(boundingbox[,2])
    plot.new()
    plot.window(xlim, ylim, asp=asp)
    if(axes) {
      box()
      axis(2, tck=tck)
      axis(1, tck=tck)
    }
  }
  if(PlotPoints){
    points(x$points, col=col, pch=20, ...)
  }
    
  if(PlotVectors){
    if(length(col)>1 && length(col)==nrow(x$points)){
      col=rep(col,rep.int(2,length(col)))
    }
    halfvect=x$vect/2*scalevecs
    if(UseAlpha) halfvect=halfvect*x$alpha
    starts=x$points-halfvect
    stops=x$points+halfvect
    # interleaved=matrix(t(cbind(starts[,1:2],stops[,1:2])),ncol=2,byrow=T)
    segments(starts[,1], starts[,2], stops[,1], stops[,2], col=col, ...)
  }
  invisible()
}

#' Subset points in dotprops object that match given conditions
#' 
#' @details \code{subset} defines either logical or numeric indices, in which
#'   case these are simply applied to the matrices that define the \code{points}, \code{vect} fields of the \code{dotprops} object
#'   etc OR a function (which is called with the 3D points array and returns T/F.
#'   OR an expression 
#'   vector).
#' @param x A dotprops object
#' @param subset A subset of points defined by indices, an expression or a function (see Details)
#' @param ... Additional parameters (currently ignored)
#' @inheritParams subset.neuron
#' @return subsetted dotprops object
#' @method subset dotprops
#' @export
#' @seealso \code{prune.dotprops}, \code{subset.neuron}
#' @examples
#' ## subset using indices ...
#' dp=kcs20[[10]]
#' dp1=subset(dp, 1:50)
#' 
#' # ... or an expression
#' dp2=subset(dp, alpha>0.7)
#' front=subset(dp, points[,'Z']<40)
#' # use a helper function
#' between=function(x, lower, upper) x>=lower & x<=upper
#' middle=middle=subset(dp, between(points[,'Z'], 40, 60))
#' 
#' # plot results in 3D
#' \donttest{
#' plot3d(front, col='red')
#' plot3d(middle, col='green')
#' plot3d(dp, col='blue')
#' }
#' 
#' \dontrun{
#' 
#' ## subset using an selection function
#' s3d=select3d()
#' dp1=subset(dp, s3d(points))
#' # special case of previous version
#' dp2=subset(dp, s3d)
#' # keep the points that were removed from dp2
#' dp2.not=subset(dp, s3d, invert=TRUE)
#' # (another way of doing the same thing)
#' dp2.not=subset(dp, Negate(s3d))
#' stopifnot(all.equal(dp1, dp2))
#' dp2=subset(dp, alpha>0.5 & s3d(pointd))
#' dp3=subset(dp, 1:10)
#' 
#' ## subset each dotprops object in a whole neuronlist
#' plot3d(kcs20)
#' s3d=select3d()
#' kcs20.partial = nlapply(kcs20, subset, s3d)
#' clear3d()
#' plot3d(kcs20.partial, col='red')
#' plot3d(kcs20, col='grey')
#' }
subset.dotprops<-function(x, subset, invert=FALSE, ...){
  e <- substitute(subset)
  r <- eval(e, x, parent.frame())
  if (!is.logical(r) && !is.numeric(r)) {
    # a function that tells us whether a point is in or out
    if(is.function(r)) r=r(x$points)
    else stop("Cannot evaluate subset")
  }
  if(is.logical(r)) {
    r <- r & !is.na(r)
    if(invert) r <- !r
  } else if(is.numeric(r)) {
    if(invert) r <- setdiff(seq_len(nvertices(x)), r)
  } else {
    stop("Subset must evaluate to a logical or numeric index")
  }
  
  x$points=x$points[r,,drop=F]
  x$alpha=x$alpha[r]
  x$vect=x$vect[r,,drop=F]
  if(!is.null(x$labels)) x$labels=x$labels[r]
  x
}

#' prune an object by removing points near (or far) from a target object
#' @export
#' @param x The object to prune. (e.g. \code{dotprops} object, see details)
#' @param target Another object with 3D points that will determine which points 
#'   in x are kept.
#' @param ... Additional arguments for methods (eventually passed to 
#'   \code{prune.default})
#' @seealso \code{\link{prune_strahler}}, \code{\link{spine}},
#'   \code{\link{prune_vertices}}
#' @examples
#' ## prune single neurons
#' \donttest{
#' plot3d(kcs20[[1]],col='blue')
#' plot3d(kcs20[[2]],col='red')
#' }
#' # prune neuron 2 down to points that are close to neuron 1
#' neuron2_close=prune(kcs20[[2]], target=kcs20[[1]], maxdist=10)
#' \donttest{
#' plot3d(neuron2_close, col='cyan', lwd=3)
#' }
#' neuron2_far=prune(kcs20[[2]], target=kcs20[[1]], maxdist=10, keep='far')
#' \donttest{
#' plot3d(neuron2_far, col='magenta', lwd=3)
#' }
#' 
#' ## Prune a neuron with a neuronlist
#' pruned=prune(kcs20[[11]], kcs20[setdiff(1:20, 11)], maxdist=8)
#' \donttest{
#' plot3d(pruned, col='red', lwd=3)
#' plot3d(kcs20[[11]], col='green', lwd=3)
#' plot3d(kcs20,col='grey')
#' }
prune<-function(x, target, ...) UseMethod("prune")

#' @export
#' @details \code{prune.neuron} depends on a more basic function 
#'   \code{\link{prune_vertices}} and is also related to
#'   \code{\link{subset.neuron}}.
#' @rdname prune
#' @seealso \code{\link{subset.neuron}}
#' @family neuron
prune.neuron<-function(x, target, ...){
  indstokeep=NextMethod(return.indices=TRUE)
  indstodrop=setdiff(seq(nrow(x$d)), which(indstokeep))
  prune_vertices(x, indstodrop, ...)
}

#' @export
#' @method prune dotprops
#' @rdname prune
#' @seealso \code{\link{subset.dotprops}}
prune.dotprops<-function(x, target, ...){
  inds=NextMethod(return.indices=TRUE)
  subset(x, inds)
}

#' @export
#' @method prune neuronlist
#' @rdname prune
prune.neuronlist<-function(x, target, ...){
  nlapply(x, prune, target=target, ...)
}

#' @export
#' @method prune default
#' @rdname prune
#' @param maxdist The threshold distance for keeping points
#' @param keep Whether to keep points in x that are near or far from the target
#' @param return.indices Whether to return the indices that pass the test rather
#'   than the 3D object/points (default \code{FALSE})
#' @importFrom nabor knn
prune.default<-function(x, target, maxdist, keep=c("near","far"), 
                        return.indices=FALSE, ...){
  keep=match.arg(keep, c("near", "far"))
  xyzx=xyzmatrix(x)
  nn_dists=drop(knn(xyzmatrix(target), xyzx, k=1)$nn.dists)
  inds=if(keep=="near") nn_dists<=maxdist else nn_dists>maxdist
  if(return.indices) inds else xyzx[inds, , drop=FALSE]
}
jefferis/nat documentation built on Aug. 26, 2018, 12:13 p.m.