# R/ell.R In p3d: Plot 3d data and fitted surfaces

#### Documented in centercenter.ellConjCompdelldellellellboxell.conjellptellptcelltanelltancuvuv.defaultuv.ell

```##
## p3d:ell.R
## 2011-12-22
##

## Consolidated ell.R, dell.R, ell.conj.R, center.R and center.ell.F
## Added function ellpt, ellptc, elltan, elltanc, ellbox to generate points on ellipses,
## axes of ellipses and tangents.

ell <-
function( center, shape , radius  = 1, n =100) {
fac <- function( x )  {
# fac(M) is a 'right factor' of a positive semidefinite M
# i.e. M = t( fac(M) ) %*% fac(M)
# similar to chol(M) but does not require M to be PD.
xx <- svd(x)
t(xx\$v) * sqrt(pmax( xx\$d,0))
}
angles = (0:n) * 2 * pi/n
if ( length(radius) > 1) {
ret <- lapply( radius, function(r) rbind(r*cbind( cos(angles), sin(angles)),NA))
circle <- do.call( rbind, ret)
}
else circle = radius * cbind( cos(angles), sin(angles))
ret <- t( c(center) + t( circle %*% fac(shape)))
attr(ret,"parms") <- list ( center = rbind( center), shape = shape, radius = radius)
class(ret) <- "ell"
ret
}

dell <-
function( x, y, radius = 1, ...) {
if ( (is.matrix(x) && (ncol(x) > 1))|| is.data.frame(x)) mat <- as.matrix(x[,1:2])
else if (is.list(x)) mat <- cbind(x\$x, x\$y)
else mat <- cbind( x,y)
}

ell.conj <-
function( center, shape, dir, radius = 1, len = 1) {
# returns conjugate axes or tangent lines to ellipse
vecs <- uv( shape, dir, radius)
list( u = list( center-len*vecs\$u, center+len*vecs\$u),
v = list( center-len*vecs\$v, center+len*vecs\$v),
tan1 = list( center + vecs\$u - len*vecs\$v,center + vecs\$u+ len*vecs\$v),
tan2 = list( center - vecs\$u - len*vecs\$v,center - vecs\$u+ len*vecs\$v),
tan3 = list( center + vecs\$v - len*vecs\$u,center + vecs\$v+ len*vecs\$u),
tan4 = list( center - vecs\$v - len*vecs\$u,center - vecs\$v+ len*vecs\$u),
center = center)
}

center <-
function( obj, ... ) UseMethod("center")

center.ell <-
function( obj, ...) attr(obj, 'parms') \$ center

ConjComp <- function( X , Z = diag( nrow(X)) , ip = diag( nrow(X)), tol = 1e-07 ) {
# also in package spida:
# keep versions consistent
# ConjComp returns a basis for the conjugate complement of the
# conjugate projection of X into span(Z) with respect to inner product with
# matrix ip.
# Note: Z is assumed to be of full column rank but not necessarily X.
xq <- qr(t(Z) %*% ip %*% X, tol = tol)
if ( xq\$rank == 0 ) return( Z )
a <- qr.Q( xq, complete = TRUE ) [ ,-(1:xq\$rank)]
Z %*% a
}

uv <- function(object,...) UseMethod('uv')

uv.ell <- function( object, u, radius = 1, ...){
p <- attr(object,"parms")
}
uv.default <-
function( object, u , radius = 1, ...) {
# returns 'unit' u and conjugate v
u <- u / sqrt( sum( u*solve(object,u)))   # 'unit' vector in direction of dir
v <- c(ConjComp( u, diag(2) , solve(object)))  # conjugate
v <- v / sqrt( sum( v * solve( object, v)))
}

ellpt <- function(ell, dir = c(0,1) , radius = 1 ) {
# point on an ellipse in a particular direction
p <- attr(ell,'parms')
dir <- cbind(dir)
dn <- sum(dir* solve(p\$radius[1]^2 * p\$shape, dir))
ax <- dir / sqrt(dn)
#disp(ax)
#disp( p\$center)
t( ax %*% rbind(radius) + as.vector(p\$center))               # returns a row for plotting
}

ellptc <- function(ell, dir = c(0,1) , radius = 1 ) {
# point on an ellipse in a conjugate direction
p <- attr(ell,'parms')
dir <- cbind(dir)
ax <- Null( solve(p\$shape, dir) )
dn <- sum(ax* solve( p\$shape, ax))
t( ax %*% rbind(radius) + as.vector(p\$center))               # returns a row for plotting
}

elltanc <- function( ell, dir = c(0,1), radius = 1, len = 1, v = c(-1,1)) {
p <- attr(ell,'parms')
ax <- ellptc( ell, dir = dir, radius = len * v)
if ( is.null(radius) ) return( t( t(ax) + dir))
t( t(ax) -as.vector(p\$center) + as.vector(pt))
}

elltan <- function( ell, dir = c(0,1), radius = 1, len = 1, v = c(-1,1)) {
p <- attr(ell,'parms')
ax <- ellpt( ell, dir = dir, radius =  len * v)
t( t(ax) -as.vector(p\$center) + as.vector(pt))
}

ellbox <- function( ell, dir = c(0,1) , radius = 1 ){
rbind(
}

ellplus <-
function (center = rep(0, 2), shape = diag(2), radius = 1, n = 100,
angles = (0:n) * 2 * pi/n, fac = chol, ellipse = all, diameters = all,
box = all, all = FALSE)
{
help <- "\n        ellplus can produce, in addition to the points of an ellipse, the\n        conjugate axes corresponding to a chol or other decomposition\n        and the surrounding parallelogram.\n        "
rbindna <- function(x, ...) {
if (nargs() == 0)
return(NULL)
if (nargs() == 1)
return(x)
rbind(x, NA, rbindna(...))
}
if (missing(ellipse) && missing(diameters) && missing(box))
all <- TRUE
circle <- function(angle) cbind(cos(angle), sin(angle))
Tr <- fac(shape)
ret <- list(t(c(center) + t(radius * circle(angles) %*% Tr)),
t(c(center) + t(radius * circle(c(0, pi)) %*% Tr)), t(c(center) +
t(radius * circle(c(pi/2, 3 * pi/2)) %*% Tr)), t(c(center) +
t(radius * rbind(c(1, 1), c(-1, 1), c(-1, -1), c(1,
-1), c(1, 1)) %*% Tr)))
do.call("rbindna", ret[c(ellipse, diameters, diameters, box)])
}
```

## Try the p3d package in your browser

Any scripts or data that you put into this service are public.

p3d documentation built on May 2, 2019, 5:25 p.m.