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

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 

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]^2 / 4 center <- c((f[3] * f[4] / 2 - b2 * f[5]), (f[1] * f[5] / 2 - f[2] * f[4] / 4)) / (b2 - f[1] * f[3]) names(center) <- c("x", "y") num <- 2 * (f[1] * f[5]^2 / 4 + f[3] * f[4]^2 / 4 + f[6] * b2 - f[2]*f[4]*f[5]/4 - f[1]*f[3]*f[6]) den1 <- (b2 - f[1]*f[3]) den2 <- sqrt((f[1] - f[3])^2 + 4*b2) den3 <- f[1] + f[3] semi.axes <- sqrt(c( num / (den1 * (den2 - den3)), num / (den1 * (-den2 - den3)) )) # calculate the angle of rotation term <- (f[1] - f[3]) / f[2] 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[1] + fit$maj * ct * ca - fit$min * st * sa y <- fit$center[2] + 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 Hrlimann > Doctorant-PhD > > Universit de Lausanne > Institut de Minralogie et Gochimie > 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 > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >

______________________________________________ R-help_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

Received on Thu 30 Sep 2010 - 14:57:08 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 30 Sep 2010 - 15:00:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive