513 lines
17 KiB
Python
Executable File
513 lines
17 KiB
Python
Executable File
#!/usr/bin/env python
|
|
#*-----------------------------------------------------------------------*
|
|
#| |
|
|
#| Copyright (c) 2018 by Paul Scherrer Institute (http://www.psi.ch) |
|
|
#| |
|
|
#| Author Thierry Zamofing (thierry.zamofing@psi.ch) |
|
|
#*-----------------------------------------------------------------------*
|
|
'''
|
|
implements an image alalyser for ESB MX
|
|
'''
|
|
import os
|
|
from scipy import fftpack, ndimage
|
|
import scipy.ndimage as ndi
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib as mpl
|
|
import numpy as np
|
|
import imgStack
|
|
|
|
def findGrid(image,numPeak=2,minFrq=2,maxFrq=None,debug=0):
|
|
d2r=np.pi/180.
|
|
s=image.shape
|
|
w1=np.hamming(s[0]).reshape((-1,1))
|
|
w2=np.hamming(s[1]).reshape((1,-1))
|
|
wnd=w1*w2
|
|
#plt.figure(num='hamming window')
|
|
#plt.imshow(wnd, interpolation="nearest")
|
|
image=wnd*image
|
|
if debug&1:
|
|
plt.figure(num='hamming window*img')
|
|
plt.imshow(image, interpolation="nearest")
|
|
|
|
fft2 = np.fft.fft2(image)
|
|
fft2=np.fft.fftshift(fft2)
|
|
fa =abs(fft2)
|
|
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]-minFrq+1:ctr[0]+minFrq, ctr[1]-minFrq+1:ctr[1]+minFrq]=0 # set dc to 0
|
|
# limit to maximal frequency
|
|
if maxFrq is not None:
|
|
fa[:ctr[0]-maxFrq,:]=0;fa[ctr[0]+maxFrq:,:]=0
|
|
fa[:,:ctr[1]-maxFrq]=0;fa[:,ctr[1]+maxFrq:]=0
|
|
x=np.arange(s[1])/float(s[1])*2.*np.pi
|
|
y=np.arange(s[0])/float(s[0])*2.*np.pi
|
|
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)
|
|
|
|
res=[] #list of tuples (freq_x,freq_y, phase)
|
|
for i in range(numPeak):
|
|
maxAmpIdx=fa.argmax()
|
|
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]
|
|
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]])
|
|
freq_y=peakPos[0]+(vn-vp)/(2.*(vp+vn-2*v0))
|
|
#calculate phase
|
|
sumCos=0.
|
|
sumSin=0.
|
|
sumAmp=0.
|
|
n=np.prod(s)
|
|
for iy in (-1,0,1):#(0,):
|
|
for ix in (-1,0,1):
|
|
v=peak[iy+1,ix+1]
|
|
fx=peakPos[1]+ix
|
|
fy=peakPos[0]+iy
|
|
amp=np.abs(v)/n; ang=np.angle(v)
|
|
sumAmp+=amp
|
|
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.
|
|
#sumAmp=np.abs(peak).sum()/n*2 #double because of conjugate part
|
|
w=np.arccos(sumCos/sumAmp)
|
|
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
|
|
if debug&1:
|
|
hi.set_data(fa)
|
|
|
|
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
|
|
|
|
|
|
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
|
|
|
|
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 ShowImage(img,title=None,cmap='gray',vmin=None, vmax=None):
|
|
plt.figure(title)
|
|
plt.imshow(img, interpolation="nearest", cmap=cmap,vmin=vmin,vmax=vmax) # ,vmin=m-3*s, vmax=m+3*s)
|
|
plt.colorbar()
|
|
|
|
def imgEqualize(img,num_bins=256):
|
|
# get image histogram
|
|
hist, bins = np.histogram(img.flatten(), num_bins, normed=True)
|
|
cdf = hist.cumsum() # cumulative distribution function
|
|
cdf = 255 * cdf / cdf[-1] # normalize
|
|
|
|
# use linear interpolation of cdf to find new pixel values
|
|
imgEqu=np.interp(img.flatten(), bins[:-1], cdf)
|
|
imgEqu=np.uint8(imgEqu.reshape(img.shape))
|
|
|
|
plt.figure();plt.plot(hist)
|
|
h2,b2=np.histogram(imgEqu.flatten(), num_bins, normed=True)
|
|
plt.figure();plt.plot(h2)
|
|
return imgEqu.reshape(img.shape)#, cdf
|
|
|
|
|
|
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)
|
|
#tolShape = roudness tolerance in object roundness (not yet implemented)
|
|
|
|
from scipy.signal import convolve2d
|
|
plt.ion()
|
|
image=image[500:2500,1000:2500]
|
|
s=image.shape
|
|
f=np.array([0.9595264, 0.9600567, 0.9608751, 0.9620137, 0.9634765, 0.9652363, 0.9672352, 0.9693891, 0.9715959, 0.9737464, 0.9757344, 0.9774676, 0.9788761, 0.9799176, 0.9805792, 0.9808776, 0.9808528, 0.9805624, 0.9800734, 0.9794550])
|
|
f=f*1000
|
|
f=f-f.mean()
|
|
|
|
m=image.mean();s=image.std()
|
|
ShowImage(image,vmin=-.85,vmax=-.99)
|
|
ShowImage(image,vmin=m-3*s,vmax=m+3*s)
|
|
ieq=imgEqualize(image)
|
|
ShowImage(ieq)
|
|
m=image.mean();s=image.std()
|
|
i2=(image-m)*(256./(3*s))+128.
|
|
i2[i2>255.]=255.
|
|
i2[i2<0.]=0.
|
|
i2=np.uint8(i2)
|
|
ShowImage(i2)
|
|
image=ndi.filters.gaussian_filter1d(image,sigma=5./3,truncate=3.,axis=0)
|
|
image=ndi.filters.gaussian_filter1d(image,sigma=5./3,truncate=3.,axis=1)
|
|
ShowImage(image,vmin=-.85,vmax=-.99)
|
|
|
|
img2=ndi.filters.convolve1d(image,f,0)
|
|
ShowImage(img2,vmin=-3,vmax=3)
|
|
|
|
img3=ndi.filters.convolve1d(image,f,1)
|
|
ShowImage(img3,vmin=-3,vmax=3)
|
|
img4=np.maximum(abs(img2),abs(img3))
|
|
ShowImage(img4,vmin=0,vmax=3)
|
|
|
|
m=image.mean();s=image.std()
|
|
img5=np.zeros(image.shape)
|
|
w=np.where(img4>1.5)
|
|
img5[w]=1
|
|
ShowImage(img5)
|
|
#imgStack.Run([image, img2,img3,img4,img5])
|
|
#w=np.where(img2>image*1.01)
|
|
#img2[:]=0
|
|
#img2[w]=1
|
|
l=int(objSize/30);l=5
|
|
img6=ndi.binary_fill_holes(img5, structure=np.ones((l,l)))
|
|
|
|
#img6 = ndi.binary_dilation(img5, iterations=2)
|
|
#img6 = ndi.binary_erosion(img5, iterations=2)
|
|
ShowImage(img6)
|
|
|
|
l=int(objSize/5)#=int(objSize/10)
|
|
if l>=3:
|
|
#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 debug&4:
|
|
plt.figure()
|
|
plt.imshow(img5, interpolation="nearest", cmap='gray')
|
|
#plt.show()
|
|
else:
|
|
img5=img4
|
|
|
|
import cv2
|
|
#cvi1=np.zeros(shape=img5.shape+(3,), dtype=np.uint8)
|
|
#cvi1[:,:,0]=image
|
|
#image*np.array([1,1,1]).reshape(-1,1,1)
|
|
s=image.shape+(1,)
|
|
cvi1=image.reshape(s)*np.ones((1,1,3),dtype=np.uint8)
|
|
#cvi1=np.ones((3,1,1),dtype=np.uint8)image
|
|
contours, hierarchy=cv2.findContours(np.uint8(img5),cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
|
|
#contours,hierarchy=cv2.findContours(np.uint8(img5),1,2)
|
|
cv2.drawContours(cvi1, contours, -1, (0,255,0), 3)
|
|
plt.figure()
|
|
plt.imshow(cvi1 , interpolation="nearest", cmap='gray')
|
|
m=cv2.moments(contours[0])
|
|
|
|
lbl = ndi.label(img5)
|
|
if debug&2:
|
|
plt.figure()
|
|
plt.imshow(lbl[0], interpolation="nearest")
|
|
#plt.show()
|
|
ctr=ndi.measurements.center_of_mass(image, lbl[0],range(lbl[1]))
|
|
ctr=np.array(ctr)
|
|
ctr2=np.ndarray(shape=(len(contours),2),dtype=np.uint16)
|
|
i=0
|
|
for c in contours:
|
|
m=cv2.moments(c)
|
|
try:
|
|
m00=m['m00']
|
|
m10=m['m10']
|
|
m01=m['m01']
|
|
#print m00
|
|
if m00>1000 and m00<7000:
|
|
ctr2[i,:]=(m10/m00,m01/m00)
|
|
i+=1
|
|
except ZeroDivisionError:
|
|
pass
|
|
#ctr2[i, :]=c[0,0]
|
|
|
|
if debug&1:
|
|
plt.figure()
|
|
plt.imshow(image, interpolation="nearest", cmap='gray')
|
|
plt.plot(ctr[:,1],ctr[:,0],'or',markeredgewidth=2, markersize=10)
|
|
plt.plot(ctr2[:i,0],ctr2[:i,1],'+y',markeredgewidth=2, markersize=10)
|
|
plt.show()
|
|
return ctr
|
|
|
|
def genImg(shape,*args):
|
|
'''args is a list of tuples (freq_x,freq_y, phase) multiple args can be added'''
|
|
image=np.ndarray(shape)#,dtype=np.uint8)
|
|
x=np.linspace(0,2*np.pi,shape[1],endpoint=False)
|
|
y=np.linspace(0,2*np.pi,shape[0],endpoint=False)
|
|
#xx, yy = np.meshgrid(x, y, sparse=True)
|
|
xx, yy = np.meshgrid(x, y)
|
|
dc=0
|
|
for i,f in enumerate(args):
|
|
(freq_x, freq_y, phase)=f
|
|
if i==0:
|
|
image = dc+np.cos(freq_x*xx + freq_y*yy - phase)
|
|
else:
|
|
image += np.cos(freq_x * xx + freq_y * yy - phase)
|
|
return image
|
|
|
|
def phase_retrieval_intensity_tranport(image, mu, delta, I_in, M, pixel_size, R2):
|
|
'''calculates the projected thickness as in eq. 12 of Paganin et al, 2002, J. Microscopy,
|
|
from an image, the absoption coefficient mu (um-1), the real part of the deviation of the refractive index
|
|
from unity delta, the uniform intensity of the incident radiation I_in (ph/s/um2), the magnification of the
|
|
image from the point source illumination M, the pixel size (um), the propagation distance R2 (um).'''
|
|
F = np.fft.rfft2(image)
|
|
k_x = 1/(pixel_size*image.shape[0]) # is there a 2pi factor?
|
|
k_y = 1/(pixel_size*image.shape[1]) # is there a 2pi factor?
|
|
#print k_x,k_y
|
|
k_array = np.zeros((F.shape[0],F.shape[1]))
|
|
A = np.zeros((F.shape[0],F.shape[1]), dtype=np.complex128)
|
|
for i_k in range(k_array.shape[0]):
|
|
for j_k in range(k_array.shape[1]):
|
|
k_array[i_k,j_k] = np.sqrt((i_k * k_x)**2 + (j_k * k_y)**2)
|
|
A[i_k,j_k] = mu * F[i_k,j_k] / (I_in * (R2 * delta * k_array[i_k,j_k]**2 / M + mu))
|
|
Fm1 = np.fft.irfft2(A)
|
|
T = np.multiply(-1/mu, np.log(Fm1))
|
|
return T
|
|
|
|
if __name__ == '__main__':
|
|
def testfftLoop():
|
|
plt.ioff()
|
|
for phi in np.arange(0.,180.,10.):
|
|
testfft(phi=phi)
|
|
plt.show()
|
|
pass
|
|
|
|
def testfft(phi=45.,frq=4.2,amp=1.,n=256.):
|
|
#find the main frequency and phase in 1-D
|
|
plt.figure(1)
|
|
d2r=np.pi/180.
|
|
x=np.arange(n)
|
|
phi=phi*d2r
|
|
y=amp*np.cos(frq*x/n*2.*np.pi-phi)
|
|
plt.plot(x,y,'y')
|
|
#y[np.where(y<-.5)]=0
|
|
#y*=np.hamming(n)
|
|
w=np.hanning(n)
|
|
plt.plot(x,w,'y')
|
|
y*=w
|
|
#y*=1.-np.cos(x/(n-1.)*2.*np.pi)
|
|
|
|
#y=[1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1]
|
|
plt.stem(x,y)
|
|
|
|
fy=np.fft.fft(y)
|
|
fya=np.abs(fy)
|
|
plt.figure(2)
|
|
plt.subplot(211)
|
|
plt.stem(x,fya)
|
|
plt.subplot(212)
|
|
plt.stem(x,np.angle(fy)/d2r)
|
|
print(np.angle(fy[frq])/d2r)
|
|
|
|
fya[0]=0
|
|
i=(fya.reshape(2,-1)[0,:]).argmax()
|
|
(vn,v0,vp)=fya[i-1:i+2]
|
|
frq_=i+(vn-vp)/(2.*(vp+vn-2*v0))
|
|
print('freq: %g phase %g %g %g'%((frq_,)+tuple(np.angle(fy[i-1:i+2])/d2r)))
|
|
|
|
#PHASE CALCULATION
|
|
plt.figure(1)
|
|
y=np.zeros(x.shape)
|
|
z=np.zeros(x.shape)
|
|
for ii in (i,i-1,i+1):
|
|
y+=np.abs(fy[ii])/n*np.cos(ii*x/n*2.*np.pi+np.angle(fy[ii]))
|
|
z+=np.abs(fy[ii])/n*np.sin(ii*x/n*2.*np.pi+np.angle(fy[ii]))
|
|
y*=2 #double because of conjugate part
|
|
z*=2 #double because of conjugate part
|
|
plt.plot(x,y,'r')
|
|
amp_=fya[i-1:i+2].sum()/n*2.
|
|
t=int(n/2)-1
|
|
#->phase: find maximum or where the sin is 0
|
|
w=np.arccos(y[t]/amp_)
|
|
if(z[t]<0): w=-w
|
|
print('amplitude %g, value at middle (%d) cos->%g sin->%g -> acos %g deg'%(amp_,t,y[t],z[t],w/d2r))
|
|
#rot angle at middle x=127 =frq_*t/n*2.*np.pi-phi_=w '%(amp_,t,y[t],phi_/d2r))
|
|
phi_=frq_*t/n*2.*np.pi-w
|
|
y=amp_*np.cos(frq_*x/n*2.*np.pi-phi_)
|
|
print('y[%d] %g'%(t,y[t]))
|
|
|
|
plt.plot(x,y,'g')
|
|
pass
|
|
|
|
def testFindGrid():
|
|
plt.ioff()
|
|
imggrid=(
|
|
('grid_20180409_115332_45deg.png', 2,None),
|
|
('grid_20180409_115332.png', 2,None),
|
|
('honeycomb.png', 2,25),
|
|
('MS01_20180411_110544.png', 10,None),
|
|
('MS02_20180411_120354.png', 10,None),
|
|
('MS03_20180411_135524.png', 10,None),
|
|
('MS04_20180411_143045.png', 10,None),
|
|
('MS04_20180411_144239.png', 10,None),
|
|
)
|
|
for (file,minFrq,maxFrq) in imggrid:
|
|
image = ndimage.imread(os.path.join(basePath,file))
|
|
image=-image
|
|
grid=findGrid(image,minFrq=minFrq,maxFrq=maxFrq,numPeak=2)
|
|
plt.figure('grid');plt.imshow(image, interpolation="nearest", cmap='gray')
|
|
plotGrid(grid, image.shape)
|
|
plt.show()
|
|
|
|
def testFindObj():
|
|
plt.ioff()
|
|
imggrid=(
|
|
('grid_20180409_115332_45deg.png', 2,None),
|
|
)
|
|
for (file,p0,p1) in imggrid:
|
|
image = ndimage.imread(os.path.join(basePath,file))
|
|
image=-image
|
|
objPos=findObj(image,debug=0)
|
|
plt.figure('findObj');plt.imshow(image, interpolation="nearest", cmap='gray')
|
|
plt.plot(objPos[:,1],objPos[:,0],'r+',markeredgewidth=2, markersize=10)
|
|
plt.axis('image')
|
|
plt.show()
|
|
|
|
def testPhaseEnhance():
|
|
|
|
mu=0.001 # um-1, optimal value found is 0.001
|
|
phase_shift=-0.4 # rad/um, optimal value found is -1
|
|
E=18 # keV
|
|
wavelength=12.4/E/10000 # um
|
|
delta=-wavelength/(2*np.pi)*phase_shift
|
|
delta=4.38560287631e-06
|
|
M=1
|
|
R2=19*1e4 # um
|
|
pixel_size=0.325 # um
|
|
flux=1e12 # ph/s
|
|
for fn in ('lyso1_scan_18keV_190mm_MosaicJ_noblend_crop.tif',
|
|
'PepT2_scan_18keV_190mm_MosaicJ_noblend_crop_FFT40.tif',
|
|
'SiN_lysoS1_scan_18keV_190mm_MosaicJ_noblend_crop_FFT40.tif',
|
|
'SiN_PepT2_scan_18keV_190mm_MosaicJ_noblend_crop.tif',):
|
|
img = ndimage.imread(os.path.join(basePath,fn))
|
|
FOV_pix = img.shape
|
|
FOV_um = (FOV_pix[0] * pixel_size, FOV_pix[1] * pixel_size)
|
|
I_in = flux / (FOV_um[0] * FOV_um[1]) # ph/s/um2
|
|
phase_image = phase_retrieval_intensity_tranport(img, mu, delta, I_in, M, pixel_size, R2)
|
|
phase_image_flipped = phase_retrieval_intensity_tranport(np.fliplr(img), mu, delta, I_in, M, pixel_size, R2)
|
|
added_image = np.add(phase_image,np.fliplr(phase_image_flipped))
|
|
|
|
|
|
m=img.mean();s=img.std()
|
|
ShowImage(img, title=fn+' raw', vmin=m-3*s, vmax=m+3*s)
|
|
m=added_image.mean();s=added_image.std()
|
|
ShowImage(added_image, title=fn+' phase retrival', vmin=m-3*s, vmax=m+3*s)
|
|
plt.show()
|
|
|
|
basePath='/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/images/'
|
|
#testfftLoop()
|
|
#testFindGrid()
|
|
testFindObj()
|
|
#testPhaseEnhance()
|
|
exit(0)
|
|
|
|
plt.ion()
|
|
#image = ndimage.imread(os.path.join(basePath, 'lyso1_scan_18keV_190mm_MosaicJ_noblend_crop.tif'))
|
|
image = ndimage.imread(os.path.join(basePath, 'lyso1_scan_18keV_190mm_MosaicJ_noblend_crop.tif'))
|
|
#image = ndimage.imread(os.path.join(basePath, 'grid_20180409_115332.png'))
|
|
image = -image
|
|
#plt.figure('input');
|
|
#plt.imshow(image, interpolation="nearest", cmap='gray')
|
|
#sbl=ndi.sobel(image)
|
|
#plt.figure('sobel');
|
|
#plt.imshow(sbl, interpolation="nearest", cmap='gray')
|
|
objPos = findObj(image, objSize=50,debug=255)
|
|
plt.figure('findObj');
|
|
plt.imshow(image, interpolation="nearest", cmap='gray')
|
|
plt.plot(objPos[:, 1], objPos[:, 0], 'r+', markeredgewidth=2, markersize=10)
|
|
plt.axis('image')
|
|
plt.show()
|
|
|
|
#d2r=np.pi/180.
|
|
#image=genImg((600,800),(4.,1.0,10.*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)
|
|
|
|
#for v in np.arange(0,2,.3):
|
|
# image=genImg((600,800),(8,.2+v,40.*d2r),(.4,5.2,30.*d2r))
|
|
# #image = genImg((600, 800), (.4,5.2,30.))
|
|
# plt.figure(10);plt.cla()
|
|
# plt.imshow(image)
|
|
# grid=findGrid(image,numPeak=2,debug=255)
|
|
# 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))
|
|
|