将数据点拟合到一个以其原点为中心的椭圆上使用R

Question

我有一个关于椭圆中心位于原点的数据拟合椭圆的问题。我已经探索了两种适合椭圆的方法，但是会产生一个任意的中心，除非我用一些假想的镜像点来处理数据。

方法＃01

脚本的这部分直接来自this有用的帖子。我直接在这里复制代码以方便。

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) 
    # [email protected] 
    # 
    # 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)) 
} 

让我们极性点的例子分布图的目的

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 

为了使椭圆的原点在中心，我们可以修改如下（肯定不是最好的的方式）

#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 

呃......这种工作可以做，但没有人会满意计算中考虑的所有虚拟点。

方法＃02

这是拟合数据的另一种间接的方式，但再次椭圆的中心位置不在原点。任何解决方法？ （a）是否有一个强大的替代方法来拟合椭圆中心在原点（0,0）的这些点？ （b）衡量椭圆拟合的好处吗？先谢谢你。

Answer 1

我不和的形式给出了我concieved真的很开心，应该有一个封闭的形式解决方案，但仍：

# 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)