Nothing
R/spheroid.R
In unfoldr: Stereological Unfolding for Spheroidal Particles
Defines functions
digitizeProfiles
em.saltykov
coefficientMatrixSpheres
binning1d
planarSection
.checkArgs
drawSpheroidIntersection
cylinders3d
spheroids3d
intersectSystem
verticalSection
coefficientMatrixSpheroids
simPoissonSystem
sectionProfiles
updateIntersections
.getAngle
Documented in
binning1d
coefficientMatrixSpheres
coefficientMatrixSpheroids
cylinders3d
digitizeProfiles
drawSpheroidIntersection
em.saltykov
intersectSystem
planarSection
sectionProfiles
simPoissonSystem
spheroids3d
updateIntersections
verticalSection
## Comment: simulation, intersections and visualization of
## spheroid, sphere and spherocylinder 3D systems
## angle in the section plane relative to z axis in 3D
.getAngle <- function(phi) {
## this is slower!
# return (abs(asin(sin(s))))
if(phi<=pi/2) { phi }
else {
if( phi <= pi) pi-phi
else if(phi<1.5*pi) phi %% pi
else 2*pi-phi
}
}
#' Check intersection
#'
#' Check itersection of either spheres, spheroids or spherocylinders
#'
#' For a given list of spheres, spheroids or spherocylinders the function tests whether an object intersects
#' one of the lateral planes (bottom, top plane included) of the simulation box. The vector returned is of length
#' equal to the number of objects in \code{S} and has entries either \code{1} for an object which is intersected by
#' and \code{0} otherwise.
#'
#' @param S list of spheres, spheroids or spherocylinders, see \code{\link{simPoissonSystem}}
#'
#' @return binary integer vector of length equal to the length of \code{S}
#'
#' @author M. Baaske
#'
#' @rdname updateIntersections
#' @export
updateIntersections <- function(S) {
.Call(C_UpdateIntersections, S, environment())
}
#' Construct section profiles
#'
#' Set up section profiles of spheroids for unfolding
#'
#' The function aggregates the necessary information for trivariate unfolding of spheroids' joint size-shape orientation distribution
#' of type either "\code{prolate}" or "\code{oblate}". The argument \code{size} is a numeric matrix of semi-axis lengths where the first column
#' corresponds to the major semi-axis and the second one to minor semi-axis. The orientation of an ellipse is assumed to be measured as the angle
#' between its major axis and vertical axis of the coordinate system in the intersection plane ('z' axis in 3D). For values in \eqn{[0,2\pi]} these
#' angles are automatically transformed to \eqn{[0,\pi/2]} as required by the unfolding procedure.
#'
#' @param size matrix of lengths of the semi-axes
#' @param alpha angle of section profiles in the plane (see details)
#' @param type name of the spheroid type, either "\code{prolate}" or "\code{oblate}" from which the
#' section profiles are assumed to come from
#'
#' @return The function returns a list which consists of either the longer or shorter
#' semi-axis length named \code{A} of section profiles corresponding to the type of spheroids used before and whose joint
#' joint distribution is to be estimated (by unfolding), the shape factor \code{S} of both semi-axes as the shape factor
#' between \eqn{(0,1]} and the orientation angle \code{alpha}, either of class "\code{prolate}" or "\code{oblate}".
#'
#' @examples
#' # load data set
#' data(data15p)
#'
#' # matrix of semi-axes lengths (major,minor)
#' AC <- data.matrix(data15p[c("A","C")])/1000
#'
#' # selecting the minor semi-axis for prolate type of spheroids:
#' # independent of nomenclature (always named \code{A})
#' sp <- sectionProfiles(AC,unlist(data15p["alpha"]))
#'
#' summary(sp$A) # here minor semi-axis because of prolate
#' summary(sp$S) # shape factor
#' summary(sp$alpha) # angle assumed to be w.r.t. (vertical) 'z' axis
#'
#' @author M. Baaske
#' @rdname sectionProfiles
#' @export
sectionProfiles <- function(size,alpha,type=c("prolate","oblate")) {
type <- match.arg(type)
stopifnot(is.matrix(size))
if(anyNA(size) || any(size<0))
stop("'size' must have non-negative values.")
if(anyNA(alpha) || !is.numeric(alpha) || any(alpha<0))
stop(paste("'alpha' must have non-negative values."))
if(max(alpha) > pi/2)
alpha <- try(sapply(alpha,.getAngle),silent=TRUE)
if(inherits(alpha,"try-error"))
warning("Could not compute 'alpha'. See its returned list element.")
structure(list("A"=if(type=="prolate") size[,2] else size[,1],
"S"=size[,2]/size[,1],
"alpha"=alpha),
class=type)
}
#' Poisson germ-grain process
#'
#' Simulation of Poisson germ-grain processes with either spheres, spheroids or spherocylinders as grains
#'
#' The function can simulate a Poisson germ-grain process according to the parameter \code{theta} within a predefined (3D) box.
#' The positions of the germs follow a uniform distribution according to a Poisson process with mean intensity parameter \code{lam}.
#' The function either randomly generates \code{type="prolate"} or \code{type="oblate"} spheroids, spheres or spherocylinders.
#' The argument \code{size} sets the name of the distribution function for the size/length of the objects, i.e. the major semi-axis
#' lengths in case of spheroids, radii for spheres or the lengths of the main axis of rotation for spherocylinders including the end caps.
#'
#' The following direction (orientation) distributions of the spheroids' major-axis, respectively, cylinders' main axis are available:
#' a uniform distribution ("\code{runifdir}"), distribution ("\code{rbetaiso}") and the "\emph{von Mises-Fisher}" ("\code{rvMisesFisher}")
#' distribution. The two last ones depend on the concentration parameter \code{kappa} which is set as part of the parameter list \code{theta}, see examples below.
#' The direction distributions generate random spherical coordinates \eqn{(\vartheta, \phi)} w.r.t. a fixed main orientation axis \code{mu}
#' with polar angle \eqn{\vartheta\in[0,\pi/2)} and azimuthal angle \eqn{\phi\in[0,2\pi)}. The simulations are always performed within a bounding 3D box which consists of a list specifying the ranges of each dimension corresponding to the lower and upper limits of the box in each direction. If the
#' argument \code{box} contains only a single range, i.e. \code{box=list(c(0,1))}, this limit isassumed for the remaining dimensions which is then simply replicated.
#' The optional argument \code{rjoint} defines a (joint) distribution function which can be any function provided by the user in order to generate
#' the required distributional parameters for the spheroids or spherocylinders. For an in-depth example of usage please see the workflow in 'simSpheroids.R'
#' and 'simCylinders.R'.
#'
#' In addition, the function supports an exact simulation type [2] of the grains. In case of spheroids and spherocylinders setting \code{size="rbinorm"}
#' declares a bivariate size-shape distribution for which the exact simulation is available. More specifically, for a bivariate normal random vector \eqn{[X,Y]}
#' with correlation parameter \eqn{\rho}, the length of the major semi-axis of a spheroid is given by \eqn{a=exp(x)} with a (logit-transformed) shape parameter
#' as \eqn{s=1/(1+exp(-y))} and thus a scaled minor semi-axis length \eqn{c=a*s}. The modification leads to a log-normally distributed length of the
#' major semi-axis. Consequently, in case of spherocylinders, the log-normally distributed length is \eqn{len=h+2*r} where \code{h} is the height and
#' \eqn{r=len/2*s} the radius. The main direction \code{u} of a spheroid or spherocylinder is determined by the major axis independent of size and shape.
#' Further, the following univariate distributions of the major semi-axis \code{a}, respectively, length \code{len} and shape \code{s} are available:
#' '\code{rbeta}', '\code{rgamma}', '\code{rlnorm}' and '\code{runif}'. One can also use '\code{const}' for simulations with constant lengths or shapes.
#' Note that only simulations with size distributions '\code{rbinorm}' or '\code{rlnorm}' can use the exact type of simulation.
#'
#' For spheres any distribution of the radii can be specified as a name of a user-defined function in the argument \code{size} as long as the formal named
#' function parameters match the actual names of the parameters exactly as defined in \code{theta}. Besides this, all other distributions given above are
#' also available. Using '\code{const}' simulates spheres of constant radii.
#'
#' The argument \code{pl>=0} denotes both the print level of intermediate output and by the same time the type of the return value. If \code{pl=10},
#' then an abbreviated list of spheroids or spheres is returned to speed up computation. Note that, the current implementation does not include routines
#' for unfolding the joint size-shape-orientation distribution of spherocylinders so far.
#'
#'
#' @param theta list of simulation parameters which must consist of elements: \code{size}, \code{shape} and \code{orientation}
#' @param lam mean number of objects per unit volume
#' @param size name of the size distribution function
#' @param shape name of the shape distribution function
#' @param orientation name of direction distribution function
#' @param type type of grain, either "\code{prolate}" or "\code{oblate}", "\code{spheres}", "\code{cylinders}"
#' @param rjoint user-defined function, which specifies the (joint) distribution of the size, shape and orientation
#' @param box simulation box
#' @param mu main orientation axis, \code{mu=c(0,0,1)} (default)
#' @param dz distance of the intersecting plane to the origin
#' @param n normal vector defining the intersecting plane
#' @param intersect options for type of return values: "\code{full}" for the simulated system together with section profiles as two lists named \code{S} and \code{sp} respectively,
#' choose "\code{only}" for section profiles only, "\code{original}" for the 3D system only and "\code{digi}" for a (binary) integer matrix \code{W} as a discretized
#' version of section profiles whose resolution depends on the chosen lattice constant \code{delta}
#' @param delta lattice constant for discretization, set to \code{0.01} (default)
#' @param intern logical, \code{FALSE} (default), whether to return only section profiles with centers inside the simulation window
#' @param perfect logical, \code{FALSE} (default), whether to simulate exactly (also called perfect)
#' @param pl integer, print level and return value type, see details
#' @param label character, a label set to each generated object, set to '\code{N}' (default)
#'
#' @return list of 3D objects depending on the chosen return type defined by the argument \code{intersect}
#'
#' @examples
#' # intensity parameter
#' lam <- 100
#'
#' # simulation bounding box
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # log normal size distribution with a constant shape factor and
#' # concentration parameter (\code{kappa=1}) for the orientation, see reference [1]
#' theta <- list("size"=list("meanlog"=-2.5,"sdlog"=0.5),
#' "shape"=list("s"=0.5),
#' "orientation"=list("kappa"=1))
#'
#' S <- simPoissonSystem(theta,lam,size="rlnorm",box=box,type="oblate",pl=1)
#' length(S)
#'
#' @references
#' \itemize{
#' \item{} {Ohser, J. and Schladitz, K. 3D images of materials structures Wiley-VCH, 2009}
#' \item{} {C. Lantu\eqn{\acute{\textrm{e}}}joul. Geostatistical simulation. Models and algorithms.
#' Springer, Berlin, 2002. Zbl 0990.86007}
#' }
#' @author M. Baaske
#' @rdname simPoissonSystem
#' @export
simPoissonSystem <- function(theta, lam, size="const", shape="const", orientation="rbetaiso",
type=c("prolate","oblate","spheres","cylinders"), rjoint=NULL,
box=list(c(0,1)), mu=c(0,0,1), dz=0, n=c(0,1,0),
intersect=c("original","full","only","digi"),
delta=0.01, intern=FALSE, perfect=FALSE, pl=0, label="N")
{
it <- pmatch(type,c("prolate","oblate","spheres","cylinders"))
if(length(it)==0 || is.na(it))
stop(paste("`type` is one of: ", paste0("`",c("prolate","oblate","sphere","cylinder"),"`",collapse=","),sep=""))
if(!is.numeric(lam) || !(lam>0) )
stop("Expected 'lam' as non-negative numeric argument")
if(length(box)==0 || !is.list(box))
stop("Expected argument 'box' as list type.")
if(length(box)==1)
box <- rep(box[1],3)
else if(length(box)!=3)
stop("Simulation box has wrong dimensions.")
names(box) <- c("xrange","yrange","zrange")
if(dz < min( box[[which(n == 1)]]) || dz > max(box[[which(n == 1)]]))
stop("The plane has to intersect the box.")
# digitization (resolution) factor
stopifnot(is.numeric(delta))
# spheroid type
type <- match.arg(type)
intersect <- match.arg(intersect)
cond <-
if(!is.null(rjoint))
{
if(!is.function(rjoint))
stop("Unknown multivariate distribution function.")
#largs <- theta[-(which(it==1))]
if(!is.list(theta) || length(theta) == 0L)
stop("Expected 'theta' as list of named arguments matching exactly with formal arguments of function `rjoint`.")
it <- match(names(theta),names(formals(rjoint)))
funname <- as.character(substitute(rjoint))
if(length(it) == 0L || anyNA(it))
stop(paste("Arguments must match formal arguments of function `",funname,"`.",sep=""))
# check function
funret <- try(do.call(rjoint,theta))
if(inherits(funret,"try-error"))
stop(paste("Error in user defined function `",funname,"`.",sep=""))
fargs <-
if(type %in% c("prolate","oblate")) {
if(!is.list(funret))
stop("Expected a list as return type in user defined function.")
c("a","b","c","u","shape","theta","phi")
} else if(type %in% c("cylinders")) {
if(!is.list(funret))
stop("Expected a list as return type in user defined function.")
c("h","r","u","theta","phi")
} else if(type %in% c("spheres")) {
if(!is.numeric(funret))
stop("Expected a numeric vector as return type in user defined function.")
c("r")
} else { stop("Undefined type!")}
if(any(!(fargs %in% names(funret)))){
stop("Argument names of return value list does not match required arguments.")
}
list("type"=type,"rdist"=funname,"box"=box,
"lam"=lam,"pl"=pl,"mu"=mu,"rho"= environment(),"label"=label,
"dz"=dz, "nsect"=n, "delta"=delta,"intern"=as.integer(intern),
"perfect"=as.integer(perfect), "intersect"=intersect)
} else {
# set defaults if not provided
if(!is.list(theta))
stop("Expected 'theta' as list of named arguments `size`, `shape`, `orientation`.")
nms <- c("size","shape","orientation")
it <- match(names(theta), nms)
if(length(it)>0L && anyNA(it))
names(theta)[it] <- nms[it]
if(shape == "const")
{
# set defaults
if(!is.list(theta$shape) || length(theta$shape) == 0L)
theta$shape <- list(1,1)
} else if(exists(shape, mode="function")) {
# check arguments of supported distribution functions
fargs <- names(formals(shape))
if(shape %in% c("rbeta","rgamma","runif"))
fargs <- fargs[-1]
if(!is.list(theta$shape) || length(theta$shape) == 0L)
stop("`shape` parameters must be given in `theta`.")
it <- match(names(theta$shape),fargs)
if(length(it)==0 || anyNA(it))
stop(paste0("Arguments of 'shape' must match formal arguments of function ",shape,sep=""))
} else {
stop(paste0("Undefined `", shape, "` distribution function."))
}
if(!is.list(theta$size) || length(theta$size)==0L)
stop("Arguments for size distribution must be given in `theta`.")
if(size == "rbinorm") {
fargs <- c("mx","my","sdx","sdy","rho","kappa")
it <- match(names(theta$size),fargs)
if(length(it)==0 || anyNA(it))
stop(paste("Arguments of 'size' must match formal arguments: ", paste0("`",fargs,"`",collapse=",")))
} else if(exists(size, mode="function")) {
# check arguments of supported distribution functions
fargs <- names(formals(size))
if(size %in% c("rlnorm","rbeta","rgamma","runif"))
fargs <- fargs[-1]
it <- match(names(theta$size),fargs)
if(length(it)==0 || anyNA(it))
stop(paste0("Arguments of 'size' must match formal arguments of function ",size,sep=""))
} else if(size == "const"){
# set defaults
if(!is.list(theta$size) || length(theta$size) == 0L)
theta$size <- list(1,1)
} else {
stop(paste0("Undefined `", size, "` distribution function."))
}
it <- pmatch(orientation,c("runifdir","rbetaiso","rvMisesFisher","const"))
if(is.na(it) && !exists(orientation, mode="function"))
stop("Undefined distribution function for the orientation/direction.")
# set defaults
if(!is.list(theta$orientation) || length(theta$orientation) == 0L)
theta$orientation <- list("kappa"=1)
list("type"=type, "lam"=lam,
"rdist"=list("size"=size, "shape"=shape, "orientation"=orientation),
"box"=box,"pl"=pl,"mu"=mu,"rho"= environment(), "dz"=dz, "nsect"=n, "delta"=delta,
"intern"=as.integer(intern),"label"=label, "perfect"=as.integer(perfect),
"intersect"=intersect)
}
.Call(C_PoissonSystem, theta, cond)
}
#' Coefficients for trivariate unfolding
#'
#' Compute coefficients of the discretized integral equation for unfolding
#'
#' In order to apply the Expectation Maximization (EM) algorithm to the stereological unfolding procedure for the joint
#' size-shape-orientation distribution in 3D one first has to compute the coefficients of a discretized integral equation.
#' This step is the most time consuming part of unfolding and therefore has been separated in its own function and can be called
#' separately if needed. The number of bin classes for the size, shape and orientation do not need to be the same, but the
#' given class limits are also used for binning the spatial values. One can define the number of cpu cores by the global option
#' \code{par.unfoldr} or passing the number of cores \code{nCores} directly to the function.
#'
#' @param breaks list of bin vectors, see \code{\link{setbreaks}}
#' @param stype type of spheroid, either "\code{prolate}" or "\code{oblate}"
#' @param check logical, whether to run some input checks
#' @param nCores number of cpu cores used to compute the coefficients
#'
#' @return numeric 6D array of coefficients
#'
#' @example inst/examples/coeffarray.R
#'
#' @author M. Baaske
#' @rdname coefficientMatrixSpheroids
#' @export
coefficientMatrixSpheroids <- function(breaks, stype=c("prolate","oblate"),
check=TRUE,nCores=getOption("par.unfoldr",1))
{
stype <- match.arg(stype)
it <- match(names(breaks), c("size","angle","shape"))
if (length(it)==0 || anyNA(it))
stop("Expected 'breaks' as named list of: 'size','angle','shape' ")
if (is.unsorted(breaks$size) || is.unsorted(breaks$angle) || is.unsorted(breaks$shape))
stop("'breaks' list must contain non-decreasingly sorted classes")
if(check) {
if(any(breaks$size<0))
stop("Breaks vector 'size' must have non-negative values.")
if(min(breaks$angle)<0 || max(breaks$angle)>pi/2)
stop(paste0("Breaks vector 'angle' must have values between zero and ",pi/2))
if(min(breaks$shape)<0 || max(breaks$shape)>1)
stop("Breaks vector 'shape' must have values between 0 and 1.")
}
.Call(C_CoefficientMatrixSpheroids,
breaks$size,breaks$angle,breaks$shape,
breaks$size,breaks$angle,breaks$shape, list(stype,nCores))
}
#' Vertical sections
#'
#' Compute vertical section profiles of a spheroid system
#'
#' The function intersects a spheroid system by a plane defined by the normal vector \code{n} either
#' equal to \code{c(0,1,0)} (default) or \code{c(1,0,0)}, which is called a vertical section. Depending on
#' the type of spheroid (either "\code{prolate} or "\code{oblate}") the returned semi-axis lengths are those
#' corresponding to the minor semi-axis or, respectively, major semi-axis in the way these are required for unfolding.
#'
#' @param S list of spheroids, see \code{\link{simPoissonSystem}}
#' @param d distance of the intersecting plane from the origin of the box
#' @param n normal vector which defines the interecting vertical plane
#' @param intern logical, \code{FALSE} (default), return all section profiles otherwise
#' only those which have their centers inside the correspondig intersection window
#'
#' @return list of sizes \code{A}, shape factors \code{S} and (vertical) angles \code{alpha}
#' of section profiles in the plane w.r.t the 'z' axis between \eqn{[0,\pi/2]}.
#'
#'
#' @examples
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # (exact) bivariate size-shape (isotropic) orientation distribution (spheroids)
#' theta <- list("size"=list("mx"=-2.5,"my"=0.5, "sdx"=0.35,"sdy"=0.25,"rho"=0.15),
#' "orientation"=list("kappa"=1))
#'
#' S <- simPoissonSystem(theta,lam=100,size="rbinorm",box=box,
#' type="prolate",perfect=TRUE,pl=1)
#'
#' sp <- verticalSection(S,d=2.5,n=c(0,1,0),intern=TRUE)
#' summary(sp$alpha)
#'
#' @author M. Baaske
#' @rdname verticalSection
#' @export
verticalSection <- function(S,d,n=c(0,1,0),intern=FALSE) {
if(!(class(S) %in% c("prolate","oblate")))
stop("Only spheroids of class 'prolate' or 'oblate' can be used.")
stopifnot(is.numeric(d))
stopifnot(is.logical(intern))
if(sum(n) != 1 || which(n == 1) > 2L)
stop("Normal vector is one of c(0,1,0) or c(1,0,0).")
ss <- .Call(C_IntersectPoissonSystem, S, list("nsect"=n,"dz"=d,"intern"=intern,"pl"=10), environment())
A <- if(class(S)=="prolate")
sapply(ss,"[[",2)
else sapply(ss,"[[",1)
## convert angle 'alpha' in the intersecting plane
## which is always between [0,pi/2] and w.r.t 'z' axis
alpha <- sapply(ss,"[[",4)
alpha <- try(0.5*pi-sapply(alpha,.getAngle),silent=TRUE) # alpha in [0,pi/2]
if(inherits(alpha,"try-error") || !is.numeric(alpha) || anyNA(alpha) )
warning("Could not compute angle 'alpha'. See the returned list element.")
structure(list("A"=A,
"S"=sapply(ss,"[[",3),
"alpha"=alpha), # w.r.t z axis in 3D
"win"=attr(ss,"win"),"plane"=attr(ss,"plane"),
class=class(S) )
}
#' Intersection in 3D
#'
#' Intersect a system of spheres, spheroids or spherocylinders by a plane
#'
#' The function intersects a given (Poisson) system made of spheres, spheroids or cylinders as grains by a plane defined by a
#' normal vector, e.g. \code{n=c(0,1,0)}, perpendicular to one of the bounding planes of the simulation box.
#' For a print level \code{pl>=0} some verbose output is given. Also it sets the type of return value. In case of spheroid
#' intersections, setting \code{pl=10}, leads to a short version of the full specification return list of section profiles
#' with elements named \code{A} (major semi-axis), \code{C} (minor semi-axis), \code{S} (the shape factor as the ratio of these two)
#' and \code{phi} as the angle in the intersecting plane between \eqn{[0,2\pi]} w.r.t. the 'x' axis. Otherwise additional components
#' are returned such as the ellipse`s rotation matrix, also named \code{A}, the center point \code{center} and a constant \code{type=10}
#' (defining the object of full ellipses among other types of intersection objects) of the section profiles. For sphere intersections only
#' a numeric vector of disc radii are returned as a short version of return values.
#'
#' @param S list of spheres, spheroids or spherocylinders, see \code{\link{simPoissonSystem}}
#' @param d distance of the the box-aligned intersecting plane from the origin
#' @param n normal vector which defines the intersecting plane
#' @param intern logical, \code{FALSE} (default), return all section profiles otherwise
#' only those which have their centers inside the corresponding intersection window
#' (if the intersected system had been simulated using exact simulation this makes sense)
#' @param pl integer, \code{pl=0} (default), for no verbose output and otherwise a full specification list of
#' section profiles in case of sphere and spheroid intersections
#'
#' @return For spheroid intersections the function returns a list of size, shape and angle of section profiles or a short version
#' of it; for sphere intersections either radii or lists containing the centers of discs and the object number.
#'
#'
#' @examples
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # constant size-shape orientation distribution (spheroids)
#' theta <- list("size"=list(0.1),"shape"=list(0.5), "orientation"=list("kappa"=10))
#'
#' S <- simPoissonSystem(theta,lam=100,box=box,type="prolate",pl=1)
#'
#' # return short version of section profiles
#' sp <- intersectSystem(S, 2.5, pl=10)
#'
#' @author M. Baaske
#' @rdname intersectSystem
#' @export
intersectSystem <- function(S, d, n=c(0,1,0), intern=FALSE, pl=0) {
if(!(class(S) %in% c("prolate","oblate","spheres","cylinders")))
stop("Unknown object class.")
stopifnot(is.numeric(d))
stopifnot(is.numeric(pl))
stopifnot(is.logical(intern))
if(sum(n)>1)
stop("Normal vector is like c(0,1,0). ")
box <- attr(S,"box")
if(is.null(box))
stop("'box' must be given as an attribute.")
if(d < min( box[[which(n == 1)]]) || d > max(box[[which(n == 1)]]))
stop("The plane has to intersect the box.")
.Call(C_IntersectPoissonSystem,S, list("nsect"=n,"dz"=d,"intern"=as.integer(intern),"pl"=as.integer(pl)), environment())
}
#' Spheroid system 3D
#'
#' Draw a spheroid system in 3D
#'
#' The function requires the package '\code{rgl}' to draw spheroids. For an example
#' please see the file 'simSpheroids.R'.
#'
#' @param S list of spheroids, see \code{\link{simPoissonSystem}}
#' @param box simulation box
#' @param draw.axes logical, whether to show the axes
#' @param draw.box logical, whether to show the bounding box
#' @param draw.bg logical, whether to show a gray background box correpsonding to the simulation box
#' @param bg.col background color used showw the box background
#' @param clipping logical, whether to clip all lateral planes of the simulation box
#' @param ... further material properties passed to 3D plotting functions
#'
#' @return NULL
#' @author M. Baaske
#' @rdname spheroids3d
#' @export
spheroids3d <- function(S, box, draw.axes=FALSE, draw.box=TRUE, draw.bg=TRUE,
bg.col="white", clipping=FALSE, ...)
{
if (!requireNamespace("rgl", quietly=TRUE))
stop("Please install 'rgl' package from CRAN repositories before running this function.")
ellipsoid3d <- function(rx=1,ry=1,rz=1,n=50,ctr=c(0,0,0), qmesh=FALSE,trans = rgl::par3d("userMatrix"),...) {
if (missing(trans) && !rgl::rgl.cur())
trans <- diag(4)
degvec <- seq(0,2*pi,length=n)
ecoord2 <- function(p) {
c(rx*cos(p[1])*sin(p[2]),ry*sin(p[1])*sin(p[2]),rz*cos(p[2]))
}
v <- apply(expand.grid(degvec,degvec),1,ecoord2)
if (qmesh)
v <- rbind(v,rep(1,ncol(v))) ## homogeneous
e <- expand.grid(1:(n-1),1:n)
i1 <- apply(e,1,function(z)z[1]+n*(z[2]-1))
i2 <- i1+1
i3 <- (i1+n-1) %% n^2 + 1
i4 <- (i2+n-1) %% n^2 + 1
i <- rbind(i1,i2,i4,i3)
if (!qmesh)
rgl::quads3d(v[1,i],v[2,i],v[3,i],...)
else
return(rgl::rotate3d(rgl::qmesh3d(v,i,material=...),matrix=trans))
}
sphere <- ellipsoid3d(qmesh=TRUE,trans=diag(4))
spheroid3d <- function (x=0,y=0,z=0, a=1, b=1,c=1, rotM, subdivide = 3, smooth = TRUE){
result <- rgl::scale3d(sphere, a,b,c)
result <- rgl::rotate3d(result,matrix=rotM)
result <- rgl::translate3d(result, x,y,z)
invisible(result)
}
N <- length(S)
ll <- lapply(S, function(x)
spheroid3d(x$center[1],x$center[2],x$center[3],
x$acb[1],x$acb[2],x$acb[3],rotM=x$rotM))
rgl::shapelist3d(ll,...)
x <- box$xrange[2]
y <- box$yrange[2]
z <- box$zrange[2]
## draw gray box
if(draw.bg) {
c3d.origin <- rgl::translate3d(rgl::scale3d(rgl::cube3d(col=bg.col, alpha=0.1),x/2,y/2,z/2),x/2,y/2,z/2)
rgl::shade3d(c3d.origin)
}
if(clipping) {
rgl::clipplanes3d(-1,0,0,x)
rgl::clipplanes3d(0,-1,0,y)
rgl::clipplanes3d(0,0,-1,z)
rgl::clipplanes3d(1,0,0,0)
rgl::clipplanes3d(0,1,0,0)
rgl::clipplanes3d(0,0,1,0)
}
if(draw.axes) {
rgl::axes3d(c('x','y','z'), pos=c(0,0,0))
rgl::title3d('','','x','y','z')
}
## draw box
if(draw.box) {
rgl::axes3d(edges = "bbox",labels=TRUE,tick=FALSE,box=TRUE,nticks=0,
expand=1.0,xlen=0,xunit=0,ylen=0,yunit=0,zlen=0,zunit=0)
}
}
#' Cylinder system 3D
#'
#' Draw spherocylinders in 3D
#'
#' The function requires the package '\code{rgl}' to draw spherocylinders into a 3D '\code{rgl}' image.
#' For an example please see the file 'simCylinders.R'.
#'
#' @param S list of cylinders, see \code{\link{simPoissonSystem}}
#' @param box simulation box
#' @param draw.axes logical, whether to show the axes
#' @param draw.box logical, whether to show the bounding box
#' @param clipping logical, whether to clip all lateral planes of the simulation box
#' @param ... further material properties passed to 3d plotting functions
#'
#' @return NULL
#' @author M. Baaske
#' @rdname cylinders3d
#' @export
cylinders3d <- function(S, box, draw.axes=FALSE, draw.box=TRUE, clipping=FALSE,...)
{
if (!requireNamespace("rgl", quietly=TRUE))
stop("Please install 'rgl' package from CRAN repositories before running this function.")
cylinder <- function(m, radius=1, h=1, rotM=diag(3), u=c(0,0,1)) {
cyl <- rgl::cylinder3d(rbind(c(0,0,0),c(0,0,1)), radius=radius,
e1=cbind(0, 0, 1), e2=cbind(1, 0, 0), sides=25,closed=-2 )
result <-rgl::scale3d(cyl,1,1,h)
result <- rgl::rotate3d(result,matrix=rotM)
## translate to midpoint of cylinder, center of mass
m <- m-0.5*h*u
result <- rgl::translate3d(result, m[1],m[2],m[3])
invisible(result)
}
args <- list(...)
ok <- sapply(S, function(x) (!is.null(x$h) && x$h > 0.0) ) # only true cylinders not 'spheres'
cyls <- lapply(S[ok], function(x) { cylinder(x$center, x$r, x$h, x$rotM, x$u) })
rgl::shapelist3d(cyls,...)
if("col" %in% names(args)) {
cols <- rep(rep(args$col,length.out=length(S)),each=2)
args$col <- NULL
} else cols <- "black"
# spheres
is.sphere <- sapply(S,function(x) !(is.null(x$h) || x$h == 0.0) )
Xc <- do.call(rbind,lapply(S[ok],function(x) rbind( c(x$origin0,x$r),c(x$origin1,x$r))))
rgl::spheres3d(Xc,radius=Xc[,4],col=cols,unlist(args))
if(!all(ok)) {
Xc <- do.call(rbind,lapply(S[!ok],function(x) rbind(c(x$center,x$r))))
rgl::spheres3d(Xc,radius=Xc[,4],col="darkgray",unlist(args))
}
x <- box$xrange[2]; y <- box$yrange[2]; z <- box$zrange[2]
c3d.origin <- rgl::translate3d(rgl::scale3d(rgl::cube3d(col="darkgray", alpha=0.1),x/2,y/2,z/2),x/2,y/2,z/2)
rgl::shade3d(c3d.origin)
if(clipping) {
rgl::clipplanes3d(-1,0,0,box$xrange[1])
rgl::clipplanes3d(0,-1,0,box$yrange[1])
rgl::clipplanes3d(0,0,-1,box$zrange[1])
rgl::clipplanes3d(1,0,0,box$xrange[2])
rgl::clipplanes3d(0,1,0,box$yrange[2])
rgl::clipplanes3d(0,0,1,box$zrange[2])
}
if(draw.axes) {
rgl::axes3d(c('x','y','z'), pos=c(0,0,0))
rgl::title3d('','','x','y','z')
}
if(draw.box) {
rgl::axes3d(edges = "bbox",labels=TRUE,tick=FALSE,box=TRUE,nticks=0,
expand=1.0,xlen=0,xunit=0,ylen=0,yunit=0,zlen=0,zunit=0)
}
}
#' Spheroid intersections 3D
#'
#' Draw section profiles of spheroids in 3D
#'
#' The function requires the package '\code{rgl}' for drawing section profiles (ellipses)
#' into a 3D '\code{rgl}' image. For a full example please see the file 'simSpheroids.R'.
#'
#' @param E a list of spheroid intersections, see \code{\link{intersectSystem}}
#' @param n the normal vector of the intersecting plane
#' @param np number of points for a polygon approximation of the ellipses
#'
#' @return NULL
#'
#' @author M. Baaske
#' @rdname drawSpheroidIntersection
#' @export
drawSpheroidIntersection <- function(E, n=c(0,1,0), np=25) {
ind <- which(n == 0)
.pointsOnEllipse <- function(E,t) {
E$phi <- E$phi + 0.5*pi
xt <- E$center[1] + E$ab[1]*cos(t)*cos(E$phi)-E$ab[2]*sin(t)*sin(E$phi)
zt <- E$center[2] + E$ab[1]*cos(t)*sin(E$phi)+E$ab[2]*sin(t)*cos(E$phi)
yt <- rep(0,length(t))
m <- matrix(0,nrow=length(xt),ncol=3)
m[,ind[1]] <- xt
m[,ind[2]] <- zt
m
}
.plotEllipse <- function(x) {
M <- .pointsOnEllipse(x,t)
rgl::polygon3d(M[,1],M[,2],M[,3],fill=TRUE,coords=ind)
}
s <- 2*pi/np
t <- seq(from=0,to=2*pi,by=s)
invisible(lapply(E,function(x) .plotEllipse(x) ))
}
.checkArgs <- function(optlist, options) {
if(is.null(names(options)))
stop("Options should be a list of named arguments.")
if (!is.list(options) || "" %in% names(options))
stop("Argument ",as.character(substitute(options)), " must be a list of named (character) elents.")
optnames <- (names(options) %in% names(optlist))
if (!all(optnames)) {
unames <- as.list(names(options)[!(optnames)])
stop(paste(c("Unknown arguments in ",as.character(substitute(options))," : ",do.call("paste", c(unames, sep = ", "))), collapse=" "))
}
return (0)
}
#' Sphere planar sections
#'
#' Planar intersection of spheres
#'
#' The function computes the planar intersection of a sphere system, i.e. an intersection with the plane whose normal
#' vector is given by \code{c(0,0,1)} and returns the diameters of the resulting discs.
#'
#' @param S list of spheres of class \code{sphere}, see \code{\link{simPoissonSystem}}
#' @param d distance of the (planar) xy-plane to the origin
#' @param intern logical, \code{FALSE} (default), whether to return all discs or
#' only those which have their centers inside the intersecting window
#' @param pl print level, \code{pl>0} for some verbose output
#'
#' @return numeric vector of disc diameters
#'
#' @examples
#' lam <- 100
#'
#' # parameter rlnorm distribution (radii)
#' theta <- list("size"=list("meanlog"=-2.5,"sdlog"=0.5))
#'
#' # simulation bounding box
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # simulate only 3D system
#' S <- simPoissonSystem(theta,lam,size="rlnorm",box=box,type="spheres",
#' intersect="original", pl=1)
#'
#' # return only objects whose centers are within
#' # the intersection window
#' sp <- planarSection(S,d=2.5,intern=TRUE,pl=1)
#'
#' # histogram of diameters
#' hist(sp)
#' summary(sp)
#'
#' # distribution of radii
#' mean(log(sp/2))
#' sd(log(sp/2))
#'
#'
#' @author M. Baaske
#' @rdname planarSection
#' @export
planarSection <- function(S, d, intern=FALSE, pl=0) {
if(!(c("spheres") %in% class(S) ))
stop("Expected spheres as a list.")
if(!is.numeric(d) || d<0)
stop("'d' must be numeric/positive.")
stopifnot(is.logical(intern))
sp <- .Call(C_IntersectPoissonSystem,S,list("nsect"=c(0,0,1),"dz"=d,"intern"=intern,"pl"=pl),environment())
if(is.list(sp)) # full list of section profiles
return (sapply(sp,function(x) 2.0*x$r))
else return (2*sp) # only diameters are returned
}
#' Binning numeric values
#'
#' Vector of count data
#'
#' The function provides basic binning (grouping) of numeric values into
#' classes defined by the breaks vector \code{bin}. The values are binned
#' according to \eqn{bin[i[j]]<x[j]\leq bin[i[j]+1]} for intervals \eqn{i=1,...,N-1}
#' and \code{length(bin)=N} of values \eqn{x[j]}, \eqn{j=1,...,|x|}. If \eqn{x[j] > bin[N]} or \eqn{x[j] < bin[1]} then \eqn{x[j]}
#' is not counted at all.
#'
#' @param x numeric values to be binned
#' @param bin non-decreasingly sorted breaks vector
#' @param na.rm logical, default \code{FALSE}, whether to remove missing values, including \code{NaN} in \code{x}
#'
#' @return Vector of count data
#'
#' @examples
#' x <- runif(100,0,1)
#' bin <- seq(0,1,by=0.1)
#' binning1d(x,bin)
#'
#' @author M. Baaske
#' @rdname binning1d
#' @export
binning1d <- function(x,bin, na.rm = FALSE) {
if (anyNA(x)) {
if(na.rm) x <- x[which (!is.na(x))]
else stop("Vectors contains missing values or NAs.")
}
if(anyNA(bin) )
stop("'bin' vector contains missing values or NAs.")
if (is.unsorted(bin))
stop("'bin' must be sorted non-decreasingly")
.Call(C_Binning1d,x,bin)
}
#' Coefficients for Expectation Maximization algorithm
#'
#' Matrix of coefficients for Wicksell's corpuscle problem
#'
#' The function computes the matrix of coefficients for solving the discretized integral equation of
#' the Wicksell's corpuscle problem.
#'
#' @param bin non-decreasing vector of class limits
#' @return array of coefficients
#'
#' @references
#' Ohser, J. and Muecklich, F. Statistical analysis of microstructures in materials science J. Wiley & Sons, 2000
#'
#' @author M. Baaske
#' @rdname coefficientMatrixSpheres
#' @export
coefficientMatrixSpheres <- function(bin) {
## count of bin classes
n <- length(bin)-1
if (n<=0)
stop("'bin' has length zero")
if (anyNA(bin))
stop("Vectors contain NA values")
if (is.unsorted(bin))
stop("'bin' must be sorted non-decreasingly")
if(any(bin<0))
stop("Breaks vector must have non-negative values.")
p <- .C(C_em_saltykov_p, as.integer(n),as.numeric(bin),
p=as.numeric(matrix(0,n,n)))$p
dim(p) <- c(n,n)
return (p)
}
#' Expectation Maximization (EM) algorithm
#'
#' Estimation of the empirical sphere diameter distribution
#'
#' The function performs the EM algorithm, see reference below.
#'
#' @param y vector of observed absolute frequencies of circle diameters
#' @param bin non-decreasing vector of class limits
#' @param maxIt maximum number of iterations to be used
#' @return vector of count data of absolute frequenties of sphere diameters
#'
#' @example inst/examples/simSpheres.R
#'
#' @references
#' Ohser, J. and Muecklich, F. Statistical analysis of microstructures in materials science J. Wiley & Sons, 2000
#'
#' @author M. Baaske
#' @rdname em.saltykov
#' @export
em.saltykov <- function(y,bin,maxIt=32) {
if (length(y)==0) stop("input array 'y' has length zero")
if (anyNA(y) || anyNA(bin) )
stop("Vectors contain NA values")
if (is.unsorted(bin))
stop("'bin' must be sorted non-decreasingly")
n <- length(y)
p <- .C(C_em_saltykov_p, as.integer(n),as.numeric(bin),
p=as.vector(matrix(0,n,n)))$p
theta <- y+1.0e-6
.C(C_em_saltykov,as.integer(n), as.integer(maxIt),
as.numeric(p),as.numeric(y),theta=as.numeric(theta))$theta
}
#' Digitization
#'
#' Digitize 2D section profiles
#'
#' A list of section profiles, e.g. discs, ellipses or segments from intersected spherocylinders, is digitized according
#' to a given resolution according to the lattice constant \code{delta} such that the result is an integer matrix
#' which can be interpretated as a binary image. An intersection window can be either provided by the user w.r.t. \eqn{[l,u]^2}
#' where \code{l,u} are lower, respectively, upper bounds of corresponding to a simulation box, or, it is taken from the section
#' profiles \code{sp} which stores the intersection window as an attribute.
#'
#' @param sp list of section profiles, see \code{\link{intersectSystem}}
#' @param delta lattice constant for discretization of section profiles
#' @param win list of length two, the intersection window, default \code{NULL}
#'
#' @examples
#' # simulation box
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#' # (exact) bivariate size-shape (isotropic) orientation distribution (spheroids)
#' theta <- list("size"=list("mx"=-2.5,"my"=0.5, "sdx"=0.35,"sdy"=0.25,"rho"=0.15),
#' "orientation"=list("kappa"=1))
#'
#' # return only 3D system
#' S <- simPoissonSystem(theta,lam=100,size="rbinorm",box=box,type="prolate",
#' intersect="original",n=c(0,1,0),mu=c(0,0,1),perfect=TRUE,pl=1)
#'
#' # vertical intersection w.r.t. 'mu' (z axis, see above)
#' sp <- intersectSystem(S, 2.5)
#'
#' # show intersecting window
#' win <- attr(sp,"win")
#'
#' # digitize (could also pass some 'win' as an argument)
#' W <- digitizeProfiles(sp, delta=0.01, win = NULL)
#' image(1:nrow(W),1:ncol(W),W,col=gray(1:0))
#'
#' @return binary (integer) matrix as an image
#'
#' @author M. Baaske
#' @rdname digitizeProfiles
#' @export
digitizeProfiles <- function(sp, delta = 0.01, win = NULL ) {
stopifnot(is.numeric(delta))
if(!is.null(win)){
if(!is.list(win) || length(win) == 0L )
stop("Expected argument 'win' as list type.")
if(length(win) == 1L)
win <- rep(win[1],2)
else if(length(win) != 2L)
stop("Digitization window 'win' must be two-dimensional.")
}
if(delta < 0 || delta > 1.0)
stop("'delta' should be positive and equal or less than 1.")
.Call(C_DigitizeProfiles, sp, delta, win, environment())
}
Try the unfoldr package in your browser
Any scripts or data that you put into this service are public.
unfoldr documentation built on Sept. 25, 2021, 3:01 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions digitizeProfiles em.saltykov coefficientMatrixSpheres binning1d planarSection .checkArgs drawSpheroidIntersection cylinders3d spheroids3d intersectSystem verticalSection coefficientMatrixSpheroids simPoissonSystem sectionProfiles updateIntersections .getAngle
Documented in binning1d coefficientMatrixSpheres coefficientMatrixSpheroids cylinders3d digitizeProfiles drawSpheroidIntersection em.saltykov intersectSystem planarSection sectionProfiles simPoissonSystem spheroids3d updateIntersections verticalSection
## Comment: simulation, intersections and visualization of
## spheroid, sphere and spherocylinder 3D systems
## angle in the section plane relative to z axis in 3D
.getAngle <- function(phi) {
## this is slower!
# return (abs(asin(sin(s))))
if(phi<=pi/2) { phi }
else {
if( phi <= pi) pi-phi
else if(phi<1.5*pi) phi %% pi
else 2*pi-phi
}
}
#' Check intersection
#'
#' Check itersection of either spheres, spheroids or spherocylinders
#'
#' For a given list of spheres, spheroids or spherocylinders the function tests whether an object intersects
#' one of the lateral planes (bottom, top plane included) of the simulation box. The vector returned is of length
#' equal to the number of objects in \code{S} and has entries either \code{1} for an object which is intersected by
#' and \code{0} otherwise.
#'
#' @param S list of spheres, spheroids or spherocylinders, see \code{\link{simPoissonSystem}}
#'
#' @return binary integer vector of length equal to the length of \code{S}
#'
#' @author M. Baaske
#'
#' @rdname updateIntersections
#' @export
updateIntersections <- function(S) {
.Call(C_UpdateIntersections, S, environment())
}
#' Construct section profiles
#'
#' Set up section profiles of spheroids for unfolding
#'
#' The function aggregates the necessary information for trivariate unfolding of spheroids' joint size-shape orientation distribution
#' of type either "\code{prolate}" or "\code{oblate}". The argument \code{size} is a numeric matrix of semi-axis lengths where the first column
#' corresponds to the major semi-axis and the second one to minor semi-axis. The orientation of an ellipse is assumed to be measured as the angle
#' between its major axis and vertical axis of the coordinate system in the intersection plane ('z' axis in 3D). For values in \eqn{[0,2\pi]} these
#' angles are automatically transformed to \eqn{[0,\pi/2]} as required by the unfolding procedure.
#'
#' @param size matrix of lengths of the semi-axes
#' @param alpha angle of section profiles in the plane (see details)
#' @param type name of the spheroid type, either "\code{prolate}" or "\code{oblate}" from which the
#' section profiles are assumed to come from
#'
#' @return The function returns a list which consists of either the longer or shorter
#' semi-axis length named \code{A} of section profiles corresponding to the type of spheroids used before and whose joint
#' joint distribution is to be estimated (by unfolding), the shape factor \code{S} of both semi-axes as the shape factor
#' between \eqn{(0,1]} and the orientation angle \code{alpha}, either of class "\code{prolate}" or "\code{oblate}".
#'
#' @examples
#' # load data set
#' data(data15p)
#'
#' # matrix of semi-axes lengths (major,minor)
#' AC <- data.matrix(data15p[c("A","C")])/1000
#'
#' # selecting the minor semi-axis for prolate type of spheroids:
#' # independent of nomenclature (always named \code{A})
#' sp <- sectionProfiles(AC,unlist(data15p["alpha"]))
#'
#' summary(sp$A) # here minor semi-axis because of prolate
#' summary(sp$S) # shape factor
#' summary(sp$alpha) # angle assumed to be w.r.t. (vertical) 'z' axis
#'
#' @author M. Baaske
#' @rdname sectionProfiles
#' @export
sectionProfiles <- function(size,alpha,type=c("prolate","oblate")) {
type <- match.arg(type)
stopifnot(is.matrix(size))
if(anyNA(size) || any(size<0))
stop("'size' must have non-negative values.")
if(anyNA(alpha) || !is.numeric(alpha) || any(alpha<0))
stop(paste("'alpha' must have non-negative values."))
if(max(alpha) > pi/2)
alpha <- try(sapply(alpha,.getAngle),silent=TRUE)
if(inherits(alpha,"try-error"))
warning("Could not compute 'alpha'. See its returned list element.")
structure(list("A"=if(type=="prolate") size[,2] else size[,1],
"S"=size[,2]/size[,1],
"alpha"=alpha),
class=type)
}
#' Poisson germ-grain process
#'
#' Simulation of Poisson germ-grain processes with either spheres, spheroids or spherocylinders as grains
#'
#' The function can simulate a Poisson germ-grain process according to the parameter \code{theta} within a predefined (3D) box.
#' The positions of the germs follow a uniform distribution according to a Poisson process with mean intensity parameter \code{lam}.
#' The function either randomly generates \code{type="prolate"} or \code{type="oblate"} spheroids, spheres or spherocylinders.
#' The argument \code{size} sets the name of the distribution function for the size/length of the objects, i.e. the major semi-axis
#' lengths in case of spheroids, radii for spheres or the lengths of the main axis of rotation for spherocylinders including the end caps.
#'
#' The following direction (orientation) distributions of the spheroids' major-axis, respectively, cylinders' main axis are available:
#' a uniform distribution ("\code{runifdir}"), distribution ("\code{rbetaiso}") and the "\emph{von Mises-Fisher}" ("\code{rvMisesFisher}")
#' distribution. The two last ones depend on the concentration parameter \code{kappa} which is set as part of the parameter list \code{theta}, see examples below.
#' The direction distributions generate random spherical coordinates \eqn{(\vartheta, \phi)} w.r.t. a fixed main orientation axis \code{mu}
#' with polar angle \eqn{\vartheta\in[0,\pi/2)} and azimuthal angle \eqn{\phi\in[0,2\pi)}. The simulations are always performed within a bounding 3D box which consists of a list specifying the ranges of each dimension corresponding to the lower and upper limits of the box in each direction. If the
#' argument \code{box} contains only a single range, i.e. \code{box=list(c(0,1))}, this limit isassumed for the remaining dimensions which is then simply replicated.
#' The optional argument \code{rjoint} defines a (joint) distribution function which can be any function provided by the user in order to generate
#' the required distributional parameters for the spheroids or spherocylinders. For an in-depth example of usage please see the workflow in 'simSpheroids.R'
#' and 'simCylinders.R'.
#'
#' In addition, the function supports an exact simulation type [2] of the grains. In case of spheroids and spherocylinders setting \code{size="rbinorm"}
#' declares a bivariate size-shape distribution for which the exact simulation is available. More specifically, for a bivariate normal random vector \eqn{[X,Y]}
#' with correlation parameter \eqn{\rho}, the length of the major semi-axis of a spheroid is given by \eqn{a=exp(x)} with a (logit-transformed) shape parameter
#' as \eqn{s=1/(1+exp(-y))} and thus a scaled minor semi-axis length \eqn{c=a*s}. The modification leads to a log-normally distributed length of the
#' major semi-axis. Consequently, in case of spherocylinders, the log-normally distributed length is \eqn{len=h+2*r} where \code{h} is the height and
#' \eqn{r=len/2*s} the radius. The main direction \code{u} of a spheroid or spherocylinder is determined by the major axis independent of size and shape.
#' Further, the following univariate distributions of the major semi-axis \code{a}, respectively, length \code{len} and shape \code{s} are available:
#' '\code{rbeta}', '\code{rgamma}', '\code{rlnorm}' and '\code{runif}'. One can also use '\code{const}' for simulations with constant lengths or shapes.
#' Note that only simulations with size distributions '\code{rbinorm}' or '\code{rlnorm}' can use the exact type of simulation.
#'
#' For spheres any distribution of the radii can be specified as a name of a user-defined function in the argument \code{size} as long as the formal named
#' function parameters match the actual names of the parameters exactly as defined in \code{theta}. Besides this, all other distributions given above are
#' also available. Using '\code{const}' simulates spheres of constant radii.
#'
#' The argument \code{pl>=0} denotes both the print level of intermediate output and by the same time the type of the return value. If \code{pl=10},
#' then an abbreviated list of spheroids or spheres is returned to speed up computation. Note that, the current implementation does not include routines
#' for unfolding the joint size-shape-orientation distribution of spherocylinders so far.
#'
#'
#' @param theta list of simulation parameters which must consist of elements: \code{size}, \code{shape} and \code{orientation}
#' @param lam mean number of objects per unit volume
#' @param size name of the size distribution function
#' @param shape name of the shape distribution function
#' @param orientation name of direction distribution function
#' @param type type of grain, either "\code{prolate}" or "\code{oblate}", "\code{spheres}", "\code{cylinders}"
#' @param rjoint user-defined function, which specifies the (joint) distribution of the size, shape and orientation
#' @param box simulation box
#' @param mu main orientation axis, \code{mu=c(0,0,1)} (default)
#' @param dz distance of the intersecting plane to the origin
#' @param n normal vector defining the intersecting plane
#' @param intersect options for type of return values: "\code{full}" for the simulated system together with section profiles as two lists named \code{S} and \code{sp} respectively,
#' choose "\code{only}" for section profiles only, "\code{original}" for the 3D system only and "\code{digi}" for a (binary) integer matrix \code{W} as a discretized
#' version of section profiles whose resolution depends on the chosen lattice constant \code{delta}
#' @param delta lattice constant for discretization, set to \code{0.01} (default)
#' @param intern logical, \code{FALSE} (default), whether to return only section profiles with centers inside the simulation window
#' @param perfect logical, \code{FALSE} (default), whether to simulate exactly (also called perfect)
#' @param pl integer, print level and return value type, see details
#' @param label character, a label set to each generated object, set to '\code{N}' (default)
#'
#' @return list of 3D objects depending on the chosen return type defined by the argument \code{intersect}
#'
#' @examples
#' # intensity parameter
#' lam <- 100
#'
#' # simulation bounding box
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # log normal size distribution with a constant shape factor and
#' # concentration parameter (\code{kappa=1}) for the orientation, see reference [1]
#' theta <- list("size"=list("meanlog"=-2.5,"sdlog"=0.5),
#' "shape"=list("s"=0.5),
#' "orientation"=list("kappa"=1))
#'
#' S <- simPoissonSystem(theta,lam,size="rlnorm",box=box,type="oblate",pl=1)
#' length(S)
#'
#' @references
#' \itemize{
#' \item{} {Ohser, J. and Schladitz, K. 3D images of materials structures Wiley-VCH, 2009}
#' \item{} {C. Lantu\eqn{\acute{\textrm{e}}}joul. Geostatistical simulation. Models and algorithms.
#' Springer, Berlin, 2002. Zbl 0990.86007}
#' }
#' @author M. Baaske
#' @rdname simPoissonSystem
#' @export
simPoissonSystem <- function(theta, lam, size="const", shape="const", orientation="rbetaiso",
type=c("prolate","oblate","spheres","cylinders"), rjoint=NULL,
box=list(c(0,1)), mu=c(0,0,1), dz=0, n=c(0,1,0),
intersect=c("original","full","only","digi"),
delta=0.01, intern=FALSE, perfect=FALSE, pl=0, label="N")
{
it <- pmatch(type,c("prolate","oblate","spheres","cylinders"))
if(length(it)==0 || is.na(it))
stop(paste("`type` is one of: ", paste0("`",c("prolate","oblate","sphere","cylinder"),"`",collapse=","),sep=""))
if(!is.numeric(lam) || !(lam>0) )
stop("Expected 'lam' as non-negative numeric argument")
if(length(box)==0 || !is.list(box))
stop("Expected argument 'box' as list type.")
if(length(box)==1)
box <- rep(box[1],3)
else if(length(box)!=3)
stop("Simulation box has wrong dimensions.")
names(box) <- c("xrange","yrange","zrange")
if(dz < min( box[[which(n == 1)]]) || dz > max(box[[which(n == 1)]]))
stop("The plane has to intersect the box.")
# digitization (resolution) factor
stopifnot(is.numeric(delta))
# spheroid type
type <- match.arg(type)
intersect <- match.arg(intersect)
cond <-
if(!is.null(rjoint))
{
if(!is.function(rjoint))
stop("Unknown multivariate distribution function.")
#largs <- theta[-(which(it==1))]
if(!is.list(theta) || length(theta) == 0L)
stop("Expected 'theta' as list of named arguments matching exactly with formal arguments of function `rjoint`.")
it <- match(names(theta),names(formals(rjoint)))
funname <- as.character(substitute(rjoint))
if(length(it) == 0L || anyNA(it))
stop(paste("Arguments must match formal arguments of function `",funname,"`.",sep=""))
# check function
funret <- try(do.call(rjoint,theta))
if(inherits(funret,"try-error"))
stop(paste("Error in user defined function `",funname,"`.",sep=""))
fargs <-
if(type %in% c("prolate","oblate")) {
if(!is.list(funret))
stop("Expected a list as return type in user defined function.")
c("a","b","c","u","shape","theta","phi")
} else if(type %in% c("cylinders")) {
if(!is.list(funret))
stop("Expected a list as return type in user defined function.")
c("h","r","u","theta","phi")
} else if(type %in% c("spheres")) {
if(!is.numeric(funret))
stop("Expected a numeric vector as return type in user defined function.")
c("r")
} else { stop("Undefined type!")}
if(any(!(fargs %in% names(funret)))){
stop("Argument names of return value list does not match required arguments.")
}
list("type"=type,"rdist"=funname,"box"=box,
"lam"=lam,"pl"=pl,"mu"=mu,"rho"= environment(),"label"=label,
"dz"=dz, "nsect"=n, "delta"=delta,"intern"=as.integer(intern),
"perfect"=as.integer(perfect), "intersect"=intersect)
} else {
# set defaults if not provided
if(!is.list(theta))
stop("Expected 'theta' as list of named arguments `size`, `shape`, `orientation`.")
nms <- c("size","shape","orientation")
it <- match(names(theta), nms)
if(length(it)>0L && anyNA(it))
names(theta)[it] <- nms[it]
if(shape == "const")
{
# set defaults
if(!is.list(theta$shape) || length(theta$shape) == 0L)
theta$shape <- list(1,1)
} else if(exists(shape, mode="function")) {
# check arguments of supported distribution functions
fargs <- names(formals(shape))
if(shape %in% c("rbeta","rgamma","runif"))
fargs <- fargs[-1]
if(!is.list(theta$shape) || length(theta$shape) == 0L)
stop("`shape` parameters must be given in `theta`.")
it <- match(names(theta$shape),fargs)
if(length(it)==0 || anyNA(it))
stop(paste0("Arguments of 'shape' must match formal arguments of function ",shape,sep=""))
} else {
stop(paste0("Undefined `", shape, "` distribution function."))
}
if(!is.list(theta$size) || length(theta$size)==0L)
stop("Arguments for size distribution must be given in `theta`.")
if(size == "rbinorm") {
fargs <- c("mx","my","sdx","sdy","rho","kappa")
it <- match(names(theta$size),fargs)
if(length(it)==0 || anyNA(it))
stop(paste("Arguments of 'size' must match formal arguments: ", paste0("`",fargs,"`",collapse=",")))
} else if(exists(size, mode="function")) {
# check arguments of supported distribution functions
fargs <- names(formals(size))
if(size %in% c("rlnorm","rbeta","rgamma","runif"))
fargs <- fargs[-1]
it <- match(names(theta$size),fargs)
if(length(it)==0 || anyNA(it))
stop(paste0("Arguments of 'size' must match formal arguments of function ",size,sep=""))
} else if(size == "const"){
# set defaults
if(!is.list(theta$size) || length(theta$size) == 0L)
theta$size <- list(1,1)
} else {
stop(paste0("Undefined `", size, "` distribution function."))
}
it <- pmatch(orientation,c("runifdir","rbetaiso","rvMisesFisher","const"))
if(is.na(it) && !exists(orientation, mode="function"))
stop("Undefined distribution function for the orientation/direction.")
# set defaults
if(!is.list(theta$orientation) || length(theta$orientation) == 0L)
theta$orientation <- list("kappa"=1)
list("type"=type, "lam"=lam,
"rdist"=list("size"=size, "shape"=shape, "orientation"=orientation),
"box"=box,"pl"=pl,"mu"=mu,"rho"= environment(), "dz"=dz, "nsect"=n, "delta"=delta,
"intern"=as.integer(intern),"label"=label, "perfect"=as.integer(perfect),
"intersect"=intersect)
}
.Call(C_PoissonSystem, theta, cond)
}
#' Coefficients for trivariate unfolding
#'
#' Compute coefficients of the discretized integral equation for unfolding
#'
#' In order to apply the Expectation Maximization (EM) algorithm to the stereological unfolding procedure for the joint
#' size-shape-orientation distribution in 3D one first has to compute the coefficients of a discretized integral equation.
#' This step is the most time consuming part of unfolding and therefore has been separated in its own function and can be called
#' separately if needed. The number of bin classes for the size, shape and orientation do not need to be the same, but the
#' given class limits are also used for binning the spatial values. One can define the number of cpu cores by the global option
#' \code{par.unfoldr} or passing the number of cores \code{nCores} directly to the function.
#'
#' @param breaks list of bin vectors, see \code{\link{setbreaks}}
#' @param stype type of spheroid, either "\code{prolate}" or "\code{oblate}"
#' @param check logical, whether to run some input checks
#' @param nCores number of cpu cores used to compute the coefficients
#'
#' @return numeric 6D array of coefficients
#'
#' @example inst/examples/coeffarray.R
#'
#' @author M. Baaske
#' @rdname coefficientMatrixSpheroids
#' @export
coefficientMatrixSpheroids <- function(breaks, stype=c("prolate","oblate"),
check=TRUE,nCores=getOption("par.unfoldr",1))
{
stype <- match.arg(stype)
it <- match(names(breaks), c("size","angle","shape"))
if (length(it)==0 || anyNA(it))
stop("Expected 'breaks' as named list of: 'size','angle','shape' ")
if (is.unsorted(breaks$size) || is.unsorted(breaks$angle) || is.unsorted(breaks$shape))
stop("'breaks' list must contain non-decreasingly sorted classes")
if(check) {
if(any(breaks$size<0))
stop("Breaks vector 'size' must have non-negative values.")
if(min(breaks$angle)<0 || max(breaks$angle)>pi/2)
stop(paste0("Breaks vector 'angle' must have values between zero and ",pi/2))
if(min(breaks$shape)<0 || max(breaks$shape)>1)
stop("Breaks vector 'shape' must have values between 0 and 1.")
}
.Call(C_CoefficientMatrixSpheroids,
breaks$size,breaks$angle,breaks$shape,
breaks$size,breaks$angle,breaks$shape, list(stype,nCores))
}
#' Vertical sections
#'
#' Compute vertical section profiles of a spheroid system
#'
#' The function intersects a spheroid system by a plane defined by the normal vector \code{n} either
#' equal to \code{c(0,1,0)} (default) or \code{c(1,0,0)}, which is called a vertical section. Depending on
#' the type of spheroid (either "\code{prolate} or "\code{oblate}") the returned semi-axis lengths are those
#' corresponding to the minor semi-axis or, respectively, major semi-axis in the way these are required for unfolding.
#'
#' @param S list of spheroids, see \code{\link{simPoissonSystem}}
#' @param d distance of the intersecting plane from the origin of the box
#' @param n normal vector which defines the interecting vertical plane
#' @param intern logical, \code{FALSE} (default), return all section profiles otherwise
#' only those which have their centers inside the correspondig intersection window
#'
#' @return list of sizes \code{A}, shape factors \code{S} and (vertical) angles \code{alpha}
#' of section profiles in the plane w.r.t the 'z' axis between \eqn{[0,\pi/2]}.
#'
#'
#' @examples
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # (exact) bivariate size-shape (isotropic) orientation distribution (spheroids)
#' theta <- list("size"=list("mx"=-2.5,"my"=0.5, "sdx"=0.35,"sdy"=0.25,"rho"=0.15),
#' "orientation"=list("kappa"=1))
#'
#' S <- simPoissonSystem(theta,lam=100,size="rbinorm",box=box,
#' type="prolate",perfect=TRUE,pl=1)
#'
#' sp <- verticalSection(S,d=2.5,n=c(0,1,0),intern=TRUE)
#' summary(sp$alpha)
#'
#' @author M. Baaske
#' @rdname verticalSection
#' @export
verticalSection <- function(S,d,n=c(0,1,0),intern=FALSE) {
if(!(class(S) %in% c("prolate","oblate")))
stop("Only spheroids of class 'prolate' or 'oblate' can be used.")
stopifnot(is.numeric(d))
stopifnot(is.logical(intern))
if(sum(n) != 1 || which(n == 1) > 2L)
stop("Normal vector is one of c(0,1,0) or c(1,0,0).")
ss <- .Call(C_IntersectPoissonSystem, S, list("nsect"=n,"dz"=d,"intern"=intern,"pl"=10), environment())
A <- if(class(S)=="prolate")
sapply(ss,"[[",2)
else sapply(ss,"[[",1)
## convert angle 'alpha' in the intersecting plane
## which is always between [0,pi/2] and w.r.t 'z' axis
alpha <- sapply(ss,"[[",4)
alpha <- try(0.5*pi-sapply(alpha,.getAngle),silent=TRUE) # alpha in [0,pi/2]
if(inherits(alpha,"try-error") || !is.numeric(alpha) || anyNA(alpha) )
warning("Could not compute angle 'alpha'. See the returned list element.")
structure(list("A"=A,
"S"=sapply(ss,"[[",3),
"alpha"=alpha), # w.r.t z axis in 3D
"win"=attr(ss,"win"),"plane"=attr(ss,"plane"),
class=class(S) )
}
#' Intersection in 3D
#'
#' Intersect a system of spheres, spheroids or spherocylinders by a plane
#'
#' The function intersects a given (Poisson) system made of spheres, spheroids or cylinders as grains by a plane defined by a
#' normal vector, e.g. \code{n=c(0,1,0)}, perpendicular to one of the bounding planes of the simulation box.
#' For a print level \code{pl>=0} some verbose output is given. Also it sets the type of return value. In case of spheroid
#' intersections, setting \code{pl=10}, leads to a short version of the full specification return list of section profiles
#' with elements named \code{A} (major semi-axis), \code{C} (minor semi-axis), \code{S} (the shape factor as the ratio of these two)
#' and \code{phi} as the angle in the intersecting plane between \eqn{[0,2\pi]} w.r.t. the 'x' axis. Otherwise additional components
#' are returned such as the ellipse`s rotation matrix, also named \code{A}, the center point \code{center} and a constant \code{type=10}
#' (defining the object of full ellipses among other types of intersection objects) of the section profiles. For sphere intersections only
#' a numeric vector of disc radii are returned as a short version of return values.
#'
#' @param S list of spheres, spheroids or spherocylinders, see \code{\link{simPoissonSystem}}
#' @param d distance of the the box-aligned intersecting plane from the origin
#' @param n normal vector which defines the intersecting plane
#' @param intern logical, \code{FALSE} (default), return all section profiles otherwise
#' only those which have their centers inside the corresponding intersection window
#' (if the intersected system had been simulated using exact simulation this makes sense)
#' @param pl integer, \code{pl=0} (default), for no verbose output and otherwise a full specification list of
#' section profiles in case of sphere and spheroid intersections
#'
#' @return For spheroid intersections the function returns a list of size, shape and angle of section profiles or a short version
#' of it; for sphere intersections either radii or lists containing the centers of discs and the object number.
#'
#'
#' @examples
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # constant size-shape orientation distribution (spheroids)
#' theta <- list("size"=list(0.1),"shape"=list(0.5), "orientation"=list("kappa"=10))
#'
#' S <- simPoissonSystem(theta,lam=100,box=box,type="prolate",pl=1)
#'
#' # return short version of section profiles
#' sp <- intersectSystem(S, 2.5, pl=10)
#'
#' @author M. Baaske
#' @rdname intersectSystem
#' @export
intersectSystem <- function(S, d, n=c(0,1,0), intern=FALSE, pl=0) {
if(!(class(S) %in% c("prolate","oblate","spheres","cylinders")))
stop("Unknown object class.")
stopifnot(is.numeric(d))
stopifnot(is.numeric(pl))
stopifnot(is.logical(intern))
if(sum(n)>1)
stop("Normal vector is like c(0,1,0). ")
box <- attr(S,"box")
if(is.null(box))
stop("'box' must be given as an attribute.")
if(d < min( box[[which(n == 1)]]) || d > max(box[[which(n == 1)]]))
stop("The plane has to intersect the box.")
.Call(C_IntersectPoissonSystem,S, list("nsect"=n,"dz"=d,"intern"=as.integer(intern),"pl"=as.integer(pl)), environment())
}
#' Spheroid system 3D
#'
#' Draw a spheroid system in 3D
#'
#' The function requires the package '\code{rgl}' to draw spheroids. For an example
#' please see the file 'simSpheroids.R'.
#'
#' @param S list of spheroids, see \code{\link{simPoissonSystem}}
#' @param box simulation box
#' @param draw.axes logical, whether to show the axes
#' @param draw.box logical, whether to show the bounding box
#' @param draw.bg logical, whether to show a gray background box correpsonding to the simulation box
#' @param bg.col background color used showw the box background
#' @param clipping logical, whether to clip all lateral planes of the simulation box
#' @param ... further material properties passed to 3D plotting functions
#'
#' @return NULL
#' @author M. Baaske
#' @rdname spheroids3d
#' @export
spheroids3d <- function(S, box, draw.axes=FALSE, draw.box=TRUE, draw.bg=TRUE,
bg.col="white", clipping=FALSE, ...)
{
if (!requireNamespace("rgl", quietly=TRUE))
stop("Please install 'rgl' package from CRAN repositories before running this function.")
ellipsoid3d <- function(rx=1,ry=1,rz=1,n=50,ctr=c(0,0,0), qmesh=FALSE,trans = rgl::par3d("userMatrix"),...) {
if (missing(trans) && !rgl::rgl.cur())
trans <- diag(4)
degvec <- seq(0,2*pi,length=n)
ecoord2 <- function(p) {
c(rx*cos(p[1])*sin(p[2]),ry*sin(p[1])*sin(p[2]),rz*cos(p[2]))
}
v <- apply(expand.grid(degvec,degvec),1,ecoord2)
if (qmesh)
v <- rbind(v,rep(1,ncol(v))) ## homogeneous
e <- expand.grid(1:(n-1),1:n)
i1 <- apply(e,1,function(z)z[1]+n*(z[2]-1))
i2 <- i1+1
i3 <- (i1+n-1) %% n^2 + 1
i4 <- (i2+n-1) %% n^2 + 1
i <- rbind(i1,i2,i4,i3)
if (!qmesh)
rgl::quads3d(v[1,i],v[2,i],v[3,i],...)
else
return(rgl::rotate3d(rgl::qmesh3d(v,i,material=...),matrix=trans))
}
sphere <- ellipsoid3d(qmesh=TRUE,trans=diag(4))
spheroid3d <- function (x=0,y=0,z=0, a=1, b=1,c=1, rotM, subdivide = 3, smooth = TRUE){
result <- rgl::scale3d(sphere, a,b,c)
result <- rgl::rotate3d(result,matrix=rotM)
result <- rgl::translate3d(result, x,y,z)
invisible(result)
}
N <- length(S)
ll <- lapply(S, function(x)
spheroid3d(x$center[1],x$center[2],x$center[3],
x$acb[1],x$acb[2],x$acb[3],rotM=x$rotM))
rgl::shapelist3d(ll,...)
x <- box$xrange[2]
y <- box$yrange[2]
z <- box$zrange[2]
## draw gray box
if(draw.bg) {
c3d.origin <- rgl::translate3d(rgl::scale3d(rgl::cube3d(col=bg.col, alpha=0.1),x/2,y/2,z/2),x/2,y/2,z/2)
rgl::shade3d(c3d.origin)
}
if(clipping) {
rgl::clipplanes3d(-1,0,0,x)
rgl::clipplanes3d(0,-1,0,y)
rgl::clipplanes3d(0,0,-1,z)
rgl::clipplanes3d(1,0,0,0)
rgl::clipplanes3d(0,1,0,0)
rgl::clipplanes3d(0,0,1,0)
}
if(draw.axes) {
rgl::axes3d(c('x','y','z'), pos=c(0,0,0))
rgl::title3d('','','x','y','z')
}
## draw box
if(draw.box) {
rgl::axes3d(edges = "bbox",labels=TRUE,tick=FALSE,box=TRUE,nticks=0,
expand=1.0,xlen=0,xunit=0,ylen=0,yunit=0,zlen=0,zunit=0)
}
}
#' Cylinder system 3D
#'
#' Draw spherocylinders in 3D
#'
#' The function requires the package '\code{rgl}' to draw spherocylinders into a 3D '\code{rgl}' image.
#' For an example please see the file 'simCylinders.R'.
#'
#' @param S list of cylinders, see \code{\link{simPoissonSystem}}
#' @param box simulation box
#' @param draw.axes logical, whether to show the axes
#' @param draw.box logical, whether to show the bounding box
#' @param clipping logical, whether to clip all lateral planes of the simulation box
#' @param ... further material properties passed to 3d plotting functions
#'
#' @return NULL
#' @author M. Baaske
#' @rdname cylinders3d
#' @export
cylinders3d <- function(S, box, draw.axes=FALSE, draw.box=TRUE, clipping=FALSE,...)
{
if (!requireNamespace("rgl", quietly=TRUE))
stop("Please install 'rgl' package from CRAN repositories before running this function.")
cylinder <- function(m, radius=1, h=1, rotM=diag(3), u=c(0,0,1)) {
cyl <- rgl::cylinder3d(rbind(c(0,0,0),c(0,0,1)), radius=radius,
e1=cbind(0, 0, 1), e2=cbind(1, 0, 0), sides=25,closed=-2 )
result <-rgl::scale3d(cyl,1,1,h)
result <- rgl::rotate3d(result,matrix=rotM)
## translate to midpoint of cylinder, center of mass
m <- m-0.5*h*u
result <- rgl::translate3d(result, m[1],m[2],m[3])
invisible(result)
}
args <- list(...)
ok <- sapply(S, function(x) (!is.null(x$h) && x$h > 0.0) ) # only true cylinders not 'spheres'
cyls <- lapply(S[ok], function(x) { cylinder(x$center, x$r, x$h, x$rotM, x$u) })
rgl::shapelist3d(cyls,...)
if("col" %in% names(args)) {
cols <- rep(rep(args$col,length.out=length(S)),each=2)
args$col <- NULL
} else cols <- "black"
# spheres
is.sphere <- sapply(S,function(x) !(is.null(x$h) || x$h == 0.0) )
Xc <- do.call(rbind,lapply(S[ok],function(x) rbind( c(x$origin0,x$r),c(x$origin1,x$r))))
rgl::spheres3d(Xc,radius=Xc[,4],col=cols,unlist(args))
if(!all(ok)) {
Xc <- do.call(rbind,lapply(S[!ok],function(x) rbind(c(x$center,x$r))))
rgl::spheres3d(Xc,radius=Xc[,4],col="darkgray",unlist(args))
}
x <- box$xrange[2]; y <- box$yrange[2]; z <- box$zrange[2]
c3d.origin <- rgl::translate3d(rgl::scale3d(rgl::cube3d(col="darkgray", alpha=0.1),x/2,y/2,z/2),x/2,y/2,z/2)
rgl::shade3d(c3d.origin)
if(clipping) {
rgl::clipplanes3d(-1,0,0,box$xrange[1])
rgl::clipplanes3d(0,-1,0,box$yrange[1])
rgl::clipplanes3d(0,0,-1,box$zrange[1])
rgl::clipplanes3d(1,0,0,box$xrange[2])
rgl::clipplanes3d(0,1,0,box$yrange[2])
rgl::clipplanes3d(0,0,1,box$zrange[2])
}
if(draw.axes) {
rgl::axes3d(c('x','y','z'), pos=c(0,0,0))
rgl::title3d('','','x','y','z')
}
if(draw.box) {
rgl::axes3d(edges = "bbox",labels=TRUE,tick=FALSE,box=TRUE,nticks=0,
expand=1.0,xlen=0,xunit=0,ylen=0,yunit=0,zlen=0,zunit=0)
}
}
#' Spheroid intersections 3D
#'
#' Draw section profiles of spheroids in 3D
#'
#' The function requires the package '\code{rgl}' for drawing section profiles (ellipses)
#' into a 3D '\code{rgl}' image. For a full example please see the file 'simSpheroids.R'.
#'
#' @param E a list of spheroid intersections, see \code{\link{intersectSystem}}
#' @param n the normal vector of the intersecting plane
#' @param np number of points for a polygon approximation of the ellipses
#'
#' @return NULL
#'
#' @author M. Baaske
#' @rdname drawSpheroidIntersection
#' @export
drawSpheroidIntersection <- function(E, n=c(0,1,0), np=25) {
ind <- which(n == 0)
.pointsOnEllipse <- function(E,t) {
E$phi <- E$phi + 0.5*pi
xt <- E$center[1] + E$ab[1]*cos(t)*cos(E$phi)-E$ab[2]*sin(t)*sin(E$phi)
zt <- E$center[2] + E$ab[1]*cos(t)*sin(E$phi)+E$ab[2]*sin(t)*cos(E$phi)
yt <- rep(0,length(t))
m <- matrix(0,nrow=length(xt),ncol=3)
m[,ind[1]] <- xt
m[,ind[2]] <- zt
m
}
.plotEllipse <- function(x) {
M <- .pointsOnEllipse(x,t)
rgl::polygon3d(M[,1],M[,2],M[,3],fill=TRUE,coords=ind)
}
s <- 2*pi/np
t <- seq(from=0,to=2*pi,by=s)
invisible(lapply(E,function(x) .plotEllipse(x) ))
}
.checkArgs <- function(optlist, options) {
if(is.null(names(options)))
stop("Options should be a list of named arguments.")
if (!is.list(options) || "" %in% names(options))
stop("Argument ",as.character(substitute(options)), " must be a list of named (character) elents.")
optnames <- (names(options) %in% names(optlist))
if (!all(optnames)) {
unames <- as.list(names(options)[!(optnames)])
stop(paste(c("Unknown arguments in ",as.character(substitute(options))," : ",do.call("paste", c(unames, sep = ", "))), collapse=" "))
}
return (0)
}
#' Sphere planar sections
#'
#' Planar intersection of spheres
#'
#' The function computes the planar intersection of a sphere system, i.e. an intersection with the plane whose normal
#' vector is given by \code{c(0,0,1)} and returns the diameters of the resulting discs.
#'
#' @param S list of spheres of class \code{sphere}, see \code{\link{simPoissonSystem}}
#' @param d distance of the (planar) xy-plane to the origin
#' @param intern logical, \code{FALSE} (default), whether to return all discs or
#' only those which have their centers inside the intersecting window
#' @param pl print level, \code{pl>0} for some verbose output
#'
#' @return numeric vector of disc diameters
#'
#' @examples
#' lam <- 100
#'
#' # parameter rlnorm distribution (radii)
#' theta <- list("size"=list("meanlog"=-2.5,"sdlog"=0.5))
#'
#' # simulation bounding box
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#'
#' # simulate only 3D system
#' S <- simPoissonSystem(theta,lam,size="rlnorm",box=box,type="spheres",
#' intersect="original", pl=1)
#'
#' # return only objects whose centers are within
#' # the intersection window
#' sp <- planarSection(S,d=2.5,intern=TRUE,pl=1)
#'
#' # histogram of diameters
#' hist(sp)
#' summary(sp)
#'
#' # distribution of radii
#' mean(log(sp/2))
#' sd(log(sp/2))
#'
#'
#' @author M. Baaske
#' @rdname planarSection
#' @export
planarSection <- function(S, d, intern=FALSE, pl=0) {
if(!(c("spheres") %in% class(S) ))
stop("Expected spheres as a list.")
if(!is.numeric(d) || d<0)
stop("'d' must be numeric/positive.")
stopifnot(is.logical(intern))
sp <- .Call(C_IntersectPoissonSystem,S,list("nsect"=c(0,0,1),"dz"=d,"intern"=intern,"pl"=pl),environment())
if(is.list(sp)) # full list of section profiles
return (sapply(sp,function(x) 2.0*x$r))
else return (2*sp) # only diameters are returned
}
#' Binning numeric values
#'
#' Vector of count data
#'
#' The function provides basic binning (grouping) of numeric values into
#' classes defined by the breaks vector \code{bin}. The values are binned
#' according to \eqn{bin[i[j]]<x[j]\leq bin[i[j]+1]} for intervals \eqn{i=1,...,N-1}
#' and \code{length(bin)=N} of values \eqn{x[j]}, \eqn{j=1,...,|x|}. If \eqn{x[j] > bin[N]} or \eqn{x[j] < bin[1]} then \eqn{x[j]}
#' is not counted at all.
#'
#' @param x numeric values to be binned
#' @param bin non-decreasingly sorted breaks vector
#' @param na.rm logical, default \code{FALSE}, whether to remove missing values, including \code{NaN} in \code{x}
#'
#' @return Vector of count data
#'
#' @examples
#' x <- runif(100,0,1)
#' bin <- seq(0,1,by=0.1)
#' binning1d(x,bin)
#'
#' @author M. Baaske
#' @rdname binning1d
#' @export
binning1d <- function(x,bin, na.rm = FALSE) {
if (anyNA(x)) {
if(na.rm) x <- x[which (!is.na(x))]
else stop("Vectors contains missing values or NAs.")
}
if(anyNA(bin) )
stop("'bin' vector contains missing values or NAs.")
if (is.unsorted(bin))
stop("'bin' must be sorted non-decreasingly")
.Call(C_Binning1d,x,bin)
}
#' Coefficients for Expectation Maximization algorithm
#'
#' Matrix of coefficients for Wicksell's corpuscle problem
#'
#' The function computes the matrix of coefficients for solving the discretized integral equation of
#' the Wicksell's corpuscle problem.
#'
#' @param bin non-decreasing vector of class limits
#' @return array of coefficients
#'
#' @references
#' Ohser, J. and Muecklich, F. Statistical analysis of microstructures in materials science J. Wiley & Sons, 2000
#'
#' @author M. Baaske
#' @rdname coefficientMatrixSpheres
#' @export
coefficientMatrixSpheres <- function(bin) {
## count of bin classes
n <- length(bin)-1
if (n<=0)
stop("'bin' has length zero")
if (anyNA(bin))
stop("Vectors contain NA values")
if (is.unsorted(bin))
stop("'bin' must be sorted non-decreasingly")
if(any(bin<0))
stop("Breaks vector must have non-negative values.")
p <- .C(C_em_saltykov_p, as.integer(n),as.numeric(bin),
p=as.numeric(matrix(0,n,n)))$p
dim(p) <- c(n,n)
return (p)
}
#' Expectation Maximization (EM) algorithm
#'
#' Estimation of the empirical sphere diameter distribution
#'
#' The function performs the EM algorithm, see reference below.
#'
#' @param y vector of observed absolute frequencies of circle diameters
#' @param bin non-decreasing vector of class limits
#' @param maxIt maximum number of iterations to be used
#' @return vector of count data of absolute frequenties of sphere diameters
#'
#' @example inst/examples/simSpheres.R
#'
#' @references
#' Ohser, J. and Muecklich, F. Statistical analysis of microstructures in materials science J. Wiley & Sons, 2000
#'
#' @author M. Baaske
#' @rdname em.saltykov
#' @export
em.saltykov <- function(y,bin,maxIt=32) {
if (length(y)==0) stop("input array 'y' has length zero")
if (anyNA(y) || anyNA(bin) )
stop("Vectors contain NA values")
if (is.unsorted(bin))
stop("'bin' must be sorted non-decreasingly")
n <- length(y)
p <- .C(C_em_saltykov_p, as.integer(n),as.numeric(bin),
p=as.vector(matrix(0,n,n)))$p
theta <- y+1.0e-6
.C(C_em_saltykov,as.integer(n), as.integer(maxIt),
as.numeric(p),as.numeric(y),theta=as.numeric(theta))$theta
}
#' Digitization
#'
#' Digitize 2D section profiles
#'
#' A list of section profiles, e.g. discs, ellipses or segments from intersected spherocylinders, is digitized according
#' to a given resolution according to the lattice constant \code{delta} such that the result is an integer matrix
#' which can be interpretated as a binary image. An intersection window can be either provided by the user w.r.t. \eqn{[l,u]^2}
#' where \code{l,u} are lower, respectively, upper bounds of corresponding to a simulation box, or, it is taken from the section
#' profiles \code{sp} which stores the intersection window as an attribute.
#'
#' @param sp list of section profiles, see \code{\link{intersectSystem}}
#' @param delta lattice constant for discretization of section profiles
#' @param win list of length two, the intersection window, default \code{NULL}
#'
#' @examples
#' # simulation box
#' box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
#' # (exact) bivariate size-shape (isotropic) orientation distribution (spheroids)
#' theta <- list("size"=list("mx"=-2.5,"my"=0.5, "sdx"=0.35,"sdy"=0.25,"rho"=0.15),
#' "orientation"=list("kappa"=1))
#'
#' # return only 3D system
#' S <- simPoissonSystem(theta,lam=100,size="rbinorm",box=box,type="prolate",
#' intersect="original",n=c(0,1,0),mu=c(0,0,1),perfect=TRUE,pl=1)
#'
#' # vertical intersection w.r.t. 'mu' (z axis, see above)
#' sp <- intersectSystem(S, 2.5)
#'
#' # show intersecting window
#' win <- attr(sp,"win")
#'
#' # digitize (could also pass some 'win' as an argument)
#' W <- digitizeProfiles(sp, delta=0.01, win = NULL)
#' image(1:nrow(W),1:ncol(W),W,col=gray(1:0))
#'
#' @return binary (integer) matrix as an image
#'
#' @author M. Baaske
#' @rdname digitizeProfiles
#' @export
digitizeProfiles <- function(sp, delta = 0.01, win = NULL ) {
stopifnot(is.numeric(delta))
if(!is.null(win)){
if(!is.list(win) || length(win) == 0L )
stop("Expected argument 'win' as list type.")
if(length(win) == 1L)
win <- rep(win[1],2)
else if(length(win) != 2L)
stop("Digitization window 'win' must be two-dimensional.")
}
if(delta < 0 || delta > 1.0)
stop("'delta' should be positive and equal or less than 1.")
.Call(C_DigitizeProfiles, sp, delta, win, environment())
}
Try the unfoldr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.