Merge pull request pyxem#1058 from viljarjf/azimuthal-integral-off-by…

…-one-bugfix

Azimuthal integral off-by-one bugfix

CHANGELOG.rst

@@ -18,6 +18,7 @@ Fixed
 - Documentation fixes and improvement. (#1028)
 - Fixed bug with flattening diffraction Vectors when there are different scales (#1024)
 - Fixed intersphinx links and improved api documentation (#1056)
+- Fix an off-by-one error in the :meth:`pyxem.signals.Diffraction2D.get_azimuthal_integral2d` (#1058)
 
 Added
 -----

examples/standards/README.rst

@@ -0,0 +1,3 @@
+Standards
+=========
+Below is a gallery of examples showing different standards in pyxem

examples/standards/pixel_coodinates.py

@@ -0,0 +1,65 @@
+"""
+====================
+Coordinates in Pyxem
+====================
+
+Pyxem is flexible in how it handles coordinates for a diffraction pattern.
+
+There are three main ways to handle coordinates in Pyxem:
+
+1. Pixel coordinates
+2. Calibrated Coordinates with evenly spaced axes
+3. Calibrated Coordinates with unevenly spaced axes (e.g. corrected for the Ewald sphere)
+"""
+
+import pyxem as pxm
+from skimage.morphology import disk
+
+
+s = pxm.signals.Diffraction2D(disk((10)))
+s.calibrate.center = None
+print(s.calibrate.center)
+
+# %%
+
+s.plot(axes_ticks=True)
+# %%
+
+# From the plot above you can see that hyperspy automatically sets the axes ticks to be centered
+# on each pixel. This means that for a 21x21 pixel image, the center is at (-10, -10) in pixel coordinates.
+# if we change the scale using the calibrate function it will automatically adjust the center.  Here it is
+# now (-1, -1)
+
+s.calibrate.scale = 0.1
+s.calibrate.units = "nm$^{-1}$"
+s.plot(axes_ticks=True)
+print(s.calibrate.center)
+
+
+# %%
+
+# Azimuthal Integration
+# ---------------------
+#
+# Now if we do integrate this dataset it will choose the appropriate center based on the center pixel.
+
+az = s.get_azimuthal_integral2d(npt=30)
+az.plot()
+
+# %%
+
+# Non-Linear Axes
+# ---------------
+#
+# Now consider the case where we have non-linear axes. In this case the center is still (10,10)
+# but things are streatched based on the effects of the Ewald Sphere.
+
+s.calibrate.beam_energy = 200
+s.calibrate.detector(pixel_size=0.1, detector_distance=3)
+print(s.calibrate.center)
+s.plot()
+
+az = s.get_azimuthal_integral2d(npt=30)
+az.plot()
+
+# %%

pyxem/tests/signals/test_diffraction2d.py

@@ -570,6 +570,43 @@ def test_internal_azimuthal_integration(self, ring):
         ring_sum = np.sum(az.data, axis=1)
         assert ring_sum.shape == (40,)
 
+    @pytest.mark.parametrize(
+        "corner",
+        [
+            (0, 0),
+            (0, -1),
+            (-1, 0),
+            (-1, -1),
+        ],
+    )
+    @pytest.mark.parametrize(
+        "shape",
+        [
+            (10, 10),  # Even square
+            (9, 9),  # Odd square
+            (4, 6),  # Even
+            (5, 9),  # Odd
+            (5, 10),  # Odd and even
+            (10, 5),  # Even and odd
+        ],
+    )
+    def test_internal_azimuthal_integration_data_range(self, corner, shape):
+        # Test that the edges of the cartesian data are included in the polar transform
+
+        max_val = 10
+        data = np.zeros(shape)
+        data[corner] = max_val
+
+        signal = Diffraction2D(data)
+
+        # Reality check
+        assert np.allclose(np.nanmax(signal.data), max_val)
+
+        # Use mean=True to conserve values
+        pol = signal.get_azimuthal_integral2d(npt=20, mean=True)
+
+        assert np.allclose(np.nanmax(pol.data), max_val)
+
 
 class TestPyFAIIntegration:
     @pytest.fixture

pyxem/tests/utils/test_calibration_utils.py

@@ -90,16 +90,20 @@ def test_set_failure(self, calibration):
             calibration.detector(pixel_size=0.1, detector_distance=1)
         calibration.beam_energy = 200
         calibration.detector(pixel_size=0.1, detector_distance=1)
+        # The center, in pixel coordinates, is (4.5, 4.5)
+        # When using a detector, this gets rounded down to 4
+        assert calibration.center == [4, 4]
         calibration.detector(
             pixel_size=0.1, detector_distance=1, beam_energy=200, units="k_nm^-1"
         )
         assert calibration.flat_ewald is False
+        assert calibration.center == [4, 4]
         with pytest.raises(ValueError):
             calibration(scale=0.01)
         assert calibration.scale is None
         with pytest.raises(ValueError):
             calibration(center=(5, 5))
-        assert calibration.center == [5, 5]
+        assert calibration.center == [4, 4]
 
         with pytest.raises(ValueError):
             calibration.detector(pixel_size=0.1, detector_distance=1, units="nm^-1")

pyxem/utils/_azimuthal_integrations.py

@@ -190,7 +190,7 @@ def _slice_radial_integrate1d(
     return ans
 
 
-def _get_factors(control_points, slices, axes):
+def _get_factors(control_points, slices, pixel_extents):
     """This function takes a set of control points (vertices of bounding polygons) and
     slices (min and max indices for each control point) and returns the factors for
     each slice. The factors are the area of the intersection of the polygon and the
@@ -199,14 +199,23 @@ def _get_factors(control_points, slices, axes):
     factors = []
     factors_slice = []
     start = 0
+    # unpack extents
+    x_extent, y_extent = pixel_extents
+    x_ext_left, x_ext_right = x_extent
+    y_ext_left, y_ext_right = y_extent
     for cp, sl in zip(control_points, slices):
         p = Polygon(cp)
         x_edges = list(range(sl[0], sl[2]))
         y_edges = list(range(sl[1], sl[3]))
         boxes = []
         for i, x in enumerate(x_edges):
             for j, y in enumerate(y_edges):
-                b = box(axes[0][x], axes[1][y], axes[0][x + 1], axes[1][y + 1])
+                b = box(
+                    x_ext_left[x],
+                    y_ext_left[y],
+                    x_ext_right[x],
+                    y_ext_right[y],
+                )
                 boxes.append(b)
         factors += list(
             shapely.area(shapely.intersection(boxes, p)) / shapely.area(boxes)

pyxem/utils/calibration.py

@@ -18,7 +18,6 @@
 
 """Utils for calibrating Diffraction Patterns."""
 
-
 import numpy as np
 import json
 
@@ -228,19 +227,20 @@ def detector(
             raise ValueError(
                 f"Unit {units} not recognized. Must be one of {unit_factors.keys()}"
             )
+
         if center is None:
-            center = np.array(self.signal.axes_manager.signal_shape) / 2
-        x_pixels = np.arange(-center[0], self.shape[0] - center[0])
-        y_pixels = np.arange(-center[1], self.shape[1] - center[1])
+            center = [(ax.size - 1) / 2 for ax in self.axes]
 
-        x_pixels *= pixel_size
-        y_pixels *= pixel_size
+        def translate_pixel_coords(px: np.ndarray) -> np.ndarray:
+            coord = pixel_size * px
+            angle = np.arctan(coord / detector_distance)
+            return np.sin(angle) * (1 / wavelength) * unit_factors[units]
 
-        x_angles = np.arctan(x_pixels / detector_distance)
-        y_angles = np.arctan(y_pixels / detector_distance)
+        x_pixels = np.arange(self.shape[0]) - center[0]
+        y_pixels = np.arange(self.shape[1]) - center[1]
 
-        x_axes = np.sin(x_angles) * (1 / wavelength) * unit_factors[units]
-        y_axes = np.sin(y_angles) * (1 / wavelength) * unit_factors[units]
+        x_axes = translate_pixel_coords(x_pixels)
+        y_axes = translate_pixel_coords(y_pixels)
 
         for ax, axis in zip(self.signal.axes_manager.signal_axes, [x_axes, y_axes]):
             if isinstance(ax, UniformDataAxis):
@@ -271,6 +271,51 @@ def __repr__(self):
     def axes(self):
         return [ax.axis for ax in self.signal.axes_manager.signal_axes][::-1]
 
+    @property
+    def pixel_extent(self):
+        """Return an array with axes [x/y, left/right, pixel_extent], as follows:
+        [
+            # x axis
+            [
+                # left
+                [boundary 1, boundary 2 ...],
+                # right
+                [boundary 1, boundary 2 ...],
+            ],
+            # y axis
+            [
+                # left
+                [boundary 1, boundary 2 ...],
+                # right
+                [boundary 1, boundary 2 ...],
+            ],
+        ]
+        """
+        if self.flat_ewald:
+            left_scales = self.scale
+            right_scales = self.scale
+        else:
+            x_scales = self.axes[0][1:] - self.axes[0][:-1]
+            x_scales = np.pad(x_scales, 1, mode="edge")
+            y_scales = self.axes[1][1:] - self.axes[1][:-1]
+            y_scales = np.pad(y_scales, 1, mode="edge")
+            left_scales = [
+                x_scales[:-1],
+                y_scales[:-1],
+            ]
+            right_scales = [
+                x_scales[1:],
+                y_scales[1:],
+            ]
+
+        extents = []
+        for ax, left_scale, right_scale in zip(self.axes, left_scales, right_scales):
+            left = ax - left_scale / 2
+            right = ax + right_scale / 2
+            extent = np.stack((left, right))
+            extents.append(extent)
+        return extents
+
     def get_slices2d(self, npt, npt_azim, radial_range=None):
         """Get the slices and factors for some image that can be used to
         slice the image for 2d integration.
@@ -293,7 +338,13 @@ def _get_radial_range(self, radial_range=None):
         if radial_range is None:
             from itertools import combinations
 
-            edges = np.reshape([[ax[0] ** 2, ax[-1] ** 2] for ax in self.axes], -1)
+            edges = np.reshape(
+                [
+                    [ax_ext[0][0] ** 2, ax_ext[1][-1] ** 2]
+                    for ax_ext in self.pixel_extent
+                ],
+                -1,
+            )
             max_range = np.max(
                 np.power(np.sum(list(combinations(edges, 2)), axis=1), 0.5)
             )
@@ -359,35 +410,42 @@ def _get_slices_and_factors(self, npt, npt_azim, radial_range):
         # get the points which bound each azimuthal pixel
         control_points = _get_control_points(npt, npt_azim, radial_range, self.affine)
 
-        # get the min and max indices for each control point using the axes
+        # get the min and max indices for each control point using the
+        pixel_ext_x, pixel_ext_y = self.pixel_extent
         min_x = (
             np.min(
-                np.searchsorted(self.axes[0], control_points[:, :, 0], side="left"),
+                np.searchsorted(
+                    pixel_ext_x[1, :], control_points[:, :, 0], side="left"
+                ),
                 axis=1,
             ).astype(int)
             - 1
         )
         min_y = (
             np.min(
-                np.searchsorted(self.axes[1], control_points[:, :, 1], side="left"),
+                np.searchsorted(
+                    pixel_ext_y[1, :], control_points[:, :, 1], side="left"
+                ),
                 axis=1,
             ).astype(int)
             - 1
         )
 
         max_x = np.max(
-            np.searchsorted(self.axes[0], control_points[:, :, 0], side="right"), axis=1
+            np.searchsorted(pixel_ext_x[0, :], control_points[:, :, 0], side="right"),
+            axis=1,
         ).astype(int)
         max_y = np.max(
-            np.searchsorted(self.axes[1], control_points[:, :, 1], side="right"), axis=1
+            np.searchsorted(pixel_ext_y[0, :], control_points[:, :, 1], side="right"),
+            axis=1,
         ).astype(int)
         # Note that if a point is outside the range of the axes it will be set to the
         # maximum value of the axis+1 and 0 if it is below the minimum value of the axis.
         slices = np.array(
             [[mx, my, mxx, myy] for mx, my, mxx, myy in zip(min_x, min_y, max_x, max_y)]
         )
-        max_y_ind = len(self.axes[1]) - 1
-        max_x_ind = len(self.axes[0]) - 1
+        max_y_ind = len(self.axes[1])
+        max_x_ind = len(self.axes[0])
 
         # set the slices to be within the range of the axes.  If the entire slice is outside
         # the range of the axes then set the slice to be the maximum value of the axis
@@ -401,7 +459,7 @@ def _get_slices_and_factors(self, npt, npt_azim, radial_range):
         slices[slices[:, 2] < 0, 2] = 0
         slices[slices[:, 3] < 0, 3] = 0
 
-        factors, factors_slice = _get_factors(control_points, slices, self.axes)
+        factors, factors_slice = _get_factors(control_points, slices, self.pixel_extent)
         return slices, factors, factors_slice
 
     @property
@@ -424,7 +482,7 @@ def center(self, center=None):
             )
         if center is None:
             for ax in self.signal.axes_manager.signal_axes:
-                ax.offset = -ax.scale * (ax.size / 2)
+                ax.offset = -ax.scale * ((ax.size - 1) / 2)
         else:
             for ax, off in zip(self.signal.axes_manager.signal_axes, center):
                 ax.offset = -off * ax.scale