Skip to content

Commit a25c3c8

Browse files
committed
Merge branch 'develop' into 'main'
Version 0.3.1 See merge request e040/e0404/pyRadPlan!120
2 parents 8899ea2 + f06528d commit a25c3c8

File tree

15 files changed

+820
-323
lines changed

15 files changed

+820
-323
lines changed

.readthedocs.yaml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ version: "2"
33
build:
44
os: "ubuntu-lts-latest"
55
tools:
6-
python: "latest"
6+
python: "3.13"
77

88
python:
99
install:

CHANGELOG.md

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,12 @@
11
# Changelog
22

3+
## Version 0.3.1
4+
Small update to drop `numba` as a mandatory dependency. `numba` accelerations can still be provided but should only compliment Array API conform base implementations and be appropriately checked.
5+
- Performant and Array API compatible candidate ray matrix setup alternative for cube raytracing
6+
- Fix readthedocs-YAML to fix Python version
7+
- Add convenience plotting function to display multiple slices `plot_distributions`
8+
- Small code quality fixes
9+
310
## Version 0.3.0 - Pre-Release
411
This update features Array API compatibility, a FRED interface, VHEE model, documentation, examples in jupytext style, some refactoring and bugfixes.
512
We further are removing python 3.9 support due to array api compatibility.

