Merge pull request #15 from LaurentRDC/develop

Develop

docs/source/api.rst

@@ -20,6 +20,10 @@ Array Utilities
 ===============
 .. automodule:: skued.array_utils
 
+Iteration/Generator Utilities
+=============================
+.. automodule:: skued.iter_utils
+
 Quantities
 ==========
 .. automodule:: skued.quantities

docs/source/tutorials/image.rst

@@ -57,6 +57,8 @@ but will not compute anything until it is requested. We can use the function
 :code:`last` to get at the final average, but we could also look at the average
 step-by-step by calling :code:`next`::
 
+	from skued import last
+
 	avg = next(averaged) # only one images is loaded, aligned and added to the average
 	total = last(averaged) # average of the entire stream
 

skued/__init__.py

@@ -2,13 +2,14 @@
 __author__ = 'Laurent P. René de Cotret'
 __email__ = 'laurent.renedecotret@mail.mcgill.ca'
 __license__ = 'MIT'
-__version__ = '0.4.6' # TODO: automatic versioning?
+__version__ = '0.4.7' # TODO: automatic versioning?
 
 from .affine import (affine_map, change_basis_mesh, change_of_basis, is_basis,
                      is_rotation_matrix, minimum_image_distance,
                      rotation_matrix, transform, translation_matrix,
                      translation_rotation_matrix)
 from .array_utils import mirror, repeated_array
+from .iter_utils import chunked, last, linspace, multilinspace
 from .parallel import pmap, preduce
 from .plot_utils import spectrum_colors, rgb_sweep
 from .quantities import electron_wavelength, interaction_parameter, lorentz

skued/image/__init__.py

@@ -3,6 +3,6 @@
 
 from .powder import angular_average, powder_center
 from .alignment import align, shift_image, diff_register
-from .symmetry import nfold_symmetry
+from .symmetry import nfold_symmetry, nfold
 from .correlation import mnxc2
-from .streaming import ialign, iaverage, isem, last
+from .streaming import ialign, iaverage, isem

skued/image/streaming.py

@@ -10,13 +10,6 @@
 
 from . import align
 
-# TODO: move into base package, e.g. iter_utils.py?
-def last(stream):
-    """ Returns the last item from a stream. """
-    # Wonderful idea from itertools recipes
-    # https://docs.python.org/3.6/library/itertools.html#itertools-recipes
-    return deque(stream, maxlen = 1)[0]
-
 def ialign(images, reference = None, fill_value = 0.0):
 	"""
 	Generator of aligned diffraction images.

skued/image/symmetry.py

@@ -3,12 +3,13 @@
 Image manipulation involving symmetry
 =====================================
 """
-from functools import partial
+from functools import partial, wraps
 from skimage.transform import rotate
 import numpy as np
+from warnings import warn
 
 # TODO: out parameter?
-def nfold_symmetry(im, center, mod, **kwargs):
+def nfold(im, mod, center = None, mask = None, **kwargs):
     """ 
     Returns an images averaged according to n-fold rotational symmetry.
     Keyword arguments are passed to skimage.transform.rotate()
@@ -17,25 +18,48 @@ def nfold_symmetry(im, center, mod, **kwargs):
     ----------
     im : array_like, ndim 2
         Image to be averaged.
-	center : array_like, shape (2,)
-		coordinates of the center (in pixels).
+    center : array_like, shape (2,) or None, optional
+        coordinates of the center (in pixels). If ``center=None``, the image is rotated around
+        its center, i.e. ``center=(rows / 2 - 0.5, cols / 2 - 0.5)``.
     mod : int
         Fold symmetry number. Valid numbers must be a divisor of 360.
+    mask : `~numpy.ndarray` or None, optional
+        Mask of `image`. The mask should evaluate to `True`
+        (or 1) on invalid pixels. If None (default), no mask
+        is used.
+    
 
     Returns
     -------
-    out : `~numpy.ndarray`
+    out : `~numpy.ndarray`, dtype float
+        Averaged image.
 
     Raises
     ------
     ValueError
         If `mod` is not a divisor of 360 deg.
+    
+    See also
+    --------
+    skimage.transform.rotate : Rotate images by interpolation.
     """
     if (360 % mod) != 0:
