Data fit to a circle using python

Question

I am trying to fit a set of data to an off-center-circle using python. However, the displacement of the center of the circle (to the center of the axis) is evaluated negative, when it should be positive.

Consider a circle of radius a with center at (θ_0, r_0). Therefore,

From this, I developed the following equations (I suppose everything is right):

Point (θ, r) on the circle:

Partial derivative of r with respect to a:

Partial derivative of r with respect to r_0:

Partial derivative of r with respect to θ_0:

This is a MWE:

import numpy as np               #... for mathematical operations and constants
from scipy import optimize       #... for data fitting
import matplotlib.pyplot as plt  #... for plots


dsts = [1.0238822745042868, 1.0179614005977726, 1.0288142236514832, 1.014240255062302, 1.0297548330833965]
angs = [1.2118767129907884, 0.46576418683129606, -0.2945355605442553, 0.7614519513889376, -0.7442328029712463]


#--> DEFINE FUNCTIONS TO FIT THE DATA
#... see: https://scipython.com/book/chapter-8-scipy/examples/non-linear-fitting-to-an-ellipse/

def circle(theta, p): #... defines the function for circle
    a, r0, theta_0 = p
    rc = r0*np.cos(theta - theta_0)
    return ( rc + (rc**2 + a**2 - r0**2)**0.5 )

def residuals(p, r, theta):
    """ Return the observed - calculated residuals using f(theta, p). """
    return r - circle(theta, p)

def jac(p, r, theta):
    """ Calculate and return the Jacobian of residuals. """
    tmp = circle(theta, p)
    a, r0, theta_0 = p
    da  = a/( tmp - r0*np.cos(theta - theta_0) )
    dr  = ( r*np.cos(theta - theta_0) - r0 )/( r - r0*np.cos(theta - theta_0) )
    dt  = tmp*r0*np.sin(theta - theta_0)/( tmp - r0*np.cos(theta - theta_0) )
    
    return -da,  -dr, -dt
    return np.array((-da, -dr, -dt)).T


#--> FIT THE DATA
p0 = (1, 0.05, 0) #... initial guess
plsq = optimize.leastsq(residuals, p0, Dfun=jac, args=(dsts, angs), col_deriv=True)


#--> CREATE FITTING CURVE ARRAYS
t_fit = np.linspace(0, 2*np.pi, 100)
r_fit = circle(t_fit, plsq[0])


#--> PLOT
cm = 1/2.54  #... centimeters-to-inches conversor
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, figsize=(8.5*cm, 8.5*cm), dpi=300)

#--> Write the title
fit_a = str(round(plsq[0][0], 4))
fit_r = str(round(plsq[0][1], 4))
fit_t = str(round(plsq[0][2]*180/np.pi, 1))
pltTitle = "$a =$" + fit_a + "; $r_0 =$ " + fit_r + "; $\\theta_0 =$ " + fit_t + "$^\circ$"
ax.set_title(pltTitle, va='bottom', fontsize=8)

#--> Plot the measured data
ax.scatter(angs, dsts)

#--> Plot the fitted orbit
ax.plot(t_fit, r_fit, linestyle='dotted')

#--> Configure the axis
ax.set_rmax(1.2)
ax.set_rticks(np.linspace(0.2, 1, 5))
ax.set_xticklabels([])      #... remove angular tick labels
ax.set_yticklabels([])      #... remove radial tick labels
ax.grid(True)

#--> Display the plot
plt.show()

Any suggestion to get positive r_0?

Answer 1

After examining your problem a bit closer I realized that there is in fact nothing wrong with your results!

Your problem just has an ambiguity with respect to r_0 and theta_0.

r( r0, theta_0 )  = r( -r0, theta_0 + [2n+1]pi )

Proper boundary conditions (or close enough start-values) should do the job just fine!

Alternatively you can also just "sanitize" the obtained results by applying something like:

def sanitize_results(a, r0, theta0):
    if r0 < 0:
        return a, -r0, (theta0 + np.pi)%(2*np.pi)
    else:
        return a, r0, theta0%(2*np.pi)

---

UPDATE

To use explicit boundary-conditions in the least-squares solver, you have to switch to the newer scipy interface, e.g.: scipy.optimize.least_squares

plsq = least_squares(
    residuals, 
    p0, 
    jac=jac, 
    args=(dsts, angs),
    bounds=((0, 0, 0), (10, 10, 2*np.pi))
    )

the obtained parameter values are then accessible via

plsq.x

... and a final warning: Your sketch is actually a bit misleading... your equation will only work if the center of the coordinate-system is inside the circle! ...otherwise your equation is not valid for all 360° (e.g. there is more than 1 possible solution for r given a set of parameters (a, r0, thata0)