diff --git a/bin/desi_fit_circles b/bin/desi_fit_circles index e11f0d01..8081f3ac 100755 --- a/bin/desi_fit_circles +++ b/bin/desi_fit_circles @@ -3,100 +3,12 @@ import sys from astropy.table import Table import numpy as np -from scipy import optimize import matplotlib.pyplot as plt import argparse - +from desimeter.circles import fit_circle from desimeter.transform.xy2qs import xy2uv, uv2xy - -#- Least squares fit to circle; adapted from "Method 2b" in -#- https://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html - -def fit_circle(xin, yin): - - #- Transform to curved focal surface which is closer to a real circle - x, y = xy2uv(xin, yin) - x_m, y_m, r = fast_fit_circle(x, y) - - #- If r is too small or too big then either this positioner wasn't moving - #- or the points are mismatched for a bad fit. - if (r < 1.0) or (r > 5.0): - raise ValueError('Bad circle fit') - - def calc_R(xc, yc): - """ calculate the distance of each data points from the center (xc, yc) """ - return np.sqrt((x-xc)**2 + (y-yc)**2) - - def f_2b(c): - """ calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """ - Ri = calc_R(*c) - return Ri - Ri.mean() - - def Df_2b(c): - """ Jacobian of f_2b - The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq""" - xc, yc = c - df2b_dc = np.empty((len(c), x.size)) - - Ri = calc_R(xc, yc) - df2b_dc[0] = (xc - x)/Ri # dR/dxc - df2b_dc[1] = (yc - y)/Ri # dR/dyc - df2b_dc = df2b_dc - df2b_dc.mean(axis=1)[:, np.newaxis] - - return df2b_dc - - - center_estimate = x_m, y_m - center_2b, ier = optimize.leastsq(f_2b, center_estimate, Dfun=Df_2b, col_deriv=True) - - xc_2b, yc_2b = center_2b - Ri_2b = calc_R(*center_2b) - R_2b = Ri_2b.mean() - residu_2b = sum((Ri_2b - R_2b)**2) - - if (R_2b < 1.0) or (R_2b > 5.0): - raise ValueError('Bad circle fit') - - #- Convert center back into CS5 x,y - xc, yc = uv2xy(xc_2b, yc_2b) - # xc, yc = xc_2b, yc_2b - - return xc, yc, R_2b - -def fast_fit_circle(x,y) : - # init - nn=len(x) - i1=np.arange(nn) - ### i2=(i1+1)%nn - i2=(i1+nn//2-1)%nn - - # midpoints - mx=((x[i1]+x[i2])/2.) - my=((y[i1]+y[i2])/2.) - nx=(y[i2]-y[i1]) - ny=-(x[i2]-x[i1]) - - # solve for intersection of perpendicular bisectors - # with s1,s2 are affine parameters of 2 adjacent perpendicular bisectors - # 2 equations: - # mx1 + nx1*s1 = mx2 + nx2*s2 - # my1 + ny1*s1 = my2 + ny2*s2 - s1 = (ny[i2]*mx[i2]-nx[i2]*my[i2]-ny[i2]*mx[i1]+nx[i2]*my[i1])/(ny[i2]*nx[i1]-nx[i2]*ny[i1]) - - # coordinates of intersections are estimates of center of circle - xc=mx[i1]+nx[i1]*s1[i1] - yc=my[i1]+ny[i1]*s1[i1] - - # first estimate of center is mean of all intersections - xc=np.mean(xc) - yc=np.mean(yc) - r=np.mean(np.sqrt((x-xc)**2+(y-yc)**2)) - - return xc,yc,r - - parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter, description="""FVC image processing""") parser.add_argument('-i','--infile', type = str, default = None, required = True, nargs="*", @@ -188,7 +100,12 @@ for iloc,loc in enumerate(x.keys()) : if pinhole[iloc] == 0 : # it's a positioner # here is the fit try: - xc,yc,r = fit_circle(x[loc],y[loc]) + #- Transform to curved focal surface which is closer to a real circle + x_cfs, y_cfs = xy2uv(x[loc], y[loc]) + #- Do the fit + xc_cfs,yc_cfs,r = fit_circle(x_cfs, y_cfs) + #- Convert center back into CS5 x,y + xc, yc = uv2xy(xc_cfs, yc_cfs) except ValueError: print("fit circle failed for loc={} x={} y={}".format(loc,xc,yc)) continue diff --git a/py/desimeter/circles.py b/py/desimeter/circles.py new file mode 100644 index 00000000..5bc28af5 --- /dev/null +++ b/py/desimeter/circles.py @@ -0,0 +1,109 @@ +""" +Utility functions to fit circles from a set of cartesian coordinates +""" + +import numpy as np +from scipy import optimize + +#- Least squares fit to circle; adapted from "Method 2b" in +#- https://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html + +def fit_circle(x, y): + """ + fit a circle from a set of 2D cartesian coordinates + Args: + x : float numpy array of coordinates along first axis of cartesian coordinate system + y : float numpy array of coordinates along second axis in same system + + returns: + xc : float, coordinates along first axis of center of circle + yc : float, coordinates along second axis of center of circle + r : float, radius of circle + """ + + x_m, y_m, r = _fast_fit_circle(x, y) + + #- If r is too small or too big then either this positioner wasn't moving + #- or the points are mismatched for a bad fit. + if (r < 1.0) or (r > 5.0): + raise ValueError('Bad circle fit') + + def calc_R(xc, yc): + """ calculate the distance of each data points from the center (xc, yc) """ + return np.sqrt((x-xc)**2 + (y-yc)**2) + + def f_2b(c): + """ calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """ + Ri = calc_R(*c) + return Ri - Ri.mean() + + def Df_2b(c): + """ Jacobian of f_2b + The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq""" + xc, yc = c + df2b_dc = np.empty((len(c), x.size)) + + Ri = calc_R(xc, yc) + df2b_dc[0] = (xc - x)/Ri # dR/dxc + df2b_dc[1] = (yc - y)/Ri # dR/dyc + df2b_dc = df2b_dc - df2b_dc.mean(axis=1)[:, np.newaxis] + + return df2b_dc + + + center_estimate = x_m, y_m + center_2b, ier = optimize.leastsq(f_2b, center_estimate, Dfun=Df_2b, col_deriv=True) + + xc_2b, yc_2b = center_2b + Ri_2b = calc_R(*center_2b) + R_2b = Ri_2b.mean() + residu_2b = sum((Ri_2b - R_2b)**2) + + if (R_2b < 1.0) or (R_2b > 5.0): + raise ValueError('Bad circle fit') + + return xc_2b, yc_2b, R_2b + +def _fast_fit_circle(x,y) : + """ + fast and approximate method to fit a circle from a set of 2D cartesian coordinates + (better use fit_circle) + Args: + xin : float numpy array of coordinates along first axis of cartesian coordinate system + yin : float numpy array of coordinates along second axis in same system + + returns: + xc : float, coordinates along first axis of center of circle + yc : float, coordinates along second axis of center of circle + r : float, radius of circle + """ + # init + nn=len(x) + i1=np.arange(nn) + ### i2=(i1+1)%nn + i2=(i1+nn//2-1)%nn + + # midpoints + mx=((x[i1]+x[i2])/2.) + my=((y[i1]+y[i2])/2.) + nx=(y[i2]-y[i1]) + ny=-(x[i2]-x[i1]) + + # solve for intersection of perpendicular bisectors + # with s1,s2 are affine parameters of 2 adjacent perpendicular bisectors + # 2 equations: + # mx1 + nx1*s1 = mx2 + nx2*s2 + # my1 + ny1*s1 = my2 + ny2*s2 + s1 = (ny[i2]*mx[i2]-nx[i2]*my[i2]-ny[i2]*mx[i1]+nx[i2]*my[i1])/(ny[i2]*nx[i1]-nx[i2]*ny[i1]) + + # coordinates of intersections are estimates of center of circle + xc=mx[i1]+nx[i1]*s1[i1] + yc=my[i1]+ny[i1]*s1[i1] + + # first estimate of center is mean of all intersections + xc=np.mean(xc) + yc=np.mean(yc) + r=np.mean(np.sqrt((x-xc)**2+(y-yc)**2)) + + return xc,yc,r +