X-ray Speckle Visibility Spectroscopy

The analysis module "skxray/core/speckle" https://github.com/scikit-xray/scikit-xray/blob/master/skxray/core/speckle.py



In [1]:

    
import xray_vision
import xray_vision.mpl_plotting as mpl_plot

import skxray.core.speckle as xsvs
import skxray.core.roi as roi
import skxray.core.correlation as corr
import skxray.core.utils as utils

import numpy as np
import os, sys

import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib.colors import LogNorm
from mpl_toolkits.axes_grid1.axes_grid import ImageGrid









    



:0: FutureWarning: IPython widgets are experimental and may change in the future.

Easily switch between interactive and static matplotlib plots



In [16]:

    
interactive_mode = False

if interactive_mode:
    %matplotlib notebook
else:
    %matplotlib inline

backend = mpl.get_backend()

This data provided by Dr. Andrei Fluerasu

L. Li, P. Kwasniewski, D. Oris, L Wiegart, L. Cristofolini, C. Carona and A. Fluerasu , "Photon statistics and speckle visibility spectroscopy with partially coherent x-rays" J. Synchrotron Rad., vol 21, p 1288-1295, 2014.



In [3]:

    
%run download.py









    



__file__
current_folder = /home/sameera/NSLs2/scikit-xray-examples/demos/speckle
current_folder = /home/sameera/NSLs2/scikit-xray-examples/demos/speckle



In [4]:

    
data_dir = "Duke_data/"
duke_rdata = np.load(data_dir+"duke_img_1_5000.npy")
duke_dark = np.load(data_dir+"duke_dark.npy")

duke_data = []
for i in range(duke_rdata.shape[0]):
    duke_data.append(duke_rdata[i] - duke_dark)
    
duke_ndata=np.asarray(duke_data)

# load the mask(s) and mask the data
mask1 = np.load(data_dir+"new_mask4.npy")
mask2 = np.load(data_dir+"Luxi_duke_mask.npy")

N_mask = ~(mask1 + mask2)



In [5]:

    
mask_data = N_mask*duke_ndata



In [6]:

    
#  get the average image
avg_img = np.average(duke_ndata, axis=0)

# if matplotlib version 1.5 or later
if float('.'.join(mpl.__version__.split('.')[:2])) >= 1.5:
    cmap = 'viridis'
else:
    cmap = 'CMRmap'

# plot the average image data after masking
plt.figure()
plt.imshow(N_mask*avg_img, vmax=1e0, cmap=cmap)
plt.title("Averaged masked data for Duke Silica Gel ")
plt.colorbar()
plt.show()

Create the Rings Mask

