# Fitting an arc from samples taken at evenly spaced angles

I have an application where I am rotating an object and measuring its position in X,Y as it rotates. There is some measurement error in both X and Y at each point (haven't yet been able to characterize, but I'm hoping roughly normally distributed.) I want to fit these points to a circle so that I can determine the center of rotation.

My first attempt was to use a least squares fit based on the data moments, like that described in this paper. The problem with this approach is that it tends to greatly underestimate the radius in my qualification experiments.

I've come up with a 'hackier' approach that seems to be working okay, by taking into account the fixed angular spacing. Here is what I'm currently doing:

1. Let's index the samples, spaced at angles $\Delta w$, so that we have $i \in [-n,n]$, with coordinates $(x_i,y_i)$.
2. First transform the coordinates, setting $(u_i,v_i) = (x_i - x_0, y_i - y_0)$.
3. Take the presumed chord formed by each $(u_i,v_i)$ and compute the radius $r$ using $\theta_i = i \Delta w$ and $r_i = \frac{\sqrt {u_i^2 + v_i^2}} {2\sin(\theta_i/2)}$.
4. Take the mean of the radius estimates and use this as our radius.
5. For each $k \in [1,n]$, compute the vector perpendicular to $(u_k,v_k)-(u_{-k},v_{-k})$. Take the average of all these vectors to yield a vector $p$.
6. Check if $p$ is facing the 'wrong direction' by summing the dot product of $p$ with the various $(u_i,v_i)$. If so, negate it.
7. The center of the circle is then $c = (x_0,y_0) + r\frac{p}{|p|}$.

The above doesn't strike me as very robust, however, as I don't believe it's doing any sort of minimization in a least squares sense. Is there a better approach to this? I also plan to repeat this process by sampling multiple times at each angle, but is there a better way to incorporate that data into a single fit?

In response to comment, some additional details:

1. Measurements are taken by rotating a mechanism to known angles, then imaging a feature and recording the feature's position at each angle.
2. There is error both in the physical positioning of the mechanism and the vision-based measurement of position.
3. I'm typically taking 21 samples ranging +/- 10 degrees.

Here is a plot of an example data set vs the circle fit from my hacky algorithm. Note that the fit here is particularly poor due to the middle $(x_0,y_0)$ point being offset to the left from the bulk of the data, since I currently always include that center point in the arc. This still works better than fitting with the moments, though, as it imposes that the points cover a roughly 20 degree portion of the circle.

Here is the raw data for the above case, spaced at 1 degree intervals:

ErrXCal ErrYCal
-0.04   0.1964
-0.0309 0.1764
-0.0099 0.1484
-0.0226 0.135
-0.0233 0.1332
-0.0261 0.1125
0.002   0.0633
0.0048  0.0404
-0.0011 0.0474
-0.0149 0.0353
-0.0073 -0.0017
0.0253  -0.019
0.0317  -0.0369
0.065   -0.069
0.0652  -0.1065
0.098   -0.1035
0.0837  -0.1064
0.1072  -0.1292
0.1179  -0.1717
0.1839  -0.175
0.1691  -0.178

• If the method in that paper fails you, then the measurement error must be very far from any kind of Normal (or even symmetric) distribution. This leads me to suspect it will be important to explain your data in more detail and perhaps provide a small dataset that exemplifies the problems you have having. The title suggests there is something special about your measurements, too: exactly how do you obtain "equally spaced angles"? Angles of what? – whuber Jul 15 '15 at 14:59
• @whuber, added some additional detail; hopefully not too much – Dan Bryant Jul 15 '15 at 15:23

These data are interesting in that they hide a subtle interdependence of errors. It can be seen only after fitting the model correctly. Choosing the right model makes a large difference in the estimates.

Let each observation be a point $Z_i$ in the plane, written with a capital letter to remind us it will be modeled as a random variable. Specifically, let $\zeta$ be the circle's center, $\rho$ its radius, and let $\theta$ be a reference angle. Then

$$Z_i = \zeta + \rho(\cos(\delta i + \theta), \sin(\delta i+\theta)) + E_i$$