-        raise ValueError('Rotational symmetry of {} is not valid.'.format(mod))
-
-    im = np.asarray(im)
+        raise ValueError('{}-fold rotational symmetry is not valid (not a divisor of 360).'.format(mod))
     angles = range(0, 360, int(360/mod))
 
-    kwargs.update({'preserve_range': True})
-    return sum(rotate(im, angle, center = center, **kwargs) for angle in angles)/len(angles)
+    im = np.array(im, dtype = np.float, copy = True)
+
+    if mask is not None:
+       im[mask] = np.nan
+
+    rotate_kwargs = {'mode': 'constant', 'preserve_range': True}
+    rotate_kwargs.update(kwargs)
+
+    stack = np.dstack([rotate(im, angle, center = center, **rotate_kwargs) for angle in angles])
+    avg = np.nanmean(stack, axis = 2)
+    return np.nan_to_num(avg)
+
+def nfold_symmetry(*args, **kwargs):
+    warn('nfold_symmetry() is deprecated. Please use nfold in \
+          the future, as it supports more features.', DeprecationWarning)
+    return nfold(*args, **kwargs)

skued/image/tests/test_streaming.py

@@ -7,7 +7,8 @@
 from scipy.stats import sem as scipy_sem
 from skimage import data
 
-from .. import ialign, iaverage, isem, shift_image, last
+from .. import ialign, iaverage, isem, shift_image
+from ... import last
 
 
 class TestIAlign(unittest.TestCase):

skued/image/tests/test_symmetry.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 import numpy as np
-from .. import nfold_symmetry
+from .. import nfold
 import unittest
 
 np.random.seed(23)
@@ -10,7 +10,7 @@ class TestNFoldSymmetry(unittest.TestCase):
 	def test_trivial(self):
 		""" Test nfold_symmetry averaging on trivial array """
 		im = np.zeros( (256, 256) )
-		rot = nfold_symmetry(im, center = (128, 128), mod = 3)
+		rot = nfold(im, mod = 3)
 		self.assertTrue(np.allclose(rot, im))
 
 	# TODO: test preserve_range = False has not effect
@@ -19,7 +19,48 @@ def test_valid_mod(self):
 		""" Test the the N-fold symmetry argument is valid """
 		im = np.empty( (128, 128) )
 		with self.assertRaises(ValueError):
-			nfold_symmetry(im, center = (64,64), mod = 1.7)
+			nfold(im, mod = 1.7)
+
+	def test_mask(self):
+		""" Test that nfold_symmetry() works correctly with a mask """
+		im = np.zeros((128, 128), dtype = np.int)
+		mask = np.zeros_like(im, dtype = np.bool)
+
+		im[0:20] = 1
+		mask[0:20] = True
+
+		rot = nfold(im, mod = 2, mask = mask)
+		self.assertTrue(np.allclose(rot, np.zeros_like(rot)))
+
+	def test_no_side_effects(self):
+		""" Test that nfold() does not modify the input image and mask """
+		im = np.empty((128, 128), dtype = np.float)
+		mask = np.zeros_like(im, dtype = np.bool)
+
+		im.setflags(write = False)
+		mask.setflags(write = False)
+
+		rot = nfold(im, center = (67, 93),mod = 3, mask = mask)
+
+	def test_output_range(self):
+		""" Test that nfold() does not modify the value range """
+		im = 1000*np.random.random(size = (256, 256))
+		mask = np.random.choice([True, False], size = im.shape)
+
+		rot = nfold(im, center = (100, 150), mod = 5, mask = mask)
+
+		self.assertLessEqual(rot.max(), im.max())
+		# In the case of a mask that overlaps with itself when rotated,
+		# the average will be zero due to nan_to_num
+		self.assertGreaterEqual(rot.min(), min(im.min(), 0))
+
+	def test_mod_1(self):
+		""" Test that nfold(mod = 1) returns an unchanged image, except
+		perhaps for a cast to float """
+		im = 1000*np.random.random(size = (256, 256))
+		rot = nfold(im, mod = 1)
+		self.assertTrue(np.allclose(im, rot))
+
 
 if __name__ == '__main__':
 	unittest.main()

skued/iter_utils.py

