# The Clear BSD License
#
# Copyright (c) 2026 Tobias Heibges
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted (subject to the limitations in the disclaimer
# below) provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""
Planar Detector Implementation
Provides a finite rectangular planar detector for imaging applications.
Examples
--------
>>> from lsurf.detectors.small import PlanarDetector
>>>
>>> detector = PlanarDetector(
... center=(0, 0, 100),
... normal=(0, 0, -1), # Facing back toward source
... width=0.1, # 10 cm wide
... height=0.1, # 10 cm tall
... name="Imaging detector"
... )
>>> result = detector.detect(rays)
"""
import numpy as np
from numpy.typing import NDArray
from ...utilities.ray_data import RayBatch
from ..base import DetectionEvent
from ..results import DetectorResult
[docs]
class PlanarDetector:
"""
Planar detector with finite rectangular size.
Detects rays that intersect a rectangular plane at a specific position
and orientation. Useful for imaging applications and beam profiling.
Parameters
----------
center : tuple of float
Center position of detector plane (x, y, z) in meters.
normal : tuple of float
Normal vector (defines which direction detector faces).
width : float
Width of detector (in local u direction) in meters.
height : float
Height of detector (in local v direction) in meters.
name : str, optional
Detector name. Default is "Planar Detector".
Attributes
----------
center : ndarray, shape (3,)
Detector center position.
normal : ndarray, shape (3,)
Unit normal vector.
width : float
Detector width.
height : float
Detector height.
u : ndarray, shape (3,)
Local x-axis direction (width direction).
v : ndarray, shape (3,)
Local y-axis direction (height direction).
name : str
Detector name.
accumulated_result : DetectorResult
All accumulated detections since last clear().
Notes
-----
The local coordinate system is constructed with:
- `normal`: direction the detector faces
- `u`: perpendicular to normal (width direction)
- `v`: perpendicular to both normal and u (height direction)
Only rays traveling toward the front face of the detector (opposite
to the normal direction) are detected.
Examples
--------
>>> detector = PlanarDetector(
... center=(0, 0, 50),
... normal=(0, 0, -1), # Facing -z direction
... width=0.1,
... height=0.1
... )
>>> result = detector.detect(rays)
"""
[docs]
def __init__(
self,
center: tuple[float, float, float],
normal: tuple[float, float, float],
width: float,
height: float,
name: str = "Planar Detector",
):
"""
Initialize planar detector.
Parameters
----------
center : tuple of float
Center position of detector plane.
normal : tuple of float
Normal vector (detector faces in this direction).
width : float
Detector width in meters.
height : float
Detector height in meters.
name : str, optional
Detector name.
"""
self.name = name
self.center = np.array(center, dtype=np.float32)
self.normal = np.array(normal, dtype=np.float32)
self.normal = self.normal / np.linalg.norm(self.normal)
self.width = width
self.height = height
# Create local coordinate system
if abs(self.normal[2]) < 0.9:
self.u = np.cross(self.normal, np.array([0, 0, 1], dtype=np.float32))
else:
self.u = np.cross(self.normal, np.array([1, 0, 0], dtype=np.float32))
self.u = self.u / np.linalg.norm(self.u)
self.v = np.cross(self.normal, self.u)
self.v = self.v / np.linalg.norm(self.v)
self._accumulated_result = DetectorResult.empty(name)
self._events: list[DetectionEvent] = []
@property
def accumulated_result(self) -> DetectorResult:
"""All accumulated detections since last clear()."""
return self._accumulated_result
@property
def events(self) -> list[DetectionEvent]:
"""
Backward compatibility: list of DetectionEvent objects.
For new code, use accumulated_result instead.
"""
return self._events
[docs]
def detect(
self, rays: RayBatch, current_time: float = 0.0, accumulate: bool = True
) -> DetectorResult:
"""
Detect rays that intersect the detector plane within bounds.
Parameters
----------
rays : RayBatch
Ray batch to test.
current_time : float, optional
Current simulation time. Default is 0.0.
accumulate : bool, optional
Whether to accumulate results. Default is True.
Returns
-------
DetectorResult
Newly detected rays.
Notes
-----
Only detects rays that:
- Are traveling toward the front face (opposite to normal)
- Intersect the plane in the forward direction
- Hit within the width x height bounds
"""
if rays is None or rays.num_rays == 0:
return DetectorResult.empty(self.name)
c = 299792458.0 # Speed of light in m/s
n = 1.0 # Refractive index of air
# Vectorized computation
active_mask = rays.active
if not np.any(active_mask):
return DetectorResult.empty(self.name)
origins = rays.positions[active_mask]
directions = rays.directions[active_mask]
active_indices = np.where(active_mask)[0]
# Check if rays intersect plane
denoms = np.dot(directions, self.normal)
# Parallel to plane or pointing away?
valid_mask = denoms < -1e-10
if not np.any(valid_mask):
return DetectorResult.empty(self.name)
# Find intersection parameter
to_center = self.center - origins[valid_mask]
t = np.dot(to_center, self.normal) / denoms[valid_mask]
# Only forward propagation
forward_mask = t > 0
if not np.any(forward_mask):
return DetectorResult.empty(self.name)
# Update indices
valid_active_indices = active_indices[valid_mask][forward_mask]
valid_origins = origins[valid_mask][forward_mask]
valid_directions = directions[valid_mask][forward_mask]
valid_t = t[forward_mask]
# Intersection points
intersections = valid_origins + valid_t[:, np.newaxis] * valid_directions
# Check if within detector bounds
offsets = intersections - self.center
local_x = np.dot(offsets, self.u)
local_y = np.dot(offsets, self.v)
in_bounds = (np.abs(local_x) <= self.width / 2) & (
np.abs(local_y) <= self.height / 2
)
if not np.any(in_bounds):
return DetectorResult.empty(self.name)
# Get final results
final_indices = valid_active_indices[in_bounds]
final_intersections = intersections[in_bounds]
final_t = valid_t[in_bounds]
# Compute arrival times
additional_times = final_t * n / c
arrival_times = rays.accumulated_time[final_indices] + additional_times
result = DetectorResult(
positions=final_intersections.astype(np.float32),
directions=rays.directions[final_indices].astype(np.float32),
times=arrival_times.astype(np.float32),
intensities=rays.intensities[final_indices].astype(np.float32),
wavelengths=rays.wavelengths[final_indices].astype(np.float32),
ray_indices=final_indices.astype(np.int32),
generations=(
rays.generations[final_indices].astype(np.int32)
if rays.generations is not None
else None
),
polarization_vectors=(
rays.polarization_vector[final_indices].astype(np.float32)
if rays.polarization_vector is not None
else None
),
detector_name=self.name,
)
if accumulate:
self._accumulated_result = DetectorResult.merge(
[self._accumulated_result, result]
)
self._events.extend(result.to_detection_events())
return result
[docs]
def clear(self) -> None:
"""
Clear all recorded detections.
Resets the detector to its initial state with no recorded events.
"""
self._accumulated_result = DetectorResult.empty(self.name)
self._events = []
[docs]
def get_image(
self, bins_u: int = 100, bins_v: int = 100
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
"""
Generate intensity image from detection events.
Parameters
----------
bins_u : int, optional
Number of bins in u direction. Default is 100.
bins_v : int, optional
Number of bins in v direction. Default is 100.
Returns
-------
u_centers : ndarray, shape (bins_u,)
Bin centers in u direction (meters).
v_centers : ndarray, shape (bins_v,)
Bin centers in v direction (meters).
image : ndarray, shape (bins_v, bins_u)
Intensity image (sum of intensities per bin).
"""
if self._accumulated_result.is_empty:
u_centers = np.linspace(-self.width / 2, self.width / 2, bins_u)
v_centers = np.linspace(-self.height / 2, self.height / 2, bins_v)
return u_centers, v_centers, np.zeros((bins_v, bins_u))
# Extract local coordinates
positions = self._accumulated_result.positions
offsets = positions - self.center
local_u = np.dot(offsets, self.u)
local_v = np.dot(offsets, self.v)
intensities = self._accumulated_result.intensities
# Create 2D histogram weighted by intensity
image, u_edges, v_edges = np.histogram2d(
local_u,
local_v,
bins=[bins_u, bins_v],
range=[
[-self.width / 2, self.width / 2],
[-self.height / 2, self.height / 2],
],
weights=intensities,
)
u_centers = (u_edges[:-1] + u_edges[1:]) / 2
v_centers = (v_edges[:-1] + v_edges[1:]) / 2
return u_centers, v_centers, image.T
[docs]
def __repr__(self) -> str:
"""Return string representation."""
return (
f"PlanarDetector(center={self.center.tolist()}, "
f"size=({self.width}, {self.height}), rays={self._accumulated_result.num_rays})"
)
[docs]
def __len__(self) -> int:
"""Return number of detected rays."""
return self._accumulated_result.num_rays
# Backward compatibility methods from old Detector base class
[docs]
def get_arrival_times(self) -> NDArray[np.float64]:
"""Get array of all arrival times."""
return self._accumulated_result.times.astype(np.float64)
[docs]
def get_arrival_angles(
self, reference_direction: NDArray[np.float32]
) -> NDArray[np.float64]:
"""Get angles between ray directions and reference direction."""
if self._accumulated_result.is_empty:
return np.array([], dtype=np.float64)
ref = reference_direction / np.linalg.norm(reference_direction)
dir_norms = self._accumulated_result.directions / np.linalg.norm(
self._accumulated_result.directions, axis=1, keepdims=True
)
cos_angles = np.dot(dir_norms, ref)
cos_angles = np.clip(cos_angles, -1.0, 1.0)
return np.arccos(cos_angles).astype(np.float64)
[docs]
def get_intensities(self) -> NDArray[np.float64]:
"""Get array of all detected intensities."""
return self._accumulated_result.intensities.astype(np.float64)
[docs]
def get_wavelengths(self) -> NDArray[np.float64]:
"""Get array of all detected wavelengths."""
return self._accumulated_result.wavelengths.astype(np.float64)
[docs]
def get_positions(self) -> NDArray[np.float32]:
"""Get array of all detection positions."""
return self._accumulated_result.positions.astype(np.float32)
[docs]
def get_total_intensity(self) -> float:
"""Get sum of all detected intensities."""
return self._accumulated_result.total_intensity