Use the skxray.core.roi module to create Ring ROIs (ROI Mask).¶ (https://github.com/scikit-xray/scikit-xray/blob/master/skxray/core/roi.py)



In [7]:

    
inner_radius = 26  # radius of the first ring
width = 1        # width of each ring
spacing = 0      # no spacing between rings
num_rings = 4    # number of rings
center = (133, 143)   # center of the spckle pattern

#  find the edges of the required rings
edges = roi.ring_edges(inner_radius, width, spacing, num_rings)
edges









    Out[7]:





array([[ 26.,  27.],
       [ 27.,  28.],
       [ 28.,  29.],
       [ 29.,  30.]])

Convert the edge values of the rings to q ( reciprocal space)



In [8]:

    
dpix = 0.055  # The physical size of the pixels

lambda_ = 1.5498  # wavelength of the X-rays
Ldet = 2200.   #   # detector to sample distance

two_theta = utils.radius_to_twotheta(Ldet, edges*dpix)
q_val = utils.twotheta_to_q(two_theta, lambda_)

q_val









    Out[8]:





array([[ 0.00263522,  0.00273658],
       [ 0.00273658,  0.00283793],
       [ 0.00283793,  0.00293929],
       [ 0.00293929,  0.00304064]])



In [9]:

    
q_ring = np.mean(q_val, axis=1)
q_ring









    Out[9]:





array([ 0.0026859 ,  0.00278726,  0.00288861,  0.00298997])

Create a labeled array using roi.rings



In [11]:

    
rings = roi.rings(edges, center, avg_img.shape)

images_sets = (mask_data, )

ring_mask = rings*N_mask

# plot the figure
fig, axes = plt.subplots(figsize=(5, 5))
axes.set_title("Ring Mask")
im = mpl_plot.show_label_array(axes, ring_mask, cmap=cmap)
plt.show()

Find the brightest pixel in any ROI in any image in the image set.

Using roi_max_counts function from skxray.core.roi module



In [12]:

    
max_cts = roi.roi_max_counts(images_sets, ring_mask)
max_cts









    Out[12]:





24

Everything looks good, next X-ray speckle visibilty spectroscopy

This function will provide the probability density of detecting photons for different integration time. Using skxray.core.speckle module



In [13]:

    
spe_cts_all, std_dev = xsvs.xsvs(images_sets, ring_mask, timebin_num=2,
                             number_of_img=30, max_cts=max_cts)

Find the integration times

using skxray.core.utils.geometric_series



In [14]:

    
time_steps = utils.geometric_series(2, 30)
time_steps









    Out[14]:





[1, 2, 4, 8, 16]

Get the mean intensity of each ring



In [15]:

    
mean_int_sets, index_list = roi.mean_intensity(mask_data, ring_mask)



In [17]:

    
plt.figure(figsize=(8, 8))
plt.title("Mean intensity of each ring")
for i in range(num_rings):
    plt.plot(mean_int_sets[:,i], label="Ring "+str(i+1))
plt.legend()    
plt.show()



In [18]:

    
mean_int_ring = np.mean(mean_int_sets, axis=0)
mean_int_ring









    Out[18]:





array([ 3.17503051,  2.03199322,  1.22062687,  0.78536136])

Get the normalized bin edges and bin centers for each integration time.

using skxray.core.speckle.normalize_bin_edges



In [19]:

    
num_times = 6
num_rois=num_rings
norm_bin_edges, norm_bin_centers = xsvs.normalize_bin_edges(num_times,
                                                            num_rois, mean_int_ring, max_cts)

1st q ring 0.0026859 (1/Angstroms)



In [20]:

    
fig, axes = plt.subplots(figsize=(6, 6))
axes.set_xlabel("K/<K>")
axes.set_ylabel("P(K)")
for i in range(4):
    art, = axes.plot(norm_bin_edges[i, 0][:-1], spe_cts_all[i, 0], '-o', label=str(time_steps[i])+" ms")
    axes.set_xlim(0, 4)
    axes.legend()
plt.title("1st q ring 0.0026859 (1/Angstroms)")
plt.show()

2nd q ring 0.00278726 (1/Angstroms)



In [21]:

    
fig, axes = plt.subplots(figsize=(6, 6))
axes.set_xlabel("K/<K>")
axes.set_ylabel("P(K)")
for i in range(4):
    art, = axes.plot(norm_bin_edges[i, 1][:-1], spe_cts_all[i, 1], '-o', label=str(time_steps[i])+" ms")
    axes.legend()
    axes.set_xlim(0, 4)
plt.title("2nd q ring 0.00278726 (1/Angstroms)")
plt.show()

X-ray speckle visibilty spectroscopy(XSVS) for differnt time steps

This function will provide the probability density of detecting photons for different integration time. Using skxray.core.speckle module

Find the new integration times

using skxray.core.utils.geometric_series



In [22]:

    
time_steps_5 = utils.geometric_series(5, 50)
time_steps_5









    Out[22]:





[1, 5, 25]

XSVS results for new integartion times 1ms, 5ms and 25ms



In [23]:

    
p_K, std_dev_5 = xsvs.xsvs(images_sets, ring_mask, timebin_num=5,
                             number_of_img=50, max_cts=max_cts)

Plot the results for each Q ring

1st q ring 0.0026859 (1/Angstroms)



In [24]:

    
fig, axes = plt.subplots(figsize=(6, 6))
axes.set_xlabel("K/<K>")
axes.set_ylabel("P(K)")
for i in range(3):
    art, = axes.plot(norm_bin_edges[i, 0][:-1], p_K[i, 0], '-o', label=str(time_steps_5[i])+" ms")
    axes.set_xlim(0, 4)
    axes.legend()
plt.title("1st q ring 0.0026859 (1/Angstroms)")
plt.show()

2nd q ring 0.00278726 (1/Angstroms)



In [25]:

    
fig, axes = plt.subplots(figsize=(6, 6))
axes.set_xlabel("K/<K>")
axes.set_ylabel("P(K)")
for i in range(3):
    art, = axes.plot(norm_bin_edges[i, 1][:-1], p_K[i, 1], '-o', label=str(time_steps_5[i])+" ms")
    axes.legend()
    axes.set_xlim(0, 4)
plt.title("2nd q ring 0.00278726 (1/Angstroms)")
plt.show()

3rd q ring 0.00288861 (1/ Angstroms)



In [26]:

    
fig, axes = plt.subplots(figsize=(6, 6))
axes.set_xlabel("K/<K>")
axes.set_ylabel("P(K)")
for i in range(3):
    art, = axes.plot(norm_bin_edges[i, 2][:-1], p_K[i, 2], '-o', label=str(time_steps_5[i])+" ms" )
    axes.set_xlim(0, 4)
    axes.legend()
plt.title("3rd q ring 0.00288861 (1/ Angstroms)")
plt.show()

4th q ring 0.0298997 (1/ Angstroms)



In [27]:

    
fig, axes = plt.subplots(figsize=(6, 6))
axes.set_xlabel("K/<K>")
axes.set_ylabel("P(K)")
for i in range(3):
    art, = axes.plot(norm_bin_edges[i, 3][:-1], p_K[i, 3], '-o', label=str(time_steps_5[i])+" ms")
    axes.set_xlim(0, 4)
    axes.legend()
plt.title("4th q ring 0.0298997 (1/ Angstroms)")
plt.show()



In [28]:

    
import skxray
print(skxray.__version__)









    



0.0.5+13.gbd24fbb



In [ ]: