"This notebook corresponds to version {{ version }} of the pipeline tool: https://github.com/NSLS-II/pipelines"
In [1]:
from databroker import DataBroker as db, get_images, get_table, get_events
from filestore.api import register_handler, deregister_handler
from filestore.retrieve import _h_registry, _HANDLER_CACHE
In [2]:
#hdr = db[{{ uid }}]
In [3]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import time
from ipywidgets import interact
In [4]:
cd /XF11ID/analysis/Analysis_Pipelines/Develop/
In [6]:
%run develop.py
%run Time_Correlation_Functions.py
%run speckle.py
In [7]:
%matplotlib notebook
#%matplotlib inline
In [8]:
BlueScan = True
DirectAcq = False
detector = 'eiger_4M_cam_img_image_lightfield' #for 4M
In [9]:
if BlueScan:
uid = '54614d43'
#uid = '95782687'
uid = '95782687'
uid= 'ff9f20c0'
uid='71720966'
uid='1663d34a'
uid = 'f505e052-3baa-47d4-bdc4-61c2eb1bcc7a' #sid= 551, 1%PEG,
uid='ee6975a1-9161' #1% wt PEG
else:
uid = '/XF11ID/data/2015/11/23/d01ab510-3cf3-4719-bee3_795_master.h5'
In [12]:
if BlueScan:
hdr = db[uid]
ev, = get_events( hdr, [detector] )
imgs = ev['data'][detector]
else:
imgs = Images(uid)
print (imgs)
Nimg=len(imgs)
In [13]:
if BlueScan:
from datetime import datetime
dt = datetime.fromtimestamp(hdr['start'].time)
path ='/XF11ID/analysis' + '/%s/%s/%s/' % (dt.year, dt.month, dt.day)
else:
path ='/XF11ID/analysis/2015/11/23/'
path
Out[13]:
In [14]:
imgs.md
Out[14]:
In [15]:
# The physical size of the pixels
dpix = imgs.md['x_pixel_size'] * 1000.
lambda_ = imgs.md['incident_wavelength'] # wavelegth of the X-rays in Angstroms
Ldet = 4812. # detector to sample distance (mm)
exposuretime= imgs.md['count_time']
acquisition_period = imgs.md['frame_time']
# deadtime= 0 # 60e-6
# timeperframe = exposuretime + deadtime
timeperframe = acquisition_period
timeperframe, exposuretime
Out[15]:
In [16]:
mask = np.load( path + str(uid)+ "_mask.npy")
In [17]:
maskr = mask[::-1,:]
In [18]:
fig, ax = plt.subplots()
im=ax.imshow(maskr, origin='lower' ,vmin=0,vmax=1,cmap='viridis')
fig.colorbar(im)
plt.show()
In [19]:
def view_image(i):
fig, ax = plt.subplots()
ax.imshow(imgs[i]*mask, interpolation='nearest', cmap='viridis',
origin='lower', norm= LogNorm(vmin=0.001, vmax=1e1) )
ax.set_title("Browse the Image Stack")
plt.show()
In [20]:
#interact(view_image, i=(0, Nimg-1))
In [21]:
def view_image(sleeps=1, ims=0, ime = 1):
fig, ax = plt.subplots()
for i in range( ims, ime ):
im=ax.imshow(imgs[i]*mask, interpolation='nearest', cmap='viridis',
origin='lower', norm= LogNorm( vmin=0.01, vmax=10 ) )
ax.set_title("images_%s"%i)
time.sleep( sleeps )
plt.draw()
#fig.colorbar(im)
#view_image(.2, 0, 2)
In [22]:
kymo_sum = np.load( path + str(uid)+"_kymo_sum.npy" )
In [23]:
bad_frames = np.where( kymo_sum > 1e5)[0]
bad_frames
Out[23]:
In [24]:
fig, axes = plt.subplots( )
axes.plot( kymo_sum, '-go' )
ax.set_ylabel('Intensity')
ax.set_xlabel('Frame')
ax.set_title('Kymograph_sum')
plt.show()
In [25]:
avg_img = np.load( path + str(uid)+"_avg_img.npy" )
avg_imgm = avg_img * mask
In [26]:
avg_imgr = avg_img[::-1,:]
avg_imgmr = avg_imgm[::-1,:]
In [27]:
fig, ax = plt.subplots()
im = ax.imshow(avg_imgmr, cmap='viridis',origin='lower',
norm= LogNorm(vmin=0.001, vmax=1e1))
ax.set_title("Masked Averaged Image_Reversed")
fig.colorbar(im)
plt.show()
In [28]:
imgs.md['beam_center_x'], imgs.md['beam_center_y']
Out[28]:
In [29]:
#center = (imgs.md['beam_center_x'], imgs.md['beam_center_y'])
center = [ 2167 - 336, 849] #for not reversed
center = [ 336, 849] #for reversed
center = [ 2167- 1830, 846]
center
Out[29]:
In [31]:
fig, ax = plt.subplots()
im = ax.imshow(avg_imgr, cmap='viridis',origin='lower', norm= LogNorm(vmin=0.001, vmax=1e1))
radius = 54
circle=plt.Circle( [center[1], center[0]], radius, color='b', alpha=1.0, lw=2, edgecolor='r',fill=False)
plt.gcf().gca().add_artist(circle)
ax.set_title("Masked Averaged Image_Reversed")
fig.colorbar(im)
rwidth = 100
x1,x2 = [center[1] - rwidth, center[1] + rwidth]
y1,y2 = [center[0] - rwidth, center[0] + rwidth]
ax.set_xlim( [x1,x2])
ax.set_ylim( [y1,y2])
plt.show()
In [32]:
bin_centers, ring_averages= circular_average(avg_imgr, center, pixel_size=(dpix, dpix), mask= maskr)
# convert to q (reciprocal space)
two_theta = utils.radius_to_twotheta(Ldet, bin_centers)
q_val = utils.twotheta_to_q(two_theta, lambda_)
In [33]:
fig,axes = plt.subplots(figsize=(8, 6))
axes.semilogy(q_val, ring_averages, '-o')
axes.set_title('Circular Average')
axes.set_ylabel('Ring Avearge')
axes.set_xlabel('Q ('r'$\AA^{-1}$ )')
axes.set_xlim(0.001, 0.02)
axes.set_ylim(0.001, 10.0)
plt.show()
In [34]:
fig,axes = plt.subplots(figsize=(8, 6))
axes.semilogy(bin_centers/dpix, ring_averages, '-o')
axes.set_title('Circular Average')
axes.set_ylabel('Ring Avearge')
axes.set_xlabel('Bin Centers, (pixel)')
axes.set_xlim(30, 250)
axes.set_ylim(0.001, 10.0)
plt.show()
In [35]:
inner_radius = 58 # radius of the first ring
width = 2 # width of each ring
spacing = (166 - 58)/9 - 2 # spacing between rings
num_rings = 6 # number of rings
# find the edges of the required rings
edges = roi.ring_edges(inner_radius, width, spacing, num_rings)
edges
Out[35]:
In [36]:
two_theta = utils.radius_to_twotheta(Ldet, edges*dpix)
q_ring_val = utils.twotheta_to_q(two_theta, lambda_)
q_ring_center = np.average(q_ring_val, axis=1)
q_ring_center
Out[36]:
In [37]:
rings = roi.rings(edges, center, avg_imgmr.shape)
ring_mask = rings*maskr
In [38]:
ring_mask
Out[38]:
In [39]:
qind, pixelist = roi.extract_label_indices( ring_mask )
noqs = len( np.unique(qind) )
nopr = np.bincount(qind, minlength=(noqs+1))[1:]
In [40]:
qind
Out[40]:
In [41]:
nopr
Out[41]:
In [42]:
pixel = roi.roi_pixel_values(avg_imgmr, ring_mask, [2] )
fig,ax=plt.subplots()
ax.plot( pixel[0][0] ,'bo', ls='-' )
Out[42]:
In [43]:
fig, axes = plt.subplots( figsize=(8, 6))
#axes.semilogy(q_val, ring_averages, '-o')
axes.plot(q_val, ring_averages, '-o')
axes.set_title('Circular Average with the Q ring values')
axes.set_ylabel('Ring Avearge')
axes.set_xlabel('Bin Centers 'r'$\AA^{-1}$')
axes.set_xlim(0.00, 0.02)
axes.set_ylim(0, 6)
for i in range(num_rings):
axes.axvline(q_ring_center[i])
plt.show()
In [47]:
# plot the figure
fig, axes = plt.subplots(figsize=(8,8))
axes.set_title("Labeled Array on Averaged Data")
im,im_label = show_label_array_on_image(axes, avg_imgmr, ring_mask, imshow_cmap='viridis',
cmap='Paired',
vmin=0.01, vmax=5, origin="lower")
#rwidth = 200
#x1,x2 = [center[1] - rwidth, center[1] + rwidth]
#y1,y2 = [center[0] - rwidth, center[0] + rwidth]
#axes.set_xlim( [x1,x2])
#axes.set_ylim( [y1,y2])
#fig.colorbar(im)
rwidth = 200
x1,x2 = [center[1] - rwidth, center[1] + rwidth]
y1,y2 = [center[0] - rwidth, center[0] + rwidth]
axes.set_xlim( [x1,x2])
axes.set_ylim( [y1,y2])
fig.colorbar(im_label)
plt.show()
In [48]:
imgs_ =imgs
imgsr = Reverse_Coordinate(imgs_, mask)
In [49]:
max_inten_ring =2
In [44]:
#plt.close('all')
In [49]:
#kymo = roi.kymograph(imgsr, ring_mask, num = max_inten_ring)
In [52]:
t0 = time.time()
data_pixel = Get_Pixel_Array( imgsr, pixelist).get_data()
run_time(t0)
In [53]:
data_pixel.shape
Out[53]:
In [54]:
pixelist.shape
Out[54]:
In [50]:
#np.save( path + 'uid_%s_data_pixel'%uid, data_pixel )
#data_pixel = np.load( path + 'uid_%s_data_pixel.npy'%(uid) )
In [51]:
pixelist_qi = np.where( qind == max_inten_ring)[0]
data_pixel_qi = data_pixel[:,pixelist_qi]
In [52]:
fig, ax = plt.subplots(figsize=(8,6))
ax.set_ylabel('Pixel')
ax.set_xlabel('Frame')
ax.set_title('Kymograph')
im = ax.imshow(data_pixel_qi.T, cmap='viridis', vmax=5.0)
fig.colorbar( im )
ax.set_aspect(30.)
plt.show()
In [53]:
mean_inten = get_mean_intensity( data_pixel, qind)
In [54]:
times = np.arange( mean_inten[1].shape[0] ) #*timeperframe # get the time for each frame
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_title("Mean intensity of each ring")
for i in range(num_rings):
ax.plot(times, mean_inten[i+1], '--o', label="Ring "+str(i+1))
ax.set_xlabel("Frame")
ax.set_ylabel("Mean Intensity")
ax.legend(loc='best')
plt.show()
Note : Enter the number of levels and number of buffers for Muliti tau one time correlation number of buffers has to be even. More details in https://github.com/scikit-xray/scikit-xray/blob/master/skbeam/core/correlation.py
In [55]:
nopr
Out[55]:
In [56]:
edges
Out[56]:
In [57]:
num_buf =8
In [58]:
ring_mask = np.array( ring_mask, dtype=int)
In [59]:
good_start = 0
good_end = Nimg #all the images
#good_start = 7600
#good_end = 7650
imgs_ =imgs[good_start: good_end-1]
imgsr = Reverse_Coordinate(imgs_, mask)
In [67]:
g2_t1, lag_steps_t1 = autocor_one_time( num_buf, ring_mask, imgsr, num_lev=None,
bad_images=None, threshold= 65500 )
In [68]:
g2, lag_steps = g2_t1, lag_steps_t1
lags = lag_steps*timeperframe
In [69]:
tg2 = np.hstack( [ lags.reshape( len(lags),1), g2] )
#np.save( path + 'g2_%s-%s--%s'%(uid,good_start, good_end), tg2)
np.savetxt( path + 'g2_%s-%s--%s.txt'%(uid,good_start, good_end), tg2)
In [62]:
#tg2 = np.loadtxt( path + 'g2_%s-%s--%s.txt'%(uid,good_start, good_end) )
#lags=tg2[:,0]
#g2=tg2[:,1:]
In [103]:
#plt.close('all')
In [63]:
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0:
sy = int(num_rings/sx)
else:
sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)
plt.axis('off')
for sn in range(num_rings):
ax = fig.add_subplot(sx,sy,sn+1 )
ax.set_ylabel("g2")
ax.set_title(" Q= " + '%.5f '%(q_ring_center[sn]) + r'$\AA^{-1}$')
y=g2[:, sn]
ax.semilogx(lags, y, '-o', markersize=6)
#ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
ax.set_ylim([ min(y) , max(y[1:]) ])
ax.set_xlim([ min(lags)+0*1e-6, max(lags)])
plt.show()
fig.tight_layout()
In [64]:
from lmfit import Model
mod = Model(corr.auto_corr_scat_factor)
In [65]:
rate = []
sx = int( round (np.sqrt(num_rings)) )
if num_rings%sx==0:
sy = int(num_rings/sx)
else:
sy = int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid, fontsize=20, y =1.08)
plt.axis('off')
for sn in range(num_rings):
ax = fig.add_subplot(sx, sy, sn+1 )
y=g2[1:, sn]
result1 = mod.fit(y, lags=lags[1:], beta=.1,
relaxation_rate =.5, baseline=1.0)
rate.append(result1.best_values['relaxation_rate'])
ax.semilogx(lags[1:], y, 'ro')
ax.semilogx(lags[1:], result1.best_fit, '-b')
ax.set_title(" Q= " + '%.5f '%(q_ring_center[sn]) + r'$\AA^{-1}$')
ax.set_ylim([min(y)*.95, max(y[1:]) *1.05])
txts = r'$\gamma$' + r'$ = %.3f$'%(rate[sn]) + r'$ s^{-1}$'
ax.text(x =0.015, y=.55, s=txts, fontsize=14, transform=ax.transAxes)
fig.tight_layout()
In [66]:
fig, ax=plt.subplots()
ax.plot(q_ring_center**2, rate, 'ro', ls='--')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()
In [67]:
D0 = np.polyfit(q_ring_center**2, rate, 1)
gmfit = np.poly1d(D0)
print ('The fitted diffusion coefficient D0 is: %.2E A^2S-1'%D0[0])
In [68]:
fig,ax = plt.subplots()
ax.plot(q_ring_center**2, rate, 'ro', ls='')
ax.plot(q_ring_center**2, gmfit(q_ring_center**2), ls='-')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()
In [69]:
good_start= 0
good_end = 10000
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
In [70]:
g12b = auto_two_Array( imgsr, ring_mask, data_pixel = data_pixel )
In [71]:
g12_num = 0 #0: the firs ring
data = g12b[:,:,g12_num]
fig, ax = plt.subplots()
im=ax.imshow( data, origin='lower' , cmap='viridis',
norm= LogNorm( vmin= 1, vmax= 1.6 ),
extent=[0, data.shape[0]*timeperframe, 0, data.shape[0]*timeperframe ] )
ax.set_title('0-%s frames--Qth= %s'%(Nimg,g12_num))
ax.set_xlabel( r'$t_1$ $(s)$', fontsize = 18)
ax.set_ylabel( r'$t_2$ $(s)$', fontsize = 18)
fig.colorbar(im)
plt.show()
In [79]:
#np.where( g12b[:,:,0] == g12b[:,:,0].min() )
In [72]:
g12b_mask = make_g12_mask(bad_frames, g12b[:,:,0].shape)
In [73]:
g12bm = masked_g12( g12b, bad_frames)
In [74]:
g2b = get_one_time_from_two_time(g12bm)
In [75]:
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0:
sy = int(num_rings/sx)
else:
sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)
plt.axis('off')
for sn in range(num_rings):
ax = fig.add_subplot(sx,sy,sn+1 )
ax.set_ylabel("g2")
ax.set_title(" Q= " + '%.5f '%(q_ring_center[sn]) + r'$\AA^{-1}$')
y=g2b[:, sn]
ax.semilogx( np.arange(len(y))*timeperframe, y, '-o', markersize=6)
#ax.semilogx(lags, y, '-o', markersize=6)
#ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
ax.set_ylim( [1.0, max(y[1:])*1.05 ] )
plt.show()
fig.tight_layout()
In [76]:
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0:
sy = int(num_rings/sx)
else:
sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)
plt.axis('off')
for sn in range(num_rings):
ax = fig.add_subplot(sx,sy,sn+1 )
ax.set_ylabel("g2")
ax.set_title(" Q= " + '%.5f '%(q_ring_center[sn]) + r'$\AA^{-1}$')
y=g2b[:, sn]
ax.semilogx( np.arange(len(y))*timeperframe, y, '--r', markersize=6)
y2=g2[:, sn]
ax.semilogx(lags, y2, 'o', markersize=6)
#ax.semilogx(lags, y, '-o', markersize=6)
#ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
ax.set_ylim( [1.0, max(y2[1:])*1.05 ] )
plt.show()
fig.tight_layout()
In [77]:
g12b_norm, g12b_not_norm, norms = auto_two_Array_g1_norm( imgsr, ring_mask, data_pixel = data_pixel )
In [78]:
g12b_not_norm_m = masked_g12( g12b_not_norm, bad_frames )
In [79]:
g2b_norm = g2b
In [80]:
g2b_not_norm = get_one_time_from_two_time(g12b_not_norm_m, norms=norms, nopr=nopr )
In [81]:
#num_rings =6
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0:
sy = int(num_rings/sx)
else:
sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)
plt.axis('off')
for sn in range(num_rings):
ax = fig.add_subplot(sx,sy,sn+1 )
ax.set_ylabel("g2")
ax.set_title(" Q= " + '%.5f '%(q_ring_center[sn]) + r'$\AA^{-1}$')
y=g2b_norm[:, sn]
#ax.semilogx( np.arange(len(y))*timeperframe, y, '--r', markersize=3)
y3=g2b_not_norm[:, sn]
ax.semilogx( np.arange(len(y))*timeperframe, y3, '-g', lw=4, markersize=3)
y2=g2[:, sn]
#ax.semilogx(lags, y2, '-ko', markersize= 3)
ax.semilogx(lags, y2, 'o', markersize=6)
#ax.semilogx(lags, y, '-o', markersize=6)
#ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
ax.set_ylim( [1.0, max(y2[1:])*1.05 ] )
plt.show()
fig.tight_layout()
In [89]:
#plt.close('all')
In [82]:
g2_aged = get_aged_g2_from_g12q( g12bm[:,:,0], slice_num =3, slice_width= 500,
slice_start=4000, slice_end= 20000-4000 )
In [83]:
show_g12q_aged_g2( g12bm[:,:,0], g2_aged,slice_width= 500, timeperframe=1,vmin= 1, vmax= 1.22 )
In [85]:
taus = get_tau_from_g12q( g12b_norm[:,:,0], slice_num = 5, slice_width=1,
slice_start=3, slice_end= 5000-1 )
In [108]:
show_g12q_taus( g12b_norm[:,:,0], taus, slice_width=50, timeperframe=1,vmin=1.01,vmax=1.55 )
In [109]:
his = histogram_taus(taus, hisbin=30, plot=True, timeperframe=timeperframe)
In [110]:
bad_frames
Out[110]:
In [89]:
s1,s2 = 0,2000
In [90]:
g4 = get_four_time_from_two_time( g12bm, g2b, rois=[s1,s2,s1,s2] )
In [111]:
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0:
sy = int(num_rings/sx)
else:
sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)
plt.axis('off')
for sn in range(num_rings):
ax = fig.add_subplot(sx,sy,sn+1 )
ax.set_ylabel("g4")
ax.set_title(" Q= " + '%.5f '%(q_ring_center[sn]) + r'$\AA^{-1}$')
y=g4[:, sn]
ax.semilogx( np.arange(len(y))*timeperframe, y, '-o', markersize=6)
#ax.semilogx(lags, y, '-o', markersize=6)
ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
#ax.set_ylim( [1.0, max(y[1:])*1.05 ] )
plt.show()
fig.tight_layout()