diff --git a/geometry.py b/geometry.py index a54783f..b1b65f2 100755 --- a/geometry.py +++ b/geometry.py @@ -248,7 +248,7 @@ class geometry: _log.debug('least square data:\nK:{}\nAA:{}'.format(K, AA)) @staticmethod - def autofocus(cam,mot,rng=(-1,1),n=30): + def autofocus(cam,mot,rng=(-1,1),n=30,saveImg=False): # cam camera object # mot motor object # rng region (min max relative to current position) to seek @@ -256,33 +256,11 @@ class geometry: # roi region of interrest to calculate sharpness # mode mode to calculate sharpness (sum/max-min/hist? of edge detection in roi) import PIL.Image - from scipy import ndimage - v=np.ndarray(shape=(len(cam),2)) - if type(cam) == list: - for i, fn in enumerate(cam): - img=PIL.Image.open(fn) - img=np.asarray(img) - s=ndimage.sobel(img) - v[i,0]=s.sum() - v[i,1]=s.std() - #fft=np.log(np.abs(np.fft.fft2(img))) - #fft=np.fft.fftshift(fft) - #s=np.array(fft.shape,dtype=np.uint16)/2 - #fft[300:700,400:800]=0 - #v[i,1]=fft.sum() - #if i&0x3==0: - # plt.figure() - # plt.imshow(fft) - fig, ax=plt.subplots() - mx=v.max(0) - mn=v.min(0) - v=(v-mn)/(mx-mn) - #ax.plot(v[:,0]) - ax.plot(v) - plt.show() - pass - else: + if mot is not None: p0=mot.get_rbv() + else: + p0=0 + if saveImg: for i,p in enumerate(np.linspace(p0+rng[0],p0+rng[1],n)): mot.move_abs(p,wait=True) pic=cam._pic# get_image() @@ -292,11 +270,55 @@ class geometry: pic=np.array(pic/scl, dtype=np.uint8) elif pic.dtype!=np.uint8: pic=np.array(pic, dtype=np.uint8) - img=PIL.Image.fromarray(pic) fn=f'/tmp/image{i:03d}.png' img.save(fn) _log.debug(f'{fn} {pic.dtype} {pic.min()} {pic.max()}') + mot.move_abs(p0) + return p0 + else: + from scipy import ndimage, signal + if type(cam) == list: + imgLst=cam + n=len(imgLst) + mtr=np.ndarray(shape=(n,)) + posLst=np.linspace(p0+rng[0], p0+rng[1], n) + for i,p in enumerate(posLst): + if type(cam)==list: + img=PIL.Image.open(imgLst[i]) + img=np.asarray(img) + else: + mot.move_abs(p, wait=True) + img=cam._pic # get_image() + img16=np.array(img, np.int16) + msk=np.array(((1, 0, -1), (2, 0, -2), (1, 0, -1)), np.int16) + sb1=signal.convolve2d(img16, msk, mode='same', boundary='fill', fillvalue=0) + sb2=signal.convolve2d(img16, msk.T, mode='same', boundary='fill', fillvalue=0) + sb=np.abs(sb1)+np.abs(sb2) + mtr[i]=sb.sum() + _log.debug(f'{i}/{p:.4g} -> {mtr[i]:.4g}') + + mx=mtr.argmax() + _log.debug(f'best focus at idx:{mx}= pos:{posLst[mx]} = metric:{mtr[mx]:.6g}') + if mx>0 and mx 255: scl=2**int(round(np.log2(mx)-8)) pic=np.array(pic/scl,dtype=np.uint8) + elif pic.dtype!=np.uint8: + pic=np.array(pic, dtype=np.uint8) except AttributeError: sim=app._camera._sim pic=cam._sim['imgSeq'][sim['imgIdx']] @@ -1582,6 +1587,11 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch) QMessageBox.about(self, "SwissMX", txt) pass + def cb_autofocus(self): + app=QApplication.instance() + geo=app._geometry + #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_testcode(self): try: @@ -1598,18 +1608,13 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch) plt.stem(x, y) plt.show(block=False) - step=2 + step=5 if step==0: vb=self.vb vb.autoRange(items=(self._goImg,)) elif step==1: testMatplotlib() elif step==2: - app=QApplication.instance() - cfg=app._cfg - geo=app._geometry - geo.autofocus(app._camera, self.tweakers['base_z']) - elif step==3: grp=pg.ItemGroup() vb.addItem(grp) obj=UsrGO.Marker((100, 100), (100, 100), mode=1) @@ -1620,7 +1625,7 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch) grp.addItem(obj) tc['grp']=grp vb.autoRange(items=(obj,)) - elif step==4: + elif step==3: grp=tc['grp'] tr=grp.transform() # UsrGO.obj_info(tr) @@ -1628,7 +1633,16 @@ Author Thierry Zamofing (thierry.zamofing@psi.ch) -.2, 1, 0, 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) #print(vb.childGroup.childItems()) pass diff --git a/workbench/autofocus.py b/workbench/autofocus.py index 4f553a0..c903fa1 100644 --- a/workbench/autofocus.py +++ b/workbench/autofocus.py @@ -4,91 +4,135 @@ Demonstrate a simple data-slicing task: given 3D data (displayed at top), select a 2D plane and interpolate data along that plane to generate a slice image (displayed at bottom). +https://www.hindawi.com/journals/js/2021/5643054/ +https://www.hindawi.com/journals/mpe/2021/8243072/ """ -## Add path to library (just for examples; you do not need this) +import logging +import numpy as np +_log=logging.getLogger(__name__) import numpy as np + +import PIL.Image +from scipy import ndimage,signal +import glob from pyqtgraph.Qt import QtCore, QtGui import pyqtgraph as pg -app=QtGui.QApplication([]) -## Create window with two ImageView widgets -win=QtGui.QMainWindow() -win.resize(800, 800) -win.setWindowTitle('pyqtgraph example: DataSlicing') -cw=QtGui.QWidget() -win.setCentralWidget(cw) -l=QtGui.QGridLayout() -cw.setLayout(l) -imv1=pg.ImageView() -imv2=pg.ImageView() -l.addWidget(imv1, 0, 0) -l.addWidget(imv2, 1, 0) -sld=QtGui.QSlider(QtCore.Qt.Horizontal) -sld.setMinimum(10) -sld.setMaximum(30) -sld.setValue(20) -sld.setTickPosition(QtGui.QSlider.TicksBelow) -sld.setTickInterval(5) +class autofocus(QtGui.QMainWindow): + def __init__(self, parent = None): + super(autofocus, self).__init__(parent) + self.resize(800, 1500) + self.setWindowTitle('pyqtgraph example: DataSlicing') + cw=QtGui.QWidget() + self.setCentralWidget(cw) + l=QtGui.QGridLayout() + cw.setLayout(l) -l.addWidget(sld, 2, 0) -win.show() + #self._imgLst=imgLst=sorted(glob.glob("../scratch/autofocus1/image*.png")) + self._imgLst=imgLst=sorted(glob.glob("../scratch/autofocus2/image*.png")) + self._metrics=mtr=np.ndarray(shape=(len(imgLst), 5)) + mtr[:]=0 + self._sld=sld=QtGui.QSlider(QtCore.Qt.Horizontal) + sld.setMinimum(0) + sld.setMaximum(len(imgLst)-1) + sld.setValue(0) + sld.setTickPosition(QtGui.QSlider.TicksBelow) + sld.setTickInterval(1) + sld.valueChanged.connect(self.cb_sld_change) -roi=pg.LineSegmentROI([[10, 64], [120, 64]], pen='r') -imv1.addItem(roi) + self._imv1=imv1=pg.ImageView() + self._imv2=imv2=pg.ImageView() -x1=np.linspace(-30, 10, 128)[:, np.newaxis, np.newaxis] -x2=np.linspace(-20, 20, 128)[:, np.newaxis, np.newaxis] -y=np.linspace(-30, 10, 128)[np.newaxis, :, np.newaxis] -z=np.linspace(-20, 20, 128)[np.newaxis, np.newaxis, :] -d1=np.sqrt(x1**2+y**2+z**2) -d2=2*np.sqrt(x1[::-1]**2+y**2+z**2) -d3=4*np.sqrt(x2**2+y[:, ::-1]**2+z**2) -data=(np.sin(d1)/d1**2)+(np.sin(d2)/d2**2)+(np.sin(d3)/d3**2) + #spl=QtGui.QSplitter(QtCore.Qt.Horizontal) + self._pw=pw=pg.PlotWidget(name='Plot1') ## giving the plots names allows us to link -import PIL.Image -from scipy import ndimage -import glob -imgLst=sorted(glob.glob("image*.png")) -v=np.ndarray(shape=(len(imgLst), 2)) - -#for i, fn in enumerate(imgLst): -# img=PIL.Image.open(fn) -# img=np.asarray(img) -# s=ndimage.sobel(img) -# v[i, 0]=s.sum() -# v[i, 1]=s.std() - #fig, ax=plt.subplots() - #mx=v.max(0) - #mn=v.min(0) - #v=(v-mn)/(mx-mn) - # ax.plot(v[:,0]) - #ax.plot(v) - #plt.show() -# pass + self._plt=plt=[] + for c in ('rgbcy'): + plt.append(pw.plot(pen=c)) -def update(): - global data, imv1, imv2, imgLst - d2=roi.getArrayRegion(data, imv1.imageItem, axes=(1, 2)) - imv2.setImage(d2) + + pw.resize(100,100) + pw.setMaximumSize(2000,200) + l.addWidget(sld, 0, 0) + l.addWidget(imv1, 1, 0) + l.addWidget(imv2, 2, 0) + l.addWidget(pw, 3, 0) + + ## Display the data + self.cb_sld_change(0,True) + mtr[1:,:]=mtr[0,:] + + self._imv1.setHistogramRange(0, 100) + self._imv1.setLevels(0, 40) + self._imv2.setHistogramRange(0, 100) + self._imv2.setLevels(0, 40) + + def cb_sld_change(self,val,auto=False): + i=self._sld.value() + _log.debug(f'{i}') + fn= self._imgLst[i] + img=PIL.Image.open(fn) + img=np.asarray(img) + #slb=ndimage.sobel(img) + img16=np.array(img,np.int16) + 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 + idx=sb[:]<20 + sbLut=sb*1 + sbLut[idx]=0 + #import numpy as np + #import matplotlib.pyplot as plt + 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) -roi.sigRegionChanged.connect(update) + #fft=np.log(np.abs(np.fft.fft2(sb))) + #fft=np.fft.fftshift(fft) + # fft[300:700,400:800]=0 + # v[i,1]=fft.sum() -## Display the data -imv1.setImage(data) -imv1.setHistogramRange(-0.01, 0.01) -imv1.setLevels(-0.003, 0.003) + self._imv1.setImage(sb,autoRange=auto,autoLevels=auto,autoHistogramRange=auto) + #self._imv2.setImage(sb,autoRange=auto,autoLevels=auto,autoHistogramRange=auto) + self._imv2.setImage(sbLut,autoRange=auto,autoLevels=auto,autoHistogramRange=auto) + mtr=self._metrics -update() + 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__': import sys + logging.basicConfig(level=logging.DEBUG, format='%(levelname)s:%(module)s:%(lineno)d:%(funcName)s:%(message)s ') + app=QtGui.QApplication([]) + af=autofocus() + af.show() if (sys.flags.interactive!=1) or not hasattr(QtCore, 'PYQT_VERSION'): QtGui.QApplication.instance().exec_() diff --git a/workbench/fiducialDetection.py b/workbench/fiducialDetection.py new file mode 100644 index 0000000..687bb6b --- /dev/null +++ b/workbench/fiducialDetection.py @@ -0,0 +1,138 @@ +# -*- coding: utf-8 -*- +""" +Demonstrate a simple data-slicing task: given 3D data (displayed at top), select +a 2D plane and interpolate data along that plane to generate a slice image +(displayed at bottom). + +https://www.hindawi.com/journals/js/2021/5643054/ +https://www.hindawi.com/journals/mpe/2021/8243072/ + +""" + +import logging +import numpy as np +_log=logging.getLogger(__name__) + +import numpy as np + +import PIL.Image +from scipy import ndimage,signal +import glob +from pyqtgraph.Qt import QtCore, QtGui +import pyqtgraph as pg + + +class fiducial(QtGui.QMainWindow): + def __init__(self, parent = None): + super(autofocus, self).__init__(parent) + self.resize(800, 1500) + self.setWindowTitle('pyqtgraph example: DataSlicing') + cw=QtGui.QWidget() + self.setCentralWidget(cw) + l=QtGui.QGridLayout() + cw.setLayout(l) + + #self._imgLst=imgLst=sorted(glob.glob("../scratch/autofocus1/image*.png")) + self._imgLst=imgLst=sorted(glob.glob("../scratch/autofocus2/image*.png")) + self._metrics=mtr=np.ndarray(shape=(len(imgLst), 5)) + mtr[:]=0 + self._sld=sld=QtGui.QSlider(QtCore.Qt.Horizontal) + sld.setMinimum(0) + sld.setMaximum(len(imgLst)-1) + sld.setValue(0) + sld.setTickPosition(QtGui.QSlider.TicksBelow) + sld.setTickInterval(1) + sld.valueChanged.connect(self.cb_sld_change) + + self._imv1=imv1=pg.ImageView() + self._imv2=imv2=pg.ImageView() + + #spl=QtGui.QSplitter(QtCore.Qt.Horizontal) + self._pw=pw=pg.PlotWidget(name='Plot1') ## giving the plots names allows us to link + + self._plt=plt=[] + for c in ('rgbcy'): + plt.append(pw.plot(pen=c)) + + + + pw.resize(100,100) + pw.setMaximumSize(2000,200) + l.addWidget(sld, 0, 0) + l.addWidget(imv1, 1, 0) + l.addWidget(imv2, 2, 0) + l.addWidget(pw, 3, 0) + + ## Display the data + self.cb_sld_change(0,True) + mtr[1:,:]=mtr[0,:] + + self._imv1.setHistogramRange(0, 100) + self._imv1.setLevels(0, 40) + self._imv2.setHistogramRange(0, 100) + self._imv2.setLevels(0, 40) + + def cb_sld_change(self,val,auto=False): + i=self._sld.value() + _log.debug(f'{i}') + fn= self._imgLst[i] + img=PIL.Image.open(fn) + img=np.asarray(img) + #slb=ndimage.sobel(img) + img16=np.array(img,np.int16) + 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 + idx=sb[:]<20 + sbLut=sb*1 + sbLut[idx]=0 + #import numpy as np + #import matplotlib.pyplot as plt + 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) + + + #fft=np.log(np.abs(np.fft.fft2(sb))) + #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) + #self._imv2.setImage(sb,autoRange=auto,autoLevels=auto,autoHistogramRange=auto) + self._imv2.setImage(sbLut,autoRange=auto,autoLevels=auto,autoHistogramRange=auto) + 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__': + import sys + logging.basicConfig(level=logging.DEBUG, format='%(levelname)s:%(module)s:%(lineno)d:%(funcName)s:%(message)s ') + + app=QtGui.QApplication([]) + fd=fiducial() + fd.show() + if (sys.flags.interactive!=1) or not hasattr(QtCore, 'PYQT_VERSION'): + QtGui.QApplication.instance().exec_()