Better macro rendering for point-focused beams

AxelHenningsson · AxelHenningsson · commit 08b6a32ec845 · 2026-03-09T14:50:15.000+01:00
diff --git a/docs/_static/diffraction_pattern.png b/docs/_static/diffraction_pattern.png
diff --git a/docs/index.html b/docs/index.html
@@ -393,7 +393,8 @@ <h1>Introduction<a class="headerlink" href="#introduction" title="Link to this h
 </pre></div>
 </div>
 </div></blockquote>
-<img alt="https://github.com/FABLE-3DXRD/xrd_simulator/blob/main/docs/source/images/diffraction_pattern.png?raw=true" class="align-center" src="https://github.com/FABLE-3DXRD/xrd_simulator/blob/main/docs/source/images/diffraction_pattern.png?raw=true" />
+<a class="reference internal image-reference" href="https://github.com/FABLE-3DXRD/xrd_simulator/blob/main/docs/source/images/diffraction_pattern.png?raw=true"><img alt="https://github.com/FABLE-3DXRD/xrd_simulator/blob/main/docs/source/images/diffraction_pattern.png?raw=true" class="align-center" src="https://github.com/FABLE-3DXRD/xrd_simulator/blob/main/docs/source/images/diffraction_pattern.png?raw=true" style="width: 95%;" />
+</a>
 <p>To compute several frames simply change the motion and collect the diffraction again. The sample may be moved before
 each computation using the same or another motion.</p>
 <blockquote>
diff --git a/tests/test_detector.py b/tests/test_detector.py
@@ -3,8 +3,11 @@
 
 import numpy as np
 import torch
+from scipy.ndimage import label
 from scipy.spatial import ConvexHull
+from scipy.spatial.transform import Rotation
 
+from xrd_simulator import polycrystal
 from xrd_simulator.beam import Beam
 from xrd_simulator.detector import Detector
 from xrd_simulator.mesh import TetraMesh
@@ -704,6 +707,107 @@ def test_lorentz_infinity_handling_nano(self):
             msg="Pattern should be empty when all peaks have infinite Lorentz",
         )
 
