Merge branch 'master' into master-reorg

Author: ZackAttack614

hexrd/core/gridutil.py

@@ -146,6 +146,10 @@ def cellCentroids(crd, con):
 
 @numba.njit(nogil=True, cache=True)
 def compute_areas(xy_eval_vtx, conn):
+    # NOTE: this function may return negative areas if the vertices
+    # are passed in the opposite order to the function. This happens
+    # if the beam vector is in the opposite direction (positive Z
+    # instead of the usual negative Z)
     areas = np.empty(len(conn))
     for i in range(len(conn)):
         vtx0x, vtx0y = xy_eval_vtx[conn[i, 0]]

hexrd/core/instrument/detector.py

@@ -9,6 +9,7 @@
     PHOSPHOR_DEFAULT,
 )
 from hexrd.core.instrument.physics_package import AbstractPhysicsPackage
+import numba
 import numpy as np
 
 from hexrd.core import constants as ct
@@ -1122,11 +1123,9 @@ def interpolate_nearest(self, xy, img, pad_with_nans=True):
 
     def interpolate_bilinear(
         self,
-        xy,
-        img,
-        pad_with_nans=True,
-        clip_to_panel=True,
-        on_panel: Optional[np.ndarray] = None,
+        xy: np.ndarray,
+        img: np.ndarray,
+        pad_with_nans: bool = True,
     ):
         """
         Interpolate an image array at the specified cartesian points.
@@ -1137,13 +1136,10 @@ def interpolate_bilinear(
             Array of cartesian coordinates in the image plane at which
             to evaluate intensity.
         img : array_like
-            2-dimensional image array.
+            2-dimensional image array. The shape must match (rows, cols).
         pad_with_nans : bool, optional
             Toggle for assigning NaN to points that fall off the detector.
             The default is True.
-        on_panel : np.ndarray, optional
-            If you want to skip clip_to_panel() for performance reasons,
-            just provide an array of which pixels are on the panel.
 
         Returns
         -------
@@ -1155,28 +1151,30 @@ def interpolate_bilinear(
         -----
         TODO: revisit normalization in here?
         """
+        fill_value = np.nan if pad_with_nans else 0
+        int_xy = np.full(len(xy), fill_value, dtype=float)
 
-        is_2d = img.ndim == 2
-        right_shape = img.shape[0] == self.rows and img.shape[1] == self.cols
-        assert (
-            is_2d and right_shape
-        ), "input image must be 2-d with shape (%d, %d)" % (
-            self.rows,
-            self.cols,
-        )
+        # clip away points too close to or off the detector edges
+        xy_clip, on_panel = self.clip_to_panel(xy, buffer_edges=True)
 
-        # initialize output with nans
-        if pad_with_nans:
-            int_xy = np.nan * np.ones(len(xy))
-        else:
-            int_xy = np.zeros(len(xy))
+        # Generate the interpolation dict
+        interp_dict = self._generate_bilinear_interp_dict(xy_clip)
 
-        if on_panel is None:
-            # clip away points too close to or off the edges of the detector
-            xy_clip, on_panel = self.clip_to_panel(xy, buffer_edges=True)
-        else:
-            xy_clip = xy[on_panel]
+        # Set the output and return
+        int_xy[on_panel] = _interpolate_bilinear(img, **interp_dict)
+        return int_xy
 
+    def _generate_bilinear_interp_dict(
+        self,
+        xy_clip: np.ndarray,
+    ) -> dict[str, np.ndarray]:
+        """Compute bilinear interpolation multipliers and indices for the panel
+
+        If you are going to be using the same panel settings and performing
+        interpolation on multiple images, it is advised to run this beforehand
+        to precompute the interpolation parameters, so you can use them
+        repeatedly.
+        """
         # grab fractional pixel indices of clipped points
         ij_frac = self.cartToPixel(xy_clip)
 
