Merge pull request #152 from desihub/spotmatch

Update to spotmatch usage

bin/desi_fvc_proc

@@ -75,6 +75,7 @@ parser.add_argument('--turbulence-correction', action = 'store_true',
                     help = 'turbulence correction assuming offsets between measured and expected coordinates are not spatially correlated. Only used in conjunction with option --expected-positions.')
 parser.add_argument('--turbulence-correction-with-pol', action = 'store_true',
                     help = 'turbulence correction using local polynomial fit instead of Gaussian processes.')
+parser.add_argument('--spotmatch-match-radius-pixels', type = float, default = 70, required = False, help = "spotmatch matching radius")
 parser.add_argument('--make-directory', action='store_true',
                     help='Make directory for output if needed.')
 
@@ -177,7 +178,7 @@ for seqid,spots in enumerate(spots_list) :
         continue
 
     if args.use_spotmatch :
-        spots = spotmatch(spots["XPIX"],spots["YPIX"])
+        spots = spotmatch(spots["XPIX"],spots["YPIX"],match_radius_pixels=args.spotmatch_match_radius_pixels)
         # drop pinholes, keep only center
         spots = spots[(spots["PINHOLE_ID"]==0)|(spots["PINHOLE_ID"]==99)]
         n_matched_fiducials = np.sum(spots['PINHOLE_ID'] == 99)
@@ -215,6 +216,13 @@ for seqid,spots in enumerate(spots_list) :
     if args.expected_positions is not None :
         log.info("reading expected positions in {}".format(args.expected_positions))
         expected_pos = Table.read(args.expected_positions) # auto-magically guess format
+
+        if (not "X_FP" in expected_pos.keys()) and "X_FP_EXP" in expected_pos.keys() :
+            log.warning("Rename X_FP_EXP,Y_FP_EXP -> X_FP,Y_FP")
+            expected_pos.rename_column("X_FP_EXP","X_FP")
+            expected_pos.rename_column("Y_FP_EXP","Y_FP")
+
+
         if not "X_FP" in expected_pos.keys() :
             if "EXP_Q_0" in  expected_pos.keys() :
                 log.info("EXP_Q_0,EXP_S_0 -> X_FP,Y_FP")
@@ -240,31 +248,70 @@ for seqid,spots in enumerate(spots_list) :
         ii = np.where(expected_pos["PINHOLE_ID"]==0)[0]
         expected_pos = expected_pos[:][ii]
 
+
     # select spots that are not already matched
     selection  = (spots["LOCATION"]==-1)
 
