! For a given set of 3 or more points in 2 dimensions,
measure the circularity of the points by
finding two concentric circles,
a) a smallest enclosing circle, and
b) a largest enclosed circle.
The difference in squared radius of the two circles is a measure of circularity.
If the points are precisely on a circle,
then the difference in radius will be zero;
! Keywords: Chart, Chebychev, Circularity, Gass, Harary, LINGO, MinMax,
Quality control, Roundness, Scatter plot, Sphere, Witzgall;
! Ref: Gass, S., C. Witzgall, and H.H. Harary (1998),
"Fitting Circles and Spheres to Coordinate Measuring Machine Data,"
International J. of Flexible Manufacturing Systems, v. 10, pp 525.;
SETS:
POINT: X, Y, R; ! The data points;
POINTG: XG, YG; ! Some reference points to graph;
ENDSETS DATA:
!Case01;
X, Y =
0 0
4 1
2 3
0 1
;
!Case02
X, Y =
100 100
104 101
102 103
100 101
;
!Case03
X, Y =
100 100
96 101
98 103
100 101
;
!Case04
X, Y =
0.8660254 0.5
0 1
0.8660254 0.5
0.8660254 0.5
0 1
0.8660254 0.5
;
!Case05
X, Y =
7.7949 4.5
0 8.996
7.7953 4.4987
7.7963 4.4998
0 9.0031
7.7939 4.5001
;
!Case06
X, Y =
2 9
6 7
7 6
;
ENDDATA
! Given a set of points ( X(k), Y( k), we want to find a reference circle,
with center ( X0, Y0) and radius R0 such that the given points are close
to being on this reference circle.
The difference between
the squared radius of point k from (X0, Y0) and
the squared radius of the reference circle is:
( X(k)  X0)^2 + (Y( k)  Y0)^2  R0^2
= X( k)^2 + Y( k)^2  2* X( k)* X0  2* Y( k)* Y0 + X0^2 + Y0^2  R0^2
= X( k)^2 + Y( k)^2  2* X( k)* X0  2* Y( k)* Y0  RHO,
where
RHO = R0^2  X0^2  Y0^2;
! We want to find the X0, Y0, and R0 that minimize the maximum of the difference
between the squared nominal radius, the squared radius of enclosing circle,
and squared radius of the enclosed circle;
SUBMODEL FitCirc:
MIN = T;
@FOR( POINT( k):
R( k) = X( k)^2 + Y( k)^2  2* X( k)* X0  2* Y( k)* Y0  RHO;
[OUTR] T >= R( k); ! At least 1 of these must be binding;
[INNR] T >=  R( k); ! At least 1 of these must be binding, but not same k;
@FREE( R( k)); ! Unconstrained in sign;
);
@FREE( X0);
@FREE( Y0);
@FREE( RHO);
! Where we should have additionally the nonlinear constraint;
! [NONL] RHO = R0^2  X0^2  Y0^2;
! We want to show that if the nonlinear constraint [NONL] is dropped,
then at an optimal solution there is a feasible solution to [NONL],
so dropping it is justified.
The argument is as follows.
At an optimum, at least one of the [INNR] constraints must hold as an =, so
T =  X( k)^2 Y( k)^2 + 2* X( k)* X0 + 2* Y( k)* Y0 + RHO
Now T >= 0, so
0 <=  X( k)^2 Y( k)^2 + 2* X( k)* X0 + 2* Y( k)* Y0 + RHO
or
X( k)^2 + Y( k)^2  2* X( k)* X0  2* Y( k)* Y0 <= RHO
or
X( k)^2 + Y( k)^2  2* X( k)* X0  2* Y( k)* Y0 + X0^2 + Y0^2 <= RHO + X0^2 + Y0^2,
or
( X( k)  X0)^2 + ( Y( k)  Y0) ^2 <= RHO + X0^2 + Y0^2,
or
0 <= RHO + X0^2 + Y0^2,
thus we can solve:
R0^2 = RHO + X0^2 + Y0^2,
Because the RHS >= 0;
ENDSUBMODEL
CALC:
@SOLVE( FitCIrc);
R0 = ( RHO + X0^2 + Y0^2)^0.5; ! Get fitted radius;
! Display a graph;
! Generate a set of points on a circle;
NPTG = 64;
POINTG = 1..NPTG;
! Generate some points on the reference circle;
@FOR( POINTG( j):
XG( j) = X0 + R0* @COS( 2*@PI()*(j1)/NPTG);
YG( j) = Y0 + R0* @SIN( 2*@PI()*(j1)/NPTG);
);
! Generate a scatter plot of the two sets of points;
@CHARTSCATTER( 'Fit a Circle to a Set of Points, Tau= '+@FORMAT( T,"5.2f"), 'xaxis', 'yaxis',
'Data', X, Y,
'Fitted', XG, YG);
ENDCALC