@@ -0,0 +1,121 @@
+# -*- coding: utf-8 -*-
+"""
+Iterator/Generator utilities 
+============================
+"""
+from collections import deque
+from itertools import islice, count
+
+def chunked(iterable, chunksize = 1):
+	"""
+	Generator yielding multiple iterables of length 'chunksize'.
+
+	Parameters
+	----------
+	iterable : iterable
+		Iterable to be chunked. 
+	chunksize : int, optional
+		Chunk size. 
+
+	Yields
+	------
+	chunk : iterable
+		Iterable of size `chunksize`. In special case of iterable not being
+		divisible by `chunksize`, the last `chunk` might be smaller.
+	"""
+	# This looks ridiculously simple now, 
+	# but I didn't always know about itertools
+	iterable = iter(iterable)
+
+	next_chunk = tuple(islice(iterable, chunksize))
+	while next_chunk:	
+		yield next_chunk
+		next_chunk = tuple(islice(iterable, chunksize))
+
+def linspace(start, stop, num, endpoint = True):
+	""" 
+	Generate linear space. This is sometimes more appropriate than
+	using `range`.
+
+	Parameters
+	----------
+	start : float
+		The starting value of the sequence.
+	stop : float
+		The end value of the sequence.
+	num : int
+		Number of samples to generate.
+	endpoint : bool, optional
+		If True (default), the endpoint is included in the linear space.
+
+	Yields
+	------
+	val : float
+
+	See also
+	--------
+	numpy.linspace : generate linear space as a dense array.
+	"""
+	# If endpoint are to be counted in,
+	# step does not count the last yield
+	if endpoint:
+		num -= 1
+
+	step = (stop - start)/num
+
+	val = start
+	for _ in range(num):
+		yield val
+		val += step
+
+	if endpoint:
+		yield stop
+
+def multilinspace(start, stop, num, endpoint = True):
+	""" 
+	Generate multilinear space, for joining the values in two iterables.
+
+	Parameters
+	----------
+	start : iterable of floats
+		The starting value. This iterable will be consumed.
+	stop : iterable of floats
+		The end value. This iterable will be consumed.
+	num : int
+		Number of samples to generate.
+	endpoint : bool, optional
+		If True (default), the endpoint is included in the linear space.
+
+	Yields
+	------
+	val : tuple
+		Tuple of the same length as start and stop
+
+	Examples
+	--------
+	>>> multispace = multilinspaces(start = (0, 0), stop = (1, 1), num = 4, endpoint = False)
+	>>> print(list(multispace))
+	[(0, 0), (0.25, 0.25), (0.5, 0.5), (0.75, 0.75)]
+
+	See also
+	--------
+	linspace : generate a linear space between two numbers
+	"""	
+	start, stop = tuple(start), tuple(stop)
+	if len(start) != len(stop):
+		raise ValueError('start and stop must have the same length')
+
+	spaces = tuple(linspace(a, b, num = num, endpoint = endpoint) for a, b in zip(start, stop))
+	yield from zip(*spaces)
+
+def last(stream):
+	""" 
+	Retrieve the last item from a stream/iterator. Generators are consumed. 
+	If empty stream, returns None. 
+	"""
+	# Wonderful idea from itertools recipes
+	# https://docs.python.org/3.6/library/itertools.html#itertools-recipes
+	try:
+		return deque(stream, maxlen = 1)[0]
+	except IndexError:	# Empty stream
+		return None

skued/parallel.py

@@ -1,30 +1,17 @@
 # -*- coding: utf-8 -*-
 """
-Parallelization utilities.
+Parallelization utilities
+=========================
 
 Functional programming-style `map` and `reduce` procedures are easily
 parallelizable. The speed gain of parallelization can offset the
 cost of spawning multiple processes for large iterables.
 """
+import multiprocessing as mp
 from collections.abc import Sized
 from functools import partial, reduce
-import multiprocessing as mp
-
-def chunked(iterable, chunksize = 1):
-	"""
-	Generator yielding multiple iterables of length 'chunksize'.
 
-	Parameters
-	----------
-	iterable : iterable
-		Must be a sized iterable, otherwise the iterable is consumed.
-	chunksize : int, optional
-	"""
-	if not isinstance(iterable, Sized):
-		iterable = tuple(iterable)
-	length = len(iterable)
-	for ndx in range(0, length, chunksize):
-		yield iterable[ndx:min(ndx + chunksize, length)]
+from .iter_utils import chunked
 
 def preduce(func, iterable, args = tuple(), kwargs = dict(), processes = None):
 	"""