Merge branch 'develop' into 'main'

Update main to version 0.2.3

See merge request e040/e0404/pyRadPlan!73

Author: wahln

examples/proton_let_optimization.py

@@ -0,0 +1,61 @@
+import logging
+
+import numpy as np
+
+from pyRadPlan import (
+    IonPlan,
+    generate_stf,
+    calc_dose_influence,
+    fluence_optimization,
+    plot_slice,
+    load_tg119,
+)
+
+from pyRadPlan.optimization.objectives import SquaredDeviation, SquaredOverdosing, MeanDose
+
+logging.basicConfig(level=logging.INFO)
+
+#  Read patient from provided TG119.mat file and validate data
+ct, cst = load_tg119()
+
+# Create a plan object
+pln = IonPlan(radiation_mode="protons", machine="Generic")
+pln.prop_opt = {"solver": "scipy"}
+
+# Generate Steering Geometry ("stf")
+stf = generate_stf(ct, cst, pln)
+
+# Calculate Dose Influence Matrix ("dij")
+dij = calc_dose_influence(ct, cst, stf, pln)
+
+# Optimization
+cst.vois[1].objectives = [SquaredDeviation(priority=100.0, d_ref=3.0)]  # Target
+cst.vois[0].objectives = [
+    SquaredOverdosing(priority=10.0, d_max=1.0),
+    SquaredOverdosing(priority=5.0, d_max=1.0, quantity="let_dose"),
+]  # OAR
+cst.vois[2].objectives = [
+    MeanDose(priority=1.0, d_ref=0.0),
+    SquaredOverdosing(priority=10.0, d_max=2.0),
+]  # BODY
+fluence = fluence_optimization(ct, cst, stf, dij, pln)
+
+# Result
+result = dij.compute_result_ct_grid(fluence)
+
+# Choose a slice to visualize
+view_slice = int(np.round(ct.size[2] / 2))
+
+# Visualize
+plot_slice(
+    ct=ct,
+    cst=cst,
+    overlay=result["physical_dose"],
+    view_slice=view_slice,
+)
+plot_slice(
+    ct=ct,
+    cst=cst,
+    overlay=result["let"],
+    view_slice=view_slice,
+)

pyRadPlan/core/_grids.py

@@ -76,9 +76,9 @@ def _check_resolution(cls, value: dict[str, float]) -> dict[str, float]:
             raise ValueError("resolution must have keys 'x', 'y', 'z'")
 
         # Check if all resolution values are positive floats. If not, try to cast them to float
-        for key in value:
-            if value[key] <= 0:
-                raise ValueError(f"resolution values must be positive floats, got {value[key]}")
+        for _, v in value.items():
+            if v <= 0:
+                raise ValueError(f"resolution values must be positive floats, got {v}")
 
         return value
 
@@ -93,11 +93,12 @@ def _check_dimensions(cls, value: tuple[int, int, int]) -> tuple[int, int, int]:
         # Check if all dimensions values are positive integers. If not, try to cast them to int
         for dim in value:
             try:
-                dim = int(dim)
+                tmpdim = int(dim)
             except ValueError:
                 raise ValueError(f"dimension value could not be casted into int, got {dim}")
-            if dim <= 0:
-                raise ValueError(f"dimension values must be positive integers, got {dim}")
+
+            if tmpdim <= 0:
+                raise ValueError(f"dimension values must be positive integers, got {tmpdim}")
 
         return value
 
@@ -162,8 +163,7 @@ def resample(
         ],
     ) -> Self:
         """
-        Create a resampled grid covering the original grid in a new
-        resolution.
+        Create a resampled grid covering the original grid in new resolution.
 
         Parameters
         ----------

pyRadPlan/core/np2sitk.py

@@ -1,6 +1,4 @@
-"""Helper Functions relating to conversion between numpy and SimpleITK data
-types.
-"""
+"""Helpers for conversion between numpy and SimpleITK."""
 
 from typing import Literal
 import SimpleITK as sitk

pyRadPlan/core/resample.py

@@ -1,3 +1,5 @@
+"""Image / Grid Resampling."""
+
 import numpy as np
 import SimpleITK as sitk
 from ._grids import Grid
@@ -11,8 +13,9 @@ def resample_image(
     target_grid: Grid = None,
 ) -> sitk.Image:  # also accept sitk images and grid points. Can be dealt with with kwargs
     """
-    Resamples an sitk Image to a target resolution, reference image
-    dimensions, or Grid.
+    Resample an sitk Image.
+
+    Use target resolution, reference image, dimensions, or Grid as input.
 
     Parameters
     ----------