+    def test_macro_scanning_3dxrd_setup(self):
+
+        radius = 100.0
+        mesh = TetraMesh.generate_mesh_from_levelset(
+            level_set=lambda x: np.linalg.norm(x) - radius,
+            bounding_radius=radius + 5.0,
+            max_cell_circumradius=50.0,
+        )
+
+        mesh.translate(-mesh.ecentroids.mean(axis=0))
+        data = os.path.join(
+            os.path.join(os.path.dirname(__file__), "data"),
+            "Fe_mp-150_conventional_standard.cif",
+        )
+        unit_cell = [3.64570000, 3.64570000, 3.64570000, 90.0, 90.0, 90.0]
+        sgname = "Fm-3m"
+        phase = Phase(unit_cell, sgname, path_to_cif_file=data)
+
+        detector_distance = 98932.41043969788
+        pixel_size = 75  # microns
+        detector_size = 2048 * pixel_size
+        energy = 55.5  # kev
+        wavelength = 12.3984198742 / energy  # angstrom
+
+        phase._setup_diffracting_planes(wavelength, 0, 20 * torch.pi / 180)
+
+        pc = polycrystal.Polycrystal(
+            mesh=mesh,
+            orientation=Rotation.random(mesh.number_of_elements).as_matrix(),
+            strain=np.zeros((3, 3)),
+            phases=[Phase(unit_cell, sgname, path_to_cif_file=None)],
+            element_phase_map=np.zeros((mesh.number_of_elements,)).astype(int),
+        )
+
+        det_corner_0 = np.array(
+            [detector_distance, -detector_size / 2.0, -detector_size / 2.0]
+        )
+        det_corner_1 = np.array(
+            [detector_distance, detector_size / 2.0, -detector_size / 2.0]
+        )
+        det_corner_2 = np.array(
+            [detector_distance, -detector_size / 2.0, detector_size / 2.0]
+        )
+
+        beam_y = 2.0  # um
+        t = beam_y / 2.0
+        dz = beam_y / 2.0
+        beam_vertices = np.array(
+            [
+                [-detector_distance, -t, -dz],
+                [-detector_distance, t, -dz],
+                [-detector_distance, t, dz],
+                [-detector_distance, -t, dz],
+                [detector_distance, -t, -dz],
+                [detector_distance, t, -dz],
+                [detector_distance, t, dz],
+                [detector_distance, -t, dz],
+            ]
+        )
+        xray_propagation_direction = np.array([1, 0, 0]) * 2 * np.pi / wavelength
+        polarization_vector = np.array([0, 1, 0])
+        xray_beam = Beam(
+            beam_vertices, xray_propagation_direction, wavelength, polarization_vector
+        )
+        rotation_angle = np.radians(20.0)
+        rotation_axis = np.array([0, 0, 1])
+        single_frame_motion = RigidBodyMotion(
+            rotation_axis,
+            rotation_angle,
+            translation=np.array([0, 0, 0]),
+        )
+
+        eiger = Detector(
+            det_corner_0,
+            det_corner_1,
+            det_corner_2,
+            pixel_size=(pixel_size, pixel_size),
+            gaussian_sigma=0.4,
+            kernel_threshold=0.001,
+            max_gaussian_kernel_radius=9,
+        )
+        peaks = pc.diffract(
+            xray_beam, single_frame_motion, detector=eiger, verbose=True
+        )
+        diffraction_pattern, peaks = eiger.render(
+            peaks,
+            frames_to_render=1,
+            method="macro",
+        )
+
+        im = diffraction_pattern[0].numpy()
+        im[np.isnan(im)] = 0
+
+        self.assertGreaterEqual(im.max(), 0, msg="Diffraction pattern is empty")
+
+        _, num_labels = label(im)
+
+        self.assertGreaterEqual(
+            num_labels, 20, msg="There appears to be very few diffraction spots"
+        )
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_utils.py b/tests/test_utils.py
@@ -1,16 +1,37 @@
 import unittest
 import warnings
+
 import numpy as np
+from scipy.spatial import ConvexHull
+from scipy.spatial.transform import Rotation
 from xfab import tools
+
 from xrd_simulator import utils
-from scipy.spatial.transform import Rotation
 
+rng = np.random.default_rng(0)
 
-class TestUtils(unittest.TestCase):
 
+class TestUtils(unittest.TestCase):
     def setUp(self):
         np.random.seed(10)  # changes all randomisation in the test
 
+    def test_sample_convex_hull_3d(self):
+        hull = ConvexHull(rng.random((100, 3)))
+        points = hull.points
+        r = points.mean(axis=0)
+        tris = points[hull.simplices]
+        a = tris[:, 0] - r
+        b = tris[:, 1] - r
+        c = tris[:, 2] - r
+        vols = np.abs(np.einsum("ij,ij->i", a, np.cross(b, c))) / 6.0
+        cents = (r + tris[:, 0] + tris[:, 1] + tris[:, 2]) / 4.0
+        cents = np.sum(cents * vols[:, None], axis=0) / vols.sum()
+        X = utils._sample_convex_hull_3d(hull, 200000)
+        A = hull.equations[:, :-1]
+        b = hull.equations[:, -1]
+        assert np.all(X @ A.T + b <= 1e-10)
+        assert np.allclose(X.mean(axis=0), cents, atol=1e-2)
+
     def test_clip_line_with_convex_polyhedron(self):
         line_points = np.ascontiguousarray([[-1.0, 0.2, 0.2], [-1.0, 0.4, 0.6]])
         line_direction = np.ascontiguousarray([1.0, 0.0, 0.0])
