From fda9750710f04d6f82db31b3f2c4252674cd0481 Mon Sep 17 00:00:00 2001 From: Thierry Zamofing Date: Thu, 19 Apr 2018 11:20:16 +0200 Subject: [PATCH] good working findxtal.py --- python/findxtal.py | 179 +++++++++++++++++++++++++++++++++------------ 1 file changed, 131 insertions(+), 48 deletions(-) diff --git a/python/findxtal.py b/python/findxtal.py index 800a037..ea2f8bb 100755 --- a/python/findxtal.py +++ b/python/findxtal.py @@ -81,7 +81,7 @@ def ffttest2(phi=45.,frq=4.2,amp=1.,n=256.): pass -def findGrid(image,numPeak=2): +def findGrid(image,numPeak=2,limFrq=None,debug=255): d2r=np.pi/180. #image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/grid_20180409_115332.png', flatten=True) # flatten=True gives a greyscale image s=image.shape @@ -97,21 +97,28 @@ def findGrid(image,numPeak=2): fft2 = np.fft.fft2(image) fft2=np.fft.fftshift(fft2) fa =abs(fft2) - plt.figure(num='log of fft: hamming wnd*image') - hi=plt.imshow(fa, interpolation="nearest",norm=mpl.colors.LogNorm(vmin=.1, vmax=fa.max())) - #plt.xlim(s[1]/2-50, s[1]/2+50);plt.ylim(s[0]/2-50, s[0]/2+50) + if debug&1: + plt.figure(num='log of fft: hamming wnd*image') + hi=plt.imshow(fa, interpolation="nearest",norm=mpl.colors.LogNorm(vmin=.1, vmax=fa.max())) + #plt.xlim(s[1]/2-50, s[1]/2+50);plt.ylim(s[0]/2-50, s[0]/2+50) ctr=np.array(image.shape,dtype=np.int16)/2 fa[ctr[0] - 1:ctr[0] + 2, ctr[1] - 1:ctr[1] + 2]=0 # set dc to 0 - hi.set_data(fa) - gen = np.zeros(fft2.shape) + # limit to maximal frequency + if limFrq is not None: + fa[:ctr[0]-limFrq,:]=0;fa[ctr[0]+limFrq:,:]=0 + fa[:,:ctr[1]-limFrq]=0;fa[:,ctr[1]+limFrq:]=0 x=np.arange(s[1])/float(s[1])*2.*np.pi y=np.arange(s[0])/float(s[0])*2.*np.pi - #x=np.linspace(0,2*np.pi,s[1],endpoint=False) - #y=np.linspace(0,2*np.pi,s[0],endpoint=False) + if debug&1: + hi.set_data(fa) + gen = np.zeros(fft2.shape) + #x=np.linspace(0,2*np.pi,s[1],endpoint=False) + #y=np.linspace(0,2*np.pi,s[0],endpoint=False) + xx, yy = np.meshgrid(x, y) + my=int(s[0]/2) mx=int(s[1]/2) - xx, yy = np.meshgrid(x, y) res=[] #list of tuples (freq_x,freq_y, phase) for i in range(numPeak): @@ -119,8 +126,9 @@ def findGrid(image,numPeak=2): maxAmpPos=np.array(divmod(maxAmpIdx,fa.shape[1]),dtype=np.int16) peakPos=maxAmpPos-ctr peak=fft2[maxAmpPos[0] - 1:maxAmpPos[0] + 2, maxAmpPos[1] - 1:maxAmpPos[1] + 2] - print(peakPos) - print(abs(peak)) + if debug&2: + print(peakPos) + print(abs(peak)) (vn, v0, vp)=np.log(fa[maxAmpPos[0], maxAmpPos[1] - 1:maxAmpPos[1] + 2]) #using log for interpolation is more precise freq_x=peakPos[1]+(vn-vp)/(2.*(vp+vn-2*v0)) (vn, v0, vp)=np.log(fa[maxAmpPos[0]-1:maxAmpPos[0]+2,maxAmpPos[1]]) @@ -137,9 +145,10 @@ def findGrid(image,numPeak=2): fy=peakPos[0]+iy amp=np.abs(v)/n; ang=np.angle(v) sumAmp+=amp - sumCos+=amp*np.cos(fx*x[mx] + fy*y[my] + ang) - sumSin+=amp*np.sin(fx*x[mx] + fy*y[my] + ang) - gen+=amp*np.cos(fx*xx + fy*yy + ang) + sumCos+=amp*np.cos(fx*x[mx+ix] + fy*y[my+iy] + ang) + sumSin+=amp*np.sin(fx*x[mx+ix] + fy*y[my+iy] + ang) + if debug&1: + gen+=amp*np.cos(fx*xx + fy*yy + ang) sumAmp*=2. #double because of conjugate part sumCos*=2. sumSin*=2. @@ -148,32 +157,82 @@ def findGrid(image,numPeak=2): if sumSin<0: w=-w phi_= freq_x*x[mx]+freq_y*y[my]-w phi_%=(np.pi*2) + #have main frequency positive and phase positive (for convinient) + if (freq_x<0 and abs(freq_x)> abs(freq_y)) or \ + (freq_y<0 and abs(freq_y)> abs(freq_x)): + freq_x = -freq_x; freq_y = -freq_y; phi_=-phi_ + if phi_<0: phi_+=2*np.pi + res.append((freq_x,freq_y,phi_)) fa[maxAmpPos[0]-1:maxAmpPos[0]+2,maxAmpPos[1]-1:maxAmpPos[1]+2]=0 # clear peak maxAmpPos_=2*ctr-maxAmpPos fa[maxAmpPos_[0]-1:maxAmpPos_[0]+2,maxAmpPos_[1]-1:maxAmpPos_[1]+2]=0 # clear conjugated peak - hi.set_data(fa) + if debug&1: + hi.set_data(fa) - gen*=2. # double because of conjugate part - for fx,fy,phase in res: - print('fx: %g fy: %g phase: %g deg'%(fx,fy,phase/d2r)) + if debug&1: + gen*=2. # double because of conjugate part + if debug&2: + for fx,fy,phase in res: + print('fx: %g fy: %g phase: %g deg'%(fx,fy,phase/d2r)) + + if debug&1: + plt.xlim(s[1]/2-50, s[1]/2+50) + plt.ylim(s[0]/2-50, s[0]/2+50) + plt.figure('image*wnd') + plt.imshow(image,interpolation="nearest") + plt.figure('reconstruct') + plt.imshow(gen,interpolation="nearest") + + plt.figure() + x=range(s[1]) + y=int(s[0]/2)-1 + plt.plot(x,image[y,:],'r') + plt.plot(x,gen[y,:],'g') + return res + +def plotGrid(grid,shape): + #x=np.linspace(0,2*np.pi,shape[1],endpoint=False) + #y=np.linspace(0,2*np.pi,shape[0],endpoint=False) + #for (freq_x, freq_y, phase) in grid: + #find points were: np.cos(freq_x*xx + freq_y*yy - phase) is max + #freq_x*xx + freq_y*yy - phase = 0 2pi, 4pi + # grid should have 2 entries + # -> 2 equations to solve + # points for: + # entry1 entry2 = (0 0), (0, 2*pi), (2*pi 0) + (fx0,fy0,p0)=grid[0] + (fx1,fy1,p1)=grid[1] + A=np.array([[fx0,fy0],[fx1,fy1]]) + A*=np.array([2*np.pi/shape[1],2*np.pi/shape[0]]) + Ai=np.asmatrix(A).I - plt.xlim(s[1]/2-50, s[1]/2+50) - plt.ylim(s[0]/2-50, s[0]/2+50) - plt.figure('image*wnd') - plt.imshow(image,interpolation="nearest") - plt.figure('reconstruct') - plt.imshow(gen,interpolation="nearest") + na=int(max(abs(fx0),abs(fy0)))+1 + nb=int(max(abs(fx1),abs(fy1)))+1 + p=np.ndarray((na*nb,2)) + #p=np.ndarray((3,2)) + #i=0 + #for x,y in ((0,0),(1,0),(0,1)): + # v = np.array([p0+x*2.*np.pi, p1+y*2.*np.pi]).reshape(-1, 1) + # p[i,:]=(Ai*v).reshape(-1) + # i+=1 - plt.figure() - x=range(s[1]) - y=int(s[0]/2)-1 - plt.plot(x,image[y,:],'r') - plt.plot(x,gen[y,:],'g') + i=0 + for b in range(nb): + for a in range(na): + v = np.array([p0+a*2.*np.pi, p1+b*2.*np.pi]).reshape(-1, 1) + p[i,:]=(Ai*v).reshape(-1) + i+=1 + + + + + #plt.plot([400,500],[400,500],'r+',markeredgewidth=2, markersize=10) + plt.plot(p[:,0],p[:,1],'r+',markeredgewidth=2, markersize=10) + plt.axis('image') pass - -def findObj(image,objSize=150,tol=0,viz=0): +def findObj(image,objSize=150,tol=0,debug=0): #objSiz is the rough diameter of the searched features in pixels #tol = tolerance in object size (not yet implemented) @@ -188,7 +247,7 @@ def findObj(image,objSize=150,tol=0,viz=0): #img2=ndi.filters.convolve1d(img2,box.reshape(-1),1) img2=ndi.filters.uniform_filter(np.float32(image),objSize*2) - if viz&32: + if debug&32: plt.imshow(image, interpolation="nearest", cmap='gray') plt.figure() plt.imshow(img2, interpolation="nearest", cmap='gray') @@ -196,14 +255,14 @@ def findObj(image,objSize=150,tol=0,viz=0): w=np.where(img2>image) img2[w]=image[w] img3=image-img2 - if viz&16: + if debug&16: plt.figure() plt.imshow(img3, interpolation="nearest", cmap='gray') #plt.show() l=int(objSize/30) if l>0: img4=ndi.binary_fill_holes(img3, structure=np.ones((l,l))) - if viz&8: + if debug&8: plt.figure() plt.imshow(img4, interpolation="nearest", cmap='gray') #plt.show() @@ -215,7 +274,7 @@ def findObj(image,objSize=150,tol=0,viz=0): #img5=ndi.binary_opening(img4, structure=np.ones((l,l))) #img5=ndi.binary_erosion(img4, structure=np.ones((l,l))) img5=ndi.binary_erosion(img4, iterations=l) - if viz&4: + if debug&4: plt.figure() plt.imshow(img5, interpolation="nearest", cmap='gray') #plt.show() @@ -237,7 +296,7 @@ def findObj(image,objSize=150,tol=0,viz=0): m=cv2.moments(contours[0]) lbl = ndi.label(img5) - if viz&2: + if debug&2: plt.figure() plt.imshow(lbl[0], interpolation="nearest") #plt.show() @@ -259,7 +318,7 @@ def findObj(image,objSize=150,tol=0,viz=0): pass #ctr2[i, :]=c[0,0] - if viz&1: + if debug&1: plt.figure() plt.imshow(image, interpolation="nearest", cmap='gray') plt.plot(ctr[:,1],ctr[:,0],'or',markeredgewidth=2, markersize=10) @@ -281,32 +340,56 @@ def genImg(shape,*args): image = dc+np.cos(freq_x*xx + freq_y*yy - phase) else: image += np.cos(freq_x * xx + freq_y * yy - phase) - - plt.imshow(image, interpolation="nearest") return image if __name__ == '__main__': #plt.ion() #ffttest() - image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/grid_20180409_115332.png') - #image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/grid_20180409_115332_45deg.png') + + image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/honeycomb.png') image=-image - #image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/honeycomb.png') + grid=findGrid(image,numPeak=2,limFrq=25,debug=2) + plt.imshow(image, interpolation="nearest", cmap='gray') + plotGrid(grid,image.shape) + plt.show() + + image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/grid_20180409_115332.png') + image=-image + grid=findGrid(image,numPeak=2,debug=2) + plt.imshow(image, interpolation="nearest", cmap='gray') + plotGrid(grid,image.shape) + plt.show() + + image = ndimage.imread('/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/grid_20180409_115332_45deg.png') + image=-image + grid=findGrid(image,numPeak=2,debug=2) + plt.imshow(image, interpolation="nearest", cmap='gray') + plotGrid(grid,image.shape) + plt.show() + + objPos=findObj(image,debug=1) + plt.show() + d2r=np.pi/180. #for phi in np.arange(0.,180.,10.): # image = genImg((600, 800), (4.5, .2, phi*d2r)) # plt.show() #image=genImg((600,800),(4.,1.0,10.*d2r)) - #image=genImg((600,800),(4.5,.2,20.*d2r)) + #image=genImg((600,800),(4.5,-3.2,70.*d2r)) + #image=genImg((600,800),(-4.5,3.2,290.*d2r)) #same image #findGrid(image,numPeak=1) + #findGrid(image,numPeak=1,debug=2) - #image=genImg((600,800),(9.5,.2,20.*d2r),(.4,5.2,60.*d2r)) - #image=genImg((600,800),(3.5,0.,0.*d2r),(0,5.,0.*d2r)) - findGrid(image,numPeak=2) - + #for v in np.arange(0,2,.3): + # image=genImg((600,800),(8,.2+v,40.*d2r),(.4,5.2,30.*d2r)) + # grid=findGrid(image,numPeak=2,debug=2) + # plt.figure(1);plt.cla() + # plt.imshow(image, interpolation="nearest", cmap='gray') + # plotGrid(grid,image.shape) + # plt.show() #image=genImg((600,800),(9.5,.2,0),(.4,5.2,0),(4,8,0)) #findObj(image,viz=1) #findObj(image,viz=255) #print(findObj(image)) - plt.show() +