How to calculate best center of an imperfect circle from N measurements of its radius

It would appear that the results you're quoting for the circularity only hold under some approximations.

One of the definitions of circularity is via fitting the data to a circle, by optimizing the following error function, which measures deviations from a perfect circle, with respect to it's parameters:

$$E=\sum_{i=1}^{N}(\sqrt{(x_i-a)^2+(y_i-b)^2}-R)^2$$

Indeed, this loss function is zero if all observations fall exactly on a circle of some center $a,b$ and some radius $R$. Minimizing this loss function with respect to $a,b,R$ will give us estimators for the best values of the position of the center and the radius of the circle.

This problem however is non-linear, and thus analytical expressions for the estimators cannot be derived. If we make the assumption that the center of the circle doesn't differ much from the origin $(0,0)$ we can successively approximate

$$\sqrt{(x_i-a)^2+(y_i-b)^2}\approx\sqrt{x_i^2+y_i^2-2(ax_i+by_i)}=r_i\sqrt{1-2\frac{a\cos\theta_i+b\sin\theta_i}{r_i}}\\\approx r_i-a\cos\theta_i-b\sin\theta_i$$

where $r_i^2=x_i^2+y_i^2$.

Applying that to the loss function we recover a linear least squares estimation problem:

$$E=\sum_{i=1}^N(r_i-R-a\cos\theta_i-b\sin\theta_i)^2$$

which we can minimize and obtain the equations

$$\frac{\partial E}{\partial a}=2\sum_i\cos\theta_i(r_i-R-a\cos\theta_i-b\sin\theta_i)=0\\\frac{\partial E}{\partial b}=2\sum_i\sin\theta_i(r_i-R-a\cos\theta_i-b\sin\theta_i)=0\\\frac{\partial E}{\partial R}=2\sum_i(r_i-R-a\cos\theta_i-b\sin\theta_i)=0\\$$

which constitutes a system of linear equations that can be written in the form:

$$\begin{pmatrix}\sum\cos^2\theta_i&\sum\cos\theta_i\sin\theta_i&\sum\cos\theta_i\\\sum\cos\theta_i\sin\theta_i&\sum\sin^2\theta_i&\sum\sin\theta_i\\\sum\cos\theta_i&\sum\sin\theta_i&\sum 1\end{pmatrix}\begin{pmatrix}a\\b\\R\end{pmatrix}=\begin{pmatrix}\sum r_i\cos\theta_i\\\sum r_i\sin\theta_{i}\\\sum r_i\end{pmatrix}$$

This is easy enough to solve explicitly but, if we would like a really simple formula, then we can also assume that the angles at which the measurements are taken are equally spaced (which is a reasonable assumption to make for engineering purposes). In this case we can show that for $\theta_i=\theta_0+\frac{2\pi (i-1)}{N}$ the following identities hold:

$$\sum_i\cos\theta_i=\sum_i\sin\theta_i=\sum_i \cos\theta_i\sin\theta_i=0\\\sum_i{\cos^2\theta_i}=\sum\sin^2\theta_i=\frac{N}{2}$$

and the linear system becomes diagonal and thus

$$a=\frac{2}{N}\sum_ir_i\cos\theta_i~,~b=\frac{2}{N}\sum_ir_i\sin\theta_i~,~r=\frac{1}{N}\sum_i r_i$$

Given this result, there are a few points to be addressed here:

1)This explains why in this confusing thread of simulations you've been finding constantly contradicting results. The factor of $2$ appearing in front of the weighted average for the center of mass is not describing an exact solution, but a solution to an approximated version of the minimization problem. As one moves the true center away from (0,0) more and more, one should find that the textbook estimators provided above are less and less accurate.

2)It is easy to see that the estimator derived above does NOT minimize the minimum circumscribed-maximum inscribed radius recipe, but the loss function shown above. We shouldn't be mixing the two, and each time we estimate roundness we should be referring to the method used to estimate it. (see for example here and here) In many of the above examples, it is not clear which method is used. Not all methods give the same answers, but they do give similar answers for many regular examples.

3)The reason for the constant being $K=2$ and not $K=1$ for example, is explained above. One could obtain a constant of $K=1$ , for example if the following loss function is considered, with angularly equidistributed measurements:

$$E^*=\sum_i(a+(R-R_i)\cos\theta_i)^2+(b+(R-R_i)\sin\theta_i)^2$$

Personally, I don't think that there is any practical value (there always is mathematical value in looking into a non-linear problem, however) in extracting a best fitting value of $K$ for centers that are significantly far away from the origin, since it does not accurately represent corrections to the approximate solutions quoted and it varies a lot with the individual set of observations obtained.