very promissing fiducial detection

geometry.py

@@ -16,10 +16,18 @@ modes:
   0x02: update_optical_center
 '''
 import logging
-import numpy as np
-
 _log=logging.getLogger(__name__)
 
+import numpy as np
+import PIL.Image
+try:
+  from scipy import ndimage, signal
+except ImportError as e:
+  _log.warning(e)
+try:
+  import cv2 as cv
+except ImportError as e:
+  _log.warning(e)
 
 class geometry:
 
@@ -255,7 +263,6 @@ class geometry:
     # n     number of images to take in region
     # roi   region of interrest to calculate sharpness
     # mode  mode to calculate sharpness (sum/max-min/hist?  of edge detection in roi)
-    import PIL.Image
     if mot is not None:
       p0=mot.get_rbv()
     else:
@@ -277,7 +284,6 @@ class geometry:
         mot.move_abs(p0)
         return p0
     else:
-      from scipy import ndimage, signal
       if type(cam) == list:
         imgLst=cam
         n=len(imgLst)
@@ -321,6 +327,45 @@ class geometry:
       return p
     pass
 
+  @staticmethod
+  def find_fiducial(cam,sz=(210,210),brd=(20,20)):
+    if type(cam)==str:
+      img=PIL.Image.open(cam)
+      img=np.asarray(img)
+    else:
+      img=cam._pic  # get_image()
+      img16=np.array(img, np.int16)
+
+    fid=np.ones((sz[1]+2*brd[1],sz[0]+2*brd[0]),dtype=np.uint8)*255
+    fid[brd[1]:sz[1]+brd[1],brd[0]:sz[0]+brd[0]]=0
+    mask=np.ones((sz[1]+2*brd[1],sz[0]+2*brd[0]),dtype=np.uint8)*255
+    mask[2*brd[1]:sz[1],2*brd[0]:sz[0]]=0
+    #https://docs.opencv.org/4.5.2/d4/dc6/tutorial_py_template_matching.html
+    #res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED )
+    res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED,mask=mask)
+
+    h,w=img.shape
+    fh2,fw2=fid.shape
+    fw2//=2
+    fh2//=2
+    mtr=np.ndarray((5,2),np.uint16)
+    corr=np.ndarray((5,),np.float32)
+    for i in range(5):
+      p=np.unravel_index(res.argmax(), res.shape)
+      corr[i]=res[p]
+      mtr[i,:]=p
+      y0=max(p[0]-fh2,0)
+      y1=min(p[0]+fh2,h)
+      x0=max(p[1]-fw2,0)
+      x1=min(p[1]+fw2,w)
+      res[y0:y1,x0:x1]*=.5
+    pos=mtr.mean(0)[::-1]+(fw2,fh2)
+    crm=corr.mean()
+    _log.debug(f'position: {pos} correlation:{crm}')
+    return (pos,crm)
+
+
+
   def pix2pos(self, p, zoom=None):
     # returns the position m(x,y) in meter relative to the optical center at a given zoom level of the pixel p(x,y)
     # if zoom=None, the last zoom value is used
@@ -627,6 +672,8 @@ if __name__=="__main__":
     import glob
     imgLst=sorted(glob.glob("scratch/autofocus2/image*.png"))
     geometry.autofocus(imgLst,None)
+  if args.mode&0x10:
+    geometry.find_fiducial("scratch/fiducial/image001.png")
 
 
 # pix2pos="[[1.0, 200.0, 400.0, 600.0, 800.0], [[[0.001182928623952055, -2.6941995127711305e-05], [-4.043716694634124e-05, -0.0011894314084263825]], [[0.0007955995220142541, -3.175003727901119e-05], [-2.0896601103372113e-05, -0.0008100805094631365]], [[0.00048302539335378367, -1.1661121407652543e-05], [-2.0673746995751222e-05, -0.0004950857899461772]], [[0.00028775903460576395, -1.3762555219494508e-05], [-9.319936861519268e-06, -0.0002889214488565999]], [[0.0001788819256630411, -6.470841493681516e-06], [-2.0336605088889967e-06, -0.0001831131753499113]]]]"

pyqtUsrObj.py

@@ -527,6 +527,17 @@ class FixTargetFrame(pg.ROI):
 
 ## Start Qt event loop unless running in interactive mode or using pyside.
 if __name__=='__main__':
+  def set_fiducial(pic):
+    # fiducial test
+    f=np.array(((0, 0, 0, 0, 0),
+                (0, 1, 1, 1, 0),
+                (0, 1, 0, 0, 0),
+                (0, 1, 1, 0, 0),
+                (0, 1, 0, 0, 0),
+                (0, 0, 0, 0, 0),), pic.dtype)
+    pic[0:6, 0:5]=f*pic.max()
+
+
   # TODO: pg.ItemGroup does not support bounding box and therefore vb.autoRange() does not work
   class ItemGroup(pg.ItemGroup):
     # own item group that supports bounding rect
@@ -667,7 +678,7 @@ if __name__=='__main__':
   arr[:, 50]=10
   arr+=np.sin(np.linspace(0, 20, 100)).reshape(1, 100)
   arr+=np.random.normal(size=(100, 100))
-
+  set_fiducial(arr)
   # add an arrow for asymmetry
   arr[10, :50]=10
   arr[9:12, 44:48]=10

swissmx.py

@@ -496,6 +496,10 @@ class WndSwissMx(QMainWindow, Ui_MainWindow):
     action.triggered.connect(self.cb_autofocus)
     self.toolBar.addAction(action)
 
+    action = QAction(icon, "Find\nFiducial", self)
+    action.triggered.connect(self.cb_find_fiducial)
+    self.toolBar.addAction(action)
+
     action = QAction(icon, "Test\nCode", self)
     action.triggered.connect(self.cb_testcode)
     self.toolBar.addAction(action)
@@ -1593,7 +1597,17 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch)
       #geo.autofocus(app._camera, self.tweakers['base_z'],rng=(-1, 1), n=30,saveImg=True)
       geo.autofocus(app._camera, self.tweakers['base_z'],rng=(-1, 1), n=10)
 
+  def cb_find_fiducial(self):
+      app=QApplication.instance()
+      geo=app._geometry
+      #geo.autofocus(app._camera, self.tweakers['base_z'],rng=(-1, 1), n=30,saveImg=True)
+      pos,corr=geo.find_fiducial(app._camera, sz=(210,210),brd=(20,20))
+      pass
+
   def cb_testcode(self):
+    app=QApplication.instance()
+    cfg=app._cfg
+    geo=app._geometry
     try:
       tc=self._testCode
       tc['idx']+=1
@@ -1608,13 +1622,14 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch)
       plt.stem(x, y)
       plt.show(block=False)
 
-    step=5
+    step=6
     if step==0:
       vb=self.vb
       vb.autoRange(items=(self._goImg,))
     elif step==1:
       testMatplotlib()
     elif step==2:
+      vb=self.vb
       grp=pg.ItemGroup()
       vb.addItem(grp)
       obj=UsrGO.Marker((100, 100), (100, 100), mode=1)
@@ -1634,15 +1649,11 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch)
                    0, 0, 1)
       grp.setTransform(tr)
     elif step==4:
-      app=QApplication.instance()
-      cfg=app._cfg
-      geo=app._geometry
       geo.autofocus(app._camera, self.tweakers['base_z'],rng=(-1, 1), n=30,saveImg=True)
     elif step==5:
-      app=QApplication.instance()
-      cfg=app._cfg
-      geo=app._geometry
       geo.autofocus(app._camera, self.tweakers['base_z'],rng=(-1, 1), n=10)
+    elif step==6:
+      geo.find_fiducial(app._camera)
     #print(vb.childGroup.childItems())
     pass
 

workbench/fiducialDetection.py

@@ -14,17 +14,50 @@ import numpy as np
 _log=logging.getLogger(__name__)
 
 import numpy as np
-
+import cv2 as cv
 import PIL.Image
 from scipy import ndimage,signal
-import glob
+import glob, os
 from pyqtgraph.Qt import QtCore, QtGui
 import pyqtgraph as pg
+pg.setConfigOptions(imageAxisOrder='row-major')
+
+
+def set_fiducial(pic):
+  # fiducial test
+  f=np.array(((0, 0, 0, 0, 0),
+              (0, 1, 1, 1, 0),
+              (0, 1, 0, 0, 0),
+              (0, 1, 1, 0, 0),
+              (0, 1, 0, 0, 0),
+              (0, 0, 0, 0, 0),), pic.dtype)
+  pic[0:6, 0:5]=f*pic.max()
+
+def wtestimages():
+  fn=os.path.expanduser('~/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/SwissMX/scratch/fiducial/image{idx:03d}.png')
+  img=np.kron(np.array([[1, 0] * 5, [0, 1] * 5] * 6,dtype=np.uint8), np.ones((100, 100),dtype=np.uint8)*255)
+  set_fiducial(img)
+  pil_img=PIL.Image.fromarray(img)
+  pil_img.save(fn.format(idx=2))
+
+  img=np.ndarray((1000,1200),dtype=np.uint8)
+  img[:]=255
+  #img=np.zeros((1200,1000),dtype=np.uint8)
+  fid=np.array([[0, 1] * 2, [1, 0] * 2] * 2,dtype=np.uint8)[:3,:3]
+  fid=np.kron(fid, np.ones((210, 210),dtype=np.uint8)*255)
+  img[100:100+fid.shape[0],200:200+fid.shape[1]]=fid
+  set_fiducial(img)
+  pil_img=PIL.Image.fromarray(img)
+  pil_img.save(fn.format(idx=3))
+
+  img=np.ndarray((1200,1000))
+
+wtestimages()
 
 
 class fiducial(QtGui.QMainWindow):
   def __init__(self, parent = None):
-    super(autofocus, self).__init__(parent)
+    super(fiducial, self).__init__(parent)
     self.resize(800, 1500)
     self.setWindowTitle('pyqtgraph example: DataSlicing')
     cw=QtGui.QWidget()
@@ -33,7 +66,7 @@ class fiducial(QtGui.QMainWindow):
     cw.setLayout(l)
 
     #self._imgLst=imgLst=sorted(glob.glob("../scratch/autofocus1/image*.png"))
-    self._imgLst=imgLst=sorted(glob.glob("../scratch/autofocus2/image*.png"))
+    self._imgLst=imgLst=sorted(glob.glob("../scratch/fiducial/*.png"))
     self._metrics=mtr=np.ndarray(shape=(len(imgLst), 5))
     mtr[:]=0
     self._sld=sld=QtGui.QSlider(QtCore.Qt.Horizontal)
@@ -69,8 +102,10 @@ class fiducial(QtGui.QMainWindow):
 
     self._imv1.setHistogramRange(0, 100)
     self._imv1.setLevels(0, 40)
-    self._imv2.setHistogramRange(0, 100)
+    #self._imv2.setHistogramRange(0, 100)
-    self._imv2.setLevels(0, 40)
+    #self._imv2.setLevels(0, 40)
+    self._imv2.setHistogramRange(0, 1)
+    self._imv2.setLevels(0, 1)
 
   def cb_sld_change(self,val,auto=False):
     i=self._sld.value()
@@ -78,53 +113,88 @@ class fiducial(QtGui.QMainWindow):
     fn= self._imgLst[i]
     img=PIL.Image.open(fn)
     img=np.asarray(img)
-    #slb=ndimage.sobel(img)
+    #pil_img=PIL.Image.fromarray(img)
-    img16=np.array(img,np.int16)
+    #pil_img.save('../scratch/fiducial/image002.png')
-    s1=np.array(((1,0,-1),(2,0,-2),(1,0,-1)),np.int16)
-    sb1=signal.convolve2d(img16, s1, mode='same', boundary='fill', fillvalue=0)
-    sb2=signal.convolve2d(img16, s1.T, mode='same', boundary='fill', fillvalue=0)
-    sb=np.abs(sb1)+np.abs(sb2)
 
-    #remove irrelevant low values
+    #fn='/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/SwissMX/scratch/fiducial/zoom0400.png'
-    idx=sb[:]<20
+    #pil_img=PIL.Image.open(fn+'_')
-    sbLut=sb*1
+    #img=np.asarray(pil_img)
-    sbLut[idx]=0
+    #img=np.array(img,dtype=np.uint8)
-    #import numpy as np
+    #pil_img=PIL.Image.fromarray(img)
-    #import matplotlib.pyplot as plt
+    #pil_img.save(fn)
-    mx=sb.max()+1
-    lut=(np.sin(np.arange(mx)/mx*np.pi-np.pi/2)+1)*128
-    #lut=((np.arcsin(np.arange(mx)/(mx-1)*2-1)/np.pi)+.5)*255
-    #lut=np.array(lut*255,np.uint16)
-    #sbLut=lut[sb]
-    #import matplotlib.pyplot as plt
-    #plt.plot(lut);plt.show()
-    #img2=ndimage.grey_dilation(sb,size=(5,5)) #, size=None, footprint=None, structure=None, output=None, mode='reflect', cval=0.0, origin=0)
-    img2=ndimage.grey_closing(sb,size=(25,25)) #, size=None, footprint=None, structure=None, output=None, mode='reflect', cval=0.0, origin=0)
 
-
+    #img16=np.array(img,np.int16)
-    #fft=np.log(np.abs(np.fft.fft2(sb)))
+    #fft=np.log(np.abs(np.fft.fft2(img)))
     #fft=np.fft.fftshift(fft)
-    # fft[300:700,400:800]=0
-    # v[i,1]=fft.sum()
 
-    self._imv1.setImage(sb,autoRange=auto,autoLevels=auto,autoHistogramRange=auto)
+    fid=np.array([[0, 1]*2, [1, 0]*2]*2, dtype=np.uint8)[:3, :3]
-    #self._imv2.setImage(sb,autoRange=auto,autoLevels=auto,autoHistogramRange=auto)
+    fid=np.kron(fid, np.ones((100, 100), dtype=np.uint8)*255)
-    self._imv2.setImage(sbLut,autoRange=auto,autoLevels=auto,autoHistogramRange=auto)
+    #res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED )
-    mtr=self._metrics
 
-    mtr[i, 0]=img.std()
+    #fid=np.ones((250,250),dtype=np.uint8)*255
-    mtr[i, 1]=sb.sum()
+    #fid[20:230,20:230]=0
-    mtr[i, 2]=sb.std()
-    mtr[i, 3]=sbLut.sum()#sb.std()
-    mtr[i, 4]=sbLut.std()#img2.sum()
 
-    mx=mtr.max(0)
+    sz=(90, 90); brd=(20, 20) # zoom 001
-    mn=mtr.min(0)
+    #sz=(130, 130); brd=(20, 20) # zoom 200
-    mtr=(mtr-mn)/(mx-mn)
+    #sz=(210, 210); brd=(20, 20) # zoom 400
+    fid=np.ones((sz[1]+2*brd[1],sz[0]+2*brd[0]),dtype=np.uint8)*255
+    fid[brd[1]:sz[1]+brd[1],brd[0]:sz[0]+brd[0]]=0
+    mask=np.ones((sz[1]+2*brd[1],sz[0]+2*brd[0]),dtype=np.uint8)*255
+    mask[2*brd[1]:sz[1],2*brd[0]:sz[0]]=0
 
-    _log.debug(f'{i} {mtr[i,:]}')
+
-    for i in range(mtr.shape[1]):
+    #https://docs.opencv.org/4.5.2/d4/dc6/tutorial_py_template_matching.html
-      self._plt[i].setData(mtr[:,i])
+    #res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED )
+    res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED,mask=mask)
+    corr_u8=np.array(res*255,np.uint8)
+    h,w=img.shape
+    fh2,fw2=fid.shape
+    fw2//=2
+    fh2//=2
+    mtr=np.ndarray((5,2),np.uint16)
+    corr=np.ndarray((5,),np.float32)
+    for i in range(5):
+      p=np.unravel_index(res.argmax(), res.shape)
+      corr[i]=res[p]
+      mtr[i,:]=p
+      y0=max(p[0]-fh2,0)
+      y1=min(p[0]+fh2,h)
+      x0=max(p[1]-fw2,0)
+      x1=min(p[1]+fw2,w)
+      res[y0:y1,x0:x1]*=.5
+      res[p[0],p[1]]=0
+
+    corr_u8=np.array(res*255,np.uint8)
+
+    ctr=mtr.mean(0,dtype=np.int32)+(fw2,fh2)
+    try:
+      img=img*1
+      img[ctr[0]-5:ctr[0]+5,ctr[1]-5:ctr[1]+5]=255
+      img[ctr[0],:]=0
+      img[:,ctr[1]]=0
+    except IndexError as e:
+      pass
+    self._imv1.setImage(img,autoRange=auto,autoLevels=auto,autoHistogramRange=auto)
+    self._imv2.setImage(res,autoRange=auto,autoLevels=auto,autoHistogramRange=auto)
+    pos=mtr.mean(0)[::-1]+(fw2,fh2)
+    _log.debug(f'position: {pos} correlation:{corr.mean()}')
+
+
+    #mtr=self._metrics
+
+    #mtr[i, 0]=img.std()
+    #mtr[i, 1]=sb.sum()
+    #mtr[i, 2]=sb.std()
+    #mtr[i, 3]=sbLut.sum()#sb.std()
+    #mtr[i, 4]=sbLut.std()#img2.sum()
+
+    #mx=mtr.max(0)
+    #mn=mtr.min(0)
+    #mtr=(mtr-mn)/(mx-mn)
+
+    #_log.debug(f'{i} {mtr[i,:]}')
+    #for i in range(mtr.shape[1]):
+    #  self._plt[i].setData(mtr[:,i])
 
 ## Start Qt event loop unless running in interactive mode.
 if __name__=='__main__':