Lord Bubbacub
Reputation: 145
least squares regression to fit circle with constrained centre point
I have a set of points describing a curve. I need to find the best fit circle for this curve whilst keeping the x coordinate of centre of the circle fixed. I.e the only variables for this circle will be the radius and the y coordinate of the centre of the circle. I've tried using some examples from the scipy cookbook which give me a very nice best fit circle (http://www.scipy.org/Cookbook/Least_Squares_Circle) but I'm not sure how I need to change this to set the x coordinate of the centre of the circle to a fixed value?
def fit_circ(xy, method):
xy = array(xy)
u = array(xy[:,0])
v = array(xy[:,1])
if method == 1:
#method_1 = 'algebraic'
# linear system defining the center in reduced coordinates (uc, vc):
# Suu * uc + Suv * vc = (Suuu + Suvv)/2
# Suv * uc + Svv * vc = (Suuv + Svvv)/2
# coordinates of the barycenter
x_m = mean(u)
y_m = mean(v)
# calculation of the reduced coordinates
u = u - x_m
v = v - y_m
Suv = sum(u*v)
Suu = sum(u**2)
Svv = sum(v**2)
Suuv = sum(u**2 * v)
Suvv = sum(u * v**2)
Suuu = sum(u**3)
Svvv = sum(v**3)
# Solving the linear system
A = array([ [ Suu, Suv ], [Suv, Svv]])
B = array([ Suuu + Suvv, Svvv + Suuv ])/2.0
xc_1, yc_1 = linalg.solve(A, B)
# Calculation of all distances from the center (xc_1, yc_1)
Ri_1 = sqrt((u-xc_1)**2 + (v-yc_1)**2)
R_1 = mean(Ri_1)
residu_1 = sum((Ri_1-R_1)**2)
residu2_1 = sum((Ri_1**2-R_1**2)**2)
xc_1 = x_m + xc_1
yc_1 = y_m + yc_1
return xc_1, yc_1, R_1, 0, residu_1
I have given some example data below. Very grateful for any help!
array([[ 0, 128],
[ 1, 128],
[ 2, 128],
[ 3, 128],
[ 4, 127],
[ 5, 127],
[ 6, 127],
[ 7, 127],
[ 8, 126],
[ 9, 126],
[ 10, 126],
[ 11, 125],
[ 12, 125],
[ 13, 125],
[ 14, 124],
[ 15, 124],
[ 16, 124],
[ 17, 124],
[ 18, 124],
[ 19, 124],
[ 20, 124],
[ 21, 124],
[ 22, 124],
[ 23, 124],
[ 24, 125],
[ 25, 125],
[ 25, 126],
[ 26, 126],
[ 26, 127],
[ 27, 127],
[ 28, 128],
[ 29, 129],
[ 29, 130],
[ 30, 130],
[ 30, 131],
[ 31, 131],
[ 32, 132],
[ 33, 133],
[ 34, 134],
[ 35, 134],
[ 35, 135],
[ 36, 135],
[ 37, 135],
[ 38, 136],
[ 46, 136],
[ 47, 135],
[ 48, 135],
[ 49, 135],
[ 50, 135],
[ 51, 135],
[ 52, 134],
[ 53, 134],
[ 54, 134],
[ 55, 133],
[ 56, 133],
[ 57, 133],
[ 58, 132],
[ 59, 132],
[ 60, 131],
[ 61, 131],
[ 62, 130],
[ 63, 130],
[ 64, 129],
[ 65, 129],
[ 66, 129],
[ 67, 128],
[ 68, 128],
[ 69, 128],
[ 70, 127],
[ 71, 127],
[ 72, 127],
[ 73, 126],
[ 74, 126],
[ 75, 126],
[ 76, 126],
[ 77, 125],
[ 78, 125],
[ 79, 125],
[ 80, 124],
[ 81, 124],
[ 82, 124],
[ 83, 123],
[ 84, 123],
[ 85, 123],
[ 86, 122],
[ 87, 122],
[ 88, 121],
[ 89, 121],
[ 90, 120],
[ 91, 120],
[ 92, 119],
[ 93, 119],
[ 94, 118],
[ 95, 118],
[ 96, 117],
[ 97, 117],
[ 98, 116],
[ 99, 115],
[ 99, 116],
[100, 115],
[101, 114],
[102, 114],
[103, 113],
[104, 113],
[105, 112],
[106, 111],
[106, 112],
[107, 111],
[108, 110],
[109, 110],
[110, 109],
[111, 109],
[112, 108],
[113, 108],
[114, 107],
[115, 106],
[115, 107],
[116, 106],
[117, 105],
[118, 105],
[119, 104],
[120, 104],
[121, 103],
[122, 103],
[123, 102],
[124, 102],
[125, 102],
[126, 101],
[127, 101],
[128, 100],
[129, 100],
[130, 99],
[131, 99],
[132, 98],
[133, 98],
[134, 97],
[135, 97],
[136, 96],
[137, 96],
[138, 95],
[139, 95],
[140, 95],
[141, 94],
[142, 94],
[143, 93],
[144, 93],
[145, 92],
[146, 92],
[147, 91],
[148, 91],
[149, 90],
[150, 90],
[151, 89],
[152, 89],
[153, 88],
[154, 88],
[155, 87],
[156, 87],
[157, 86],
[158, 85],
[158, 86],
[159, 85],
[160, 84],
[161, 83],
[161, 84],
[162, 83],
[163, 82],
[164, 81],
[164, 82],
[165, 80],
[165, 81],
[166, 80],
[167, 79],
[168, 78],
[168, 79],
[169, 77],
[169, 78],
[170, 77],
[171, 76],
[172, 75],
[172, 76],
[173, 74],
[173, 75],
[174, 73],
[174, 74],
[175, 72],
[175, 73],
[176, 70],
[176, 71],
[177, 69],
[177, 70],
[178, 68],
[178, 69],
[179, 67],
[179, 68],
[180, 66],
[180, 67],
[181, 64],
[181, 65],
[182, 56],
[182, 57],
[182, 58],
[182, 59],
[182, 60],
[182, 61],
[182, 62],
[182, 63],
[183, 54],
[183, 55],
[184, 53],
[184, 54],
[185, 52],
[185, 53],
[186, 52],
[187, 51],
[188, 51],
[189, 50],
[190, 50],
[191, 49],
[192, 49],
[193, 49],
[194, 48],
[195, 47],
[195, 48],
[196, 47],
[197, 46],
[198, 45],
[198, 46],
[199, 44],
[199, 45],
[200, 43],
[200, 44],
[201, 42],
[201, 43],
[202, 41],
[202, 42],
[203, 39],
[203, 40],
[204, 37],
[204, 38],
[205, 36],
[205, 37],
[206, 34],
[206, 35],
[207, 32],
[207, 33],
[208, 30],
[208, 31],
[209, 29],
[209, 30],
[210, 27],
[210, 28],
[211, 26],
[211, 27],
[212, 24],
[212, 25],
[213, 23],
[213, 24],
[214, 22],
[214, 23],
[215, 21],
[215, 22],
[216, 19],
[216, 20],
[217, 18],
[217, 19],
[218, 17],
[218, 18],
[219, 15],
[219, 16],
[220, 14],
[220, 15],
[221, 12],
[221, 13],
[222, 10],
[222, 11],
[223, 7],
[223, 8],
[223, 9],
[224, 6]])
Upvotes: 0
Views: 1092
Answers (1)
Pascal Getreuer
Reputation: 3256
We can keep the x coordinate of the circle fixed by replacing the linear system solve
xc_1, yc_1 = linalg.solve(A, B)
with
xc_1 = x_fixed - x_m
yc_1 = A[:, 1].dot(B - A[:, 0] * xc_1) / A[:, 1].dot(A[:, 1])
in which "x_fixed" is the circle's x coordinate that we specify, and all else is the same as before.
Explanation: rather than solving for xc_1, we're setting it according to x_fixed and finding the least squares fit for yc_1 on the reduced system. For further details on what this system means, see the derivation in "Least-Squares Circle Fit" by Randy Bullock.
With your test data, here is the resulting fit for a couple different values for x_fixed:
xc_1 yc_1 R_1 residu_1
-------------------------------
10.0 -116.02 247.66 2710.54
35.0 -71.07 199.45 2866.73
Upvotes: 0
Related Questions
- How to find the center of circle using the least square fit in python?
- Least square method in python?
- Python: Circular regression to determine center and radius
- python optimize.leastsq: fitting a circle to 3d set of points
- How to fit a circle with a given radius to sample data points
- How to fit a circle to a set of points with a constrained radius?
- Best fit circle with fix center in MATLAB
- Least squares for circle detection
- Fitting a 3d points of an arc to a circle (regression in Python)
- least squares curve fitting