diff --git a/skbeam_copied/.coveragerc b/skbeam_copied/.coveragerc new file mode 100644 index 00000000..906b7275 --- /dev/null +++ b/skbeam_copied/.coveragerc @@ -0,0 +1,13 @@ +[run] +source=skbeam +[report] +omit = + */python?.?/* + */site-packages/nose/* + skbeam/__init__.py + # ignore _version.py and versioneer.py + .*version.* + *_version.py + +exclude_lines = + def set_default diff --git a/skbeam_copied/.flake8 b/skbeam_copied/.flake8 new file mode 100644 index 00000000..96d71e3f --- /dev/null +++ b/skbeam_copied/.flake8 @@ -0,0 +1,11 @@ +[flake8] +exclude = + .git, + __pycache__, + build, + dist, + versioneer.py, + pyxrf/_version.py, + doc/conf.py, + doc/sphinxext/* +max-line-length = 115 diff --git a/skbeam_copied/.gitattributes b/skbeam_copied/.gitattributes new file mode 100644 index 00000000..4951c496 --- /dev/null +++ b/skbeam_copied/.gitattributes @@ -0,0 +1 @@ +skbeam/_version.py export-subst diff --git a/skbeam_copied/.gitignore b/skbeam_copied/.gitignore new file mode 100644 index 00000000..f9e50e75 --- /dev/null +++ b/skbeam_copied/.gitignore @@ -0,0 +1,82 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +venv/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +coverage.xml +*,cover +.hypothesis/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log + +# Sphinx documentation +doc/_build/ +doc/resource/api/core/_as_gen/ + +# pytest +.pytest_cache/ + +# PyBuilder +target/ + +# Editor files +#mac +.DS_Store +*~ + +#vim +*.swp +*.swo + +#pycharm +.idea/* + + +#Ipython Notebook +.ipynb_checkpoints + +# Project-specific ignores: +skbeam/core/accumulators/histogram.c diff --git a/skbeam_copied/.travis.yml b/skbeam_copied/.travis.yml new file mode 100644 index 00000000..ebd95583 --- /dev/null +++ b/skbeam_copied/.travis.yml @@ -0,0 +1,145 @@ +language: python +sudo: false + +addons: + apt_packages: + - pandoc +env: + global: + - SUBMIT_CODECOV: false + - BUILD_DOCS: false + - GH_REF: github.com/scikit-beam/scikit-beam.git + - secure: "KzntGlmLDUZlAB8ZiSRBtONJypv4atrnFxl+THCW8kg6YbiODcCxG9cez7mfmh63wpJ54+zeGvgGljeVvjb7bjB2p+4hZy+mLoP9xoE1w5qF4XaQ5YTOYaSQqCeWa5cEIi+854gFX7gOfW5qL+EJf7h3HtmndI15zaU5gOPjSDU=" + +matrix: + include: + - python: 3.7 + env: BUILD_DOCS=true SUBMIT_CODECOV=true NUMPY=1.17 + - os: osx + language: generic + env: TRAVIS_PYTHON_VERSION=3.6 NUMPY=1.17 + - os: osx + language: generic + env: TRAVIS_PYTHON_VERSION=3.7 NUMPY=1.17 + +python: + - 3.6 + - 3.7 + +env: + - NUMPY=1.16 + - NUMPY=1.17 + - NUMPY=1.18 + +before_install: + - | + set -e + if [ "$TRAVIS_OS_NAME" == "linux" ]; then + arch="Linux" + elif [ "$TRAVIS_OS_NAME" == "osx" ]; then + arch="MacOSX" + else + echo "Unknown arch $TRAVIS_OS_NAME" + exit 1 + fi + wget https://repo.continuum.io/miniconda/Miniconda3-latest-${arch}-x86_64.sh -O miniconda.sh + chmod +x miniconda.sh + ./miniconda.sh -b -p ~/mc + source ~/mc/etc/profile.d/conda.sh + conda update conda --yes + export CONDARC=ci/condarc + +install: + - export GIT_FULL_HASH=`git rev-parse HEAD` + - conda create -n testenv python=$TRAVIS_PYTHON_VERSION numpy=$NUMPY + - conda activate testenv + - conda install xraylib # this is the only package not available on PyPI + - pip install -r requirements.txt + - pip install -r requirements-dev.txt + # # need to build_ext -i for the tests so that the .so is local to the source + # # code. We could also setup.py develop, but I'm not sure if that is any + # # better + - python setup.py install build_ext -i + - pip install codecov + - conda list + - pip list + +before_script: + # define a merge size check function that ensures no huge file was committed + # to the repo by accident + - | + size_check() { + GIT_TARGET_EXTRA="+refs/heads/${TRAVIS_BRANCH}" + GIT_SOURCE_EXTRA="+refs/pull/${TRAVIS_PULL_REQUEST}/merge" + echo TRAVIS_BRANCH=$TRAVIS_BRANCH + echo TRAVIS_PULL_REQUEST=$TRAVIS_PULL_REQUEST + echo TRAVIS_REPO_SLUG=$TRAVIS_REPO_SLUG + echo GIT_TARGET_EXTRA=$GIT_TARGET_EXTRA + echo GIT_SOURCE_EXTRA=$GIT_SOURCE_EXTRA + mkdir ~/repo-clone + pushd ~/repo-clone + git init && git remote add -t ${TRAVIS_BRANCH} origin git://github.com/${TRAVIS_REPO_SLUG}.git + git fetch origin ${GIT_TARGET_EXTRA} + git checkout -qf FETCH_HEAD + git tag travis-merge-target + git gc --aggressive + TARGET_SIZE=`du -s . | sed -e "s/\t.*//"` + echo TARGET_SIZE=$TARGET_SIZE + git pull origin ${GIT_SOURCE_EXTRA} + git gc --aggressive + MERGE_SIZE=`du -s . | sed -e "s/\t.*//"` + popd + echo MERGE_SIZE=$MERGE_SIZE + if [ "${MERGE_SIZE}" != "${TARGET_SIZE}" ]; then + SIZE_DIFF=`expr \( ${MERGE_SIZE} - ${TARGET_SIZE} \)` + else + SIZE_DIFF=0 + fi + + echo SIZE_DIFF=$SIZE_DIFF + + echo -e "Estimated content size difference = ${SIZE_DIFF} kB" + + if test ${SIZE_DIFF} -lt $1; then + echo "Size check passed" + return 0 + else + echo "Size check failed" + return 1 + fi + } + +script: + - flake8 + - coverage run run_tests.py + - coverage report -m + # Check the merge size to make sure nothing huge was committed, but only do + # it on one of the branches + - if [ $BUILD_DOCS == 'true' ]; then + size_check 100; + fi; + - | + set -e + if [ $BUILD_DOCS == 'true' ]; then + pip install -r requirements-doc.txt + pushd ../ + git clone https://github.com/scikit-beam/scikit-beam-examples.git + popd + cd doc + bash build_docs.sh + cd .. + fi + +after_success: + - if [ $SUBMIT_CODECOV == 'true' ]; then codecov; fi; + - | + set -e + if [ $TRAVIS_BRANCH == 'master' ] && [ $TRAVIS_PULL_REQUEST == 'false' ] && [ $BUILD_DOCS ]; then + openssl aes-256-cbc -K $encrypted_b4ed46da4197_key -iv $encrypted_b4ed46da4197_iv -in skbeam-docs-deploy.enc -out skbeam-docs-deploy -d + eval `ssh-agent -s` + chmod 600 skbeam-docs-deploy + ssh-add skbeam-docs-deploy + cd doc + bash push_docs.sh + cd .. + fi diff --git a/skbeam_copied/skbeam/core/__init__.py b/skbeam_copied/skbeam/core/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/skbeam_copied/skbeam/core/correlation.py b/skbeam_copied/skbeam/core/correlation.py new file mode 100644 index 00000000..836e5568 --- /dev/null +++ b/skbeam_copied/skbeam/core/correlation.py @@ -0,0 +1,1174 @@ +# coding=utf-8 +# ###################################################################### +# Developed at the NSLS-II, Brookhaven National Laboratory # +# Developed by Sameera K. Abeykoon, February 2014 # +# # +# 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. # +######################################################################## + +""" + +This module is for functions specific to time correlation + +""" +from __future__ import absolute_import, division, print_function +from .utils import multi_tau_lags +from .roi import extract_label_indices +from collections import namedtuple +import numpy as np +from scipy.signal import fftconvolve +# for a convenient status bar +try: + from tqdm import tqdm +except ImportError: + def tqdm(iterator): + return iterator + + +import logging +logger = logging.getLogger(__name__) + + +def _one_time_process(buf, G, past_intensity_norm, future_intensity_norm, + label_array, num_bufs, num_pixels, img_per_level, + level, buf_no, norm, lev_len): + """Reference implementation of the inner loop of multi-tau one time + correlation + + This helper function calculates G, past_intensity_norm and + future_intensity_norm at each level, symmetric normalization is used. + + .. warning :: This modifies inputs in place. + + Parameters + ---------- + buf : array + image data array to use for correlation + G : array + matrix of auto-correlation function without normalizations + past_intensity_norm : array + matrix of past intensity normalizations + future_intensity_norm : array + matrix of future intensity normalizations + label_array : array + labeled array where all nonzero values are ROIs + num_bufs : int, even + number of buffers(channels) + num_pixels : array + number of pixels in certain ROI's + ROI's, dimensions are : [number of ROI's]X1 + img_per_level : array + to track how many images processed in each level + level : int + the current multi-tau level + buf_no : int + the current buffer number + norm : dict + to track bad images + lev_len : array + length of each level + + Notes + ----- + .. math:: + G = + .. math:: + past_intensity_norm = + .. math:: + future_intensity_norm = + """ + img_per_level[level] += 1 + # in multi-tau correlation, the subsequent levels have half as many + # buffers as the first + i_min = num_bufs // 2 if level else 0 + for i in range(i_min, min(img_per_level[level], num_bufs)): + # compute the index into the autocorrelation matrix + t_index = level * num_bufs // 2 + i + delay_no = (buf_no - i) % num_bufs + + # get the images for correlating + past_img = buf[level, delay_no] + future_img = buf[level, buf_no] + + # find the normalization that can work both for bad_images + # and good_images + ind = int(t_index - lev_len[:level].sum()) + normalize = img_per_level[level] - i - norm[level+1][ind] + + # take out the past_ing and future_img created using bad images + # (bad images are converted to np.nan array) + if np.isnan(past_img).any() or np.isnan(future_img).any(): + norm[level + 1][ind] += 1 + else: + for w, arr in zip([past_img*future_img, past_img, future_img], + [G, past_intensity_norm, future_intensity_norm]): + binned = np.bincount(label_array, weights=w)[1:] + arr[t_index] += ((binned / num_pixels - + arr[t_index]) / normalize) + return None # modifies arguments in place! + + +results = namedtuple( + 'correlation_results', + ['g2', 'lag_steps', 'internal_state'] +) + +_internal_state = namedtuple( + 'correlation_state', + ['buf', + 'G', + 'past_intensity', + 'future_intensity', + 'img_per_level', + 'label_array', + 'track_level', + 'cur', + 'pixel_list', + 'num_pixels', + 'lag_steps', + 'norm', + 'lev_len'] +) + +_two_time_internal_state = namedtuple( + 'two_time_correlation_state', + ['buf', + 'img_per_level', + 'label_array', + 'track_level', + 'cur', + 'pixel_list', + 'num_pixels', + 'lag_steps', + 'g2', + 'count_level', + 'current_img_time', + 'time_ind', + 'norm', + 'lev_len'] +) + + +def _init_state_one_time(num_levels, num_bufs, labels): + """Initialize a stateful namedtuple for the generator-based multi-tau + for one time correlation + + Parameters + ---------- + num_levels : int + num_bufs : int + labels : array + Two dimensional labeled array that contains ROI information + + Returns + ------- + internal_state : namedtuple + The namedtuple that contains all the state information that + `lazy_one_time` requires so that it can be used to pick up + processing after it was interrupted + """ + (label_array, pixel_list, num_rois, num_pixels, lag_steps, buf, + img_per_level, track_level, cur, norm, + lev_len) = _validate_and_transform_inputs(num_bufs, num_levels, labels) + + # G holds the un normalized auto- correlation result. We + # accumulate computations into G as the algorithm proceeds. + G = np.zeros(((num_levels + 1) * num_bufs // 2, num_rois), + dtype=np.float64) + # matrix for normalizing G into g2 + past_intensity = np.zeros_like(G) + # matrix for normalizing G into g2 + future_intensity = np.zeros_like(G) + + return _internal_state( + buf, + G, + past_intensity, + future_intensity, + img_per_level, + label_array, + track_level, + cur, + pixel_list, + num_pixels, + lag_steps, + norm, + lev_len, + ) + + +def lazy_one_time(image_iterable, num_levels, num_bufs, labels, + internal_state=None): + """Generator implementation of 1-time multi-tau correlation + + If you do not want multi-tau correlation, set num_levels to 1 and + num_bufs to the number of images you wish to correlate + + Parameters + ---------- + image_iterable : iterable of 2D arrays + num_levels : int + how many generations of downsampling to perform, i.e., the depth of + the binomial tree of averaged frames + num_bufs : int, must be even + maximum lag step to compute in each generation of downsampling + labels : array + Labeled array of the same shape as the image stack. + Each ROI is represented by sequential integers starting at one. For + example, if you have four ROIs, they must be labeled 1, 2, 3, + 4. Background is labeled as 0 + internal_state : namedtuple, optional + internal_state is a bucket for all of the internal state of the + generator. It is part of the `results` object that is yielded from + this generator + + Yields + ------ + namedtuple + A `results` object is yielded after every image has been processed. + This `reults` object contains, in this order: + + - `g2`: the normalized correlation + shape is (len(lag_steps), num_rois) + - `lag_steps`: the times at which the correlation was computed + - `_internal_state`: all of the internal state. Can be passed back in + to `lazy_one_time` as the `internal_state` parameter + + Notes + ----- + The normalized intensity-intensity time-autocorrelation function + is defined as + + .. math:: + + g_2(q, t') = \\frac{ }{^2} + + t' > 0 + + Here, ``I(q, t)`` refers to the scattering strength at the momentum + transfer vector ``q`` in reciprocal space at time ``t``, and the brackets + ``<...>`` refer to averages over time ``t``. The quantity ``t'`` denotes + the delay time + + This implementation is based on published work. [1]_ + + References + ---------- + + .. [1] D. Lumma, L. B. Lurio, S. G. J. Mochrie and M. Sutton, + "Area detector based photon correlation in the regime of + short data batches: Data reduction for dynamic x-ray + scattering," Rev. Sci. Instrum., vol 71, p 3274-3289, 2000. + + """ + + if internal_state is None: + internal_state = _init_state_one_time(num_levels, num_bufs, labels) + # create a shorthand reference to the results and state named tuple + s = internal_state + + # iterate over the images to compute multi-tau correlation + for image in image_iterable: + # Compute the correlations for all higher levels. + level = 0 + + # increment buffer + s.cur[0] = (1 + s.cur[0]) % num_bufs + + # Put the ROI pixels into the ring buffer. + s.buf[0, s.cur[0] - 1] = np.ravel(image)[s.pixel_list] + buf_no = s.cur[0] - 1 + # Compute the correlations between the first level + # (undownsampled) frames. This modifies G, + # past_intensity, future_intensity, + # and img_per_level in place! + _one_time_process(s.buf, s.G, s.past_intensity, s.future_intensity, + s.label_array, num_bufs, s.num_pixels, + s.img_per_level, level, buf_no, s.norm, s.lev_len) + + # check whether the number of levels is one, otherwise + # continue processing the next level + processing = num_levels > 1 + + level = 1 + while processing: + if not s.track_level[level]: + s.track_level[level] = True + processing = False + else: + prev = (1 + (s.cur[level - 1] - 2) % num_bufs) + s.cur[level] = ( + 1 + s.cur[level] % num_bufs) + + s.buf[level, s.cur[level] - 1] = (( + s.buf[level - 1, prev - 1] + + s.buf[level - 1, s.cur[level - 1] - 1]) / 2) + + # make the track_level zero once that level is processed + s.track_level[level] = False + + # call processing_func for each multi-tau level greater + # than one. This is modifying things in place. See comment + # on previous call above. + buf_no = s.cur[level] - 1 + _one_time_process(s.buf, s.G, s.past_intensity, + s.future_intensity, s.label_array, num_bufs, + s.num_pixels, s.img_per_level, level, buf_no, + s.norm, s.lev_len) + level += 1 + + # Checking whether there is next level for processing + processing = level < num_levels + + # If any past intensities are zero, then g2 cannot be normalized at + # those levels. This if/else code block is basically preventing + # divide-by-zero errors. + if len(np.where(s.past_intensity == 0)[0]) != 0: + g_max = np.where(s.past_intensity == 0)[0][0] + else: + g_max = s.past_intensity.shape[0] + + g2 = (s.G[:g_max] / (s.past_intensity[:g_max] * + s.future_intensity[:g_max])) + yield results(g2, s.lag_steps[:g_max], s) + + +def multi_tau_auto_corr(num_levels, num_bufs, labels, images): + """Wraps generator implementation of multi-tau + + Original code(in Yorick) for multi tau auto correlation + author: Mark Sutton + + For parameter description, please reference the docstring for + lazy_one_time. Note that there is an API difference between this function + and `lazy_one_time`. The `images` argument is at the end of this function + signature here for backwards compatibility, but is the first argument in + the `lazy_one_time()` function. The semantics of the variables remain + unchanged. + """ + gen = lazy_one_time(images, num_levels, num_bufs, labels) + for result in gen: + pass + return result.g2, result.lag_steps + + +def auto_corr_scat_factor(lags, beta, relaxation_rate, baseline=1): + """ + This model will provide normalized intensity-intensity time + correlation data to be minimized. + + Parameters + ---------- + lags : array + delay time + beta : float + optical contrast (speckle contrast), a sample-independent + beamline parameter + relaxation_rate : float + relaxation time associated with the samples dynamics. + baseline : float, optional + baseline of one time correlation + equal to one for ergodic samples + + Returns + ------- + g2 : array + normalized intensity-intensity time autocorreltion + + Notes : + ------- + The intensity-intensity autocorrelation g2 is connected to the intermediate + scattering factor(ISF) g1 + + .. math:: + g_2(q, \\tau) = \\beta_1[g_1(q, \\tau)]^{2} + g_\\infty + + For a system undergoing diffusive dynamics, + .. math:: + g_1(q, \\tau) = e^{-\\gamma(q) \\tau} + .. math:: + g_2(q, \\tau) = \\beta_1 e^{-2\\gamma(q) \\tau} + g_\\infty + + These implementation are based on published work. [1]_ + + References + ---------- + .. [1] L. Li, P. Kwasniewski, D. Orsi, L. Wiegart, L. Cristofolini, + C. Caronna and A. Fluerasu, " Photon statistics and speckle + visibility spectroscopy with partially coherent X-rays," + J. Synchrotron Rad. vol 21, p 1288-1295, 2014 + + """ + return beta * np.exp(-2 * relaxation_rate * lags) + baseline + + +def two_time_corr(labels, images, num_frames, num_bufs, num_levels=1): + """Wraps generator implementation of multi-tau two time correlation + + This function computes two-time correlation + Original code : author: Yugang Zhang + + Returns + ------- + results : namedtuple + + For parameter definition, see the docstring for the `lazy_two_time()` + function in this module + """ + gen = lazy_two_time(labels, images, num_frames, num_bufs, num_levels) + for result in gen: + pass + return two_time_state_to_results(result) + + +def lazy_two_time(labels, images, num_frames, num_bufs, num_levels=1, + two_time_internal_state=None): + """Generator implementation of two-time correlation + + If you do not want multi-tau correlation, set num_levels to 1 and + num_bufs to the number of images you wish to correlate + + Multi-tau correlation uses a scheme to achieve long-time correlations + inexpensively by downsampling the data, iteratively combining successive + frames. + + The longest lag time computed is ``num_levels * num_bufs``. + + See Also: ``multi_tau_auto_corr`` for a non-generator implementation + + Parameters + ---------- + labels : array + labeled array of the same shape as the image stack; + each ROI is represented by a distinct label (i.e., integer) + images : iterable of 2D arrays + dimensions are: (rr, cc), iterable of 2D arrays + num_frames : int + number of images to use + default is number of images + num_bufs : int, must be even + maximum lag step to compute in each generation of + downsampling + num_levels : int, optional + how many generations of downsampling to perform, i.e., + the depth of the binomial tree of averaged frames + default is one + + Yields + ------ + namedtuple + A ``results`` object is yielded after every image has been processed. + This `reults` object contains, in this order: + + - ``g2``: the normalized correlation + shape is (num_rois, len(lag_steps), len(lag_steps)) + - ``lag_steps``: the times at which the correlation was computed + - ``_internal_state``: all of the internal state. Can be passed back in + to ``lazy_one_time`` as the ``internal_state`` parameter + + Notes + ----- + The two-time correlation function is defined as + + .. math:: + C(q,t_1,t_2) = \\frac{}{} + + Here, the ensemble averages are performed over many pixels of detector, + all having the same ``q`` value. The average time or age is equal to + ``(t1+t2)/2``, measured by the distance along the ``t1 = t2`` diagonal. + The time difference ``t = |t1 - t2|``, with is distance from the + ``t1 = t2`` diagonal in the perpendicular direction. + In the equilibrium system, the two-time correlation functions depend only + on the time difference ``t``, and hence the two-time correlation contour + lines are parallel. + + [1]_ + + References + ---------- + + .. [1] A. Fluerasu, A. Moussaid, A. Mandsen and A. Schofield, + "Slow dynamics and aging in collodial gels studied by x-ray + photon correlation spectroscopy," Phys. Rev. E., vol 76, p + 010401(1-4), 2007. + + """ + if two_time_internal_state is None: + two_time_internal_state = _init_state_two_time(num_levels, num_bufs, + labels, num_frames) + # create a shorthand reference to the results and state named tuple + s = two_time_internal_state + + for img in images: + s.cur[0] = (1 + s.cur[0]) % num_bufs # increment buffer + + s.count_level[0] = 1 + s.count_level[0] + + # get the current image time + s = s._replace(current_img_time=(s.current_img_time + 1)) + + # Put the image into the ring buffer. + s.buf[0, s.cur[0] - 1] = (np.ravel(img))[s.pixel_list] + + # Compute the two time correlations between the first level + # (undownsampled) frames. two_time and img_per_level in place! + _two_time_process(s.buf, s.g2, s.label_array, num_bufs, + s.num_pixels, s.img_per_level, s.lag_steps, + s.current_img_time, + level=0, buf_no=s.cur[0] - 1) + + # time frame for each level + s.time_ind[0].append(s.current_img_time) + + # check whether the number of levels is one, otherwise + # continue processing the next level + processing = num_levels > 1 + + # Compute the correlations for all higher levels. + level = 1 + while processing: + if not s.track_level[level]: + s.track_level[level] = 1 + processing = False + else: + prev = 1 + (s.cur[level - 1] - 2) % num_bufs + s.cur[level] = 1 + s.cur[level] % num_bufs + s.count_level[level] = 1 + s.count_level[level] + + s.buf[level, s.cur[level] - 1] = (s.buf[level - 1, prev - 1] + + s.buf[level - 1, + s.cur[level - 1] - 1])/2 + + t1_idx = (s.count_level[level] - 1) * 2 + + current_img_time = ((s.time_ind[level - 1])[t1_idx] + + (s.time_ind[level - 1])[t1_idx + 1])/2. + + # time frame for each level + s.time_ind[level].append(current_img_time) + + # make the track_level zero once that level is processed + s.track_level[level] = 0 + + # call the _two_time_process function for each multi-tau level + # for multi-tau levels greater than one + # Again, this is modifying things in place. See comment + # on previous call above. + _two_time_process(s.buf, s.g2, s.label_array, num_bufs, + s.num_pixels, s.img_per_level, s.lag_steps, + current_img_time, + level=level, buf_no=s.cur[level]-1) + level += 1 + + # Checking whether there is next level for processing + processing = level < num_levels + yield s + + +def two_time_state_to_results(state): + """Convert the internal state of the two time generator into usable results + + Parameters + ---------- + state : namedtuple + The internal state that is yielded from `lazy_two_time` + + Returns + ------- + results : namedtuple + A results object that contains the two time correlation results + and the lag steps + """ + for q in range(np.max(state.label_array)): + x0 = (state.g2)[q, :, :] + (state.g2)[q, :, :] = (np.tril(x0) + np.tril(x0).T - + np.diag(np.diag(x0))) + return results(state.g2, state.lag_steps, state) + + +def _two_time_process(buf, g2, label_array, num_bufs, num_pixels, + img_per_level, lag_steps, current_img_time, + level, buf_no): + """ + Parameters + ---------- + buf: array + image data array to use for two time correlation + g2: array + two time correlation matrix + shape (number of labels(ROI), number of frames, number of frames) + label_array: array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, etc. corresponding to the order they are specified + in edges and segments + num_bufs: int, even + number of buffers(channels) + num_pixels : array + number of pixels in certain ROI's + ROI's, dimensions are len(np.unique(label_array)) + img_per_level: array + to track how many images processed in each level + lag_steps : array + delay or lag steps for the multiple tau analysis + shape num_levels + current_img_time : int + the current image number + level : int + the current multi-tau level + buf_no : int + the current buffer number + """ + img_per_level[level] += 1 + + # in multi-tau correlation other than first level all other levels + # have to do the half of the correlation + if level == 0: + i_min = 0 + else: + i_min = num_bufs//2 + + for i in range(i_min, min(img_per_level[level], num_bufs)): + t_index = level*num_bufs//2 + i + + delay_no = (buf_no - i) % num_bufs + + past_img = buf[level, delay_no] + future_img = buf[level, buf_no] + + # get the matrix of correlation function without normalizations + tmp_binned = (np.bincount(label_array, + weights=past_img*future_img)[1:]) + # get the matrix of past intensity normalizations + pi_binned = (np.bincount(label_array, + weights=past_img)[1:]) + + # get the matrix of future intensity normalizations + fi_binned = (np.bincount(label_array, + weights=future_img)[1:]) + + tind1 = (current_img_time - 1) + + tind2 = (current_img_time - lag_steps[t_index] - 1) + + if not isinstance(current_img_time, int): + nshift = 2**(level-1) + for i in range(-nshift+1, nshift+1): + g2[:, int(tind1+i), + int(tind2+i)] = (tmp_binned/(pi_binned * + fi_binned))*num_pixels + else: + g2[:, int(tind1), int(tind2)] = tmp_binned/(pi_binned * fi_binned)*num_pixels + + +def _init_state_two_time(num_levels, num_bufs, labels, num_frames): + """Initialize a stateful namedtuple for two time correlation + + Parameters + ---------- + num_levels : int + num_bufs : int + labels : array + Two dimensional labeled array that contains ROI information + num_frames : int + number of images to use + default is number of images + Returns + ------- + internal_state : namedtuple + The namedtuple that contains all the state information that + `lazy_two_time` requires so that it can be used to pick up processing + after it was interrupted + """ + (label_array, pixel_list, num_rois, num_pixels, lag_steps, + buf, img_per_level, track_level, cur, norm, + lev_len) = _validate_and_transform_inputs(num_bufs, num_levels, labels) + + # to count images in each level + count_level = np.zeros(num_levels, dtype=np.int64) + + # current image time + current_img_time = 0 + + # generate a time frame for each level + time_ind = {key: [] for key in range(num_levels)} + + # two time correlation results (array) + g2 = np.zeros((num_rois, num_frames, num_frames), dtype=np.float64) + + return _two_time_internal_state( + buf, + img_per_level, + label_array, + track_level, + cur, + pixel_list, + num_pixels, + lag_steps, + g2, + count_level, + current_img_time, + time_ind, + norm, + lev_len, + ) + + +def _validate_and_transform_inputs(num_bufs, num_levels, labels): + """ + This is a helper function to validate inputs and create initial state + inputs for both one time and two time correlation + + Parameters + ---------- + num_bufs : int + num_levels : int + labels : array + labeled array of the same shape as the image stack; + each ROI is represented by a distinct label (i.e., integer) + + Returns + ------- + label_array : array + labels of the required region of interests(ROI's) + pixel_list : array + 1D array of indices into the raveled image for all + foreground pixels (labeled nonzero) + e.g., [5, 6, 7, 8, 14, 15, 21, 22] + num_rois : int + number of region of interests (ROI) + num_pixels : array + number of pixels in each ROI + lag_steps : array + the times at which the correlation was computed + buf : array + image data for correlation + img_per_level : array + to track how many images processed in each level + track_level : array + to track processing each level + cur : array + to increment the buffer + norm : dict + to track bad images + lev_len : array + length of each levels + """ + if num_bufs % 2 != 0: + raise ValueError("There must be an even number of `num_bufs`. You " + "provided %s" % num_bufs) + label_array, pixel_list = extract_label_indices(labels) + + # map the indices onto a sequential list of integers starting at 1 + label_mapping = {label: n+1 + for n, label in enumerate(np.unique(label_array))} + # remap the label array to go from 1 -> max(_labels) + for label, n in label_mapping.items(): + label_array[label_array == label] = n + + # number of ROI's + num_rois = len(label_mapping) + + # stash the number of pixels in the mask + num_pixels = np.bincount(label_array)[1:] + + # Convert from num_levels, num_bufs to lag frames. + tot_channels, lag_steps, dict_lag = multi_tau_lags(num_levels, num_bufs) + + # these norm and lev_len will help to find the one time correlation + # normalization norm will updated when there is a bad image + norm = {key: [0] * len(dict_lag[key]) for key in (dict_lag.keys())} + lev_len = np.array([len(dict_lag[i]) for i in (dict_lag.keys())]) + + # Ring buffer, a buffer with periodic boundary conditions. + # Images must be keep for up to maximum delay in buf. + buf = np.zeros((num_levels, num_bufs, len(pixel_list)), + dtype=np.float64) + # to track how many images processed in each level + img_per_level = np.zeros(num_levels, dtype=np.int64) + # to track which levels have already been processed + track_level = np.zeros(num_levels, dtype=bool) + # to increment buffer + cur = np.ones(num_levels, dtype=np.int64) + + return (label_array, pixel_list, num_rois, num_pixels, + lag_steps, buf, img_per_level, track_level, cur, + norm, lev_len) + + +def one_time_from_two_time(two_time_corr): + """ + This will provide the one-time correlation data from two-time + correlation data. + + Parameters + ---------- + two_time_corr : array + matrix of two time correlation + shape (number of labels(ROI's), number of frames, number of frames) + + Returns + ------- + one_time_corr : array + matrix of one time correlation + shape (number of labels(ROI's), number of frames) + """ + + one_time_corr = np.zeros((two_time_corr.shape[0], two_time_corr.shape[2])) + for g in two_time_corr: + for j in range(two_time_corr.shape[2]): + one_time_corr[:, j] = np.trace(g, offset=j)/two_time_corr.shape[2] + return one_time_corr + + +class CrossCorrelator: + ''' + Compute a 1D or 2D cross-correlation on data. + + This uses a mask, which may be binary (array of 0's and 1's), + or a list of non-negative integer id's to compute cross-correlations + separately on. + + The symmetric averaging scheme introduced here is inspired by a paper + from Schätzel, although the implementation is novel in that it + allows for the usage of arbitrary masks. [1]_ + + Examples + -------- + + >> ccorr = CrossCorrelator(mask.shape, mask=mask) + >> # correlated image + >> cimg = cc(img1) + or, mask may may be ids + >> cc = CrossCorrelator(ids) + #(where ids is same shape as img1) + >> cc1 = cc(img1) + >> cc12 = cc(img1, img2) + # if img2 shifts right of img1, point of maximum correlation is shifted + # right from correlation center + + References + ---------- + .. [1] Schätzel, Klaus, Martin Drewel, and Sven Stimac. “Photon + correlation measurements at large lag times: improving + statistical accuracy.” Journal of Modern Optics 35.4 (1988): + 711-718. + + + ''' + def __init__(self, shape, mask=None, normalization=None): + ''' + Prepare the spatial correlator for various regions specified by the + id's in the image. + + Parameters + ---------- + shape : 1 or 2-tuple + The shape of the incoming images or curves. May specify 1D or + 2D shapes by inputting a 1 or 2-tuple + + mask : 1D or 2D np.ndarray of int, optional + Each non-zero integer represents unique bin. Zero integers are + assumed to be ignored regions. If None, creates a mask with + all points set to 1 + + normalization: string or list of strings, optional + These specify the normalization and may be any of the + following: + 'regular' : divide by pixel number + 'symavg' : use symmetric averaging + Defaults to ['regular'] normalization + ''' + if normalization is None: + normalization = ['regular'] + elif not isinstance(normalization, list): + normalization = list([normalization]) + self.normalization = normalization + + if mask is None: + mask = np.ones(shape) + + self.shape = shape + self.ndim = len(shape) + if self.ndim == 1: + mask = mask.reshape((1, shape[0])) + + if self.ndim > 2: + raise ValueError("Dimensions other than 1 or 2 not supported") + + # A terse nomenclature was chosen for ease of reading code. + # Here is the Nomenclature (for both init and call sections): + # 1. pi, pj: the list of yindex (rows), xindex (cols) into original + # image + # 2. pis, pjs : a subselection of the above + # 3. ppii, ppjj : the list of yindex (rows), xindex (cols) into + # smaller sub images + # 4. ppiis, ppjjs : a subselection of the above + # 5. bind : the bin id per pixel (defined by mask) + # 6. idpos : a list of pointers into the position arrays for each id + + # initialize subregions information for the correlation + # first find indices of subregions and sort them by subregion id + pi, pj = np.where(mask != 0) + # typing not necessary, but it may be better to be explicit + # (especially if we cythonize in the future) + bind = mask[pi, pj].astype(np.uint32) + + # sort them so that we can easily index into subregions with one number + order = np.argsort(bind) + bind = bind[order] + pi = pi[order] + pj = pj[order] + + # make array of pointers into these subregions + idpos = 0 + inds_tmp = 1 + np.where(np.not_equal(bind[1:], bind[:-1]))[0] + idpos = np.append(idpos, inds_tmp) + idpos = np.append(idpos, len(bind)).astype(np.uint32) + self.idpos = idpos + + self.nids = len(idpos[:-1]) + + # finds the shape of the subregions + shapes = np.zeros((self.nids, 4), dtype=np.uint32) + for i in range(self.nids): + index_start, index_end = idpos[i], idpos[i+1] + row_min = pi[index_start:index_end].min() + row_max = pi[index_start:index_end].max() + col_min = pj[index_start:index_end].min() + col_max = pj[index_start:index_end].max() + shapes[i] = [row_min, row_max, col_min, col_max] + + # make indices for subregions arrays and their shapes + ppii = pi.copy() + ppjj = pj.copy() + # subtract the minimum index, this trick + # will then index into some other array properly + for i in range(self.nids): + index_start, index_stop = idpos[i], idpos[i+1] + ppii[index_start:index_stop] -= shapes[i, 0] + ppjj[index_start:index_stop] -= shapes[i, 2] + # now save to object + self.ppii = ppii + self.ppjj = ppjj + self.pi = pi + self.pj = pj + + # redefine shapes to be the shapes of the subregions + # make shapes be for regions + shapes = np.array(1 + (np.diff(shapes, axis=-1)[:, [0, 2]])) + self.shapes = shapes + # We now have two sets of positions of the subregions (pi,pj) in + # subregion and (pii,pjj) in images. pos is a pointer such that + # (pos[i]:pos[i+1]) is the indices in the position arrays of subregion + # i. + + # Making a list of arrays holding the masks for each id. Ideally, mask + # is binary so this is one element to quickly index original images + self.submasks = list() + self.centers = list() + # the positions of each axes of each correlation + self.positions = list() + self.maskcorrs = list() + # regions where the correlations are not zero + self.pxlst_maskcorrs = list() + + # basically saving bunch of mask related stuff like indexing etc, just + # to save some time when actually computing the cross correlations + for i in range(self.nids): + submask = np.zeros(self.shapes[i, :]) + index_start, index_stop = idpos[i], idpos[i+1] + ppiis = ppii[index_start:index_stop] + ppjjs = ppjj[index_start:index_stop] + submask[ppiis, ppjjs] = 1 + self.submasks.append(submask) + + maskcorr = _cross_corr(submask) + # choose some small value to threshold + maskcorr *= maskcorr > .5 + # This following line was originally placed two lines below. It was moved here in order + # to avoid runtime warnings during the execution of '>' (zeros are replaced by 'nan' values) + self.pxlst_maskcorrs.append(maskcorr > 0) + maskcorr[np.where(maskcorr == 0)] = np.nan + self.maskcorrs.append(maskcorr) + # centers are shape//2 as performed by fftshift + center = np.array(maskcorr.shape)//2 + self.centers.append(np.array(maskcorr.shape)//2) + if self.ndim == 1: + self.positions.append(np.arange(maskcorr.shape[1]) - center[1]) + elif self.ndim == 2: + self.positions.append([np.arange(maskcorr.shape[0]) - + center[0], + np.arange(maskcorr.shape[1]) - + center[1]]) + + if self.nids == 1: + self.positions = self.positions[0] + self.centers = self.centers[0] + + def __call__(self, img1, img2=None, normalization=None): + ''' Run the cross correlation on an image/curve or against two + images/curves + + Parameters + ---------- + img1 : 1D or 2D np.ndarray + The image (or curve) to run the cross correlation on + + img2 : 1D or 2D np.ndarray + If not set to None, run cross correlation of this image (or + curve) against img1. Default is None. + + normalization : string or list of strings + normalization types. If not set, use internally saved + normalization parameters + + Returns + ------- + ccorrs : 1d or 2d np.ndarray + An image of the correlation. The zero correlation is + located at shape//2 where shape is the 1 or 2-tuple + shape of the array + + ''' + if normalization is None: + normalization = self.normalization + + if img1.shape != self.shape: + raise ValueError("Image not expected shape." + + "Got {}, ".format(img1.shape) + + "expected {}".format(self.shape) + ) + # reshape for 1D case + if self.ndim == 1: + img1 = img1.reshape((1, self.shape[0])) + + if img2 is None: + self_correlation = True + else: + self_correlation = False + if img2.shape != self.shape: + raise ValueError("Second image not expected shape. " + + "Got {}".format(img2.shape) + + " expected {}".format(self.shape)) + if self.ndim == 1: + img2 = img2.reshape((1, self.shape[0])) + + ccorrs = list() + rngiter = tqdm(range(self.nids)) + + idpos = self.idpos + for i in rngiter: + index_start, index_stop = idpos[i], idpos[i+1] + ppiis = self.ppii[index_start:index_stop] + ppjjs = self.ppjj[index_start:index_stop] + pis = self.pi[index_start:index_stop] + pjs = self.pj[index_start:index_stop] + tmpimg = np.zeros(self.shapes[i, :]) + tmpimg[ppiis, ppjjs] = img1[pis, pjs] + + if not self_correlation: + tmpimg2 = np.zeros_like(tmpimg) + tmpimg2[ppiis, ppjjs] = img2[pis, pjs] + tmpimg2[ppiis, ppjjs] = img2[pis, pjs] + + if self_correlation: + ccorr = _cross_corr(tmpimg) + else: + ccorr = _cross_corr(tmpimg, tmpimg2) + + # Note, in this code, non-overlapping regions will now get np.nan + # also, for sym averaging, if Icorr*Icorr2==0, then we also get + # np.nan + if 'symavg' in normalization: + # do symmetric averaging + Icorr = _cross_corr(tmpimg * self.submasks[i], + self.submasks[i]) + if self_correlation: + Icorr2 = _cross_corr(self.submasks[i], tmpimg * + self.submasks[i]) + else: + Icorr2 = _cross_corr(self.submasks[i], tmpimg2 * + self.submasks[i]) + ccorr *= self.maskcorrs[i]/Icorr/Icorr2 + + if 'regular' in normalization: + if self_correlation: + ccorr /= self.maskcorrs[i] * \ + np.average(tmpimg[ppiis, ppjjs])**2 + else: + ccorr /= self.maskcorrs[i] *\ + np.average(tmpimg[ppiis, ppjjs]) *\ + np.average(tmpimg2[ppiis, ppjjs]) + if self.ndim == 1: + ccorr = ccorr.reshape(-1) + ccorrs.append(ccorr) + + if len(ccorrs) == 1: + ccorrs = ccorrs[0] + + return ccorrs + + +def _cross_corr(img1, img2=None): + ''' Compute the cross correlation of one (or two) images. + + Parameters + ---------- + img1 : np.ndarray + the image or curve to cross correlate + + img2 : 1d or 2d np.ndarray, optional + If set, cross correlate img1 against img2. A shift of img2 + to the right of img1 will lead to a shift of the point of + highest correlation to the right. + Default is set to None + ''' + ndim = img1.ndim + + if img2 is None: + img2 = img1 + + if img1.shape != img2.shape: + errorstr = "Image shapes don't match. " + errorstr += "(img1 : {}; img2 : {})"\ + .format(img1.shape, img2.shape) + raise ValueError(errorstr) + + # need to reverse indices for second image + # fftconvolve(A,B) = FFT^(-1)(FFT(A)*FFT(B)) + # but need FFT^(-1)(FFT(A(x))*conj(FFT(B(x)))) = FFT^(-1)(A(x)*B(-x)) + reverse_index = [slice(None, None, -1) for i in range(ndim)] + imgc = fftconvolve(img1, img2[tuple(reverse_index)], mode='full') + + return imgc diff --git a/skbeam_copied/skbeam/core/mask.py b/skbeam_copied/skbeam/core/mask.py new file mode 100644 index 00000000..7c6f9569 --- /dev/null +++ b/skbeam_copied/skbeam/core/mask.py @@ -0,0 +1,187 @@ +# ###################################################################### +# 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, January 2016 # +# # +# 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. # +######################################################################## + +""" +This module is for functions specific to mask or threshold an image +basically to clean images +""" + +from __future__ import absolute_import, division, print_function +import numpy as np + +import logging +import scipy.stats as sts +logger = logging.getLogger(__name__) + + +def bad_to_nan_gen(images, bad): + """ + Convert the images marked as "bad" in `bad` by their index in + images into a np.nan array + + Parameters + ---------- + images : iterable + Iterable of 2-D arrays + bad : list + List of integer indices into the `images` parameter that mark those + images as "bad". + + Yields + ------ + img : array + if image is bad it will convert to np.nan array otherwise no + change to the array + """ + ret_val = None + for n, im in enumerate(images): + if n in bad: + if ret_val is None: + ret_val = np.empty(im.shape) + ret_val[:] = np.nan + yield ret_val + else: + yield im + + +def threshold(images, threshold, mask=None): + """ + This generator sets all pixels whose value is greater than `threshold` + to 0 and yields the thresholded images out + + Parameters + ---------- + images : iterable + Iterable of 2-D arrays + threshold : float + threshold value to remove the hot spots in the image + mask : array + array with values above the threshold marked as 0 and values + below marked as 1. + shape is (num_columns, num_rows) of the image, optional None + + Yields + ------- + mask : array + array with values above the threshold marked as 0 and values + below marked as 1. + shape is (num_columns, num_rows) of the image + """ + if mask is None: + mask = np.ones_like(images[0]) + for im in images: + bad_pixels = np.where(im >= threshold) + if len(bad_pixels[0]) != 0: + mask[bad_pixels] = 0 + yield mask + + +def margin(img_shape, edge_size): + """ + Mask the edge of an image + + Parameters + ----------- + img_shape: tuple + The shape of the image + edge_size: int + Number of pixels to mask from the edge + + Returns + -------- + 2darray: + The mask array, bad pixels are 0 + """ + mask = np.ones(img_shape, dtype=bool) + mask[edge_size:-edge_size, edge_size:-edge_size] = 0. + return ~mask + + +def binned_outlier(img, r, alpha, bins, mask=None): + """ + Generates a mask by identifying outlier pixels in bins and masks any + pixels which have a value greater or less than alpha * std away from the + mean + + Parameters + ---------- + img: 2darray + The image + r: 2darray + The array which maps pixels to bins + alpha: float or tuple or, 1darray + Then number of acceptable standard deviations, if tuple then we use + a linear distribution of alphas from alpha[0] to alpha[1], if array + then we just use that as the distribution of alphas + bins: list + The bin edges + mask: 1darray, bool + A starting flattened mask + + Returns + -------- + 2darray: + The mask + """ + + if mask is None: + working_mask = np.ones(img.shape).astype(bool) + else: + working_mask = np.copy(mask).astype(bool) + if working_mask.shape != img.shape: + working_mask = working_mask.reshape(img.shape) + msk_img = img[working_mask] + msk_r = r[working_mask] + + int_r = np.digitize(r, bins[:-1], True) - 1 + # integration + mean = sts.binned_statistic(msk_r, msk_img, bins=bins, + statistic='mean')[0] + std = sts.binned_statistic(msk_r, msk_img, bins=bins, + statistic=np.std)[0] + if type(alpha) is tuple: + alpha = np.linspace(alpha[0], alpha[1], len(std)) + threshold = alpha * std + lower = mean - threshold + upper = mean + threshold + + # single out the too low and too high pixels + working_mask *= img > lower[int_r] + working_mask *= img < upper[int_r] + + return working_mask.astype(bool) diff --git a/skbeam_copied/skbeam/core/roi.py b/skbeam_copied/skbeam/core/roi.py new file mode 100644 index 00000000..0e4960f9 --- /dev/null +++ b/skbeam_copied/skbeam/core/roi.py @@ -0,0 +1,759 @@ +#! encoding: utf-8 +# ###################################################################### +# 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. # +######################################################################## + +""" +This module contain convenience methods to generate ROI labeled arrays for +simple shapes such as rectangles and concentric circles. +""" +from __future__ import absolute_import, division, print_function + +import collections +import scipy.ndimage.measurements as ndim +from skimage.draw import line +from skimage import img_as_float, feature, color, draw +from skimage.measure import ransac, CircleModel +import numpy as np +from . import utils +import logging + +logger = logging.getLogger(__name__) + + +def rectangles(coords, shape): + """ + This function wil provide the indices array for rectangle region of + interests. + + Parameters + ---------- + coords : iterable + coordinates of the upper-left corner and width and height of each + rectangle: e.g., [(x, y, w, h), (x, y, w, h)] + + shape : tuple + Image shape which is used to determine the maximum extent of output + pixel coordinates. Order is (rr, cc). + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are specified + in coords. Order is (rr, cc). + + """ + + labels_grid = np.zeros(shape, dtype=np.int64) + + for i, (col_coor, row_coor, col_val, row_val) in enumerate(coords): + + left, right = np.max([col_coor, 0]), np.min([col_coor + col_val, + shape[0]]) + top, bottom = np.max([row_coor, 0]), np.min([row_coor + row_val, + shape[1]]) + + slc1 = slice(left, right) + slc2 = slice(top, bottom) + + if np.any(labels_grid[slc1, slc2]): + raise ValueError("overlapping ROIs") + + # assign a different scalar for each roi + labels_grid[slc1, slc2] = (i + 1) + + return labels_grid + + +def rings(edges, center, shape): + """ + Draw annual (ring-shaped) shaped regions of interest. + + Each ring will be labeled with an integer. Regions outside any ring will + be filled with zeros. + + Parameters + ---------- + edges: list + giving the inner and outer radius of each ring + e.g., [(1, 2), (11, 12), (21, 22)] + center: tuple + point in image where r=0; may be a float giving subpixel precision. + Order is (rr, cc). + shape: tuple + Image shape which is used to determine the maximum extent of output + pixel coordinates. Order is (rr, cc). + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are specified + in edges. + """ + edges = np.atleast_2d(np.asarray(edges)).ravel() + if not 0 == len(edges) % 2: + raise ValueError("edges should have an even number of elements, " + "giving inner, outer radii for each ring") + if not np.all(np.diff(edges) >= 0): + raise ValueError("edges are expected to be monotonically increasing, " + "giving inner and outer radii of each ring from " + "r=0 outward") + r_coord = utils.radial_grid(center, shape).ravel() + return _make_roi(r_coord, edges, shape) + + +def ring_edges(inner_radius, width, spacing=0, num_rings=None): + """ Calculate the inner and outer radius of a set of rings. + + The number of rings, their widths, and any spacing between rings can be + specified. They can be uniform or varied. + + Parameters + ---------- + inner_radius : float + inner radius of the inner-most ring + + width : float or list of floats + ring thickness + If a float, all rings will have the same thickness. + + spacing : float or list of floats, optional + margin between rings, 0 by default + If a float, all rings will have the same spacing. If a list, + the length of the list must be one less than the number of + rings. + + num_rings : int, optional + number of rings + Required if width and spacing are not lists and number + cannot thereby be inferred. If it is given and can also be + inferred, input is checked for consistency. + + Returns + ------- + edges : array + inner and outer radius for each ring + + Example + ------- + # Make two rings starting at r=1px, each 5px wide + >>> ring_edges(inner_radius=1, width=5, num_rings=2) + [(1, 6), (6, 11)] + # Make three rings of different widths and spacings. + # Since the width and spacings are given individually, the number of + # rings here is simply inferred. + >>> ring_edges(inner_radius=1, width=(5, 4, 3), spacing=(1, 2)) + [(1, 6), (7, 11), (13, 16)] + """ + # All of this input validation merely checks that width, spacing, and + # num_rings are self-consistent and complete. + width_is_list = isinstance(width, collections.abc.Iterable) + spacing_is_list = isinstance(spacing, collections.abc.Iterable) + if (width_is_list and spacing_is_list): + if len(width) != len(spacing) - 1: + raise ValueError("List of spacings must be one less than list " + "of widths.") + if num_rings is None: + try: + num_rings = len(width) + except TypeError: + try: + num_rings = len(spacing) + 1 + except TypeError: + raise ValueError("Since width and spacing are constant, " + "num_rings cannot be inferred and must be " + "specified.") + else: + if width_is_list: + if num_rings != len(width): + raise ValueError("num_rings does not match width list") + if spacing_is_list: + if num_rings-1 != len(spacing): + raise ValueError("num_rings does not match spacing list") + + # Now regularlize the input. + if not width_is_list: + width = np.ones(num_rings) * width + if not spacing_is_list: + spacing = np.ones(num_rings - 1) * spacing + + # The inner radius is the first "spacing." + all_spacings = np.insert(spacing, 0, inner_radius) + steps = np.array([all_spacings, width]).T.ravel() + edges = np.cumsum(steps).reshape(-1, 2) + + return edges + + +def segmented_rings(edges, segments, center, shape, offset_angle=0): + """ + Parameters + ---------- + edges : array + inner and outer radius for each ring + + segments : int or list + number of pie slices or list of angles in radians + That is, 8 produces eight equal-sized angular segments, + whereas a list can be used to produce segments of unequal size. + + center : tuple + point in image where r=0; may be a float giving subpixel precision. + Order is (rr, cc). + + shape: tuple + Image shape which is used to determine the maximum extent of output + pixel coordinates. Order is (rr, cc). + + angle_offset : float or array, optional + offset in radians from offset_angle=0 along the positive X axis + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are specified + in edges and segments + + See Also + -------- + ring_edges : Calculate the inner and outer radius of a set of rings. + + """ + edges = np.asarray(edges).ravel() + if not 0 == len(edges) % 2: + raise ValueError("edges should have an even number of elements, " + "giving inner, outer radii for each ring") + if not np.all(np.diff(edges) >= 0): + raise ValueError("edges are expected to be monotonically increasing, " + "giving inner and outer radii of each ring from " + "r=0 outward") + + agrid = utils.angle_grid(center, shape) + + agrid[agrid < 0] = 2*np.pi + agrid[agrid < 0] + + segments_is_list = isinstance(segments, collections.abc.Iterable) + if segments_is_list: + segments = np.asarray(segments) + offset_angle + else: + # N equal segments requires N+1 bin edges spanning 0 to 2pi. + segments = np.linspace(0, 2*np.pi, num=1+segments, endpoint=True) + segments += offset_angle + + # the indices of the bins(angles) to which each value in input + # array(angle_grid) belongs. + ind_grid = (np.digitize(np.ravel(agrid), segments, + right=False)).reshape(shape) + + label_array = np.zeros(shape, dtype=np.int64) + # radius grid for the image_shape + rgrid = utils.radial_grid(center, shape) + + # assign indices value according to angles then rings + len_segments = len(segments) + for i in range(len(edges) // 2): + indices = (edges[2*i] <= rgrid) & (rgrid < edges[2*i + 1]) + # Combine "segment #" and "ring #" to get unique label for each. + label_array[indices] = ind_grid[indices] + (len_segments - 1) * i + + return label_array + + +def roi_max_counts(images_sets, label_array): + """ + Return the brightest pixel in any ROI in any image in the image set. + + Parameters + ---------- + images_sets : array + iterable of 4D arrays + shapes is: (len(images_sets), ) + + label_array : array + labeled array; 0 is background. + Each ROI is represented by a distinct label (i.e., integer). + + Returns + ------- + max_counts : int + maximum pixel counts + """ + max_cts = 0 + for img_set in images_sets: + for img in img_set: + max_cts = max(max_cts, ndim.maximum(img, label_array)) + return max_cts + + +def roi_pixel_values(image, labels, index=None): + """ + This will provide intensities of the ROI's of the labeled array + according to the pixel list + eg: intensities of the rings of the labeled array + + Parameters + ---------- + image : array + image data dimensions are: (rr, cc) + + labels : array + labeled array; 0 is background. + Each ROI is represented by a distinct label (i.e., integer). + + index_list : list, optional + labels list + eg: 5 ROI's + index = [1, 2, 3, 4, 5] + + Returns + ------- + roi_pix : list + intensities of the ROI's of the labeled array according + to the pixel list + + """ + if labels.shape != image.shape: + raise ValueError("Shape of the image data should be equal to" + " shape of the labeled array") + if index is None: + index = np.arange(1, np.max(labels) + 1) + + roi_pix = [] + for n in index: + roi_pix.append(image[labels == n]) + return roi_pix, index + + +def mean_intensity(images, labeled_array, index=None): + """Compute the mean intensity for each ROI in the image list + + Parameters + ---------- + images : list + List of images + labeled_array : array + labeled array; 0 is background. + Each ROI is represented by a nonzero integer. It is not required that + the ROI labels are contiguous + index : int, list, optional + The ROI's to use. If None, this function will extract averages for all + ROIs + + Returns + ------- + mean_intensity : array + The mean intensity of each ROI for all `images` + Dimensions: + + - len(mean_intensity) == len(index) + - len(mean_intensity[0]) == len(images) + index : list + The labels for each element of the `mean_intensity` list + """ + if labeled_array.shape != images[0].shape[0:]: + raise ValueError( + "`images` shape (%s) needs to be equal to the labeled_array shape" + "(%s)" % (images[0].shape, labeled_array.shape)) + # handle various input for `index` + if index is None: + index = list(np.unique(labeled_array)) + index.remove(0) + try: + len(index) + except TypeError: + index = [index] + # pre-allocate an array for performance + # might be able to use list comprehension to make this faster + mean_intensity = np.zeros((images.shape[0], len(index))) + for n, img in enumerate(images): + # use a mean that is mask-aware + mean_intensity[n] = ndim.mean(img, labeled_array, index=index) + return mean_intensity, index + + +def circular_average(image, calibrated_center, threshold=0, nx=100, + pixel_size=(1, 1), min_x=None, max_x=None, mask=None): + """Circular average of the the image data + The circular average is also known as the radial integration + + Parameters + ---------- + image : array + Image to compute the average as a function of radius + calibrated_center : tuple + The center of the image in pixel units + argument order should be (row, col) + threshold : int, optional + Ignore counts below `threshold` + default is zero + nx : int, optional + number of bins in x + defaults is 100 bins + pixel_size : tuple, optional + The size of a pixel (in a real unit, like mm). + argument order should be (pixel_height, pixel_width) + default is (1, 1) + min_x : float, optional number of pixels + Left edge of first bin defaults to minimum value of x + max_x : float, optional number of pixels + Right edge of last bin defaults to maximum value of x + mask : mask for 2D data. Assumes 1 is non masked and 0 masked. + None defaults to no mask. + + Returns + ------- + bin_centers : array + The center of each bin in R. shape is (nx, ) + ring_averages : array + Radial average of the image. shape is (nx, ). + + See Also + -------- + bad_to_nan_gen : Create a mask with np.nan entries + bin_grid : Bin and integrate an image, given the radial array of pixels + Useful for nonlinear spacing (Ewald curvature) + """ + radial_val = utils.radial_grid(calibrated_center, image.shape, pixel_size) + + if mask is not None: + w = np.where(mask == 1) + radial_val = radial_val[w] + image = image[w] + + bin_edges, sums, counts = utils.bin_1D(np.ravel(radial_val), + np.ravel(image), nx, + min_x=min_x, + max_x=max_x) + th_mask = counts > threshold + ring_averages = sums[th_mask] / counts[th_mask] + + bin_centers = utils.bin_edges_to_centers(bin_edges)[th_mask] + + return bin_centers, ring_averages + + +def kymograph(images, labels, num): + """ + This function will provide data for graphical representation of pixels + variation over time for required ROI. + + Parameters + ---------- + images : array + Image stack. dimensions are: (num_img, num_rows, num_cols) + labels : array + labeled array; 0 is background. Each ROI is represented by an integer + num : int + The ROI to turn into a kymograph + + Returns + ------- + kymograph : array + data for graphical representation of pixels variation over time + for required ROI + + """ + kymo = [] + for n, img in enumerate(images): + kymo.append((roi_pixel_values(img, labels == num)[0])) + + return np.vstack(kymo) + + +def extract_label_indices(labels): + """ + This will find the label's required region of interests (roi's), + number of roi's count the number of pixels in each roi's and pixels + list for the required roi's. + + Parameters + ---------- + labels : array + labeled array; 0 is background. + Each ROI is represented by a distinct label (i.e., integer). + + Returns + ------- + label_mask : array + 1D array labeling each foreground pixel + e.g., [1, 1, 1, 1, 2, 2, 1, 1] + + indices : array + 1D array of indices into the raveled image for all + foreground pixels (labeled nonzero) + e.g., [5, 6, 7, 8, 14, 15, 21, 22] + """ + img_dim = labels.shape + + # TODO Make this tighter. + w = np.where(np.ravel(labels) > 0) + grid = np.indices((img_dim[0], img_dim[1])) + pixel_list = np.ravel((grid[0] * img_dim[1] + grid[1]))[w] + + # discard the zeros + label_mask = labels[labels > 0] + + return label_mask, pixel_list + + +def _make_roi(coords, edges, shape): + """ Helper function to create ring rois and bar rois + + Parameters + ---------- + coords : array + shape is image shape + edges : list + List of tuples of inner (left or top) and outer (right or bottom) + edges of each roi. + e.g., edges=[(1, 2), (11, 12), (21, 22)] + shape : tuple + Shape of the image in which to create the ROIs + e.g., shape=(512, 512) + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are + specified in `edges`. + Has shape=`image shape` + """ + label_array = np.digitize(coords, edges, right=False) + # Even elements of label_array are in the space between rings. + label_array = (np.where(label_array % 2 != 0, label_array, 0) + 1) // 2 + return label_array.reshape(shape) + + +def bar(edges, shape, horizontal=True, values=None): + """Draw bars defined by `edges` from one edge to the other of `image_shape` + + Bars will be horizontal or vertical depending on the value of `horizontal` + + Parameters + ---------- + edges : list + List of tuples of inner (left or top) and outer (right or bottom) + edges of each bar. + e.g., edges=[(1, 2), (11, 12), (21, 22)] + shape : tuple + Shape of the image in which to create the ROIs + e.g., shape=(512, 512) + horizontal : bool, optional + True: Make horizontal bars + False: Make vertical bars + Defaults to True + values : array, optional + image pixels co-ordinates + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are + specified in `edges`. + Has shape=`image shape` + + Note + ---- + The primary use case is in GISAXS. + """ + edges = np.atleast_2d(np.asarray(edges)).ravel() + if not 0 == len(edges) % 2: + raise ValueError("edges should have an even number of elements, " + "giving inner, outer edge value for each bar") + if not np.all(np.diff(edges) >= 0): + raise ValueError("edges are expected to be monotonically increasing, " + "giving inner and outer radii of each bar from " + "r=0 outward") + if values is None: + values = np.repeat(range(shape[0]), shape[1]) + if not horizontal: + values = np.tile(range(shape[1]), shape[0]) + + return _make_roi(values, edges, shape) + + +def box(shape, v_edges, h_edges=None, h_values=None, v_values=None): + """Draw box shaped rois when the horizontal and vertical edges + are provided. + + Parameters + ---------- + shape : tuple + Shape of the image in which to create the ROIs + e.g., shape=(512, 512) + v_edges : list + giving the inner and outer edges of each vertical bar + e.g., [(1, 2), (11, 12), (21, 22)] + h_edges : list, optional + giving the inner and outer edges of each horizontal bar + e.g., [(1, 2), (11, 12), (21, 22)] + h_values : array, optional + image pixels co-ordinates in horizontal direction + shape has to be image shape + v_values : array, optional + image pixels co-ordinates in vertical direction + shape has to be image shape + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are specified + in edges. + + Note + ---- + To draw boxes according to the image pixels co-ordinates has to provide + both h_values and v_values. The primary use case is in GISAXS. + e.g., v_values=gisaxs_qy, h_values=gisaxs_qx + + """ + if h_edges is None: + h_edges = v_edges + + if h_values is None and v_values is None: + v_values, h_values = np.mgrid[:shape[0], :shape[1]] + elif h_values.shape != v_values.shape: + raise ValueError("Shape of the h_values array should be equal to" + " shape of the v_values array") + for edges in (h_edges, v_edges): + edges = np.atleast_2d(np.asarray(edges)).ravel() + if not 0 == len(edges) % 2: + raise ValueError("edges should have an even number of elements, " + "giving inner, outer edges for each roi") + coords = [] + for h in h_edges: + for v in v_edges: + coords.append((h[0], v[0], h[1]-h[0], v[1] - v[0])) + + return rectangles(coords, v_values.shape) + + +def lines(end_points, shape): + """ + Parameters + ---------- + end_points : iterable + coordinates of the starting point and the ending point of each + line: e.g., [(start_x, start_y, end_x, end_y), (x1, y1, x2, y2)] + shape : tuple + Image shape which is used to determine the maximum extent of output + pixel coordinates. Order is (rr, cc). + + Returns + ------- + label_array : array + Elements not inside any ROI are zero; elements inside each + ROI are 1, 2, 3, corresponding to the order they are specified + in coords. Order is (rr, cc). + + """ + label_array = np.zeros(shape, dtype=np.int64) + label = 0 + for points in end_points: + if len(points) != 4: + raise ValueError("end points should have four number of" + " elements, giving starting co-ordinates," + " ending co-ordinates for each line") + rr, cc = line(np.max([points[0], 0]), np.max([points[1], 0]), + np.min([points[2], shape[0]-1]), + np.min([points[3], shape[1]-1])) + label += 1 + label_array[rr, cc] = label + return label_array + + +def auto_find_center_rings(avg_img, sigma=1, no_rings=4, min_samples=3, + residual_threshold=1, max_trials=1000): + """This will find the center of the speckle pattern and the radii of the + most intense rings. + + Parameters + ---------- + avg_img : 2D array + shape of the image + sigma : float, optional + Standard deviation of the Gaussian filter. + no_rings : int, optional + number of rings + min_sample : int, optional + The minimum number of data points to fit a model to. + residual_threshold : float, optional + Maximum distance for a data point to be classified as an inlier. + max_trials : int, optional + Maximum number of iterations for random sample selection. + + Returns + ------- + center : tuple + center co-ordinates of the speckle pattern + image : 2D array + Indices of pixels that belong to the rings, + directly index into an array + radii : list + values of the radii of the rings + + Note + ---- + scikit-image ransac + method(http://www.imagexd.org/tutorial/lessons/1_ransac.html) is used to + automatically find the center and the most intense rings. + """ + + image = img_as_float(color.rgb2gray(avg_img)) + edges = feature.canny(image, sigma) + coords = np.column_stack(np.nonzero(edges)) + edge_pts_xy = coords[:, ::-1] + radii = [] + + for i in range(no_rings): + model_robust, inliers = ransac(edge_pts_xy, CircleModel, min_samples, + residual_threshold, + max_trials=max_trials) + if i == 0: + center = int(model_robust.params[0]), int(model_robust.params[1]) + radii.append(model_robust.params[2]) + + rr, cc = draw.circle_perimeter(center[1], center[0], + int(model_robust.params[2]), + shape=image.shape) + image[rr, cc] = i + 1 + edge_pts_xy = edge_pts_xy[~inliers] + + return center, image, radii diff --git a/skbeam_copied/skbeam/core/tests/__init__.py b/skbeam_copied/skbeam/core/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/skbeam_copied/skbeam/core/tests/test_correlation.py b/skbeam_copied/skbeam/core/tests/test_correlation.py new file mode 100644 index 00000000..46c5100d --- /dev/null +++ b/skbeam_copied/skbeam/core/tests/test_correlation.py @@ -0,0 +1,385 @@ +# ###################################################################### +# 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_equal, assert_array_almost_equal +import pytest + +import skbeam.core.utils as utils +from skbeam.core.correlation import (multi_tau_auto_corr, + auto_corr_scat_factor, + lazy_one_time, + lazy_two_time, two_time_corr, + two_time_state_to_results, + one_time_from_two_time, + CrossCorrelator) +from skbeam.core.mask import bad_to_nan_gen +from skbeam.core.roi import ring_edges, segmented_rings + + +logger = logging.getLogger(__name__) + + +def setup(): + global num_levels, num_bufs, xdim, ydim, stack_size, img_stack, rois + num_levels = 6 + num_bufs = 4 # must be even + xdim = 256 + ydim = 512 + stack_size = 100 + img_stack = np.random.randint(1, 3, (stack_size, xdim, ydim)) + rois = np.zeros_like(img_stack[0]) + # make sure that the ROIs can be any integers greater than 1. + # They do not have to start at 1 and be continuous + rois[0:xdim//10, 0:ydim//10] = 5 + rois[xdim//10:xdim//5, ydim//10:ydim//5] = 3 + + +def test_lazy_vs_original(): + setup() + # run the correlation on the full stack + full_gen_one = lazy_one_time( + img_stack, num_levels, num_bufs, rois) + for gen_state_one in full_gen_one: + pass + + g2, lag_steps = multi_tau_auto_corr(num_levels, num_bufs, + rois, img_stack) + + assert np.all(g2 == gen_state_one.g2) + assert np.all(lag_steps == gen_state_one.lag_steps) + + full_gen_two = lazy_two_time(rois, img_stack, stack_size, + num_bufs, num_levels) + for gen_state_two in full_gen_two: + pass + final_gen_result_two = two_time_state_to_results(gen_state_two) + + two_time = two_time_corr(rois, img_stack, stack_size, + num_bufs, num_levels) + + assert np.all(two_time[0] == final_gen_result_two.g2) + assert np.all(two_time[1] == final_gen_result_two.lag_steps) + + +def test_lazy_two_time(): + setup() + # run the correlation on the full stack + full_gen = lazy_two_time(rois, img_stack, stack_size, + stack_size, 1) + for full_state in full_gen: + pass + final_result = two_time_state_to_results(full_state) + + # make sure we have essentially zero correlation in the images, + # since they are random integers + assert np.average(final_result.g2-1) < 0.01 + + # run the correlation on the first half + gen_first_half = lazy_two_time(rois, img_stack[:stack_size//2], stack_size, + num_bufs=stack_size, num_levels=1) + for first_half_state in gen_first_half: + pass + # run the correlation on the second half by passing in the state from the + # first half + gen_second_half = lazy_two_time(rois, img_stack[stack_size//2:], + stack_size, num_bufs=stack_size, + num_levels=1, + two_time_internal_state=first_half_state) + + for second_half_state in gen_second_half: + pass + result = two_time_state_to_results(second_half_state) + + assert np.all(full_state.g2 == result.g2) + + +def test_lazy_one_time(): + setup() + # run the correlation on the full stack + full_gen = lazy_one_time(img_stack, num_levels, num_bufs, rois) + for full_result in full_gen: + pass + + # make sure we have essentially zero correlation in the images, + # since they are random integers + assert np.average(full_result.g2-1) < 0.01 + + # run the correlation on the first half + gen_first_half = lazy_one_time( + img_stack[:stack_size//2], num_levels, num_bufs, rois) + for first_half_result in gen_first_half: + pass + # run the correlation on the second half by passing in the state from the + # first half + gen_second_half = lazy_one_time( + img_stack[stack_size//2:], num_levels, num_bufs, rois, + internal_state=first_half_result.internal_state + ) + + for second_half_result in gen_second_half: + pass + + assert np.all(full_result.g2 == + second_half_result.g2) + + +def test_two_time_corr(): + setup() + y = [] + for i in range(50): + y.append(img_stack[0]) + two_time = two_time_corr(rois, np.asarray(y), 50, + num_bufs=50, num_levels=1) + assert np.all(two_time[0]) + + # check the number of buffers are even + with pytest.raises(ValueError): + two_time_corr(rois, np.asarray(y), 50, num_bufs=25, num_levels=1) + + +def test_auto_corr_scat_factor(): + num_levels, num_bufs = 3, 4 + tot_channels, lags, dict_lags = utils.multi_tau_lags(num_levels, num_bufs) + beta = 0.5 + relaxation_rate = 10.0 + baseline = 1.0 + + g2 = auto_corr_scat_factor(lags, beta, relaxation_rate, baseline) + + assert_array_almost_equal(g2, np.array([1.5, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0]), decimal=8) + + +def test_bad_images(): + setup() + g2, lag_steps = multi_tau_auto_corr(4, num_bufs, + rois, img_stack) + # introduce bad images + bad_img_list = [3, 21, 35, 48] + # convert each bad image to np.nan array + images = bad_to_nan_gen(img_stack, bad_img_list) + + # then use new images (including bad images) + g2_n, lag_steps_n = multi_tau_auto_corr(4, num_bufs, + rois, images) + + assert_array_almost_equal(g2[:, 0], g2_n[:, 0], decimal=3) + assert_array_almost_equal(g2[:, 1], g2_n[:, 1], decimal=3) + + +def test_one_time_from_two_time(): + num_lev = 1 + num_buf = 10 # must be even + x_dim = 10 + y_dim = 10 + stack = 10 + imgs = np.random.randint(1, 3, (stack, x_dim, y_dim)) + roi = np.zeros_like(imgs[0]) + # make sure that the ROIs can be any integers greater than 1. + # They do not have to start at 1 and be continuous + roi[0:x_dim//10, 0:y_dim//10] = 5 + roi[x_dim//10:x_dim//5, y_dim//10:y_dim//5] = 3 + + g2, lag_steps, _state = two_time_corr(roi, imgs, stack, + num_buf, num_lev) + + one_time = one_time_from_two_time(g2) + assert_array_almost_equal(one_time[0, :], np.array([1.0, 0.9, 0.8, 0.7, + 0.6, 0.5, 0.4, 0.3, + 0.2, 0.1])) + + +@pytest.mark.skipif(int(np.__version__.split('.')[1]) > 14, reason="Test is numerically unstable") +def test_CrossCorrelator1d(): + ''' Test the 1d version of the cross correlator with these methods: + -method='regular', no mask + -method='regular', masked + -method='symavg', no mask + -method='symavg', masked + ''' + np.random.seed(123) + # test 1D data + sigma = .1 + Npoints = 100 + x = np.linspace(-10, 10, Npoints) + + sigma = .2 + # purposely have sparsely filled values (with lots of zeros) + peak_positions = (np.random.random(10)-.5)*20 + y = np.zeros_like(x) + for peak_position in peak_positions: + y += np.exp(-(x-peak_position)**2/2./sigma**2) + + mask_1D = np.ones_like(y) + mask_1D[10:20] = 0 + mask_1D[60:90] = 0 + mask_1D[111:137] = 0 + mask_1D[211:237] = 0 + mask_1D[411:537] = 0 + + mask_1D *= mask_1D[::-1] + + cc1D = CrossCorrelator(mask_1D.shape) + cc1D_symavg = CrossCorrelator(mask_1D.shape, normalization='symavg') + cc1D_masked = CrossCorrelator(mask_1D.shape, mask=mask_1D) + cc1D_masked_symavg = CrossCorrelator(mask_1D.shape, mask=mask_1D, + normalization='symavg') + + assert_equal(cc1D.nids, 1) + + ycorr_1D = cc1D(y) + ycorr_1D_masked = cc1D_masked(y*mask_1D) + ycorr_1D_symavg = cc1D_symavg(y) + ycorr_1D_masked_symavg = cc1D_masked_symavg(y*mask_1D) + + assert_array_almost_equal(ycorr_1D[::20], + np.array([-1.155123e-14, 6.750373e-03, + 6.221636e-01, 7.105527e-01, + 1.187275e+00, 2.984563e+00, + 1.092725e+00, 1.198341e+00, + 1.045922e-01, 5.451511e-06])) + assert_array_almost_equal(ycorr_1D_masked[::20], + np.array([-5.172377e-16, np.nan, 7.481473e-01, + 6.066887e-02, 4.470989e-04, + 2.330335e+00, np.nan, 7.109758e-01, + np.nan, 2.275846e-14])) + + assert_array_almost_equal(ycorr_1D_symavg[::20], + np.array([-5.3002753, 1.54268227, 0.86220476, + 0.57715207, 0.86503802, 2.94383202, + 0.7587901, 0.99763715, 0.16800951, + 1.23506293])) + + assert_array_almost_equal(ycorr_1D_masked_symavg[::20][:-1], + np.array([-5.30027530e-01, np.nan, + 1.99940257e+00, 7.33127871e-02, + 1.00000000e+00, 2.15887870e+00, + np.nan, 9.12832602e-01, + np.nan])) + + +def test_CrossCorrelator2d(): + ''' Test the 2D case of the cross correlator. + With non-binary labels. + ''' + np.random.seed(123) + # test 2D data + Npoints2 = 10 + x2 = np.linspace(-10, 10, Npoints2) + X, Y = np.meshgrid(x2, x2) + Z = np.random.random((Npoints2, Npoints2)) + + np.random.seed(123) + sigma = .2 + # purposely have sparsely filled values (with lots of zeros) + # place peaks in random positions + peak_positions = (np.random.random((2, 10))-.5)*20 + for peak_position in peak_positions: + Z += np.exp(-((X - peak_position[0])**2 + + (Y - peak_position[1])**2)/2./sigma**2) + + mask_2D = np.ones_like(Z) + mask_2D[1:2, 1:2] = 0 + mask_2D[7:9, 4:6] = 0 + mask_2D[1:2, 9:] = 0 + + # Compute with segmented rings + edges = ring_edges(1, 3, num_rings=2) + segments = 5 + x0, y0 = np.array(mask_2D.shape)//2 + + maskids = segmented_rings(edges, segments, (y0, x0), mask_2D.shape) + + cc2D_ids = CrossCorrelator(mask_2D.shape, mask=maskids) + cc2D_ids_symavg = CrossCorrelator(mask_2D.shape, mask=maskids, + normalization='symavg') + + # 10 ids + assert_equal(cc2D_ids.nids, 10) + + ycorr_ids_2D = cc2D_ids(Z) + ycorr_ids_2D_symavg = cc2D_ids_symavg(Z) + index = 0 + ycorr_ids_2D[index][ycorr_ids_2D[index].shape[0]//2] + assert_array_almost_equal(ycorr_ids_2D[index] + [ycorr_ids_2D[index].shape[0]//2], + np.array([1.22195059, 1.08685771, + 1.43246508, 1.08685771, 1.22195059 + ]) + ) + + index = 1 + ycorr_ids_2D[index][ycorr_ids_2D[index].shape[0]//2] + assert_array_almost_equal(ycorr_ids_2D[index] + [ycorr_ids_2D[index].shape[0]//2], + np.array([1.24324268, 0.80748997, + 1.35790022, 0.80748997, 1.24324268 + ]) + ) + + index = 0 + ycorr_ids_2D_symavg[index][ycorr_ids_2D[index].shape[0]//2] + assert_array_almost_equal(ycorr_ids_2D_symavg[index] + [ycorr_ids_2D[index].shape[0]//2], + np.array([0.84532695, 1.16405848, 1.43246508, + 1.16405848, 0.84532695]) + ) + + index = 1 + ycorr_ids_2D_symavg[index][ycorr_ids_2D[index].shape[0]//2] + assert_array_almost_equal(ycorr_ids_2D_symavg[index] + [ycorr_ids_2D[index].shape[0]//2], + np.array([0.94823482, 0.8629459, 1.35790022, + 0.8629459, 0.94823482]) + ) + + +def test_CrossCorrelator_badinputs(): + with pytest.raises(ValueError): + CrossCorrelator((1, 1, 1)) + + with pytest.raises(ValueError): + cc = CrossCorrelator((10, 10)) + a = np.ones((10, 11)) + cc(a) + + with pytest.raises(ValueError): + cc = CrossCorrelator((10, 10)) + a = np.ones((10, 10)) + a2 = np.ones((10, 11)) + cc(a, a2) diff --git a/skbeam_copied/skbeam/core/tests/test_mask.py b/skbeam_copied/skbeam/core/tests/test_mask.py new file mode 100644 index 00000000..a9258a52 --- /dev/null +++ b/skbeam_copied/skbeam/core/tests/test_mask.py @@ -0,0 +1,139 @@ +# ###################################################################### +# 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_equal + +import skbeam.core.mask as mask + +logger = logging.getLogger(__name__) + + +def test_threshold_mask(): + xdim = 10 + ydim = 10 + stack_size = 10 + img_stack = np.random.randint(1, 3, (stack_size, xdim, ydim)) + + img_stack[0][0, 1] = 100 + img_stack[0][9, 1] = 98 + img_stack[6][8, 8] = 75 + img_stack[7][6, 6] = 80 + + th = mask.threshold(img_stack, 75) + + for final in th: + pass + + y = np.ones_like(img_stack[0]) + y[0, 1] = 0 + y[9, 1] = 0 + y[8, 8] = 0 + y[6, 6] = 0 + + assert_array_equal(final, y) + + +def test_bad_to_nan_gen(): + xdim = 2 + ydim = 2 + stack_size = 5 + img_stack = np.random.randint(1, 3, (stack_size, xdim, ydim)) + + bad_list = [1, 3] + + img = mask.bad_to_nan_gen(img_stack, bad_list) + y = [] + for im in img: + y.append(im) + + assert np.isnan(np.asarray(y)[1]).all() + assert np.isnan(np.asarray(y)[3]).all() + assert not np.isnan(np.asarray(y)[4]).all() + + +def test_margin(): + size = (10, 10) + edge = 1 + mask1 = mask.margin(size, edge) + mask2 = np.zeros(size) + mask2[:, :edge] = 1 + mask2[:, -edge:] = 1 + mask2[:edge, :] = 1 + mask2[-edge:, :] = 1 + mask2 = mask2.astype(bool) + assert_array_equal(mask1, ~mask2) + + +def test_ring_blur_mask(): + from skbeam.core import recip + g = recip.geo.Geometry( + detector='Perkin', pixel1=.0002, pixel2=.0002, + dist=.23, + poni1=.209, poni2=.207, + # rot1=.0128, rot2=-.015, rot3=-5.2e-8, + wavelength=1.43e-11 + ) + r = g.rArray((2048, 2048)) + # make some sample data + Z = 100 * np.cos(50 * r) ** 2 + 150 + + np.random.seed(10) + pixels = [] + for i in range(0, 100): + a, b = np.random.randint(low=0, high=2048), \ + np.random.randint(low=0, high=2048) + if np.random.random() > .5: + # Add some hot pixels + Z[a, b] = np.random.randint(low=200, high=255) + else: + # and dead pixels + Z[a, b] = np.random.randint(low=0, high=10) + pixels.append((a, b)) + pixel_size = [getattr(g, k) for k in ['pixel1', 'pixel2']] + rres = np.hypot(*pixel_size) + bins = np.arange(np.min(r) - rres/2., np.max(r) + rres / 2., rres) + msk = mask.binned_outlier(Z, r, (3., 3), bins, mask=None) + a = set(zip(*np.nonzero(~msk))) + b = set(pixels) + a_not_in_b = a - b + b_not_in_a = b - a + + # We have not over masked 10% of the number of bad pixels + assert len(a_not_in_b) / len(b) < .1 + # Make certain that we have masked over 90% of the bad pixels + assert len(b_not_in_a) / len(b) < .1 diff --git a/skbeam_copied/skbeam/core/tests/test_roi.py b/skbeam_copied/skbeam/core/tests/test_roi.py new file mode 100644 index 00000000..54ccc227 --- /dev/null +++ b/skbeam_copied/skbeam/core/tests/test_roi.py @@ -0,0 +1,410 @@ +# ###################################################################### +# 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 skbeam.core import roi +from skbeam.core import utils +import itertools +from skimage import morphology + +from numpy.testing import (assert_array_equal, assert_array_almost_equal, + assert_equal, assert_raises, assert_almost_equal) + +logger = logging.getLogger(__name__) + + +def test_rectangles(): + shape = (15, 26) + roi_data = np.array(([2, 2, 6, 3], [6, 7, 8, 5], [8, 18, 5, 10]), + dtype=np.int64) + + all_roi_inds = roi.rectangles(roi_data, shape) + + roi_inds, pixel_list = roi.extract_label_indices(all_roi_inds) + + ty = np.zeros(shape).ravel() + ty[pixel_list] = roi_inds + + re_mesh = ty.reshape(*shape) + for i, (col_coor, row_coor, col_val, row_val) in enumerate(roi_data, 0): + ind_co = np.column_stack(np.where(re_mesh == i + 1)) + + left, right = np.max([col_coor, 0]), np.min([col_coor + col_val, + shape[0]]) + top, bottom = np.max([row_coor, 0]), np.min([row_coor + row_val, + shape[1]]) + assert_almost_equal(left, ind_co[0][0]) + assert_almost_equal(right-1, ind_co[-1][0]) + assert_almost_equal(top, ind_co[0][1]) + assert_almost_equal(bottom-1, ind_co[-1][-1]) + + +def test_rings(): + center = (100., 100.) + img_dim = (200, 205) + first_q = 10. + delta_q = 5. + num_rings = 7 # number of Q rings + one_step_q = 5.0 + step_q = [2.5, 3.0, 5.8] + + # test when there is same spacing between rings + edges = roi.ring_edges(first_q, width=delta_q, spacing=one_step_q, + num_rings=num_rings) + print("edges there is same spacing between rings ", edges) + label_array = roi.rings(edges, center, img_dim) + print("label_array there is same spacing between rings", label_array) + label_mask, pixel_list = roi.extract_label_indices(label_array) + # number of pixels per ROI + num_pixels = np.bincount(label_mask, minlength=(np.max(label_mask)+1)) + num_pixels = num_pixels[1:] + + # test when there is same spacing between rings + edges = roi.ring_edges(first_q, width=delta_q, spacing=2.5, + num_rings=num_rings) + print("edges there is same spacing between rings ", edges) + label_array = roi.rings(edges, center, img_dim) + print("label_array there is same spacing between rings", label_array) + label_mask, pixel_list = roi.extract_label_indices(label_array) + # number of pixels per ROI + num_pixels = np.bincount(label_mask, minlength=(np.max(label_mask)+1)) + num_pixels = num_pixels[1:] + + # test when there is different spacing between rings + edges = roi.ring_edges(first_q, width=delta_q, spacing=step_q, + num_rings=4) + print("edges when there is different spacing between rings", edges) + label_array = roi.rings(edges, center, img_dim) + print("label_array there is different spacing between rings", label_array) + label_mask, pixel_list = roi.extract_label_indices(label_array) + # number of pixels per ROI + num_pixels = np.bincount(label_mask, minlength=(np.max(label_mask)+1)) + num_pixels = num_pixels[1:] + + # test when there is no spacing between rings + edges = roi.ring_edges(first_q, width=delta_q, num_rings=num_rings) + print("edges", edges) + label_array = roi.rings(edges, center, img_dim) + print("label_array", label_array) + label_mask, pixel_list = roi.extract_label_indices(label_array) + # number of pixels per ROI + num_pixels = np.bincount(label_mask, minlength=(np.max(label_mask)+1)) + num_pixels = num_pixels[1:] + + # Did we draw the right number of rings? + print(np.unique(label_array)) + actual_num_rings = len(np.unique(label_array)) - 1 + assert_equal(actual_num_rings, num_rings) + + # Does each ring have more pixels than the last, being larger? + ring_areas = np.bincount(label_array.ravel())[1:] + area_comparison = np.diff(ring_areas) + print(area_comparison) + areas_monotonically_increasing = np.all(area_comparison > 0) + assert areas_monotonically_increasing + + # Test various illegal inputs + assert_raises(ValueError, + lambda: roi.ring_edges(1, 2)) # need num_rings + # width incompatible with num_rings + assert_raises(ValueError, + lambda: roi.ring_edges(1, [1, 2, 3], num_rings=2)) + # too few spacings + assert_raises(ValueError, + lambda: roi.ring_edges(1, [1, 2, 3], [1])) + # too many spacings + assert_raises(ValueError, + lambda: roi.ring_edges(1, [1, 2, 3], [1, 2, 3])) + # num_rings conflicts with width, spacing + assert_raises(ValueError, + lambda: roi.ring_edges(1, [1, 2, 3], [1, 2], 5)) + w_edges = [[5, 7], [1, 2]] + assert_raises(ValueError, roi.rings, w_edges, center=(4, 4), + shape=(20, 20)) + + +def _helper_check(pixel_list, inds, num_pix, edges, center, + img_dim, num_qs): + # recreate the indices using pixel_list and inds values + ty = np.zeros(img_dim).ravel() + ty[pixel_list] = inds + # get the grid values from the center + grid_values = utils.radial_grid(img_dim, center) + + # get the indices into a grid + zero_grid = np.zeros((img_dim[0], img_dim[1])) + for r in range(num_qs): + vl = (edges[r][0] <= grid_values) & (grid_values < edges[r][1]) + zero_grid[vl] = r + 1 + + # check the num_pixels + num_pixels = [] + for r in range(num_qs): + num_pixels.append(int((np.histogramdd(np.ravel(grid_values), bins=1, + range=[[edges[r][0], + (edges[r][1] - + 0.000001)]]))[0][0])) + assert_array_equal(num_pix, num_pixels) + + +def test_segmented_rings(): + center = (75, 75) + img_dim = (150, 140) + first_q = 5 + delta_q = 5 + num_rings = 4 # number of Q rings + slicing = 4 + + edges = roi.ring_edges(first_q, width=delta_q, spacing=4, + num_rings=num_rings) + print("edges", edges) + + label_array = roi.segmented_rings(edges, slicing, center, + img_dim, offset_angle=0) + print("label_array for segmented_rings", label_array) + + # Did we draw the right number of ROIs? + label_list = np.unique(label_array.ravel()) + actual_num_labels = len(label_list) - 1 + num_labels = num_rings * slicing + assert_equal(actual_num_labels, num_labels) + + # Did we draw the right ROIs? (1-16 with some zeros around too) + assert_array_equal(label_list, np.arange(num_labels + 1)) + + # A brittle test to make sure the exactly number of pixels per label + # is never accidentally changed: + # number of pixels per ROI + num_pixels = np.bincount(label_array.ravel()) + expected_num_pixels = [18372, 59, 59, 59, 59, 129, 129, 129, + 129, 200, 200, 200, 200, 269, 269, 269, 269] + assert_array_equal(num_pixels, expected_num_pixels) + + +def test_roi_pixel_values(): + images = morphology.diamond(8) + # width incompatible with num_rings + + label_array = np.zeros((256, 256)) + + # different shapes for the images and labels + assert_raises(ValueError, + lambda: roi.roi_pixel_values(images, label_array)) + # create a label mask + center = (8., 8.) + inner_radius = 2. + width = 1 + spacing = 1 + edges = roi.ring_edges(inner_radius, width, spacing, num_rings=5) + rings = roi.rings(edges, center, images.shape) + + intensity_data, index = roi.roi_pixel_values(images, rings) + assert_array_equal(intensity_data[0], ([1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1])) + assert_array_equal([1, 2, 3, 4, 5], index) + + +def test_roi_max_counts(): + img_stack1 = np.random.randint(0, 60, size=(50, ) + (50, 50)) + img_stack2 = np.random.randint(0, 60, size=(100, ) + (50, 50)) + + img_stack1[0][20, 20] = 60 + + samples = (img_stack1, img_stack2) + + label_array = np.zeros((img_stack1[0].shape)) + + label_array[img_stack1[0] < 20] = 1 + label_array[img_stack1[0] > 40] = 2 + + assert_array_equal(60, roi.roi_max_counts(samples, label_array)) + + +def test_static_test_sets(): + images1 = [] + for i in range(10): + int_array = np.tril(i*np.ones(50)) + int_array[int_array == 0] = i*100 + images1.append(int_array) + + images2 = [] + for i in range(20): + int_array = np.triu(i*np.ones(50)) + int_array[int_array == 0] = i*100 + images2.append(int_array) + + samples = {'sample1': np.asarray(images1), 'sample2': np.asarray(images2)} + + roi_data = np.array(([2, 30, 12, 15], [40, 20, 15, 10]), dtype=np.int64) + + label_array = roi.rectangles(roi_data, shape=(50, 50)) + + # get the mean intensities of image sets given as a dictionary + roi_data = [] + for k, v in sorted(samples.items()): + intensity, index_list = roi.mean_intensity(v, label_array) + roi_data.append(intensity) + + return_values = [roi_data[0][:, 0], roi_data[0][:, 1], + roi_data[1][:, 0], roi_data[1][:, 1], ] + expected_values = [ + np.asarray([float(x) for x in range(0, 1000, 100)]), + np.asarray([float(x) for x in range(0, 10, 1)]), + np.asarray([float(x) for x in range(0, 20, 1)]), + np.asarray([float(x) for x in range(0, 2000, 100)]) + ] + err_msg = ['roi%s of sample%s is incorrect' % (i, j) + for i, j in itertools.product((1, 2), (1, 2))] + for returned, expected, err in zip(return_values, + expected_values, err_msg): + assert_array_equal(returned, expected, + err_msg=err, verbose=True) + + +def test_circular_average(): + image = np.zeros((12, 12)) + calib_center = (5, 5) + inner_radius = 1 + + edges = roi.ring_edges(inner_radius, width=1, spacing=1, num_rings=2) + labels = roi.rings(edges, calib_center, image.shape) + image[labels == 1] = 10 + image[labels == 2] = 10 + bin_cen, ring_avg = roi.circular_average(image, calib_center, nx=6) + + assert_array_almost_equal(bin_cen, [0.70710678, 2.12132034, + 3.53553391, 4.94974747, 6.36396103, + 7.77817459], decimal=6) + assert_array_almost_equal(ring_avg, [8., 2.5, 5.55555556, 0., + 0., 0.], decimal=6) + + bin_cen1, ring_avg1 = roi.circular_average(image, calib_center, min_x=0, + max_x=10, nx=None) + assert_array_almost_equal(bin_cen1, [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, + 7.5, 8.5]) + mask = np.ones_like(image) + mask[4:6, 2:3] = 0 + + bin_cen_masked, ring_avg_masked = roi.circular_average(image, + calib_center, + min_x=0, + max_x=10, + nx=6, + mask=mask) + + assert_array_almost_equal(bin_cen_masked, [0.83333333, 2.5, 4.16666667, + 5.83333333, 7.5, 9.16666667]) + + assert_array_almost_equal(ring_avg_masked, [8.88888889, 3.84615385, 2.5, + 0., 0., 0.]) + + +def test_kymograph(): + calib_center = (25, 25) + inner_radius = 5 + + edges = roi.ring_edges(inner_radius, width=2, num_rings=1) + labels = roi.rings(edges, calib_center, (50, 50)) + + images = [] + num_images = 100 + for i in range(num_images): + int_array = i*np.ones(labels.shape) + images.append(int_array) + + kymograph_data = roi.kymograph(np.asarray(images), labels, num=1) + # make sure the the return array has the expected dimensions + expected_shape = (num_images, np.sum(labels[labels == 1])) + assert kymograph_data.shape[0] == expected_shape[0] + assert kymograph_data.shape[1] == expected_shape[1] + # make sure we got one element from each image + assert np.all(kymograph_data[:, 0] == np.arange(num_images)) + # given the input data, every row of kymograph_data should be the same + # number + for row in kymograph_data: + assert np.all(row == row[0]) + + +def test_bars_boxes(): + edges = [[3, 4], [5, 7]] + shape = (10, 10) + h_label_array = roi.bar(edges, shape) + v_label_array = roi.bar(edges, shape, horizontal=False) + w_edges = [[5, 7], [1, 2]] + + assert_array_equal(h_label_array[3, :], np.ones((10,), dtype=np.int)) + assert_array_equal(v_label_array[:, 3], np.ones((10,), dtype=np.int)) + assert_array_equal(np.unique(h_label_array)[1:], np.array([1, 2])) + assert_array_equal(np.unique(v_label_array)[1:], np.array([1, 2])) + assert_raises(ValueError, roi.bar, w_edges, shape) + + box_array = roi.box(shape, edges) + + assert_array_equal(box_array[3, 3], 1) + assert_array_equal(box_array[5, 5], 4) + assert_array_equal(np.unique(box_array)[1:], np.array([1, 2, 3, 4])) + + n_edges = edges = [[3, 4], [5, 7], [8]] + h_values, v_values = np.mgrid[0:10, 0:10] + h_val, v_val = np.mgrid[0:20, 0:20] + + assert_raises(ValueError, roi.box, shape, n_edges) + assert_raises(ValueError, roi.box, shape, edges, h_values=h_val, + v_values=v_values) + + +def test_lines(): + points = ([30, 45, 50, 256], [56, 60, 80, 150]) + shape = (256, 150) + label_array = roi.lines(points, shape) + assert_array_equal(np.array([1, 2]), np.unique(label_array)[1:]) + + assert_raises(ValueError, roi.lines, + ([10, 12, 30], [30, 45, 50, 256]), shape) + + +def test_auto_find_center_rings(): + + x = np.linspace(-5, 5, 200) + X, Y = np.meshgrid(x, x) + image = 100 * np.cos(np.sqrt(x**2 + Y**2))**2 + 50 + + center, image, radii = roi.auto_find_center_rings(image, sigma=20, + no_rings=2) + + assert_equal((99, 99), center) + assert_array_equal(41., np.round(radii[0])) diff --git a/skbeam_copied/skbeam/core/tests/test_utils.py b/skbeam_copied/skbeam/core/tests/test_utils.py new file mode 100644 index 00000000..48971d44 --- /dev/null +++ b/skbeam_copied/skbeam/core/tests/test_utils.py @@ -0,0 +1,572 @@ +# ###################################################################### +# 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 six +import numpy as np +import sys +import pytest + +import numpy.testing as npt +from numpy.testing import (assert_array_equal, assert_array_almost_equal, + assert_equal, assert_almost_equal) + +import skbeam.core.utils as core + +import logging + +try: + from pyFAI.geometry import Geometry + pf = True +except ImportError: + pf = False + +logger = logging.getLogger(__name__) + + +def test_bin_1D(): + # set up simple data + x = np.linspace(0, 1, 100) + y = np.arange(100) + nx = 10 + # make call + edges, val, count = core.bin_1D(x, y, nx) + # check that values are as expected + assert_array_almost_equal(edges, + np.linspace(0, 1, nx + 1, endpoint=True)) + assert_array_almost_equal(val, + np.sum(y.reshape(nx, -1), axis=1)) + assert_array_equal(count, + np.ones(nx) * 10) + + +def test_bin_1D_2(): + """ + Test for appropriate default value handling + """ + # set up simple data + x = np.linspace(0, 1, 100) + y = np.arange(100) + nx = core._defaults["bins"] + min_x = None + max_x = None + # make call + edges, val, count = core.bin_1D(x=x, y=y, nx=nx, min_x=min_x, max_x=max_x) + # check that values are as expected + assert_array_almost_equal(edges, + np.linspace(0, 1, nx + 1, endpoint=True)) + assert_array_almost_equal(val, + np.sum(y.reshape(nx, -1), axis=1)) + assert_array_equal(count, + np.ones(nx)) + + +def test_bin_1D_limits(): + # set up simple data + x = np.linspace(0, 1, 100) + y = np.arange(100) + nx = 10 + min_x, max_x = .25, .75 + # make call + edges, val, count = core.bin_1D(x, y, nx, min_x, max_x) + # check that values are as expected + assert_array_almost_equal(edges, + np.linspace(min_x, max_x, nx + 1, endpoint=True)) + assert_array_almost_equal(val, + np.sum(y[25:75].reshape(nx, -1), axis=1)) + assert_array_equal(count, + np.ones(nx) * 5) + + +def _bin_edges_helper(p_dict): + bin_edges = core.bin_edges(**p_dict) + assert_almost_equal(0, np.ptp(np.diff(bin_edges))) + if 'nbins' in p_dict: + nbins = p_dict['nbins'] + assert_equal(nbins + 1, len(bin_edges)) + if 'step' in p_dict: + step = p_dict['step'] + assert_almost_equal(step, np.diff(bin_edges)) + if 'range_max' in p_dict: + range_max = p_dict['range_max'] + assert np.all(bin_edges <= range_max) + if 'range_min' in p_dict: + range_min = p_dict['range_min'] + assert np.all(bin_edges >= range_min) + if 'range_max' in p_dict and 'step' in p_dict: + step = p_dict['step'] + range_max = p_dict['range_max'] + assert (range_max - bin_edges[-1]) < step + + +def _bin_edges_exceptions(param_dict): + with pytest.raises(ValueError): + core.bin_edges(**param_dict) + + +param_test_bin_edges = ['range_min', 'range_max', 'step', 'nbins'] + + +@pytest.mark.parametrize("drop_key", param_test_bin_edges) +def test_bin_edges(drop_key): + test_dict = {'range_min': 1.234, + 'range_max': 5.678, + 'nbins': 42, + 'step': np.pi / 10} + test_dict.pop(drop_key) + _bin_edges_helper(test_dict) + + +param_test_bin_edges_exceptions = [ + # no entries + {}, + # 4 entries + {'range_min': 1.234, 'range_max': 5.678, 'nbins': 42, + 'step': np.pi / 10}, + # 2 entries + {'range_min': 1.234, 'step': np.pi / 10}, + # 1 entry + {'range_min': 1.234, }, + # max < min + {'range_max': 1.234, 'range_min': 5.678, 'step': np.pi / 10}, + # step > max - min + {'range_min': 1.234, 'range_max': 5.678, 'step': np.pi * 10}, + # nbins == 0 + {'range_min': 1.234, 'range_max': 5.678, 'nbins': 0}] + + +@pytest.mark.parametrize("fail_dict", param_test_bin_edges_exceptions) +def test_bin_edges_exceptions(fail_dict): + _bin_edges_exceptions(fail_dict) + + +def test_grid3d(): + size = 10 + q_max = np.array([1.0, 1.0, 1.0]) + q_min = np.array([-1.0, -1.0, -1.0]) + dqn = np.array([size, size, size]) + param_dict = {'nx': dqn[0], + 'ny': dqn[1], + 'nz': dqn[2], + 'xmin': q_min[0], + 'ymin': q_min[1], + 'zmin': q_min[2], + 'xmax': q_max[0], + 'ymax': q_max[1], + 'zmax': q_max[2]} + # slice tricks + # this make a list of slices, the imaginary value in the + # step is interpreted as meaning 'this many values' + slc = [slice(_min + (_max - _min) / (s * 2), + _max - (_max - _min) / (s * 2), + 1j * s) + for _min, _max, s in zip(q_min, q_max, dqn)] + # use the numpy slice magic to make X, Y, Z these are dense meshes with + # points in the center of each bin + X, Y, Z = np.mgrid[slc] + + # make and ravel the image data (which is all ones) + II = np.ones_like(X).ravel() + + # make input data (Nx3 + data = np.array([np.ravel(X), + np.ravel(Y), + np.ravel(Z)]).T + + (mean, occupancy, + std_err, bounds) = core.grid3d(data, II, **param_dict) + + # check the values are as expected + npt.assert_array_equal(mean.ravel(), II) + npt.assert_array_equal(occupancy, np.ones_like(occupancy)) + npt.assert_array_equal(std_err, 0) + + +def test_process_grid_std_err(): + size = 10 + q_max = np.array([1.0, 1.0, 1.0]) + q_min = np.array([-1.0, -1.0, -1.0]) + dqn = np.array([size, size, size]) + param_dict = {'nx': dqn[0], + 'ny': dqn[1], + 'nz': dqn[2], + 'xmin': q_min[0], + 'ymin': q_min[1], + 'zmin': q_min[2], + 'xmax': q_max[0], + 'ymax': q_max[1], + 'zmax': q_max[2]} + # slice tricks + # this make a list of slices, the imaginary value in the + # step is interpreted as meaning 'this many values' + slc = [slice(_min + (_max - _min) / (s * 2), + _max - (_max - _min) / (s * 2), + 1j * s) + for _min, _max, s in zip(q_min, q_max, dqn)] + # use the numpy slice magic to make X, Y, Z these are dense meshes with + # points in the center of each bin + X, Y, Z = np.mgrid[slc] + + # make and ravel the image data (which is all ones) + I = np.hstack([j * np.ones_like(X).ravel() for j in range(1, 101)]) # noqa: E741 + + # make input data (N*5x3) + data = np.vstack([np.tile(_, 100) + for _ in (np.ravel(X), np.ravel(Y), np.ravel(Z))]).T + (mean, occupancy, + std_err, bounds) = core.grid3d(data, I, **param_dict) + + # check the values are as expected + npt.assert_array_equal(mean, + np.ones_like(X) * np.mean(np.arange(1, 101))) + npt.assert_array_equal(occupancy, np.ones_like(occupancy) * 100) + # need to convert std -> ste (standard error) + # according to wikipedia ste = std/sqrt(n) + npt.assert_array_almost_equal(std_err, + (np.ones_like(occupancy) * + np.std(np.arange(1, 101)) / np.sqrt(100))) + + +def test_bin_edge2center(): + test_edges = np.arange(11) + centers = core.bin_edges_to_centers(test_edges) + assert_array_almost_equal(.5, centers % 1) + assert_equal(10, len(centers)) + + +def test_small_verbosedict(): + expected_string = ("You tried to access the key 'b' " + "which does not exist. " + "The extant keys are: ['a']") + dd = core.verbosedict() + dd['a'] = 1 + assert_equal(dd['a'], 1) + try: + dd['b'] + except KeyError as e: + assert_equal(eval(six.text_type(e)), expected_string) + else: + # did not raise a KeyError + assert (False) + + +def test_large_verbosedict(): + expected_sting = ("You tried to access the key 'a' " + "which does not exist. There are 100 " + "extant keys, which is too many to show you") + + dd = core.verbosedict() + for j in range(100): + dd[j] = j + # test success + for j in range(100): + assert_equal(dd[j], j) + # test failure + try: + dd['a'] + except KeyError as e: + assert_equal(eval(six.text_type(e)), expected_sting) + else: + # did not raise a KeyError + assert (False) + + +def test_d_q_conversion(): + assert_equal(2 * np.pi, core.d_to_q(1)) + assert_equal(2 * np.pi, core.q_to_d(1)) + test_data = np.linspace(.1, 5, 100) + assert_array_almost_equal(test_data, core.d_to_q(core.q_to_d(test_data)), + decimal=12) + assert_array_almost_equal(test_data, core.q_to_d(core.d_to_q(test_data)), + decimal=12) + + +def test_q_twotheta_conversion(): + wavelength = 1 + q = np.linspace(0, 4 * np.pi, 100) + assert_array_almost_equal(q, + core.twotheta_to_q( + core.q_to_twotheta(q, wavelength), + wavelength), + decimal=12) + two_theta = np.linspace(0, np.pi, 100) + assert_array_almost_equal(two_theta, + core.q_to_twotheta( + core.twotheta_to_q(two_theta, + wavelength), + wavelength), + decimal=12) + + +def test_radius_to_twotheta(): + dist_sample = 100 + radius = np.linspace(50, 100) + + two_theta = np.array( + [0.46364761, 0.47177751, 0.47984053, 0.48783644, 0.49576508, + 0.5036263, 0.51142, 0.51914611, 0.52680461, 0.53439548, + 0.54191875, 0.54937448, 0.55676277, 0.56408372, 0.57133748, + 0.57852421, 0.58564412, 0.5926974, 0.59968432, 0.60660511, + 0.61346007, 0.62024949, 0.62697369, 0.63363301, 0.6402278, + 0.64675843, 0.65322528, 0.65962874, 0.66596924, 0.67224718, + 0.67846301, 0.68461716, 0.6907101, 0.69674228, 0.70271418, + 0.70862627, 0.71447905, 0.720273, 0.72600863, 0.73168643, + 0.73730693, 0.74287063, 0.74837805, 0.75382971, 0.75922613, + 0.76456784, 0.76985537, 0.77508925, 0.78027, 0.78539816]) + + assert_array_almost_equal(two_theta, + core.radius_to_twotheta(dist_sample, + radius), decimal=8) + + +def test_multi_tau_lags(): + levels = 3 + channels = 8 + + delay_steps = [0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24, 28] + dict_dly = {} + dict_dly[1] = (0, 1, 2, 3, 4, 5, 6, 7) + dict_dly[3] = (16, 20, 24, 28) + dict_dly[2] = (8, 10, 12, 14) + tot_channels, lag_steps, dict_lags = core.multi_tau_lags(levels, channels) + + assert_array_equal(16, tot_channels) + assert_array_equal(delay_steps, lag_steps) + assert_array_equal(dict_dly[1], dict_lags[1]) + assert_array_equal(dict_dly[3], dict_lags[3]) + + +def test_wedge_integration(): + with pytest.raises(NotImplementedError): + core.wedge_integration(src_data=None, center=None, theta_start=None, + delta_theta=None, r_inner=None, delta_r=None) + + +def test_subtract_reference_images(): + num_images = 10 + img_dims = 200 + ones = np.ones((img_dims, img_dims)) + img_lst = [ones * _ for _ in range(num_images)] + img_arr = np.asarray(img_lst) + is_dark_lst = [True] + is_dark = False + was_dark = True + while len(is_dark_lst) < num_images: + if was_dark: + is_dark = False + else: + is_dark = np.random.rand() > 0.5 + was_dark = is_dark + is_dark_lst.append(is_dark) + + is_dark_arr = np.asarray(is_dark_lst) + # make sure that a list of 2d images can be passed in + core.subtract_reference_images(imgs=img_lst, is_reference=is_dark_arr) + # make sure that the reference arr can actually be a list + core.subtract_reference_images(imgs=img_arr, is_reference=is_dark_lst) + # make sure that both input arrays can actually be lists + core.subtract_reference_images(imgs=img_arr, is_reference=is_dark_lst) + + # test that the number of returned images is equal to the expected number + # of returned images + num_expected_images = is_dark_lst.count(False) + # subtract an additional value if the last image is a reference image + # num_expected_images -= is_dark_lst[len(is_dark_lst)-1] + subtracted = core.subtract_reference_images(img_lst, is_dark_lst) + try: + assert_equal(num_expected_images, len(subtracted)) + except AssertionError as ae: + print('is_dark_lst: {0}'.format(is_dark_lst)) + print('num_expected_images: {0}'.format(num_expected_images)) + print('len(subtracted): {0}'.format(len(subtracted))) + six.reraise(AssertionError, ae, sys.exc_info()[2]) + # test that the image subtraction values are behaving as expected + img_sum_lst = [img_dims * img_dims * val for val in range(num_images)] + expected_return_val = 0 + dark_val = 0 + for idx, (is_dark, img_val) in enumerate(zip(is_dark_lst, img_sum_lst)): + if is_dark: + dark_val = img_val + else: + expected_return_val = expected_return_val - dark_val + img_val + # test that the image subtraction was actually processed correctly + return_sum = sum(subtracted) + try: + while True: + return_sum = sum(return_sum) + except TypeError: + # thrown when return_sum is a single number + pass + + try: + assert_equal(expected_return_val, return_sum) + except AssertionError as ae: + print('is_dark_lst: {0}'.format(is_dark_lst)) + print('expected_return_val: {0}'.format(expected_return_val)) + print('return_sum: {0}'.format(return_sum)) + six.reraise(AssertionError, ae, sys.exc_info()[2]) + + +def _fail_img_to_relative_xyi_helper(input_dict): + with pytest.raises(ValueError): + core.img_to_relative_xyi(**input_dict) + + +param_test_img_to_relative_fails = [ + # invalid values of x and y + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_x': -1, + 'pixel_size_y': -1}, + # valid value of x, no value for y + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_x': 1}, + # valid value of y, no value for x + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_y': 1}, + # valid value of y, invalid value for x + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_x': -1, + 'pixel_size_y': 1}, + # valid value of x, invalid value for y + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_x': 1, + 'pixel_size_y': -1}, + # invalid value of x, no value for y + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_x': -1}, + # invalid value of y, no value for x + {'img': np.ones((100, 100)), 'cx': 50, 'cy': 50, 'pixel_size_y': -1}] + + +@pytest.mark.parametrize("fail_dict", param_test_img_to_relative_fails) +def test_img_to_relative_fails(fail_dict): + _fail_img_to_relative_xyi_helper(fail_dict) + + +def test_img_to_relative_xyi(random_seed=None): + from skbeam.core.utils import img_to_relative_xyi + # make the RNG deterministic + if random_seed is not None: + np.random.seed(42) + # set the maximum image dims + maxx = 2000 + maxy = 2000 + # create a randomly sized image + nx = int(np.random.rand() * maxx) + ny = int(np.random.rand() * maxy) + # create a randomly located center + cx = np.random.rand() * nx + cy = np.random.rand() * ny + # generate the image + img = np.ones((nx, ny)) + # generate options for the x center to test edge conditions + cx_lst = [0, cx, nx] + # generate options for the y center to test edge conditions + cy_lst = [0, cy, ny] + for cx, cy in zip(cx_lst, cy_lst): + # call the function + x, y, i = img_to_relative_xyi(img=img, cx=cx, cy=cy) + logger.debug('y {0}'.format(y)) + logger.debug('sum(y) {0}'.format(sum(y))) + expected_total_y = sum(np.arange(ny, dtype=np.int64) - cy) * nx + logger.debug('expected_total_y {0}'.format(expected_total_y)) + logger.debug('x {0}'.format(x)) + logger.debug('sum(x) {0}'.format(sum(x))) + expected_total_x = sum(np.arange(nx, dtype=np.int64) - cx) * ny + logger.debug('expected_total_x {0}'.format(expected_total_x)) + expected_total_intensity = nx * ny + try: + assert_almost_equal(sum(x), expected_total_x, decimal=0) + assert_almost_equal(sum(y), expected_total_y, decimal=0) + assert_equal(sum(i), expected_total_intensity) + except AssertionError as ae: + logger.error('img dims: ({0}, {1})'.format(nx, ny)) + logger.error('img center: ({0}, {1})'.format(cx, cy)) + logger.error('sum(returned_x): {0}'.format(sum(x))) + logger.error('expected_x: {0}'.format(expected_total_x)) + logger.error('sum(returned_y): {0}'.format(sum(y))) + logger.error('expected_y: {0}'.format(expected_total_y)) + logger.error('sum(returned_i): {0}'.format(sum(i))) + logger.error('expected_x: {0}'.format(expected_total_intensity)) + six.reraise(AssertionError, ae, sys.exc_info()[2]) + + +def run_image_to_relative_xyi_repeatedly(): + level = logging.ERROR + ch = logging.StreamHandler() + ch.setLevel(level) + logger.addHandler(ch) + logger.setLevel(level) + + num_calls = 0 + while True: + test_img_to_relative_xyi() + num_calls += 1 + if num_calls % 10 == 0: + print('{0} calls successful'.format(num_calls)) + + +def test_angle_grid(): + a = core.angle_grid((3, 3), (7, 7)) + assert_equal(a[3, -1], 0) + assert_almost_equal(a[3, 0], np.pi) + assert_almost_equal(a[4, 4], np.pi / 4) # (1, 1) should be 45 degrees + # The documented domain is [-pi, pi]. + correct_domain = np.all((a < np.pi + 0.1) & (a > -np.pi - 0.1)) + assert correct_domain + + +def test_radial_grid(): + a = core.radial_grid((3, 3), (7, 7)) + assert_equal(a[3, 3], 0) + assert_equal(a[3, 4], 1) + + +def test_geometric_series(): + time_series = core.geometric_series(common_ratio=5, number_of_images=150) + + assert_array_equal(time_series, [1, 5, 25, 125]) + + +def test_bin_grid(): + if not pf: + pytest.skip("'Geometry' can not be imported from 'pyFAI.geometry'.") + geo = Geometry( + detector='Perkin', pixel1=.0002, pixel2=.0002, + dist=.23, + poni1=.209, poni2=.207, + # poni1=0, poni2=0, + # rot1=.0128, rot2=-.015, rot3=-5.2e-8, + wavelength=1.43e-11 + ) + r_array = geo.rArray((2048, 2048)) + img = r_array.copy() + x, y = core.bin_grid(img, r_array, (geo.pixel1, geo.pixel2)) + + assert_array_almost_equal(y, x, decimal=2) diff --git a/skbeam_copied/skbeam/core/tests/utils.py b/skbeam_copied/skbeam/core/tests/utils.py new file mode 100644 index 00000000..5cd84316 --- /dev/null +++ b/skbeam_copied/skbeam/core/tests/utils.py @@ -0,0 +1,11 @@ +from __future__ import print_function, absolute_import, division + +import numpy as np + + +def parabola_gen(x, center, height, width): + return width * (x-center)**2 + height + + +def gauss_gen(x, center, height, width): + return height * np.exp(-((x-center) / width)**2) diff --git a/skbeam_copied/skbeam/core/utils.py b/skbeam_copied/skbeam/core/utils.py new file mode 100644 index 00000000..ed061954 --- /dev/null +++ b/skbeam_copied/skbeam/core/utils.py @@ -0,0 +1,1276 @@ +#! encoding: utf-8 +# ###################################################################### +# 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. # +######################################################################## +""" +This module is for the 'core' data types. +""" +from __future__ import absolute_import, division, print_function +import six + +from six.moves import zip +from six import string_types + +import time +import sys + +from collections import namedtuple, defaultdict, deque +from collections.abc import MutableMapping +import numpy as np +from itertools import tee + +import logging +import scipy.stats as sts + +logger = logging.getLogger(__name__) + +md_value = namedtuple("md_value", ['value', 'units']) + +_defaults = { + "bins": 100, + 'nx': 100, + 'ny': 100, + 'nz': 100 +} + + +class NotInstalledError(ImportError): + ''' + Custom exception that should be subclassed to handle + specific missing libraries + + ''' + pass + + +class MD_dict(MutableMapping): + """ + A class to make dealing with the meta-data scheme for DataExchange easier + + Examples + -------- + Getting and setting data by path is possible + + >>> tt = MD_dict() + >>> tt['name'] = 'test' + >>> tt['nested.a'] = 2 + >>> tt['nested.b'] = (5, 'm') + >>> tt['nested.a'].value + 2 + >>> tt['nested.a'].units is None + True + >>> tt['name'].value + 'test' + >>> tt['nested.b'].units + 'm' + """ + + def __init__(self, md_dict=None): + # TODO properly walk the input on upgrade dicts -> MD_dict + if md_dict is None: + md_dict = dict() + + self._dict = md_dict + self._split = '.' + + def __repr__(self): + return self._dict.__repr__() + + # overload __setitem__ so dotted paths work + def __setitem__(self, key, val): + + key_split = key.split(self._split) + tmp = self._dict + for k in key_split[:-1]: + try: + tmp = tmp[k]._dict + except Exception: + tmp[k] = type(self)() + tmp = tmp[k]._dict + if isinstance(tmp, md_value): + # TODO make message better + raise KeyError("trying to use a leaf node as a branch") + + # if passed in an md_value, set it and return + if isinstance(val, md_value): + tmp[key_split[-1]] = val + return + # catch the case of a bare string + elif isinstance(val, string_types): + # a value with out units + tmp[key_split[-1]] = md_value(val, 'text') + return + # not something easy, try to guess what to do instead + try: + # if the second element is a string or None, cast to named tuple + if isinstance(val[1], string_types) or val[1] is None: + print('here') + tmp[key_split[-1]] = md_value(*val) + # else, assume whole thing is the value with no units + else: + tmp[key_split[-1]] = md_value(val, None) + # catch any type errors from trying to index into non-indexable things + # or from trying to use iterables longer than 2 + except TypeError: + tmp[key_split[-1]] = md_value(val, None) + + def __getitem__(self, key): + key_split = key.split(self._split) + tmp = self._dict + for k in key_split[:-1]: + try: + tmp = tmp[k]._dict + except Exception: + tmp[k] = type(self)() + tmp = tmp[k]._dict + + if isinstance(tmp, md_value): + # TODO make message better + raise KeyError("trying to use a leaf node as a branch") + + return tmp.get(key_split[-1], None) + + def __delitem__(self, key): + # pass one delete the entry + # TODO make robust to non-keys + key_split = key.split(self._split) + tmp = self._dict + for k in key_split[:-1]: + # make sure we are grabbing the internal dict + tmp = tmp[k]._dict + del tmp[key_split[-1]] + # TODO pass 2 remove empty branches + + def __len__(self): + return len(list(iter(self))) + + def __iter__(self): + return _iter_helper([], self._split, self._dict) + + +def _iter_helper(path_list, split, md_dict): + """ + Recursively walk the tree and return the names of the leaves + """ + for k, v in six.iteritems(md_dict): + if isinstance(v, md_value): + yield split.join(path_list + [k]) + else: + for inner_v in _iter_helper(path_list + [k], split, v._dict): + yield inner_v + + +class verbosedict(dict): + """ + A sub-class of dict which raises more verbose errors if + a key is not found. + """ + + def __getitem__(self, key): + try: + v = dict.__getitem__(self, key) + except KeyError: + if len(self) < 25: + new_msg = ("You tried to access the key '{key}' " + "which does not exist. The " + "extant keys are: {valid_keys}").format( + key=key, valid_keys=list(self)) + else: + new_msg = ("You tried to access the key '{key}' " + "which does not exist. There " + "are {num} extant keys, which is too many to " + "show you").format( + key=key, num=len(self)) + six.reraise(KeyError, KeyError(new_msg), sys.exc_info()[2]) + return v + + +class RCParamDict(MutableMapping): + """A class to make dealing with storing default values easier. + + RC params is a hold- over from the UNIX days where configuration + files are 'rc' files. See + http://en.wikipedia.org/wiki/Configuration_file + + Examples + -------- + Getting and setting data by path is possible + + >>> tt = RCParamDict() + >>> tt['name'] = 'test' + >>> tt['nested.a'] = 2 + """ + _delim = '.' + + def __init__(self): + # the dict to hold the keys at this level + self._dict = dict() + # the defaultdict (defaults to just accepting it) of + # validator functions + self._validators = defaultdict(lambda: lambda x: True) + + # overload __setitem__ so dotted paths work + def __setitem__(self, key, val): + # try to split the key + splt_key = key.split(self._delim, 1) + # if more than one part, recurse + if len(splt_key) > 1: + try: + tmp = self._dict[splt_key[0]] + except KeyError: + tmp = RCParamDict() + self._dict[splt_key[0]] = tmp + + if not isinstance(tmp, RCParamDict): + raise KeyError("name space is borked") + + tmp[splt_key[1]] = val + else: + if not self._validators[key]: + # TODO improve the validation error + raise ValueError("fails to validate, improve this") + self._dict[key] = val + + def __getitem__(self, key): + # try to split the key + splt_key = key.split(self._delim, 1) + if len(splt_key) > 1: + return self._dict[splt_key[0]][splt_key[1]] + else: + return self._dict[key] + + def __delitem__(self, key): + splt_key = key.split(self._delim, 1) + if len(splt_key) > 1: + self._dict[splt_key[0]].__delitem__(splt_key[1]) + else: + del self._dict[key] + + def __len__(self): + return len(list(iter(self))) + + def __iter__(self): + return self._iter_helper([]) + + def _iter_helper(self, path_list): + """ + Recursively walk the tree and return the names of the leaves + """ + for key, val in six.iteritems(self._dict): + if isinstance(val, RCParamDict): + for k in val._iter_helper(path_list + [key, ]): + yield k + else: + yield self._delim.join(path_list + [key, ]) + + def __repr__(self): + # recursively get the formatted list of strings + str_list = self._repr_helper(0) + # return as a single string + return '\n'.join(str_list) + + def _repr_helper(self, tab_level): + # to accumulate the strings into + str_list = [] + # list of the elements at this level + elm_list = [] + # list of sub-levels + nested_list = [] + # loop over the local _dict and sort out which + # keys are nested and which are this level + for key, val in six.iteritems(self._dict): + if isinstance(val, RCParamDict): + nested_list.append(key) + else: + elm_list.append(key) + + # sort the keys in both lists + elm_list.sort() + nested_list.sort() + + # loop over and format the keys/vals at this level + for elm in elm_list: + str_list.append(" " * tab_level + + "{key}: {val}".format( + key=elm, val=self._dict[elm])) + # deal with the nested groups + for nested in nested_list: + # add the label for the group name + str_list.append(" " * tab_level + + "{key}:".format(key=nested)) + # add the strings from _all_ the nested groups + str_list.extend( + self._dict[nested]._repr_helper(tab_level + 1)) + return str_list + + +keys_core = { + "pixel_size": { + "description": ("2 element tuple defining the (x y) dimensions of the " + "pixel"), + "type": tuple, + "units": "um", + }, + "voxel_size": { + "description": ("3 element tuple defining the (x y z) dimensions " + "of the voxel"), + "type": tuple, + "units": "um" + }, + "calibrated_center": { + "description": ("2 element tuple defining the (x y) center of the " + "detector in pixels"), + "type": tuple, + "units": "pixel", + }, + "detector_size": { + "description": ("2 element tuple defining no. of pixels(size) in the " + "detector X and Y direction"), + "type": tuple, + "units": "pixel", + }, + "detector_tilt_angles": { + "description": "Detector tilt angle", + "type": tuple, + "units": " degrees", + }, + "dist_sample": { + "description": "distance from the sample to the detector (mm)", + "type": float, + "units": "mm", + }, + "wavelength": { + "description": "wavelength of incident radiation (Angstroms)", + "type": float, + "units": "angstrom", + }, + "ub_mat": { + "description": "UB matrix(orientation matrix) 3x3 array", + "type": "ndarray", + }, + "energy": { + "description": "scanning energy for data collection", + "type": float, + "units": "keV", + }, + "array_dimensions": { + "description": "axial lengths of the array (Pixels)", + "x_dimension": { + "description": "x-axis array length as int", + "type": int, + "units": "pixels" + }, + "y_dimension": { + "description": "y-axis array length as int", + "type": int, + "units": "pixels" + }, + "z_dimension": { + "description": "z-axis array length as int", + "type": int, + "units": "pixels" + } + }, + "bounding_box": { + "description": ("physical extents of the array: useful for ", + "volume alignment, transformation, merge and ", + "spatial comparison of multiple volumes"), + "x_min": { + "description": "minimum spatial coordinate along the x-axis", + "type": float, + "units": "um" + }, + "x_max": { + "description": "maximum spatial coordinate along the x-axis", + "type": float, + "units": "um" + }, + "y_min": { + "description": "minimum spatial coordinate along the y-axis", + "type": float, + "units": "um" + }, + "y_max": { + "description": "maximum spatial coordinate along the y-axis", + "type": float, + "units": "um" + }, + "z_min": { + "description": "minimum spatial coordinate along the z-axis", + "type": float, + "units": "um" + }, + "z_max": { + "description": "maximum spatial coordinate along the z-axis", + "type": float, + "units": "um" + }, + }, +} + + +def subtract_reference_images(imgs, is_reference): + """ + Function to subtract a series of measured images from + background/dark current/reference images. The nearest reference + image in the reverse temporal direction is subtracted from each + measured image. + + Parameters + ---------- + imgs : numpy.ndarray + Array of 2-D images + + is_reference : 1-D boolean array + true : image is reference image + false : image is measured image + + Returns + ------- + img_corr : numpy.ndarray + len(img_corr) == len(img_arr) - len(is_reference_img == true) + img_corr is the array of measured images minus the reference + images. + + Raises + ------ + ValueError + Possible causes: + is_reference contains no true values + Raised when the first image in the array is not a reference image. + + """ + # an array of 1, 0, 1,.. should work too + if not is_reference[0]: + # use ValueError because the user passed in invalid data + raise ValueError("The first image is not a reference image") + # grab the first image + ref_imge = imgs[0] + # make an array of zeros of the correct type + corrected_image = deque() + # zip together (lazy like this is really izip), images and flags + for imgs, ref in zip(imgs[1:], is_reference[1:]): + # if this is a ref image, save it and move on + if ref: + ref_imge = imgs + continue + # else, do the subtraction + corrected_image.append(imgs - ref_imge) + + # return the output as a list + return list(corrected_image) + + +def img_to_relative_xyi(img, cx, cy, pixel_size_x=None, pixel_size_y=None): + """ + Convert the 2D image to a list of x y I coordinates where + x == x_img - detector_center[0] and + y == y_img - detector_center[1] + + Parameters + ---------- + img: `ndarray` + 2D image + cx : float + Image center in the x direction + cy : float + Image center in the y direction + pixel_size_x : float, optional + Pixel size in x + pixel_size_y : float, optional + Pixel size in y + **kwargs: dict + Bucket for extra parameters in an unpacked dictionary + + Returns + ------- + x : `ndarray` + x-coordinate of pixel. shape (N, ) + y : `ndarray` + y-coordinate of pixel. shape (N, ) + I : `ndarray` + intensity of pixel. shape (N, ) + """ + if pixel_size_x is not None and pixel_size_y is not None: + if pixel_size_x <= 0: + raise ValueError('Input parameter pixel_size_x must be greater ' + 'than 0. Your value was ' + + six.text_type(pixel_size_x)) + if pixel_size_y <= 0: + raise ValueError('Input parameter pixel_size_y must be greater ' + 'than 0. Your value was ' + + six.text_type(pixel_size_y)) + elif pixel_size_x is None and pixel_size_y is None: + pixel_size_x = 1 + pixel_size_y = 1 + else: + raise ValueError('pixel_size_x and pixel_size_y must both be None or ' + 'greater than zero. You passed in values for ' + 'pixel_size_x of {0} and pixel_size_y of {1}' + ''.format(pixel_size_x, pixel_size_y)) + + # Caswell's incredible terse rewrite + x, y = np.meshgrid(pixel_size_x * (np.arange(img.shape[0]) - cx), + pixel_size_y * (np.arange(img.shape[1]) - cy)) + + # return x, y and intensity as 1D arrays + return x.ravel(), y.ravel(), img.ravel() + + +def bin_1D(x, y, nx=None, min_x=None, max_x=None): + """ + Bin the values in y based on their x-coordinates + + Parameters + ---------- + x : array + position + y : array + intensity + nx : integer, optional + number of bins to use defaults to default bin value + min_x : float, optional + Left edge of first bin defaults to minimum value of x + max_x : float, optional + Right edge of last bin defaults to maximum value of x + + Returns + ------- + edges : array + edges of bins, length nx + 1 + + val : array + sum of values in each bin, length nx + + count : array + The number of counts in each bin, length nx + """ + + # handle default values + if min_x is None: + min_x = np.min(x) + if max_x is None: + max_x = np.max(x) + if nx is None: + nx = int(max_x - min_x) + + # use a weighted histogram to get the bin sum + bins = np.linspace(start=min_x, stop=max_x, num=nx + 1, endpoint=True) + val, _ = np.histogram(a=x, bins=bins, weights=y) + # use an un-weighted histogram to get the counts + count, _ = np.histogram(a=x, bins=bins) + # return the three arrays + return bins, val, count + + +def radial_grid(center, shape, pixel_size=None): + """Convert a cartesian grid (x,y) to the radius relative to some center + + Parameters + ---------- + center : tuple + point in image where r=0; may be a float giving subpixel precision. + Order is (rr, cc). + shape : tuple + Image shape which is used to determine the maximum extent of output + pixel coordinates. + Order is (rr, cc). + pixel_size : sequence, optional + The physical size of the pixels. + len(pixel_size) should be the same as len(shape) + defaults to (1,1) + + Returns + ------- + r : array + The distance of each pixel from `center` + Shape of the return value is equal to the `shape` input parameter + """ + + if pixel_size is None: + pixel_size = (1, 1) + + X, Y = np.meshgrid(pixel_size[1] * (np.arange(shape[1]) - center[1]), + pixel_size[0] * (np.arange(shape[0]) - center[0])) + return np.sqrt(X * X + Y * Y) + + +def angle_grid(center, shape, pixel_size=None): + """ + Make a grid of angular positions. + + Read note for our conventions here -- there be dragons! + + Parameters + ---------- + center : tuple + point in image where r=0; may be a float giving subpixel precision. + Order is (rr, cc). + + shape: tuple + Image shape which is used to determine the maximum extent of output + pixel coordinates. Order is (rr, cc). + + Returns + ------- + agrid : array + angular position (in radians) of each array element in range [-pi, pi] + + Note + ---- + :math:`\\theta`, the counter-clockwise angle from the positive x axis, + assuming the positive y-axis points upward. + :math:`\\theta \\el [-\\pi, \\pi]`. In array indexing and the conventional + axes for images (origin in upper left), positive y is downward. + """ + + if pixel_size is None: + pixel_size = (1, 1) + + # row is y, column is x. "so say we all. amen." + x, y = np.meshgrid(pixel_size[1] * (np.arange(shape[1]) - + center[1]), + pixel_size[0] * (np.arange(shape[0]) - + center[0])) + return np.arctan2(y, x) + + +def radius_to_twotheta(dist_sample, radius): + """ + Converts radius from the calibrated center to scattering angle + (2:math:`2\\theta`) with known detector to sample distance. + + Parameters + ---------- + dist_sample : float + distance from the sample to the detector (mm) + + radius : array + The L2 norm of the distance of each pixel from the calibrated center. + + Returns + ------- + two_theta : array + An array of :math:`2\\theta` values + """ + return np.arctan(radius / dist_sample) + + +def wedge_integration(src_data, center, theta_start, + delta_theta, r_inner, delta_r): + """ + Implementation of caking. + + Parameters + ---------- + scr_data : ndarray + The source-data to be integrated + + center : ndarray + The center of the ring in pixels + + theta_start : float + The angle of the start of the wedge from the + image y-axis in degrees + + delta_theta : float + The angular width of the wedge in degrees. Positive + angles go clockwise, negative go counter-clockwise. + + r_inner : float + The inner radius in pixel units, Must be non-negative + + delta_r : float + The length of the wedge in the radial direction + in pixel units. Must be non-negative + + Returns + ------- + float + The integrated intensity under the wedge + """ + raise NotImplementedError() + + +def bin_edges(range_min=None, range_max=None, nbins=None, step=None): + """ + Generate bin edges. The last value in the returned array is + the right edge of the last bin, the rest of the values are the + left edges of each bin. + + If `range_max` is specified all bin edges will be less than or + equal to it's value. + + If `range_min` is specified all bin edges will be greater than + or equal to it's value + + If `nbins` is specified then there will be than number of bins and + the returned array will have length `nbins + 1` (as the right most + edge is included) + + If `step` is specified then bin width is approximately `step` (It is + not exact due to the nature of floats). The arrays generated by + `np.cumsum(np.ones(nbins) * step)` and `np.arange(nbins) * step` are + not identical. This function uses the second method in all cases + where `step` is specified. + + .. warning :: If the set :code:`(range_min, range_max, step)` is + given there is no guarantee that :code:`range_max - range_min` + is an integer multiple of :code:`step`. In this case the left + most bin edge is :code:`range_min` and the right most bin edge + is less than :code:`range_max` and the distance between the + right most edge and :code:`range_max` is not greater than + :code:`step` (this is the same behavior as the built-in + :code:`range()`). It is not recommended to specify bins in this + manner. + + Parameters + ---------- + range_min : float, optional + The minimum value that may be included as a bin edge + + range_max : float, optional + The maximum value that may be included as a bin edge + + nbins : int, optional + The number of bins, if specified the length of the returned + value will be nbins + 1 + + step : float, optional + The step between the bins + + Returns + ------- + edges : np.array + An array of floats for the bin edges. The last value is the + right edge of the last bin. + """ + num_valid_args = sum((range_min is not None, range_max is not None, + step is not None, nbins is not None)) + if num_valid_args != 3: + raise ValueError("Exactly three of the arguments must be non-None " + "not {}.".format(num_valid_args)) + + if range_min is not None and range_max is not None: + if range_max <= range_min: + raise ValueError("The minimum must be less than the maximum") + + if nbins is not None: + if nbins <= 0: + raise ValueError("The number of bins must be positive") + + # The easy case + if step is None: + return np.linspace(range_min, range_max, nbins + 1, endpoint=True) + + # in this case, the user gave use min, max, and step + if nbins is None: + if step > (range_max - range_min): + raise ValueError("The step can not be greater than the difference " + "between min and max") + nbins = int((range_max - range_min) // step) + ret = range_min + np.arange(nbins + 1) * step + # if the last value is greater than the max (should never happen) + if ret[-1] > range_max: + return ret[:-1] + if range_max - ret[-1] > 1e-10 * step: + logger.debug("Inconsistent " + "(range_min, range_max, step) " + "and step does not evenly divide " + "(range_min - range_max). " + "The bins has been truncated.\n" + "min: %f max: %f step: %f gap: %f", + range_min, range_max, + step, range_max - ret[-1]) + return ret + + # in this case we got range_min, nbins, step + if range_max is None: + return range_min + np.arange(nbins + 1) * step + + # in this case we got range_max, nbins, step + if range_min is None: + return range_max - np.arange(nbins + 1)[::-1] * step + + +def grid3d(q, img_stack, + nx=None, ny=None, nz=None, + xmin=None, xmax=None, ymin=None, + ymax=None, zmin=None, zmax=None, + binary_mask=None): + """Grid irregularly spaced data points onto a regular grid via histogramming + + This function will process the set of reciprocal space values (q), the + image stack (img_stack) and grid the image data based on the bounds + provided, using defaults if none are provided. + + Parameters + ---------- + q : ndarray + (Qx, Qy, Qz) - HKL values - Nx3 array + img_stack : ndarray + Intensity array of the images + dimensions are: [num_img][num_rows][num_cols] + nx : int, optional + Number of voxels along x + ny : int, optional + Number of voxels along y + nz : int, optional + Number of voxels along z + xmin : float, optional + Minimum value along x. Defaults to smallest x value in q + ymin : float, optional + Minimum value along y. Defaults to smallest y value in q + zmin : float, optional + Minimum value along z. Defaults to smallest z value in q + xmax : float, optional + Maximum value along x. Defaults to largest x value in q + ymax : float, optional + Maximum value along y. Defaults to largest y value in q + zmax : float, optional + Maximum value along z. Defaults to largest z value in q + binary_mask : ndarray, optional + The binary mask provides a mechanism to remove unwanted pixels + from the images. + Binary mask can be two different shapes. + - 1: 2-D with binary_mask.shape == np.asarray(img_stack[0]).shape + - 2: 3-D with binary_mask.shape == np.asarray(img_stack).shape + + Returns + ------- + mean : ndarray + intensity grid. The values in this grid are the + mean of the values that fill with in the grid. + occupancy : ndarray + The number of data points that fell in the grid. + std_err : ndarray + This is the standard error of the value in the + grid box. + bounds : list + tuple of (min, max, step) for x, y, z in order: [x_bounds, + y_bounds, z_bounds] + + """ + try: + from ..ext import ctrans + except ImportError: + raise NotImplementedError( + "ctrans is not available on your platform. See " + "https://github.com/scikit-beam/scikit-beam/issues/418 " + "to follow updates to this problem.") + + # validate input + img_stack = np.asarray(img_stack) + # todo determine if we're going to support masked arrays + # todo masked arrays seemed to have been punted to `process_to_q` + + # check to see if the binary mask and the image stack are identical shapes + if binary_mask is None or binary_mask.shape == img_stack.shape: + # do a dance :) + pass + elif binary_mask.shape == img_stack[0].shape: + # this is still a valid mask, so make it the same dimensions + # as img_stack. + # should probably change this to use something similar to: + # todo http://stackoverflow.com/questions/5564098/ + binary_mask = np.tile(np.ravel(binary_mask), img_stack.shape[0]) + + else: + raise ValueError("The binary mask must be the same shape as the" + "img_stack ({0}) or a single image in the image " + "stack ({1}). The input binary mask is shaped ({2})" + "".format(img_stack.shape, img_stack[0].shape, + binary_mask.shape)) + + q = np.atleast_2d(q) + if q.ndim != 2: + raise ValueError("q.ndim must be a 2-D array of shape Nx3 array. " + "You provided an array with {0} dimensions." + "".format(q.ndim)) + if q.shape[1] != 3: + raise ValueError("The shape of q must be an Nx3 array, not {0}X{1}" + " which you provided.".format(*q.shape)) + + # set defaults for qmin, qmax, dq + qmin = np.min(q, axis=0) + qmax = np.max(q, axis=0) + dqn = [_defaults['nx'], _defaults['ny'], _defaults['nz']] + + # pad the upper edge by just enough to ensure that all of the + # points are in-bounds with the binning rules: lo <= val < hi + qmax += np.spacing(qmax) + + # check for non-default input + for target, input_vals in ((dqn, (nx, ny, nz)), + (qmin, (xmin, ymin, zmin)), + (qmax, (xmax, ymax, zmax))): + for j, in_val in enumerate(input_vals): + if in_val is not None: + target[j] = in_val + + # format bounds + bounds = np.array([qmin, qmax, dqn]).T + + # creating (Qx, Qy, Qz, I) Nx4 array - HKL values and Intensity + # getting the intensity value for each pixel + q = np.insert(q, 3, np.ravel(img_stack), axis=1) + if binary_mask is not None: + q = q[np.ravel(binary_mask)] + + # 3D grid of the data set + # starting time for gridding + t1 = time.time() + + # call the c library + + total, total2, occupancy, std_err = ctrans.grid3d(q, qmin, qmax, dqn) + mean = total / occupancy + + # ending time for the gridding + t2 = time.time() + logger.info("Done processed in {0} seconds".format(t2 - t1)) + + # No. of values zero in the grid + empt_nb = (occupancy == 0).sum() + + # log some information about the grid at the debug level + logger.debug("There are %2e bins in the grid", mean.size) + if empt_nb: + logger.debug("There are %.2e values zero in the grid", empt_nb) + + return mean, occupancy, std_err, bounds + + +def bin_edges_to_centers(input_edges): + """ + Helper function for turning a array of bin edges into + an array of bin centers + + Parameters + ---------- + input_edges : array-like + N + 1 values which are the left edges of N bins + and the right edge of the last bin + + Returns + ------- + centers : ndarray + A length N array giving the centers of the bins + """ + input_edges = np.asarray(input_edges) + return (input_edges[:-1] + input_edges[1:]) * 0.5 + + +# https://docs.python.org/2/library/itertools.html#recipes +def pairwise(iterable): + "s -> (s0,s1), (s1,s2), (s2, s3), ..." + a, b = tee(iterable) + next(b, None) + return zip(a, b) + + +def q_to_d(q): + """ + Helper function to convert :math:`d` to :math:`q`. The point + of this function is to prevent fat-fingered typos. + + By definition the relationship is: + + .. math:: + + q = \\frac{2 \\pi}{d} + + + Parameters + ---------- + q : array + An array of q values + + Returns + ------- + d : array + An array of d (plane) spacing in the inverse of the units of ``q`` + + """ + return (2 * np.pi) / np.asarray(q) + + +def d_to_q(d): + """ + Helper function to convert :math:`d` to :math:`q`. + The point of this function is to prevent fat-fingered typos. + + By definition the relationship is: + + .. math:: + + d = \\frac{2 \\pi}{q} + + Parameters + ---------- + d : array + An array of d (plane) spacing + + Returns + ------- + q : array + An array of q values in the inverse of the units of ``d`` + + + """ + return (2 * np.pi) / np.asarray(d) + + +def q_to_twotheta(q, wavelength): + r""" + Helper function to convert q to two-theta. + + By definition the relationship is: + + .. math:: + + \sin\left(\frac{2\theta}{2}\right) = \frac{\lambda q}{4 \pi} + + thus + + .. math:: + + 2\theta_n = 2 \arcsin\left(\frac{\lambda q}{4 \pi}\right) + + Parameters + ---------- + q : array + An array of :math:`q` values + + wavelength : float + Wavelength of the incoming x-rays + + Returns + ------- + two_theta : array + An array of :math:`2\theta` values in radians + """ + q = np.asarray(q) + wavelength = float(wavelength) + pre_factor = wavelength / (4 * np.pi) + return 2 * np.arcsin(q * pre_factor) + + +def twotheta_to_q(two_theta, wavelength): + r""" + Helper function to convert two-theta to q + + By definition the relationship is + + .. math:: + + \sin\left(\frac{2\theta}{2}\right) = \frac{\lambda q}{4 \pi} + + thus + + .. math:: + + q = \frac{4 \pi \sin\left(\frac{2\theta}{2}\right)}{\lambda} + + + + Parameters + ---------- + two_theta : array + An array of :math:`2\theta` values + + wavelength : float + Wavelength of the incoming x-rays + + Returns + ------- + q : array + An array of :math:`q` values in the inverse of the units + of ``wavelength`` + """ + two_theta = np.asarray(two_theta) + wavelength = float(wavelength) + pre_factor = ((4 * np.pi) / wavelength) + return pre_factor * np.sin(two_theta / 2) + + +def multi_tau_lags(multitau_levels, multitau_channels): + """ + Standard multiple-tau algorithm for finding the lag times (delay + times). + + Parameters + ---------- + multitau_levels : int + number of levels of multiple-taus + multitau_channels : int + number of channels or number of buffers in auto-correlators + normalizations (must be even) + + Returns + ------- + total_channels : int + total number of channels ( or total number of delay times) + lag_steps : ndarray + delay or lag steps for the multiple tau analysis + dict_lags : dict + dictionary of delays for each multitau_levels + + Notes + ----- + The multi-tau correlation scheme was used for finding the lag times + (delay times). + + References: text [1]_ + + .. [1] K. Schätzela, M. Drewela and S. Stimaca, "Photon correlation + measurements at large lag times: Improving statistical accuracy," + J. Mod. Opt., vol 35, p 711–718, 1988. + """ + + if (multitau_channels % 2 != 0): + raise ValueError("Number of multiple tau channels(buffers)" + " must be even. You provided {0} " + .format(multitau_channels)) + + # total number of channels ( or total number of delay times) + tot_channels = (multitau_levels + 1) * multitau_channels // 2 + + lag = [] + dict_lags = {} + lag_steps = np.arange(0, multitau_channels) + dict_lags[1] = lag_steps + for i in range(2, multitau_levels + 1): + y = [] + for j in range(0, multitau_channels // 2): + value = (multitau_channels // 2 + j) * (2 ** (i - 1)) + lag.append(value) + y.append(value) + dict_lags[i] = y + + lag_steps = np.append(lag_steps, np.array(lag)) + return tot_channels, lag_steps, dict_lags + + +def geometric_series(common_ratio, number_of_images, first_term=1): + """ + This will provide the geometric series for the integration. + Last values of the series has to be less than or equal to number + of images + ex: number_of_images = 100, first_term =1 + common_ratio = 2, geometric_series = 1, 2, 4, 8, 16, 32, 64 + common_ratio = 3, geometric_series = 1, 3, 9, 27, 81 + + Parameters + ---------- + common_ratio : float + common ratio of the series + + number_of_images : int + number of images + + first_term : float, optional + first term in the series + + Return + ------ + geometric_series : list + time series + + Notes + ----- + .. math:: + a + ar + ar^2 + ar^3 + ar^4 + ... + + a - first term in the series + + r - is the common ratio + """ + + geometric_series = [first_term] + + while geometric_series[-1] * common_ratio < number_of_images: + geometric_series.append(geometric_series[-1] * common_ratio) + return geometric_series + + +def bin_grid(image, r_array, pixel_sizes, statistic='mean', mask=None, + bins=None): + """ + Bin and integrate an image, given the radial array of pixels + + Parameters + ---------- + image: np.array + The image in quesion + r_array: np.array + The array which maps pixel positions to tilt/rotation corrected radii + pixel_sizes: tuple + The size of the pixels in the same units as the r_array + statistic: str or func, optional + The statistic to compute over the integration, defaults to mean + mask: bool array, optional + The array of pixels to be removed from the image before integration + bins: array, optional + The bins to use in the integration, if none given the function will + give its best assessment based on the pixel_size and r_array + + Returns + ------- + bin_centers : array + The center of each bin in R + int_stat : array + Radial integrated statistic of the image. + + See Also + -------- + circular_average : circularly average an image, assuming linear radial + spacing (less general) + + """ + if mask is None: + mask = np.ones(image.shape, dtype=int).astype(bool) + if bins is None: + res = np.hypot(*pixel_sizes) + bins = np.arange(np.min(r_array) - res * .5, + np.max(r_array) + res * .5, res) + + int_stat, bin_edge, bin_num = sts.binned_statistic(r_array[mask], + image[mask], + statistic=statistic, + bins=bins) + + bin_centers = bin_edges_to_centers(bin_edge) + + return bin_centers, int_stat