Merge pull request #9 from LaurentRDC/develop

Version 0.4.2

appveyor.yml

@@ -12,7 +12,7 @@ environment:
       PYTHON_ARCH: "64" # needs to be set for CMD_IN_ENV to succeed. If a mix
                         # of 32 bit and 64 bit builds are needed, move this
                         # to the matrix section.
-      TEST_CMD: "python -m unittest discover"
+      TEST_CMD: "python -m unittest discover --verbose"
 
   matrix:
 

setup.py

@@ -21,9 +21,6 @@
 with open('README.rst') as f:
     readme = f.read()
 
-with open('CHANGELOG.rst') as f:
-    changes = f.read()
-
 with open('requirements.txt') as f:
     requirements = [line for line in f.read().split('\n') if len(line.strip())]
 
@@ -35,7 +32,7 @@
     setup(
         name = 'scikit-ued',
         description = 'Collection of algorithms and functions for ultrafast electron diffraction',
-        long_description = '\n\n'.join([readme, changes]),
+        long_description = '\n\n'.join(readme),
         license = LICENSE,
         url = 'http://scikit-ued.readthedocs.io',
         download_url = 'http://github.com/LaurentRDC/scikit-ued',
@@ -52,7 +49,8 @@
         zip_safe = False,
         classifiers = ['Environment :: Console',
                        'Intended Audience :: Science/Research',
-                       'Topic :: Scientific'
+                       'Topic :: Scientific/Engineering',
+                       'Topic :: Scientific/Engineering :: Physics',
                        'License :: OSI Approved :: MIT License',
                        'Natural Language :: English',
                        'Operating System :: OS Independent',

skued/__init__.py

@@ -2,7 +2,7 @@
 __author__ = 'Laurent P. René de Cotret'
 __email__ = 'laurent.renedecotret@mail.mcgill.ca'
 __license__ = 'MIT'
-__version__ = '0.4.1'
+__version__ = '0.4.2' # TODO: automatic versioning
 
 from .array_utils import repeated_array
 from .parallel import pmap, preduce

skued/baseline/algorithms.py

@@ -53,8 +53,8 @@ def _iterative_baseline(array, max_iter, mask, background_regions, axes, approx_
 
 	# Readjust size for odd input signals
 	# If axis was padded, trim
-	if background.shape != original_shape:
-		for axis in axes:
+	for axis in axes:
+		while background.shape[axis] > original_shape[axis]:
 			background = np.swapaxes(np.swapaxes(background, 0, axis)[:-1], 0, axis)
 	return background
 
@@ -157,7 +157,7 @@ def _dwt_approx_rec(array, level, wavelet, mode, axis):
 	reconstructed = pywt.waverec([app_coeffs] + [None,]*len(det_coeffs), wavelet = wavelet, mode = 'constant', axis = axis)
 
 	# Sometimes pywt.waverec returns a signal that is longer than the original signal
-	if reconstructed.shape[axis] > array.shape[axis]:
+	while reconstructed.shape[axis] > array.shape[axis]:
 		reconstructed = np.swapaxes(np.swapaxes(reconstructed, 0, axis)[:array.shape[axis]], 0, axis)
 	return  reconstructed
 
@@ -225,9 +225,9 @@ def _dwt_approx_rec2(array, level, wavelet, mode, axis):
 								  wavelet = wavelet, mode = mode, axes = axis)
 
 	# Sometimes pywt.waverec returns a signal that is longer than the original signal
-	for ax in axis:
-		if reconstructed.shape[ax] > array.shape[ax]:
-			reconstructed = np.swapaxes(np.swapaxes(reconstructed, 0, ax)[:-1], 0, ax)
+	for ax in range(reconstructed.ndim):
+		while reconstructed.shape[ax] > array.shape[ax]:
+			reconstructed = np.swapaxes(np.swapaxes(reconstructed, 0, ax)[:-1], ax, 0)
 	return reconstructed
 
 def baseline_dt(array, max_iter, level = None, first_stage = 'sym6', wavelet = 'qshift1', 

skued/baseline/tests/test_algorithms.py

@@ -23,17 +23,17 @@ def test_approx_rec(self):
 			self.assertSequenceEqual(rec_arr.shape, arr.shape)
 
 		with self.subTest('1D along axis'):
-			arr2 = np.random.random(size = (21, 52))
+			arr2 = np.random.random(size = (101, 104))
 			rec_arr2 = _dwt_approx_rec(arr2, level = 2, wavelet = 'sym4', mode = 'constant', axis = 1)
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
 		with self.subTest('2D shape'):
-			arr2 = np.random.random(size = (22,52))
+			arr2 = np.random.random(size = (102,104))
 			rec_arr2 = _dwt_approx_rec2(arr2, level = 2, wavelet = 'sym6', mode = 'constant', axis = (-2, -1))
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
 		with self.subTest('2D along axes'):
-			arr2 = np.random.random(size = (22,52, 5))
+			arr2 = np.random.random(size = (102,94,25))
 			rec_arr2 = _dwt_approx_rec2(arr2, level = 2, wavelet = 'sym6', mode = 'constant', axis = (0, 1))
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
@@ -88,6 +88,29 @@ def test_trivial_case(self):
 		""" The baseline for an array of zeros should be zeros """
 		arr = np.zeros_like(self.arr)
 		self.assertTrue(np.allclose(arr, baseline_dt(arr, max_iter = 10, level = None)))
+
+	def test_final_shape(self):
+		""" Test that baseline_dt returns an array of the same shape as input """
+
+		with self.subTest('1D baseline of odd length'):
+			arr = np.random.random(size = (101,))
+			b = baseline_dt(arr, max_iter = 10, level = None)
+			self.assertSequenceEqual(arr.shape, b.shape)
+
+		with self.subTest('1D baseline of even length'):
+			arr = np.random.random(size = (100,))
+			b = baseline_dt(arr, max_iter = 10, level = None)
+			self.assertSequenceEqual(arr.shape, b.shape)
+
+		with self.subTest('baseline along axis of odd length'):
+			arr = np.random.random(size = (101, 113))
+			b = baseline_dt(arr, max_iter = 10, level = None, axis = 0)
+			self.assertSequenceEqual(arr.shape, b.shape)
+
+		with self.subTest('baseline along axis of even length'):
+			arr = np.random.random(size = (102, 104))
+			b = baseline_dt(arr, max_iter = 10, level = None, axis = 1)
+			self.assertSequenceEqual(arr.shape, b.shape)
 
 	def test_2d_along_axis(self):
 		""" Test that iterating over array rows and baseline_dt along axis are equivalent """

skued/image_analysis/__init__.py

@@ -2,4 +2,5 @@
 """ Diffraction image analysis """
 
 from .powder import angular_average, powder_center
-from .alignment import align, shift_image
+from .alignment import align, shift_image
+from .symmetry import nfold_symmetry

skued/image_analysis/powder.py

@@ -1,13 +1,10 @@
 # -*- coding: utf-8 -*-
 """
-Image manipulation involving symmetry
+Image manipulation of powder diffraction
 """
-import matplotlib.pyplot as plt
 import numpy as np
-from scipy.signal import fftconvolve, find_peaks_cwt
-from skimage.filters import gaussian, threshold_local
-
-import matplotlib.pyplot as plt
+from scipy.signal import fftconvolve
+from skimage.filters import threshold_local
 
 _PC_CACHE = dict()
 def powder_center(image):
@@ -54,7 +51,7 @@ def powder_center(image):
 	full_center = np.array(np.unravel_index(np.argmax(corr), dims = corr.shape))
 	return tuple(full_center/2)
 
-def angular_average(image, center, mask = None, extras = None):
+def angular_average(image, center, mask = None, extras = None, angular_bounds = None):
 	"""
 	This function returns an angularly-averaged pattern computed from a diffraction pattern, 
 	e.g. polycrystalline diffraction.
@@ -69,9 +66,13 @@ def angular_average(image, center, mask = None, extras = None):
 		Evaluates to True on invalid elements of array.
 	extras : dict-like or None, optional
 		if not None, this dict-like object will be updated with: 
-
 			extras['error'] : `~numpy.ndarray`
 				standard error in mean across radii.
+	angular_bounds : 2-tuple or None, optional
+		If not None, the angles between first and second elements of `angular_bounds`
+		(inclusively) will be used for the average. Angle bounds are specified in degrees.
+		0 degrees is defined as the positive x-axis. Angle bounds outside [0, 360) are mapped back
+		to [0, 360).
 
 	Returns
 	-------
@@ -80,8 +81,6 @@ def angular_average(image, center, mask = None, extras = None):
 	average : `~numpy.ndarray`, shape (N,)
 		Angular-average of the array.
 	"""
-	#TODO: N dimensions. The only problem is the generation of an arbitrary
-	#	   meshgrid on the image.
 	#TODO: axis parameter? so batches of angular averages can be computed
 	#	   at the same time.
 	image = np.asarray(image, dtype = np.float)
@@ -96,6 +95,12 @@ def angular_average(image, center, mask = None, extras = None):
 					   np.arange(0, image.shape[1]))
 	R = np.rint(np.sqrt( (X - xc)**2 + (Y - yc)**2 )).astype(np.int)
 
+	if angular_bounds:
+		mi, ma = angular_bounds
+		mi, ma = mi % 360, ma % 360
+		angles = np.rad2deg(np.arctan2(Y - yc, X - xc))
+		mask[np.logical_not(np.logical_and(mi <= angles, angles <= ma))] = True
+
 	valid = np.logical_not(mask)
 	R_v, image_v = R[valid].ravel(), image[valid].ravel()
 	px_bin = np.bincount(R_v, weights = image_v)

skued/image_analysis/symmetry.py

@@ -0,0 +1,40 @@
+# -*- coding: utf-8 -*-
+"""
+Image manipulation involving symmetry
+"""
+from functools import partial
+from skimage.transform import rotate
+import numpy as np
+
+# TODO: out parameter?
+def nfold_symmetry(im, center, mod, **kwargs):
+    """ 
+    Returns an images averaged according to n-fold rotational symmetry.
+    Keyword arguments are passed to skimage.transform.rotate()
+
+    Parameters
+    ----------
+    im : array_like, ndim 2
+        Image to be averaged.
+	center : array_like, shape (2,)
+		coordinates of the center (in pixels).
+    mod : int
+        Fold symmetry number. Valid numbers must be a divisor of 360.
+
+    Returns
+    -------
+    out : `~numpy.ndarray`
+
+    Raises
+    ------
+    ValueError
+        If `mod` is not a divisor of 360 deg.
+    """
+    if (360 % mod) != 0:
+        raise ValueError('Rotational symmetry of {} is not valid.'.format(mod))
+
+    im = np.asarray(im)
+    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)

skued/image_analysis/tests/test_alignment.py

@@ -6,8 +6,6 @@
 from skimage.filters import gaussian
 from skimage.transform import rotate
 
-import matplotlib.pyplot as plt
-
 from .. import shift_image, align
 
 class TestAlign(unittest.TestCase):

skued/image_analysis/tests/test_powder.py

@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-import matplotlib.pyplot as plt
+
 import numpy as np
 from .. import angular_average, powder_center
 import unittest
@@ -79,6 +79,35 @@ def test_ring(self):
 		self.assertEqual(intensity.max(), image.max())
 		self.assertSequenceEqual(radius.shape, intensity.shape)
 
+	def test_angular_bounds(self):
+		""" Test angular_average with a restrictive angular_bounds argument """
+		image = np.zeros(shape = (256, 256), dtype = np.float)
+		center = (image.shape[0]/2, image.shape[1]/2)
+		xc, yc = center
+
+		# Create an image with a wide ring
+		extent = np.arange(0, image.shape[0])
+		xx, yy = np.meshgrid(extent, extent)
+		rr = np.sqrt((xx - xc)**2 + (yy - yc)**2)
+		angles = np.rad2deg(np.arctan2(yy - yc, xx - xc))
+		image[np.logical_and(0 <= angles, angles <= 60)] = 1
+
+		with self.subTest('Outside angle bounds'):
+			radius, intensity = angular_average(image, center, angular_bounds = (0, 60))
+			self.assertTrue(np.allclose(intensity, np.ones_like(intensity)))
+
+		with self.subTest('Overlapping bounds'):
+			radius, intensity = angular_average(image, center, angular_bounds = (15, 75))
+			self.assertFalse(np.allclose(intensity, np.ones_like(intensity)))
+
+		with self.subTest('Inside angle bounds'):
+			radius, intensity = angular_average(image, center, angular_bounds = (60, 360))
+			self.assertTrue(np.allclose(intensity, np.zeros_like(intensity)))
+
+		with self.subTest('Inside angle bounds with 360deg rollover'):
+			radius, intensity = angular_average(image, center, angular_bounds = (60 + 360, 360 + 360))
+			self.assertTrue(np.allclose(intensity, np.zeros_like(intensity)))
+
 	def test_ring_with_mask(self):
 		""" Test angular_average on an image with a wide ring """
 		image = np.zeros(shape = (256, 256), dtype = np.float)

skued/image_analysis/tests/test_symmetry.py

@@ -0,0 +1,25 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+from .. import nfold_symmetry
+import unittest
+
+np.random.seed(23)
+
+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)
+		self.assertTrue(np.allclose(rot, im))
+
+	# TODO: test preserve_range = False has not effect
+
+	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)
+
+if __name__ == '__main__':
+	unittest.main()