[R] drawing ellipses

[Martin Maechler]

>>>>> "Yves" == Yves Brostaux <brostaux.y at fsagx.ac.be> writes:

    Yves> Hello again,

    Yves> First I want to thank all the people who answered my
    Yves> question about line width in graphs. I promise I will
    Yves> learn the 'par' help page by heart for the end of the
    Yves> month !


    Yves> I now want to trace some ellipses to emphasize groups
    Yves> of data. I found how to trace circles with
    Yves> 'symbols()', but no ellipse. I'm planning on writing
    Yves> my own function based on 'polygon()' and the ellipse
    Yves> equation. Does anybody already have done similar work
    Yves> or have another solution ?

yes. several.

1) The ellipse  package  does some.
2) The cluster  package (which is recommended, hence should be
		        available everywhere)
    has an ellipsoidhull() function and a plot method for its results
    (in the 2D case, i.e when the ellipsoid is an ellipse)

3) A `client' (researcher in another department here) needed
   something  based on geometric parameters rather than
   covariance matrices.
   So I ended up doing the following, I hope it helps.

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><

-------------------------------------------------------------------

## Martin's Solution : radially equispaced ("radial gleichmässig") :

ellipsePoints <- function(a,b, alpha = 0, loc = c(0,0), n = 201)
{
    ## Purpose: ellipse points,radially equispaced, given geometric par.s
    ## -------------------------------------------------------------------------
    ## Arguments: a, b : length of half axes in (x,y) direction
    ##            alpha: angle (in degrees) for rotation
    ##            loc  : center of ellipse
    ##            n    : number of points
    ## -------------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 19 Mar 2002, 16:26

    B <- min(a,b)
    A <- max(a,b)
    ## B <= A
    d2 <- (A-B)*(A+B)                   #= A^2 - B^2
    phi <- 2*pi*seq(0,1, len = n)
    sp <- sin(phi)
    cp <- cos(phi)
    r <- a*b / sqrt(B^2 + d2 * sp^2)
    xy <- r * cbind(cp, sp)
    ## xy are the ellipse points for alpha = 0 and loc = (0,0)
    al <- alpha * pi/180
    ca <- cos(al)
    sa <- sin(al)
    xy %*% rbind(c(ca, sa), c(-sa, ca)) + cbind(rep(loc[1],n),
                                                rep(loc[2],n))
}

## Example:

ep <- ellipsePoints(2,5)
str(ep)
plot(ep, type="n",asp=1) ; polygon(ep, col = 2)

## Movie : rotating ellipse  :
nTurns <- 4 # #{full 360 deg turns}
for(al in 1:(nTurns*360)) {
    ep <- ellipsePoints(3,6, alpha=al, loc = c(5,2))
    plot(ep,type="l",xlim=c(-1,11),ylim=c(-4,8),
         asp=1, axes = FALSE, xlab="", ylab="")
}

## Movie : rotating _filled_ ellipse  :
for(al in 1:180) {
    ep <- ellipsePoints(3,6, alpha=al, loc = c(5,2))
    plot(ep,type="n",xlim=c(-1,11),ylim=c(-4,8),
         asp=1, axes = FALSE, xlab="", ylab="")
    polygon(ep,col=2,border=3,lwd=2.5)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._