From b0d24a05cd7487090f3bbb186a6f6cd5d6ba73e8 Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Wed, 17 Jun 2026 17:28:42 +0200 Subject: [PATCH 1/5] Add `IndexedRayTransferIntegrator` and `IndexedRayTransferEmitter` classes --- cherab/tools/raytransfer/emitters.pxd | 19 +++ cherab/tools/raytransfer/emitters.pyx | 194 ++++++++++++++++++++++++++ 2 files changed, 213 insertions(+) diff --git a/cherab/tools/raytransfer/emitters.pxd b/cherab/tools/raytransfer/emitters.pxd index 49bb527a..724203b0 100644 --- a/cherab/tools/raytransfer/emitters.pxd +++ b/cherab/tools/raytransfer/emitters.pxd @@ -19,6 +19,7 @@ # under the Licence from numpy cimport ndarray +from raysect.core.math.function.float.function3d cimport Function3D from raysect.optical cimport World, Primitive, Ray, Spectrum, Point3D, Vector3D, AffineMatrix3D from raysect.optical.material cimport VolumeIntegrator, InhomogeneousVolumeEmitter @@ -44,6 +45,13 @@ cdef class CartesianRayTransferIntegrator(RayTransferIntegrator): AffineMatrix3D world_to_primitive, AffineMatrix3D primitive_to_world) +cdef class IndexedRayTransferIntegrator(RayTransferIntegrator): + + cpdef Spectrum integrate(self, Spectrum spectrum, World world, Ray ray, Primitive primitive, + InhomogeneousVolumeEmitter material, Point3D start_point, Point3D end_point, + AffineMatrix3D world_to_primitive, AffineMatrix3D primitive_to_world) + + cdef class RayTransferEmitter(InhomogeneousVolumeEmitter): cdef: @@ -72,6 +80,17 @@ cdef class CartesianRayTransferEmitter(RayTransferEmitter): cdef: double _dx, _dy, _dz + cpdef Spectrum emission_function(self, Point3D point, Vector3D direction, Spectrum spectrum, + World world, Ray ray, Primitive primitive, + AffineMatrix3D world_to_primitive, AffineMatrix3D primitive_to_world) + + +cdef class IndexedRayTransferEmitter(InhomogeneousVolumeEmitter): + + cdef: + Function3D index_function + int _bins + cpdef Spectrum emission_function(self, Point3D point, Vector3D direction, Spectrum spectrum, World world, Ray ray, Primitive primitive, AffineMatrix3D world_to_primitive, AffineMatrix3D primitive_to_world) \ No newline at end of file diff --git a/cherab/tools/raytransfer/emitters.pyx b/cherab/tools/raytransfer/emitters.pyx index c0883baa..d94201d5 100644 --- a/cherab/tools/raytransfer/emitters.pyx +++ b/cherab/tools/raytransfer/emitters.pyx @@ -25,6 +25,7 @@ with other integrators is not guaranteed. """ import numpy as np +from raysect.core.math.function.float.function3d cimport Function3D, autowrap_function3d from raysect.optical cimport World, Primitive, Ray, Spectrum, Point3D, Vector3D, AffineMatrix3D from raysect.optical.material cimport VolumeIntegrator, InhomogeneousVolumeEmitter from libc.math cimport sqrt, atan2, M_PI as pi @@ -224,6 +225,89 @@ cdef class CartesianRayTransferIntegrator(RayTransferIntegrator): return spectrum + +cdef class IndexedRayTransferIntegrator(RayTransferIntegrator): + """Calculates the distances traveled by the ray through indexed regions in 3D coordinate system: :math:`(X, Y, Z)`. + + This integrator is used with the `IndexedRayTransferEmitter` material class and + an index function to calculate ray transfer matrices (geometry matrices). + The value for each voxel is stored in respective bin of the spectral array. + The distances traveled by the ray through the voxel is calculated + approximately and the accuracy depends on the integration step. + + Parameters + ---------- + step : float + Integration step, by default 0.001. + min_samples : int + Number of minimum samples of integration, by default 2. + """ + def __init__(self, double step=0.001, int min_samples=2): + super().__init__(step=step, min_samples=min_samples) + + @cython.wraparound(False) + @cython.cdivision(True) + @cython.initializedcheck(False) + @cython.nonecheck(False) + cpdef Spectrum integrate( + self, + Spectrum spectrum, + World world, + Ray ray, + Primitive primitive, + InhomogeneousVolumeEmitter material, + Point3D start_point, + Point3D end_point, + AffineMatrix3D world_to_primitive, + AffineMatrix3D primitive_to_world + ): + + cdef: + Point3D start, end + Vector3D direction + int isource_pre, isource_current, it, n + double length, t, dt, x, y, z, res + + if not isinstance(material, IndexedRayTransferEmitter): + raise TypeError( + 'Only IndexedRayTransferEmitter material is supported ' + 'by IndexedRayTransferIntegrator' + ) + start = start_point # start point in world coordinates + end = end_point # end point in world coordinates + direction = start.vector_to(end) # direction of integration + length = direction.get_length() # integration length + if length < 0.1 * self._step: # return if ray's path is too short + return spectrum + direction = direction.normalise() # normalized direction + # number of points along ray's trajectory + n = max(self._min_samples, (length / self._step)) + dt = length / n # integration step + # cython performs checks on attributes of external class, + # so it's better to do the checks before the loop + isource_current = -1 + isource_pre = -1 + res = 0 + for it in range(n): + t = (it + 0.5) * dt + x = start.x + direction.x * t # x coordinates of the points + y = start.y + direction.y * t # y coordinates of the points + z = start.z + direction.z * t # z coordinates of the points + isource_current = material.index_function(x, y, z) # get geometry grid index + if isource_current != isource_pre: # we moved to the next cell + if isource_pre > -1: + # writing results for the current source + spectrum.samples_mv[isource_pre] += res + isource_pre = isource_current + res = 0 + if isource_current > -1: + res += dt + if isource_current > -1: + spectrum.samples_mv[isource_current] += res + + return spectrum + + cdef class RayTransferEmitter(InhomogeneousVolumeEmitter): """ Basic class for ray transfer emitters defined on a regular 3D grid. Ray transfer emitters @@ -569,3 +653,113 @@ cdef class CartesianRayTransferEmitter(RayTransferEmitter): return spectrum spectrum.samples_mv[isource] += 1. # unit emissivity return spectrum + + +cdef class IndexedRayTransferEmitter(InhomogeneousVolumeEmitter): + """A unit emitter defined by an index function, which can be used + to calculate ray transfer matrices (geometry matrices). + + Note that for performance reason there are no boundary checks in `emission_function()`, + or in `IndexedRayTransferIntegrator`, + so this emitter must be placed inside a bounding box. + + Parameters + ---------- + index_function : callable + Callable objects taking 3 positional arguments :math:`(X, Y, Z)`. + bins : int + Number of bins for the spectral array, by default 0. + integration_step : float + The length of line integration step, by default 0.01. + integrator : :obj:`~raysect.optical.material.emitter.inhomogeneous.VolumeIntegrator` + Volume integrator, by default `.IndexedRayTransferIntegrator(integration_step)`. + + Examples + -------- + .. code-block:: python + + from numpy import hypot + from raysect.optical import World, translate, Point3D + from raysect.primitive import Cylinder + from cherab.tools.raytransfer import IndexedRayTransferEmitter + + def index_func(x, y, z): + if hypot(x, y) <= 1: + return 0 + elif 1 < hypot(x, y) <= 2: + return 1 + else: + return -1 # must be set -1 outside meshes. + + world = World() + bins = 2 # Note that bins must be same as the number of meshes. + # Here thinking of two meshes like bins.shape = (2, ). + material = IndexedRayTransferEmitter(index_func, bins, integration_step=0.001) + eps = 1.e-6 # ray must never leave the grid when passing through the volume + radius = 2.0 - eps + height = 10.0 + cylinder = Cylinder( + radius, + height, + material=material, + parent=world, + transform=translate(0, 0, -0.5 * height) + ) + + camera.spectral_bins = material.bins + # ray transfer matrix will be calculated for 600.5 nm + camera.min_wavelength = 600. + camera.max_wavelength = 601. + """ + + cdef: + readonly Function3D index_function + readonly int _bins + + def __init__( + self, + object index_function not None, + int bins=0, + double integration_step=0.01, + VolumeIntegrator integrator=None + ): + + integrator = integrator or IndexedRayTransferIntegrator(step=integration_step) + super().__init__(integrator=integrator) + + self.index_function = autowrap_function3d(index_function) + self._bins = bins + + @property + def bins(self): + """ + Number of indexed regions, must not exceed the maximum of `index_function`. + + :rtype: int + """ + return self._bins + + @cython.wraparound(False) + @cython.cdivision(True) + @cython.initializedcheck(False) + @cython.nonecheck(False) + cpdef Spectrum emission_function( + self, Point3D point, + Vector3D direction, + Spectrum spectrum, + World world, + Ray ray, + Primitive primitive, + AffineMatrix3D world_to_primitive, + AffineMatrix3D primitive_to_world + ): + + cdef: + int isource + + # index of the light source for the current region + isource = self.index_function(point.x, point.y, point.z) + if isource < 0: # grid cell is not mapped to any light source + return spectrum + spectrum.samples_mv[isource] += 1. # unit emissivity + return spectrum From fc94bc2ce37974a88fb96505b1cebd0dce948a5f Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Wed, 17 Jun 2026 17:33:46 +0200 Subject: [PATCH 2/5] Refactor `IndexedRayTransferEmitter` by removing unused variables and cleaning imports --- cherab/tools/raytransfer/emitters.pyx | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/cherab/tools/raytransfer/emitters.pyx b/cherab/tools/raytransfer/emitters.pyx index d94201d5..5d6a6fdf 100644 --- a/cherab/tools/raytransfer/emitters.pyx +++ b/cherab/tools/raytransfer/emitters.pyx @@ -25,7 +25,7 @@ with other integrators is not guaranteed. """ import numpy as np -from raysect.core.math.function.float.function3d cimport Function3D, autowrap_function3d +from raysect.core.math.function.float.function3d cimport autowrap_function3d from raysect.optical cimport World, Primitive, Ray, Spectrum, Point3D, Vector3D, AffineMatrix3D from raysect.optical.material cimport VolumeIntegrator, InhomogeneousVolumeEmitter from libc.math cimport sqrt, atan2, M_PI as pi @@ -712,10 +712,6 @@ cdef class IndexedRayTransferEmitter(InhomogeneousVolumeEmitter): camera.max_wavelength = 601. """ - cdef: - readonly Function3D index_function - readonly int _bins - def __init__( self, object index_function not None, From d86532973fbace49ee3998e313e913f9c6d5c0cd Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Thu, 18 Jun 2026 10:00:51 +0200 Subject: [PATCH 3/5] Make `index_function` public in `IndexedRayTransferEmitter` class --- cherab/tools/raytransfer/emitters.pxd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cherab/tools/raytransfer/emitters.pxd b/cherab/tools/raytransfer/emitters.pxd index 724203b0..d7f5508e 100644 --- a/cherab/tools/raytransfer/emitters.pxd +++ b/cherab/tools/raytransfer/emitters.pxd @@ -88,7 +88,7 @@ cdef class CartesianRayTransferEmitter(RayTransferEmitter): cdef class IndexedRayTransferEmitter(InhomogeneousVolumeEmitter): cdef: - Function3D index_function + public Function3D index_function int _bins cpdef Spectrum emission_function(self, Point3D point, Vector3D direction, Spectrum spectrum, From 2d03077cc4d28a0632f98c625e3873a05ea4687e Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Thu, 18 Jun 2026 10:01:26 +0200 Subject: [PATCH 4/5] Use sphinx docstrings format --- cherab/tools/raytransfer/emitters.pyx | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/cherab/tools/raytransfer/emitters.pyx b/cherab/tools/raytransfer/emitters.pyx index 5d6a6fdf..6e82c806 100644 --- a/cherab/tools/raytransfer/emitters.pyx +++ b/cherab/tools/raytransfer/emitters.pyx @@ -235,12 +235,8 @@ cdef class IndexedRayTransferIntegrator(RayTransferIntegrator): The distances traveled by the ray through the voxel is calculated approximately and the accuracy depends on the integration step. - Parameters - ---------- - step : float - Integration step, by default 0.001. - min_samples : int - Number of minimum samples of integration, by default 2. + :param float step: Integration step, by default 0.001. + :param int min_samples: Number of minimum samples of integration, by default 2. """ def __init__(self, double step=0.001, int min_samples=2): super().__init__(step=step, min_samples=min_samples) @@ -663,19 +659,15 @@ cdef class IndexedRayTransferEmitter(InhomogeneousVolumeEmitter): or in `IndexedRayTransferIntegrator`, so this emitter must be placed inside a bounding box. - Parameters - ---------- - index_function : callable - Callable objects taking 3 positional arguments :math:`(X, Y, Z)`. - bins : int - Number of bins for the spectral array, by default 0. - integration_step : float - The length of line integration step, by default 0.01. - integrator : :obj:`~raysect.optical.material.emitter.inhomogeneous.VolumeIntegrator` + :param callable index_function: Callable objects taking 3 positional arguments :math:`(X, Y, Z)`. + :param int bins: Number of bins for the spectral array, by default 0. + :param float integration_step: The length of line integration step, by default 0.01. + :param VolumeIntegrator integrator: Volume integrator, by default `.IndexedRayTransferIntegrator(integration_step)`. Volume integrator, by default `.IndexedRayTransferIntegrator(integration_step)`. - Examples - -------- + :ivar callable index_function: Callable objects taking 3 positional arguments :math:`(X, Y, Z)`. + :ivar int bins: Number of bins for the spectral array. + .. code-block:: python from numpy import hypot From 6311a982070c9e9499b47ea0c4f18d42a48bc985 Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Thu, 18 Jun 2026 10:39:01 +0200 Subject: [PATCH 5/5] Add tests for IndexedRayTransferEmitter class and its functionality --- cherab/tools/tests/test_raytransfer.py | 139 ++++++++++++++++++++++++- 1 file changed, 135 insertions(+), 4 deletions(-) diff --git a/cherab/tools/tests/test_raytransfer.py b/cherab/tools/tests/test_raytransfer.py index 76ecd4d4..f6aa5390 100644 --- a/cherab/tools/tests/test_raytransfer.py +++ b/cherab/tools/tests/test_raytransfer.py @@ -17,12 +17,24 @@ # under the Licence. import unittest + import numpy as np -from raysect.optical import World, Ray, Point3D, Point2D, Vector3D, NumericalIntegrator, Spectrum +from raysect.core.math.function.float.function3d.interpolate import Discrete3DMesh +from raysect.optical import NumericalIntegrator, Point2D, Point3D, Ray, Spectrum, Vector3D, World from raysect.primitive import Box, Cylinder, Subtract -from cherab.tools.raytransfer import RayTransferBox, RayTransferCylinder, CartesianRayTransferEmitter, CylindricalRayTransferEmitter -from cherab.tools.raytransfer import RayTransferPipeline0D, RayTransferPipeline1D, RayTransferPipeline2D + from cherab.tools.inversions import ToroidalVoxelGrid +from cherab.tools.raytransfer import ( + CartesianRayTransferEmitter, + CylindricalRayTransferEmitter, + IndexedRayTransferEmitter, + IndexedRayTransferIntegrator, + RayTransferBox, + RayTransferCylinder, + RayTransferPipeline0D, + RayTransferPipeline1D, + RayTransferPipeline2D, +) class TestRayTransferCylinder(unittest.TestCase): @@ -89,7 +101,7 @@ def test_voxel_map_3d(self): self.assertTrue(np.all(mask_ref == rtc.mask) and np.all(np.array(inv_vmap_ref) == np.array(rtc.invert_voxel_map())) and rtc.bins == 8) def test_integration_2d(self): - """ Testing against ToroidalVoxelGrid""" + """Testing against ToroidalVoxelGrid""" world = World() rtc = RayTransferCylinder(radius_outer=4., height=2., n_radius=2, n_height=2, radius_inner=2., parent=world) rtc.step = 0.001 * rtc.step @@ -242,6 +254,125 @@ def test_evaluate_function_3d(self): self.assertTrue(np.allclose(spectrum_test, spectrum.samples, atol=0.001)) +class TestIndexedRayTransferEmitter(unittest.TestCase): + """ + Test cases for IndexedRayTransferEmitter class. + """ + + @staticmethod + def index_func(x, y, z): + if 0 <= x < 3 and 0 <= y < 3 and 0 <= z < 3: + return int(x) + 3 * int(y) + 9 * int(z) + return -1 + + @staticmethod + def _vertex_index(i, j, k): + return i + 4 * j + 16 * k + + @classmethod + def _build_discrete3dmesh_index_function(cls): + # 4x4x4 grid vertices covering a 3x3x3 set of unit cells. + vertex_coords = np.array( + [[float(i), float(j), float(k)] for k in range(4) for j in range(4) for i in range(4)], + dtype=np.float64, + ) + + tetrahedra = [] + tetrahedra_data = [] + for i, j, k in np.ndindex(3, 3, 3): + cell_value = float(i + 3 * j + 9 * k) + v000 = cls._vertex_index(i, j, k) + v100 = cls._vertex_index(i + 1, j, k) + v010 = cls._vertex_index(i, j + 1, k) + v110 = cls._vertex_index(i + 1, j + 1, k) + v001 = cls._vertex_index(i, j, k + 1) + v101 = cls._vertex_index(i + 1, j, k + 1) + v011 = cls._vertex_index(i, j + 1, k + 1) + v111 = cls._vertex_index(i + 1, j + 1, k + 1) + + # Split each cube into 6 tetrahedra using the 000 -> 111 body diagonal. + tetrahedra.extend( + [ + [v000, v100, v110, v111], + [v000, v100, v101, v111], + [v000, v001, v101, v111], + [v000, v001, v011, v111], + [v000, v010, v011, v111], + [v000, v010, v110, v111], + ] + ) + tetrahedra_data.extend([cell_value] * 6) + + return Discrete3DMesh( + vertex_coords, + np.array(tetrahedra, dtype=np.int32), + np.array(tetrahedra_data, dtype=np.float64), + limit=False, + default_value=-1.0, + ) + + def test_evaluate_function(self): + """ + Test IndexedRayTransferEmitter with NumericalIntegrator. + """ + world = World() + material = IndexedRayTransferEmitter(self.index_func, bins=27, integrator=NumericalIntegrator(0.0001)) + box = Box(lower=Point3D(0, 0, 0), upper=Point3D(2.99999, 2.99999, 2.99999), material=material, parent=world) + ray = Ray(origin=Point3D(4.0, 4.0, 4.0), direction=Vector3D(-1.0, -1.0, -1.0) / np.sqrt(3), min_wavelength=500.0, max_wavelength=501.0, bins=material.bins) + spectrum = ray.trace(world) + spectrum_test = np.zeros(material.bins) + spectrum_test[0] = spectrum_test[13] = spectrum_test[26] = np.sqrt(3.0) + self.assertTrue(np.allclose(spectrum_test, spectrum.samples, atol=0.001)) + + def test_default_integrator(self): + """ + Test IndexedRayTransferEmitter with IndexedRayTransferIntegrator. + """ + world = World() + material = IndexedRayTransferEmitter(self.index_func, bins=27) + self.assertTrue(isinstance(material.integrator, IndexedRayTransferIntegrator)) + box = Box(lower=Point3D(0, 0, 0), upper=Point3D(2.99999, 2.99999, 2.99999), material=material, parent=world) + ray = Ray(origin=Point3D(4.0, 4.0, 4.0), direction=Vector3D(-1.0, -1.0, -1.0) / np.sqrt(3), min_wavelength=500.0, max_wavelength=501.0, bins=material.bins) + spectrum = ray.trace(world) + spectrum_test = np.zeros(material.bins) + spectrum_test[0] = spectrum_test[13] = spectrum_test[26] = np.sqrt(3.0) + self.assertTrue(np.allclose(spectrum_test, spectrum.samples, atol=0.001)) + + def test_discrete3dmesh_as_index_function(self): + """ + Test IndexedRayTransferEmitter with a Discrete3DMesh index function　equivalent to index_func. + """ + mesh_index_func = self._build_discrete3dmesh_index_function() + + # Check representative points inside all voxels and outside the domain. + for i, j, k in np.ndindex(3, 3, 3): + x, y, z = i + 0.25, j + 0.25, k + 0.25 + self.assertEqual(int(mesh_index_func(x, y, z)), self.index_func(x, y, z)) + self.assertEqual(int(mesh_index_func(-0.1, 1.0, 1.0)), -1) + self.assertEqual(int(mesh_index_func(3.1, 1.0, 1.0)), -1) + + # Check equivalent ray transfer matrix row from both index functions. + ray = Ray( + origin=Point3D(4.0, 4.0, 4.0), + direction=Vector3D(-1.0, -1.0, -1.0) / np.sqrt(3), + min_wavelength=500.0, + max_wavelength=501.0, + bins=27, + ) + + world_ref = World() + material_ref = IndexedRayTransferEmitter(self.index_func, bins=27, integrator=NumericalIntegrator(0.0001)) + box_ref = Box(lower=Point3D(0, 0, 0), upper=Point3D(2.99999, 2.99999, 2.99999), material=material_ref, parent=world_ref) + spectrum_ref = ray.trace(world_ref) + + world_mesh = World() + material_mesh = IndexedRayTransferEmitter(mesh_index_func, bins=27, integrator=NumericalIntegrator(0.0001)) + box_mesh = Box(lower=Point3D(0, 0, 0), upper=Point3D(2.99999, 2.99999, 2.99999), material=material_mesh, parent=world_mesh) + spectrum_mesh = ray.trace(world_mesh) + + self.assertTrue(np.allclose(spectrum_ref.samples, spectrum_mesh.samples, atol=0.001)) + + class TestRayTransferPipeline0D(unittest.TestCase): """ Test cases for RayTransferPipeline0D class.