Merge pull request #16 from LaurentRDC/develop

Develop

docs/source/tutorials/image.rst

@@ -6,8 +6,8 @@
 Image Analysis Tutorial
 ***********************
 
-Due to the high electron cross-section, background signals (or baseline) are
-much more of a problem for electron diffraction than equivalent X-ray experiments.
+Diffraction patterns analysis is essentially specialized image processing. This tutorial
+will show some of the image processing and analysis techniques that are part of the :mod:`skued.image` module.
 
 Contents
 ========
@@ -78,6 +78,10 @@ memory usage; this allows the use of multiple processes in parallel::
 	    # write to disk of display
 	    pass
 
+Scikit-ued provides a few streaming statistical functions (:func:`ivar`, :func:`istd`, 
+:func:`isem`, :func:`iaverage`) which are tested to have exact same results are the familiar
+:mod:`numpy` and :mod:`scipy` equivalents. You can also find more specialized streaming functions (e.g. :func:`ialign`).
+
 Example: averaging with error
 ------------------------------
 

docs/source/tutorials/index.rst

@@ -23,4 +23,5 @@ We currently have the following tutorials:
    structure
    simulation
    baseline
-   image
+   image
+   plotting

docs/source/tutorials/plotting.rst

@@ -0,0 +1,61 @@
+.. include:: ../references.txt
+
+.. _plotting_tutorial:
+
+*****************
+Plotting Tutorial
+*****************
+
+Time-resolved diffraction measurements can benefit from careful plotting, especially
+for time-series of data. This tutorial will go over some of the functions that make this easier
+in :mod:`skued`.
+
+Contents
+========
+
+* :ref:`colors`
+
+.. _colors:
+
+Colors
+======
+
+Time-order
+----------
+
+Time-resolved data possesses obvious ordering, and as such we can use colors to
+enhance visualization of such data.
+
+My favourite way to plot time-series data is to make use of rainbow colors to indicate time:
+
+.. plot::
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    from skued import spectrum_colors
+
+    t = np.linspace(0, 10, num = 1000)
+    y = np.exp(-t/0.4) + 10 *(1 - np.exp(-t/2)) + np.sin(2*np.pi*t)
+
+    fig, ax = plt.subplots(1,1)
+    ax.scatter(t, y, c = list(spectrum_colors(y.size)))
+
+    ax.set_title('Some realistic dynamics')
+    ax.set_xlabel('Time (ps)')
+    ax.set_ylabel('Intensity (a.u.)')
+
+    plt.show()
+
+This functionality can be accessed through the :func:`spectrum_colors` generator::
+
+    from skued import spectrum_colors, last
+
+    colors = spectrum_colors(5)
+    
+    # Colors always go from purple to red (increasing wavelength)
+    purple = next(colors)
+    red = last(colors)
+
+You can see some examples of uses of :func:`spectrum_colors` in the :ref:`baseline tutorial <baseline_tutorial>`.
+
+:ref:`Return to Top <plotting_tutorial>`

skued/image/__init__.py

@@ -5,4 +5,4 @@
 from .alignment import align, shift_image, diff_register
 from .symmetry import nfold_symmetry, nfold
 from .correlation import mnxc2
-from .streaming import ialign, iaverage, isem
+from .streaming import ialign, iaverage, isem, istd, ivar

skued/image/streaming.py

@@ -2,58 +2,64 @@
 Streaming operations on arrays/images
 =====================================
 """
-from collections import deque
+
 from functools import partial
 from itertools import repeat
+from math import sqrt
 
 import numpy as np
 
 from . import align
 
-def ialign(images, reference = None, fill_value = 0.0):
+def ialign(images, reference = None, mask = None, fill_value = 0.0):
 	"""
