Added new Bragg peak functional and auto Brillouin zone determination (…

…#42)

* pyqtgraph moved home of PyQT5 functionality; added catch all for diff versions

* Added new functionality to namespace

* Added auto Brilluoin zone generation

* New Bragg peak determination with 2D persistence

* added correct behaviour for masking

* added test for new Bragg peak detection method

* added test for BZ generation via determined BZ consistency

* updated documentation to reflect new functionalities

* typo in cl argument

* added tutorial for new Bragg peak finding and BZ generation

* Updated documentation examples for BZ generation

* Return all scattering vectors not in the center BZ

* added documentation for equiv scat vec meth

* fixed typo in name Brilluoin

* Formatting

* Formatting

CHANGELOG.rst

@@ -1,6 +1,11 @@
 Changelog
 =========
 
+Release 2.1.11
+--------------
+
+* Added new Bragg peak determination functionality to combat new datasets with sub-optimal signal-to-noise ratio. Further added Brillouin zone determination based on Bragg peak locations. 
+
 Release 2.1.10
 --------------
 

docs/api.rst

@@ -220,6 +220,8 @@ Single-crystal indexing
     :nosignatures:
 
     bragg_peaks
+    bragg_peaks_persistence
+    brillouin_zones
 
 ==========
 Simulation

docs/tutorials/image.rst

@@ -369,6 +369,88 @@ We can plot the result:
 
 How cool is that!
 
+We can also use the method of 2D peak topological persistence to find Bragg peaks. Using the 
+:func:`bragg_peaks_persistence` function, we determine the peaks via:
+
+	>>> from skued import diffread, bragg_peaks_persistence
+	>>> 
+	>>> im = diffread('docs/tutorials/data/ab_initio_monolayer_mos2.tif')
+	>>>
+	>>> peaks, bd, bd_indices, persistence = bragg_peaks_persistence(im, prominence=0.1)
+
+We plot the results. The left plot shows that the standard :func:`bragg_peaks` finds no peaks,
+while :func:`bragg_peaks_persistence` finds all of them, given on the right :
+
+.. plot:: 
+
+	from skued import diffread, bragg_peaks_persistence, bragg_peaks
+	import matplotlib.pyplot as plt
+	import numpy as np
+	from matplotlib.patches import Circle
+
+	im = diffread("data/ab_initio_monolayer_mos2.tif")
+
+	fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3))
+
+	peaks = np.array(bragg_peaks(im, mask=np.ones_like(im, dtype=bool))).reshape(-1,2)
+	peaks_persistence, _, _, _ = bragg_peaks_persistence(im, prominence=0.1)
+	for pp, image, ax in zip([peaks, peaks_persistence], [im, im], [ax1, ax2]):
+		ax.imshow(im, origin='lower', interpolation='bicubic', vmin=-1, vmax=1)
+		for r, c in pp:
+			ax.add_patch(Circle(xy=(r, c), radius=5, ec="r", fc="none"))
+		ax.get_xaxis().set_visible(False)
+		ax.get_yaxis().set_visible(False)
+	plt.tight_layout()
+	plt.show()
+
+Finding Brillouin zones
+=======================
+
+We can also determine Brillouin zones based on the location of Bragg peaks.
+For normal incidence diffraction, such patterns are regular polygons for 
+simple crystal structures, but for tilted-incidence diffraction or macromolecules,
+the geometry becomes unclear, necessitating a different approach. The Brillouin zone,
+the Wigner-Seitz unit cell in reciprocal space, is also the Voronoi regions determined
+by the Bragg peak locations. This functionality is implemented in this package 
+by the class :class:`brillouin_zones`. 
+
+	>>> from skued import diffread, autocenter, bragg_peaks_persistence, brillouin_zones
+	>>> 
+	>>> im = diffread('docs/tutorials/data/ab_initio_monolayer_mos2.tif')
+	>>>
+	>>> peaks, bd, bd_indices, persistence = bragg_peaks_persistence(im, prominence=0.1)
+	>>> center = np.array([im.shape[0]//2, im.shape[1]//2])
+	>>> BZ = brillouin_zones(im, mask=np.ones(im.shape), peaks=peaks.astype(int), center=center.astype(int))
+
+We can even determine equivalent scattering vectors at all visible
+Brillouin zones via the :meth:`getEquivalentScatteringVectors` method, useful for averaging the signal at a given reduced wavevector
+for low SNR single-crystal datasets, visualized explicitly on the right.
+
+.. plot:: 
+
+	from skued import diffread, bragg_peaks_persistence, brillouin_zones
+	import matplotlib.pyplot as plt
+	import numpy as np
+
+	im = diffread("data/ab_initio_monolayer_mos2.tif")
+	center = np.array([im.shape[0]//2, im.shape[1]//2])
+	peaks_persistence, _, _, _ = bragg_peaks_persistence(im, prominence=0.1)
+	BZ = brillouin_zones(im, mask=np.ones(im.shape), peaks=peaks_persistence.astype(int), center=center.astype(int))
+	BZ.getVisibleBZs(symmetry=6)
+	EQUIV = BZ.getEquivalentScatteringVectors(np.array([119, 174]), use_visible=True)
+	fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3))
+	for ax in [ax1, ax2]:
+		ax.imshow(im, origin='lower', interpolation='bicubic', vmin=-1, vmax=1)
+		BZ.renderVisibleBZs(ax, **dict(edgecolor='white'))
+
+		ax.axis('off')
+
+	ax2.scatter(EQUIV[:, 0], EQUIV[:, 1], marker="x", color="r", s=20)
+
+	plt.tight_layout()
+	plt.show()
+
+
 .. _powder:
 
 Image analysis on polycrystalline diffraction patterns

skued/__init__.py

@@ -45,12 +45,14 @@
     align,
     autocenter,
     azimuthal_average,
+    brillouin_zones,
     combine_masks,
     detector_scattvectors,
     ialign,
     isnr,
     itrack_peak,
     bragg_peaks,
+    bragg_peaks_persistence,
     mask_from_collection,
     mask_image,
     nfold,

skued/array_utils.py

@@ -180,7 +180,7 @@ def cart2spherical(x, y, z):
     t : `~numpy.ndarray`
         Azimuthal coordinate in radians.
     """
-    r = np.sqrt(x**2 + y**2 + z**2)
+    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
     p = np.arctan2(y, x)
     t = np.arccos(z / r)
     return r, p, t

skued/eproperties.py

@@ -22,7 +22,7 @@ def lorentz(keV):
     ----------
     .. Kirkland 2010 Eq. 2.2
     """
-    return 1 + (elementary_charge * keV * 1e3) / (electron_mass * speed_of_light**2)
+    return 1 + (elementary_charge * keV * 1e3) / (electron_mass * speed_of_light ** 2)
 
 
 def electron_wavelength(keV):
@@ -53,7 +53,7 @@ def electron_wavelength(keV):
     wavelength_meters = (
         Planck
         * speed_of_light
-        / np.sqrt(eV * (2 * electron_mass * speed_of_light**2 + eV))
+        / np.sqrt(eV * (2 * electron_mass * speed_of_light ** 2 + eV))
     )
     return wavelength_meters * 1e10  # wavelength in angstroms
 
@@ -81,7 +81,7 @@ def electron_velocity(keV):
     .. Kirkland 2010 Eq. 2.3
     """
     eV = elementary_charge * keV * 1e3
-    m0c2 = electron_mass * speed_of_light**2
+    m0c2 = electron_mass * speed_of_light ** 2
     v_over_c = np.sqrt(eV * (eV + 2 * m0c2)) / (m0c2 + eV)
     return (speed_of_light * v_over_c) * 1e10  # speed in Angstroms
 
@@ -110,6 +110,6 @@ def interaction_parameter(keV):
     return (
         (2 * np.pi)
         / (electron_wavelength(keV) * V)
-        * (electron_mass * speed_of_light**2 + elementary_charge * V)
-        / (2 * electron_mass * speed_of_light**2 + elementary_charge * V)
+        * (electron_mass * speed_of_light ** 2 + elementary_charge * V)
+        / (2 * electron_mass * speed_of_light ** 2 + elementary_charge * V)
     )

skued/image/__init__.py

@@ -2,9 +2,10 @@
 """ Diffraction image analysis """
 
 from .alignment import align, ialign, itrack_peak
+from .brillouin import brillouin_zones
 from .calibration import detector_scattvectors, powder_calq
 from .center import autocenter
-from .indexing import bragg_peaks
+from .indexing import bragg_peaks, bragg_peaks_persistence
 from .metrics import (
     combine_masks,
     isnr,