@@ -1196,20 +1194,24 @@ def interpolate_bilinear(
         j_ceil = j_floor + 1
         j_ceil_img = _fix_indices(j_ceil, 0, self.cols - 1)
 
-        # first interpolate at top/bottom rows
-        row_floor_int = (j_ceil - ij_frac[:, 1]) * img[
-            i_floor_img, j_floor_img
-        ] + (ij_frac[:, 1] - j_floor) * img[i_floor_img, j_ceil_img]
-        row_ceil_int = (j_ceil - ij_frac[:, 1]) * img[
-            i_ceil_img, j_floor_img
-        ] + (ij_frac[:, 1] - j_floor) * img[i_ceil_img, j_ceil_img]
-
-        # next interpolate across cols
-        int_vals = (i_ceil - ij_frac[:, 0]) * row_floor_int + (
-            ij_frac[:, 0] - i_floor
-        ) * row_ceil_int
-        int_xy[on_panel] = int_vals
-        return int_xy
+        # Compute differences between raw coordinates to use for interpolation
+        j_ceil_sub = j_ceil - ij_frac[:, 1]
+        j_floor_sub = ij_frac[:, 1] - j_floor
+        i_ceil_sub = i_ceil - ij_frac[:, 0]
+        i_floor_sub = ij_frac[:, 0] - i_floor
+
+        return {
+            # Compute interpolation multipliers for every pixel
+            'cc': j_ceil_sub * i_ceil_sub,
+            'fc': j_floor_sub * i_ceil_sub,
+            'cf': j_ceil_sub * i_floor_sub,
+            'ff': j_floor_sub * i_floor_sub,
+            # Store needed pixel indices
+            'i_floor_img': i_floor_img,
+            'j_floor_img': j_floor_img,
+            'i_ceil_img': i_ceil_img,
+            'j_ceil_img': j_ceil_img,
+        }
 
     def make_powder_rings(
         self,
@@ -2138,3 +2140,63 @@ def _row_edge_vec(rows, pixel_size_row):
 
 def _col_edge_vec(cols, pixel_size_col):
     return pixel_size_col * (np.arange(cols + 1) - 0.5 * cols)
+
+
+@numba.njit(nogil=True, cache=True)
+def _interpolate_bilinear(
+    img: np.ndarray,
+    cc: np.ndarray,
+    fc: np.ndarray,
+    cf: np.ndarray,
+    ff: np.ndarray,
+    i_floor_img: np.ndarray,
+    j_floor_img: np.ndarray,
+    i_ceil_img: np.ndarray,
+    j_ceil_img: np.ndarray,
+) -> np.ndarray:
+    # The math is faster and uses the GIL less (which is more
+    # multi-threading friendly) when we run this code in numba.
+    result = np.zeros(i_floor_img.shape[0], dtype=img.dtype)
+    on_panel_idx = np.arange(i_floor_img.shape[0])
+    _interpolate_bilinear_in_place(
+        img,
+        cc,
+        fc,
+        cf,
+        ff,
+        i_floor_img,
+        j_floor_img,
+        i_ceil_img,
+        j_ceil_img,
+        on_panel_idx,
+        result,
+    )
+    return result
+
+
+@numba.njit(nogil=True, cache=True)
+def _interpolate_bilinear_in_place(
+    img: np.ndarray,
+    cc: np.ndarray,
+    fc: np.ndarray,
+    cf: np.ndarray,
+    ff: np.ndarray,
+    i_floor_img: np.ndarray,
+    j_floor_img: np.ndarray,
+    i_ceil_img: np.ndarray,
+    j_ceil_img: np.ndarray,
+    on_panel_idx: np.ndarray,
+    output_img: np.ndarray,
+):
+    # The math is faster and uses the GIL less (which is more
+    # multi-threading friendly) when we run this code in numba.
+    # Running in-place eliminates some intermediary arrays for
+    # even faster performance.
+    for i in range(on_panel_idx.shape[0]):
+        idx = on_panel_idx[i]
+        output_img[idx] += (
+            cc[i] * img[i_floor_img[i], j_floor_img[i]] +
+            fc[i] * img[i_floor_img[i], j_ceil_img[i]] +
+            cf[i] * img[i_ceil_img[i], j_floor_img[i]] +
+            ff[i] * img[i_ceil_img[i], j_ceil_img[i]]
+        )

hexrd/core/instrument/hedm_instrument.py

@@ -2941,6 +2941,12 @@ def _extract_ring_line_positions(
         # strip relevant objects out of current patch
         vtx_angs, vtx_xys, conn, areas, xys_eval, ijs = patch
 
+        # These areas can be negative if the beam vector is in
+        # the opposite direction than it normally is in (positive
+        # Z instead of the usual negative Z). Take the absolute
+        # value of the areas to ensure they are positive.
+        areas = np.abs(areas)
+
         # need to reshape eval pts for interpolation
         xy_eval = np.vstack([xys_eval[0].flatten(), xys_eval[1].flatten()]).T
 

hexrd/core/material/unitcell.py

@@ -592,28 +592,19 @@ def CalcStar(self, v, space, applyLaue=False):
         this function calculates the symmetrically equivalent hkls (or uvws)
         for the reciprocal (or direct) point group symmetry.
         '''
-        if space == 'd':
-            mat = self.dmt.astype(np.float64)
-            if applyLaue:
-                sym = self.SYM_PG_d_laue.astype(np.float64)
-            else:
-                sym = self.SYM_PG_d.astype(np.float64)
-        elif space == 'r':
-            mat = self.rmt.astype(np.float64)
-            if applyLaue:
-                sym = self.SYM_PG_r_laue.astype(np.float64)
-            else:
-                sym = self.SYM_PG_r.astype(np.float64)
-        elif space == 'c':
+        if space not in ('d', 'r', 'c'):
+            raise ValueError('CalcStar: unrecognized space.')
+
+        if space == 'c':
             mat = np.eye(3)
-            if applyLaue:
-                sym = self.SYM_PG_c_laue.astype(np.float64)
-            else:
-                sym = self.SYM_PG_c.astype(np.float64)
         else:
-            raise ValueError('CalcStar: unrecognized space.')
+            mat = (self.dmt if space == 'd' else self.rmt).astype(float)
+
+        suffix = '_laue' if applyLaue else ''
+        name = f'SYM_PG_{space}{suffix}'
+        sym = getattr(self, name).astype(float)
 
-        vv = np.array(v).astype(np.float64)
+        vv = np.asarray(v).astype(float)
         return _calcstar(vv, sym, mat)
 
     def CalcPositions(self):