-	Generator of aligned diffraction images.
-
-	Parameters
-	----------
-	images : iterable
-		Iterable of ndarrays of shape (N,M)
-	reference : `~numpy.ndarray` or None, optional
-		If not None, this is the reference image to which all images will be aligned. Otherwise,
-		images will be aligned to the first element of the iterable 'images'. 
-	fill_value : float, optional
-		Edges will be filled with `fill_value` after shifting.
-    
-	Yields
-	------
-	aligned : ndarray, ndim 2
-		Aligned image
-
-	Notes
-	-----
-	Diffraction images exhibit high symmetry in most cases, therefore images
-	are cropped to a quarter of their size before alignment.
+    Generator of aligned diffraction images.
+
+    Parameters
+    ----------
+    images : iterable
+        Iterable of ndarrays of shape (N,M)
+    reference : `~numpy.ndarray`, shape (M,N)
+        Images in `images` will be aligned onto the `reference` image. If
+        'reference' is None (default), the first image in the 'images' stream
+        is used as a reference
+    mask : `~numpy.ndarray` or None, optional
+        Mask that evaluates to True on invalid pixels.
+    fill_value : float, optional
+        Edges will be filled with `fill_value` after alignment.
+
+    Yields
+    ------
+    aligned : `~numpy.ndarray`
+        Aligned image. If `reference` is None, the first aligned image is the reference.
+
+    See Also
+    --------
+    skued.image.align : align a single diffraction pattern onto a reference.
 	"""
 	images = iter(images)
 
 	if reference is None:
 		reference = next(images)
 		yield reference
 
-	yield from map(partial(align, reference = reference), images)
+	yield from map(partial(align, reference = reference, mask = mask, fill_value =  fill_value), images)
 
 def iaverage(images, weights = None):
     """ 
-    Streaming average of diffraction images. This generator can be used to 
-    observe a live averaging.
+    Streaming average of diffraction images. This is equivalent
+    to `numpy.average(axis = 2)` for a stack of images.
 
     Parameters
     ----------
     images : iterable of ndarrays
         Images to be averaged. This iterable can also a generator.
     weights : iterable of ndarray, iterable of floats, or None, optional
-        Array of weights. See `numpy.average` for further information. If None (default), 
-        total picture intensity of valid pixels is used to weight each picture.
+        Iterable of weights associated with the values in each item of `images`. 
+        Each value in an element of `images` contributes to the average 
+        according to its associated weight. The weights array can either be a float
+        or an array of the same shape as any element of `images`. If weights=None, 
+        then all data in each element of `images` are assumed to have a weight equal to one.
     
     Yields
     ------
@@ -62,7 +68,7 @@ def iaverage(images, weights = None):
     
     See Also
     --------
-    numpy.average : average for dense arrays
+    numpy.average : (weighted) average for dense arrays
     """
     images = iter(images)
 
@@ -81,43 +87,124 @@ def iaverage(images, weights = None):
         #print(sum_of_weights)
         yield weighted_sum/sum_of_weights
 
-def isem(images):
+def ivar(images, ddof = 1, weights = None):
     """ 
-    Streaming standard error in the mean (SEM) of images. This is equivalent to
-    calling `scipy.mstats.sem` with `ddof = 1`.
-
+    Streaming variance of a set of images, per pixel. Weights are also supported.
+    This is equivalent to calling `numpy.var(axis = 2)` on a stack of images.
+    
     Parameters
     ----------
     images : iterable of ndarrays
-        Images to be averaged. This iterable can also a generator.
+        Images to be combined. This iterable can also a generator.
+    ddof : int, optional
+        Means Delta Degrees of Freedom.  The divisor used in calculations
+        is ``N - ddof``, where ``N`` represents the number of elements.
+        By default `ddof` is one.
+    weights : iterable of ndarray, iterable of floats, or None, optional
+        Iterable of weights associated with the values in each item of `images`. 
+        Each value in an element of `images` contributes to the variance 
+        according to its associated weight. The weights array can either be a float
+        or an array of the same shape as any element of `images`. If weights=None, 
+        then all data in each element of `images` are assumed to have a weight equal to one.
     
     Yields
     ------
-    sem: `~numpy.ndarray`
-        Standard error in the mean. 
+    var: `~numpy.ndarray`
+        Variance on a per-pixel basis. 
     
-    See also
+    See Also
     --------