describes points $Z_i$ on the circle that have been displaced randomly by vectors $E_i$ ($i=1,2,\ldots, 21$ for these data). The angular increment $\delta$ is one degree for these data. The data themselves are considered realizations of the $Z_i$, written $z_i$.

A least-squares solution is appropriate when the $E_i$ are independent and each has a circularly symmetric distribution around $(0,0)$ (generalizing the assumption of approximate normality). It finds values of $\zeta = (\xi,\eta)$, $\rho$, and $\theta$ minimizing the sum of squared residuals

$$f(\zeta,\rho,\theta) = \sum_i |z_i - (\zeta + \rho(\cos(\delta i + \theta), \sin(\delta i+\theta)))|^2.$$

Although the gradient of $f$ is a nonlinear function of $\theta$, $f$ is readily minimized provided we can supply a reasonable starting estimate for the parameters. One such estimate can be obtained by selecting three widespread points that typify the data locations. I have taken the centroids of three groups of five points: the first five, the last five, and the middle five. The starting parameters describe the circumcircle of these three points. Let's call this the "geometric fit."

The figure displays the geometric fit and the subsequent least squares fit obtained by starting with the geometric fit and applying a nonlinear solver to improve the value of $f$.

Visually, the geometric fit is excellent: the circular arc seems to split the data points (shown as hollow dots) pretty evenly. The three centroids used to find this fit are shown as solid dots.

The least squares fit does not look quite as good. To understand why not, I have drawn it a little differently. Recall that the model prescribes a location for each point along known angles. Thus, to each of those (21) known angles corresponds a definite predicted location. These are shown as small solid red dots. The correspondence between each observation $z_i$ and its fit $\hat z_i$ is indicated by drawing a light arrow from one towards the other.

For the first time we begin to see the role played by the sequence of points. They aren't just some geometric jumble, as shown in the left hand plot, but really are a two-dimensional time series with an associated structure. And an interesting structure it is! Look at how the gray arrows seem to spin around as they track along the fitted curve. The residuals are not independent. In other words, the locational error is changing systematically over time, and not entirely randomly as assumed.

This is evident in a plot of the residuals:

The x-coordinates are shown with circles (and a Loess smooth in black) while the y-coordinates are shown with triangles (and a Loess smooth in red). There is a systematic variation in the x-coordinates: they start high, drop low, then go back high. (Geometrically, the points start to the right of the fitted circle, move to its left, and then back to its right. They wobble up and down, too, but those wobbles look random.)

Because of this, the least-squares fit is not entirely trustworthy. However, because the x-residuals have at least cycled through one phase of high-low-high, the least-squares fit is probably pretty good compared to the geometric fit (which pays no attention to the underlying angle associated with each point and is confused by the systematic variation in the x-direction). Here are the two sets of parameter estimates:

                zeta(1)   zeta(2)       rho     theta
Geometric     0.3683045 0.1630800 0.3936789 -3.145837
Least squares 1.1618824 0.5750505 1.2705868 -2.860703


The least-squares radius of 1.27 is three times the geometric radius of 0.40. That large difference attests to the importance of using the right model to make the fit.

For an improved fit, look into multidimensional time-series methods. A simple one would be to adopt a distributional assumption and correlation structure for the residuals and perform the fitting using maximum likelihood.

This R code did the work.

