-
Notifications
You must be signed in to change notification settings - Fork 60
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
X-Ray Speckle Visibility Spectroscopy #293
Changes from all commits
fb94bb8
6baa470
558270d
a735748
147caf8
d26e986
61341f8
95470d6
d59d965
62a7100
75fb1bf
882a3a7
9212960
aae0f9b
1be549c
f81441b
a5468a8
c33dd73
d7dffee
78791fa
bb3aeb3
4fb539e
3e3b995
8bc9d6d
2be429d
7584f08
48cd057
1d28ed1
c98f0a5
180db89
4b783d4
7f75022
be553a8
5811f1e
de067b9
64a31a9
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,290 @@ | ||
# ###################################################################### | ||
# Copyright (c) 2014, Brookhaven Science Associates, Brookhaven # | ||
# National Laboratory. All rights reserved. # | ||
# # | ||
# Developed at the NSLS-II, Brookhaven National Laboratory # | ||
# Developed by Sameera K. Abeykoon and Yugang Zhang, June 2015 # | ||
# # | ||
# Redistribution and use in source and binary forms, with or without # | ||
# modification, are permitted provided that the following conditions # | ||
# are met: # | ||
# # | ||
# * Redistributions of source code must retain the above copyright # | ||
# notice, this list of conditions and the following disclaimer. # | ||
# # | ||
# * Redistributions in binary form must reproduce the above copyright # | ||
# notice this list of conditions and the following disclaimer in # | ||
# the documentation and/or other materials provided with the # | ||
# distribution. # | ||
# # | ||
# * Neither the name of the Brookhaven Science Associates, Brookhaven # | ||
# National Laboratory nor the names of its contributors may be used # | ||
# to endorse or promote products derived from this software without # | ||
# specific prior written permission. # | ||
# # | ||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS # | ||
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # | ||
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS # | ||
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE # | ||
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, # | ||
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES # | ||
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR # | ||
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) # | ||
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, # | ||
# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OTHERWISE) ARISING # | ||
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # | ||
# POSSIBILITY OF SUCH DAMAGE. # | ||
######################################################################## | ||
|
||
""" | ||
X-ray speckle visibility spectroscopy(XSVS) - Dynamic information of | ||
the speckle patterns are obtained by analyzing the speckle statistics | ||
and calculating the speckle contrast in single scattering patterns. | ||
|
||
This module will provide XSVS analysis tools | ||
""" | ||
|
||
from __future__ import (absolute_import, division, print_function) | ||
import six | ||
import numpy as np | ||
import time | ||
|
||
from . import correlation as corr | ||
from . import roi | ||
from .utils import bin_edges_to_centers, geometric_series | ||
|
||
import logging | ||
logger = logging.getLogger(__name__) | ||
|
||
|
||
def xsvs(image_sets, label_array, number_of_img, timebin_num=2, | ||
max_cts=None): | ||
""" | ||
This function will provide the probability density of detecting photons | ||
for different integration times. | ||
|
||
The experimental probability density P(K) of detecting photons K is | ||
obtained by histogramming the speckle counts over an ensemble of | ||
equivalent pixels and over a number of speckle patterns recorded | ||
with the same integration time T under the same condition. | ||
|
||
Parameters | ||
---------- | ||
image_sets : array | ||
sets of images | ||
label_array : array | ||
labeled array; 0 is background. | ||
Each ROI is represented by a distinct label (i.e., integer). | ||
number_of_img : int | ||
number of images (how far to go with integration times when finding | ||
the time_bin, using skxray.utils.geometric function) | ||
timebin_num : int, optional | ||
integration time; default is 2 | ||
max_cts : int, optional | ||
the brightest pixel in any ROI in any image in the image set. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. What does this default to if no value is passed in? |
||
defaults to using skxray.core.roi.roi_max_counts to determine the brightest | ||
pixel in any of the ROIs | ||
|
||
Returns | ||
------- | ||
prob_k_all : array | ||
probability density of detecting photons | ||
prob_k_std_dev : array | ||
standard deviation of probability density of detecting photons | ||
|
||
Notes | ||
----- | ||
These implementation is based on following references | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Can you put in a PR with your speckle notebook against scikit-xray-examples and add another bullet to this section noting that there is an example in https://github.com/scikit-xray/scikit-xray-examples with that notebook? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Sure I will put a PR with a speckle notebook against scikit-xray-examples |
||
References: text [1]_, text [2]_ | ||
|
||
.. [1] 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. | ||
|
||
.. [2] R. Bandyopadhyay, A. S. Gittings, S. S. Suh, P.K. Dixon and | ||
D.J. Durian "Speckle-visibilty Spectroscopy: A tool to study | ||
time-varying dynamics" Rev. Sci. Instrum. vol 76, p 093110, 2005. | ||
|
||
There is an example in https://github.com/scikit-xray/scikit-xray-examples | ||
It will demonstrate the use of these functions in this module for | ||
experimental data. | ||
|
||
""" | ||
if max_cts is None: | ||
max_cts = roi.roi_max_counts(image_sets, label_array) | ||
|
||
# find the label's and pixel indices for ROI's | ||
labels, indices = corr.extract_label_indices(label_array) | ||
|
||
# number of ROI's | ||
u_labels = list(np.unique(labels)) | ||
num_roi = len(u_labels) | ||
|
||
# create integration times | ||
time_bin = geometric_series(timebin_num, number_of_img) | ||
|
||
# number of times in the time bin | ||
num_times = len(time_bin) | ||
|
||
# number of pixels per ROI | ||
num_pixels = np.bincount(labels, minlength=(num_roi+1))[1:] | ||
|
||
# probability density of detecting photons | ||
prob_k_all = np.zeros([num_times, num_roi], dtype=np.object) | ||
|
||
# square of probability density of detecting photons | ||
prob_k_pow_all = np.zeros_like(prob_k_all) | ||
|
||
# standard deviation of probability density of detecting photons | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The docstring says that this is the |
||
prob_k_std_dev = np.zeros_like(prob_k_all) | ||
|
||
# get the bin edges for each time bin for each ROI | ||
bin_edges = np.zeros(prob_k_all.shape[0], dtype=prob_k_all.dtype) | ||
for i in range(num_times): | ||
bin_edges[i] = np.arange(max_cts*2**i) | ||
|
||
start_time = time.time() # used to log the computation time (optionally) | ||
|
||
for i, images in enumerate(image_sets): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Why are we accepting an iterable of image_sets as the input for this function? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @ericdill Input can be different image_sets collected at different points of the sample as we discussed in your xsvs_plotting tools PR in X-ray-vision |
||
# Ring buffer, a buffer with periodic boundary conditions. | ||
# Images must be keep for up to maximum delay in buf. | ||
buf = np.zeros([num_times, timebin_num], | ||
dtype=np.object) # matrix of buffers | ||
|
||
# to track processing each time level | ||
track_level = np.zeros(num_times) | ||
|
||
# to increment buffer | ||
cur = np.full(num_times, timebin_num) | ||
|
||
# to track how many images processed in each level | ||
img_per_level = np.zeros(num_times, dtype=np.int64) | ||
|
||
prob_k = np.zeros_like(prob_k_all) | ||
prob_k_pow = np.zeros_like(prob_k_all) | ||
|
||
for n, img in enumerate(images): | ||
cur[0] = (1 + cur[0]) % timebin_num | ||
# read each frame | ||
# Put the image into the ring buffer. | ||
buf[0, cur[0] - 1] = (np.ravel(img))[indices] | ||
|
||
_process(num_roi, 0, cur[0] - 1, buf, img_per_level, labels, | ||
max_cts, bin_edges[0], prob_k, prob_k_pow) | ||
|
||
# check whether the number of levels is one, otherwise | ||
# continue processing the next level | ||
level = 1 | ||
|
||
while level < num_times: | ||
if not track_level[level]: | ||
track_level[level] = 1 | ||
else: | ||
prev = 1 + (cur[level - 1] - 2) % timebin_num | ||
cur[level] = 1 + cur[level] % timebin_num | ||
|
||
buf[level, cur[level]-1] = (buf[level-1, | ||
prev-1] + | ||
buf[level-1, | ||
cur[level - 1] - 1]) | ||
track_level[level] = 0 | ||
|
||
_process(num_roi, level, cur[level]-1, buf, img_per_level, | ||
labels, max_cts, bin_edges[level], prob_k, | ||
prob_k_pow) | ||
level += 1 | ||
|
||
prob_k_all += (prob_k - prob_k_all)/(i + 1) | ||
prob_k_pow_all += (prob_k_pow - prob_k_pow_all)/(i + 1) | ||
|
||
prob_k_std_dev = np.power((prob_k_pow_all - | ||
np.power(prob_k_all, 2)), .5) | ||
|
||
logger.info("Processing time for XSVS took %s seconds." | ||
"", (time.time() - start_time)) | ||
return prob_k_all, prob_k_std_dev | ||
|
||
|
||
def _process(num_roi, level, buf_no, buf, img_per_level, labels, | ||
max_cts, bin_edges, prob_k, prob_k_pow): | ||
""" | ||
Internal helper function. This modifies inputs in place. | ||
|
||
This helper function calculate probability of detecting photons for | ||
each integration time. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Can you copy this pattern and add a There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Great, thanks |
||
|
||
.. warning :: This function mutates the input values. | ||
|
||
Parameters | ||
---------- | ||
num_roi : int | ||
number of ROI's | ||
level : int | ||
current time level(integration time) | ||
buf_no : int | ||
current buffer number | ||
buf : array | ||
image data array to use for XSVS | ||
img_per_level : int | ||
to track how many images processed in each level | ||
labels : array | ||
labels of the required region of interests(ROI's) | ||
max_cts: int | ||
maximum pixel count | ||
bin_edges : array | ||
bin edges for each integration times and each ROI | ||
prob_k : array | ||
probability density of detecting photons | ||
prob_k_pow : array | ||
squares of probability density of detecting photons | ||
""" | ||
img_per_level[level] += 1 | ||
u_labels = list(np.unique(labels)) | ||
|
||
for j, label in enumerate(u_labels): | ||
roi_data = buf[level, buf_no][labels == label] | ||
spe_hist, bin_edges = np.histogram(roi_data, bins=bin_edges, | ||
density=True) | ||
spe_hist = np.nan_to_num(spe_hist) | ||
prob_k[level, j] += (spe_hist - | ||
prob_k[level, j])/(img_per_level[level]) | ||
|
||
prob_k_pow[level, j] += (np.power(spe_hist, 2) - | ||
prob_k_pow[level, j])/(img_per_level[level]) | ||
|
||
|
||
def normalize_bin_edges(num_times, num_rois, mean_roi, max_cts): | ||
""" | ||
This will provide the normalized bin edges and bin centers for each | ||
integration time. | ||
|
||
Parameters | ||
---------- | ||
num_times : int | ||
number of integration times for XSVS | ||
num_rois : int | ||
number of ROI's | ||
mean_roi : array | ||
mean intensity of each ROI | ||
shape (number of ROI's) | ||
max_cts : int | ||
maximum pixel counts | ||
|
||
Returns | ||
------- | ||
norm_bin_edges : array | ||
normalized speckle count bin edges | ||
shape (num_times, num_rois) | ||
norm_bin_centers :array | ||
normalized speckle count bin centers | ||
shape (num_times, num_rois) | ||
""" | ||
norm_bin_edges = np.zeros((num_times, num_rois), dtype=object) | ||
norm_bin_centers = np.zeros_like(norm_bin_edges) | ||
for i in range(num_times): | ||
for j in range(num_rois): | ||
norm_bin_edges[i, j] = np.arange(max_cts*2**i)/(mean_roi[j]*2**i) | ||
norm_bin_centers[i, j] = bin_edges_to_centers(norm_bin_edges[i, j]) | ||
|
||
return norm_bin_edges, norm_bin_centers |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,91 @@ | ||
# ###################################################################### | ||
# Copyright (c) 2014, Brookhaven Science Associates, Brookhaven # | ||
# National Laboratory. All rights reserved. # | ||
# # | ||
# Redistribution and use in source and binary forms, with or without # | ||
# modification, are permitted provided that the following conditions # | ||
# are met: # | ||
# # | ||
# * Redistributions of source code must retain the above copyright # | ||
# notice, this list of conditions and the following disclaimer. # | ||
# # | ||
# * Redistributions in binary form must reproduce the above copyright # | ||
# notice this list of conditions and the following disclaimer in # | ||
# the documentation and/or other materials provided with the # | ||
# distribution. # | ||
# # | ||
# * Neither the name of the Brookhaven Science Associates, Brookhaven # | ||
# National Laboratory nor the names of its contributors may be used # | ||
# to endorse or promote products derived from this software without # | ||
# specific prior written permission. # | ||
# # | ||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS # | ||
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # | ||
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS # | ||
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE # | ||
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, # | ||
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES # | ||
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR # | ||
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) # | ||
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, # | ||
# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OTHERWISE) ARISING # | ||
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # | ||
# POSSIBILITY OF SUCH DAMAGE. # | ||
######################################################################## | ||
from __future__ import absolute_import, division, print_function | ||
import logging | ||
|
||
import numpy as np | ||
from numpy.testing import (assert_array_almost_equal, | ||
assert_almost_equal) | ||
|
||
from skimage.morphology import convex_hull_image | ||
|
||
from .. import speckle as xsvs | ||
from .. import roi as roi | ||
from skxray.testing.decorators import skip_if | ||
|
||
logger = logging.getLogger(__name__) | ||
|
||
|
||
def test_xsvs(): | ||
images = [] | ||
for i in range(5): | ||
int_array = np.tril((i + 2) * np.ones(10)) | ||
int_array[int_array == 0] = (i + 1) | ||
images.append(int_array) | ||
|
||
images_sets = [np.asarray(images), ] | ||
roi_data = np.array(([4, 2, 2, 2], [0, 5, 4, 4]), dtype=np.int64) | ||
label_array = roi.rectangles(roi_data, shape=images[0].shape) | ||
|
||
prob_k_all, std = xsvs.xsvs(images_sets, label_array, timebin_num=2, | ||
number_of_img=5, max_cts=None) | ||
|
||
assert_array_almost_equal(prob_k_all[0, 0], | ||
np.array([0., 0., 0.2, 0.2, 0.4])) | ||
assert_array_almost_equal(prob_k_all[0, 1], | ||
np.array([0., 0.2, 0.2, 0.2, 0.4])) | ||
|
||
|
||
def test_normalize_bin_edges(): | ||
num_times = 3 | ||
num_rois = 2 | ||
mean_roi = np.array([2.5, 4.0]) | ||
max_cts = 5 | ||
|
||
bin_edges, bin_cen = xsvs.normalize_bin_edges(num_times, num_rois, | ||
mean_roi, max_cts) | ||
|
||
assert_array_almost_equal(bin_edges[0, 0], np.array([0., 0.4, 0.8, | ||
1.2, 1.6])) | ||
|
||
assert_array_almost_equal(bin_edges[2, 1], np.array([0., 0.0625, 0.125, | ||
0.1875, 0.25, 0.3125, | ||
0.375, 0.4375, 0.5, | ||
0.5625, 0.625, 0.6875, | ||
0.75, 0.8125, 0.875, | ||
0.9375, 1., 1.0625, | ||
1.125, 1.1875])) | ||
|
||
assert_array_almost_equal(bin_cen[0, 0], np.array([0.2, 0.6, 1., 1.4])) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The docstring is out of sync with the variable name. I imagine this should be called
image_sets
. The docstring should explain that this is an iterable of iterables of images where the outer iterable represents different points in the sample and the inner iterable represents the images collected at each point in the sample.