@@ -149,7 +170,7 @@ def test_get_misorientations(self):
 
     def test_diffractogram_deprecated(self):
         """Test that _diffractogram raises a deprecation warning.
-        
+
         .. deprecated::
             This test verifies that _diffractogram is properly marked as deprecated.
             The function will be removed in a future version.
@@ -171,7 +192,7 @@ def test_diffractogram_deprecated(self):
             # Check that deprecation warning was raised
             self.assertTrue(
                 any(issubclass(warning.category, DeprecationWarning) for warning in w),
-                msg="_diffractogram should raise DeprecationWarning"
+                msg="_diffractogram should raise DeprecationWarning",
             )
 
         # Verify function still works correctly (for backward compatibility)
@@ -186,23 +207,23 @@ def test_diffractogram_deprecated(self):
 
     def test_contained_by_intervals_deprecated(self):
         """Test that _contained_by_intervals raises a deprecation warning.
-        
+
         .. deprecated::
             This test verifies that _contained_by_intervals is properly marked as deprecated.
             The function will be removed in a future version.
         """
         intervals = [[0.0, 0.5], [0.7, 1.0]]
-        
+
         # Verify deprecation warning is raised
         with warnings.catch_warnings(record=True) as w:
             warnings.simplefilter("always")
             result = utils._contained_by_intervals(0.3, intervals)
             # Check that deprecation warning was raised
             self.assertTrue(
                 any(issubclass(warning.category, DeprecationWarning) for warning in w),
-                msg="_contained_by_intervals should raise DeprecationWarning"
+                msg="_contained_by_intervals should raise DeprecationWarning",
             )
-        
+
         # Verify function still works correctly (for backward compatibility)
         self.assertTrue(result, msg="0.3 should be contained in [0.0, 0.5]")
 
diff --git a/xrd_simulator/detector.py b/xrd_simulator/detector.py
@@ -1591,6 +1591,13 @@ def _project_convex_hull(
     ) -> npt.NDArray:
         """Project convex hull onto detector region.
 
+        Swithces between a sample driven and detector driven approcimation depending on
+        the detector pixel size and the convex hull size. When the convex hull is large
+        compared to the pixel size, the detector driven approach is used. This means that
+        rays traced from pixel centroids are used as projection lines. When the convex
+        hull is small compared to the pixel size, the sample driven approach is used. This
+        means that points are sampled from the interior of the hull and used as projection lines.
+
         Parameters
         ----------
         proj_context : dict