-    # match
-    indices_of_expected_pos,distances = match_same_system(spots["X_FP"][selection],spots["Y_FP"][selection],expected_pos["X_FP"],expected_pos["Y_FP"])
-
-    is_matched = (distances<args.max_match_distance)&(indices_of_expected_pos>=0)
-    ii=np.where(selection)[0]
-    selection[ii]          &=  is_matched
-    indices_of_expected_pos = indices_of_expected_pos[is_matched]
-    distances               = distances[is_matched]
-
-    # add columns after matching fibers
-    for k1,k2 in zip(["X_FP","Y_FP"],["X_FP_EXP","Y_FP_EXP"]) :
-        if k2 not in spots.keys() : spots[k2] = np.zeros(len(spots))
-        spots[k2][selection]=expected_pos[k1][indices_of_expected_pos]
-    for k in ["EXP_Q_0","EXP_S_0","PETAL_LOC","DEVICE_LOC","DEVICE_ID","DEVICE_TYPE","LOCATION"] :
-        if k in expected_pos.keys() :
-            if k not in spots.keys() :
-                if k in ["DEVICE_ID","DEVICE_TYPE"] :
-                    spots[k] = np.repeat("None           ",len(spots))
-                else :
-                    spots[k] = np.zeros(len(spots))
-            spots[k][selection]=expected_pos[k][indices_of_expected_pos]
-
+    if args.use_spotmatch :
+        spots = spotmatch(spots["XPIX"],spots["YPIX"],expected_x_fp=expected_pos["X_FP"],expected_y_fp=expected_pos["Y_FP"],expected_location=expected_pos["LOCATION"],fvc2fp=tx,match_radius_pixels=args.spotmatch_match_radius_pixels)
+        #spots = spotmatch(spots["XPIX"],spots["YPIX"],expected_x_fp=expected_pos["X_FP"],expected_y_fp=expected_pos["Y_FP"],expected_location=expected_pos["LOCATION"],fvc2fp=None,match_radius_pixels=args.spotmatch_match_radius_pixels)
+
+
+        # add info
+        nspots = len(spots["LOCATION"])
+        for k in ["X_FP","Y_FP","X_FP_EXP","Y_FP_EXP","X_FP_METRO","Y_FP_METRO"] :
+            if k not in spots.dtype.names :
+                spots[k] = np.zeros(nspots,dtype=float)
+
+        spots["X_FP"],spots["Y_FP"] = tx.fvc2fp(spots["XPIX"],spots["YPIX"])
+
+        is_matched = (spots["LOCATION"]>=0)
+        loc2i = {loc:i for i,loc in enumerate(spots["LOCATION"])}
+
+        ii=[]
+        jj=[]
+        for j,loc in enumerate(expected_pos["LOCATION"]) :
+            if loc in loc2i :
+                ii.append(loc2i[loc])
+                jj.append(j)
+        spots["X_FP_EXP"][ii]=expected_pos["X_FP"][jj]
+        spots["Y_FP_EXP"][ii]=expected_pos["Y_FP"][jj]
+
+        metrology = load_metrology()
+        ii=[]
+        jj=[]
+        for j,loc in enumerate(metrology["LOCATION"]) :
+            if loc in loc2i :
+                ii.append(loc2i[loc])
+                jj.append(j)
+        spots["X_FP_METRO"][ii]=metrology["X_FP"][jj]
+        spots["Y_FP_METRO"][ii]=metrology["Y_FP"][jj]
+
+
+    else :
+
+        # match
+        indices_of_expected_pos,distances = match_same_system(spots["X_FP"][selection],spots["Y_FP"][selection],expected_pos["X_FP"],expected_pos["Y_FP"])
+        is_matched = (distances<args.max_match_distance)&(indices_of_expected_pos>=0)
+        ii=np.where(selection)[0]
+        selection[ii]          &=  is_matched
+        indices_of_expected_pos = indices_of_expected_pos[is_matched]
+        distances               = distances[is_matched]
+
+        # add columns after matching fibers
+        for k1,k2 in zip(["X_FP","Y_FP"],["X_FP_EXP","Y_FP_EXP"]) :
+            if k2 not in spots.keys() : spots[k2] = np.zeros(len(spots))
+            spots[k2][selection]=expected_pos[k1][indices_of_expected_pos]
+        for k in ["EXP_Q_0","EXP_S_0","PETAL_LOC","DEVICE_LOC","DEVICE_ID","DEVICE_TYPE","LOCATION"] :
+            if k in expected_pos.keys() :
+                if k not in spots.keys() :
+                    if k in ["DEVICE_ID","DEVICE_TYPE"] :
+                        spots[k] = np.repeat("None           ",len(spots))
+                    else :
+                        spots[k] = np.zeros(len(spots))
+                spots[k][selection]=expected_pos[k][indices_of_expected_pos]
+
+
     if args.expected_positions is not None and ( args.turbulence_correction or args.turbulence_correction_with_pol ):
 
         if args.turbulence_correction_with_pol :
@@ -285,12 +332,13 @@ for seqid,spots in enumerate(spots_list) :
             # apply it anyway because there might be a good reason for this
         spots["X_FP"][selection] = new_x
         spots["Y_FP"][selection] = new_y
-
-    # for spots with metrology X_FP_EXP=X_FP_METRO
-    selection = (spots["X_FP_METRO"]!=0)
-    spots["X_FP_EXP"][selection]=spots["X_FP_METRO"][selection]
-    selection = (spots["Y_FP_METRO"]!=0)
-    spots["Y_FP_EXP"][selection]=spots["Y_FP_METRO"][selection]
+
+    if "X_FP_METRO" in spots.dtype.names :
+        # for spots with metrology X_FP_EXP=X_FP_METRO
+        selection = (spots["X_FP_METRO"]!=0)
+        spots["X_FP_EXP"][selection]=spots["X_FP_METRO"][selection]
+        selection = (spots["Y_FP_METRO"]!=0)
+        spots["Y_FP_EXP"][selection]=spots["Y_FP_METRO"][selection]
 
     # write transfo
     if args.output_transform is not None :
@@ -300,7 +348,6 @@ for seqid,spots in enumerate(spots_list) :
             tx.write_jsonfile(args.output_transform)
             print("wrote transform in {}".format(args.output_transform))
 
-
     if args.field_model is not None :
         log.info("Reading field model in {}".format(args.field_model))
         with open(args.field_model) as file :
@@ -312,14 +359,11 @@ for seqid,spots in enumerate(spots_list) :
         spots["RA"][ii]  = ra
         spots["DEC"][ii] = dec
 