README.md

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,8 @@ Development is lead by the [Radiotherapy Optimization group](https://www.dkfz.de
1616
## Concept and Goals
1717
pyRadPlan is a multi-modality treatment planning toolkit in python born from the established Matlab-based toolkit [matRad](http://www.matRad.org). As such, pyRadPlan aims to provide a framework as well as tools for combining dose calculation with optimization with focus on ion planning.
1818

19+
One of pyRadPlan's main goals is for its algorithms to be [Python Array API](https://data-apis.org/array-api/) conform as much as possible, to allow free choice of compute backends and devices (i.e., [`numpy`](https://numpy.org/), [`cupy`](https://cupy.dev/), [`pytorch`](https://pytorch.org/)). Further, its data structures (more precisely, their metadata) should be easily serializable for exchange with novel LLMs. This facilitate performant treatment planning research while enabling the integration of AI - such as deep neural networks implemented in [`pytorch`](https://pytorch.org/) or interaction with an LLM Agent via [`pydantic-ai`](https://ai.pydantic.dev/) - at every step of the treatment planning workflow.
20+
1921
### Data Structures
2022
pyRadPlan uses a similar datastructure and workflow concept as in matRad, while trying to ensure that the corresponding datastructures can be easily imported and exported from/to matRad. This facilitates the application of either algortihms from matRad or native pyRadPlan at any stage of the treatment planning workflow.
2123

examples/utils_plotting.py

Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
# %% [markdown]
2+
"""# Multiple examples on how to use the visualization tools provided by pyRadPlan."""
3+
# %% [markdown]
4+
# This example demonstrates the usage of `plot_slice()`, `plot_distributions()` and `plot_3d()`.
5+
6+
# To display this script in a Jupyter Notebook, you need to install jupytext via pip and run the following command.
7+
# This will create a .ipynb file in the same directory:
8+
9+
# ```bash
10+
# pip install jupytext
11+
# jupytext --to notebook path/to/this/file/utils_plotting.py
12+
13+
# %%
14+
# Import necessary libraries
15+
import logging
16+
17+
import matplotlib.pyplot as plt
18+
19+
from pyRadPlan import (
20+
IonPlan,
21+
generate_stf,
22+
calc_dose_influence,
23+
fluence_optimization,
24+
plot_slice,
25+
plot_multiple_slices,
26+
load_tg119,
27+
)
28+
29+
logging.basicConfig(level=logging.INFO)
30+
31+
# %% [markdown]
32+
# ### Just like in other examples, we need to generate some distributions/quantities to visualize.
33+
# %%
34+
# load TG119
35+
ct, cst = load_tg119()
36+
37+
# Create a plan object
38+
pln = IonPlan(radiation_mode="carbon", machine="Generic")
39+
pln.prop_opt = {"solver": "scipy"}
40+
# Lets calc the biological dose too
41+
pln.prop_dose_calc = {"calc_bio_dose": True}
42+
43+
# Generate Steering Geometry ("stf")
44+
stf = generate_stf(ct, cst, pln)
45+
46+
# Calculate Dose Influence Matrix ("dij")
47+
dij = calc_dose_influence(ct, cst, stf, pln)
48+
49+
# Run fluence optimization and compute the result
50+
fluence = fluence_optimization(ct, cst, stf, dij, pln)
51+
result = dij.compute_result_ct_grid(fluence)
52+
53+
# %% [markdown]
54+
# ### Visualizing a single slice with `plot_slice()`
55+
# %%
56+
57+
# Visualize only ct and choosen quantity
58+
plot_slice(ct=ct, overlay=result["physical_dose"])
59+
60+
# Visualize ct, cst and choosen quantity
61+
plot_slice(ct=ct, cst=cst, overlay=result["physical_dose"])
62+
63+
# %% [markdown]
64+
# ## Visualze more abstract settings using `plot_slice()` <br>
65+
66+
# `plot_slice()` has multiple tweakable parameters:
67+
# - **ct**: The CT data to visualize
68+
# - **cst**: The structure set to visualize (optional)
69+
# - **overlay**: The quantity to visualize
70+
# - **view_slice**: Which slice to visualize (default: middle slice)
71+
# - **plane**: Which plane to visualize (default: "axial")
72+
# - **overlay_alpha**: Transparency of the overlay (default: 0.5)
73+
# - **overlay_rel_threshold**: Relative threshold for the overlay (default: 0.01)
74+
# - **overlay_unit**: The unit of the overlay quantity (default: "Gy")
75+
# - **save_filename**: If provided, saves the plot to a file
76+
# - **show_plot**: Whether to show the plot (default: True)
77+
# - **use_global_max**: If True, uses the global maximum of the overlay for scaling (default: False)
78+
# - **ax**: If provided, plots on the given axes (useful for subplots) - use of 'plot_distributions()' is recommended
79+
# %%
80+
# Feel free to change the parameters to see how they affect the plot.
81+
plot_slice(
82+
ct=ct,
83+
cst=cst,
84+
overlay=result["physical_dose"],
85+
view_slice=64, # Visualize slice 10
86+
plane="coronal", # axial, coronal or sagittal
87+
overlay_alpha=0.5, # Transparency of the overlay
88+
overlay_unit="Gy", # Gy, dimensionless, etc.
89+
overlay_rel_threshold=0.01, # Relative threshold for the overlay
90+
contour_line_width=1.0, # Width of the contour lines
91+
save_filename=None, # alt: path/to/save/plot.png
92+
show_plot=True, # Show the plot
93+
use_global_max=False, # Do not use global max for scaling
94+
)
95+
96+
# %% [markdown]
97+
# ### Visualizing multiple overlays/quantites with `plot_distributions()`
98+
# This function might be useful in cases you want to compare multiple overlays/quantities/beams side by side.
99+
100+
# `plot_distributions()` has similar parameters to `plot_slice()`:
101+
# - **ct**: The CT data to visualize
102+
# - **cst**: The structure set to visualize (optional)
103+
# - **overlays**: List of overlay images to visualize
104+
# - **view_slice**: Which slices to visualize (default: middle slice)
105+
# - **plane**: Which plane to visualize (default: "axial")
106+
# - **overlay_alpha**: Transparency of the overlay (default: 0.5)
107+
# - **overlay_unit**: List of units for each overlay (default: "Gy")
108+
# - **overlay_rel_threshold**: Relative threshold for the overlay (default: 0.01)
109+
# - **contour_line_width**: Line width for the contour lines (default: 1.0)
110+
# - **save_filename**: If provided, saves the plot to a file
111+
# - **show_plot**: Whether to show the plot (default: True)
112+
# - **use_global_max**: If True, uses the global maximum of the overlay for scaling
113+
# - **overlay_titles**: Custom titles for each overlay type (optional)
114+
# %%
115+
# Feel free to change the parameters to see how they affect the plot.
116+
plot_multiple_slices(
117+
ct=ct,
118+
cst=cst,
119+
overlays=[result["physical_dose"], result["effect"]],
120+
view_slice=[64, 65], # Visualize slices 64 and 65
121+
plane="axial", # axial, coronal or sagittal
122+
overlay_alpha=0.5, # Transparency of the overlay
123+
overlay_unit=["Gy", "dimensionless"], # Units for each overlay
124+
overlay_rel_threshold=0.01, # Relative threshold for the overlay
125+
contour_line_width=1.0, # Width of the contour lines
126+
save_filename=None, # alt: path/to/save/plot.png
127+
show_plot=True, # Show the plot
128+
use_global_max=False, # Do not use global max for scaling
129+
overlay_titles=["Physical Dose", "Biological Effect"], # Titles for each overlay
130+
)
131+
132+
# %% [markdown]
133+
# ### Being more flexible with only `plot_slice()` is also possible:
134+
# %%
135+
# Create a figure with subplots for side-by-side comparison
136+
# 2 rows (physical dose, biological effect) x 2 cols (slice1, slice2)
137+
slices = [64, 65] # Slices to visualize
138+
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
139+
140+
# Plot physical dose for both slices
141+
for i, slice_idx in enumerate(slices):
142+
plot_slice(
143+
ct=ct,
144+
cst=cst,
145+
overlay=result["physical_dose"],
146+
view_slice=slice_idx, # Single slice
147+
plane="axial",
148+
overlay_unit="Gy",
149+
show_plot=False,
150+
ax=axes[0, i], # Top row
151+
)
152+
axes[0, i].set_title(f"Physical Dose - Slice {slice_idx}")
153+
154+
# Plot biological effect for both slices
155+
for i, slice_idx in enumerate(slices):
156+
plot_slice(
157+
ct=ct,
158+
cst=cst,
159+
overlay=result["effect"],
160+
view_slice=slice_idx, # Single slice
161+
plane="axial",
162+
overlay_unit="dimensionless",
163+
show_plot=False,
164+
ax=axes[1, i], # Bottom row
165+
)
166+
axes[1, i].set_title(f"Biological Effect - Slice {slice_idx}")
167+
168+
plt.tight_layout()
169+
plt.show()

pyRadPlan/__init__.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
from .dose._calc_dose import calc_dose_influence, calc_dose_forward
3636
from .optimization._fluence_optimization import fluence_optimization
3737
from .analysis._dvh import DVH, DVHCollection
38-
from .visualization import plot_slice
38+
from .visualization import plot_slice, plot_multiple_slices
3939
from .io import load_patient, load_tg119
4040
from .core import xp_utils
4141

@@ -68,6 +68,7 @@
6868
"SteeringInformation",
6969
"validate_stf",
7070
"plot_slice",
71+
"plot_multiple_slices",
7172
"load_patient",
7273
"load_tg119",
7374
]

pyRadPlan/dose/engines/_factory.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,8 @@ def get_available_engines(pln: Union[Plan, dict[str]]) -> dict[str, Type[DoseEng
4848
4949
Returns
5050
-------
51-
list
52-
A list of available engines.
51+
dict[str, Type[DoseEngineBase]]
52+
A dictionary containing the available engines.
5353
"""
5454
pln = validate_pln(pln)
5555
return {

pyRadPlan/raytracer/_base.py

Lines changed: 68 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,17 @@
77

88
import numpy as np
99
import SimpleITK as sitk
10+
import array_api_compat
1011

11-
from pyRadPlan.core.np2sitk import linear_indices_to_image_coordinates
12-
from pyRadPlan.geometry import lps
13-
from pyRadPlan.stf._beam import Beam
12+
from ..core.xp_utils.typing import Array
13+
from ..core.np2sitk import linear_indices_to_image_coordinates
14+
from ..geometry import lps
15+
from ..stf._beam import Beam
1416

15-
from ._perf import fast_spatial_circle_lookup
17+
try:
18+
from ._numba_perf import fast_spatial_circle_lookup
19+
except ImportError:
20+
fast_spatial_circle_lookup = None
1621

1722
logger = logging.getLogger(__name__)
1823

@@ -173,40 +178,24 @@ def trace_cubes(self, beam: Union[dict[str, Any], Beam]) -> list[sitk.Image]:
173178
)
174179
ray_matrix_scale = 1 + ray_matrix_bev_y / beam.sad
175180

176-
# num_candidate_rays = 2 * np.ceil(500.0 / ray_spacing).astype(np.int64) + 1
177-
178-
spacing_range = ray_spacing * np.arange(
179-
np.floor(-500.0 / ray_spacing), np.ceil(500.0 / ray_spacing) + 1, dtype=self.precision
180-
)
181-
candidate_ray_coords_x, candidate_ray_coords_z = np.meshgrid(spacing_range, spacing_range)
182-
183181
# If we have reference positions, we use them to restrict the raytracing region
184182
reference_positions_bev = ray_matrix_scale * np.array(
185183
[ray.ray_pos_bev for ray in beam.rays]
186184
)
187185

188-
# use a precompiled numba function to speed up the spatial lookup
189-
candidate_ray_mx = fast_spatial_circle_lookup(
190-
candidate_ray_coords_x,
191-
candidate_ray_coords_z,
192-
reference_positions_bev,
193-
self.lateral_cut_off,
186+
spacing_range = ray_spacing * np.arange(
187+
np.floor(-500.0 / ray_spacing), np.ceil(500.0 / ray_spacing) + 1, dtype=self.precision
194188
)
195189

196-
# candidate_ray_mx = np.full(candidate_ray_coords_x.shape, False, dtype=np.bool)
190+
candidate_ray_mx = self._get_candidate_ray_matrix(spacing_range, reference_positions_bev)
197191

198-
# for i in range(reference_positions_bev.shape[0]):
199-
# notix = (candidate_ray_coords_x - reference_positions_bev[i, 0]) ** 2 + (
200-
# candidate_ray_coords_z - reference_positions_bev[i, 2]
201-
# ) ** 2 <= self.lateral_cut_off**2
202-
# candidate_ray_mx[notix] = True
192+
ray_idx_z, ray_idx_x = np.nonzero(candidate_ray_mx)
203193

204-
ray_matrix_bev = np.hstack(
194+
ray_matrix_bev = np.column_stack(
205195
(
206-
candidate_ray_coords_x[candidate_ray_mx].reshape(-1, 1),
207-
ray_matrix_bev_y
208-
* np.ones(np.sum(candidate_ray_mx), dtype=self.precision).reshape(-1, 1),
209-
candidate_ray_coords_z[candidate_ray_mx].reshape(-1, 1),
196+
spacing_range[ray_idx_x],
197+
np.full(ray_idx_x.shape[0], ray_matrix_bev_y, dtype=self.precision),
198+
spacing_range[ray_idx_z],
210199
)
211200
)
212201

@@ -301,6 +290,57 @@ def trace_cubes(self, beam: Union[dict[str, Any], Beam]) -> list[sitk.Image]:
301290
return rad_depth_cubes
302291
# scale_factor[valid_ix] = lengths[valid_ix] / d12[valid_ix]
303292

293+
def _get_candidate_ray_matrix(self, spacing_range: Array, ref_pos_bev: Array) -> Array:
294+
"""Get candidate ray matrix for given ray spacing and reference positions."""
295+
296+
xp = array_api_compat.array_namespace(spacing_range, ref_pos_bev)
297+
298+
# Use numba accelerated code if possible
299+
if array_api_compat.is_numpy_namespace(xp) and fast_spatial_circle_lookup is not None:
300+
candidate_ray_coords_x, candidate_ray_coords_z = np.meshgrid(
301+
spacing_range, spacing_range
302+
)
303+
return fast_spatial_circle_lookup(
304+
candidate_ray_coords_x, candidate_ray_coords_z, ref_pos_bev, self.lateral_cut_off
305+
)
306+
307+
# Array API compliant code
308+
n_candidates = array_api_compat.size(spacing_range)
309+
310+
candidate_ray_mx = xp.zeros((n_candidates, n_candidates), dtype=xp.bool)
311+
312+
r2 = self.lateral_cut_off**2
313+
314+
# use buffers to avoid repeated allocations in the loop
315+
buffer_x = xp.empty_like(spacing_range)
316+
buffer_z = xp.empty_like(spacing_range)
317+
318+
for i in range(ref_pos_bev.shape[0]):
319+
# Use in-place operations as much as possible
320+
buffer_x[:] = spacing_range
321+
buffer_z[:] = spacing_range
322+
323+
buffer_x -= ref_pos_bev[i, 0]
324+
buffer_z -= ref_pos_bev[i, 2]
325+
326+
buffer_x *= buffer_x
327+
buffer_z *= buffer_z
328+
329+
# simple full boolean or
330+
# candidate_ray_mx |= (buffer_x[:, None] + buffer_z[None, :) <= r2
331+
332+
# reduce update region by finding z ranges that are valid
333+
z_ok = xp.astype(buffer_z <= r2, xp.uint8) # argmax later only works on numeric types
334+
if not xp.any(z_ok):
335+
continue
336+
337+
z0 = xp.argmax(z_ok)
338+
z1 = array_api_compat.size(z_ok) - xp.argmax(z_ok[::-1])
339+
340+
candidate_ray_mx[z0:z1, :] |= (buffer_z[z0:z1, None] + buffer_x[None, :]) <= r2
341+
342+
return candidate_ray_mx
343+
304344
@abstractmethod
305345
def _initialize_geometry(self):
306346
"""

pyRadPlan/raytracer/_numba_perf.py

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
from numba import njit, prange
2+
import numpy as np
3+
4+
5+
@njit(parallel=True, nogil=True, cache=True)
6+
def fast_spatial_circle_lookup(
7+
mesh_x: np.ndarray, mesh_z: np.ndarray, lookup_pos: np.ndarray, radius: float
8+
) -> np.ndarray:
9+
"""Lookup all points within a circle in a meshgrid.
10+
11+
Parameters
12+
----------
13+
mesh_x : np.ndarray
14+
The x-coordinates of the meshgrid (NxN)
15+
mesh_z : np.ndarray
16+
The y-coordinates of the meshgrid (NxN)
17+
lookup_pos : np.ndarray
18+
the reference positions to lookup (Mx3)
19+
radius : float
20+
The radius of the circle.
21+
22+
Returns
23+
-------
24+
np.ndarray
25+
The indices of the points within the circle.
26+
"""
27+
candidate_ray_mx = np.full(mesh_x.shape, False, dtype=np.bool_)
28+
for i in prange(lookup_pos.shape[0]):
29+
ix = (mesh_x.ravel() - lookup_pos[i, 0]) ** 2 + (
30+
mesh_z.ravel() - lookup_pos[i, 2]
31+
) ** 2 <= radius**2
32+
candidate_ray_mx.ravel()[ix] = 1
33+
return candidate_ray_mx

0 commit comments

Comments
 (0)