# Re: [R] Fitting a half-ellipse curve

From: Niklaus Hurlimann <niklaus.hurlimann_at_unil.ch>
Date: Thu, 30 Sep 2010 16:54:57 +0200

Hello Michael,

Thanks very much for the hint to my problem this helps really a lot.

Cheers
Niklaus

```--
Niklaus H√ľrlimann
Doctorant-PhD

L'Anthropole
CH-1015 Lausanne
Suisse

E-mail: Niklaus.Hurlimann_at_unil.ch
Tel:+41(0)21 692 4452

attached mail follows:

Hello Niklaus,

I'm not sure if the following is the sort of thing you are looking for (?)

You can fit an ellipse to your data using a deterministic least
squares method. The following is a function that I use to do this...

fit.ellipse <- function (x, y = NULL)
{

# Least squares fitting of an ellipse to point data

# using the algorithm described in:
#   Radim Halir & Jan Flusser. 1998.
#   Numerically stable direct least squares fitting of ellipses.
#   Proceedings of the 6th International Conference in Central Europe
#   on Computer Graphics and Visualization. WSCG '98, p. 125-132
#
# Adapted from the original Matlab code by Michael Bedward
# <michael.bedward_at_gmail.com>
#
# Arguments:
# x, y - the x and y coordinates of the data points
#
# Returns a list with the following elements:
#
# coef - coefficients of the ellipse as described by the general
#        quadratic:  ax^2 + bxy + cy^2 + dx + ey + f = 0
#
# center - center x and y
#
# major - major semi-axis length
#
# minor - minor semi-axis length
#

EPS <- 1.0e-8
dat <- xy.coords(x, y)

D1 <- cbind(dat\$x * dat\$x, dat\$x * dat\$y, dat\$y * dat\$y)
D2 <- cbind(dat\$x, dat\$y, 1)
S1 <- t(D1) %*% D1
S2 <- t(D1) %*% D2
S3 <- t(D2) %*% D2
T <- -solve(S3) %*% t(S2)
M <- S1 + S2 %*% T
M <- rbind(M[3,] / 2, -M[2,], M[1,] / 2)
evec <- eigen(M)\$vec
cond <- 4 * evec[1,] * evec[3,] - evec[2,]^2
a1 <- evec[, which(cond > 0)]
f <- c(a1, T %*% a1)
names(f) <- letters[1:6]

# calculate the center and lengths of the semi-axes
b2 <- f^2 / 4
center <- c((f * f / 2 - b2 * f), (f * f / 2 - f *
f / 4)) / (b2 - f * f)
names(center) <- c("x", "y")

num <- 2 * (f * f^2 / 4 + f * f^2 / 4 + f * b2 -
f*f*f/4 - f*f*f)
den1 <- (b2 - f*f)
den2 <- sqrt((f - f)^2 + 4*b2)
den3 <- f + f

semi.axes <- sqrt(c( num / (den1 * (den2 - den3)),  num / (den1 *
(-den2 - den3)) ))

# calculate the angle of rotation
term <- (f - f) / f
angle <- atan(1 / term) / 2

list(coef=f, center = center, major = max(semi.axes), minor =
min(semi.axes), angle = unname(angle))
}

There are quite a few functions / packages to draw ellipses in R, but
the following is will work directly with the output of the above
function...

get.ellipse <- function ( fit, n=360 )
{
# Calculate points on an ellipse described by
# the fit argument as returned by fit.ellipse
#
# n is the number of points to render

tt <- seq(0, 2*pi, length=n)
sa <- sin(fit\$angle)
ca <- cos(fit\$angle)
ct <- cos(tt)
st <- sin(tt)

x <- fit\$center + fit\$maj * ct * ca - fit\$min * st * sa
y <- fit\$center + fit\$maj * ct * sa + fit\$min * st * ca
cbind(x=x, y=y)
}

So if your data were in a matrix or data.frame X...

efit <- fit.ellipse( X )
e <- get.ellipse( efit )
plot(X)
lines(e, col="red")

Hope this helps,
Michael

On 29 September 2010 23:45, Niklaus Hurlimann <niklaus.hurlimann_at_unil.ch> wrote:
> Dear mailing list,
>
> I have following array:
>
> † † † X2 † † † † † † † † Y2
> [1,] 422.7900 † † †6.0
> [2,] 469.8007 † † †10.5
> [3,] 483.9428 † † †11.0
> [4,] 532.4917 † † †25.5
> [5,] 596.1942 † † †33.5
> [6,] 630.8496 † † †40.5
> [7,] 733.2996 † † †45.0
> [8,] 946.4779 † † †32.0
> [9,] 996.8068 † † †35.5
> [10,] 1074.3310 †23.0
>
> I do afterwards the following:
>
> plot.new()
>
> plot.window(xlim=c(min(X1)-50,max(X1)+50),
> ylim=c(min(Y1)-25,max(Y1)+25))
>
> axis(2, cex.axis=1.2)
> axis(1, cex.axis=1.2)
>
> points(X2, Y2, type="p", pch=0, cex=1.2, col="black")
>
> title(main="Dyke Geometry Along Strike", cex.main=1.2, font.main=4)
> title(xlab="distance [m]", cex.lab=1.2)
> title(ylab="half-thickness [m]", cex.lab=1.2)
>
> box()
>
>
> I would like curve fitting where I proceeded with a non
> linear-regression with the according function below:
>
> nls(formula = Y2 ~ abs(b*abs(1-X2^2/a^2)^(1/2)), start = list( a=282,
> b=40))
>
> The formula should give the y-positive part of an ellipse in my opinion
> or I might be completely wrong.
>
> In the end I would like to fit a curve of half an ellipse through the
> data. †I might could do this as well with a 2nd order polynomial fit. I
> am grateful for any suggestions and comments to my problem.
>
>
> Cheers
>
>
>
> --
> Niklaus HŁrlimann
> Doctorant-PhD
>
> Universitť de Lausanne
> Institut de Minťralogie et Gťochimie
> L'Anthropole
> CH-1015 Lausanne
> Suisse
>
> E-mail: Niklaus.Hurlimann_at_unil.ch
> Tel:+41(0)21 692 4452
>
> † † † †[[alternative HTML version deleted]]
>
>
> ______________________________________________
> R-help_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> and provide commented, minimal, self-contained, reproducible code.
>
>

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help