Nothing
#' 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 dotprops objects
#' @param x A dotprops object
#' @param y A scalar or 3-vector that will be applied to the dotprops object
#' @return A new dotprops object
#' @export
#' @method * dotprops
#' @rdname dotprops-arithmetic
`*.dotprops` <- function(x,y) {
ly=length(y)
if(!ly%in%c(1,3)) stop("expects a numeric vector of length 1 or 3")
xyzmatrix(x)<-t(t(xyzmatrix(x))*y)
x
}
#' @export
#' @method + dotprops
#' @rdname dotprops-arithmetic
`+.dotprops` <- function(x,y) {
ly=length(y)
if(!ly%in%c(1,3)) stop("expects a numeric vector of length 1 or 3")
xyzmatrix(x)<-t(t(xyzmatrix(x))+y)
x
}
#' @method - dotprops
#' @export
#' @rdname dotprops-arithmetic
`-.dotprops` <- function(x,y) x+(-y)
#' @method / dotprops
#' @export
#' @rdname dotprops-arithmetic
`/.dotprops` <- function(x,y) x*(1/y)
#' @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).
#' \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)
}
#' 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]
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.