"This notebook corresponds to version {{ version }} of the pipeline tool: https://github.com/NSLS-II/pipelines"
In [612]:
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 [4]:
hdr = db[{{ uid }}]
In [6]:
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 [5]:
cd /home/yuzhang/chx-pipelines/Develops/
In [ ]:
cd /XF11ID/analysis/Analysis_Pipelines/Develop/
In [7]:
%run develop.py
In [8]:
%matplotlib notebook
#%matplotlib inline
In [9]:
BlueScan = True
DirectAcq = False
detector = 'eiger_4M_cam_img_image_lightfield' #for 4M
In [10]:
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 [16]:
if BlueScan:
hdr = db[uid]
ev, = get_events( hdr, [detector] )
imgs = ev['data'][detector]
else:
imgs = Images(uid)
print (imgs)
Nimg=len(imgs)
In [17]:
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[17]:
In [18]:
imgs.md
Out[18]:
In [19]:
# 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[19]:
In [20]:
mask = np.load( path + str(uid)+ "_mask.npy")
In [21]:
maskr = mask[::-1,:]
In [613]:
fig, ax = plt.subplots()
im=ax.imshow(maskr, origin='lower' ,vmin=0,vmax=1,cmap='viridis')
fig.colorbar(im)
plt.show()
In [23]:
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 [24]:
#interact(view_image, i=(0, Nimg-1))
In [25]:
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 [26]:
kymo_sum = np.load( path + str(uid)+"_kymo_sum.npy" )
In [614]:
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 [615]:
avg_img = np.load( path + str(uid)+"_avg_img.npy" )
avg_imgm = avg_img * mask
In [616]:
avg_imgr = avg_img[::-1,:]
avg_imgmr = avg_imgm[::-1,:]
In [617]:
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 [618]:
imgs.md['beam_center_x'], imgs.md['beam_center_y']
Out[618]:
In [619]:
#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[619]:
In [620]:
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 [622]:
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 [623]:
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 [624]:
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 [625]:
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[625]:
In [626]:
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[626]:
In [627]:
rings = roi.rings(edges, center, avg_imgmr.shape)
ring_mask = rings*maskr
In [628]:
ring_mask
Out[628]:
In [629]:
qind, pixelist = roi.extract_label_indices( ring_mask )
noqs = len( np.unique(qind) )
nopr = np.bincount(qind, minlength=(noqs+1))[1:]
In [630]:
qind
Out[630]:
In [50]:
nopr
Out[50]:
In [51]:
pixel = roi.roi_pixel_values(avg_imgmr, ring_mask, [2] )
fig,ax=plt.subplots()
ax.plot( pixel[0][0] ,'bo', ls='-' )
Out[51]:
In [631]:
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 [638]:
#plt.close('all')
In [632]:
# 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 [56]:
imgs_ =imgs
imgsr = Reverse_Coordinate(imgs_, mask)
In [58]:
max_inten_ring =2
In [59]:
#kymo = roi.kymograph(imgsr, ring_mask, num = max_inten_ring)
In [60]:
t0 = time.time()
data_pixel = Get_Pixel_Array( imgsr, pixelist).get_data()
run_time(t0)
In [61]:
pixelist_qi = np.where( qind == 2)[0]
data_pixel_qi = data_pixel[:,pixelist_qi]
In [639]:
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 [63]:
mean_inten = get_mean_intensity( data_pixel, qind)
In [634]:
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()
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/skxray/core/correlation.py
In [65]:
nopr
Out[65]:
In [66]:
edges
Out[66]:
In [67]:
qind
Out[67]:
In [68]:
ring_mask = np.array( ring_mask, dtype=int)
In [636]:
good_start = 0
good_end = 7600
imgs_ =imgs[good_start: good_end-1]
imgsr2 = Reverse_Coordinate(imgs_, mask)
In [637]:
num_lev = 15
num_buf = 8
t0 = time.time()
g2_0_7600, lag_steps = corr.multi_tau_auto_corr(num_lev, num_buf,ring_mask, imgsr2)
run_time(t0)
In [645]:
g2_0_7600 = g2
In [643]:
lags = lag_steps*timeperframe
In [644]:
tg2 = np.hstack( [ lags.reshape( len(lags),1), g2_0_7600] )
#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 [646]:
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 [77]:
from lmfit import Model
mod = Model(corr.auto_corr_scat_factor)
In [647]:
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.02)
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 [648]:
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 [649]:
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 [650]:
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 [82]:
%run two_time.py
In [657]:
good_start= 0
good_end = 10000
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
In [658]:
g12b = auto_two_Array( imgsr, ring_mask, data_pixel = data_pixel )
In [659]:
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 [ ]:
#np.where( g12b[:,:,0] == g12b[:,:,0].min() )
In [660]:
g12b_mask = make_g12_mask(g12b,[7629])
g12bm = masked_g12(g12b,[7629])
In [661]:
g2b = get_one_time_from_two_time(g12bm, timeperframe= timeperframe)
In [662]:
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 [663]:
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 [664]:
g12b_norm, g12b_not_norm, norms = auto_two_Array_g1_norm( imgsr, ring_mask, data_pixel = data_pixel )
In [667]:
g2b_norm = get_one_time_from_two_time(g12b_norm,
norms=None, nopr=nopr, timeperframe= timeperframe)
In [665]:
g2b_not_norm = get_one_time_from_two_time(g12b_not_norm,
norms=norms, nopr=nopr, timeperframe= timeperframe)
In [670]:
#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 [681]:
good_start= 0
good_end = 100
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start
num_buf = 4
In [682]:
g2_t2, lag_steps_t = autocor_one_time( num_buf, ring_mask, imgsr,
num_lev=None, bad_images=[ 94,95,96,97,98,99,100] )
In [675]:
good_start= 0
good_end = 95
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start
num_buf = 4
In [676]:
g2_t3, lag_steps_t3 =autocor_one_time( num_buf, ring_mask, imgsr,
num_lev=None, bad_images=[ ] )
In [677]:
g2_t2[:,0]
Out[677]:
In [678]:
g2_t3[:,0]
Out[678]:
In [ ]:
In [683]:
g12_num =0
g2b = show_g12_age_cuts( g12bm, g12_num=g12_num, slice_num =3, slice_width= 500,
slice_start=4000, slice_end= 20000-4000,
timeperframe= timeperframe,vmin= 1.0, vmax = 1.22) #vmax=g12b[:,:,g12_num].max() )
In [684]:
taus = show_g12_tau_cuts(g12b_norm, g12_num=0, slice_num = 5, slice_width= 1.0,
slice_start=3, slice_end= 5000-1,draw_scale_tau =1.0, vmin=1.01,vmax=1.55 )
In [686]:
his = his_tau(taus, hisbin = 10, plot=True)
In [ ]: