Merged in feature/RAM-3022_wl_bb_location (pull request #277)

Feature/RAM-3022 WL BB location plot

Approved-by: Randy Taylor

docs/source/changelog.rst

@@ -13,6 +13,12 @@ Planar Imaging
   ROIs. This result will also be included in the ``results_data`` structure.
   This is not done for light/rad phantoms.
 
+Winston-Lutz
+^^^^^^^^^^^^
+
+* Normal Winston-Lutz analyses (not multi-target/multi-field) can now plot a visualization of the BB position
+  relative to the determined isocenter. After analyzing a WL set, call ``plot_location()``. See :ref:`wl_visualizing_bb`.
+
 v 3.16.0
 --------
 

docs/source/winston_lutz.rst

@@ -166,7 +166,25 @@ You can also pass in your couch coordinates and the new values will be generated
     print(wl.bb_shift_instructions(couch_vrt=0.41, couch_lng=96.23, couch_lat=0.12))
     # New couch coordinates (mm): VRT: 0.32; LNG: 96.11; LAT: 0.11
 
-.. _using_file_names_wl:
+.. _wl_visualizing_bb:
+
+Visualizing the Isocenter-to-BB
+-------------------------------
+
+.. versionadded:: 3.17
+
+The isocenter and BB can be visualized together after analysis by calling ``plot_location``:
+
+.. code-block:: python
+
+    wl = WinstonLutz(...)
+    wl.analyze(...)
+    wl.plot_location()
+
+This will result in a 3D plot visualizing the BB (true physical size) and the isocenter (true physical size) in the room coordinates
+like so:
+
+.. image:: images/iso_visualization.png
 
 
 Accessing data
@@ -239,6 +257,8 @@ You may optionally analyze a single image if that is your preference. Obviously,
 
 This class does not have all the methods that ``WinstonLutz`` has for mostly obvious reasons and lower likelihood of being used directly.
 
+.. _using_file_names_wl:
+
 .. _passing-in-axis-values:
 
 Passing in Axis values

pylinac/core/geometry.py

@@ -10,6 +10,7 @@
 import numpy as np
 from matplotlib.patches import Circle as mpl_Circle
 from matplotlib.patches import Rectangle as mpl_Rectangle
+from mpl_toolkits.mplot3d.art3d import Line3D
 
 from .utilities import is_iterable
 
@@ -399,7 +400,7 @@ def distance_to(self, point: Point) -> float:
 
     def plot2axes(
         self, axes: plt.Axes, width: float = 1, color: str = "w", **kwargs
-    ) -> None:
+    ) -> Line3D:
         """Plot the line to an axes.
 
         Parameters
@@ -409,14 +410,15 @@ def plot2axes(
         color : str
             The color of the line.
         """
-        axes.plot(
+        lines = axes.plot(
             (self.point1.x, self.point2.x),
             (self.point1.y, self.point2.y),
             (self.point1.z, self.point2.z),
             linewidth=width,
             color=color,
             **kwargs,
         )
+        return lines[0]
 
 
 class Rectangle:

pylinac/core/image_generator/utils.py

@@ -213,7 +213,10 @@ def generate_winstonlutz(
                 alpha=field_alpha,
             )
         )