#
# The point sequence.
# All arrays of points will be maintained as 2 x n matrices, one point
# per column.
#
z <- matrix(scan(text="-0.04   0.1964
-0.0309 0.1764
-0.0099 0.1484
-0.0226 0.135
-0.0233 0.1332
-0.0261 0.1125
0.002   0.0633
0.0048  0.0404
-0.0011 0.0474
-0.0149 0.0353
-0.0073 -0.0017
0.0253  -0.019
0.0317  -0.0369
0.065   -0.069
0.0652  -0.1065
0.098   -0.1035
0.0837  -0.1064
0.1072  -0.1292
0.1179  -0.1717
0.1839  -0.175
0.1691  -0.178"), nrow=2, dimnames=list(c("x", "y")))
#
# The nominal angles, in increments of one degree.
#
theta <- 1/180 * pi * (1:dim(z)[2])
#
# The objective function.
#
f <- function(parms) {
zeta <- parms[1:2]   # Center
theta.0 <- parms[4]  # Start angle, radians

theta.d <- theta + theta.0
w <- z - (zeta + rho*rbind(cos(theta.d), sin(theta.d)))
sum(w^2)
}
#
# Rough estimate of starting values: use a triangle to approximate the data.
#
n <- dim(z)[2]                  # Number of data points
k <- 2                          # Will use clusters of 2*k+1 points
pts <- sapply(c(k+1, floor(n/2), n-k),
function(i) c(mean(z[1, (i-k):(i+k)]), mean(z[2, (i-k):(i+k)])))
#
# Find the center and radius of the circumcircle.
#
v <- apply(pts, 1, diff)        # Direction vectors for the triangle's sides
x <- solve(cbind(c(-v[1,2], v[1,1]), c(v[2,2], -v[2,1])), apply(v, 2, mean))
zeta <- v[1, ]/2 - x[1] * c(v[1,2], -v[1,1]) + pts[, 1]
rho <- sqrt(sum((p <- pts[, 1]-zeta)^2))
theta.0 <- atan2(p[2], p[1]) - theta[1] # Guess at the angle of origin
#
# Plot this fit.
#
par(mfrow=c(1,2))
a <- seq(0, 360, 1) / 180 * pi            # Angular coordinates on a circle
plot(t(z), asp=1, main="Geometric fit")   # The points
lines(t(rho*rbind(cos(a), sin(a)) + zeta))# The fitted circle (clipped)
points(t(pts), pch=16)                    # The triangle's vertices
points(matrix(zeta, 1, 2), pch=16)        # The circle's center (clipped)
#
# Refine this estimate.
#
fit <- nlm(f, c(zeta, rho, theta.0))      # Newton-Raphson improvement
rbind(c(zeta, rho, theta.0), fit$estimate)# Compare the two fits prediction <- fit$estimate[3]*rbind(cos(theta+fit$estimate[4]), sin(theta+fit$estimate[4])) + fit\$estimate[1:2]
#
# Plot it.
#
plot(rbind(t(z), t(prediction)), type="n", asp=1,
main="Least squares fit")
invisible(sapply(1:n, function(i)
arrows(z[1,i], z[2,i], prediction[1,i], prediction[2,i],
length=0.1, angle=20, col="#a0a0a0")))
lines(t(prediction), col="Red", lty=3)
points(t(z))
points(t(prediction), col="Red", pch=16, cex=0.7)
#
# Plot the residuals.
#
par(mfrow=c(1,1))
residuals <- z - prediction
k <- max(abs(residuals[2,]))
plot(c(0,n+1)*k, range(residuals[2,]), type="n",
asp = diff(range(residuals[2,])) / (k*1.1),
xaxt="n", bty="n",
xlab="", ylab="", main="Residual Vectors")
abline(h=0, col="Gray")
invisible(sapply(1:n, function(i)
arrows(i*k, 0, i*k+residuals[1,i], residuals[2,i],  length=0.1, angle=20)))
#
# Plot the residual components.
#
plot(residuals[1,], type="n", ylab="", main="Residuals")
lines(lowess(residuals[2,], f=0.4), col="Red", lwd=2)
lines(lowess(residuals[1,], f=0.4), lwd=2)
points(residuals[1,])
points(residuals[2,], col="Red", pch=2)

• Very interesting; the data is probably revealing subtle errors inside the kinematics of the mechanism that wind up being expressed in systematic (yet non-linear) ways. I'll have to try some experiments with revisiting the same points repeatedly to see how consistent the error is over the series. – Dan Bryant Jul 15 '15 at 19:52
• On a side note, did you generate that 'vector-error' plot in Matlab or Octave? That's a really handy tool for 2d data like this. – Dan Bryant Jul 15 '15 at 19:52
• It's in R. I'll post the code if you want it. – whuber Jul 15 '15 at 20:02
• Ah, very cool; R looks like a very interesting tool. If it's not too much trouble, please do post the code; I may be able to use it to delve further into some of the details of what's happening on the system. – Dan Bryant Jul 15 '15 at 20:05