-    scipy.stats.sem : standard error in the mean of dense arrays.
+    numpy.var : variance calculation for dense arrays. Weights are not supported.
     
     References
     ----------
-    .. [#] D. Knuth, The Art of Computer Programming 3rd Edition, Vol. 2, p. 232
+    .. [#] D. H. D. West, Updating the mean and variance estimates: an improved method.
+        Communications of the ACM Vol. 22, Issue 9, pp. 532 - 535 (1979)
     """
     images = iter(images)
 
+    if weights is None:
+        weights = repeat(1.0)
+    weights = iter(weights)
+    sum_of_weights = np.array(next(weights), copy = True)
+
     first = next(images)
-    old_M = new_M = np.array(first, copy = True)
+    old_mean = new_mean = np.array(first, copy = True)
     old_S = new_S = np.zeros_like(first, dtype = np.float)
     yield np.zeros_like(first)  # No error if no averaging
 
-    for k, image in enumerate(images, start = 2):
+    for image, weight in zip(images, weights):
+
+        sum_of_weights += weight
+
+        _sub = weight * (image - old_mean)
+        new_mean[:] = old_mean + _sub/sum_of_weights
+        new_S[:] = old_S + _sub*(image - new_mean)
+
+        # In case there hasn't been enough measurements yet,
+        # yield zeros.
+        if np.any(sum_of_weights - ddof <= 0):
+            yield np.zeros_like(first)
+        else:
+            yield new_S/(sum_of_weights - ddof) # variance = S / k-1, sem = std / sqrt(k)    
+
+        old_mean[:] = new_mean
+        old_S[:] = new_S
+
+def istd(images, ddof = 1, weights = None):
+    """ 
+    Streaming standard deviation of images. Weights are also supported.
+    This is equivalent to calling `numpy.std(axis = 2)` on a stack of images.
+
+    Parameters
+    ----------
+    images : iterable of ndarrays
+        Images to be combined. This iterable can also a generator.
+    ddof : int, optional
+        Means Delta Degrees of Freedom.  The divisor used in calculations
+        is ``N - ddof``, where ``N`` represents the number of elements.
+        By default `ddof` is one.
+    weights : iterable of ndarray, iterable of floats, or None, optional
+        Iterable of weights associated with the values in each item of `images`. 
+        Each value in an element of `images` contributes to the standard deviation
+        according to its associated weight. The weights array can either be a float
+        or an array of the same shape as any element of `images`. If weights=None, 
+        then all data in each element of `images` are assumed to have a weight equal to one.
+    
+    Yields
+    ------
+    std: `~numpy.ndarray`
+        Standard deviation on a per-pixel basis.
+
+    See Also
+    --------
+    numpy.std : standard deviation calculation of dense arrays. Weights are not supported.
+    """
+    yield from map(np.sqrt, ivar(images, ddof = ddof, weights = weights))
 
-        _sub = image - old_M
-        new_M[:] = old_M + _sub/k
-        new_S[:] = old_S + _sub*(image - new_M)
-
-        yield np.sqrt(new_S/(k*(k-1))) # variance = S / k-1, sem = std / sqrt(k)    
+def isem(images, ddof = 1):
+    """ 
+    Streaming standard error in the mean (SEM) of images. This is equivalent to
+    calling `scipy.stats.sem(axis = 2)` on a stack of images.
 
-        old_M[:] = new_M
-        old_S[:] = new_S
+    Parameters
+    ----------
+    images : iterable of ndarrays
+        Images to be combined. This iterable can also a generator.
+    ddof : int, optional
+        Means Delta Degrees of Freedom.  The divisor used in calculations
+        is ``N - ddof``, where ``N`` represents the number of elements.
+        By default `ddof` is one.
+    
+    Yields
+    ------
+    sem: `~numpy.ndarray`
+        Standard error in the mean. 
+    
+    See Also
+    --------
+    scipy.stats.sem : standard error in the mean of dense arrays.
+    """
+    # TODO: include weights
+    for k, std in enumerate(istd(images, ddof = ddof), start = 1):
+        yield std / sqrt(k)