Merged in refactor/finally_explain_low (pull request #365)

finally fully deconvolve the Low paper.

Approved-by: Randy Taylor

docs/source/winston_lutz.rst

@@ -535,7 +535,6 @@ Couch shift
 This method is used to determine the shift instructions. Specifically, equations 6, 7, and 9.
 
 
-
 .. note::
 
     If doing research, it is very important to note that Low implicitly used the "Varian" coordinate system.
@@ -546,28 +545,82 @@ This method is used to determine the shift instructions. Specifically, equations
     To use a different scale use the ``machine_scale`` parameter, shown here :ref:`passing-a-coordinate-system`.
     Also see :ref:`scale`.
 
-Implementation for the couch shift is as follows:
+As a prologue, we should explain the differences in our coordinate system vs Low's. In the paper
+the gantry coordinate system is defined as such:
 
-For each image we determine (equation 6a):
+"Let :math:`(X_{G},Y_{G},Z_{G})` represent a gantry coordinate system defined as follows: the :math:`X_{G}` axis
+corresponds to the gantry rotation axis, with its positive direction pointing away from the gantry (toward the couch); the :math:`Z_{G}`
+axis coincides with the central axis of the radiation beam, pointing into the collimator, and the :math:`Y_{G}`
+axis is assigned such that :math:`(X_{G},Y_{G},Z_{G})` forms a right-handed coordinate system."
+
+Later in equation 2 they defined a couch coordinate system as:
 
 .. math::
 
-    \textbf{A}(\phi, \theta) = \begin{pmatrix} -\cos(\phi) & \sin(\phi) & 0 \\ \cos(\theta)\sin(\phi) & \cos(\theta)\cos(\phi) & -\sin(\theta) \end{pmatrix}
+    X_{c} = VERT
 
-where :math:`\theta` is the gantry angle and :math:`\phi` is the couch angle.
+    Y_{c} = LAT
 
-.. warning::
+    Z_{c} = -AP
+
+but this is inconsistent because AP is the same as VERT in normal conventions.
+Using the sentence after eqn 2 as a guide: "A couch angle of :math:`\phi` degrees corresponds
+to a signed rotation around the :math:`Z_{C}` axis of the couch coordinate system."
+
+Just from this sentence it would appear :math:`Z_{C}` should be VERT. The only way this makes
+sense is if the couch coordinate system is from a HFS patient perspective where VERT is actually
+Superior-Inferior and AP is VERT. If this is true, everything falls into place.
+
+Finally, it is worth nothing that in the Varian "Standard" coordinate system (the one assumed by Low), couch vertical
+*increases* as the couch is lowered. This may explain the negative sign in the Z axis of equation 2, whereas
+for the IEC 61217 scale, vertical increases as the couch is raised.
+
+We thus can generate the following table for transforming the equations from Low to IEC 61217/pylinac coordinates:
+
+
++-----------------------+---+----------------+---+------------------------+---+-----------------+---+----------------------------+
+| Low Couch coordinates |   | Low Couch Axes |   | Low Gantry coordinates |   | Low Gantry axes |   | Pylinac Gantry coordinates |
++=======================+===+================+===+========================+===+=================+===+============================+
+| Xc                    | = | VERT           | = | Xg                     | = | LONG            | = | -Y                         |
++-----------------------+---+----------------+---+------------------------+---+-----------------+---+----------------------------+
+| Yc                    | = | LAT            | = | Yg                     | = | LAT             | = | X                          |
++-----------------------+---+----------------+---+------------------------+---+-----------------+---+----------------------------+
+| Zc                    | = | -AP            | = | Zg                     | = | VERT            | = | Z                          |
++-----------------------+---+----------------+---+------------------------+---+-----------------+---+----------------------------+
+
+
+We can now address the implementation for the couch shift as follows:
+
+For each image we determine (equation 6a):
 
-  The Low paper appears to have incorrect signs for some of the matrix. Using synthetic images with known shifts
-  can prove the correctness of the algorithm.
+.. math::
+
+    \textbf{A}(\phi, \theta) = \begin{pmatrix} -\cos(\phi) & -\sin(\phi) & 0 \\ -\cos(\theta)\sin(\phi) & \cos(\theta)\cos(\phi) & -\sin(\theta) \end{pmatrix}
 
-The :math:`\xi` matrix is then calculated (equation 7):
+where :math:`\theta` is the gantry angle and :math:`\phi` is the couch angle (Remember, these must be in Varian Standard axes values).
+
+The :math:`\xi` matrix is then calculated with an altered equation 7:
 
 .. math::
 
-    \xi = (x_{1},y_{1}, ..., x_{i},y_{i},..., x_{n},y_{n})^{T}
+    \xi = (y_{1},-x_{1}, ..., y_{i},-x_{i},..., y_{n},-x_{n})^{T}
+
+where :math:`x_{i}` and :math:`y_{i}` are the scalar shifts from the field CAX to the BB for :math:`n` images *in our coordinate system*.
+
+.. note::
+
+   Equation 7 has :math:`(x_{i},y_{i})` convention, but we use :math:`(y_{i},-x_{i})`. Using Figure 1
+   as a guide, it says the x-axis is pointing toward the gantry (and points to the right in the figure).
+   Both the x and y axes appear inverted here. In the figure both :math:`x_{i}` and :math:`y_{i}`
+   values are presumably positive in the example, but are both negative when using the right-hand rule of Low's gantry coordinate
+   definition ("...with [:math:`X_{G}`] positive direction pointing away from the gantry").
+   Thus :math:`x_{i}` becomes :math:`-x_{i}` when using Low's own coordinate convention. Using the above table conversion this becomes :math:`y_{i}`. A similar conversion is true for the y-axis.
 
-where :math:`x_{i}` and :math:`y_{i}` are the x and y scalar shifts from the field CAX to the BB for :math:`n` images.
+   .. math::
+
+         +x_{Fig1} = -x_{Low} = y_{Pylinac}
+
+         +y_{Fig1} = -y_{Low} = -x_{Pylinac}
 
 From equation 9 we can calculate :math:`\textbf{B}(\phi_{1},\theta_{1},..., \phi_{n},\theta_{n})`:
 
@@ -583,6 +636,14 @@ We can then solve for the shift vector :math:`\boldsymbol{\Delta}`:
 
     \boldsymbol{\Delta} = \textbf{B}(\phi_{1},\theta_{1},..., \phi_{n},\theta_{n}) \cdot \xi
 
+Using our conversion table above, :math:`\boldsymbol{\Delta}` resolves to:
+
+.. math::
+
+    \boldsymbol{\Delta} = \begin{pmatrix} \Delta \text{VERT} \\ \Delta \text{LAT} \\ \Delta \text{-AP} \end{pmatrix}_{Low} = \begin{pmatrix} \Delta \text{-LONG} \\ \Delta \text{LAT} \\ \Delta \text{VERT} \end{pmatrix}_{Pylinac}
+
+where :math:`(...)_{Low}` and :math:`(...)_{Pylinac}` are in their respective coordinate systems.
+
 Implementation
 ^^^^^^^^^^^^^^
 

pylinac/winston_lutz.py

@@ -223,13 +223,21 @@ class BBFieldMatch:
 
     @property
     def bb_field_vector_mm(self) -> Vector:
-        """The vector from the BB to the field CAX"""
-        return (self.bb - self.field) / self.dpmm
+        """The vector from the BB to the field CAX *IN COORDINATE SPACE*."""
+        v = (self.bb - self.field) / self.dpmm
+        v.y = (
+            -v.y
+        )  # invert the y-axis; positive is down in image space but negative in coordinate space
+        return v
 
     @property
     def bb_epid_vector_mm(self) -> Vector:
-        """The vector from the BB to the field CAX"""
-        return (self.bb - self.epid) / self.dpmm
+        """The vector from the BB to the field CAX *IN COORDINATE SPACE*."""
+        v = (self.bb - self.epid) / self.dpmm
+        v.y = (
+            -v.y
+        )  # invert the y-axis; positive is down in image space but negative in coordinate space
+        return v
 
     @property
     def bb_field_distance_mm(self) -> float:
@@ -483,7 +491,6 @@ def __init__(
             self.detection_conditions = conditions
         super().__init__(file, use_filenames=use_filenames, **kwargs)
         self._is_analyzed = False
-        self.flipud()
 
     def analyze(
         self,
@@ -529,7 +536,9 @@ def analyze(
             # convert from mm to pixels and add to the detected points
             for p in detected_bb_points:
                 p.x += lat * self.dpmm
-                p.y += sup_inf * self.dpmm
+                p.y -= (
+                    sup_inf * self.dpmm
+                )  # we subtract because the detected point is in image space, not coordinate space so we convert the shift from coordinate to image space
         bb_matches = self.find_bb_matches(detected_points=detected_bb_points)
         if len(bb_matches) != len(field_matches):
             raise ValueError("The number of detected fields and BBs do not match")
@@ -658,7 +667,7 @@ def plot(
         epid_handle = ax.axhline(y=self.epid.y, color="b")
         # show the field CAXs
         for match in self.arrangement_matches.values():
-            (field_handle,) = ax.plot(match.field.x, match.bb.y, "gs", ms=8)
+            (field_handle,) = ax.plot(match.field.x, match.field.y, "gs", ms=8)
             (bb_handle,) = ax.plot(match.bb.x, match.bb.y, "co", ms=10)
         if legend:
             ax.legend(
@@ -686,10 +695,10 @@ def plot(
                 max([match.bb.y for match in self.arrangement_matches.values()])
                 + 20 * self.dpmm
             )
-            ax.set_ylim([min_y, max_y])
+            ax.set_ylim([max_y, min_y])
             ax.set_xlim([min_x, max_x])
-        ax.set_yticklabels([])
-        ax.set_xticklabels([])
+        # ax.set_yticklabels([])
+        # ax.set_xticklabels([])
         ax.set_title("\n".join(wrap(Path(self.path).name, 30)), fontsize=10)
         ax.set_xlabel(
             f"G={self.gantry_angle:.0f}, B={self.collimator_angle:.0f}, P={self.couch_angle:.0f}"
@@ -1209,11 +1218,14 @@ def analyze(
                 img.analyze(bb_size_mm, low_density_bb, open_field, shift_vector=shift)
         bb_config = BBArrangement.ISO[0]
         bb_config.bb_size_mm = bb_size_mm
+        # in the vanilla WL case, the BB can only be represented by non-couch-kick images
+        # the ray trace cannot handle the kick currently
         self.bb = BB3D(
             bb_config=bb_config,
             ray_lines=[
                 img.arrangement_matches["Iso"].bb_to_field_projection
                 for img in self.images
+                if img.variable_axis in (Axis.GANTRY, Axis.REFERENCE)
             ],
         )
         self._is_analyzed = True
@@ -1323,25 +1335,27 @@ def bb_shift_vector(self) -> Vector:
                 rotation=img.couch_angle,
             )
             A[2 * idx : 2 * idx + 2, :] = np.array(
-                # note the signs are different than the paper; based on
-                # synthetic data we can prove this. See docs
                 [
-                    [-cos(couch), sin(couch), 0],
+                    [-cos(couch), -sin(couch), 0],
                     [
-                        cos(gantry) * sin(couch),
+                        -cos(gantry) * sin(couch),
                         cos(gantry) * cos(couch),
                         -sin(gantry),
                     ],
                 ]
-            )  # equation 6 (minus delta)
+            )  # equation 6 (minus delta portion)
             epsilon[2 * idx : 2 * idx + 2] = np.array(
-                [[img.cax2bb_vector.y], [img.cax2bb_vector.x]]
-            )  # equation 7
+                [
+                    [img.arrangement_matches["Iso"].bb_field_vector_mm.y],
+                    [-img.arrangement_matches["Iso"].bb_field_vector_mm.x],
+                ]
+            )  # equation 7; note that we have (y, -x) instead of (x, y) as in the paper. This is because
+            # of coordinate system convention differences. See docs.
 
         B = linalg.pinv(A)
-        delta = B.dot(epsilon)  # equation 9
-        # we use the negative for all values because it's from the iso POV -> linac not the linac -> iso POV
-        return Vector(x=-delta[1][0], y=-delta[0][0], z=-delta[2][0])
+        delta = B.dot(epsilon).squeeze()  # equation 9
+        # y is negative because their coordinate system is out=positive, in ours it's negative; see convestion table in docs
+        return Vector(y=-delta[0], x=delta[1], z=delta[2])
 
     def bb_shift_instructions(
         self,
@@ -1599,11 +1613,7 @@ def plot_location(
             + 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
         # plot the x,y,z origin lines
         x_line = Line(Point(-limit, 0, 0), Point(limit, 0, 0))
         x_line.plot2axes(ax, color="green", alpha=0.5)
@@ -2011,10 +2021,6 @@ class WinstonLutzMultiTargetMultiFieldImage(WLBaseImage):
     detection_conditions = [is_round, is_symmetric, is_modest_size]
     field_conditions = [is_square, is_right_square_size]
 
-    def __init__(self, *args, **kwargs):
-        super().__init__(*args, **kwargs)
-        self.flipud()  # restore to original view; vanilla WL may need to revert
-
     def find_field_centroids(self, is_open_field: bool) -> list[Point]:
         """Find the centroid of the radiation field based on a 50% height threshold.
         This applies the field detection conditions and also a nearness condition.
@@ -2168,7 +2174,7 @@ def plot_location(
         y_line.plot2axes(ax, color="green", alpha=0.5)
         z_line = Line(Point(0, 0, -100), Point(0, 0, 100))
         z_line.plot2axes(
-            ax, color="green", alpha=0.5, label="Determined isocenter (x,y,z)"
+            ax, color="green", alpha=0.5, label="Nominal isocenter (x,y,z)"
         )
         if plot_bb:
             for bb in self.bbs:

tests_basic/test_winstonlutz.py

@@ -824,6 +824,34 @@ def test_bb_shift_vector(self):
             self.wl.bb_shift_vector.z, -self.offset_mm_up, delta=0.05
         )
 
+    def test_virtual_shift_has_zero_remaining_shift(self):
+        """The virtual shift should have zero remaining shift after applying it"""
+        wl = self.new_instance()
+        wl.analyze(
+            bb_size_mm=self.bb_size,
+            machine_scale=self.machine_scale,
+            low_density_bb=self.low_density_bb,
+            open_field=self.open_field,
+            apply_virtual_shift=True,
+        )
+        self.assertAlmostEqual(wl.bb_shift_vector.as_scalar(), 0, delta=0.05)
+
+    def test_bb3d_measured_position(self):
+        self.assertAlmostEqual(
+            self.wl.bb.measured_position.x, -self.offset_mm_left, delta=0.03
+        )
+        self.assertAlmostEqual(
+            self.wl.bb.measured_position.y, self.offset_mm_in, delta=0.03
+        )
+        self.assertAlmostEqual(
+            self.wl.bb.measured_position.z, self.offset_mm_up, delta=0.03
+        )
+
+    def test_bb3d_nominal_position(self):
+        self.assertAlmostEqual(self.wl.bb.nominal_position.x, 0, delta=0.01)
+        self.assertAlmostEqual(self.wl.bb.nominal_position.y, 0, delta=0.01)
+        self.assertAlmostEqual(self.wl.bb.nominal_position.z, 0, delta=0.01)
+
 
 class Synthetic1mmLeftNoCouch(SyntheticWLMixin, TestCase):
     images_axes = (
@@ -939,35 +967,6 @@ class Synthetic1mmOut1SidedCouch(SyntheticWLMixin, TestCase):
     )
 
 
-class SyntheticVirtualShift(SyntheticWLMixin, TestCase):
-    """When virtually shifting, the BB should be at the isocenter and with these
-    perfect synthetic images the error is ~0.0"""
-
-    offset_mm_in = 1
-    offset_mm_left = 1
-    cax2bb_max_distance = 0
-    cax2bb_mean_distance = 0
-    cax2bb_median_distance = 0
-    couch_iso_size = 0
-    images_axes = (
-        (0, 0, 0),
-        (90, 0, 0),
-        (180, 0, 0),
-        (270, 0, 0),
-        (0, 0, 45),
-        (0, 0, 90),
-        (0, 0, 270),
-        (0, 0, 315),
-    )
-    apply_virtual_shift = True
-
-    def test_bb_shift_vector(self):
-        """The vector is ~0.0 when virtually shifting"""
-        self.assertAlmostEqual(self.wl.bb_shift_vector.as_scalar(), 0, delta=0.05)
-
-    # def test_bb_size_doesnt_change_result(self):
-
-
 class WLDemo(WinstonLutzMixin, TestCase):
     num_images = 17
     gantry_iso_size = 1