-
-
     if not args.sequence :
         outfile=args.outfile
     else :
         outfile=args.outfile.replace(".csv","-F{:04d}.csv".format(seqid))
 
-
     # write spots
     spots.write(outfile,format="csv",overwrite=True)
     print("wrote {}".format(outfile))

bin/plot_fvc_residuals

@@ -47,7 +47,7 @@ else :
     ym = table["Y_FP_METRO"][ii]
     pid = table["PINHOLE_ID"][ii]
 
-selection = (table["LOCATION"]>0)&(table["X_FP_EXP"]!=0)&(table["PINHOLE_ID"]>0)
+selection = (table["LOCATION"]>=0)&(table["X_FP_EXP"]!=0)&(table["PINHOLE_ID"]>0)
 number_of_fiducials = np.unique(table["LOCATION"][selection]).size
 print("Number of fiducials that are ON=",number_of_fiducials)
 
@@ -58,8 +58,12 @@ fig = plt.figure(figsize=(6,6))
 a = plt.subplot(1,1,1)
 a.set_title(os.path.basename(filename))
 
-a.plot(table["X_FP"],table["Y_FP"],".",alpha=0.5,label="all spots")
+a.plot(table["X_FP"],table["Y_FP"],".",alpha=0.5,label="all measured spots")
+#a.plot(table["X_FP_EXP"],table["Y_FP_EXP"],"x",alpha=0.5,label="X,Y_FP_EXP")
+#a.plot(table["X_FP_METRO"],table["Y_FP_METRO"],"+",alpha=0.5,label="X,Y_FP_METRO")
+#plt.legend(loc="upper left")
 
+#plt.show()
 if not args.expected : # plotting match to fiducials
     # plotting all of FIF and GIF
     metrology = load_metrology()
@@ -75,40 +79,40 @@ if args.expected :
 else :
     label = "matched metrology"
     marker = "o"
+
 a.scatter(xm,ym,marker=marker,edgecolors="orange",facecolors="none",label=label)
+
 if args.expected :
-    from desimeter.transform.fvc2fp import FVC2FP
-    from desimeter.io import fvc2fp_filename,load_metrology
-    metrology = load_metrology()
-    input_tx = FVC2FP.read_jsonfile(fvc2fp_filename())
-    xpix=np.array([2000.,]) ; ypix=np.array([0.,])
-    xfp1,yfp1 = input_tx.fvc2fp(xpix,ypix)
-    xfp2,yfp2 = input_tx.fvc2fp(xpix+1,ypix)
-    pixel2fp  = np.hypot(xfp2-xfp1, yfp2-yfp1)[0] # mm
-
-    match_radius_pixels = 70
-    match_radius_mm = pixel2fp*match_radius_pixels
-    print("draw match radius = {} pixels = {} mm".format(match_radius_pixels,match_radius_mm))
+    radius_mm = 6. #pixel2fp*match_radius_pixels
     angle=np.linspace(0,2*np.pi,50)
-    ca=match_radius_mm*np.cos(angle)
-    sa=match_radius_mm*np.sin(angle)
-    for xxm,yym,xx,yy,ppid in zip(xm,ym,x,y,pid) :
-        if ppid==0 : # a positioner
-            a.plot(xxm+ca,yym+sa,"-",color="orange")
-            a.plot([xxm,xx],[yym,yy],"-",color="orange")
+    ca=radius_mm*np.cos(angle)
+    sa=radius_mm*np.sin(angle)
+    if "X_FP_METRO" in table.dtype.names and "DEVICE_TYPE" in table.dtype.names :
+        ok=table["DEVICE_TYPE"]=="POS"
+        for xx,yy in zip(table["X_FP_METRO"][ok],table["Y_FP_METRO"][ok]) :
+            a.plot(xx+ca,yy+sa,"-",color="gray",alpha=0.4)
+    # also display unmatched
+    metrology = load_metrology()
+    unmatched = ~np.in1d(metrology["LOCATION"],table["LOCATION"])
+    unmatched_positioners = unmatched & (metrology["DEVICE_TYPE"]=="POS")
+    print("unmatched positioners = {}".format(list(metrology["DEVICE_ID"][unmatched_positioners])))
+    for xx,yy in zip(metrology["X_FP"][unmatched_positioners],metrology["Y_FP"][unmatched_positioners]):
+        a.plot(xx+ca,yy+sa,"-",color="red",alpha=0.4)
+
 
 dx=xm-x
 dy=ym-y
 
 if args.posids is not None :
     posids=args.posids.split(",")
+    metrology = load_metrology()
     for posid in posids :
