Port solution analysis functionality from matlab (OpenwaterHealth#105)

* [WIP] finishing up solution analysis

* Fix default simulation parameters

* update masking

* fix TIC calculation

* Fix ruff errors

* remove stale references

* Style changes to fix CI

* Fix broken imports

* remove spurious file

* fix pulse sequences in unit tests

* skip the unit tests that still need to be revamped

* document and test get_effective_origin

* add a few more docstrings

* add a few more docstrings

* remove unused function mask_focus

* add more docstrings

* make get_beam_bounds return tuple

Makes type annotation more clear

* add more docstrings

* migrate the offset grid test to the new get_offset_grid

unit test rescued from the graveyard!

---------

Co-authored-by: Ebrahim Ebrahim <ebrahim.ebrahim@kitware.com>

Author: peterhollender

notebooks/test_first.py

@@ -92,7 +92,7 @@
 ds['p_min'].sel(ele=0).plot.imshow()
 
 # %% [markdown]
-# We can examine the output object, which is an `xarray.DataSet` object with 3 data variables: `p_max` (Peak Positive Pressure), `p_min` (Peak Negative Pressure), and `ita` (Time Averaged Intensity). It's attributes also contain the `source` pulse (an `xarray.DataArray`), and `output`, the raw K-Wave output structure.
+# We can examine the output object, which is an `xarray.DataSet` object with 3 data variables: `p_max` (Peak Positive Pressure), `p_min` (Peak Negative Pressure), and `intensity` (Time Averaged Intensity). It's attributes also contain the `source` pulse (an `xarray.DataArray`), and `output`, the raw K-Wave output structure.
 
 # %%
 ds

src/openlifu/__init__.py

@@ -15,9 +15,6 @@
     apod_methods,
     delay_methods,
     focal_patterns,
-    get_beamwidth,
-    mask_focus,
-    offset_grid,
 )
 from openlifu.db import Database, User
 
@@ -63,9 +60,6 @@
     "focal_patterns",
     "delay_methods",
     "apod_methods",
-    "get_beamwidth",
-    "mask_focus",
-    "offset_grid",
     "SimSetup",
     "Database",
     "User",

src/openlifu/bf/__init__.py

@@ -1,9 +1,6 @@
 from .apod_methods import ApodizationMethod
 from .delay_methods import DelayMethod
 from .focal_patterns import FocalPattern, SinglePoint, Wheel
-from .get_beamwidth import get_beamwidth
-from .mask_focus import mask_focus
-from .offset_grid import offset_grid
 from .pulse import Pulse
 from .sequence import Sequence
 
@@ -14,8 +11,5 @@
     "FocalPattern",
     "SinglePoint",
     "Pulse",
-    "Sequence",
-    "offset_grid",
-    "mask_focus",
-    "get_beamwidth"
+    "Sequence"
 ]

src/openlifu/bf/get_beamwidth.py

@@ -20,6 +20,18 @@ class Sequence(DictMixin):
     pulse_train_interval: float = 1.0 # s
     pulse_train_count: int = 1
 
+    def __post_init__(self):
+        if self.pulse_interval <= 0:
+            raise ValueError("Pulse interval must be positive")
+        if self.pulse_count <= 0:
+            raise ValueError("Pulse count must be positive")
+        if self.pulse_train_interval < 0:
+            raise ValueError("Pulse train interval must be non-negative")
+        elif (self.pulse_train_interval > 0) and (self.pulse_train_interval < (self.pulse_interval * self.pulse_count)):
+            raise ValueError("Pulse train interval must be greater than or equal to the total pulse interval")
+        if self.pulse_train_count <= 0:
+            raise ValueError("Pulse train count must be positive")
+
     def get_table(self) -> pd.DataFrame:
         """
         Get a table of the sequence parameters

src/openlifu/bf/mask_focus.py

@@ -11,11 +11,11 @@
 
 @dataclass
 class Point:
+    position: np.ndarray = field(default_factory=lambda: np.array([0.0, 0.0, 0.0])) # mm
     id: str = "point"
     name: str = "Point"
     color: Any = (1.0, 0.0, 0.0)
     radius: float = 1. # mm
-    position: np.ndarray = field(default_factory=lambda: np.array([0.0, 0.0, 0.0])) # mm
     dims: Tuple[str, str, str] = ("x","y","z")
     units: str = "mm"
 

src/openlifu/bf/offset_grid.py

@@ -252,7 +252,7 @@ def calc_solution(
         simulation_result_aggregated: xa.Dataset = xa.Dataset()
         scaled_solution_analysis: SolutionAnalysis = SolutionAnalysis()
         foci: List[Point] = self.focal_pattern.get_targets(target)
-        simulation_cycles = np.max([np.round(self.pulse.duration * self.pulse.frequency), 20])
+        simulation_cycles = np.min([np.round(self.pulse.duration * self.pulse.frequency), 20])
 
         # updating solution sequence if pulse mismatch
         if (self.sequence.pulse_count % len(foci)) != 0:
@@ -273,7 +273,7 @@ def calc_solution(
                     cycles = simulation_cycles,
                     dt=sim_options.dt,
                     t_end=sim_options.t_end,
-                    amplitude = 1,
+                    amplitude = self.pulse.amplitude,
                     gpu = use_gpu
                 )
             delays_to_stack.append(delays)
@@ -325,11 +325,11 @@ def calc_solution(
             pnp_aggregated = solution.simulation_result['p_min'].max(dim="focal_point_index")
             ppp_aggregated = solution.simulation_result['p_max'].max(dim="focal_point_index")
             # TODO: Ensure this mean is weighted by the number of times each point is focused on, once openlifu supports hitting points different numbers of times
-            intensity_aggregated = solution.simulation_result['ita'].mean(dim="focal_point_index")
+            intensity_aggregated = solution.simulation_result['intensity'].mean(dim="focal_point_index")
             simulation_result_aggregated = deepcopy(solution.simulation_result)
             simulation_result_aggregated = simulation_result_aggregated.drop_dims("focal_point_index")
             simulation_result_aggregated['p_min'] = pnp_aggregated
             simulation_result_aggregated['p_max'] = ppp_aggregated
-            simulation_result_aggregated['ita'] = intensity_aggregated
+            simulation_result_aggregated['intensity'] = intensity_aggregated
 
         return solution, simulation_result_aggregated, scaled_solution_analysis