@@ -1605,12 +1612,10 @@ def _project_convex_hull(
         Returns
         -------
         torch.Tensor
-            Array of clip lengths (path lengths) through the hull for each
-            detector ray.
+            Array of projected fractional contribution to each detector pixel
+            in the box. The sum total contribution is proportional to the
+            volume of the convex hull.
         """
-        # Get pixel coordinates - keep in torch
-        pixel_coords = self.pixel_coordinates[box[0] : box[1], box[2] : box[3], :]
-        ray_points = pixel_coords.reshape((box[1] - box[0]) * (box[3] - box[2]), 3)
 
         # ConvexHull from scipy returns numpy, convert to torch
         hull = proj_context["convex_hull"]
@@ -1620,14 +1625,62 @@ def _project_convex_hull(
         )
         plane_points = -plane_offsets * plane_normals
 
+        n_pixels_in_box = (box[1] - box[0]) * (box[3] - box[2])
+
         # Ray direction in torch
         scattered_vec = proj_context["scattered_wave_vector"]
         ray_direction = scattered_vec / torch.linalg.norm(scattered_vec)
 
-        clip_lengths = utils._clip_line_with_convex_polyhedron(
-            ray_points, ray_direction, plane_points, plane_normals
-        )
-        clip_lengths = clip_lengths.reshape(box[1] - box[0], box[3] - box[2])
+        if n_pixels_in_box == 1:
+            # tetra is fully contained in a detector pixel, use the volume of the tetra as intensity
+            return torch.ones(1, 1) * proj_context["convex_hull"].volume
+        elif (
+            n_pixels_in_box <= 4
+        ):  # tetra is small compared to the pixel size, adapt to sampling approach.
+            # This is a decent approximation when the tetras are small in relation to the pixel size
+            n_samples = (
+                7 * n_pixels_in_box
+            )  # TODO the ad-hoc number 7 here should perhaps be a parameter to give controll over the accuracy of the statistics.
+            ray_points = utils._sample_convex_hull_3d(
+                proj_context["convex_hull"], n_samples
+            )
+            ray_points.reshape(n_samples, 3)
+            ray_points = ensure_torch(ray_points)
+
+            clip_lengths = utils._clip_line_with_convex_polyhedron(
+                ray_points, ray_direction, plane_points, plane_normals
+            )
+            # now map the clip lengths to the detector pixels
+            ray_directions = torch.tile(ray_direction, (n_samples, 1))
+            out = self._get_intersection(ray_directions, ray_points)
+            zd = out[:, 0]
+            yd = out[:, 1]
+
+            row_index = ensure_numpy(zd / self.pixel_size_z).astype(int) - box[0]
+            col_index = ensure_numpy(yd / self.pixel_size_y).astype(int) - box[2]
+
+            box_intensities = np.zeros((box[1] - box[0], box[3] - box[2]))
+            np.add.at(box_intensities, (row_index, col_index), clip_lengths)
+            box_intensities /= box_intensities.sum()
+            box_intensities *= proj_context["convex_hull"].volume
+
+            return ensure_torch(box_intensities)
+
+        else:
+            # This is a good approximation when the tetras are large in relation to the pixel size
+            # Get pixel coordinates and project on pixel centorid ray-lines
+            pixel_coords = self.pixel_coordinates[box[0] : box[1], box[2] : box[3], :]
+            ray_points = pixel_coords.reshape((box[1] - box[0]) * (box[3] - box[2]), 3)
+
+            clip_lengths = utils._clip_line_with_convex_polyhedron(
+                ray_points, ray_direction, plane_points, plane_normals
+            )
+
+            # ensure that scattering is proportional to the material volume
+            # that is scattering.
+            clip_lengths /= clip_lengths.sum()
+            clip_lengths *= proj_context["convex_hull"].volume
+            clip_lengths = clip_lengths.reshape(box[1] - box[0], box[3] - box[2])
 
         return clip_lengths
 
diff --git a/xrd_simulator/utils.py b/xrd_simulator/utils.py
@@ -83,6 +83,8 @@
 
 from xrd_simulator.cuda import get_selected_device
 
+rng = np.random.default_rng()
+
 
 def _cif_open(cif_file):
     """Helper function to be able to use the ``.CIFread`` of ``xfab``."""
@@ -122,6 +124,47 @@ def _print_progress(progress_fraction, message):
         print("", flush=True)
 
 
+def _sample_convex_hull_3d(hull, n_samples):
+    """Uniformly randomly sample points from the interior of a convex hull.
+
+    Parameters
+    ----------
+    hull : scipy.spatial.ConvexHull
+        Convex hull to sample from.
+    n_samples : int
+        Number of samples to draw.
+
+    Returns
+    -------
+    numpy.ndarray
+        Sampled points, shape ``(n_samples, 3)``.
+    """
+
+    # Split the convex hull into tetrahedra
+    points = hull.points
+    r = points.mean(axis=0)
+    tris = points[hull.simplices]
+    tets = np.empty((len(tris), 4, 3), dtype=points.dtype)
+    tets[:, 0] = r
+    tets[:, 1:] = tris
+    a = tets[:, 1] - tets[:, 0]
+    b = tets[:, 2] - tets[:, 0]
+    c = tets[:, 3] - tets[:, 0]
+
+    # Compute the volumes of the tetrahedra and randomly select tets
+    # based on volume weights.
+    vols = np.abs(np.einsum("ij,ij->i", a, np.cross(b, c))) / 6.0
+    idx = rng.choice(len(tets), size=n_samples, p=vols / vols.sum())
+
+    # Now use barycentric coordinates to uniformly sample points from
+    # the selected tetrahedra
+    w = rng.exponential(size=(n_samples, 4))
+    w /= w.sum(axis=1, keepdims=True)
+    random_points = np.einsum("ni,nij->nj", w, tets[idx])
+
+    return random_points
+
+
 def _clip_line_with_convex_polyhedron(
     line_points, line_direction, plane_points, plane_normals
 ):