-        long_offset = bb_projection_long(
+        # we return the negative because this function
+        # will return the offset in PLOTTING space, not coordinate space
+        # which is inverted in the long direction
+        long_offset = -bb_projection_long(
             offset_in=offset_mm_in,
             offset_up=offset_mm_up,
             offset_left=offset_mm_left,

pylinac/winston_lutz.py

@@ -29,10 +29,11 @@
 import tempfile
 import webbrowser
 from dataclasses import dataclass
+from functools import cached_property
 from itertools import zip_longest
 from pathlib import Path
 from textwrap import wrap
-from typing import BinaryIO, Iterable, Sequence
+from typing import BinaryIO, Iterable, Sequence, TypedDict
 
 import argue
 import matplotlib.pyplot as plt
@@ -148,6 +149,81 @@ def to_human(arrangement: dict) -> str:
         return f"'{a['name']}': {lr} {abs(a['offset_left_mm'])}mm, {ud} {abs(a['offset_up_mm'])}mm, {io} {abs(a['offset_in_mm'])}mm"
 
 
+class NominalBB(TypedDict):
+    """Input for BB location"""
+
+    offset_left_mm: float
+    offset_up_mm: float
+    offset_in_mm: float
+    bb_diameter_mm: float
+
+
+class BB:
+    """A representation of a BB in 3D space"""
+
+    def __repr__(self):
+        return self.nominal_position
+
+    def __init__(self, nominal_bb: NominalBB, ray_lines: list[Line]):
+        self.nominal_bb = nominal_bb
+        self.ray_lines = ray_lines
+
+    @cached_property
+    def measured_position(self) -> Point:
+        """The 3D measured position of the BB based on the ray-tracing lines in MM"""
+        initial_guess = self.nominal_position.as_array()
+        bounds = [(-200, 200), (-200, 200), (-200, 200)]
+        result = optimize.minimize(
+            max_distance_to_lines, initial_guess, args=self.ray_lines, bounds=bounds
+        )
+        return Point(result.x)
+
+    @cached_property
+    def nominal_position(self) -> Point:
+        """The nominal location of the BB in MM"""
+        return Point(
+            x=-self.nominal_bb["offset_left_mm"],
+            y=-self.nominal_bb["offset_in_mm"],
+            z=self.nominal_bb["offset_up_mm"],
+        )
+
+    @cached_property
+    def delta_vector(self) -> Vector:
+        """The shift from measured BB location to nominal as a vector in MM"""
+        return self.measured_position - self.nominal_position
+
+    @cached_property
+    def delta_distance(self):
+        """The scalar distance between the measured BB location and nominal in MM"""
+        return self.measured_position.distance_to(self.nominal_position)
+
+    def plot_nominal(self, axes: plt.Axes, color: str):
+        """Plot the BB nominal position"""
+        x, y, z = create_sphere_surface(
+            radius=self.nominal_bb["bb_diameter_mm"] / 2, center=self.nominal_position
+        )
+        axes.plot_surface(x, y, z, color=color)
+
+    def plot_measured(self, axes: plt.Axes, color: str):
+        """Plot the BB measured position"""
+        x, y, z = create_sphere_surface(
+            radius=self.nominal_bb["bb_diameter_mm"] / 2, center=self.measured_position
+        )
+        axes.plot_surface(x, y, z, color=color)
+
+
+def create_sphere_surface(
+    radius: float, center: Point
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Create a sphere surface for plotting"""
+    u = np.linspace(0, 2 * np.pi, 100)
+    v = np.linspace(0, np.pi, 100)
+    x = radius * np.outer(np.cos(u), np.sin(v)) + center.x
+    y = radius * np.outer(np.sin(u), np.sin(v)) + center.y
+    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) + center.z
+    return x, y, z
+
+
 class Axis(enum.Enum):
     GANTRY = "Gantry"  #:
     COLLIMATOR = "Collimator"  #:
@@ -526,18 +602,18 @@ def cax_line_projection(self) -> Line:
         p1.x = self.cax2bb_vector.x * cos(self.gantry_angle) + 20 * sin(
             self.gantry_angle
         )
-        p1.y = self.cax2bb_vector.x * -sin(self.gantry_angle) + 20 * cos(
+        p1.z = self.cax2bb_vector.x * -sin(self.gantry_angle) + 20 * cos(
             self.gantry_angle
         )
-        p1.z = self.cax2bb_vector.y
+        p1.y = self.cax2bb_vector.y
         # point 2 - ray destination
         p2.x = self.cax2bb_vector.x * cos(self.gantry_angle) - 20 * sin(
             self.gantry_angle
         )
-        p2.y = self.cax2bb_vector.x * -sin(self.gantry_angle) - 20 * cos(
+        p2.z = self.cax2bb_vector.x * -sin(self.gantry_angle) - 20 * cos(
             self.gantry_angle
         )
-        p2.z = self.cax2bb_vector.y
+        p2.y = self.cax2bb_vector.y
         line = Line(p1, p2)
         return line
 
@@ -660,6 +736,7 @@ class WinstonLutz:
     machine_scale: MachineScale  #:
     image_type = WinstonLutz2D
     is_from_cbct: bool = False
+    _bb_diameter = float
 
     def __init__(
         self,
@@ -923,7 +1000,7 @@ def analyze(
         Parameters
         ----------
         bb_size_mm
-            The expected size of the BB in mm. The actual size of the BB can be +/-2mm from the passed value.
+            The expected diameter of the BB in mm. The actual size of the BB can be +/-2mm from the passed value.
         machine_scale
             The scale of the machine. Shift vectors depend on this value.
         low_density_bb
@@ -940,6 +1017,7 @@ def analyze(
         for img in self.images:
             img.analyze(bb_size_mm, low_density_bb, open_field)
         self._is_analyzed = True
+        self._bb_diameter = bb_size_mm
 
     @lru_cache()
     def _minimize_axis(self, axes: Axis | tuple[Axis, ...] = (Axis.GANTRY,)):
@@ -1257,6 +1335,96 @@ def plot_axis_images(
         if show:
             plt.show()
 
+    def plot_location(
+        self,
+        show: bool = True,
+        viewbox_mm: float | None = None,
+        plot_bb: bool = True,
+        plot_isocenter_sphere: bool = True,
+    ):
+        """Plot the isocenter and size as a sphere in 3D space relative to the BB. The
+        iso is at the origin.
+
+        Only images where the couch was at zero are considered.
+
+        Parameters
+        ----------
+        show : bool
+            Whether to plot the image.
+        viewbox_mm : float
+            The default size of the 3D space to plot in mm in each axis.
+        plot_bb : bool
+            Whether to plot the BB location; the size is also considered.
+        plot_isocenter_sphere : bool
+            Whether to plot the gantry + collimator isocenter size.
+        """
+        limit = (
+            viewbox_mm
+            or max(
+                np.abs(
+                    (
+                        self.bb_shift_vector.x,
+                        self.bb_shift_vector.y,
+                        self.bb_shift_vector.z,
+                    )
+                )
+            )
+            + self._bb_diameter
+        )
+        ax = plt.axes(projection="3d")
+        _, relevant_images = self._get_images(
+            axis=(Axis.REFERENCE, Axis.GB_COMBO, Axis.COLLIMATOR, Axis.GANTRY)
+        )
+        # we can represent the iso sphere as a BB object; the nominal object isn't used, just the BB size
+        # the ray lines are what we want to plot as a sphere
+        iso_sphere = BB(
+            nominal_bb={
+                "offset_left_mm": 0,
+                "offset_in_mm": 0,
+                "offset_up_mm": 0,
+                "bb_diameter_mm": self._bb_diameter,
+            },
+            ray_lines=[image.cax_line_projection for image in relevant_images],
+        )
+        # plot the x,y,z origin lines
+        x_line = Line(Point(-limit, 0, 0), Point(limit, 0, 0))
+        x_line.plot2axes(ax, color="red", label="isocenter (x)")
+        y_line = Line(Point(0, -limit, 0), Point(0, limit, 0))
+        y_line.plot2axes(ax, color="green", label="isocenter (y)")
+        z_line = Line(Point(0, 0, -limit), Point(0, 0, limit))
+        z_line.plot2axes(ax, color="blue", label="isocenter (z)")
+        if plot_bb:
+            iso_sphere.plot_measured(ax, color="cyan")
+            # create an empty, fake line so we can add a label for the legend
+            fake_line = Line(Point(0, 0, 0), Point(0, 0, 0))
+            fake_line.plot2axes(ax, color="cyan", label=f"BB ({self._bb_diameter}mm)")
+        if plot_isocenter_sphere:
+            x, y, z = create_sphere_surface(
+                radius=self.gantry_coll_iso_size / 2, center=Point(0, 0, 0)
+            )
+            ax.plot_surface(x, y, z, alpha=0.2, color="magenta")
+            # create an empty, fake line so we can add a label for the legend
+            fake_line = Line(Point(0, 0, 0), Point(0, 0, 0))
+            fake_line.plot2axes(
+                ax,
+                color="magenta",
+                label=f"Isocenter sphere ({self.gantry_coll_iso_size:3.2f}mm)",
+            )
+        ax.legend()
+        # set the limits of the 3D plot; they must be the same in all axes for equal aspect ratio
+        ax.set(
+            xlabel="X (mm), Right (+)",
+            ylabel="Y (mm), In (+)",
+            zlabel="Z (mm), Up (+)",
+            title="Isocenter Visualization",
+            ylim=[-limit, limit],
+            xlim=[-limit, limit],
+            zlim=[-limit, limit],
+        )
+
+        if show:
+            plt.show()
+
     def plot_images(
         self,
         axis: Axis = Axis.GANTRY,
@@ -1750,7 +1918,6 @@ def _find_bb(self, bb_of_interest: dict) -> Point:
                 max_thresh -= 0.03 * spread
                 if max_thresh < hmin:
                     raise ValueError(BB_ERROR_MESSAGE)
-
         # determine the center of mass of the BB
         inv_img = image.load(self.array)
         # we invert so BB intensity increases w/ attenuation
@@ -1828,7 +1995,7 @@ def run_demo():
         print(wl.results())
         wl.plot_images()
 
-    def analyze(self, bb_arrangement: Iterable[dict]):
+    def analyze(self, bb_arrangement: Sequence[dict]):
         """Analyze the WL images.
 
         Parameters
@@ -1845,7 +2012,7 @@ def analyze(self, bb_arrangement: Iterable[dict]):
                 try:
                     image_set.append(img.as_analyzed(bb))
                 except ValueError:
-                    pass  # didn't find the field and/or BB; likely occluded or not part of the plan
+                    pass
             if not image_set:
                 raise ValueError(f"Did not find any field/bb pairs for bb: {bb}")
             self.analyzed_images[BBArrangement.to_human(bb)] = image_set
@@ -2130,7 +2297,7 @@ def _bb_projection_with_rotation(
     gantry: float,
     couch: float = 0,
     sad: float = 1000,
-) -> np.ndarray:
+) -> (float, float):
     """Calculate the isoplane projection onto the panel at the given SSD.
 
     This function applies a rotation around the gantry plane (X/Z) to the
@@ -2171,4 +2338,4 @@ def _bb_projection_with_rotation(
     imager_projection = (
         np.array([rotated_positions[1], rotated_positions[2]]) * bb_magnification
     )
-    return imager_projection
+    return imager_projection[0], imager_projection[1]