Anpassen von Datenpunkten an eine Ellipse mit Mittelpunkt im Ursprung mit R

Question

Ich habe eine Frage zum Anpassen von Ellipsen an Daten mit dem Ellipsenzentrum im Ursprung. Ich habe zwei Methoden untersucht, die auf Ellipsen passen, aber einen willkürlichen Mittelpunkt erzeugen, wenn ich die Daten nicht mit einigen imaginären Spiegelpunkten manipuliere.

Methode # 01

Dieser Teil des Skripts kommt direkt aus this nützlichen Beitrag. Ich kopiere die Codes direkt hier.

fit.ellipse <- function (x, y = NULL) { 
    # from: 
    # http://r.789695.n4.nabble.com/Fitting-a-half-ellipse-curve-tp2719037p2720560.html 
    # 
    # 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 (2010) 
    # michael.bedward@gmail.com 
    # 
    # Subsequently improved by John Minter (2012) 
    # 
    # Arguments: 
    # x, y - x and y coordinates of the data points. 
    #  If a single arg is provided it is assumed to be a 
    #  two column matrix. 
    # 
    # 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 
    # 
    # see http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2288654/ 
    # J. R. Minter 
    # for the center, linear algebra to the rescue 
    # center is the solution to the pair of equations 
    # 2ax + by + d = 0 
    # bx + 2cy + e = 0 
    # or 
    # | 2a b | |x| |-d| 
    # | b 2c | * |y| = |-e| 
    # or 
    # A x = b 
    # or 
    # x = Ainv b 
    # or 
    # x = solve(A) %*% b 
    A <- matrix(c(2*f[1], f[2], f[2], 2*f[3]), nrow=2, ncol=2, byrow=T) 
    b <- matrix(c(-f[4], -f[5]), nrow=2, ncol=1, byrow=T) 
    soln <- solve(A) %*% b 

    b2 <- f[2]^2/4 

    center <- c(soln[1], soln[2]) 
    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)) 
} 

Nehmen wir ein Beispiel Verteilung von polaren Punkte für Illustrationszwecke

X<-structure(list(x_polar = c(0, 229.777200000011, 246.746099999989, 
-10.8621999999741, -60.8808999999892, 75.8904999999795, -83.938199999975, 
-62.9770000000135, 49.1650999999838, 52.3093000000226, 49.6891000000178, 
-66.4248999999836, 34.3671999999788, 242.386400000018, 343.60619999998 
), y_polar = c(0, 214.868299999973, 161.063599999994, -68.8972000000067, 
-77.0230000000447, 93.2863000000361, -16.2356000000145, 27.7828000000445, 
-17.8077000000048, 2.10540000000037, 25.6866000000155, -84.6034999999683, 
-31.1800000000512, 192.010800000047, 222.003700000001)), .Names = c("x_polar", 
"y_polar"), row.names = c(NA, -15L), class = "data.frame") 


efit <- fit.ellipse(X) 
e <- get.ellipse(efit) 
#plot 
par(bg=NA) 
plot(X, pch=3, col='gray', lwd=2, axes=F, xlab="", ylab="", type='n', 
    ylim=c(min(X$y_polar)-150, max(X$y_polar)), xlim=c(min(X$x_polar)-150, max(X$x_polar))) #blank plot 
points(X$x_polar, X$y_polar, pch=3, col='gray', lwd=2, axes=F, xlab="", ylab="") #observations 
lines(e, col="red", lwd=3, lty=2) #plotting the ellipse 
points(0,0,col=2, lwd=2, cex=2) #center/origin 

den Ursprung der Ellipse zu bringen, in der Mitte ändern wir konnten wie folgt (sicherlich nicht die beste Art und Weise es zu tun)

#generate mirror coordinates 
X$x_polar_mirror<- -X$x_polar 
X$y_polar_mirror<- -X$y_polar 

mydata<-as.matrix(data.frame(c(X$x_polar, X$x_polar_mirror), c(X$y_polar, X$y_polar_mirror))) 
#fit the data 
efit <- fit.ellipse(mydata) 
e <- get.ellipse(efit) 
par(bg=NA) 
plot(mydata, pch=3, col='gray', lwd=2, axes=F, xlab="", ylab="", type='n', 
    ylim=c(min(X$y_polar)-150, max(X$y_polar)), xlim=c(min(X$x_polar)-150, max(X$x_polar))) 
points(X$x_polar, X$y_polar, pch=3, col='gray', lwd=2, axes=F, xlab="", ylab="") 
lines(e, col="red", lwd=3, lty=2) 
points(0,0,col=2, lwd=2, cex=2) #center 

Nun ... es macht den Job, aber keiner wäre glücklich mit all diesen imaginären Punkten, die in der Berechnung berücksichtigt werden.

Methode # 02

Dies ist ein weiterer indirekter Weg, um die Daten des Anbringens aber wieder die Mittel Ellipse ist der Ursprung nicht. Irgendeine Problemumgehung?

require(car) 
dataEllipse(X$x_polar, X$y_polar, levels=c(0.15, 0.7), 
      xlim=c(-150, 400), ylim=c(-200,300)) 

Meine Fragen: (a) ist eine robuste Alternative Art und Weise, diese Punkte mit der Ellipse Mittelpunkt im Ursprung des Anpassens (0,0)? (b) Gibt es ein Maß für die Güte der Ellipsenanpassung? Vielen Dank im Voraus.

Answer 1

ich mit aproach nicht wirklich glücklich bin ich concieved habe, sollte es eine geschlossene Lösung sein, aber nach wie vor:

# Ellipse equasion with center in (0, 0) with semiaxis pars[1] and pars[2] rotated by pars[3]. 
# t and pars[3] in radians 

ellipsePoints <- function(t, pars) { 
    data.frame(x = cos(pars[3]) * pars[1] * cos(t) - sin(pars[3]) * pars[2] * sin(t), 
      y = sin(pars[3]) * pars[1] * cos(t) + cos(pars[3]) * pars[2] * sin(t)) 
} 


# Way to fit an ellipse through minimising distance to data points. 
# If weighted then points which are most remote from center will have bigger impact. 

ellipseBrute <- function(x, y, pars, weighted = FALSE) { 

    d <- sqrt(x**2 + y**2) 
    t <- asin(y/d) 
    w <- (d/sum(d))**weighted 

    t[x == 0 & y == 0] <- 0 

    ep <- ellipsePoints(t, pars) 

    sum(w*(sqrt(ep$x**2 + ep$y**2) - d)**2) 
} 

# Fit through optim. 

opt_res <- optim(c(diff(range(X$x_polar)), 
        diff(range(X$y_polar)), 
        2*pi)/2, 
       ellipseBrute, 
       x = X$x_polar, y = X$y_polar, 
       weighted = TRUE 
       ) 
# Check resulting ellipse throuh plot 

df <- ellipsePoints(seq(0, 2*pi, length.out = 1e3), opt_res$par) 

plot(y ~ x, df, col = 'blue', t = 'l', 
    xlim = range(c(X$x_polar, df$x)), 
    ylim = range(c(X$y_polar, df$y))) 
points(0, 0, pch = 3, col = 'blue') 

points(y_polar ~ x_polar, X)