-        i=np.where(table["DEVICE_ID"]==posid)[0]
+        i=np.where(metrology["DEVICE_ID"]==posid)[0]
         if len(i)==0 :
-            print(posid,"not matched")
+            print(posid,"not in metrology")
             continue
         i=i[0]
-        plt.plot(table["X_FP"][i],table["Y_FP"][i],"o",label=posid)
+        plt.plot(metrology["X_FP"][i],metrology["Y_FP"][i],"o",alpha=0.6,markersize=10,label=posid)
 
 
 
@@ -119,13 +123,16 @@ elif jj.size > 0 : # sadly
     print("Large residuals")
     for j in jj :
         i=ii[j]
-        line = "dx={:4.3f}mm dy={:4.3f}mm".format(dx[j],dy[j])
-        for k in table.dtype.names :
+        line=""
+        #for k in table.dtype.names :
+        for k in ["DEVICE_ID","X_FP","Y_FP"] :
+            if not k in table.dtype.names : continue
             if k=="XERR" or k=="YERR" : continue
             try :
-                line += " {}={:4.3f}".format(k,table[k][i])
+                line += " {}={:4.2f}".format(k,table[k][i])
             except ValueError :
                 line += " {}={}".format(k,table[k][i])
+        line += " dX={:4.2f}mm dY={:4.2f}mm".format(dx[j],dy[j])
         print(line)
         label=None
         if j==jj[0]: label="large residual"
@@ -138,19 +145,16 @@ a.set_ylabel("Y_FP (mm)")
 
 dist=np.sqrt(dx**2+dy**2)
 
-rms=np.sqrt(np.mean(dist[dist<0.2]**2))
-blabla="rms(2D) = {:3.1f} um".format(rms*1000.)
+rms1=np.sqrt(np.mean(dist[dist<0.2]**2))
+blabla="rms(2D) (<200um):{:3.1f} um".format(rms1*1000.)
+rms2=np.sqrt(np.mean(dist[dist<6.]**2))
+blabla+= ", all:{:3.1f} mm".format(rms2)
 print(blabla)
 a.legend(loc="upper left",title=blabla)
 
 a2=fig.add_subplot(666)
 a2.hist(dist[dist<0.1]*1000.)
 a2.set_xlabel("dist. (um)")
-if args.expected :
-    plt.figure("dist")
-    plt.hist(dist)
-    plt.xlabel("dist. (mm)")
-    print("max(dist)=",np.max(dist))
 
 if False :
     print("Write coordinates of missing or incorrect fiducials")

py/desimeter/spotmatch.py

@@ -118,7 +118,7 @@ def _write_spotmatch_targets_file(x_fp,y_fp,location,filename,fvc2fp=None) :
 
 
     xpix,ypix = fvc2fp.fp2fvc(x_fp,y_fp)
-
+    
     with open(filename,"w") as ofile :
         for i in range(xpix.size) :
             ofile.write("{} {:4.3f} {:4.3f} 12.000 0.001 {}\n".format(location[i],xpix[i],ypix[i],flags[i]))
@@ -233,7 +233,7 @@ def _write_spotmatch_reference_pos_file(filename,fvc2fp=None) :
 
     print("wrote",filename)
 
-def spotmatch(xpix,ypix,expected_x_fp=None,expected_y_fp=None,expected_location=None,verbose=0,match_radius_pixels=70) :
+def spotmatch(xpix,ypix,expected_x_fp=None,expected_y_fp=None,expected_location=None,verbose=0,match_radius_pixels=70,fvc2fp=None) :
     """
     Wrapper to spotmatch in desimeter. Calls the C executable 'match_positions' that has to be in the path. All inputs
     to match_positions are generated in this call using desimeter metrology table, the default fvc2fp transform and the inputs
@@ -248,6 +248,7 @@ def spotmatch(xpix,ypix,expected_x_fp=None,expected_y_fp=None,expected_location=
        expected_location :  1D numpy array with targets location index ( = petal_loc*10000+device_loc)
        verbose : 0 or 1
        match_radius_pixels : match radius in pixels
+       fvc2fp : desimeter.transform.fvc2fp.FVC2FP object
 
     returns :
        an astropy.table.Table object with at least the columns
@@ -277,7 +278,8 @@ def spotmatch(xpix,ypix,expected_x_fp=None,expected_y_fp=None,expected_location=
 
     tmp_dir=tempfile.gettempdir()
 
-    fvc2fp = FVC2FP.read(fvc2fp_filename())
+    if fvc2fp is None :
+        fvc2fp = FVC2FP.read(fvc2fp_filename())
     exp_pixel_scale = _compute_pixel_scale(fvc2fp)