@@ -70,8 +73,10 @@ def resample_numpy_array(
     target_grid: Grid = None,
 ) -> np.ndarray:
     """
-    Resamples a numpy grid to a target resolution, reference image
-    dimensions, or Grid.
+    Resample a numpy grid.
+
+    Use target resolution, reference image, dimensions, or Grid as input.
+
     We convert the numpy array to an sitk.Image, so the dimensions will be
     switched. More exactly,
     the numpy array will be index i,j,k <-> z,y,x, which will converted to

pyRadPlan/ct/_ct.py

@@ -91,6 +91,9 @@ def validate_cube_hu(
             if cube_hu is None:
                 raise ValueError("HU cube not present in dictionary.")
 
+            if isinstance(cube_hu, sitk.Image) and "origin" not in data:
+                data["origin"] = cube_hu.GetOrigin()
+
             if isinstance(cube_hu, np.ndarray):
                 # Now we do check if it is matRad data
                 if data.get("cubeDim", None) is not None:
@@ -133,16 +136,6 @@ def validate_cube_hu(
             is4d = cube_hu.GetDimension() == 4
             # Now we parse the rest of the data if additional data is available in the dictionary
 
-            # TODO: Either set up a Grid or check for x/y/z to do the data management
-            if "origin" in data:
-                data["cube_hu"].SetOrigin(data["origin"])
-
-            elif all(key in data for key in ("x", "y", "z")):
-                origin = np.array([data["x"][0], data["y"][0], data["z"][0]], dtype=float)
-                if is4d:
-                    origin = np.append(origin, [0.0])
-                data["cube_hu"].SetOrigin(origin)
-
             if "direction" in data:
                 data["cube_hu"].SetDirection(data["direction"])
 
@@ -153,6 +146,23 @@ def validate_cube_hu(
                 else:
                     cube_hu.SetSpacing([resolution["x"], resolution["y"], resolution["z"]])
 
+            # TODO: Either set up a Grid or check for x/y/z to do the data management
+            if "origin" in data:
+                data["cube_hu"].SetOrigin(data["origin"])
+
+            else:
+                data["cube_hu"].SetOrigin(
+                    -np.array(data["cube_hu"].GetSize())
+                    / 2.0
+                    * np.array(data["cube_hu"].GetSpacing())
+                )
+                if is4d:
+                    data["cube_hu"].SetOrigin(np.append(data["cube_hu"].GetOrigin(), [0.0]))
+            # elif all(key in data for key in ("x", "y", "z")):
+            #     origin = np.array([data["x"][0], data["y"][0], data["z"][0]], dtype=float)
+            #     if is4d:
+            #         origin = np.append(origin, [0.0])
+            #     data["cube_hu"].SetOrigin(origin)
         return data
 
     @computed_field

pyRadPlan/dij/_dij.py

@@ -1,6 +1,6 @@
 """Contains the dij class as a (collection of) influence matrices."""
 
-from typing import Any, Union, Annotated, cast
+from typing import Any, Union, Annotated, Optional, cast
 from pydantic import (
     Field,
     field_validator,
@@ -41,15 +41,17 @@ class Dij(PyRadPlanBaseModel):
     ct_grid: Annotated[Grid, Field(default=None)]
 
     physical_dose: Annotated[NDArray, Field(default=None)]
-    let_dose: Annotated[NDArray, Field(default=None)]
+    let_dose: Annotated[NDArray, Field(default=None, alias="mLETDose")]
     alpha_dose: Annotated[NDArray, Field(default=None)]
     sqrt_beta_dose: Annotated[NDArray, Field(default=None)]
 
     num_of_beams: Annotated[int, Field(default=None)]
 
-    bixel_num: Annotated[np.ndarray, Field(default=None)]
-    ray_num: Annotated[np.ndarray, Field(default=None)]
-    beam_num: Annotated[np.ndarray, Field(default=None)]
+    bixel_num: Annotated[NDArray, Field(default=None)]
+    ray_num: Annotated[NDArray, Field(default=None)]
+    beam_num: Annotated[NDArray, Field(default=None)]
+
+    rad_depth_cubes: Optional[list[NDArray]] = Field(default=None)
 
     @computed_field
     @property
@@ -109,11 +111,30 @@ def validate_matrices(cls, v: Any, info: ValidationInfo) -> np.ndarray:
                 if mat.shape[0] != info.data["dose_grid"].num_voxels:
                     raise ValueError(f"{info.field_name} shape inconsistent with ct grid")
 
+        if info.context and "from_matRad" in info.context and info.context["from_matRad"]:
+            if v is not None:
+                for i in range(v.size):
+                    shape = (
+                        int(info.data["dose_grid"].dimensions[2]),
+                        int(info.data["dose_grid"].dimensions[0]),
+                        int(info.data["dose_grid"].dimensions[1]),
+                    )
+                    v.flat[i] = swap_orientation_sparse_matrix(
+                        v.flat[i],
+                        shape,
+                        (1, 2),  # (65, 100, 100) example
+                    )
+                    if v.flat[i] is not None and not isinstance(v.flat[i], sp.csc_matrix):
+                        v.flat[i] = sp.csc_matrix(v.flat[i])
+            else:
+                v = np.array([0])
+            return [[v]]  # TODO
+
         return v
 
     @field_validator("dose_grid", "ct_grid", mode="before")
     @classmethod
-    def validate_grid(cls, grid: Union[Grid, dict]) -> Union[Grid, dict]:
+    def validate_grid(cls, grid: Union[Grid, dict], info: ValidationInfo) -> Union[Grid, dict]:
         """
         Validate grid dictionaries.
 
