Nothing
##' Class containing functions to reconstruct \link{StitchedOutline}s
##' and store the associated data
##'
##' @description The function \code{reconstruct} reconstructs outline
##' into spherical surface Reconstruct outline into spherical
##' surface.
##' @importFrom geometry tsearch sph2cart
##' @author David Sterratt
##' @export
ReconstructedOutline <- R6Class("ReconstructedOutline",
inherit = OutlineCommon,
private = list(
## @field ims spherical coordinates of pixel corners
ims = NULL
),
public = list(
##' @field ol Annotated outline
ol = NULL,
##' @field ol0 Original Annotated outline
ol0 = NULL,
##' @field Pt Transformed cartesian mesh points
Pt = NULL,
##' @field Tt Transformed triangulation
Tt = NULL,
##' @field Ct Transformed links
Ct = NULL,
##' @field Cut Transformed links
Cut = NULL,
##' @field Bt Transformed binary vector representation
##' of edge indices onto a binary vector representation of the
##' indices of the points linked by the edge
Bt = NULL,
##' @field Lt Transformed lengths
Lt = NULL,
##' @field ht Transformed correspondences
ht = NULL,
##' @field u Indices of unique points in untransformed space
u = NULL,
##' @field U Transformed indices of unique points in untransformed space
U = NULL,
##' @field Rsett Transformed rim set
Rsett = NULL,
##' @field i0t Transformed marker
i0t = NULL,
##' @field H mapping from edges onto corresponding edges
H = NULL,
##' @field Ht Transformed mapping from edges onto corresponding edges
Ht = NULL,
##' @field phi0 Rim angle
phi0 = NULL,
##' @field R Radius of spherical template
R = NULL,
##' @field lambda0 Longitude of pole on rim
lambda0 = NULL,
##' @field lambda Longitudes of transformed mesh points
lambda = NULL,
##' @field phi Latitudes of transformed mesh points
phi = NULL,
##' @field Ps Location of mesh point on sphere in spherical coordinates
Ps = NULL,
##' @field n Number of mesh points
n = 500,
##' @field alpha Weighting of areas in energy function
alpha = 8,
##' @field x0 Area cut-off coefficient
x0 = 0.5,
##' @field nflip0 Initial number flipped triangles
nflip0 = NULL,
##' @field nflip Final number flipped triangles
nflip = NULL,
##' @field opt Optimisation object
opt = NULL,
##' @field E.tot Energy function including area
E.tot = NULL,
##' @field E.l Energy function based on lengths alone
E.l = NULL,
##' @field mean.strain Mean strain
mean.strain = NULL,
##' @field mean.logstrain Mean log strain
mean.logstrain = NULL,
##' @field debug Debug function
debug = NULL,
##' @description Load \link{AnnotatedOutline} into ReconstructedOutline object
##' @param ol \code{\link{AnnotatedOutline}} object, containing the following information
##' @param n Number of points in triangulation.
##' @param alpha Area scaling coefficient
##' @param x0 Area cut-off coefficient
##' @param plot.3d Whether to show 3D picture during optimisation.
##' @param dev.flat Device to plot grid onto. Value of \code{NA} (default)
##' means no plotting.
##' @param dev.polar Device display projection. Value of NA
##' (default) means no plotting.
##' @param report Function to report progress.
##' @param debug If \code{TRUE} print extra debugging output
loadOutline = function(ol,
n=500, alpha=8, x0=0.5,
plot.3d=FALSE, dev.flat=NA, dev.polar=NA, report=retistruct::report,
debug=FALSE) {
self$ol0 <- ol$clone()
self$n <- n
self$alpha <- alpha
self$x0 <- x0
self$debug <- debug
ol$triangulate()
ol$stitchTears()
ol$triangulate(suppress.external.steiner=TRUE)
if (length(ol$corrs)) {
ol$triangulate(suppress.external.steiner=TRUE)
}
## Transform the rim set
## ol$orderRset()
self$ol <- ol
self$phi0 <- ol$phi0
self$lambda0 <- ol$lambda0
report("Merging points...")
self$mergePointsEdges()
report("Projecting to sphere...")
self$projectToSphere()
},
##' @description Reconstruct
##' Reconstruction proceeds in a number of stages:
##'
##' \enumerate{
##'
##' \item The flat object is triangulated with at least \code{n}
##' triangles. This can introduce new vertices in the rim.
##'
##' \item The triangulated object is stitched.
##'
##' \item The stitched object is triangulated again, but this time it
##' is not permitted to add extra vertices to the rim.
##'
##' \item The corresponding points determined by the stitching process
##' are merged to form a new set of merged points and a new
##' triangulation.
##'
##' \item The merged points are projected roughly to a sphere.
##'
##' \item The locations of the points on the sphere are moved so as to
##' minimise the energy function.
##' }
##'
##' @param plot.3d If \code{TRUE} make a 3D plot in an RGL window
##' @param dev.flat Device handle for plotting flatplot updates to. If
##' \code{NA} don't make any flat plots
##' @param dev.polar Device handle for plotting polar plot updates
##' to. If \code{NA} don't make any polar plots.
##' @param Control argument to pass to \code{optim}
##' @param report Function to report progress.
reconstruct = function(plot.3d=FALSE, dev.flat=NA, dev.polar=NA,
report=getOption("retistruct.report")) {
## ## Initial plot in 3D space
## if (plot.3d) {
## sphericalplot(r)
## }
## }
## Check for flipped triangles and record initial number
ft <- flipped.triangles(self$getPoints(), self$Tt, self$R)
self$nflip0 <- sum(ft$flipped)
report("Optimising mapping with no area constraint using BFGS...")
self$optimiseMapping(alpha=0, x0=0, nu=1,
plot.3d=plot.3d,
dev.flat=dev.flat, dev.polar=dev.polar)
report("Optimising mapping with area constraint using FIRE...")
## FIXME: Need to put in some better heuristics for scaling
## maxmove, and perhaps other parameters
self$optimiseMappingCart(alpha=self$alpha, x0=self$x0, nu=1,
dtmax=500, maxmove=0.002*sqrt(self$ol$A.tot),
tol=1e-5,
plot.3d=plot.3d,
dev.flat=dev.flat, dev.polar=dev.polar)
report("Optimising mapping with strong area constraint using BFGS...")
self$optimiseMapping(alpha=self$alpha, x0=self$x0, nu=1,
plot.3d=plot.3d,
dev.flat=dev.flat, dev.polar=dev.polar)
report("Optimising mapping with weak area constraint using BFGS...")
self$optimiseMapping(alpha=self$alpha, x0=self$x0, nu=0.5,
plot.3d=plot.3d,
dev.flat=dev.flat, dev.polar=dev.polar)
report(paste("Mapping optimised. Deformation energy E:", format(self$opt$value, 5),
";", self$nflip, "flipped triangles."))
},
##' @description Merge stitched points and edges.
##' Create merged and transformed versions (all suffixed with \code{t})
##' of a number of existing variables, as well
##' as a matrix \code{Bt}, which maps a binary vector representation
##' of edge indices onto a binary vector representation of the
##' indices of the points linked by the edge.
##' Sets following fields
##' \itemize{
##' \item{\code{Pt}}{Transformed point locations}
##' \item{\code{Tt}}{Transformed triangulation}
##' \item{\code{Ct}}{Transformed connection set}
##' \item{\code{Cut}}{Transformed symmetric connection set}
##' \item{\code{Bt}}{Transformed binary vector representation
##' of edge indices onto a binary vector representation of the
##' indices of the points linked by the edge}
##' \item{\code{Lt}}{Transformed edge lengths}
##' \item{\code{ht}}{Transformed correspondences}
##' \item{\code{u}}{Indices of unique points in untransformed space}
##' \item{\code{U}}{Transformed indices of unique points in untransformed space}
##' \item{\code{Rset}}{The set of points on the rim (which has been reordered)}
##' \item{\code{Rsett}}{Transformed set of points on rim}
##' \item{\code{i0t}}{Transformed index of the landmark}
##' \item{\code{H}}{mapping from edges onto corresponding edges}
##' \item{\code{Ht}}{Transformed mapping from edges onto corresponding edges}
##' }
mergePointsEdges = function() {
h <- self$ol$h
T <- self$ol$T
Cu <- self$ol$Cu
L <- self$ol$L
P <- self$ol$getPointsScaled()
gf <- self$ol$gf
## Form the mapping from a new set of consecutive indices
## the existing indices onto the existing indices
u <- unique(h)
## Transform the point set into the new indices
Pt <- P[u,]
## Transform the point correspondance mapping to the new index space
ht <- c()
for (i in 1:length(h)) {
ht[i] <- which(u == h[i])
}
## DOESN'T WORK
## Form the inverse mapping from the existing indices to the new
## set of consecutive indices
## uinv <- c()
## uinv[u] <- 1:length(u)
## ht <- uinv[h[u]]
## Transform the triangulation to the new index space
Tt <- matrix(ht[T], ncol=3)
## Tansform the forward pointer into the new indices
gft <- ht[gf]
## Determine H, the mapping from edges onto corresponding edges
Cut <- matrix(ht[Cu], ncol=2)
Cut <- t(apply(Cut, 1, sort))
M <- nrow(Cut)
H <- rep(0, M)
for (i in 1:M) {
if (!H[i]) {
H[i] <- i
for (j in i:M) {
if (identical(Cut[i,], Cut[j,])) {
H[j] <- i
}
}
}
}
## Form the mapping from a new set of consecutive edge indices
## onto the existing edge indices
U <- unique(H)
## Transform the edge set into the new indices
Cut <- Cut[U,]
## Transform the edge correspondance mapping to the new index space
Ht <- c()
for (i in 1:length(H)) {
Ht[i] <- which(U == H[i])
}
## Create the lengths of the merged edges by averaging
Lt <- c()
for (k in 1:length(U)) {
is <- which(Ht == k)
## if (length(is)>1) {
## print(L[is])
## }
Lt[k] <- mean(L[is])
}
## Transform the rim set
Rset <- order.Rset(self$ol$getRimSet(), self$ol$gf, self$ol$h)
Rsett <- unique(ht[Rset])
i0t <- ht[self$ol$i0]
## Create the symmetric connection set
Ct <- rbind(Cut, Cut[,2:1])
## Matrix to map line segments onto the points they link
## Bt <- Matrix(0, nrow(Pt), nrow(Ct), sparse=TRUE)
Bt <- matrix(0, nrow(Pt), nrow(Ct))
for (i in 1:nrow(Ct)) {
Bt[Ct[i,1],i] <- 1
}
self$Pt = Pt
self$Tt = Tt
self$Ct = Ct
self$Cut = Cut
self$Bt = Bt
self$Lt = Lt
self$ht = ht
self$u = u
self$U = U
self$Rsett = Rsett
self$i0t = i0t
self$H = H
self$Ht= Ht
},
##' @description Project mesh points in the flat outline onto a sphere
##' This takes the mesh points from the flat outline and maps them to
##' the curtailed sphere. It uses the area of the flat outline and
##' \code{phi0} to determine the radius \code{R} of the sphere. It
##' tries to get a good first approximation by using the function
##' \code{\link{stretchMesh}}.
##' The following fields are set:
##' \itemize{
##' \item{\code{phi}}{Latitude of mesh points.}
##' \item{\code{lmabda}}{Longitude of mesh points.}
##' \item{\code{R}}{Radius of sphere.}
##' }
projectToSphere = function() {
Rsett <- self$Rsett
i0t <- self$i0t
A.tot <- self$ol$A.tot
Cut <- self$Cut
Lt <- self$Lt
phi0 <- self$phi0
lambda0 <- self$lambda0
Nt <- nrow(self$Pt)
Nphi <- Nt - length(Rsett)
## From this we can infer what the radius should be from the formula
## for the area of a sphere which is cut off at a latitude of phi0
## area = 2 * PI * R^2 * (sin(phi0)+1)
R <- sqrt(A.tot/(2*pi*(sin(phi0)+1)))
## Find lengths between successive points on rim
C <- matrix(NA, nrow(self$Pt), nrow(self$Pt))
for (i in 1:nrow(Cut)) {
C[Cut[i,1],Cut[i,2]] <- Lt[i]
C[Cut[i,2],Cut[i,1]] <- Lt[i]
}
L.Rsett <- rep(NA, length(Rsett))
for (i in 1:length(Rsett)) {
L.Rsett[i] <- C[Rsett[i],Rsett[mod1(i+1, length(Rsett))]]
}
## Check that this length matches the length computed from the AnnotatedOutline
## FIXME - this doesn't work for one retina - need to check why
## if (sum(L.Rsett) != sum(getRimLengths(r))) {
## stop("Internal error: Mismatch in rim lengths")
## }
## Stretch mesh points to circle
Ps <- stretchMesh(Cut, Lt, Rsett, circle(L=L.Rsett))
x <- Ps[,1]
y <- Ps[,2]
phi <- -pi/2 + sqrt(x^2 + y^2)*(phi0+pi/2)
phi[Rsett] <- phi0
lambda <- atan2(y, x)
lambda <- lambda - lambda[i0t] + lambda0
self$phi <- phi
self$lambda <- lambda
self$R <- R
self$phi0 <- phi0
self$lambda0 <- lambda0
self$Ps <- Ps
},
##' @description Return strains edges are under in spherical retina
##' Set information about how edges on the sphere
##' have been deformed from their flat state.
##' @return A list containing two data frames \code{flat} and \code{spherical}.
##' Each data frame contains for each edge in the flat or spherical meshes:
##' \itemize{
##' \item{\code{L}}{Length of the edge in the flat outline }
##' \item{\code{l}}{Length of the corresponding edge on the sphere}
##' \item{\code{strain}}{The strain of each connection}
##' \item{\code{logstrain}}{The logarithmic strain of each connection}
##' }
getStrains = function() {
## Original lengths in flattened outline is a vector with
## M elements, the number of rows of Cu
L <- self$ol$L
## New lengths in reconstructed object is a vector with Mt < M
## elements, the number of rows of Cut
lt <- compute.lengths(self$phi, self$lambda, self$Cut, self$R)
## For each connection in the flattened object, we want the length of
## the corresponding connection in the reconstructed object
## The mapping Ht achieves this
l <- lt[self$Ht]
stretch <- l/L
strain <- stretch - 1
logstrain <- log(stretch)
## Compute quantities in spherical retina too
Lt <- self$Lt
stretcht <- lt/Lt
straint <- stretcht - 1
logstraint <- log(stretcht)
return(list(flat=
data.frame(L=L, l=l,
strain=strain, logstrain=logstrain),
spherical=
data.frame(L=Lt, l=lt,
strain=straint, logstrain=logstraint)))
},
##' @description Optimise the mapping from the flat outline to the sphere
##' @param alpha Area penalty scaling coefficient
##' @param x0 Area penalty cut-off coefficient
##' @param nu Power to which to raise area
##' @param optim.method Method to pass to \code{optim}
##' @param plot.3d If \code{TRUE} make a 3D plot in an RGL window
##' @param dev.flat Device handle for plotting flatplot updates to. If
##' \code{NA} don't make any flat plots
##' @param dev.polar Device handle for plotting polar plot updates
##' to. If \code{NA} don't make any polar plots.
##' @param control Control argument to pass to \code{optim}
optimiseMapping = function(alpha=4, x0=0.5, nu=1, optim.method="BFGS",
plot.3d=FALSE, dev.flat=NA, dev.polar=NA,
control=list()) {
phi <- self$phi
lambda <- self$lambda
R <- self$R
phi0 <- self$phi0
lambda0 <- self$lambda0
Tt <- self$Tt
A <- self$ol$A
Cut <- self$Cut
Ct <- self$Ct
Lt <- self$Lt
Bt <- self$Bt
Rsett <- self$Rsett
i0t <- self$i0t
Nt <- nrow(self$Pt)
Nphi <- Nt - length(Rsett)
## Optimisation and plotting
opt <- list()
opt$p <- c(phi[-Rsett], lambda[-i0t])
opt$conv <- 1
count <- 0
while (opt$conv) {
## Optimise
opt <- stats::optim(opt$p, E, gr=dE,
method=optim.method,
T=Tt, A=A, Cu=Cut, C=Ct, L=Lt, B=Bt, R=R,
alpha=alpha, N=Nt, x0=x0, nu=nu,
Rset=Rsett, i0=i0t, phi0=phi0, lambda0=lambda0, Nphi=Nphi,
verbose=FALSE, control=control)
## Report
E.tot <- E(opt$p, Cu=Cut, C=Ct, L=Lt, B=Bt, R=R, T=Tt, A=A,
alpha=alpha, N=Nt, x0=x0, nu=nu,
Rset=Rsett, i0=i0t, phi0=phi0, lambda0=lambda0, Nphi=Nphi)
E.l <- E(opt$p, Cu=Cut, C=Ct, L=Lt, B=Bt, R=R, T=Tt, A=A,
alpha=0, N=Nt, x0=x0, nu=nu,
Rset=Rsett, i0=i0t, phi0=phi0, lambda0=lambda0, Nphi=Nphi)
ft <- flipped.triangles(cbind(phi=phi, lambda=lambda), Tt, R)
nflip <- sum(ft$flipped)
report(sprintf("E = %8.5f | E_L = %8.5f | E_A = %8.5f | %3d flippped triangles", E.tot, E.l, E.tot - E.l, nflip))
if (nflip & self$debug) {
print(data.frame(rbind(id=which(ft$flipped),
A=A[ft$flipped],
a=ft$areas[ft$flipped])))
}
## Decode p vector
phi <- rep(phi0, Nt)
phi[-Rsett] <- opt$p[1:Nphi]
lambda <- rep(lambda0, Nt)
lambda[-i0t] <- opt$p[Nphi+1:(Nt-1)]
self$phi <- phi
self$lambda <- lambda
self$opt <- opt
self$nflip <- sum(ft$flipped)
self$E.tot <- E.tot
self$E.l <- E.l
self$mean.strain <- mean(abs(self$getStrains()$spherical$strain))
self$mean.logstrain <- mean(abs(self$getStrains()$spherical$logstrain))
## Plot
if (plot.3d) {
sphericalplot(self, datapoints=FALSE, strain=FALSE)
}
if (!is.na(dev.flat)) {
dev.set(dev.flat)
flatplot(self, grid=TRUE, strain=TRUE, mesh=FALSE, markup=FALSE,
datapoints=FALSE, landmarks=FALSE,
image=FALSE)
}
if (!is.na(dev.polar)) {
## Wipe any previous reconstruction of coordinates of pixels and feature sets
private$ims <- NULL
self$clearFeatureSets()
dev.set(dev.polar)
projection(self, mesh=TRUE,
datapoints=FALSE, landmarks=FALSE,
image=FALSE)
}
}
},
##' @description Optimise the mapping from the flat outline to the sphere
##' @param alpha Area penalty scaling coefficient
##' @param x0 Area penalty cut-off coefficient
##' @param nu Power to which to raise area
##' @param method Method to pass to \code{optim}
##' @param plot.3d If \code{TRUE} make a 3D plot in an RGL window
##' @param dev.flat Device handle for plotting grid to
##' @param dev.polar Device handle for plotting polar plot to
##' @param ... Extra arguments to pass to \code{\link{fire}}
optimiseMappingCart = function(alpha=4, x0=0.5, nu=1, method="BFGS",
plot.3d=FALSE, dev.flat=NA, dev.polar=NA, ...) {
phi <- self$phi
lambda <- self$lambda
R <- self$R
phi0 <- self$phi0
lambda0 <- self$lambda0
Tt <- self$Tt
A <- self$ol$A
Cut <- self$Cut
Ct <- self$Ct
Lt <- self$Lt
Bt <- self$Bt
Rsett <- self$Rsett
i0t <- self$i0t
Nt <- nrow(self$Pt)
Nphi <- Nt - length(Rsett)
## Optimisation and plotting
opt <- list()
opt$x <- sphere.spherical.to.sphere.cart(cbind(phi=phi, lambda=lambda), R)
opt$conv <- 1
## Compute "mass" for each node
minL <- rep(Inf, nrow(self$Pt))
for (i in 1:nrow(Cut)) {
minL[Cut[i,1]] <- min(minL[Cut[i,1]], Lt[i])
minL[Cut[i,2]] <- min(minL[Cut[i,2]], Lt[i])
}
m <- 1/minL
m <- m/mean(m)
count <- 50
while (opt$conv && count) {
## Optimise
opt <- fire(opt$x,
force=function(x) {Fcart(x, Ct, Lt, Tt, A, R, alpha, x0, nu)},
restraint=function(x) {Rcart(x, R, Rsett, i0t, phi0, lambda0)},
dt=1,
nstep=200,
m=m, verbose=TRUE, report=report, ...)
count <- count - 1
## Report
E.tot <- Ecart(opt$x, Cu=Cut, L=Lt, R=R, T=Tt, A=A,
alpha=alpha, x0=x0, nu=nu)
E.l <- Ecart(opt$x, Cu=Cut, L=Lt, R=R, T=Tt, A=A,
alpha=0, x0=x0, nu=0)
s <- sphere.cart.to.sphere.spherical(opt$x, R)
phi <- s[,"phi"]
lambda <- s[,"lambda"]
ft <- flipped.triangles(cbind(phi=phi, lambda=lambda), Tt, R)
nflip <- sum(ft$flipped)
report(sprintf("E = %8.5f | E_L = %8.5f | E_A = %8.5f | %3d flippped triangles", E.tot, E.l, E.tot - E.l, nflip))
if (nflip) {
print(data.frame(rbind(id=which(ft$flipped),
A=A[ft$flipped],
a=ft$areas[ft$flipped])))
}
## Plot
if (plot.3d) {
sphericalplot(list(phi=phi, lambda=lambda, R=R,
Tt=Tt, Rsett=Rsett, gb=self$ol$gb, ht=self$ol$ht),
datapoints=FALSE)
}
if (!is.na(dev.flat)) {
dev.set(dev.flat)
flatplot(self, grid=TRUE, strain=TRUE, mesh=FALSE, markup=FALSE,
datapoints=FALSE, landmarks=FALSE,
image=FALSE)
}
if (!is.na(dev.polar)) {
## Wipe any previous reconstruction of coordinates of pixels and feature sets
private$ims <- NULL
self$clearFeatureSets()
dev.set(dev.polar)
self$phi <- phi
self$lambda <- lambda
projection(self, mesh=TRUE,
datapoints=FALSE, landmarks=FALSE,
image=FALSE)
}
}
self$phi <- phi
self$lambda <- lambda
self$opt <- opt
self$nflip <- sum(ft$flipped)
self$E.tot <- E.tot
self$E.l <- E.l
self$mean.strain <- mean(abs(self$getStrains()$spherical$strain))
self$mean.logstrain <- mean(abs(self$getStrains()$spherical$logstrain))
},
##' @description Transform an image into the reconstructed space
##' Transform an image into the reconstructed space. The four corner
##' coordinates of each pixel are transformed into spherical
##' coordinates and a mask matrix with the same dimensions as
##' \code{im} is created. This has \code{TRUE} for pixels that should
##' be displayed and \code{FALSE} for ones that should not.
##' Sets the field
##' \itemize{
##' \item{\code{ims}}{Coordinates of corners of pixels in spherical coordinates}
##' }
transformImage = function() {
im <- self$ol$getImage()
if (!is.null(im)) {
## Need to find the *boundaries* of pixels
N <- ncol(im)
M <- nrow(im)
## Create grid coords of corners of pixels. These run from the
## top left of the image down each column of the image.
xs <- 0:N
ys <- M:0
## x-coords of pixel corners, arranged in (N+1) by (M+1) grid
## Ditto for y-coords
I <- cbind(as.vector(outer(ys*0, xs, FUN="+")),
as.vector(outer(ys, xs*0, FUN="+")))
## Find Barycentric coordinates of corners of pixels
Ib <- tsearch(self$ol$getPoints()[,"X"], self$ol$getPoints()[,"Y"],
self$ol$T, I[,1], I[,2], bary=TRUE)
rm(I)
gc()
## Find 3D coordinates of mesh points
Pc <- sph2cart(theta=self$lambda, phi=self$phi, r=1)
private$ims <- bary2sph(Ib, self$Tt, Pc)
}
},
##' @description Get coordinates of corners of pixels of image in spherical
##' coordinates
##' @return Coordinates of corners of pixels in spherical coordinates
getIms = function() {
if (is.null(private$ims)) {
report("Transforming image...")
## Force garbage collection; not great practice, but this
## procedure is imemory intensive for large images
gc()
self$transformImage()
gc()
}
return(private$ims)
},
##' @description Get location of tear coordinates in spherical coordinates
##' @return Location of tear coordinates in spherical coordinates
getTearCoords = function() {
Tss <- list()
for (TF in self$ol$TFset) {
## Convert indices to the spherical frame of reference
j <- self$ht[TF]
Tss <- c(Tss, list(cbind(phi=self$phi[j], lambda=self$lambda[j])))
}
return(Tss)
},
##' @description Get \link{ReconstructedFeatureSet}
##' @param type Base type of \link{FeatureSet} as string.
##' E.g. \code{PointSet} returns a \link{ReconstructedPointSet}
getFeatureSet = function(type) {
type <- paste0("Reconstructed", type)
fs <- super$getFeatureSet(type)
if (is.null(fs)) {
self$reconstructFeatureSets()
fs <- super$getFeatureSet(type)
}
return(fs)
},
##' @description Reconstruct any attached feature sets.
reconstructFeatureSets = function() {
self$featureSets <- lapply(self$ol$getFeatureSets(), function(x) x$reconstruct(self))
},
##' @description Get mesh points in spherical coordinates
##' @return Matrix with columns \code{phi} (latitude) and \code{lambda}
##' (longitude)
getPoints = function() {
return(cbind(phi=self$phi, lambda=self$lambda))
},
##' @description Return location of point on sphere corresponding
##' to point on the flat outline
##' @param P Cartesian coordinates on flat outline as a matrix
##' with \code{X} and \code{Y} columns
mapFlatToSpherical = function(P) {
if (!(is.numeric(P))) {
stop("P must be numeric")
}
if (!(is.matrix(P))) {
stop("P must be matrix")
}
if (!(all(c("X", "Y") %in% colnames(P)))) {
stop("P should have columns named X and Y")
}
## Meshpoints in Cartesian coordinates
Ptc <- sph2cart(theta=self$lambda, phi=self$phi, r=1)
Pb <- tsearch(self$ol$getPoints()[,"X"],
self$ol$getPoints()[,"Y"],
self$ol$T,
P[,"X"],
P[,"Y"], bary=TRUE)
oo <- is.na(Pb$idx) # Points outwith outline
if (any(oo)) {
warning(paste(sum(oo), "points outwith the outline will be ignored"))
}
Pb$p <- Pb$p[!oo,,drop=FALSE]
Pb$idx <- Pb$idx[!oo]
return(bary2sph(Pb, self$Tt, Ptc))
}
)
)
##' Plot \code{\link{ReconstructedOutline}} object. This adds a mesh
##' of gridlines from the spherical retina (described by points
##' \code{phi}, \code{lambda} and triangulation \code{Tt} and cut-off
##' point \code{phi0}) onto a flattened retina (described by points
##' \code{P} and triangulation \code{T}).
##'
##' @title Flat plot of reconstructed outline
##' @param x \code{\link{ReconstructedOutline}} object
##' @param axt whether to plot axes
##' @param xlim x-limits
##' @param ylim y-limits
##' @param grid Whether or not to show the grid lines of
##' latitude and longitude
##' @param strain Whether or not to show the strain
##' @param ... Other plotting parameters
##' @method flatplot ReconstructedOutline
##' @author David Sterratt
##' @importFrom grDevices rainbow palette
##' @export
flatplot.ReconstructedOutline <- function(x, axt="n",
xlim=NULL, ylim=NULL,
grid=TRUE,
strain=FALSE,
...) {
## NextMethod()
flatplot(x$ol, ...)
if (strain) {
o <- x$getStrains()
palette(rainbow(100))
P <- x$ol$getPoints()
scols <- strain.colours(o$flat$logstrain)
Cu <- x$ol$Cu
segments(P[Cu[,1],1], P[Cu[,1],2],
P[Cu[,2],1], P[Cu[,2],2], col=round(scols))
}
## Plot a gridline from the spherical retina (described by points phi,
## lambda and triangulation Tt) onto a flattened retina (described by
## points P and triangulation T). The gridline is described by a
## normal n to a plane and a distance to the plane. The intersection of
## the plane and the spehere is the gridline.
get.gridline.flat <- function(P, T, phi, lambda, Tt, n, d, ...) {
mu <- compute.intersections.sphere(phi, lambda, Tt, n, d)
## Take out rows that are not intersections. If a plane intersects
## one edge of a triangle and the opposing vertex, in the row
## corresponding to the triangle, there will be a 0, a 1 and a
## value between 0 and 1. We get rid of the 1 in the
## following. Triangles in which one line is in the plane have mu
## values 0, 1 and NaN; we want to include these.
tri.int <- (rowSums((mu >= 0) & (mu <= 1), na.rm=TRUE) == 2)
## | apply(mu, 1, function(x) setequal(x, c(0, 1, NaN))))
if (any(tri.int)) {
T <- T[tri.int,,drop=FALSE]
mu <- mu[tri.int,,drop=FALSE]
## Create a logical matrix of which points are involved in lines
## that interscect the plane.
line.int <- (mu >= 0) & (mu < 1)
## If any element of mu contained a NaN, due to a line being in
## the plane, this should be set to false as the point opposite
## the NaN is not in the plane
line.int[is.na(line.int)] <- FALSE
## Order rows so that the false indicator is in the third column
T[!line.int[,2] ,] <- T[!line.int[,2], c(3,1,2)]
mu[!line.int[,2],] <- mu[!line.int[,2],c(3,1,2)]
T[!line.int[,1] ,] <- T[!line.int[,1], c(2,3,1)]
mu[!line.int[,1],] <- mu[!line.int[,1],c(2,3,1)]
P <- cbind(mu[,1] * P[T[,3],] + (1-mu[,1]) * P[T[,2],],
mu[,2] * P[T[,1],] + (1-mu[,2]) * P[T[,3],])
# suppressWarnings(segments(P1[,1], P1[,2], P2[,1], P2[,2], ...))
} else {
P <- matrix(0, nrow=1, ncol=4)
}
colnames(P) <- c("X1", "Y1", "X2", "Y2")
return(P)
}
if (grid) {
grid.int.minor <- 15
grid.int.major <- 45
grid.maj.col <- getOption("grid.maj.col")
grid.min.col <- getOption("grid.min.col")
phi0d <- x$phi0 * 180/pi
P <- matrix(0, nrow=0, ncol=4)
cols <- NULL
Phis <- setdiff(seq(-90, phi0d, by=grid.int.minor), phi0d)
Lambdas <- seq(0, 180-grid.int.minor, by=grid.int.minor)
for (Phi in Phis) {
if ((!(Phi %% grid.int.major) || Phi == phi0d)) {
col <- grid.maj.col
} else {
col <- grid.min.col
}
P1 <- get.gridline.flat(x$ol$getPoints(), x$ol$T, x$phi, x$lambda, x$Tt,
c(0,0,1), sin(Phi*pi/180))
cols <- c(cols, rep(col, nrow(P1)))
P <- rbind(P, P1)
}
for (Lambda in Lambdas) {
if (!(Lambda %% grid.int.major)) {
col <- grid.maj.col
} else {
col <- grid.min.col
}
Lambda <- Lambda * pi/180
P1 <- get.gridline.flat(x$ol$getPoints(), x$ol$T, x$phi, x$lambda, x$Tt,
c(sin(Lambda),cos(Lambda),0), 0)
cols <- c(cols, rep(col, nrow(P1)))
P <- rbind(P, P1)
}
if (nrow(P) > 0) {
segments(P[,"X1"], P[,"Y1"], P[,"X2"], P[,"Y2"], col=cols)
}
}
}
##' Draw a projection of a \code{\link{ReconstructedOutline}}. This method sets up
##' the grid lines and the angular labels and draws the image.
##'
##' @title Projection of a reconstructed outline
##' @param r \code{ReconstructedOutline} object
##' @param transform Transform function to apply to spherical coordinates
##' before rotation
##' @param axisdir Direction of axis (North pole) of sphere in
##' external space as matrix with column names \code{phi} (elevation)
##' and \code{lambda} (longitude).
##' @param projection Projection in which to display object,
##' e.g. \code{\link{azimuthal.equalarea}} or \code{\link{sinusoidal}}
##' @param proj.centre Location of centre of projection as matrix with
##' column names \code{phi} (elevation) and \code{lambda} (longitude).
##' @param lambdalim Limits of longitude (in degrees) to display
##' @param philim Limits of latitude (in degrees) to display
##' @param labels Vector of 4 labels to plot at 0, 90, 180 and 270 degrees
##' @param mesh If \code{TRUE}, plot mesh
##' @param grid Whether or not to show the grid lines of
##' latitude and longitude
##' @param grid.bg Background colour of the grid
##' @param grid.int.minor Interval between minor grid lines in degrees
##' @param grid.int.major Interval between major grid lines in degrees
##' @param colatitude If \code{TRUE} have radial labels plotted with
##' respect to colatitude rather than latitude
##' @param pole If \code{TRUE} indicate the pole with a "*"
##' @param image If \code{TRUE}, show the image
##' @param markup If \code{TRUE}, plot markup, i.e. reconstructed tears
##' @param add If \code{TRUE}, don't draw axes; add to existing plot.
##' @param max.proj.dim Maximum width of the image created in pixels
##' @param ... Graphical parameters to pass to plotting functions
##' @method projection ReconstructedOutline
##' @export
projection.ReconstructedOutline <- function(r,
transform=identity.transform,
axisdir=cbind(phi=90, lambda=0),
projection=azimuthal.equalarea,
proj.centre=cbind(phi=0, lambda=0),
lambdalim=c(-180, 180),
philim=c(-90, 90),
labels=c(0, 90, 180, 270),
mesh=FALSE,
grid=TRUE,
grid.bg="transparent",
grid.int.minor=15,
grid.int.major=45,
colatitude=TRUE,
pole=FALSE,
image=TRUE,
markup=TRUE,
add=FALSE,
max.proj.dim=getOption("max.proj.dim"),
...) {
Call <- match.call(expand.dots=TRUE)
plot.image <- image
## Compute grid lines
## Lines of latitude (parallels)
## Determine the major and minor parallels
phis.maj <- seq(-90, 90, by=grid.int.major)
phis.maj <- c(philim[1],
phis.maj[(phis.maj > philim[1]) & (phis.maj < philim[2])],
philim[2])
phis.min <- seq(-90, 90, by=grid.int.minor)
phis.min <- c(philim[0],
phis.min[(phis.min > philim[1]) & (phis.min < philim[2])],
philim[2])
phis.min <- setdiff(phis.min, phis.maj)
## Longitudes at which to draw lines; the smaller the by interval,
## the smoother
lambdas <- seq(lambdalim[1], lambdalim[2], by=1)
## Compute the minor and and major parallels to draw
paras.min <- projection(pi/180*cbind(phi =as.vector(outer(c(lambdas, NA)*0, phis.min, FUN="+")),
lambda=as.vector(outer(c(lambdas, NA), phis.min*0, FUN="+"))),
proj.centre=pi/180*proj.centre)
paras.maj <- projection(pi/180*cbind(phi =as.vector(outer(c(lambdas, NA)*0, phis.maj, FUN="+")),
lambda=as.vector(outer(c(lambdas, NA), phis.maj*0, FUN="+"))),
proj.centre=pi/180*proj.centre)
## Lines of longitude (meridians)
## Determine the major and minor parallels
lambdas.maj <- c(rev(seq(0 , lambdalim[1], by=-grid.int.major)),
seq(grid.int.major, lambdalim[2], by= grid.int.major))
lambdas.min <- c(rev(seq(0 , lambdalim[1], by=-grid.int.minor)),
seq(grid.int.minor, lambdalim[2], by= grid.int.minor))
lambdas.min <- setdiff(lambdas.min, lambdas.maj)
## Latitudes at which to draw lines; the smaller the by interval,
## the smoother
phis <- seq(philim[1], philim[2], by=1)
## Compute the minor and and major meridians to draw
merids.min <- projection(pi/180*cbind(phi =as.vector(outer(c(phis, NA), lambdas.min*0, FUN="+")),
lambda=as.vector(outer(c(phis*0, NA), lambdas.min, FUN="+"))),
proj.centre=pi/180*proj.centre)
merids.maj <- projection(pi/180*cbind(phi =as.vector(outer(c(phis, NA), lambdas.maj*0, FUN="+")),
lambda=as.vector(outer(c(phis*0, NA), lambdas.maj, FUN="+"))),
proj.centre=pi/180*proj.centre)
## Set up the plot region
xlim <- range(na.omit(rbind(paras.min, paras.maj))[,"x"])
ylim <- range(na.omit(rbind(paras.min, paras.maj))[,"y"])
if (!add) {
plot(NA, NA, xlim=xlim, ylim=ylim,# xaxs="i", yaxs="i",
type = "n", axes = FALSE, xlab = "", ylab = "", asp=1,
main=Call$main, bg=Call$bg)
}
## Plot an image
## Get the spherical coordinates of the corners of pixels
if (plot.image) {
ims <- r$getIms()
gc()
if (!is.null(ims)) {
## Reconstitute image from stored values of phi and lambda
## coordinates of corners of pixels
## Get the size of the image
im <- r$ol$getImage()
M <- nrow(im)
N <- ncol(im)
## Downsample the image by first selecting rows and columns to
## look at
by <- ceiling(max(N, M)/max.proj.dim) # Number of pixels to merge
Ms <- seq(1, M - (M %% by), by=by)
Ns <- seq(1, N - (N %% by), by=by)
## Downsample the image
## This should not create a new version of the image
if (by > 1) {
report("Downsampling image by factor of ", by)
im <- im[Ms, Ns]
}
## Now need to do the more complex job of downsampling the matrix
## containing the coordinates of the corners of pixels
if (by > 1) {
report("Downsampling pixel corner spherical coordinates by factor of ", by)
imsmask <- matrix(FALSE, M+1, N+1)
imsmask[c(Ms, (max(Ms) + by)), c(Ns, (max(Ns) + by))] <- TRUE
ims <- ims[imsmask,]
}
## Convenience variables for the new image sizes
M <- nrow(im)
N <- ncol(im)
## Target size of of block, i.e. the number of pixels to
## transform in one go
S <- 200*200 # 200*8000
## Number of blocks
B <- ceiling(N*M/S)
## Numnber of columns in the first B-1 blocks
C <- floor(N/B)
## In the Bth block there will be N-(B-1)*C columns
## print(paste("M =", M, "; N =", N, "; S =", S, "; B =", B, "; C =", C))
## Number of columns in block k
Ck <- C
for (k in 1:B) {
report("Projecting block ", k, "/", B)
## Actual number of columns, since the last block may have a
## different number of chunks to C
if (k == B) {
Ck <- N - (B - 1)*C
}
## Index of leftmost column of image matrix needed
j0 <- (k - 1)*C + 1
## Index of rightmost column of image matrix needed
j1 <- (k - 1)*C + Ck
## Transform the pixel coordinates and compute x and y positions
## of corners of pixels.
## inds <- (k - 1)*(M + 1)*C + (1:((M + 1)*(Ci + 1)))
l0 <- (j0 - 1)*(M + 1) + 1
l1 <- j1*(M + 1) + (M + 1)
## print(paste("k =", k, "; Ck =", Ck, "; j0 =", j0, "; j1 =", j1,
## "; l0 =", l0, "; l1 =", l1))
inds <- l0:l1
rc <- projection(
rotate.axis(
transform(ims[inds,],
phi0=r$phi0),
axisdir*pi/180),
lambdalim=lambdalim*pi/180,
proj.centre=pi/180*proj.centre)
xpos <- matrix(rc[,"x"], M + 1, Ck + 1)
ypos <- matrix(rc[,"y"], M + 1, Ck + 1)
## Convert these to format suitable for polygon()
impx <- rbind(as.vector(xpos[1:M , 1:Ck ]),
as.vector(xpos[1:M , 2:(Ck+1)]),
as.vector(xpos[2:(M+1), 2:(Ck+1)]),
as.vector(xpos[2:(M+1), 1:Ck ]),
NA)
impy <- rbind(as.vector(ypos[1:M , 1:Ck ]),
as.vector(ypos[1:M , 2:(Ck+1)]),
as.vector(ypos[2:(M+1), 2:(Ck+1)]),
as.vector(ypos[2:(M+1), 1:Ck ]),
NA)
## Pixels outside the image should be masked. The mask is a matrix
## the same size as the image, containing TRUE for pixels that
## should be displayed and FALSE for those that should be
## masked. It is computed by finding the corners of the poly-pixel
## lie outwith the outline. These corners will have the coordinate
## NA. print("sum(!is.na(colSums(impx[1:4,])))")
## print(sum(!is.na(colSums(impx[1:4,]))))
immask <- matrix(!is.na(colSums(impx[1:4,])), M, Ck)
## We want to get rid of any poly-pixels that cross either end of
## the longitude range in a pseudocylindrical projection. A simple
## way of doing this is to say that if a pixel is very large,
## don't plot it.
## This code is chronically slow, hence replacing with the
## pmin/pmax version
## bigpx <- which(apply(impx[1:4,], 2,
## function(x) {max(x) - min(x)}) > 0.1*abs(diff(xlim)) |
## apply(impy[1:4,], 2,
## function(y) {max(y) - min(y)}) > 0.1*abs(diff(ylim)))
bigpx <- which((pmax(impx[1,], impx[2,], impx[3,], impx[4,]) -
pmin(impx[1,], impx[2,], impx[3,], impx[4,]) >
0.1*abs(diff(xlim))) |
(pmax(impy[1,], impy[2,], impy[3,], impy[4,]) -
pmin(impy[1,], impy[2,], impy[3,], impy[4,]) >
0.1*abs(diff(xlim))))
immask[bigpx] <- FALSE
## Plot the polygon, masking as we go
graphics::polygon(impx[,immask], impy[,immask],
col=im[1:M,(k-1)*C+(1:Ck)][immask],
border=im[1:M,(k-1)*C+(1:Ck)][immask])
}
}
}
## Plot the mesh
if (mesh) {
Pt <- projection(rotate.axis(transform(cbind(phi=r$phi, lambda=r$lambda),
phi0=r$phi0),
axisdir*pi/180),
lambdalim=lambdalim*pi/180,
proj.centre=pi/180*proj.centre)
trimesh(r$Tt, Pt, col="gray", add=TRUE)
}
grid.maj.col <- getOption("grid.maj.col")
grid.min.col <- getOption("grid.min.col")
## Plot the grid
if (grid) {
## Minor paralells and meridians
lines(paras.min, col=grid.min.col)
lines(merids.min, col=grid.min.col)
## Major lines of latitude on top of all minor lines
lines(paras.maj, col=grid.maj.col)
lines(merids.maj, col=grid.maj.col)
## Boundary of projection
boundary <- projection("boundary")
graphics::polygon(boundary[,"x"], boundary[,"y"], border="black")
}
## Plot rim in visutopic space
rs <- cbind(phi=r$phi0, lambda=seq(0, 2*pi, len=360))
rs.rot <- rotate.axis(transform(rs, phi0=r$phi0), axisdir*pi/180)
## "Home" position for a cyclops looking ahead
## r$axisdir = cbind(phi=0, lambda=0)
lines(projection(rs.rot, lambdalim=lambdalim*pi/180, lines=TRUE,
proj.centre=pi/180*proj.centre),
col=getOption("TF.col"))
## Projection of pole
if (pole) {
oa.rot <- rotate.axis(transform(cbind(phi=-pi/2, lambda=0), phi0=r$phi0),
axisdir*pi/180)
points(projection(oa.rot, proj.centre=pi/180*proj.centre),
pch="*", col=getOption("TF.col"), cex=2)
}
## Plot tears
if (markup) {
Tss <- r$getTearCoords()
for (Ts in Tss) {
## Plot
lines(projection(rotate.axis(transform(Ts, phi0=r$phi0),
axisdir*pi/180),
lines=TRUE,
lambdalim=lambdalim*pi/180,
proj.centre=pi/180*proj.centre),
col=getOption("TF.col"), lwd=Call$lwd, lty=Call$lty)
}
}
if (grid) {
## Longitude labels around rim - not on actual frame of reference!
if (!is.null(labels)) {
## Longitudes (meridians) at which to plot at
angles <- seq(0, by=2*pi/length(labels), len=length(labels))
## This is a nasty hack: we want to find out how far to plot the
## labels from the rim. We can't do this simply in terms of
## angles, so we have to find the right angular distance to give
## the desired fraction of the axes at which to plot the
## labels. This is done by this optimisation function.
label.fax <- 0.02 # Fraction of axis length from axes to plot labels
opt <- stats::optimise(function(a) {
rs0 <- cbind(phi=r$phi0, lambda=angles[1])
rs <- cbind(phi=r$phi0 + a, lambda=angles[1])
rc0 <- projection(rotate.axis(transform(rs0, phi0=r$phi0),
axisdir*pi/180),
proj.centre=pi/180*proj.centre)
rc <- projection(rotate.axis(transform(rs, phi0=r$phi0),
axisdir*pi/180),
proj.centre=pi/180*proj.centre)
return((vecnorm(rc - rc0) - label.fax*abs(diff(xlim)))^2)
}
,interval=c(1, 20)*pi/180)
lambda.label.off <- opt$minimum
## Now plot the labels themselves. Phew!!
rs <- cbind(phi=r$phi0 + lambda.label.off, lambda=angles)
rc <- projection(rotate.axis(transform(rs, phi0=r$phi0), axisdir*pi/180),
proj.centre=pi/180*proj.centre)
text(rc[,"x"], rc[,"y"], labels, xpd=TRUE)
}
## Latitude Labels
## rlabels <- c(seq(philim[1], philim[2], by=grid.int.major))
rlabels <- phis.maj
rs <- cbind(phi=rlabels*pi/180, lambda=proj.centre[1,"lambda"])
rc <- projection(rs, proj.centre=pi/180*proj.centre)
text(rc[,"x"], rc[,"y"], rlabels + ifelse(colatitude, 90, 0),
xpd=TRUE, adj=c(1, 1), col=grid.maj.col)
}
}
##' @export
##' @importFrom graphics abline
lvsLplot.ReconstructedOutline <- function(r, ...) {
Call <- match.call(expand.dots=TRUE)
o <- r$getStrains()$spherical
palette(rainbow(100)) ## Green is about 35; dark blue about 70
cols <- strain.colours(o$logstrain)
L <- o$L
l <- o$l
units <- r$ol$units
xlab <- expression(paste("Flat ", italic(L)))
ylab <- expression(paste("Reconstructed ", italic(l)))
if (!is.na(units)) {
xlab <- eval(substitute(expression(paste("Flat ", italic(L), " (", units1, ")")),
list(units1=units)))
ylab <- eval(substitute(expression(paste("Reconstructed ", italic(l), " (",
units1, ")")),
list(units1=units)))
}
plot(L, l, col=cols, pch=20,
xlim=c(0, max(L, l)), ylim=c(0, max(L, l)),
xlab=xlab, ylab=ylab,
asp=1, main=Call$main)
par(xpd=FALSE)
abline(0, 1)
abline(0, 0.75, col="blue")
abline(0, 1.25, col="red")
text(0.2*max(L), 0.2*max(L)*0.5, "25% compressed", col="blue",
pos=4)
text(0.75*max(L), 0.75*max(L)*1.25, "25% expanded", col="red",
pos=2)
}
##' Draw a spherical plot of reconstructed outline. This method just
##' draws the mesh.
##'
##' @title Spherical plot of reconstructed outline
##' @param r \code{\link{ReconstructedOutline}} object
##' @param strain If \code{TRUE}, plot the strain
##' @param surf If \code{TRUE}, plot the surface
##' @param ... Other graphics parameters -- not used at present
##' @method sphericalplot ReconstructedOutline
##' @author David Sterratt
##' @import rgl
##' @export
sphericalplot.ReconstructedOutline <- function(r,
strain=FALSE,
surf=TRUE, ...) {
## Obtain Cartesian coordinates of points
Ps <- r$getPoints()
P <- sphere.spherical.to.sphere.cart(Ps)
rgl.clear()
if (surf) {
## Outer triangles
fac <- 1.005
triangles3d(matrix(fac*P[t(r$Tt[,c(2,1,3)]),1], nrow=3),
matrix(fac*P[t(r$Tt[,c(2,1,3)]),2], nrow=3),
matrix(fac*P[t(r$Tt[,c(2,1,3)]),3], nrow=3),
color="darkgrey", alpha=1)
## Inner triangles
triangles3d(matrix(P[t(r$Tt),1], nrow=3),
matrix(P[t(r$Tt),2], nrow=3),
matrix(P[t(r$Tt),3], nrow=3),
color="white", alpha=1)
}
## Plot any flipped triangles
ft <- flipped.triangles(Ps, r$Tt)
with(ft, points3d(cents[flipped,1], cents[flipped,2], cents[flipped,3],
col="blue", size=5))
## Shrink so that they appear inside the hemisphere
fac <- 0.997
ht <- r$ht
gb <- r$ol$gb
rgl.lines(fac*rbind(P[ht[gb[gb]],1], P[ht[gb],1]),
fac*rbind(P[ht[gb[gb]],2], P[ht[gb],2]),
fac*rbind(P[ht[gb[gb]],3], P[ht[gb],3]),
lwd=3, color=getOption("TF.col"))
fac <- 1.006
rgl.lines(fac*rbind(P[ht[gb[gb]],1], P[ht[gb],1]),
fac*rbind(P[ht[gb[gb]],2], P[ht[gb],2]),
fac*rbind(P[ht[gb[gb]],3], P[ht[gb],3]),
lwd=3, color=getOption("TF.col"))
if (strain) {
o <- r$getStrains()
palette(rainbow(100))
scols <- strain.colours(o$spherical$logstrain)
fac <- 0.999
P1 <- fac*P[r$Cut[,1],]
P2 <- fac*P[r$Cut[,2],]
width <- 0.02
## Compute displacement vector to make sure that strips are
## parallel to surface of sphere
d <- extprod3d(P1, P2-P1)
d <- width/2*d/vecnorm(d)
PA <- P1 - d
PB <- P1 + d
PC <- P2 + d
PD <- P2 - d
## This is a ridiculously inefficient way of drawing the strain,
## but if you try presenting a color vector, it makes each line
## multi-coloured. It has taking HOURS of fiddling round to
## discover this! GRRRRRRRRRRRRRRRRR!
for (i in 1:nrow(PA)) {
quads3d(rbind(PA[i,1], PB[i,1], PC[i,1], PD[i,1]),
rbind(PA[i,2], PB[i,2], PC[i,2], PD[i,2]),
rbind(PA[i,3], PB[i,3], PC[i,3], PD[i,3]),
color=round(scols[i]), alpha=1)
}
}
}
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.