Merge pull request #53 from desihub/fit-circles

Fit circles

bin/desi_fit_circles

@@ -0,0 +1,151 @@
+#!/usr/bin/env python
+
+import sys
+from astropy.table import Table
+import numpy as np
+import matplotlib.pyplot as plt
+import argparse
+
+from desimeter.circles import fit_circle
+from desimeter.transform.xy2qs import xy2uv, uv2xy
+
+parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+                                     description="""FVC image processing""")
+parser.add_argument('-i','--infile', type = str, default = None, required = True, nargs="*",
+                    help = 'path to desimeter CSV files')
+parser.add_argument('-o','--outfile', type = str, default = None, required = True,
+                    help = 'path to output CSV ASCII file')
+parser.add_argument('--plot', action = 'store_true', 
+                    help = 'plot some circles')
+
+args  = parser.parse_args()
+
+nmaxplot=100
+
+x={}
+y={}
+xexp={}
+yexp={}
+first=True
+for filename in args.infile :
+    print(filename)
+    t=Table.read(filename)
+    #print(t.dtype.names)
+    selection=(t["LOCATION"]>0)
+    location_and_pinhole=(np.array(t["LOCATION"])*10+np.array(t["PINHOLE_ID"])).astype(int)
+    if first :
+        for loc in location_and_pinhole[selection] :
+            x[loc] = []
+            y[loc] = []
+            xexp[loc] = float(t["X_FP_EXP"][location_and_pinhole==loc][0])
+            yexp[loc] = float(t["Y_FP_EXP"][location_and_pinhole==loc][0])
+            #print(loc,xexp[loc],yexp[loc])
+        first=False
+
+    for loc in location_and_pinhole[selection] :
+        ii = np.where(location_and_pinhole==loc)[0]
+        if ii.size > 1 :
+            print("several matched for LOCATION ",loc)
+            continue
+        i=ii[0]
+        if not loc in x.keys() :
+            x[loc] = []
+            y[loc] = []
+            xexp[loc] = float(t["X_FP_EXP"][location_and_pinhole==loc][0])
+            yexp[loc] = float(t["Y_FP_EXP"][location_and_pinhole==loc][0])
+
+        x[loc].append(float(t["X_FP"][i]))
+        y[loc].append(float(t["Y_FP"][i]))
+
+location_and_pinhole=np.array(list(x.keys()),dtype=int)
+location=location_and_pinhole//10
+pinhole=location_and_pinhole%10
+print("number of positioners:",np.sum(pinhole==0))
+print("number of fiducials:",np.sum(pinhole==1))
+print("number of pinholes:",np.sum(pinhole>=1))
+
+ndots=len(location_and_pinhole)
+
+theta=np.linspace(0,2*np.pi,50)
+
+
+xfp_metro=np.zeros(ndots)
+yfp_metro=np.zeros(ndots)
+xfp_meas=np.zeros(ndots)
+yfp_meas=np.zeros(ndots)
+
+
+count=0
+for iloc,loc in enumerate(x.keys()) :
+    if len(x[loc])<6 : continue
+    x[loc]=np.array(x[loc])
+    y[loc]=np.array(y[loc])
+    ii=np.where(x[loc]!=0)[0]
+    x[loc]=x[loc][ii]
+    y[loc]=y[loc][ii]
+
+    #print("{} nmeas={} rms(x)={} rms(y)={}".format(loc,ii.size,np.std(x[loc]),np.std(y[loc])))
+
+    if pinhole[iloc] == 0 and np.std(x[loc])<1. :
+        # this is a non-moving positioner, I don't use this
+        if False and args.plot and count<nmaxplot :
+            plt.figure("circles")
+            plt.plot(np.mean(x[loc]),np.mean(y[loc]),"x",color="k",markersize=12)
+        continue
+    count += 1
+
+    xc=np.median(x[loc])
+    yc=np.median(y[loc])
+
+    if pinhole[iloc] == 0 : # it's a positioner
+        # here is the fit
+        try:
+            #- 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
+
+        if iloc%100==0 :
+            print("{}/{} loc={} x={} y={} r={}".format(iloc,len(x),loc,xc,yc,r))
+        if r<0.1 : continue
+
+        if args.plot and count<nmaxplot :
+            plt.figure("circles")
+            plt.plot(x[loc],y[loc],"o")
+            plt.plot(xexp[loc],yexp[loc],"x")
+            theta=np.linspace(0,2*np.pi,50)
+            plt.plot(xc+r*np.cos(theta),yc+r*np.sin(theta),"-",color="green")
+            plt.plot(xc,yc,"+",color="green")
+
+
+    xfp_metro[iloc]=xexp[loc]
+    yfp_metro[iloc]=yexp[loc]
+    xfp_meas[iloc]=xc
+    yfp_meas[iloc]=yc
+
+dx=xfp_meas-xfp_metro
+dy=yfp_meas-yfp_metro
+dr=np.sqrt(dx**2+dy**2)
+print("median offset = {:4.1f} um",np.median(dr[dr!=0])*1000.)
+
+ii=np.where((xfp_metro!=0)&(dr<0.2))[0]
+
+# make a table out of that
+t2=Table([location[ii],pinhole[ii],xfp_metro[ii],yfp_metro[ii],xfp_meas[ii],yfp_meas[ii]],names=["LOCATION","PINHOLE_ID","X_FP_METRO","Y_FP_METRO","X_FP","Y_FP"],dtype=[int,int,float,float,float,float])
+
+t2.write(args.outfile,format="csv",overwrite=True)
+print("wrote",args.outfile)         
+
+if args.plot :
+    plt.figure("quiver",figsize=(6.1, 6))
+    plt.quiver(xfp_meas[ii],yfp_meas[ii],dx[ii],dy[ii])
+    plt.plot(xfp_meas[ii][pinhole[ii]>0],yfp_meas[ii][pinhole[ii]>0],".",c="red")
+    plt.show()
+
+
+