@@ -123,7 +144,14 @@ def validate_grid(cls, grid: Union[Grid, dict]) -> Union[Grid, dict]:
         """
         # Check if it is a dictionary and then try to create a Grid object
         if isinstance(grid, dict):
-            grid = Grid.model_validate(grid)
+            if info.context and "from_matRad" in info.context and info.context["from_matRad"]:
+                grid["dimensions"] = np.array(
+                    [grid["dimensions"][1], grid["dimensions"][0], grid["dimensions"][2]]
+                )
+                # TODO: might swap offset and resolution
+                grid = Grid.model_validate(grid)
+            else:
+                grid = Grid.model_validate(grid)
         return grid
 
     @field_validator("beam_num", mode="before")
@@ -184,7 +212,6 @@ def grid_serializer(
         context = info.context
         if context and context.get("matRad") == "mat-file":
             return value.to_matrad(context=context["matRad"])
-
         return handler(value, info)
 
     @field_serializer("physical_dose", "let_dose", "alpha_dose", "sqrt_beta_dose")
@@ -204,6 +231,18 @@ def physical_dose_serializer(self, value: np.ndarray, info: SerializationInfo) -
                 )
                 if value.flat[i] is not None and not isinstance(value.flat[i], sp.csc_matrix):
                     value.flat[i] = sp.csc_matrix(value.flat[i])
+        # return 0 if value is None. savemat() cant handle 'None'
+        elif context and context.get("matRad") == "mat-file" and value is None:
+            value = np.array([0])
+        return value
+
+    @field_serializer("rad_depth_cubes")
+    def rad_depth_cubes_serializer(self, value: np.ndarray, info: SerializationInfo) -> np.ndarray:
+        context = info.context
+        if context and context.get("matRad") == "mat-file" and value is not None:
+            # TODO: it might be necessary to rotate the cube!
+            return value
+        # return 0 if value is None. savemat() cant handle 'None'
         elif context and context.get("matRad") == "mat-file" and value is None:
             value = np.array([0])
         return value
@@ -253,7 +292,7 @@ def get_result_arrays_from_intensity(
             if self.physical_dose is None:
                 raise ValueError("Physical dose must be calculated for dose-weighted let")
 
-            indices = out["physical_dose"] > 0.0
+            indices = out["physical_dose"] > 0.05 * np.max(out["physical_dose"])
 
             let_dose = self.let_dose.flat[scenario_index] @ intensity
             out["let"] = np.zeros_like(let_dose)

pyRadPlan/dose/engines/_base_pencilbeam_particle.py

@@ -5,6 +5,7 @@
 from typing import cast, Literal, Any
 import logging
 import time
+from copy import deepcopy
 
 import numpy as np
 from scipy.interpolate import interp1d
@@ -735,6 +736,8 @@ def _calc_lateral_particle_cut_off(self, cut_off_level, stf_element):
                 comp_fac=1.0, depths=depth_values, cut_off=np.inf * np.ones_like(depth_values)
             )
 
+            base_kernel = deepcopy(kernel)  # similar to base_data in matRad
+
             # TODO: this could probably be done for all depths with one calculation without a loop
             for j, current_depth in enumerate(depth_values):
                 # If there's no cut-off set, we do not need to find it. Set it to infinity
@@ -744,12 +747,13 @@ def _calc_lateral_particle_cut_off(self, cut_off_level, stf_element):
                 # depth
                 bixel = {
                     "energy_ix": energy_ix,
-                    "kernel": kernel,
+                    "kernel": base_kernel,
                     "radial_dist_sq": radial_dist_sq,
                     "sigma_ini_sq": largest_sigma_sq4unique_energies[
                         cnt - 1
                     ],  # TODO: check if this result is correct (rounded but may be right)
-                    "rad_depths": (current_depth + kernel.offset) * np.ones_like(radial_dist_sq),
+                    "rad_depths": (current_depth + base_kernel.offset)
+                    * np.ones_like(radial_dist_sq),
                     "v_tissue_index": np.ones_like(radial_dist_sq),
                     "v_alpha_x": 0.5 * np.ones_like(radial_dist_sq),
                     "v_beta_x": 0.05 * np.ones_like(radial_dist_sq),

pyRadPlan/dose/engines/_svdpb.py

@@ -299,7 +299,7 @@ def _compute_bixel(self, curr_ray: dict[str], k: int) -> dict[str, Any]:
         bixel_dose *= ((sad) / geo_depths) ** 2
 
         bixel["physical_dose"] = bixel_dose
-
+        bixel["weight"] = curr_ray["beamlets"][k]["weight"]
         bixel["ix"] = curr_ray["ix"]
 
         return bixel

pyRadPlan/optimization/problems/_optiprob.py

@@ -191,6 +191,8 @@ def _initialize(self):
                 )
 
         # TODO: manage scenarios
+
+        # Manage quantites by getting them from the objective quantities
         self._quantities = [q(self._dij) for q in quantities]
 
         # obtain cache info to match quantities with objectives