diff --git a/examples/verification/tst01_console_selftest.py b/examples/verification/tst01_console_selftest.py index 9f027fce..045fabd4 100644 --- a/examples/verification/tst01_console_selftest.py +++ b/examples/verification/tst01_console_selftest.py @@ -2,8 +2,6 @@ import sys -import base58 - from openlifu_sdk.io.LIFUInterface import LIFUInterface # set PYTHONPATH=%cd%\src;%PYTHONPATH% @@ -46,8 +44,7 @@ print("Get HW ID") hw_id = interface.hvcontroller.get_hardware_id() print(f"HW ID: {hw_id}") -encoded_id = base58.b58encode(bytes.fromhex(hw_id)).decode() -print(f"OW-LIFU-CON-{encoded_id}") +print(f"OW-LIFU-CON-{hw_id}") print("Get Temperature1") temp1 = interface.hvcontroller.get_temperature1() diff --git a/notebooks/get_all_versions.py b/notebooks/get_all_versions.py new file mode 100644 index 00000000..4fb7bd1e --- /dev/null +++ b/notebooks/get_all_versions.py @@ -0,0 +1,102 @@ +from __future__ import annotations + +import logging +import os +import sys +import time + +if os.name == 'nt': + pass +else: + pass + + +from openlifu.io.LIFUInterface import LIFUInterface + +# set PYTHONPATH=%cd%\src;%PYTHONPATH% +# python notebooks/test_watertank.py + +""" +Test script to automate: +1. Connect to the device. +2. Test HVController: Turn HV on/off and check voltage. +3. Test Device functionality. +""" + +# TO BE USED TO MONITOR TEMPERATURE CURVE TO SEE HOW LONG IT TAKES TO COOL DOWN + +# Configure logging +logger = logging.getLogger(__name__) +logger.setLevel(logging.INFO) + +# Prevent duplicate handlers and cluttered terminal output +if not logger.hasHandlers(): + handler = logging.StreamHandler() + handler.setFormatter(logging.Formatter("%(asctime)s - %(levelname)s - %(message)s")) + logger.addHandler(handler) + logger.propagate = False + +log_interval = 1 # seconds; you can adjust this variable as needed +num_modules = 2 # Number of modules in the system + +use_external_power_supply = False # Select whether to use console or power supply + +logger.info("Starting LIFU Test Script...") +interface = LIFUInterface(ext_power_supply=use_external_power_supply) +tx_connected, hv_connected = interface.is_device_connected() + +if not use_external_power_supply and not tx_connected: + logger.warning("TX device not connected. Attempting to turn on 12V...") + interface.hvcontroller.turn_12v_on() + + # Give time for the TX device to power up and enumerate over USB + time.sleep(2) + + # Cleanup and recreate interface to reinitialize USB devices + interface.stop_monitoring() + del interface + time.sleep(1) # Short delay before recreating + + logger.info("Reinitializing LIFU interface after powering 12V...") + interface = LIFUInterface(ext_power_supply=use_external_power_supply) + + # Re-check connection + tx_connected, hv_connected = interface.is_device_connected() + +if not use_external_power_supply: + if hv_connected: + logger.info(f" HV Connected: {hv_connected}") + else: + logger.error("❌ HV NOT fully connected.") + sys.exit(1) +else: + logger.info(" Using external power supply") + +if tx_connected: + logger.info(f" TX Connected: {tx_connected}") + logger.info("✅ LIFU Device fully connected.") +else: + logger.error("❌ TX NOT fully connected.") + sys.exit(1) + +# Verify communication with the devices +if not interface.txdevice.ping(): + logger.error("Failed to ping the transmitter device.") + sys.exit(1) + +if not use_external_power_supply and not interface.hvcontroller.ping(): + logger.error("Failed to ping the console device.") + sys.exit(1) + +print(f"console version: {interface.hvcontroller.get_version()}") + +logger.info("Enumerate TX7332 chips") +# num_tx_devices = interface.txdevice.get_tx_module_count() +num_tx_devices = 10 + +for module in range(num_tx_devices+1): + try: + tx_firmware_version = interface.txdevice.get_version(module=module) + logger.info(f"TX Firmware Version: {tx_firmware_version}") + except Exception as e: + logger.error(f"Error querying TX firmware version: {e}") diff --git a/pyproject.toml b/pyproject.toml index 4652de68..738dbf8d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,74 +29,72 @@ classifiers = [ ] dynamic = ["version"] dependencies = [ - "xarray[io]", - "numpy<2", - "pandas", - "scipy", - "requests", + "xarray[io]>=2026.2.0", + "numpy>=2.0.0", + "pandas>=3.0.2", + "scipy>=1.14.1", + "requests>=2.33.1", ] [project.optional-dependencies] jupyter = [ - "ipykernel", - "matplotlib", + "ipykernel>=7.2.0", + "matplotlib>=3.10.0", ] mesh = [ - "embreex; platform_machine=='x86_64' or platform_machine=='AMD64'", - "trimesh", - "scikit-image", - "vtk", - "Pillow", - "OpenEXR" + "embreex>=2.17.7.post7; platform_machine=='x86_64' or platform_machine=='AMD64'", + "trimesh>=4.11.5", + "scikit-image>=0.26.0", + "vtk>=9.6.1", + "Pillow>=12.2.0", + "OpenEXR>=3.4.9" ] db = [ - "nibabel", - "pydicom" + "nibabel>=5.4.2", + "pydicom>=3.0.2" ] io = [ - "base58", - "crc", - "crcmod", - "pyserial", - "openlifu-sdk>=2.0.12" + "crc>=7.1.0", + "crcmod>=1.7", + "openlifu-sdk>=2.0.14" ] sim = [ "k-wave-python==0.4.0", - "nvidia-ml-py" + "nvidia-ml-py>=13.595.45" ] cloud = [ - "watchdog", - "python-socketio[client]" + "watchdog>=6.0.0", + "python-socketio[client]>=5.16.1" ] photogrammetry = [ - "embreex; platform_machine=='x86_64' or platform_machine=='AMD64'", - "opencv-contrib-python", - "onnxruntime==1.18.0", - "trimesh", - "scikit-image", - "vtk", - "Pillow", - "OpenEXR" + "embreex>=2.17.7.post7; platform_machine=='x86_64' or platform_machine=='AMD64'", + "opencv-contrib-python>=4.11.0.86", + "onnxruntime>=1.20.0", + "trimesh>=4.11.5", + "scikit-image>=0.26.0", + "vtk>=9.6.1", + "Pillow>=12.2.0", + "OpenEXR>=3.4.9" ] dev = [ "pytest >=6", "pytest-cov >=3", - "pytest-mock", - "dvc[gdrive]", + "pytest-mock>=3.15.1", + "dvc[gdrive]>=3.67.1", ] docs = [ "openlifu[mesh, db, cloud]", "sphinx>=7.0", "myst_parser>=0.13", - "sphinx_copybutton", - "sphinx_autodoc_typehints", + "sphinx_copybutton>=0.5.2", + "sphinx_autodoc_typehints>=3.9.11", "furo>=2023.08.17", ] test = [ "openlifu[mesh, db, io, sim, cloud, photogrammetry]", "pytest >=6", "pytest-cov >=3", - "pytest-mock", + "pytest-mock>=3.15.1", ] all = [ "openlifu[jupyter, mesh, db, io, sim, cloud, photogrammetry, test, dev, docs]", diff --git a/src/openlifu/bf/apod_methods/maxangle.py b/src/openlifu/bf/apod_methods/maxangle.py index e30065df..368bfb29 100644 --- a/src/openlifu/bf/apod_methods/maxangle.py +++ b/src/openlifu/bf/apod_methods/maxangle.py @@ -10,16 +10,25 @@ from openlifu.bf.apod_methods import ApodizationMethod from openlifu.geo.point import Point from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.field_display import summarize_fields from openlifu.util.units import getunittype from openlifu.xdc import Transducer @dataclass class MaxAngle(ApodizationMethod): - max_angle: Annotated[float, OpenLIFUFieldData("Maximum acceptance angle", "Maximum acceptance angle for each element from the vector normal to the element surface")] = 30.0 + max_angle: Annotated[float, OpenLIFUFieldData( + name="Maximum acceptance angle", + description="Maximum acceptance angle for each element from the vector normal to the element surface", + units_field="units", display_units="deg", precision=1, + )] = 30.0 """Maximum acceptance angle for each element from the vector normal to the element surface""" - units: Annotated[str, OpenLIFUFieldData("Angle units", "Angle units")] = "deg" + units: Annotated[str, OpenLIFUFieldData( + name="Angle units", + description="Angle units", + unit_options=("deg", "rad"), + )] = "deg" """Angle units""" def __post_init__(self): @@ -47,3 +56,7 @@ def to_table(self) -> pd.DataFrame: records = [{"Name": "Type", "Value": "Max Angle", "Unit": ""}, {"Name": "Max Angle", "Value": self.max_angle, "Unit": self.units}] return pd.DataFrame.from_records(records) + + def get_summary(self) -> str: + """Return a one-liner summary of the apodization parameters.""" + return summarize_fields(self, ("max_angle",)) diff --git a/src/openlifu/bf/apod_methods/piecewiselinear.py b/src/openlifu/bf/apod_methods/piecewiselinear.py index 33179cba..c88e02ab 100644 --- a/src/openlifu/bf/apod_methods/piecewiselinear.py +++ b/src/openlifu/bf/apod_methods/piecewiselinear.py @@ -10,19 +10,32 @@ from openlifu.bf.apod_methods import ApodizationMethod from openlifu.geo.point import Point from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.field_display import summarize_fields from openlifu.util.units import getunittype from openlifu.xdc import Transducer @dataclass class PiecewiseLinear(ApodizationMethod): - zero_angle: Annotated[float, OpenLIFUFieldData("Zero Apodization Angle", "Angle at and beyond which the piecewise linear apodization is 0%")] = 90.0 + zero_angle: Annotated[float, OpenLIFUFieldData( + name="Zero apodization angle", + description="Angle at and beyond which the piecewise linear apodization is 0%", + units_field="units", display_units="deg", precision=1, + )] = 90.0 """Angle at and beyond which the piecewise linear apodization is 0%""" - rolloff_angle: Annotated[float, OpenLIFUFieldData("Rolloff start angle", "Angle below which the piecewise linear apodization is 100%")] = 45.0 + rolloff_angle: Annotated[float, OpenLIFUFieldData( + name="Rolloff start angle", + description="Angle below which the piecewise linear apodization is 100%", + units_field="units", display_units="deg", precision=1, + )] = 45.0 """Angle below which the piecewise linear apodization is 100%""" - units: Annotated[str, OpenLIFUFieldData("Angle units", "Angle units")] = "deg" + units: Annotated[str, OpenLIFUFieldData( + name="Angle units", + description="Angle units", + unit_options=("deg", "rad"), + )] = "deg" """Angle units""" def __post_init__(self): @@ -60,3 +73,7 @@ def to_table(self) -> pd.DataFrame: {"Name": "Rolloff Angle", "Value": self.rolloff_angle, "Unit": self.units}, ] return pd.DataFrame.from_records(records) + + def get_summary(self) -> str: + """Return a one-liner summary of the apodization parameters.""" + return summarize_fields(self, ("zero_angle", "rolloff_angle")) diff --git a/src/openlifu/bf/apod_methods/uniform.py b/src/openlifu/bf/apod_methods/uniform.py index fa2c28f3..a294db6b 100644 --- a/src/openlifu/bf/apod_methods/uniform.py +++ b/src/openlifu/bf/apod_methods/uniform.py @@ -10,12 +10,17 @@ from openlifu.bf.apod_methods import ApodizationMethod from openlifu.geo.point import Point from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.field_display import summarize_fields from openlifu.xdc import Transducer @dataclass class Uniform(ApodizationMethod): - value: Annotated[float, OpenLIFUFieldData("Value", "Uniform apodization value between 0 and 1.")] = 1.0 + value: Annotated[float, OpenLIFUFieldData( + name="Value", + description="Uniform apodization value between 0 and 1.", + precision=2, + )] = 1.0 """Uniform apodization value between 0 and 1.""" def calc_apodization(self, arr: Transducer, target: Point, params: xa.Dataset, transform:np.ndarray | None=None): @@ -30,3 +35,7 @@ def to_table(self): records = [{"Name": "Type", "Value": "Uniform", "Unit": ""}, {"Name": "Value", "Value": self.value, "Unit": ""}] return pd.DataFrame.from_records(records) + + def get_summary(self) -> str: + """Return a one-liner summary of the apodization parameters.""" + return summarize_fields(self, ("value",)) diff --git a/src/openlifu/bf/delay_methods/direct.py b/src/openlifu/bf/delay_methods/direct.py index c66a7c73..4d67e82e 100644 --- a/src/openlifu/bf/delay_methods/direct.py +++ b/src/openlifu/bf/delay_methods/direct.py @@ -10,12 +10,17 @@ from openlifu.bf.delay_methods import DelayMethod from openlifu.geo.point import Point from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.field_display import summarize_fields from openlifu.xdc import Transducer @dataclass class Direct(DelayMethod): - c0: Annotated[float, OpenLIFUFieldData("Speed of Sound (m/s)", "Speed of sound in the medium (m/s)")] = 1480.0 + c0: Annotated[float, OpenLIFUFieldData( + name="Speed of sound", + description="Speed of sound in the medium", + units="m/s", precision=0, + )] = 1480.0 """Speed of sound in the medium (m/s)""" def __post_init__(self): @@ -46,3 +51,7 @@ def to_table(self) -> pd.DataFrame: records = [{"Name": "Type", "Value": "Direct", "Unit": ""}, {"Name": "Default Sound Speed", "Value": self.c0, "Unit": "m/s"}] return pd.DataFrame.from_records(records) + + def get_summary(self) -> str: + """Return a one-liner summary of the delay method.""" + return summarize_fields(self, ("c0",)) diff --git a/src/openlifu/bf/focal_patterns/focal_pattern.py b/src/openlifu/bf/focal_patterns/focal_pattern.py index f12af4f4..b67c101c 100644 --- a/src/openlifu/bf/focal_patterns/focal_pattern.py +++ b/src/openlifu/bf/focal_patterns/focal_pattern.py @@ -9,6 +9,7 @@ from openlifu.bf import focal_patterns from openlifu.geo.point import Point from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.field_display import summarize_fields from openlifu.util.units import getunittype @@ -18,10 +19,18 @@ class FocalPattern(ABC): Abstract base class for representing a focal pattern """ - target_pressure: Annotated[float, OpenLIFUFieldData("Target pressure", "Target pressure of the focal pattern in given units")] = 1.0 + target_pressure: Annotated[float, OpenLIFUFieldData( + name="Target pressure", + description="Target pressure of the focal pattern", + units_field="units", display_units="kPa", precision=0, + )] = 1.0 """Target pressure of the focal pattern in given units""" - units: Annotated[str, OpenLIFUFieldData("Pressure units", "Pressure units")] = "Pa" + units: Annotated[str, OpenLIFUFieldData( + name="Pressure units", + description="Pressure units (Pa, kPa, MPa)", + unit_options=("Pa", "kPa", "MPa"), + )] = "Pa" """Pressure units""" def __post_init__(self): @@ -51,6 +60,15 @@ def num_foci(self): """ pass + @abstractmethod + def get_order(self): + """ + Get the order of foci in the focal pattern + + :returns: List of indices of foci in the order they are used in the pulse sequence + """ + pass + def to_dict(self): """ Convert the focal pattern to a dictionary @@ -83,3 +101,11 @@ def to_table(self) -> pd.DataFrame: :returns: Pandas DataFrame of the focal pattern parameters """ pass + + def get_summary(self) -> str: + """Return a one-liner summary of this focal pattern's parameters. + + Subclasses may override to include their own additional fields; by + default the base class summarizes ``target_pressure``. + """ + return summarize_fields(self, ("target_pressure",)) diff --git a/src/openlifu/bf/focal_patterns/single.py b/src/openlifu/bf/focal_patterns/single.py index bdfeea7d..2be7bbef 100644 --- a/src/openlifu/bf/focal_patterns/single.py +++ b/src/openlifu/bf/focal_patterns/single.py @@ -32,6 +32,14 @@ def num_foci(self): """ return 1 + def get_order(self): + """ + Get the order of foci in the focal pattern + + :returns: List of indices of foci in the order they are used in the pulse sequence + """ + return [1] + def to_table(self) -> pd.DataFrame: """ Get a table of the focal pattern parameters diff --git a/src/openlifu/bf/focal_patterns/wheel.py b/src/openlifu/bf/focal_patterns/wheel.py index c768ceeb..90599616 100644 --- a/src/openlifu/bf/focal_patterns/wheel.py +++ b/src/openlifu/bf/focal_patterns/wheel.py @@ -1,7 +1,7 @@ from __future__ import annotations from dataclasses import dataclass -from typing import Annotated +from typing import Annotated, List import numpy as np import pandas as pd @@ -9,6 +9,7 @@ from openlifu.bf.focal_patterns import FocalPattern from openlifu.geo.point import Point from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.field_display import summarize_fields @dataclass @@ -17,25 +18,66 @@ class Wheel(FocalPattern): Class for representing a wheel pattern """ - center: Annotated[bool, OpenLIFUFieldData("Include center point?", "Whether to include the center for the wheel pattern")] = True + center: Annotated[bool, OpenLIFUFieldData( + name="Include center point?", + description="Whether to include the center for the wheel pattern", + )] = True """Whether to include the center for the wheel pattern""" - num_spokes: Annotated[int, OpenLIFUFieldData("Number of spokes", "Number of spokes in the wheel pattern")] = 4 + num_spokes: Annotated[int | List[int], OpenLIFUFieldData( + name="Number of spokes", + description="Number of spokes in the wheel pattern", + )] = 4 """Number of spokes in the wheel pattern""" - spoke_radius: Annotated[float, OpenLIFUFieldData("Spoke radius", "Radius of the spokes in the wheel pattern")] = 1.0 # mm + spoke_radius: Annotated[float | List[float], OpenLIFUFieldData( + name="Spoke radius", + description="Radius of the spokes in the wheel pattern", + units_field="distance_units", display_units="mm", precision=1, + )] = 1.0 # mm """Radius of the spokes in the wheel pattern""" - distance_units: Annotated[str, OpenLIFUFieldData("Units", "Units of the wheel pattern parameters")] = "mm" + spoke_phase: Annotated[float | List[float], OpenLIFUFieldData( + name="Spoke phase", + description="Phase of the spokes in the wheel pattern", + units_field="angle_units", display_units="deg", precision=1, + )] = 0.0 # degrees + """Phase of the spokes in the wheel pattern""" + + distance_units: Annotated[str, OpenLIFUFieldData( + name="Distance units", + description="Units of the wheel pattern parameters", + unit_options=("mm", "cm", "m"), + )] = "mm" """Units of the wheel pattern parameters""" + order: Annotated[list[int] | None, OpenLIFUFieldData("Focus order", "Order of Foci (1-indexed) in the sequence")] = None + """Order of Foci (1-indexed) in the sequence. This is a list of integers that specifies the order in which the foci are used in the pulse sequence. If None, the foci are used in the order they are listed in the `foci` attribute.""" + def __post_init__(self): if not isinstance(self.center, bool): raise TypeError(f"Center must be a boolean, got {type(self.center).__name__}.") - if not isinstance(self.num_spokes, int) or self.num_spokes < 1: - raise ValueError(f"Number of spokes must be a positive integer, got {self.num_spokes}.") - if not isinstance(self.spoke_radius, int | float) or self.spoke_radius <= 0: - raise ValueError(f"Spoke radius must be a positive number, got {self.spoke_radius}.") + if isinstance(self.num_spokes, int): + if self.num_spokes < 1: + raise ValueError(f"Number of spokes must be a positive integer, got {self.num_spokes}.") + elif isinstance(self.num_spokes, list): + if not all(isinstance(n, int) and n > 0 for n in self.num_spokes): + raise ValueError(f"All elements of num_spokes must be positive integers, got {self.num_spokes}.") + else: + raise TypeError(f"num_spokes must be an int or list of ints, got {type(self.num_spokes).__name__}.") + + if isinstance(self.spoke_radius, int | float): + if self.spoke_radius <= 0: + raise ValueError(f"Spoke radius must be a positive number, got {self.spoke_radius}.") + elif isinstance(self.spoke_radius, list): + if not all(isinstance(r, int | float) and r > 0 for r in self.spoke_radius): + raise ValueError(f"All elements of spoke_radius must be positive numbers, got {self.spoke_radius}.") + if isinstance(self.spoke_phase, list) and len(self.spoke_phase) != len(self.spoke_radius): + raise ValueError(f"Length of spoke_phase list must match length of spoke_radius list, got {len(self.spoke_phase)} and {len(self.spoke_radius)}.") + if isinstance(self.num_spokes, list) and len(self.spoke_radius) != len(self.num_spokes): + raise ValueError(f"Length of spoke_radius list must match length of num_spokes list, got {len(self.spoke_radius)} and {len(self.num_spokes)}.") + else: + raise TypeError(f"spoke_radius must be a number or list of numbers, got {type(self.spoke_radius).__name__}.") super().__post_init__() def get_targets(self, target: Point): @@ -45,23 +87,41 @@ def get_targets(self, target: Point): :param target: Target point of the focal pattern :returns: List of target points """ + m = target.get_matrix(center_on_point=True) + if isinstance(self.num_spokes, int): + if isinstance(self.spoke_radius, int | float): + spoke_radius_list = [self.spoke_radius] + spoke_phase_list = [self.spoke_phase] + else: + spoke_radius_list = self.spoke_radius + spoke_phase_list = self.spoke_phase if isinstance(self.spoke_phase, list) else [self.spoke_phase]*len(spoke_radius_list) + num_spokes_list = [self.num_spokes]*len(spoke_radius_list) + else: + num_spokes_list = self.num_spokes + spoke_radius_list = self.spoke_radius + spoke_phase_list = self.spoke_phase if isinstance(self.spoke_phase, list) else [self.spoke_phase]*len(spoke_radius_list) + + n_points = sum(num_spokes_list) + int(self.center) + if self.center: targets = [target.copy()] - targets[0].id = f"{target.id}_center" - targets[0].id = f"{target.id} (Center)" + targets[0].id = f"{target.id}_01" + targets[0].id = f"{target.id} (1/{n_points}, Center)" else: targets = [] - m = target.get_matrix(center_on_point=True) - for i in range(self.num_spokes): - theta = 2*np.pi*i/self.num_spokes - local_position = self.spoke_radius * np.array([np.cos(theta), np.sin(theta), 0.0]) - position = np.dot(m, np.append(local_position, 1.0))[:3] - spoke = Point(id=f"{target.id}_{np.rad2deg(theta):.0f}deg", - name=f"{target.name} ({np.rad2deg(theta):.0f}°)", + + for (num_spokes, spoke_radius, spoke_phase) in zip(num_spokes_list, spoke_radius_list, spoke_phase_list): + for j in range(num_spokes): + point_index = len(targets) + 1 + theta = 2*np.pi*j/num_spokes + np.deg2rad(spoke_phase) + local_position = spoke_radius * np.array([np.cos(theta), np.sin(theta), 0.0]) + position = np.dot(m, np.append(local_position, 1.0))[:3] + spoke = Point(id=f"{target.id}_{point_index:02d}", + name=f"{target.name} ({point_index}/{n_points}, {spoke_radius:.1f} mm, {np.rad2deg(theta):.0f}°)", position=position, units=self.distance_units, radius=target.radius) - targets.append(spoke) + targets.append(spoke) return targets def num_foci(self) -> int: @@ -70,7 +130,21 @@ def num_foci(self) -> int: :returns: Number of foci """ - return int(self.center) + self.num_spokes + if isinstance(self.num_spokes, int): + return int(self.center) + self.num_spokes + else: + return int(self.center) + sum(self.num_spokes) + + def get_order(self): + """ + Get the order of foci in the focal pattern + + :returns: List of indices of foci in the order they are used in the pulse sequence + """ + if self.order is not None: + return self.order + else: + return list(range(1, self.num_foci() + 1)) def to_table(self) -> pd.DataFrame: """ @@ -86,3 +160,10 @@ def to_table(self) -> pd.DataFrame: {"Name": "Spoke Radius", "Value": self.spoke_radius, "Unit": self.distance_units}, ] return pd.DataFrame.from_records(records) + + def get_summary(self) -> str: + """Return a one-liner summary including base + wheel-specific fields.""" + return summarize_fields( + self, + ("target_pressure", "num_spokes", "spoke_radius"), + ) diff --git a/src/openlifu/bf/pulse.py b/src/openlifu/bf/pulse.py index de40135a..e4c47a10 100644 --- a/src/openlifu/bf/pulse.py +++ b/src/openlifu/bf/pulse.py @@ -8,6 +8,7 @@ from openlifu.util.annotations import OpenLIFUFieldData from openlifu.util.dict_conversion import DictMixin +from openlifu.util.field_display import summarize_fields @dataclass @@ -16,13 +17,25 @@ class Pulse(DictMixin): Class for representing a sinusoidal pulse """ - frequency: Annotated[float, OpenLIFUFieldData("Frequency (Hz)", "Frequency of the pulse in Hz")] = 1.0 # Hz + frequency: Annotated[float, OpenLIFUFieldData( + name="Frequency", + description="Frequency of the pulse", + units="Hz", display_units="kHz", precision=1, + )] = 1.0 # Hz """Frequency of the pulse in Hz""" - amplitude: Annotated[float, OpenLIFUFieldData("Amplitude (AU)", "Amplitude of the pulse (between 0 and 1). ")] = 1.0 # AU + amplitude: Annotated[float, OpenLIFUFieldData( + name="Amplitude", + description="Amplitude of the pulse (between 0 and 1).", + precision=2, + )] = 1.0 # AU """Amplitude of the pulse in arbitrary units (AU) between 0 and 1""" - duration: Annotated[float, OpenLIFUFieldData("Duration (s)", "Duration of the pulse in s")] = 1.0 # s + duration: Annotated[float, OpenLIFUFieldData( + name="Duration", + description="Duration of the pulse", + units="s", display_units="ms", precision=1, + )] = 1.0 # s """Duration of the pulse in s""" def __post_init__(self): @@ -61,3 +74,7 @@ def to_table(self) -> pd.DataFrame: {"Name": "Amplitude", "Value": self.amplitude, "Unit": "AU"}, {"Name": "Duration", "Value": self.duration, "Unit": "s"}] return pd.DataFrame.from_records(records) + + def get_summary(self) -> str: + """Return a one-liner summary like ``Frequency: 400 kHz, Amplitude: 1, Duration: 5 ms``.""" + return summarize_fields(self, ("frequency", "amplitude", "duration")) diff --git a/src/openlifu/bf/sequence.py b/src/openlifu/bf/sequence.py index a5601075..b86c4c64 100644 --- a/src/openlifu/bf/sequence.py +++ b/src/openlifu/bf/sequence.py @@ -15,16 +15,30 @@ class Sequence(DictMixin): Class for representing a sequence of pulses """ - pulse_interval: Annotated[float, OpenLIFUFieldData("Pulse interval (s)", "Interval between pulses in the sequence (s)")] = 1.0 # s + pulse_interval: Annotated[float, OpenLIFUFieldData( + name="Pulse interval", + description="Interval between pulses in the sequence", + units="s", display_units="ms", precision=1, + )] = 1.0 # s """Interval between pulses in the sequence (s)""" - pulse_count: Annotated[int, OpenLIFUFieldData("Pulse count", "Number of pulses in the sequence")] = 1 + pulse_count: Annotated[int, OpenLIFUFieldData( + name="Pulse count", + description="Number of pulses in the sequence", + )] = 1 """Number of pulses in the sequence""" - pulse_train_interval: Annotated[float, OpenLIFUFieldData("Pulse train interval (s)", "Interval between pulse trains in the sequence (s)")] = 1.0 # s + pulse_train_interval: Annotated[float, OpenLIFUFieldData( + name="Pulse train interval", + description="Interval between pulse trains in the sequence", + units="s", display_units="s", precision=2, + )] = 1.0 # s """Interval between pulse trains in the sequence (s)""" - pulse_train_count: Annotated[int, OpenLIFUFieldData("Pulse train count", "Number of pulse trains in the sequence")] = 1 + pulse_train_count: Annotated[int, OpenLIFUFieldData( + name="Pulse train count", + description="Number of pulse trains in the sequence", + )] = 1 """Number of pulse trains in the sequence""" def __post_init__(self): @@ -72,3 +86,17 @@ def get_sequence_duration(self) -> float: else: interval = self.pulse_train_interval return interval * self.pulse_train_count + + def get_summary(self) -> str: + """Return a one-liner summary of the sequence parameters. + + Format: ``"{pulse_count} pulses every {pulse_interval}ms, + repeated {pulse_train_count}x every {pulse_train_interval}s"``. + Numeric values use ``%g`` formatting (no trailing zeros). + """ + pulse_interval_ms = self.pulse_interval * 1000.0 + pulse_train_interval_s = self.pulse_train_interval + return ( + f"{int(self.pulse_count)} pulses every {pulse_interval_ms:g}ms, " + f"repeated {int(self.pulse_train_count)}x every {pulse_train_interval_s:g}s" + ) diff --git a/src/openlifu/cloud/api/api.py b/src/openlifu/cloud/api/api.py index 0948ad8c..956251f6 100644 --- a/src/openlifu/cloud/api/api.py +++ b/src/openlifu/cloud/api/api.py @@ -16,8 +16,8 @@ class Api: - def __init__(self): - self._request = Request() + def __init__(self, api_url: str): + self._request = Request(api_url) self._request.debug_log = True self._databases = DatabasesApi(self._request) self._protocols = ProtocolsApi(self._request) diff --git a/src/openlifu/cloud/api/request.py b/src/openlifu/cloud/api/request.py index 8ed8520e..3c55811a 100644 --- a/src/openlifu/cloud/api/request.py +++ b/src/openlifu/cloud/api/request.py @@ -7,7 +7,6 @@ import urllib3 from requests.adapters import HTTPAdapter -from openlifu.cloud.const import API_URL from openlifu.cloud.utils import logger_cloud, to_json urllib3.disable_warnings(urllib3.exceptions.InsecureRequestWarning) @@ -34,7 +33,8 @@ def init_poolmanager(self, *args, **kwargs): class Request: TIMEOUT = (5, 300) - def __init__(self): + def __init__(self, api_url: str): + self._api_url = api_url self.headers = {} self.session = requests.Session() adapter = SlicerAdapter( @@ -48,7 +48,7 @@ def _log_request(self, method: str, url: str, start_time: float, status_code: in def get(self, url: str) -> str: start = time.perf_counter() - response = self.session.get(API_URL + url, headers=self.headers, timeout=self.TIMEOUT, verify=False) + response = self.session.get(self._api_url + url, headers=self.headers, timeout=self.TIMEOUT, verify=False) self._log_request("GET", url, start, response.status_code) logger_cloud.debug(f"GET: {url}, status_code: {response.status_code}\nresponse: {response.text}") @@ -57,7 +57,7 @@ def get(self, url: str) -> str: def get_bytes(self, url: str) -> bytes: start = time.perf_counter() - response = self.session.get(API_URL + url, headers=self.headers, timeout=self.TIMEOUT, verify=False) + response = self.session.get(self._api_url + url, headers=self.headers, timeout=self.TIMEOUT, verify=False) self._log_request("GET_BYTES", url, start, response.status_code) logger_cloud.debug(f"GET bytes: {url}, status_code: {response.status_code}") @@ -66,7 +66,7 @@ def get_bytes(self, url: str) -> bytes: def post(self, url: str, dto) -> str: start = time.perf_counter() - response = self.session.post(API_URL + url, data=to_json(dto), headers=self.headers, timeout=self.TIMEOUT, verify=False) + response = self.session.post(self._api_url + url, data=to_json(dto), headers=self.headers, timeout=self.TIMEOUT, verify=False) self._log_request("POST", url, start, response.status_code) logger_cloud.debug(f"POST: {url}, body: {to_json(dto)}, status_code: {response.status_code}\nresponse: {response.text}") @@ -75,7 +75,7 @@ def post(self, url: str, dto) -> str: def post_bytes(self, url: str, data) -> str: start = time.perf_counter() - response = self.session.post(API_URL + url, data=data, headers=self.headers, timeout=self.TIMEOUT, verify=False) + response = self.session.post(self._api_url + url, data=data, headers=self.headers, timeout=self.TIMEOUT, verify=False) self._log_request("POST_BYTES", url, start, response.status_code) logger_cloud.debug(f"POST bytes: {url}, status_code: {response.status_code}\nresponse: {response.text}") @@ -84,7 +84,7 @@ def post_bytes(self, url: str, data) -> str: def put(self, url: str, dto) -> str: start = time.perf_counter() - response = self.session.put(API_URL + url, data=to_json(dto), headers=self.headers, timeout=self.TIMEOUT, verify=False) + response = self.session.put(self._api_url + url, data=to_json(dto), headers=self.headers, timeout=self.TIMEOUT, verify=False) self._log_request("PUT", url, start, response.status_code) logger_cloud.debug(f"PUT: {url}, body: {to_json(dto)}, status_code: {response.status_code}\nresponse: {response.text}") @@ -93,7 +93,7 @@ def put(self, url: str, dto) -> str: def delete(self, url: str) -> str: start = time.perf_counter() - response = self.session.delete(API_URL + url, headers=self.headers, timeout=self.TIMEOUT, verify=False) + response = self.session.delete(self._api_url + url, headers=self.headers, timeout=self.TIMEOUT, verify=False) self._log_request("DELETE", url, start, response.status_code) logger_cloud.debug(f"DELETE: {url}, status_code: {response.status_code}\nresponse: {response.text}") diff --git a/src/openlifu/cloud/cloud.py b/src/openlifu/cloud/cloud.py index 5ccddcae..1b740eac 100644 --- a/src/openlifu/cloud/cloud.py +++ b/src/openlifu/cloud/cloud.py @@ -17,10 +17,10 @@ from openlifu.cloud.components.sessions import Sessions from openlifu.cloud.components.solutions import Solutions from openlifu.cloud.components.subjects import Subjects -from openlifu.cloud.components.systems import Systems from openlifu.cloud.components.transducers import Transducers from openlifu.cloud.components.users import Users from openlifu.cloud.components.volumes import Volumes +from openlifu.cloud.const import API_URL_DEV, API_URL_PROD, ENV_DEV, ENV_PROD from openlifu.cloud.filesystem_observer import FilesystemObserver from openlifu.cloud.status import Status from openlifu.cloud.sync_thread import SyncThread @@ -30,10 +30,14 @@ class Cloud: - def __init__(self): + def __init__(self, environment: str = ENV_PROD): + if environment == ENV_DEV: + api_url = API_URL_DEV + else: + api_url = API_URL_PROD self._filesystem_observer = FilesystemObserver(self._on_file_system_update) - self._api = Api() - self._websocket = Websocket(self._on_websocket_update) + self._api = Api(api_url) + self._websocket = Websocket(api_url, self._on_websocket_update) self._components: List[AbstractComponent] = [] self._sync_thread = SyncThread(self._on_status_changed) self._db_path: Path | None = None @@ -146,7 +150,9 @@ def _create_components(self): self._components.clear() self._components.append(Users(self._api, self._db_path, self._db.id, self._sync_thread)) self._components.append(Protocols(self._api, self._db_path, self._db.id, self._sync_thread)) - self._components.append(Systems(self._api, self._db_path, self._db.id, self._sync_thread)) + # Systems are managed elsewhere and the per-database "systems" + # folder is no longer part of the local database layout, so the + # Systems sync component is intentionally not registered here. self._components.append(Transducers(self._api, self._db_path, self._db.id, self._sync_thread)) self._components.append( Subjects(self._api, self._db_path, self._db.id, self._sync_thread) @@ -172,7 +178,7 @@ def _create_components(self): logger_cloud.setLevel(logging.DEBUG) logger_cloud.addHandler(logging.StreamHandler(sys.stdout)) - cloud = Cloud() + cloud = Cloud(ENV_DEV) token = os.getenv("TOKEN") db_path = os.getenv("DB_PATH") diff --git a/src/openlifu/cloud/const.py b/src/openlifu/cloud/const.py index 731e7af6..8089fcbb 100644 --- a/src/openlifu/cloud/const.py +++ b/src/openlifu/cloud/const.py @@ -1,6 +1,10 @@ from __future__ import annotations -API_URL = "https://api.openwater.health" +API_URL_PROD = "https://api.openwater.health" +API_URL_DEV = "https://dev.api.openwater.health" + +ENV_PROD = "prod" +ENV_DEV = "dev" CONFIG_FILE = "config" DATA_FILE = "data" diff --git a/src/openlifu/cloud/ws.py b/src/openlifu/cloud/ws.py index cd41f05b..d3195c82 100644 --- a/src/openlifu/cloud/ws.py +++ b/src/openlifu/cloud/ws.py @@ -5,14 +5,14 @@ import socketio from socketio import exceptions -from openlifu.cloud.const import API_URL from openlifu.cloud.utils import logger_cloud DATABASE_UPDATES_NS = "/database_updates" class Websocket: - def __init__(self, update_callback: Callable[[dict], None]): + def __init__(self, api_url: str, update_callback: Callable[[dict], None]): + self._api_url = api_url self._sio: socketio.Client | None = None self._database_id = None self._auth = {} @@ -28,7 +28,7 @@ def authenticate(self, access_token: str): self.connect(self._database_id) def connect(self, database_id: int): - self.log(f"Attempting connection to {API_URL} for DB {database_id}") + self.log(f"Attempting connection to {self._api_url} for DB {database_id}") if self._sio is not None: self.disconnect() @@ -71,7 +71,7 @@ def on_update(data): try: self._sio.connect( - f"{API_URL}/socket.io", + f"{self._api_url}/socket.io", auth=self._auth, namespaces=[DATABASE_UPDATES_NS], transports=["websocket"], diff --git a/src/openlifu/db/__init__.py b/src/openlifu/db/__init__.py index f563a8c1..133eb9e4 100644 --- a/src/openlifu/db/__init__.py +++ b/src/openlifu/db/__init__.py @@ -1,13 +1,19 @@ from __future__ import annotations from openlifu.db.database import Database +from openlifu.db.plan import Plan +from openlifu.db.planning_session import PlanningSession from openlifu.db.session import Session +from openlifu.db.sonication_session import SonicationSession from openlifu.db.subject import Subject from openlifu.db.user import User __all__ = [ - "Subject", - "Session", "Database", + "Plan", + "PlanningSession", + "Session", + "SonicationSession", + "Subject", "User", ] diff --git a/src/openlifu/db/database.py b/src/openlifu/db/database.py index f82dfa48..b88be62f 100644 --- a/src/openlifu/db/database.py +++ b/src/openlifu/db/database.py @@ -10,7 +10,7 @@ from typing import Dict, List from openlifu.nav.photoscan import Photoscan, load_data_from_photoscan -from openlifu.plan import Protocol, Run, Solution +from openlifu.plan import Protocol, Run, Solution, SolutionAnalysis from openlifu.util.json import PYFUSEncoder from openlifu.util.types import PathLike from openlifu.util.volume_conversion import ( @@ -20,7 +20,10 @@ from openlifu.xdc import Transducer, TransducerArray from openlifu.xdc.util import load_transducer_from_file +from .plan import Plan +from .planning_session import PlanningSession from .session import Session +from .sonication_session import SonicationSession from .subject import Subject from .user import User @@ -568,6 +571,750 @@ def write_solution(self, session:Session, solution:Solution, on_conflict: OnConf self.logger.info(f"Wrote solution with ID {solution.id} to the database.") + def write_solution_analysis( + self, + session: Session, + solution_id: str, + analysis: SolutionAnalysis, + on_conflict: OnConflictOpts | str = OnConflictOpts.OVERWRITE, + ) -> None: + """Write a SolutionAnalysis next to its parent Solution. + + The analysis file lives at ``/_analysis.json``. Defaults to overwriting + because the analysis is a derived artifact: re-running ``Solution.analyze`` should always be safe to + re-persist over a stale copy. + """ + on_conflict = _normalize_on_conflict(on_conflict) + analysis_filepath = self.get_solution_analysis_filepath(session.subject_id, session.id, solution_id) + if analysis_filepath.exists(): + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"SolutionAnalysis for solution {solution_id} already exists in the database." + ) + if on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping SolutionAnalysis for solution {solution_id} as it already exists." + ) + return + if on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting SolutionAnalysis for solution {solution_id} in the database." + ) + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + analysis_filepath.parent.mkdir(parents=True, exist_ok=True) + analysis_filepath.write_text(analysis.to_json(compact=False)) + self.logger.info( + f"Wrote SolutionAnalysis for solution {solution_id} to the database." + ) + + def load_solution_analysis(self, session: Session, solution_id: str) -> SolutionAnalysis: + """Load the SolutionAnalysis associated with the given Solution within the given Session.""" + analysis_filepath = self.get_solution_analysis_filepath(session.subject_id, session.id, solution_id) + if not analysis_filepath.exists() or not analysis_filepath.is_file(): + self.logger.error( + f"SolutionAnalysis file not found for solution {solution_id}, session {session.id}" + ) + raise FileNotFoundError( + f"SolutionAnalysis file not found for solution {solution_id}, session {session.id}" + ) + analysis = SolutionAnalysis.from_json(analysis_filepath.read_text()) + self.logger.info(f"Loaded SolutionAnalysis for solution {solution_id}") + return analysis + + def delete_solution( + self, + session: Session, + solution_id: str, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Delete a solution and its associated files (analysis, .nc) from the database. + + Removes the ``{session_dir}/solutions/{solution_id}/`` directory and drops + ``solution_id`` from the session's ``solutions.json`` index. + + Args: + session: Session that owns the solution. + solution_id: ID of the solution to delete. + on_conflict: Behavior when the solution does not exist ('error' or 'skip'). + """ + on_conflict = _normalize_on_conflict(on_conflict) + solution_ids = self.get_solution_ids(session.subject_id, session.id) + + if solution_id not in solution_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Solution ID {solution_id} does not exist for session {session.id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Cannot delete solution ID {solution_id} as it does not exist for session {session.id}." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error' or 'skip'.") + + solution_dir = self.get_solution_filepath(session.subject_id, session.id, solution_id).parent + if solution_dir.is_dir(): + shutil.rmtree(solution_dir) + + solution_ids.remove(solution_id) + self.write_solution_ids(session, solution_ids) + + self.logger.info(f"Removed solution with ID {solution_id} from the database.") + + def purge_orphaned_solutions(self, session: Session) -> List[str]: + """Delete any on-disk solutions for ``session`` that are not tracked by ``session.solutions``. + + The authoritative list of solutions belonging to a session is + ``session.solutions`` (a list of ``SolutionInfo``). This method reconciles the on-disk state + against that list: any ``{session_dir}/solutions/{sid}/`` whose ``sid`` is not in + ``{si.solution_id for si in session.solutions}`` is removed, and the ``solutions.json`` index + is trimmed accordingly. + + Returns: + List of solution IDs that were purged (may be empty). + """ + tracked_ids = {si.solution_id for si in session.solutions} + on_disk_ids = self.get_solution_ids(session.subject_id, session.id) + purged = [sid for sid in on_disk_ids if sid not in tracked_ids] + for sid in purged: + self.delete_solution(session, sid, on_conflict=OnConflictOpts.SKIP) + return purged + + # ------------------------------------------------------------------ + # Subject-scoped solution read/write/delete (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def write_solution_at_subject_scope( + self, + subject_id: str, + solution: Solution, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Write a Solution to the subject-scoped ``solutions/{sid}/`` directory. + + The solution JSON and ``.nc`` files land at + ``subjects/{subject_id}/solutions/{solution.id}/``, and the solution id is + added to the subject's ``solutions/solutions.json`` index. No Session is + required or referenced. + """ + on_conflict = _normalize_on_conflict(on_conflict) + solution_ids = self.get_subject_solution_ids(subject_id) + + if solution.id in solution_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Solution with ID {solution.id} already exists at subject scope for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting solution with ID {solution.id} at subject scope for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping solution with ID {solution.id} at subject scope as it already exists." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + + solution_json_filepath = self.get_subject_solution_filepath(subject_id, solution.id) + solution_json_filepath.parent.mkdir(parents=True, exist_ok=True) + solution.to_files(solution_json_filepath) + + if solution.id not in solution_ids: + solution_ids.append(solution.id) + self.write_subject_solution_ids(subject_id, solution_ids) + + self.logger.info( + f"Wrote solution {solution.id} at subject scope for subject {subject_id}." + ) + + def load_solution_at_subject_scope(self, subject_id: str, solution_id: str) -> Solution: + """Load a subject-scoped Solution by id. + + No Session is required or referenced. Raises ``FileNotFoundError`` if the solution + does not exist under ``subjects/{subject_id}/solutions/``. + """ + solution_json_filepath = self.get_subject_solution_filepath(subject_id, solution_id) + if not (solution_json_filepath.exists() and solution_json_filepath.is_file()): + self.logger.error( + f"Solution file not found at subject scope for solution {solution_id}, subject {subject_id}" + ) + raise FileNotFoundError( + f"Solution file not found at subject scope for solution {solution_id}, subject {subject_id}" + ) + solution = Solution.from_files(solution_json_filepath) + self.logger.info(f"Loaded solution {solution_id} at subject scope for subject {subject_id}.") + return solution + + def write_solution_analysis_at_subject_scope( + self, + subject_id: str, + solution_id: str, + analysis: SolutionAnalysis, + on_conflict: OnConflictOpts | str = OnConflictOpts.OVERWRITE, + ) -> None: + """Write a SolutionAnalysis next to its parent subject-scoped Solution.""" + on_conflict = _normalize_on_conflict(on_conflict) + analysis_filepath = self.get_subject_solution_analysis_filepath(subject_id, solution_id) + if analysis_filepath.exists(): + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"SolutionAnalysis for solution {solution_id} already exists at subject scope." + ) + if on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping SolutionAnalysis for solution {solution_id} as it already exists at subject scope." + ) + return + if on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting SolutionAnalysis for solution {solution_id} at subject scope." + ) + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + analysis_filepath.parent.mkdir(parents=True, exist_ok=True) + analysis_filepath.write_text(analysis.to_json(compact=False)) + self.logger.info( + f"Wrote SolutionAnalysis for solution {solution_id} at subject scope for subject {subject_id}." + ) + + def load_solution_analysis_at_subject_scope( + self, subject_id: str, solution_id: str + ) -> SolutionAnalysis: + """Load a SolutionAnalysis for a subject-scoped Solution.""" + analysis_filepath = self.get_subject_solution_analysis_filepath(subject_id, solution_id) + if not (analysis_filepath.exists() and analysis_filepath.is_file()): + self.logger.error( + f"SolutionAnalysis file not found at subject scope for solution {solution_id}, subject {subject_id}" + ) + raise FileNotFoundError( + f"SolutionAnalysis file not found at subject scope for solution {solution_id}, subject {subject_id}" + ) + analysis = SolutionAnalysis.from_json(analysis_filepath.read_text()) + self.logger.info( + f"Loaded SolutionAnalysis for solution {solution_id} at subject scope for subject {subject_id}." + ) + return analysis + + def delete_solution_at_subject_scope( + self, + subject_id: str, + solution_id: str, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Delete a subject-scoped Solution and its associated files. + + Removes ``subjects/{subject_id}/solutions/{solution_id}/`` and drops + ``solution_id`` from the subject's ``solutions.json`` index. Does not consult + or affect any Session record. + """ + on_conflict = _normalize_on_conflict(on_conflict) + solution_ids = self.get_subject_solution_ids(subject_id) + + if solution_id not in solution_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Solution ID {solution_id} does not exist at subject scope for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Cannot delete solution ID {solution_id} at subject scope as it does not exist for subject {subject_id}." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error' or 'skip'.") + + solution_dir = self.get_subject_solution_dir(subject_id, solution_id) + if solution_dir.is_dir(): + shutil.rmtree(solution_dir) + + solution_ids.remove(solution_id) + self.write_subject_solution_ids(subject_id, solution_ids) + + self.logger.info( + f"Removed solution {solution_id} at subject scope for subject {subject_id}." + ) + + # ------------------------------------------------------------------ + # Plan read/write/delete (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def write_plan( + self, + subject_id: str, + plan: Plan, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Write a Plan to the subject's ``plans/{plan_id}/`` directory.""" + on_conflict = _normalize_on_conflict(on_conflict) + if plan.subject_id is not None and plan.subject_id != subject_id: + raise ValueError( + f"Plan.subject_id ({plan.subject_id!r}) does not match subject_id argument ({subject_id!r})." + ) + plan_ids = self.get_plan_ids(subject_id) + if plan.id in plan_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Plan with ID {plan.id} already exists for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting plan {plan.id} for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping plan {plan.id} for subject {subject_id} as it already exists." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + + plan_filename = self.get_plan_filename(subject_id, plan.id) + plan_filename.parent.parent.mkdir(parents=True, exist_ok=True) + plan_filename.parent.mkdir(exist_ok=True) + plan.to_file(plan_filename) + + if plan.id not in plan_ids: + plan_ids.append(plan.id) + self.write_plan_ids(subject_id, plan_ids) + + self.logger.info(f"Wrote plan {plan.id} for subject {subject_id}.") + + def load_plan(self, subject_id: str, plan_id: str) -> Plan: + """Load a Plan by id from the subject's ``plans/`` directory.""" + plan_filename = self.get_plan_filename(subject_id, plan_id) + if not (plan_filename.exists() and plan_filename.is_file()): + self.logger.error( + f"Plan file not found for plan {plan_id}, subject {subject_id}" + ) + raise FileNotFoundError( + f"Plan file not found for plan {plan_id}, subject {subject_id}" + ) + plan = Plan.from_file(plan_filename) + self.logger.info(f"Loaded plan {plan_id} for subject {subject_id}.") + return plan + + def delete_plan( + self, + subject_id: str, + plan_id: str, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Delete a Plan and drop it from the subject's plans index. + + Does NOT touch the pre-solutions the Plan references -- those live at subject + scope and may be referenced by the parent PlanningSession or another Plan. + """ + on_conflict = _normalize_on_conflict(on_conflict) + plan_ids = self.get_plan_ids(subject_id) + if plan_id not in plan_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Plan ID {plan_id} does not exist for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Cannot delete plan {plan_id} as it does not exist for subject {subject_id}." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error' or 'skip'.") + + plan_dir = self.get_plan_dir(subject_id, plan_id) + if plan_dir.is_dir(): + shutil.rmtree(plan_dir) + + plan_ids.remove(plan_id) + self.write_plan_ids(subject_id, plan_ids) + + self.logger.info(f"Removed plan {plan_id} for subject {subject_id}.") + + # ------------------------------------------------------------------ + # PlanningSession read/write/delete (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def write_planning_session( + self, + subject_id: str, + planning_session: PlanningSession, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Write a PlanningSession to the subject's ``planning_sessions/`` directory.""" + on_conflict = _normalize_on_conflict(on_conflict) + if planning_session.subject_id is not None and planning_session.subject_id != subject_id: + raise ValueError( + "PlanningSession.subject_id " + f"({planning_session.subject_id!r}) does not match subject_id argument ({subject_id!r})." + ) + # Validate virtual-fit results reference targets on the session. + target_ids = {t.id for t in planning_session.targets} + for target_id, list_of_transforms in planning_session.virtual_fit_results.items(): + if target_id not in target_ids: + raise ValueError( + f"PlanningSession {planning_session.id} virtual_fit_results references target " + f"{target_id} that is not in its targets list." + ) + if len(list_of_transforms) < 1: + raise ValueError( + f"PlanningSession {planning_session.id} virtual_fit_results provides no " + f"transforms for target {target_id}." + ) + + ids = self.get_planning_session_ids(subject_id) + if planning_session.id in ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"PlanningSession with ID {planning_session.id} already exists for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting PlanningSession {planning_session.id} for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping PlanningSession {planning_session.id} for subject {subject_id} as it already exists." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + + filename = self.get_planning_session_filename(subject_id, planning_session.id) + filename.parent.parent.mkdir(parents=True, exist_ok=True) + filename.parent.mkdir(exist_ok=True) + planning_session.update_modified_time() + planning_session.to_file(filename) + + if planning_session.id not in ids: + ids.append(planning_session.id) + self.write_planning_session_ids(subject_id, ids) + + self.logger.info( + f"Wrote PlanningSession {planning_session.id} for subject {subject_id}." + ) + + def load_planning_session(self, subject_id: str, planning_session_id: str) -> PlanningSession: + """Load a PlanningSession by id from the subject's ``planning_sessions/`` directory.""" + filename = self.get_planning_session_filename(subject_id, planning_session_id) + if not (filename.exists() and filename.is_file()): + self.logger.error( + f"PlanningSession file not found for planning session {planning_session_id}, subject {subject_id}" + ) + raise FileNotFoundError( + f"PlanningSession file not found for planning session {planning_session_id}, subject {subject_id}" + ) + ps = PlanningSession.from_file(filename) + self.logger.info( + f"Loaded PlanningSession {planning_session_id} for subject {subject_id}." + ) + return ps + + def delete_planning_session( + self, + subject_id: str, + planning_session_id: str, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Delete a PlanningSession and drop it from the subject's planning-sessions index. + + Does NOT touch pre-solutions on disk or any Plans that were finalized from this + PlanningSession. Those live at subject scope and are independent. + """ + on_conflict = _normalize_on_conflict(on_conflict) + ids = self.get_planning_session_ids(subject_id) + if planning_session_id not in ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"PlanningSession ID {planning_session_id} does not exist for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Cannot delete PlanningSession {planning_session_id} as it does not exist for subject {subject_id}." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error' or 'skip'.") + + session_dir = self.get_planning_session_dir(subject_id, planning_session_id) + if session_dir.is_dir(): + shutil.rmtree(session_dir) + + ids.remove(planning_session_id) + self.write_planning_session_ids(subject_id, ids) + + self.logger.info( + f"Removed PlanningSession {planning_session_id} for subject {subject_id}." + ) + + # ------------------------------------------------------------------ + # SonicationSession read/write/delete (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def write_sonication_session( + self, + subject_id: str, + sonication_session: SonicationSession, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Write a SonicationSession to the subject's ``sonication_sessions/`` directory.""" + on_conflict = _normalize_on_conflict(on_conflict) + if sonication_session.subject_id is not None and sonication_session.subject_id != subject_id: + raise ValueError( + "SonicationSession.subject_id " + f"({sonication_session.subject_id!r}) does not match subject_id argument ({subject_id!r})." + ) + ids = self.get_sonication_session_ids(subject_id) + if sonication_session.id in ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"SonicationSession with ID {sonication_session.id} already exists for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting SonicationSession {sonication_session.id} for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping SonicationSession {sonication_session.id} for subject {subject_id} as it already exists." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + + filename = self.get_sonication_session_filename(subject_id, sonication_session.id) + filename.parent.parent.mkdir(parents=True, exist_ok=True) + filename.parent.mkdir(exist_ok=True) + sonication_session.update_modified_time() + sonication_session.to_file(filename) + + if sonication_session.id not in ids: + ids.append(sonication_session.id) + self.write_sonication_session_ids(subject_id, ids) + + self.logger.info( + f"Wrote SonicationSession {sonication_session.id} for subject {subject_id}." + ) + + def load_sonication_session(self, subject_id: str, sonication_session_id: str) -> SonicationSession: + """Load a SonicationSession by id from the subject's ``sonication_sessions/`` directory.""" + filename = self.get_sonication_session_filename(subject_id, sonication_session_id) + if not (filename.exists() and filename.is_file()): + self.logger.error( + f"SonicationSession file not found for sonication session {sonication_session_id}, subject {subject_id}" + ) + raise FileNotFoundError( + f"SonicationSession file not found for sonication session {sonication_session_id}, subject {subject_id}" + ) + ss = SonicationSession.from_file(filename) + self.logger.info( + f"Loaded SonicationSession {sonication_session_id} for subject {subject_id}." + ) + return ss + + def delete_sonication_session( + self, + subject_id: str, + sonication_session_id: str, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Delete a SonicationSession and drop it from the subject's sonication-sessions index. + + Does NOT delete the referenced Plan, subject-scoped photoscans, or subject-scoped + solutions -- those are independent artifacts. + """ + on_conflict = _normalize_on_conflict(on_conflict) + ids = self.get_sonication_session_ids(subject_id) + if sonication_session_id not in ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"SonicationSession ID {sonication_session_id} does not exist for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Cannot delete SonicationSession {sonication_session_id} as it does not exist for subject {subject_id}." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error' or 'skip'.") + + session_dir = self.get_sonication_session_dir(subject_id, sonication_session_id) + if session_dir.is_dir(): + shutil.rmtree(session_dir) + + ids.remove(sonication_session_id) + self.write_sonication_session_ids(subject_id, ids) + + self.logger.info( + f"Removed SonicationSession {sonication_session_id} for subject {subject_id}." + ) + + # ------------------------------------------------------------------ + # Subject-scoped Photoscan storage (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def write_photoscan_at_subject_scope( + self, + subject_id: str, + photoscan: Photoscan, + model_data_filepath: str | None = None, + texture_data_filepath: str | None = None, + mtl_data_filepath: str | None = None, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Write a Photoscan to ``subjects/{subject_id}/photoscans/{photoscan_id}/``. + + Physical files live at subject scope; logical ownership is per-SonicationSession + (via ``SonicationSession.photoscan_ids``). The same photoscan is not duplicated + on disk when it's referenced by multiple sonication sessions. + + The model data file is required on first write; on overwrite, the caller may + omit the data filepaths to keep the existing files in place. Texture and MTL + files are optional. + """ + on_conflict = _normalize_on_conflict(on_conflict) + photoscan_ids = self.get_subject_photoscan_ids(subject_id) + if photoscan.id in photoscan_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Photoscan with ID {photoscan.id} already exists at subject scope for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.OVERWRITE: + self.logger.info( + f"Overwriting photoscan {photoscan.id} at subject scope for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Skipping photoscan {photoscan.id} at subject scope as it already exists." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error', 'overwrite', or 'skip'.") + + photoscan_metadata_filepath = self.get_subject_photoscan_metadata_filepath(subject_id, photoscan.id) + photoscan_parent_dir = photoscan_metadata_filepath.parent + photoscan_parent_dir.mkdir(parents=True, exist_ok=True) + + if model_data_filepath: + model_data_filepath = Path(model_data_filepath) + if not model_data_filepath.exists(): + raise FileNotFoundError(f"Model data filepath does not exist: {model_data_filepath}") + photoscan.model_filename = model_data_filepath.name + shutil.copy(model_data_filepath, photoscan_parent_dir) + elif not photoscan.model_filename or not (photoscan_parent_dir / photoscan.model_filename).exists(): + raise ValueError(f"Cannot find model file associated with photoscan {photoscan.id}.") + + if texture_data_filepath: + texture_data_filepath = Path(texture_data_filepath) + if not texture_data_filepath.exists(): + raise FileNotFoundError(f"Texture data filepath does not exist: {texture_data_filepath}") + photoscan.texture_filename = texture_data_filepath.name + shutil.copy(texture_data_filepath, photoscan_parent_dir) + elif photoscan.texture_filename and not (photoscan_parent_dir / photoscan.texture_filename).exists(): + raise ValueError(f"Cannot find texture file associated with photoscan {photoscan.id}.") + + if mtl_data_filepath: + mtl_data_filepath = Path(mtl_data_filepath) + if not mtl_data_filepath.exists(): + raise FileNotFoundError(f"MTL filepath does not exist: {mtl_data_filepath}") + photoscan.mtl_filename = mtl_data_filepath.name + shutil.copy(mtl_data_filepath, photoscan_parent_dir) + elif photoscan.mtl_filename and not (photoscan_parent_dir / photoscan.mtl_filename).exists(): + raise ValueError(f"Cannot find photoscan materials file associated with photoscan {photoscan.id}.") + + photoscan.to_file(photoscan_metadata_filepath) + + if photoscan.id not in photoscan_ids: + photoscan_ids.append(photoscan.id) + self.write_subject_photoscan_ids(subject_id, photoscan_ids) + + self.logger.info( + f"Wrote photoscan {photoscan.id} at subject scope for subject {subject_id}." + ) + + def load_photoscan_at_subject_scope( + self, + subject_id: str, + photoscan_id: str, + load_data: bool = False, + ): + """Load a subject-scoped Photoscan and optionally its model + texture data. + + Returns ``photoscan`` if ``load_data=False``, else + ``(photoscan, (model_data, texture_data))`` matching the legacy + ``load_photoscan`` shape. + """ + photoscan_metadata_filepath = self.get_subject_photoscan_metadata_filepath(subject_id, photoscan_id) + if not (photoscan_metadata_filepath.exists() and photoscan_metadata_filepath.is_file()): + raise FileNotFoundError( + f"Photoscan file not found at subject scope for photoscan {photoscan_id}, subject {subject_id}" + ) + photoscan = Photoscan.from_file(photoscan_metadata_filepath) + if load_data: + model_data, texture_data = load_data_from_photoscan( + photoscan, Path(photoscan_metadata_filepath.parent), + ) + return photoscan, (model_data, texture_data) + return photoscan + + def get_photoscan_absolute_filepaths_info_at_subject_scope( + self, subject_id: str, photoscan_id: str, + ) -> Dict: + """Return a dict of absolute file paths for the subject-scoped photoscan.""" + photoscan_metadata_filepath = self.get_subject_photoscan_metadata_filepath(subject_id, photoscan_id) + photoscan_metadata_directory = Path(photoscan_metadata_filepath).parent + with open(photoscan_metadata_filepath) as f: + photoscan = json.load(f) + info = { + "id": photoscan["id"], + "name": photoscan["name"], + "model_abspath": photoscan_metadata_directory / photoscan["model_filename"], + "texture_abspath": photoscan_metadata_directory / photoscan["texture_filename"], + "photoscan_approved": photoscan["photoscan_approved"], + } + if "mtl_filename" in photoscan: + info["mtl_abspath"] = photoscan_metadata_directory / photoscan["mtl_filename"] + return info + + def delete_photoscan_at_subject_scope( + self, + subject_id: str, + photoscan_id: str, + on_conflict: OnConflictOpts | str = OnConflictOpts.ERROR, + ) -> None: + """Delete a subject-scoped Photoscan and drop it from the subject's photoscans index. + + Does NOT scrub SonicationSessions that reference this photoscan. Cleaning up + dangling references is the caller's responsibility (a session-level compact + step, not a DB-level cascade). + """ + on_conflict = _normalize_on_conflict(on_conflict) + photoscan_ids = self.get_subject_photoscan_ids(subject_id) + if photoscan_id not in photoscan_ids: + if on_conflict == OnConflictOpts.ERROR: + raise ValueError( + f"Photoscan ID {photoscan_id} does not exist at subject scope for subject {subject_id}." + ) + elif on_conflict == OnConflictOpts.SKIP: + self.logger.info( + f"Cannot delete photoscan {photoscan_id} at subject scope as it does not exist for subject {subject_id}." + ) + return + else: + raise ValueError("Invalid 'on_conflict' option. Use 'error' or 'skip'.") + + photoscan_dir = self.get_subject_photoscan_dir(subject_id, photoscan_id) + if photoscan_dir.is_dir(): + shutil.rmtree(photoscan_dir) + + photoscan_ids.remove(photoscan_id) + self.write_subject_photoscan_ids(subject_id, photoscan_ids) + + self.logger.info( + f"Removed photoscan {photoscan_id} at subject scope for subject {subject_id}." + ) + def choose_session(self, subject, options=None): # Implement the logic to choose a session raise NotImplementedError("Method not yet implemented") @@ -670,6 +1417,24 @@ def get_solution_ids(self, subject_id:str, session_id:str) -> List[str]: return json.loads(solutions_filename.read_text())["solution_ids"] + def get_subject_solution_ids(self, subject_id: str) -> List[str]: + """List IDs of all solutions stored under the subject's subject-scoped ``solutions/`` directory. + + Independent of the legacy session-scoped ``get_solution_ids``. Returns ``[]`` if the + subject has no solutions index file yet. + """ + solutions_filename = self.get_subject_solutions_filename(subject_id) + if not (solutions_filename.exists() and solutions_filename.is_file()): + self.logger.info("Subject-scoped solutions file not found for subject %s.", subject_id) + return [] + return json.loads(solutions_filename.read_text())["solution_ids"] + + def write_subject_solution_ids(self, subject_id: str, solution_ids: List[str]) -> None: + """Overwrite the subject-scoped solutions index.""" + solutions_filepath = self.get_subject_solutions_filename(subject_id) + solutions_filepath.parent.mkdir(parents=True, exist_ok=True) + solutions_filepath.write_text(json.dumps({"solution_ids": solution_ids})) + def get_photocollection_reference_numbers(self, subject_id: str, session_id: str) -> List[str]: """Get a list of reference numbers of the photocollections associated with the given session""" photocollection_filename = self.get_photocollections_filename(subject_id, session_id) @@ -912,7 +1677,48 @@ def load_session(self, subject, session_id, options=None): options = {} session_filename = self.get_session_filename(subject.id, session_id) if os.path.isfile(session_filename): + # Read raw JSON first so we can detect legacy sessions that lack + # the photoscans / photocollections index fields and migrate them + # from the on-disk index files. Present-but-empty lists are + # respected (user may have intentionally cleared them). + with open(session_filename) as _legacy_fh: + _raw_session_dict = json.load(_legacy_fh) session = Session.from_file(session_filename) + if 'photoscans' not in _raw_session_dict: + try: + session.photoscans = list(self.get_photoscan_ids(subject.id, session_id) or []) + except Exception: # pylint: disable=broad-exception-caught + session.photoscans = [] + if 'photocollections' not in _raw_session_dict: + try: + session.photocollections = list(self.get_photocollection_reference_numbers(subject.id, session_id) or []) + except Exception: # pylint: disable=broad-exception-caught + session.photocollections = [] + # Drop any transducer_tracking_results whose photoscan_id no + # longer exists in this session's photoscans index. The session + # JSON and the photoscans index can drift out of sync (e.g. a + # photoscan was deleted from the index without the session JSON + # being rewritten), and ``write_session`` rejects sessions whose + # tracking results reference unknown photoscans. Sanitizing on + # load means a freshly-loaded session can round-trip through + # ``write_session`` without surprise validation errors on exit. + if session.transducer_tracking_results: + indexed_photoscan_ids = set(self.get_photoscan_ids(subject.id, session_id)) + kept: list = [] + dropped: list[str] = [] + for result in session.transducer_tracking_results: + if result.photoscan_id in indexed_photoscan_ids: + kept.append(result) + else: + dropped.append(result.photoscan_id) + if dropped: + self.logger.warning( + "Dropping %d transducer_tracking_result(s) from session " + "%s of subject %s referencing photoscan id(s) not in " + "this session's photoscans index: %s", + len(dropped), session_id, subject.id, sorted(set(dropped)), + ) + session.transducer_tracking_results = kept self.logger.info(f"Loaded session {session_id} for subject {subject.id}") return session else: @@ -1022,11 +1828,53 @@ def get_solution_filepath(self, subject_id, session_id, solution_id) -> Path: session_dir = self.get_session_dir(subject_id, session_id) return Path(session_dir) / 'solutions' / solution_id / f"{solution_id}.json" + def get_solution_analysis_filepath(self, subject_id, session_id, solution_id) -> Path: + """Get the solution-analysis json file for the solution with the given ID. + + Lives next to the solution itself so the analysis and the solution that produced it stay together. + """ + session_dir = self.get_session_dir(subject_id, session_id) + return Path(session_dir) / 'solutions' / solution_id / f"{solution_id}_analysis.json" + def get_solutions_filename(self, subject_id, session_id) -> Path: """Get the path to the overall solutions json file for the requested session""" session_dir = self.get_session_dir(subject_id, session_id) return Path(session_dir) / 'solutions' / 'solutions.json' + # ------------------------------------------------------------------ + # Subject-scoped solution storage (see SESSION_SPLIT_DESIGN.md). + # + # In the split-session model, Solutions live under the subject rather + # than under a specific Session. This lets a Solution be referenced + # (via SolutionInfo) by a PlanningSession's pre-solutions list, a + # finalized Plan's pre-solutions list, or a SonicationSession's + # final solution field, without ever being duplicated on disk. + # ------------------------------------------------------------------ + + def get_subject_solutions_filename(self, subject_id: str) -> Path: + """Path to the subject-scoped solutions index (``subjects/{sid}/solutions/solutions.json``).""" + return Path(self.get_subject_dir(subject_id)) / 'solutions' / 'solutions.json' + + def get_subject_solution_dir(self, subject_id: str, solution_id: str) -> Path: + """Directory holding a subject-scoped solution's files.""" + return Path(self.get_subject_dir(subject_id)) / 'solutions' / solution_id + + def get_subject_solution_filepath(self, subject_id: str, solution_id: str) -> Path: + """Path to a subject-scoped solution's JSON file (``{solution_id}.solution.json``). + + The ``.solution.json`` extension identifies the file type so the raw filename + (rather than an id-suffix convention) carries the discriminator across sibling + tables in a future relational-DB migration. + """ + return self.get_subject_solution_dir(subject_id, solution_id) / f"{solution_id}.solution.json" + + def get_subject_solution_analysis_filepath(self, subject_id: str, solution_id: str) -> Path: + """Path to a subject-scoped solution's analysis JSON, sitting next to the solution. + + Uses ``.solution_analysis.json`` to match the ``.solution.json`` companion. + """ + return self.get_subject_solution_dir(subject_id, solution_id) / f"{solution_id}.solution_analysis.json" + def get_photocollections_filename(self, subject_id, session_id) -> Path: """Get the path to the overall photocollections json file for the requested session""" session_dir = self.get_session_dir(subject_id, session_id) @@ -1064,6 +1912,63 @@ def get_volume_metadata_filepath(self, subject_id, volume_id): def get_photoscan_metadata_filepath(self, subject_id, session_id, photoscan_id): return Path(self.get_session_dir(subject_id, session_id)) / 'photoscans' / photoscan_id / f'{photoscan_id}.json' + # ------------------------------------------------------------------ + # Split-session path helpers (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def get_plans_filename(self, subject_id: str) -> Path: + """Path to the subject-scoped plans index (``subjects/{sid}/plans/plans.json``).""" + return Path(self.get_subject_dir(subject_id)) / 'plans' / 'plans.json' + + def get_plan_dir(self, subject_id: str, plan_id: str) -> Path: + """Directory holding a Plan's files.""" + return Path(self.get_subject_dir(subject_id)) / 'plans' / plan_id + + def get_plan_filename(self, subject_id: str, plan_id: str) -> Path: + """Path to a Plan's JSON file (``{plan_id}.plan.json``). + + The ``.plan.json`` extension identifies the file type so the id itself does + not need a ``_plan`` suffix; the same ``id`` can be reused across a + PlanningSession, a Plan finalized from it, and a SonicationSession using it. + """ + return self.get_plan_dir(subject_id, plan_id) / f'{plan_id}.plan.json' + + def get_planning_sessions_filename(self, subject_id: str) -> Path: + """Path to the subject-scoped planning-sessions index.""" + return Path(self.get_subject_dir(subject_id)) / 'planning_sessions' / 'planning_sessions.json' + + def get_planning_session_dir(self, subject_id: str, planning_session_id: str) -> Path: + """Directory holding a PlanningSession's files.""" + return Path(self.get_subject_dir(subject_id)) / 'planning_sessions' / planning_session_id + + def get_planning_session_filename(self, subject_id: str, planning_session_id: str) -> Path: + """Path to a PlanningSession's JSON file (``{planning_session_id}.planning.json``).""" + return self.get_planning_session_dir(subject_id, planning_session_id) / f'{planning_session_id}.planning.json' + + def get_sonication_sessions_filename(self, subject_id: str) -> Path: + """Path to the subject-scoped sonication-sessions index.""" + return Path(self.get_subject_dir(subject_id)) / 'sonication_sessions' / 'sonication_sessions.json' + + def get_sonication_session_dir(self, subject_id: str, sonication_session_id: str) -> Path: + """Directory holding a SonicationSession's files.""" + return Path(self.get_subject_dir(subject_id)) / 'sonication_sessions' / sonication_session_id + + def get_sonication_session_filename(self, subject_id: str, sonication_session_id: str) -> Path: + """Path to a SonicationSession's JSON file (``{sonication_session_id}.sonication.json``).""" + return self.get_sonication_session_dir(subject_id, sonication_session_id) / f'{sonication_session_id}.sonication.json' + + def get_subject_photoscans_filename(self, subject_id: str) -> Path: + """Path to the subject-scoped photoscans index.""" + return Path(self.get_subject_dir(subject_id)) / 'photoscans' / 'photoscans.json' + + def get_subject_photoscan_dir(self, subject_id: str, photoscan_id: str) -> Path: + """Directory holding a subject-scoped Photoscan's files.""" + return Path(self.get_subject_dir(subject_id)) / 'photoscans' / photoscan_id + + def get_subject_photoscan_metadata_filepath(self, subject_id: str, photoscan_id: str) -> Path: + """Path to a subject-scoped Photoscan's metadata JSON.""" + return self.get_subject_photoscan_dir(subject_id, photoscan_id) / f'{photoscan_id}.json' + def get_run_dir(self, subject_id, session_id, run_id): run_dir = self.get_session_dir(subject_id, session_id) / 'runs' / f'{run_id}' return run_dir @@ -1143,6 +2048,70 @@ def write_solution_ids(self, session:Session, solution_ids:List[str]): solutions_filepath.parent.mkdir(exist_ok=True) # Make solutions directory in case it does not exist solutions_filepath.write_text(json.dumps(solutions_data)) + # ------------------------------------------------------------------ + # Split-session id-index read/write helpers (see SESSION_SPLIT_DESIGN.md). + # ------------------------------------------------------------------ + + def get_plan_ids(self, subject_id: str) -> List[str]: + """List IDs of Plans stored under ``subjects/{sid}/plans/``.""" + plans_filename = self.get_plans_filename(subject_id) + if not (plans_filename.exists() and plans_filename.is_file()): + self.logger.info("Plans file not found for subject %s.", subject_id) + return [] + return json.loads(plans_filename.read_text()).get("plan_ids", []) + + def write_plan_ids(self, subject_id: str, plan_ids: List[str]) -> None: + """Overwrite the subject-scoped plans index.""" + plans_filename = self.get_plans_filename(subject_id) + plans_filename.parent.mkdir(parents=True, exist_ok=True) + plans_filename.write_text(json.dumps({"plan_ids": plan_ids})) + + def get_planning_session_ids(self, subject_id: str) -> List[str]: + """List IDs of PlanningSessions stored under ``subjects/{sid}/planning_sessions/``.""" + idx = self.get_planning_sessions_filename(subject_id) + if not (idx.exists() and idx.is_file()): + self.logger.info("Planning-sessions file not found for subject %s.", subject_id) + return [] + return json.loads(idx.read_text()).get("planning_session_ids", []) + + def write_planning_session_ids(self, subject_id: str, ids: List[str]) -> None: + """Overwrite the subject-scoped planning-sessions index.""" + idx = self.get_planning_sessions_filename(subject_id) + idx.parent.mkdir(parents=True, exist_ok=True) + idx.write_text(json.dumps({"planning_session_ids": ids})) + + def get_sonication_session_ids(self, subject_id: str) -> List[str]: + """List IDs of SonicationSessions stored under ``subjects/{sid}/sonication_sessions/``.""" + idx = self.get_sonication_sessions_filename(subject_id) + if not (idx.exists() and idx.is_file()): + self.logger.info("Sonication-sessions file not found for subject %s.", subject_id) + return [] + return json.loads(idx.read_text()).get("sonication_session_ids", []) + + def write_sonication_session_ids(self, subject_id: str, ids: List[str]) -> None: + """Overwrite the subject-scoped sonication-sessions index.""" + idx = self.get_sonication_sessions_filename(subject_id) + idx.parent.mkdir(parents=True, exist_ok=True) + idx.write_text(json.dumps({"sonication_session_ids": ids})) + + def get_subject_photoscan_ids(self, subject_id: str) -> List[str]: + """List IDs of photoscans stored under ``subjects/{sid}/photoscans/``. + + Independent of the legacy session-scoped ``get_photoscan_ids``. Returns ``[]`` + when the subject has no subject-scoped photoscans index yet. + """ + idx = self.get_subject_photoscans_filename(subject_id) + if not (idx.exists() and idx.is_file()): + self.logger.info("Subject-scoped photoscans file not found for subject %s.", subject_id) + return [] + return json.loads(idx.read_text()).get("photoscan_ids", []) + + def write_subject_photoscan_ids(self, subject_id: str, photoscan_ids: List[str]) -> None: + """Overwrite the subject-scoped photoscans index.""" + idx = self.get_subject_photoscans_filename(subject_id) + idx.parent.mkdir(parents=True, exist_ok=True) + idx.write_text(json.dumps({"photoscan_ids": photoscan_ids})) + @staticmethod def get_default_user_dir(): """ diff --git a/src/openlifu/db/plan.py b/src/openlifu/db/plan.py new file mode 100644 index 00000000..7cec76ea --- /dev/null +++ b/src/openlifu/db/plan.py @@ -0,0 +1,146 @@ +from __future__ import annotations + +import copy +import json +from dataclasses import asdict, dataclass, field +from datetime import datetime +from pathlib import Path +from typing import Annotated, Dict, List + +import numpy as np + +from openlifu.db.session import SolutionInfo +from openlifu.geo.point import Point +from openlifu.geo.transforms import ArrayTransform +from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.json import PYFUSEncoder +from openlifu.util.strings import sanitize + + +@dataclass +class Plan: + """Immutable finalized treatment plan produced by a :class:`PlanningSession`. + + A ``Plan`` pins down: which target, which transducer at which pose (the virtual-fit + pose that was approved on the parent PlanningSession), and against which protocol. + Optionally carries :class:`SolutionInfo` references to pre-solutions computed at + planning time so the sonication operator can review expected pressures before + treatment. + + A ``Plan`` is frozen once written. Editing means "create a new Plan from an updated + PlanningSession". A single PlanningSession may finalize multiple Plans over time. + + A :class:`SonicationSession` references a ``Plan`` by ``id`` (frozen input); the + Plan's ``target``, ``volume_id``, ``protocol_id``, ``transducer_id``, and + ``array_transform`` are the reference "what we're trying to hit" for the sonication. + """ + + id: Annotated[str | None, OpenLIFUFieldData("Plan ID", "ID of this plan")] = None + """ID of this plan""" + + name: Annotated[str | None, OpenLIFUFieldData("Plan name", "Plan name")] = None + """Plan name""" + + subject_id: Annotated[str | None, OpenLIFUFieldData("Subject ID", "ID of the parent subject of this plan")] = None + """ID of the parent subject of this plan""" + + volume_id: Annotated[str | None, OpenLIFUFieldData("Volume ID", "ID of the subject volume this plan was finalized against")] = None + """ID of the subject volume this plan was finalized against""" + + protocol_id: Annotated[str | None, OpenLIFUFieldData("Protocol ID", "ID of the protocol this plan was finalized against")] = None + """ID of the protocol this plan was finalized against""" + + transducer_id: Annotated[str | None, OpenLIFUFieldData("Transducer ID", "ID of the transducer this plan was finalized against")] = None + """ID of the transducer this plan was finalized against""" + + target: Annotated[Point | None, OpenLIFUFieldData("Target", "The Point target this plan committed to sonicating.")] = None + """The :class:`Point` target this plan committed to sonicating.""" + + array_transform: Annotated[ArrayTransform, OpenLIFUFieldData("Array transform", "The transducer affine transform matrix with units representing the approved virtual-fit pose this plan committed to.")] = field(default_factory=lambda: ArrayTransform(np.eye(4), "mm")) + """The transducer affine transform matrix with units representing the approved + virtual-fit pose this plan committed to.""" + + pre_solutions: Annotated[List[SolutionInfo], OpenLIFUFieldData("Pre-solutions", "SolutionInfo references to pre-solutions computed at planning time. The actual Solution files live at subject scope (``subjects/{sid}/solutions/``) and are shared with the parent PlanningSession's pre_solutions list. Empty if the user finalized the plan without computing pre-solutions.")] = field(default_factory=list) + """:class:`SolutionInfo` references to pre-solutions computed at planning time. + + The actual :class:`~openlifu.plan.Solution` files live at subject scope + (``subjects/{sid}/solutions/``) and are shared with the parent + :class:`PlanningSession`'s ``pre_solutions`` list. Empty if the user finalized the + plan without computing pre-solutions. + """ + + parent_planning_session_id: Annotated[str | None, OpenLIFUFieldData("Parent planning session ID", "ID of the PlanningSession that finalized this Plan. Provenance only; a Plan does not depend on its parent still existing.")] = None + """ID of the PlanningSession that finalized this Plan. Provenance only; a Plan does + not depend on its parent still existing.""" + + date_created: Annotated[datetime, OpenLIFUFieldData("Date created", "Date the plan was finalized")] = field(default_factory=datetime.now) + """Date the plan was finalized""" + + notes: Annotated[str, OpenLIFUFieldData("Plan notes", "Free-form notes recorded at finalization time.")] = "" + """Free-form notes recorded at finalization time.""" + + attrs: Annotated[dict, OpenLIFUFieldData("Custom attributes", "Dictionary of additional custom attributes to save with the plan")] = field(default_factory=dict) + """Dictionary of additional custom attributes to save with the plan""" + + def __post_init__(self): + if self.id is None and self.name is None: + self.id = "plan" + if self.id is None: + self.id = sanitize(self.name, "snake") + if self.name is None: + self.name = self.id + + @staticmethod + def from_file(filename) -> Plan: + """Load a Plan from a JSON file.""" + with open(filename) as f: + return Plan.from_dict(json.load(f)) + + @staticmethod + def from_dict(d: Dict) -> Plan: + """Reconstruct a Plan from its dictionary representation.""" + d = dict(d) # shallow copy; we mutate below + if "date_created" in d and isinstance(d["date_created"], str): + d["date_created"] = datetime.fromisoformat(d["date_created"]) + if "array_transform" in d and isinstance(d["array_transform"], dict): + d["array_transform"] = ArrayTransform.from_dict(d["array_transform"]) + if "target" in d and isinstance(d["target"], dict): + d["target"] = Point.from_dict(d["target"]) + if "pre_solutions" in d: + d["pre_solutions"] = [ + s if isinstance(s, SolutionInfo) else SolutionInfo(**s) + for s in d["pre_solutions"] + ] + return Plan(**d) + + def to_dict(self) -> Dict: + """Serialize the Plan to a dictionary.""" + d = copy.deepcopy(self.__dict__) + d["date_created"] = d["date_created"].isoformat() + d["array_transform"] = asdict(d["array_transform"]) + if d["target"] is not None: + d["target"] = d["target"].to_dict() + d["pre_solutions"] = [asdict(s) for s in d["pre_solutions"]] + return d + + @staticmethod + def from_json(json_string: str) -> Plan: + """Load a Plan from a JSON string.""" + return Plan.from_dict(json.loads(json_string)) + + def to_json(self, compact: bool) -> str: + """Serialize a Plan to a JSON string. + + Args: + compact: if enabled then the string is compact (not pretty). Disable for pretty. + """ + if compact: + return json.dumps(self.to_dict(), separators=(",", ":"), cls=PYFUSEncoder) + return json.dumps(self.to_dict(), indent=4, cls=PYFUSEncoder) + + def to_file(self, filename) -> None: + """Write the Plan to a JSON file, creating parent directories as needed.""" + Path(filename).parent.parent.mkdir(exist_ok=True) # plans directory + Path(filename).parent.mkdir(exist_ok=True) # {plan_id} directory + with open(filename, "w") as f: + f.write(self.to_json(compact=False)) diff --git a/src/openlifu/db/planning_session.py b/src/openlifu/db/planning_session.py new file mode 100644 index 00000000..7752325d --- /dev/null +++ b/src/openlifu/db/planning_session.py @@ -0,0 +1,176 @@ +from __future__ import annotations + +import copy +import json +from dataclasses import asdict, dataclass, field +from datetime import datetime +from pathlib import Path +from typing import Annotated, Dict, List, Tuple + +from openlifu.db.session import SolutionInfo +from openlifu.geo.point import Point +from openlifu.geo.transforms import ArrayTransform +from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.json import PYFUSEncoder +from openlifu.util.strings import sanitize + + +@dataclass +class PlanningSession: + """A working document for developing a treatment :class:`Plan`. + + Owns targets, virtual-fit results, and pre-solutions. Mutable throughout its + lifetime; only explicit :meth:`to_file` or a database ``write_planning_session`` + call touches disk. + + "Finalizing" a PlanningSession produces an immutable :class:`Plan` (a + ``db.finalize_plan(...)`` helper, added in a later commit). The same + PlanningSession may finalize multiple Plans over time; every Plan id it has + produced is tracked in :attr:`finalized_plan_ids`. + + ``pre_solutions`` entries are :class:`SolutionInfo` references to Solutions + stored at subject scope (``subjects/{sid}/solutions/``). The same Solution + file may be referenced by this session's ``pre_solutions`` list and one or + more Plans finalized from it, without duplication on disk. + """ + + id: Annotated[str | None, OpenLIFUFieldData("Planning session ID", "ID of this planning session")] = None + """ID of this planning session""" + + name: Annotated[str | None, OpenLIFUFieldData("Planning session name", "Planning session name")] = None + """Planning session name""" + + subject_id: Annotated[str | None, OpenLIFUFieldData("Subject ID", "ID of the parent subject of this planning session")] = None + """ID of the parent subject of this planning session""" + + volume_id: Annotated[str | None, OpenLIFUFieldData("Volume ID", "ID of the subject volume associated with this planning session")] = None + """ID of the subject volume associated with this planning session""" + + protocol_id: Annotated[str | None, OpenLIFUFieldData("Protocol ID", "ID of the protocol used for this planning session")] = None + """ID of the protocol used for this planning session""" + + transducer_id: Annotated[str | None, OpenLIFUFieldData("Transducer ID", "ID of the transducer associated with this planning session")] = None + """ID of the transducer associated with this planning session""" + + date_created: Annotated[datetime, OpenLIFUFieldData("Date created", "Date of creation of the planning session")] = field(default_factory=datetime.now) + """Date of creation of the planning session""" + + date_modified: Annotated[datetime, OpenLIFUFieldData("Date modified", "Date of last modification of the planning session")] = field(default_factory=datetime.now) + """Date of last modification of the planning session""" + + targets: Annotated[List[Point], OpenLIFUFieldData("Targets", "Targets saved to this planning session")] = field(default_factory=list) + """Targets saved to this planning session""" + + virtual_fit_results: Annotated[Dict[str, List[Tuple[bool, ArrayTransform]]], OpenLIFUFieldData("Virtual fit results", "Dictionary mapping target IDs to a list of (approval, transform) pairs, one entry per virtual-fit candidate for that target.")] = field(default_factory=dict) + """Virtual fit results. + + Dictionary mapping target IDs to a list of ``(approval, transform)`` pairs, where + ``approval`` is a boolean indicating whether the specific virtual fit ``transform`` + has been approved for sonication, and ``transform`` is a transducer array-to-volume + transform generated by the virtual fit for that target. + """ + + pre_solutions: Annotated[List[SolutionInfo], OpenLIFUFieldData("Pre-solutions", "SolutionInfo references to pre-solutions computed against virtual-fit poses on this planning session. Actual Solution files live at subject scope.")] = field(default_factory=list) + """:class:`SolutionInfo` references to pre-solutions computed against virtual-fit + poses on this planning session. + + The actual :class:`~openlifu.plan.Solution` files live at subject scope + (``subjects/{sid}/solutions/``). Solutions on disk under the subject that have no + matching ``SolutionInfo`` here (and are not referenced from any Plan or + SonicationSession) are candidates for orphan cleanup, but that cleanup is not + performed automatically by session save. + """ + + finalized_plan_ids: Annotated[List[str], OpenLIFUFieldData("Finalized plan IDs", "IDs of every Plan this planning session has finalized. Each Plan is immutable; multiple Plans per session are supported.")] = field(default_factory=list) + """IDs of every :class:`Plan` this planning session has finalized. Each Plan is + immutable; multiple Plans per session are supported.""" + + attrs: Annotated[dict, OpenLIFUFieldData("Custom attributes", "Dictionary of additional custom attributes to save to the planning session")] = field(default_factory=dict) + """Dictionary of additional custom attributes to save to the planning session""" + + def __post_init__(self): + if self.id is None and self.name is None: + self.id = "planning_session" + if self.id is None: + self.id = sanitize(self.name, "snake") + if self.name is None: + self.name = self.id + if isinstance(self.targets, Point): + self.targets = [self.targets] + else: + self.targets = list(self.targets) + + @staticmethod + def from_file(filename) -> PlanningSession: + """Load a PlanningSession from a JSON file.""" + with open(filename) as f: + return PlanningSession.from_dict(json.load(f)) + + @staticmethod + def from_dict(d: Dict) -> PlanningSession: + """Reconstruct a PlanningSession from its dictionary representation.""" + d = dict(d) # shallow copy; we mutate below + if "date_created" in d and isinstance(d["date_created"], str): + d["date_created"] = datetime.fromisoformat(d["date_created"]) + if "date_modified" in d and isinstance(d["date_modified"], str): + d["date_modified"] = datetime.fromisoformat(d["date_modified"]) + if "targets" in d: + targets = d["targets"] + if isinstance(targets, list): + d["targets"] = [ + p if isinstance(p, Point) else Point.from_dict(p) for p in targets + ] + elif isinstance(targets, dict): + d["targets"] = [Point.from_dict(targets)] + elif isinstance(targets, Point): + d["targets"] = [targets] + if "virtual_fit_results" in d: + d["virtual_fit_results"] = { + target_id: [ + (approval, t if isinstance(t, ArrayTransform) else ArrayTransform.from_dict(t)) + for approval, t in list_of_transforms + ] + for target_id, list_of_transforms in d["virtual_fit_results"].items() + } + if "pre_solutions" in d: + d["pre_solutions"] = [ + s if isinstance(s, SolutionInfo) else SolutionInfo(**s) + for s in d["pre_solutions"] + ] + return PlanningSession(**d) + + def to_dict(self) -> Dict: + """Serialize the PlanningSession to a dictionary.""" + d = copy.deepcopy(self.__dict__) + d["date_created"] = d["date_created"].isoformat() + d["date_modified"] = d["date_modified"].isoformat() + d["targets"] = [p.to_dict() for p in d["targets"]] + d["virtual_fit_results"] = { + target_id: [(approval, asdict(t)) for [approval, t] in list_of_transforms] + for target_id, list_of_transforms in d["virtual_fit_results"].items() + } + d["pre_solutions"] = [asdict(s) for s in d["pre_solutions"]] + return d + + @staticmethod + def from_json(json_string: str) -> PlanningSession: + """Load a PlanningSession from a JSON string.""" + return PlanningSession.from_dict(json.loads(json_string)) + + def to_json(self, compact: bool) -> str: + """Serialize a PlanningSession to a JSON string.""" + if compact: + return json.dumps(self.to_dict(), separators=(",", ":"), cls=PYFUSEncoder) + return json.dumps(self.to_dict(), indent=4, cls=PYFUSEncoder) + + def to_file(self, filename) -> None: + """Write the PlanningSession to a JSON file, creating parent directories as needed.""" + Path(filename).parent.parent.mkdir(exist_ok=True) # planning_sessions directory + Path(filename).parent.mkdir(exist_ok=True) # {planning_session_id} directory + with open(filename, "w") as f: + f.write(self.to_json(compact=False)) + + def update_modified_time(self, time: datetime | None = None) -> None: + if time is None: + time = datetime.now() + self.date_modified = time diff --git a/src/openlifu/db/session.py b/src/openlifu/db/session.py index aac35bcf..3de0ca1b 100644 --- a/src/openlifu/db/session.py +++ b/src/openlifu/db/session.py @@ -5,7 +5,7 @@ from dataclasses import asdict, dataclass, field from datetime import datetime from pathlib import Path -from typing import Annotated, Dict, List, Tuple +from typing import Annotated, ClassVar, Dict, FrozenSet, List, Tuple import numpy as np @@ -16,11 +16,38 @@ from openlifu.util.strings import sanitize +@dataclass +class PhotoscanRegistration: + """A registration of a photoscan model into the volume coordinate frame. + + A photoscan may have multiple registration attempts stored on the session; at most one is + expected to be approved at a time, though this is not enforced at the dataclass level. + Downstream transducer-tracking results refer back to a specific registration by ``id``. + """ + + photoscan_id: Annotated[str, OpenLIFUFieldData("Photoscan ID", "ID of the photoscan that this registration applies to")] + """ID of the photoscan that this registration applies to""" + + transform: Annotated[ArrayTransform, OpenLIFUFieldData("Photoscan-to-volume transform", "Transform that registers the photoscan model to the volume's skin segmentation")] + """Transform that registers the photoscan model to the volume's skin segmentation""" + + approval: Annotated[bool, OpenLIFUFieldData("Registration approved?", "Whether this photoscan registration has been approved by the user.")] = False + """Whether this photoscan registration has been approved by the user.""" + + id: Annotated[str | None, OpenLIFUFieldData("Registration ID", "Stable identifier for this photoscan registration, unique within the session. Survives reordering or deletion of other registrations so that downstream references (e.g. transducer tracking results) remain valid.")] = None + """Stable identifier for this photoscan registration, unique within the session.""" + + @dataclass class TransducerTrackingResult: """ Class representing the results of running the transducer tracking algorithm. + + Each result registers a transducer pose against a volume in the context of a particular + photoscan registration. The photoscan-to-volume transform that the tracking was performed + against is stored separately as a :class:`PhotoscanRegistration` and referenced here by + ``photoscan_registration_id``. """ photoscan_id: Annotated[str, OpenLIFUFieldData("Photoscan ID", "ID of the photoscan object used for transducer tracking")] @@ -29,16 +56,105 @@ class TransducerTrackingResult: transducer_to_volume_transform: Annotated[ArrayTransform, OpenLIFUFieldData("Transducer to volume transform", "Transform output by transducer tracking algorithm to register the transducer surface to the volume")] """Transform output by transducer tracking algorithm to register the transducer surface to the volume""" - photoscan_to_volume_transform: Annotated[ArrayTransform, OpenLIFUFieldData("Photoscan to volume transform", "Transform output by the transducer tracking algorithm to register the photoscan model the volume's skin segmentation")] - """Transform output by the transducer tracking algorithm to register the photoscan model the volume's skin segmentation""" + photoscan_registration_id: Annotated[str | None, OpenLIFUFieldData("Photoscan registration ID", "ID of the PhotoscanRegistration this tracking result was computed against. May be None for legacy entries imported from sessions saved before the registration concept was split out.")] = None + """ID of the :class:`PhotoscanRegistration` this tracking result was computed against.""" + + approval: Annotated[bool, OpenLIFUFieldData("Tracking approved?", "Approval state of the transducer tracking result. `True` means the user has provided some kind of confirmation that the transform result agrees with reality.")] = False + """Approval state of the transducer tracking result. ``True`` means the user has provided + some kind of confirmation that the transform result agrees with reality.""" + + id: Annotated[str | None, OpenLIFUFieldData("Result ID", "Stable identifier for this tracking result, unique within the session. Survives reordering or deletion of other results so that downstream references (e.g. solutions) remain valid.")] = None + """Stable identifier for this tracking result, unique within the session. Survives reordering or + deletion of other results so that downstream references (e.g. solutions) remain valid.""" + + target_id: Annotated[str | None, OpenLIFUFieldData("Target ID", "ID of the openlifu Point target this tracking result was computed for, if any.")] = None + """ID of the openlifu Point target this tracking result was computed for, if any.""" + +@dataclass +class SolutionInfo: + """Per-solution provenance record stored on a :class:`Session`. + + A ``Solution`` on disk is a compact, hardware-ready object with no intrinsic knowledge of which + target, transducer pose, or protocol produced it. ``SolutionInfo`` carries that provenance at the + session level so consumers can link a loaded ``Solution`` back to its context and so the session + can cascade-delete solutions when the referenced target / virtual-fit / tracking result is removed. + + Solutions on disk that have no matching ``SolutionInfo`` entry in the owning session are + considered orphaned: they are not loaded, and they are purged from the ``solutions/`` directory + on session save. + """ + + VALID_TRANSDUCER_TRANSFORM_SOURCES: ClassVar[FrozenSet[str]] = frozenset({"virtual_fit", "localization"}) + """Allowed values for :attr:`transducer_transform_source`.""" + + solution_id: Annotated[str, OpenLIFUFieldData("Solution ID", "ID of the openlifu Solution this record describes. Corresponds to a Solution stored under the session's ``solutions/`` directory in the database.")] + """ID of the openlifu Solution this record describes.""" + + protocol_id: Annotated[str, OpenLIFUFieldData("Protocol ID", "ID of the protocol that produced this solution.")] + """ID of the protocol that produced this solution.""" + + target_id: Annotated[str, OpenLIFUFieldData("Target ID", "ID of the openlifu Point target this solution was computed for. Must match an entry in the session's ``targets`` list.")] + """ID of the openlifu Point target this solution was computed for.""" + + transducer_id: Annotated[str, OpenLIFUFieldData("Transducer ID", "ID of the transducer this solution was computed for.")] + """ID of the transducer this solution was computed for.""" - transducer_to_volume_tracking_approved: Annotated[bool, OpenLIFUFieldData("Transducer tracking approved?", "Approval state of transducer to volume tracking result. `True` means the user has provided some kind of confirmation that the transform result agrees with reality.")] = False - """Approval state of transducer to volume tracking result. `True` means the user has provided some kind of - confirmation that the transform result agrees with reality.""" + transducer_transform_source: Annotated[str, OpenLIFUFieldData("Transducer transform source", "Provenance of the transducer pose used to compute this solution. One of 'virtual_fit' or 'localization'.")] + """Provenance of the transducer pose used to compute this solution. - photoscan_to_volume_tracking_approved: Annotated[bool, OpenLIFUFieldData("Photoscan tracking approved?", "Approval state of photoscan to volume tracking result. `True` means the user has provided some kind of confirmation that the transform result agrees with reality.")] = False - """Approval state of photoscan to volume tracking result. `True` means the user has provided some kind of - confirmation that the transform result agrees with reality.""" + Must be one of :attr:`VALID_TRANSDUCER_TRANSFORM_SOURCES` (``'virtual_fit'`` or + ``'localization'``). See :attr:`transducer_transform_source_id` for the specific + virtual-fit / tracking-result identifier. + """ + + transducer_transform_source_id: Annotated[str | None, OpenLIFUFieldData("Transducer transform source ID", "Opaque identifier of the specific virtual-fit or transducer-tracking result the transducer pose came from at compute time. Interpretation is source-dependent; SlicerOpenLIFU uses ':' for VF and the TransducerTrackingResult.id for TT. Downstream consumers use this to determine whether the solution's source is still 'live' (present + approved), 'revoked' (present but no longer approved), or 'missing' (removed from the session).")] = None + """Opaque identifier of the specific virtual-fit or transducer-tracking result the + transducer pose came from at compute time. + + Interpretation is source-dependent -- see the ``transducer_transform_source`` field for + the kind. SlicerOpenLIFU uses ``":"`` for VF (composite key that stays + stable across approval-only changes) and the :class:`TransducerTrackingResult.id` for TT. + Downstream consumers can use this to determine whether the solution's source is still + live (present + approved), revoked (present but no longer approved), or missing (removed + from the session), and to warn / block sending an out-of-date solution to hardware. + ``None`` for legacy entries that predate this field; consumers should degrade gracefully + (e.g. treat as "legacy" status rather than "missing"). + """ + + approved: Annotated[bool, OpenLIFUFieldData("Approved", "Whether the user has approved this solution for sonication. Approval is tracked at the session-provenance layer because generating a Solution via Python has no approval concept; approval is a user-session-time decision.")] = False + """Whether the user has approved this solution for sonication.""" + + computed_at: Annotated[datetime | None, OpenLIFUFieldData("Computed at", "Timestamp at which the Solution was first computed. Serialized as ISO 8601. Left as None for legacy entries that predate this field.")] = None + """Timestamp at which the Solution was first computed. ``None`` for legacy entries that predate the field.""" + + array_transform: Annotated[ArrayTransform | None, OpenLIFUFieldData("Array transform", "The transducer array-to-volume transform matrix (with units) that was used when this solution was computed. Downstream consumers should reproduce this exact pose when displaying the solution so the peak-negative-pressure volume renders where it was computed, independent of subsequent virtual-fit / transducer-tracking approval changes. Left as None for legacy entries that predate this field.")] = None + """Transducer array-to-volume transform that was used when this solution was computed. + + The Solution's PNP / intensity volumes are stored in transducer-local coordinates; their world + position is ``array_transform.matrix @ local_coords``. Showing a solution should snap the + transducer to *this* matrix rather than to the currently-approved virtual-fit / transducer-tracking + result, because the currently-approved result may have changed since the solution was computed + (or, in the case of a pre-solution computed off a virtual-fit result whose approval was later + revoked, may no longer be available at all). ``None`` for legacy entries that predate this field; + consumers should fall back to their previous "snap to best current approved" behavior in that case. + """ + + def __post_init__(self): + if self.transducer_transform_source not in self.VALID_TRANSDUCER_TRANSFORM_SOURCES: + raise ValueError( + f"transducer_transform_source must be one of {sorted(self.VALID_TRANSDUCER_TRANSFORM_SOURCES)}, " + f"got {self.transducer_transform_source!r}" + ) + # Accept ISO-8601 strings for ``computed_at`` transparently so callers (like + # :meth:`Session.from_dict`) that don't pre-parse the timestamp still get a + # ``datetime`` on the resulting instance. + if isinstance(self.computed_at, str): + self.computed_at = datetime.fromisoformat(self.computed_at) + # Accept dict for ``array_transform`` transparently so callers that don't pre-parse + # (like :meth:`Session.from_dict` decoding ``solutions`` entries) still get an + # ``ArrayTransform`` on the resulting instance. + if isinstance(self.array_transform, dict): + self.array_transform = ArrayTransform.from_dict(self.array_transform) @dataclass class Session: @@ -71,6 +187,15 @@ class Session: transducer_id: Annotated[str | None, OpenLIFUFieldData("Transducer ID", "ID of the transducer associated with this session")] = None """ID of the transducer associated with this session""" + solution_id: Annotated[str, OpenLIFUFieldData("Solution ID", "ID of the most recently computed sonication Solution for this session, or '' if there is none. Cleared whenever the array_transform changes because a Solution is only valid for the transducer pose it was computed against.")] = "" + """ID of the most recently computed sonication ``Solution`` for this session, or ``""`` if there is none. + + Cleared whenever the ``array_transform`` changes (manual move, virtual fit, transducer tracking), because a + ``Solution`` is only valid for the specific transducer pose it was computed against. Consumers loading a + session can use this to fetch the persisted ``Solution`` (and its analysis) from the database rather than + re-running the simulation. + """ + array_transform: Annotated[ArrayTransform, OpenLIFUFieldData("Array transform", "The transducer affine transform matrix with units, situating the transducer in space")] = field(default_factory=lambda: ArrayTransform(np.eye(4), "mm")) """The transducer affine transform matrix with units, situating the transducer in space""" @@ -80,6 +205,18 @@ class Session: markers: Annotated[List[Point], OpenLIFUFieldData("Markers", "Registration markers saved to this session")] = field(default_factory=list) """Registration markers saved to this session""" + photoscans: Annotated[List[str], OpenLIFUFieldData("Photoscan IDs", "IDs of photoscans that belong to this session. Each ID corresponds to a Photoscan stored under the session's ``photoscans/`` directory in the database.")] = field(default_factory=list) + """IDs of photoscans that belong to this session. Each ID corresponds to a Photoscan + stored under the session's ``photoscans/`` directory in the database. This is the + authoritative list used to decide which photoscans to keep on save; legacy sessions + that omit this field are auto-populated from the on-disk index on load.""" + + photocollections: Annotated[List[str], OpenLIFUFieldData("Photocollection reference numbers", "Reference numbers of photocollections that belong to this session. Each entry corresponds to a directory under the session's ``photocollections/`` directory in the database.")] = field(default_factory=list) + """Reference numbers (scan IDs) of photocollections that belong to this session. + Each entry corresponds to a directory under the session's ``photocollections/`` + directory in the database. Legacy sessions that omit this field are auto-populated + from the on-disk index on load.""" + attrs: Annotated[dict, OpenLIFUFieldData("Custom attributes", "Dictionary of additional custom attributes to save to the session")] = field(default_factory=dict) """Dictionary of additional custom attributes to save to the session""" @@ -98,6 +235,20 @@ class Session: transducer_tracking_results: Annotated[List[TransducerTrackingResult], OpenLIFUFieldData("Tracking results", "List of any transducer tracking results")] = field(default_factory=list) """List of any transducer tracking results""" + photoscan_registrations: Annotated[List[PhotoscanRegistration], OpenLIFUFieldData("Photoscan registrations", "List of photoscan-to-volume registrations stored on this session.")] = field(default_factory=list) + """List of photoscan-to-volume registrations stored on this session. Each transducer tracking + result references one of these registrations by ``photoscan_registration_id``.""" + + solutions: Annotated[List[SolutionInfo], OpenLIFUFieldData("Solutions", "Per-solution provenance records: one entry per Solution belonging to this session, carrying the target id, transducer id, protocol id, and transducer-transform source ('virtual_fit' or 'localization'). Solutions on disk that have no matching entry here are considered orphaned and are purged on save.")] = field(default_factory=list) + """Per-solution provenance records. See :class:`SolutionInfo`. + + Every :class:`~openlifu.plan.Solution` that belongs to this session should have an entry here. + On-disk solutions without a matching entry are orphaned and are purged on save. Consumers use + this list to link a loaded ``Solution`` back to its target / transducer / protocol / pose + provenance and to cascade-delete solutions when a referenced target, virtual fit, or tracking + result is removed. + """ + def __post_init__(self): if self.id is None and self.name is None: self.id = "session" @@ -141,17 +292,62 @@ def from_dict(d:Dict): raise ValueError("Sessions no longer recognize a volume attribute -- it is now volume_id.") if 'array_transform' in d: d['array_transform'] = ArrayTransform.from_dict(d['array_transform']) + + # PhotoscanRegistrations are split out of TT results as of the multi-registration refactor. + # Old session JSONs lack this key; if absent we start with an empty list and may populate it + # below when migrating legacy TT entries that still carry an embedded photoscan_to_volume_transform. + if 'photoscan_registrations' in d: + d['photoscan_registrations'] = [ + PhotoscanRegistration( + photoscan_id=p['photoscan_id'], + transform=ArrayTransform.from_dict(p['transform']), + approval=p.get('approval', False), + id=p.get('id'), + ) + for p in d['photoscan_registrations'] + ] + else: + d['photoscan_registrations'] = [] + if 'transducer_tracking_results' in d: - d['transducer_tracking_results'] = [ - TransducerTrackingResult( - t['photoscan_id'], - ArrayTransform.from_dict(t['transducer_to_volume_transform']), - ArrayTransform.from_dict(t['photoscan_to_volume_transform']), - t['transducer_to_volume_tracking_approved'], - t['photoscan_to_volume_tracking_approved'] - ) - for t in d['transducer_tracking_results'] - ] + # Per-photoscan counter for any registrations we synthesize during legacy migration; + # continues past the count of registrations already present so we don't collide. + pr_count_by_photoscan: Dict[str, int] = {} + for pr in d['photoscan_registrations']: + pr_count_by_photoscan[pr.photoscan_id] = pr_count_by_photoscan.get(pr.photoscan_id, 0) + 1 + + migrated_tt: List[TransducerTrackingResult] = [] + for t in d['transducer_tracking_results']: + if 'photoscan_to_volume_transform' in t: + # Legacy entry: split the embedded PV transform out into its own registration. + pid = t['photoscan_id'] + n = pr_count_by_photoscan.get(pid, 0) + pr_count_by_photoscan[pid] = n + 1 + synthesized_pr_id = f"{pid}__pr__{n:02d}" + d['photoscan_registrations'].append(PhotoscanRegistration( + photoscan_id=pid, + transform=ArrayTransform.from_dict(t['photoscan_to_volume_transform']), + approval=t.get('photoscan_to_volume_tracking_approved', False), + id=synthesized_pr_id, + )) + migrated_tt.append(TransducerTrackingResult( + photoscan_id=pid, + transducer_to_volume_transform=ArrayTransform.from_dict(t['transducer_to_volume_transform']), + photoscan_registration_id=synthesized_pr_id, + approval=t.get('transducer_to_volume_tracking_approved', t.get('approval', False)), + id=t.get('id'), + target_id=t.get('target_id'), + )) + else: + migrated_tt.append(TransducerTrackingResult( + photoscan_id=t['photoscan_id'], + transducer_to_volume_transform=ArrayTransform.from_dict(t['transducer_to_volume_transform']), + photoscan_registration_id=t.get('photoscan_registration_id'), + approval=t.get('approval', t.get('transducer_to_volume_tracking_approved', False)), + id=t.get('id'), + target_id=t.get('target_id'), + )) + d['transducer_tracking_results'] = migrated_tt if isinstance(d['targets'], list): if len(d['targets'])>0 and isinstance(d['targets'][0], dict): d['targets'] = [Point.from_dict(p) for p in d['targets']] @@ -171,6 +367,8 @@ def from_dict(d:Dict): d['markers'] = [Point.from_dict(d['markers'])] elif isinstance(d['markers'], Point): d['markers'] = [d['markers']] + if 'solutions' in d: + d['solutions'] = [SolutionInfo(**s) for s in d['solutions']] return Session(**d) def to_dict(self): @@ -192,6 +390,8 @@ def to_dict(self): ] d['transducer_tracking_results'] = [asdict(t) for t in d['transducer_tracking_results']] + d['photoscan_registrations'] = [asdict(r) for r in d['photoscan_registrations']] + d['solutions'] = [asdict(s) for s in d['solutions']] return d diff --git a/src/openlifu/db/sonication_session.py b/src/openlifu/db/sonication_session.py new file mode 100644 index 00000000..cc8ce7f9 --- /dev/null +++ b/src/openlifu/db/sonication_session.py @@ -0,0 +1,177 @@ +from __future__ import annotations + +import copy +import json +from dataclasses import asdict, dataclass, field +from datetime import datetime +from pathlib import Path +from typing import Annotated, Dict, List + +from openlifu.db.session import ( + PhotoscanRegistration, + SolutionInfo, + TransducerTrackingResult, +) +from openlifu.geo.transforms import ArrayTransform +from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.json import PYFUSEncoder +from openlifu.util.strings import sanitize + + +@dataclass +class SonicationSession: + """At-treatment-time session. + + References a :class:`Plan` by id (immutable input); the Plan's ``target``, + ``volume_id``, ``protocol_id``, ``transducer_id``, and ``array_transform`` are the + reference "what we're trying to hit". The SonicationSession itself owns "what we + actually did": photoscan registrations, transducer-tracking results, the final + :class:`~openlifu.plan.Solution` reference, and the list of :class:`~openlifu.plan.Run` + ids performed against this session. + + Unlike :class:`PlanningSession` there is at most one solution per SonicationSession + (``solution: Optional[SolutionInfo]``, not a list). Recomputing replaces it in + place. Multi-solution generation and selection is a plausible future feature but + is deliberately out of scope for the initial split-session refactor. + + Photoscans are physically stored at subject scope + (``subjects/{sid}/photoscans/{pid}/``) but logically owned by this session via + :attr:`photoscan_ids`. Two SonicationSessions on the same subject each own their + own list. + """ + + id: Annotated[str | None, OpenLIFUFieldData("Sonication session ID", "ID of this sonication session")] = None + """ID of this sonication session""" + + name: Annotated[str | None, OpenLIFUFieldData("Sonication session name", "Sonication session name")] = None + """Sonication session name""" + + subject_id: Annotated[str | None, OpenLIFUFieldData("Subject ID", "ID of the parent subject of this sonication session")] = None + """ID of the parent subject of this sonication session""" + + plan_id: Annotated[str | None, OpenLIFUFieldData("Plan ID", "ID of the Plan this sonication session was launched from. Immutable frozen reference; volume/target/protocol/transducer/array_transform come from the Plan.")] = None + """ID of the :class:`Plan` this sonication session was launched from. Immutable + frozen reference; ``volume_id``, ``target``, ``protocol_id``, ``transducer_id``, + and ``array_transform`` all come from the Plan.""" + + date_created: Annotated[datetime, OpenLIFUFieldData("Date created", "Date of creation of the sonication session")] = field(default_factory=datetime.now) + """Date of creation of the sonication session""" + + date_modified: Annotated[datetime, OpenLIFUFieldData("Date modified", "Date of last modification of the sonication session")] = field(default_factory=datetime.now) + """Date of last modification of the sonication session""" + + photoscan_ids: Annotated[List[str], OpenLIFUFieldData("Photoscan IDs", "IDs of photoscans owned by this sonication session. Photoscan files live at subject scope; a photoscan_id in this list is the assertion that this session owns that photoscan.")] = field(default_factory=list) + """IDs of photoscans owned by this sonication session. + + Photoscan files live at subject scope (``subjects/{sid}/photoscans/{pid}/``); + a photoscan_id in this list is the assertion that this session owns that + photoscan. Two SonicationSessions on the same subject each own their own list. + """ + + photoscan_registrations: Annotated[List[PhotoscanRegistration], OpenLIFUFieldData("Photoscan registrations", "List of photoscan-to-volume registrations stored on this sonication session. Each transducer tracking result references one of these registrations by ``photoscan_registration_id``.")] = field(default_factory=list) + """List of photoscan-to-volume registrations stored on this sonication session. + Each transducer tracking result references one of these registrations by + ``photoscan_registration_id``.""" + + transducer_tracking_results: Annotated[List[TransducerTrackingResult], OpenLIFUFieldData("Transducer tracking results", "List of transducer tracking results computed on this sonication session")] = field(default_factory=list) + """List of transducer tracking results computed on this sonication session""" + + solution: Annotated[SolutionInfo | None, OpenLIFUFieldData("Solution", "SolutionInfo reference to the ONE final solution computed for this sonication session, or None. Actual Solution files live at subject scope. Recomputing replaces this in place; the run_ids history captures what was actually delivered under previous solutions.")] = None + """:class:`SolutionInfo` reference to the ONE final solution computed for this + sonication session, or ``None``. + + Actual :class:`~openlifu.plan.Solution` files live at subject scope + (``subjects/{sid}/solutions/``). Recomputing replaces this in place; the + :attr:`run_ids` history captures what was actually delivered under previous + solutions. Multi-solution generation and selection is a plausible future + feature but is deliberately out of scope for the initial split-session refactor. + """ + + run_ids: Annotated[List[str], OpenLIFUFieldData("Run IDs", "IDs of Runs performed against this sonication session, in chronological order. Actual Run JSON files live under the session's runs/ subdirectory.")] = field(default_factory=list) + """IDs of :class:`~openlifu.plan.Run`\\ s performed against this sonication session, + in chronological order. Actual Run JSON files live under the session's ``runs/`` + subdirectory.""" + + attrs: Annotated[dict, OpenLIFUFieldData("Custom attributes", "Dictionary of additional custom attributes to save to the sonication session")] = field(default_factory=dict) + """Dictionary of additional custom attributes to save to the sonication session""" + + def __post_init__(self): + if self.id is None and self.name is None: + self.id = "sonication_session" + if self.id is None: + self.id = sanitize(self.name, "snake") + if self.name is None: + self.name = self.id + + @staticmethod + def from_file(filename) -> SonicationSession: + """Load a SonicationSession from a JSON file.""" + with open(filename) as f: + return SonicationSession.from_dict(json.load(f)) + + @staticmethod + def from_dict(d: Dict) -> SonicationSession: + """Reconstruct a SonicationSession from its dictionary representation.""" + d = dict(d) # shallow copy; we mutate below + if "date_created" in d and isinstance(d["date_created"], str): + d["date_created"] = datetime.fromisoformat(d["date_created"]) + if "date_modified" in d and isinstance(d["date_modified"], str): + d["date_modified"] = datetime.fromisoformat(d["date_modified"]) + if "photoscan_registrations" in d: + d["photoscan_registrations"] = [ + pr if isinstance(pr, PhotoscanRegistration) else PhotoscanRegistration( + photoscan_id=pr["photoscan_id"], + transform=ArrayTransform.from_dict(pr["transform"]), + approval=pr.get("approval", False), + id=pr.get("id"), + ) + for pr in d["photoscan_registrations"] + ] + if "transducer_tracking_results" in d: + d["transducer_tracking_results"] = [ + t if isinstance(t, TransducerTrackingResult) else TransducerTrackingResult( + photoscan_id=t["photoscan_id"], + transducer_to_volume_transform=ArrayTransform.from_dict(t["transducer_to_volume_transform"]), + photoscan_registration_id=t.get("photoscan_registration_id"), + approval=t.get("approval", False), + id=t.get("id"), + target_id=t.get("target_id"), + ) + for t in d["transducer_tracking_results"] + ] + if "solution" in d and d["solution"] is not None and not isinstance(d["solution"], SolutionInfo): + d["solution"] = SolutionInfo(**d["solution"]) + return SonicationSession(**d) + + def to_dict(self) -> Dict: + """Serialize the SonicationSession to a dictionary.""" + d = copy.deepcopy(self.__dict__) + d["date_created"] = d["date_created"].isoformat() + d["date_modified"] = d["date_modified"].isoformat() + d["photoscan_registrations"] = [asdict(pr) for pr in d["photoscan_registrations"]] + d["transducer_tracking_results"] = [asdict(t) for t in d["transducer_tracking_results"]] + d["solution"] = asdict(d["solution"]) if d["solution"] is not None else None + return d + + @staticmethod + def from_json(json_string: str) -> SonicationSession: + """Load a SonicationSession from a JSON string.""" + return SonicationSession.from_dict(json.loads(json_string)) + + def to_json(self, compact: bool) -> str: + """Serialize a SonicationSession to a JSON string.""" + if compact: + return json.dumps(self.to_dict(), separators=(",", ":"), cls=PYFUSEncoder) + return json.dumps(self.to_dict(), indent=4, cls=PYFUSEncoder) + + def to_file(self, filename) -> None: + """Write the SonicationSession to a JSON file, creating parent directories as needed.""" + Path(filename).parent.parent.mkdir(exist_ok=True) # sonication_sessions directory + Path(filename).parent.mkdir(exist_ok=True) # {sonication_session_id} directory + with open(filename, "w") as f: + f.write(self.to_json(compact=False)) + + def update_modified_time(self, time: datetime | None = None) -> None: + if time is None: + time = datetime.now() + self.date_modified = time diff --git a/src/openlifu/plan/__init__.py b/src/openlifu/plan/__init__.py index d3189b3d..48f49515 100644 --- a/src/openlifu/plan/__init__.py +++ b/src/openlifu/plan/__init__.py @@ -14,5 +14,5 @@ "SolutionAnalysis", "SolutionAnalysisOptions", "TargetConstraints", - "ParameterConstraint", + "ParameterConstraint" ] diff --git a/src/openlifu/plan/param_constraint.py b/src/openlifu/plan/param_constraint.py index 0364b9d7..fe08d163 100644 --- a/src/openlifu/plan/param_constraint.py +++ b/src/openlifu/plan/param_constraint.py @@ -10,8 +10,8 @@ PARAM_STATUS_SYMBOLS = { "ok": "✅", - "warning": "❗", - "error": "❌" + "warning": "⚠️", + "error": "⛔", } @dataclass diff --git a/src/openlifu/plan/protocol.py b/src/openlifu/plan/protocol.py index eb47df3f..391a03e0 100644 --- a/src/openlifu/plan/protocol.py +++ b/src/openlifu/plan/protocol.py @@ -328,6 +328,7 @@ def calc_solution( apodizations_to_stack: List[np.ndarray] = [] simulation_result_aggregated: xa.Dataset = xa.Dataset() foci: List[Point] = self.focal_pattern.get_targets(target) + order = self.focal_pattern.get_order() # updating solution sequence if pulse mismatch if (self.sequence.pulse_count % len(foci)) != 0: @@ -357,9 +358,9 @@ def calc_solution( voltage=voltage, sequence=self.sequence, foci=foci, + order=order, target=target, simulation_result=xa.Dataset(), - approved=False, description= ( f"A solution computed for the {self.name} protocol with transducer {transducer.name}" f" for target {target.id}." @@ -387,7 +388,7 @@ def calc_solution( pnp_aggregated = solution.simulation_result['p_min'].max(dim="focal_point_index", keep_attrs=True) ppp_aggregated = solution.simulation_result['p_max'].max(dim="focal_point_index", keep_attrs=True) # 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['intensity'].mean(dim="focal_point_index", keep_attrs=True) + intensity_aggregated = solution.get_ita(solution.simulation_result['intensity']) 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 diff --git a/src/openlifu/plan/solution.py b/src/openlifu/plan/solution.py index 6077a953..5cdb1881 100644 --- a/src/openlifu/plan/solution.py +++ b/src/openlifu/plan/solution.py @@ -7,7 +7,7 @@ from dataclasses import asdict, dataclass, field from datetime import datetime from pathlib import Path -from typing import Annotated, Dict, List, Tuple +from typing import Annotated, Dict, List, Literal, Tuple import numpy as np import xarray as xa @@ -25,6 +25,11 @@ model_tx_temperature_rise, ) from openlifu.sim import SimSetup, run_simulation +from openlifu.sim.thermal import ( + compute_heat_source, + generate_pulse_events, + run_thermal_simulation, +) from openlifu.util.annotations import OpenLIFUFieldData from openlifu.util.checkgpu import gpu_available from openlifu.util.json import PYFUSEncoder @@ -85,6 +90,9 @@ class Solution: what determines how many times each point will be used. """ + order: Annotated[list[int] | None, OpenLIFUFieldData("Focus order", "Order of Foci (1-indexed) in the sequence")] = None + """Order of Foci (1-indexed) in the sequence. This is a list of integers that specifies the order in which the foci are used in the pulse sequence. If None, the foci are used in the order they are listed in the `foci` attribute.""" + # there was "target_id" in the matlab software, but here we do not have the concept of a target ID. # I believe this was only needed in the matlab software because solutions were organized by target rather # than having their own unique solution ID. We do have unique solution IDs so it's possible we don't need @@ -98,10 +106,6 @@ class Solution: simulation_result: Annotated[xa.Dataset, OpenLIFUFieldData("Simulation result", "The xarray Dataset of simulation results")] = field(default_factory=xa.Dataset) """The xarray Dataset of simulation results""" - approved: Annotated[bool, OpenLIFUFieldData("Approved?", "Approval state of this solution as a sonication plan. `True` means the user has provided some kind of confirmation that the solution is safe and acceptable to be executed.")] = False - """Approval state of this solution as a sonication plan. `True` means the user has provided some - kind of confirmation that the solution is safe and acceptable to be executed.""" - def __post_init__(self): self.logger = logging.getLogger(__name__) if self.delays is not None: @@ -192,6 +196,220 @@ def simulate(self, dim='focal_point_index', ) + def simulate_thermal( + self, + params: xa.Dataset, + acoustic_result: xa.Dataset | None = None, + pressure_field: str = 'p_min', + alpha_power: float = 0.9, + dt: float = 0.1, + t_end: float | None = None, + oversample: int = 1, + T0: float = 0.0, + T_init: np.ndarray | None = None, + mode: Literal['direct', 'superpose', 'impulse'] = 'impulse', + ) -> xa.Dataset: + """Run a thermal diffusion simulation for this Solution. + + Uses the acoustic pressure field (from a prior call to :meth:`simulate` + or provided via ``acoustic_result``) together with tissue parameter + maps to compute the volumetric heat deposition rate at each focus, then + integrates the heterogeneous heat equation over the pulse sequence. + + For typical LIFU parameters (pulses of 1-100 ms with sequences of + seconds to minutes), the default ``mode='impulse'`` marches a single + temperature field forward in time and applies each pulse as an + instantaneous energy deposit. Pulse events do NOT need to be aligned + with the output-``dt`` grid: fractional sub-steps are inserted at + every event, so multiple pulses per output step -- including pulses + that would fall between samples -- all deposit their energy exactly. + This is dramatically faster than resolving pulse edges directly and, + for ``pulse_duration`` much smaller than the thermal timescale, is + essentially exact. + + .. warning:: The default openlifu ``Material`` definitions in + :mod:`openlifu.seg.material` have ``attenuation=0`` for water, + tissue, skull, and air. With those defaults the heat deposition + ``Q = alpha * p^2 / (rho * c)`` is zero and no temperature rise + will be simulated. Pass a segmentation method whose materials use + realistic ``attenuation`` values (e.g. 0.3-0.5 dB/cm/MHz for + tissue, ~4 dB/cm/MHz for skull) if you want to see thermal effects. + + Args: + params: Tissue parameter maps (an xarray Dataset). Must contain + ``density``, ``sound_speed``, ``attenuation``, ``specific_heat``, + and ``thermal_conductivity`` variables on the simulation grid. + Typically produced by + :meth:`openlifu.seg.SegmentationMethod.seg_params` or + :meth:`openlifu.seg.SegmentationMethod.ref_params`. + acoustic_result: The acoustic simulation Dataset with a + ``focal_point_index`` dimension. If None, uses + ``self.simulation_result``. + pressure_field: Name of the pressure variable in + ``acoustic_result`` to use for heat-source calculation. + Default ``'p_min'`` (peak-negative-pressure amplitude). + alpha_power: Power law exponent for attenuation. Default 0.9. + dt: Output time step (s). Default 0.1. + t_end: Simulation duration (s). If None, defaults to the full + sequence duration (as reported by ``self.sequence``). + oversample: Number of internal sub-steps per output step for the + explicit diffusion solver. + T0: Background temperature offset (degC). Only affects the + reported peak temperature; the stored ``temperature_rise`` is + always the delta from ambient. + T_init: Optional initial temperature-rise field. + mode: ``'impulse'`` (default) treats each pulse as an + instantaneous energy deposit and superposes per-focus + impulse responses -- ideal for pulses much shorter than the + thermal diffusion timescale. ``'superpose'`` precomputes a + single-pulse response with the source held on for + ``pulse_duration`` and superposes shifted copies. + ``'direct'`` walks the sequence event by event, matching the + LOFUgen1 MATLAB implementation. + + Returns: + An xarray Dataset with ``temperature_rise`` (dims ``t, x, y, z``), + ``Tmax``, and ``sim_time`` variables. See + :func:`openlifu.sim.run_thermal_simulation` for details. + """ + if acoustic_result is None: + if self.simulation_result is None or len(self.simulation_result) == 0: + raise ValueError( + "No acoustic_result provided and no simulation_result found on the Solution. " + "Call `simulate()` first or pass `acoustic_result` explicitly." + ) + acoustic_result = self.simulation_result + if pressure_field not in acoustic_result: + raise ValueError( + f"Pressure field '{pressure_field}' not found in acoustic_result. " + f"Available: {list(acoustic_result.data_vars)}" + ) + + pressure_da = acoustic_result[pressure_field] + if 'focal_point_index' not in pressure_da.dims: + raise ValueError( + f"acoustic_result['{pressure_field}'] must have a 'focal_point_index' dim." + ) + pressure_da = pressure_da.transpose('focal_point_index', 'x', 'y', 'z') + + Q_data = compute_heat_source( + params=params, + pressure=pressure_da.data, + freq=self.pulse.frequency, + alpha_power=alpha_power, + ) + heat_sources = xa.DataArray( + Q_data, + dims=['focal_point_index', 'x', 'y', 'z'], + coords={ + 'focal_point_index': pressure_da.coords['focal_point_index'], + 'x': pressure_da.coords['x'], + 'y': pressure_da.coords['y'], + 'z': pressure_da.coords['z'], + }, + attrs={'units': 'W/m^3', 'long_name': 'Volumetric heat deposition rate'}, + ) + + events = generate_pulse_events( + pulse_count=self.sequence.pulse_count, + pulse_interval=self.sequence.pulse_interval, + pulse_train_count=self.sequence.pulse_train_count, + pulse_train_interval=self.sequence.pulse_train_interval, + pulse_duration=self.pulse.duration, + num_foci=self.num_foci(), + order=self.order, + ) + + if t_end is None: + t_end = self.sequence.get_sequence_duration() + + return run_thermal_simulation( + heat_sources=heat_sources, + params=params, + events=events, + dt=dt, + t_end=t_end, + oversample=oversample, + T0=T0, + T_init=T_init, + mode=mode, + pulse_duration=self.pulse.duration, + ) + + def get_mainlobe_mask(self, + simulation_result: xa.Dataset | None = None, + options: SolutionAnalysisOptions = SolutionAnalysisOptions(), + units: str | None = None + ) -> xa.DataArray: + """Get a masked version of the simulation result where only the mainlobe is unmasked. + """ + return self.get_mask(simulation_result=simulation_result, type='mainlobe', options=options, units=units) + + def get_sidelobe_mask(self, + simulation_result: xa.Dataset | None = None, + options: SolutionAnalysisOptions = SolutionAnalysisOptions(), + units: str | None = None + ) -> xa.DataArray: + """Get a masked version of the simulation result where only the sidelobe is unmasked. + """ + return self.get_mask(simulation_result=simulation_result, type='sidelobe', options=options, units=units) + + def get_mask(self, + simulation_result: xa.Dataset | None = None, + type='mainlobe', + options: SolutionAnalysisOptions = SolutionAnalysisOptions(), + units: str | None = None) -> xa.DataArray: + """Get a masked version of the simulation result where only the mainlobe is unmasked. + """ + if simulation_result is None: + if self.simulation_result is None or len(self.simulation_result)==0: + raise ValueError("No simulation result provided for masking, and no simulation result found in the Solution.") + simulation_result = self.simulation_result + simulation_result_scaled = rescale_coords(simulation_result, options.distance_units) + masks = [] + wavelength = options.ref_sound_speed / self.pulse.frequency + for focus_index in range(self.num_foci()): + focus = self.foci[focus_index].get_position(units=options.distance_units) + apodization = self.apodizations[focus_index] + origin = self.transducer.get_effective_origin(apodizations=apodization, units=options.distance_units) + aperture_dia = self.transducer.get_effective_aperture_radius(apodizations=apodization, units=options.distance_units)*2 + focal_dist = np.linalg.norm(np.array(focus) - np.array(origin)) + fnum = np.max([focal_dist / aperture_dia, 1.0]) + if type == 'mainlobe': + if options.mainlobe_radius is None: + distance = fnum * wavelength + else: + distance = options.mainlobe_radius + aspect_ratio = options.mainlobe_aspect_ratio + operator = '<' + elif type == 'sidelobe': + if options.sidelobe_radius is None: + distance = 2.5 * fnum * wavelength + else: + distance = options.sidelobe_radius + aspect_ratio = options.mainlobe_aspect_ratio + operator = '>' + else: + raise ValueError(f"Invalid mask type {type}. Must be 'mainlobe' or 'sidelobe'.") + mask = get_mask( + simulation_result_scaled.isel(focal_point_index=focus_index), + focus = focus, + origin = origin, + distance = distance, + operator = operator, + aspect_ratio = aspect_ratio + ) + if type == 'sidelobe' and options.sidelobe_zmin is not None: + z_mask = simulation_result_scaled.isel(focal_point_index=focus_index).z > options.sidelobe_zmin + mask = mask.where(z_mask, False) + masks.append(mask) + masks = xa.concat(masks, dim='focal_point_index') + if units is None and 'units' in simulation_result.coords['x'].attrs: + units = simulation_result.coords['x'].attrs['units'] + if units is not None: + masks = rescale_coords(masks, units) + return masks + def analyze(self, simulation_result: xa.Dataset | None = None, options: SolutionAnalysisOptions = SolutionAnalysisOptions(), @@ -211,40 +429,39 @@ def analyze(self, solution_analysis = SolutionAnalysis() dt = 1 / (self.pulse.frequency * 20) - t = self.pulse.calc_time(dt) - input_signal_V = self.pulse.calc_pulse(t) * self.voltage if simulation_result is None: if self.simulation_result is None or len(self.simulation_result)==0: raise ValueError("No simulation result provided for analysis, and no simulation result found in the Solution.") simulation_result = self.simulation_result - pnp_MPa_all = rescale_data_arr(rescale_coords(simulation_result['p_min'], options.distance_units),"MPa") - ipa_Wcm2_all = rescale_data_arr(rescale_coords(simulation_result['intensity'], options.distance_units), "W/cm^2") + simulation_result_scaled = rescale_coords(simulation_result, options.distance_units) + pnp_MPa_all = rescale_data_arr(simulation_result_scaled['p_min'],"MPa") + ipa_Wcm2_all = rescale_data_arr(simulation_result_scaled['intensity'], "W/cm^2") if options.sidelobe_radius is np.nan: options.sidelobe_radius = options.mainlobe_radius - standoff_Z = options.standoff_density * 1500 c_tic = 40e-3 # W cm-1 A_cm = self.transducer.get_area("cm") d_eq_cm = np.sqrt(4*A_cm / np.pi) ele_sizes_cm2 = np.array([elem.get_area("cm") for elem in self.transducer.elements]) - # xyz = np.stack(np.meshgrid(*coords, indexing="xy"), axis=-1) #TODO: if fus.Axis is defined, coords.ndgrid(dim="z") - # z_mask = xyz[..., -1] >= options.sidelobe_zmin #TODO: probably wrong here, should be z{1}>=options.sidelobe_zmin; - solution_analysis.duty_cycle_pulse_train_pct = self.get_pulsetrain_dutycycle()*100 solution_analysis.duty_cycle_sequence_pct = self.get_sequence_dutycycle()*100 if self.sequence.pulse_train_interval == 0: solution_analysis.sequence_duration_s = float(self.sequence.pulse_interval * self.sequence.pulse_count * self.sequence.pulse_train_count) else: solution_analysis.sequence_duration_s = float(self.sequence.pulse_train_interval * self.sequence.pulse_train_count) - ita_mWcm2 = rescale_coords(self.get_ita(intensity=simulation_result['intensity'], units="mW/cm^2"), options.distance_units) + ita_mWcm2 = self.get_ita(intensity=simulation_result_scaled['intensity'], units="mW/cm^2") power_W = np.zeros(self.num_foci()) TIC = np.zeros(self.num_foci()) + + mainlobe_masks = self.get_mainlobe_mask(simulation_result=simulation_result_scaled, options=options) + sidelobe_masks = self.get_sidelobe_mask(simulation_result=simulation_result_scaled, options=options) + for focus_index in range(self.num_foci()): pnp_MPa = pnp_MPa_all.isel(focal_point_index=focus_index) ipa_Wcm2 = ipa_Wcm2_all.isel(focal_point_index=focus_index) @@ -265,30 +482,16 @@ def analyze(self, amplitude=self.pulse.amplitude * self.voltage, ), axis=1) - mainlobe_mask = get_mask( - pnp_MPa, - focus = focus, - origin = origin, - distance = options.mainlobe_radius, - operator = '<', - aspect_ratio = options.mainlobe_aspect_ratio - ) - - sidelobe_mask = get_mask( - pnp_MPa, - focus = focus, - origin = origin, - distance = options.sidelobe_radius, - operator = '>', - aspect_ratio=options.mainlobe_aspect_ratio - ) - z_dim = pnp_MPa.dims[2] - z_mask = pnp_MPa.coords[z_dim] > options.sidelobe_zmin - sidelobe_mask = sidelobe_mask.where(z_mask, False) + mainlobe_mask = mainlobe_masks.isel(focal_point_index=focus_index) + sidelobe_mask = sidelobe_masks.isel(focal_point_index=focus_index) + z_mask = pnp_MPa.z > options.sidelobe_zmin + sidelobe_mask = sidelobe_mask.where(z_mask, other=False) pnp_mainlobe = pnp_MPa.where(mainlobe_mask) + ipa_mainlobe = ipa_Wcm2.where(mainlobe_mask) pk = float(pnp_mainlobe.max()) - mainlobe_focus = find_centroid(pnp_mainlobe, pk*10**(-3/20), "mm") + ipk = float(ipa_mainlobe.max()) + mainlobe_focus = find_centroid(ipa_mainlobe, ipk*10**(-3/20), "mm") solution_analysis.focal_centroid_lat_mm += [mainlobe_focus[0]] solution_analysis.focal_centroid_ele_mm += [mainlobe_focus[1]] solution_analysis.focal_centroid_ax_mm += [mainlobe_focus[2]] @@ -296,6 +499,15 @@ def analyze(self, solution_analysis.mainlobe_pnp_MPa += [pk] solution_analysis.focal_gain += [pk*1e6/np.max(p0_Pa)] + if options.beamwidth_radius is None: + aperture_dia = self.transducer.get_effective_aperture_radius(apodizations=apodization, units=options.distance_units)*2 + focal_dist = np.linalg.norm(np.array(focus) - np.array(origin)) + fnum = np.max([focal_dist / aperture_dia, 1.0]) + wavelength = options.ref_sound_speed / self.pulse.frequency + beamwidth_radius = 2.0 * fnum * wavelength + else: + beamwidth_radius = options.beamwidth_radius + for dim, named_dim, scale in zip(pnp_MPa.dims, ("lat","ele","ax"), options.mainlobe_aspect_ratio): for threshdB in [3, 6]: attr_name = f'beamwidth_{named_dim}_{threshdB}dB_mm' @@ -307,8 +519,8 @@ def analyze(self, dim=dim, cutoff=cutoff, origin=origin, - min_offset=-scale*options.beamwidth_radius, - max_offset=scale*options.beamwidth_radius) + min_offset=-scale*beamwidth_radius, + max_offset=scale*beamwidth_radius) bw = getunitconversion(options.distance_units, "mm") * bw bw = [*bw0, bw] setattr(solution_analysis, attr_name, bw) @@ -345,7 +557,7 @@ def analyze(self, power_W[focus_index] = np.mean(np.sum(i0ta_Wcm2 * ele_sizes_cm2 * self.apodizations[focus_index, :])) TIC[focus_index] = power_W[focus_index] / (d_eq_cm * c_tic) solution_analysis.p0_MPa += [1e-6*np.max(p0_Pa)] - solution_analysis.global_ispta_mWcm2 = float((ita_mWcm2*z_mask).max()) + solution_analysis.global_ispta_mWcm2 = float(ita_mWcm2.where(z_mask).max()) solution_analysis.MI = (np.max(solution_analysis.mainlobe_pnp_MPa)/np.sqrt(self.pulse.frequency*1e-6)) solution_analysis.TIC = np.mean(TIC) solution_analysis.voltage_V = self.voltage @@ -491,14 +703,13 @@ def get_ita(self, intensity: xa.DataArray | None = None, units: str = "mW/cm^2") intensity_scaled = rescale_data_arr(self.simulation_result['intensity'], units) sequence_dutycycle = self.get_sequence_dutycycle() pulse_seq = (np.arange(self.sequence.pulse_count) - 1) % self.num_foci() + 1 - counts = np.zeros((1, 1, 1, self.num_foci())) - for i in range(self.num_foci()): - counts[0, 0, 0, i] = np.sum(pulse_seq == (i+1)) - intensity = intensity_scaled.copy(deep=True) - isppa_avg = np.sum(np.expand_dims(intensity.data, axis=-1) * counts, axis=-1) / np.sum(counts) - intensity.data = isppa_avg * sequence_dutycycle - - return intensity + if self.order is not None: + pulse_seq = np.array(list(self.order)*(self.sequence.pulse_count//len(self.order))) + else: + pulse_seq = np.arange(self.sequence.pulse_count) % self.num_foci() + 1 + counts = [np.sum(pulse_seq == (i+1)) for i in range(self.num_foci())] + counts_da = xa.DataArray(counts, dims=['focal_point_index'], coords={'focal_point_index': np.arange(self.num_foci())}) + return ((intensity_scaled * counts_da).sum(dim='focal_point_index') / counts_da.sum(dim='focal_point_index')) * sequence_dutycycle def to_dict(self, include_simulation_data: bool = False) -> dict: """Serialize a Solution to a dictionary @@ -586,6 +797,11 @@ def from_dict(solution_dict: dict) -> Solution: engine='scipy', ) + # Backward-compat: pre-migration Solution JSONs carried an ``approved`` field. + # Approval has moved to the Session's ``SolutionInfo`` records, so drop any + # legacy ``approved`` key here silently to keep old files loadable. + solution_dict.pop("approved", None) + return Solution(**solution_dict) diff --git a/src/openlifu/plan/solution_analysis.py b/src/openlifu/plan/solution_analysis.py index d3e8a680..7e3d3b0a 100644 --- a/src/openlifu/plan/solution_analysis.py +++ b/src/openlifu/plan/solution_analysis.py @@ -17,11 +17,12 @@ logger = logging.getLogger(__name__) DEFAULT_ORIGIN = np.zeros(3) +DEFAULT_SIDELOBE_ZMIN_MM = 10.0 PARAM_FORMATS = { "mainlobe_pnp_MPa": ["max", "0.3f", "MPa", "Mainlobe Peak Negative Pressure"], "mainlobe_isppa_Wcm2": ["max", "0.1f", "W/cm^2", "Mainlobe I_SPPA"], - "mainlobe_ispta_mWcm2": ["mean", "0.1f", "mW/cm^2", "Mainlobe I_SPTA"], + "mainlobe_ispta_mWcm2": ["max", "0.1f", "mW/cm^2", "Mainlobe I_SPTA"], "target_position_lat_mm": ["mean", "0.1f", "mm", "Target Position (Lateral)"], "target_position_ele_mm": ["mean", "0.1f", "mm", "Target Position (Elevation)"], "target_position_ax_mm": ["mean", "0.1f", "mm", "Target Position (Axial)"], @@ -50,7 +51,7 @@ "duty_cycle_pulse_train_pct": [None, "0.1f", "%", "Pulse Train Duty Cycle"], "duty_cycle_sequence_pct": [None, "0.1f", "%", "Sequence Duty Cycle"], "sequence_duration_s": [None, "0.0f", "s", "Sequence Duration"], - "estimated_tx_temperature_rise_C": [None, "0.2f", "°C", "Estimated TX Temperature Rise"]} + "estimated_tx_temperature_rise_C": [None, "0.2f", "°C", "Est. Transmitter Heating"]} @dataclass class SolutionAnalysis(DictMixin): @@ -237,40 +238,84 @@ def to_json(self, compact:bool) -> str: @dataclass class SolutionAnalysisOptions(DictMixin): - standoff_sound_speed: Annotated[float, OpenLIFUFieldData("Standoff sound speed (m/s)", "Speed of sound in standoff, for calculating initial impedance")] = 1500.0 + standoff_sound_speed: Annotated[float, OpenLIFUFieldData( + name="Standoff sound speed", + description="Speed of sound in standoff, for calculating initial impedance", + units="m/s", precision=0, + )] = 1500.0 """Speed of sound in standoff, for calculating initial impedance""" - standoff_density: Annotated[float, OpenLIFUFieldData("Standoff density (kg/m³)", "Density of standoff medium (kg/m³)")] = 1000.0 + standoff_density: Annotated[float, OpenLIFUFieldData( + name="Standoff density", + description="Density of standoff medium", + units="kg/m^3", precision=0, + )] = 1000.0 """Density of standoff medium (kg/m³)""" - ref_sound_speed: Annotated[float, OpenLIFUFieldData("Reference sound speed (m/s)", "Reference speed of sound in the medium (m/s)")] = 1500.0 + ref_sound_speed: Annotated[float, OpenLIFUFieldData( + name="Reference sound speed", + description="Reference speed of sound in the medium", + units="m/s", precision=0, + )] = 1500.0 """Reference speed of sound in the medium (m/s)""" - ref_density: Annotated[float, OpenLIFUFieldData("Reference density (kg/m³)", "Reference density (kg/m³)")] = 1000.0 + ref_density: Annotated[float, OpenLIFUFieldData( + name="Reference density", + description="Reference density", + units="kg/m^3", precision=0, + )] = 1000.0 """Reference density (kg/m³)""" - mainlobe_aspect_ratio: Annotated[Tuple[float, float, float], OpenLIFUFieldData("Mainlobe aspect ratio (lat,ele,ax)", "Aspect ratio of the mainlobe mask")] = (1., 1., 5.) - """Aspect ratio of the mainlobe ellipsoid mask, in the form (lat,ele,ax). (1,1,5) means an ellipsoid 5x as long as it is wide.""" - - mainlobe_radius: Annotated[float, OpenLIFUFieldData("Mainlobe mask radius", "Size of the mainlobe mask, in the units provided for Distance units (`distance_units`)")] = 2.5e-3 - """Size of the mainlobe mask, in the units provided for Distance units (`distance_units`). The mainlobe mask is an ellipsoid with this radius, scaled by the `mainlobe_aspect_ratio`.""" - - beamwidth_radius: Annotated[float, OpenLIFUFieldData("Beamwidth search radius", "Size of the beamwidth search, in the units provided for Distance units (`distance_units`)")] = 5e-3 + mainlobe_aspect_ratio: Annotated[Tuple[float, float, float], OpenLIFUFieldData( + name="Mainlobe aspect ratio (lat,ele,ax)", + description="Aspect ratio of the mainlobe mask", + precision=1, + )] = (1., 1., 7.) + """Aspect ratio of the mainlobe ellipsoid mask, in the form (lat,ele,ax). (1,1,7) means an ellipsoid 7x as long as it is wide.""" + + mainlobe_radius: Annotated[float | None, OpenLIFUFieldData( + name="Mainlobe mask radius", + description="Size of the mainlobe mask, in the units provided for Distance units (`distance_units`)", + units_field="distance_units", display_units="mm", precision=2, + )] = None + """Size of the mainlobe mask, in the units provided for Distance units (`distance_units`). The mainlobe mask is an ellipsoid with this radius, scaled by the `mainlobe_aspect_ratio`. If not provided, will be calculated from estimated beamwidth""" + + beamwidth_radius: Annotated[float | None, OpenLIFUFieldData( + name="Beamwidth search radius", + description="Size of the beamwidth search, in the units provided for Distance units (`distance_units`)", + units_field="distance_units", display_units="mm", precision=2, + )] = None """Size of the beamwidth search, in the units provided for Distance units (`distance_units`). The beamwidth is found along the lateral and elevation lines perpendicular to the focus axis.""" - sidelobe_radius: Annotated[float, OpenLIFUFieldData("Sidelobe radius", "Size of the sidelobe mask, in the units provided for Distance units (`distance_units`)")] = 3e-3 - """Size of the sidelobe mask, in the units provided for Distance units (`distance_units`). Pressure outside of this ellipsoid (scaled by `mainlobe_aspect_ratio`) is considered outside of the focal region.""" - - sidelobe_zmin: Annotated[float, OpenLIFUFieldData("Sidelobe minimum z", "Minimum z coordinate of the sidelobe mask, in the units provided for Distance units (`distance_units`)")] = 1e-3 + sidelobe_radius: Annotated[float | None, OpenLIFUFieldData( + name="Sidelobe radius", + description="Size of the sidelobe mask, in the units provided for Distance units (`distance_units`)", + units_field="distance_units", display_units="mm", precision=2, + )] = None + """Size of the sidelobe mask, in the units provided for Distance units (`distance_units`). Pressure outside of this ellipsoid (scaled by `mainlobe_aspect_ratio`) is considered outside of the focal region. If not provided, will be estimated from the beamwidth * 1.5""" + + sidelobe_zmin: Annotated[float | None, OpenLIFUFieldData( + name="Sidelobe minimum z", + description="Minimum z coordinate of the sidelobe mask, in the units provided for Distance units (`distance_units`)", + units_field="distance_units", display_units="mm", precision=2, + )] = None """Minimum z coordinate of the sidelobe mask, in the units provided for Distance units (`distance_units`). This value is used to ignore emitted pressure artifacts.""" - distance_units: Annotated[str, OpenLIFUFieldData("Distance units", "The units used for distance measurements")] = "m" + distance_units: Annotated[str, OpenLIFUFieldData( + name="Distance units", + description="The units used for distance measurements", + unit_options=("mm", "cm", "m"), + )] = "m" """The units used for distance measurements""" param_constraints: Annotated[Dict[str, ParameterConstraint], OpenLIFUFieldData("Parameter constraints", None)] = field(default_factory=dict) """TODO: Add description""" def __post_init__(self): + if not isinstance(self.distance_units, str): + raise TypeError("Distance units must be a string") + if getunittype(self.distance_units) != 'distance': + raise ValueError(f"Distance units must be a length unit, got {self.distance_units}") if self.standoff_sound_speed <= 0: raise ValueError("Standoff sound speed must be greater than 0") if self.standoff_density <= 0: @@ -284,18 +329,16 @@ def __post_init__(self): self.mainlobe_aspect_ratio = tuple(self.mainlobe_aspect_ratio) # Ensure it's a tuple if not all(isinstance(x, int | float) for x in self.mainlobe_aspect_ratio): raise TypeError("Mainlobe aspect ratio must contain only numbers") - if not isinstance(self.mainlobe_radius, int | float) or self.mainlobe_radius <= 0: + if self.mainlobe_radius is not None and (not isinstance(self.mainlobe_radius, int | float) or self.mainlobe_radius <= 0): raise ValueError("Mainlobe radius must be a positive number") - if not isinstance(self.beamwidth_radius, int | float) or self.beamwidth_radius <= 0: + if self.beamwidth_radius is not None and (not isinstance(self.beamwidth_radius, int | float) or self.beamwidth_radius <= 0): raise ValueError("Beamwidth radius must be a positive number") - if not isinstance(self.sidelobe_radius, int | float) or self.sidelobe_radius <= 0: + if self.sidelobe_radius is not None and (not isinstance(self.sidelobe_radius, int | float) or self.sidelobe_radius <= 0): raise ValueError("Sidelobe radius must be a positive number") + if self.sidelobe_zmin is None: + self.sidelobe_zmin = getunitconversion("mm", self.distance_units) * DEFAULT_SIDELOBE_ZMIN_MM if not isinstance(self.sidelobe_zmin, int | float) or self.sidelobe_zmin < 0: raise ValueError("Sidelobe minimum z must be a non-negative number") - if not isinstance(self.distance_units, str): - raise TypeError("Distance units must be a string") - if getunittype(self.distance_units) != 'distance': - raise ValueError(f"Distance units must be a length unit, got {self.distance_units}") @classmethod def from_dict(cls: Type[SolutionAnalysisOptions], parameter_dict: Dict[str, Any]) -> SolutionAnalysisOptions: @@ -313,6 +356,15 @@ def from_dict(cls: Type[SolutionAnalysisOptions], parameter_dict: Dict[str, Any] return cls(**parameter_dict) + def get_summary(self) -> str: + """Return a one-liner summary of the analysis options. + + Returns an empty string: the solution-analysis options are too + numerous to render meaningfully on a collapsible header, so callers + should fall back to showing only the section title. + """ + return "" + def find_centroid(da: xa.DataArray, cutoff:float, units:None) -> np.ndarray: """Find the centroid of a thresholded region of a DataArray""" if units is not None and getunittype(units) != 'distance': @@ -351,11 +403,11 @@ def get_focus_matrix(focus, origin=[0,0,0]) -> np.ndarray: M[3,3] = 1 return M -def get_gridded_transformed_coords(da: xa.DataArray, matrix: np.ndarray, as_dataset=True): - """Transform the coords of a DataArray using a transform matrix. +def get_gridded_transformed_coords(da: xa.DataArray | xa.Dataset, matrix: np.ndarray, as_dataset=True): + """Transform the coords of a DataArray or Dataset using a transform matrix. Args: - da: DataArray whose coordinates will be used + da: DataArray or Dataset whose coordinates will be used matrix: a 4x4 coordinate transformation matrix, transforming from the desired coordinate system to the coordinate system of `da` as_dataset: Whether to return the transformed coords as a numpy array or an xarray Dataset @@ -372,13 +424,13 @@ def get_gridded_transformed_coords(da: xa.DataArray, matrix: np.ndarray, as_data coords = xa.Dataset({f'd_{dim}': (da.dims, coords[...,i]) for i, dim in enumerate(da.dims)}, coords=da.coords) return coords -def get_offset_grid(da: xa.DataArray, focus, origin=DEFAULT_ORIGIN, as_dataset=True): +def get_offset_grid(da: xa.DataArray | xa.Dataset, focus, origin=DEFAULT_ORIGIN, as_dataset=True): """Transform the coords of a DataArray that is in transducer coordinates to focus coordinates See `get_focus_matrix` for the meaning of "focus coordinates" Args: - da: DataArray whose coordinates will be used (presumably the transducer coordinates) + da: DataArray or Dataset whose coordinates will be used (presumably the transducer coordinates) focus: A 3D point describing the focus location in the coordinates of `da` origin: A 3D point describing the "effective origin" in the coordinates of `da` (see `Transducer.get_effective_origin` for the meaning of this). @@ -391,12 +443,12 @@ def get_offset_grid(da: xa.DataArray, focus, origin=DEFAULT_ORIGIN, as_dataset=T coords = get_gridded_transformed_coords(da, M, as_dataset=as_dataset) return coords -def calc_dist_from_focus(da: xa.DataArray, focus, origin=DEFAULT_ORIGIN, aspect_ratio=[1,1,1], as_dataarray=True): +def calc_dist_from_focus(da: xa.DataArray | xa.Dataset, focus, origin=DEFAULT_ORIGIN, aspect_ratio=[1,1,1], as_dataarray=True): """Compute a distance map from a focus point in transducer space, using a possibly distorted metric that respects the symmetry of the focus shape (e.g. it could be cigar-shaped). Args: - da: DataArray that will supply the coordnate grid (presumably transducer coordinates) + da: DataArray or Dataset that will supply the coordnate grid (presumably transducer coordinates) focus: A 3D point describing the focus location in the coordinates of `da` origin: A 3D point describing the "effective origin" in the coordinates of `da` (see `Transducer.get_effective_origin` for the meaning of this). @@ -413,7 +465,7 @@ def calc_dist_from_focus(da: xa.DataArray, focus, origin=DEFAULT_ORIGIN, aspect_ return dist def get_mask( - da: xa.DataArray, + da: xa.DataArray | xa.Dataset, focus, distance:float, origin=DEFAULT_ORIGIN, @@ -425,7 +477,7 @@ def get_mask( The focus region is an ellipsoid centered at the focus point. Args: - da: DataArray that will supply the coordnate grid (presumably transducer coordinates) + da: DataArray or Dataset that will supply the coordnate grid (presumably transducer coordinates) focus: A 3D point describing the focus location in the coordinates of `da` distance: How far from the `focus` to include points in the mask. See `calc_dist_from_focus` for the distorted metric under which a ball of points becomes an ellispoid in euclidean space. @@ -504,6 +556,7 @@ def get_beam_bounds( origin=DEFAULT_ORIGIN, min_offset:float | None=None, max_offset:float | None=None, + clip_to_bounds:bool=True, ) -> Tuple[float, float]: """Determine how far along a focal coordinate system axis a DataArray's value stays above a certain cutoff. @@ -535,11 +588,15 @@ def get_beam_bounds( da_negoff = da_negoff.where(da_negoff < float(cutoff), drop=True) if da_negoff.size > 0: negoff = float(da_negoff.coords[f'offset_d{dim}'][-1]) + elif clip_to_bounds: + negoff = float(interp_da.coords[f'offset_d{dim}'][0]) else: negoff = np.nan da_posoff = da_posoff.where(da_posoff < float(cutoff), drop=True) if da_posoff.size > 0: posoff = float(da_posoff.coords[f'offset_d{dim}'][0]) + elif clip_to_bounds: + posoff = float(interp_da.coords[f'offset_d{dim}'][-1]) else: posoff = np.nan return negoff, posoff @@ -621,16 +678,16 @@ def model_tx_temperature_rise(voltage: float, T0 = T0_degC if T0 < 20 or T0 > 40: - logger.warning("Initial temperature T0 must be between 20 and 40 degrees Celsius for the model to be valid.") + logger.debug("Initial temperature T0 must be between 20 and 40 degrees Celsius for the electronics thermal model to be valid.") if P < 50 or P > 500: - logger.warning("Squared Voltage must be between 50 and 500 V^2 for the model to be valid.") + logger.debug("Squared Voltage must be between 50 and 500 V^2 for the electronics thermal model to be valid.") if t < 1 or t > 600: - logger.warning("Time t must be between 1 and 600 seconds for the model to be valid.") + logger.debug("Time t must be between 1 and 600 seconds for the electronics thermal model to be valid.") if frequency_kHz < 380 or frequency_kHz > 420: - logger.warning("Frequency must be between 380 and 420 kHz for the model to be valid.") + logger.debug("Frequency must be between 380 and 420 kHz for the electronics thermal model to be valid.") # Predict power law parameters using polynomial regression (degree 2) n = (2.131832 + -0.003475*P + -0.044916*T0 + diff --git a/src/openlifu/seg/material.py b/src/openlifu/seg/material.py index 2e79d913..19b77750 100644 --- a/src/openlifu/seg/material.py +++ b/src/openlifu/seg/material.py @@ -99,28 +99,28 @@ def from_dict(d: dict[str, Any]): TISSUE = Material(name="tissue", sound_speed=1540.0, density=1000.0, - attenuation=0.0, + attenuation=0.3, specific_heat=3600.0, thermal_conductivity=0.5) SKULL = Material(name="skull", sound_speed=4080.0, density=1900.0, - attenuation=0.0, + attenuation=0.6, specific_heat=1100.0, thermal_conductivity=0.3) AIR = Material(name="air", sound_speed=344.0, density=1.25, - attenuation=0.0, + attenuation=1.0, specific_heat=1012.0, thermal_conductivity=0.025) STANDOFF = Material(name="standoff", sound_speed=1420.0, density=1000.0, - attenuation=1.0, + attenuation=0.2, specific_heat=4182.0, thermal_conductivity=0.598) diff --git a/src/openlifu/seg/seg_method.py b/src/openlifu/seg/seg_method.py index 82c2270d..e143c89f 100644 --- a/src/openlifu/seg/seg_method.py +++ b/src/openlifu/seg/seg_method.py @@ -122,3 +122,19 @@ def to_table(self) -> pd.DataFrame: :returns: Pandas DataFrame of the segmentation method parameters """ pass + + def get_summary(self) -> str: + """Return a one-liner summary of the segmentation method. + + Default implementation returns the human-friendly form of the class name + (e.g. ``"Uniform Tissue"``); subclasses may override to provide more + detail. + """ + # Insert spaces before capital letters: "UniformTissue" -> "Uniform Tissue" + name = type(self).__name__ + result = [] + for i, ch in enumerate(name): + if i > 0 and ch.isupper() and not name[i - 1].isupper(): + result.append(" ") + result.append(ch) + return "".join(result) diff --git a/src/openlifu/seg/virtual_fit.py b/src/openlifu/seg/virtual_fit.py index dd616885..1c318ca6 100644 --- a/src/openlifu/seg/virtual_fit.py +++ b/src/openlifu/seg/virtual_fit.py @@ -52,45 +52,92 @@ class VirtualFitOptions(DictMixin): yaw: 90 degrees minus the polar spherical coordinate. """ - units: Annotated[str, OpenLIFUFieldData("Length units", "The units of length used in the length attributes of this class")] = "mm" + units: Annotated[str, OpenLIFUFieldData( + name="Length units", + description="The units of length used in the length attributes of this class", + unit_options=("mm", "cm", "m"), + )] = "mm" """The units of length used in the length attributes of this class""" - transducer_steering_center_distance: Annotated[float, OpenLIFUFieldData("Steering center distance", "Distance from the transducer origin axially to the center of the steering zone in the units `units`")] = 50. + transducer_steering_center_distance: Annotated[float, OpenLIFUFieldData( + name="Steering center distance", + description="Distance from the transducer origin axially to the center of the steering zone", + units_field="units", display_units="mm", precision=2, + )] = 50. """Distance from the transducer origin axially to the center of the steering zone in the units `units`""" steering_limits: Annotated[Tuple[Tuple[float, float], Tuple[float, float], Tuple[float, float]], - OpenLIFUFieldData("Steering limits", "Steering bounds along each axis from the transducer origin, in the units `units`")] = ((-50, 50), (-50, 50), (-50, 50)) + OpenLIFUFieldData( + name="Steering limits", + description="Steering bounds along each axis from the transducer origin", + units_field="units", display_units="mm", precision=1, + )] = ((-50, 50), (-50, 50), (-50, 50)) """Distance from the transducer origin axially to the center of the steering zone in the units `units`""" - pitch_range: Annotated[Tuple[float, float], OpenLIFUFieldData("Pitch range (deg)", "Range of pitches to include in the transducer fitting search grid, in degrees")] = (-10, 150) + pitch_range: Annotated[Tuple[float, float], OpenLIFUFieldData( + name="Pitch range", + description="Range of pitches to include in the transducer fitting search grid", + units="deg", precision=0, + )] = (-10, 150) """Range of pitches to include in the transducer fitting search grid, in degrees""" - pitch_step: Annotated[float, OpenLIFUFieldData("Pitch step size (deg)", "Pitch step size when forming the transducer fitting search grid, in degrees")] = 5 + pitch_step: Annotated[float, OpenLIFUFieldData( + name="Pitch step size", + description="Pitch step size when forming the transducer fitting search grid", + units="deg", precision=1, + )] = 5 """Pitch step size when forming the transducer fitting search grid, in degrees""" - yaw_range: Annotated[Tuple[float, float], OpenLIFUFieldData("Yaw range (deg)", "Range of yaws to include in the transducer fitting search grid, in degrees")] = (-65, 65) + yaw_range: Annotated[Tuple[float, float], OpenLIFUFieldData( + name="Yaw range", + description="Range of yaws to include in the transducer fitting search grid", + units="deg", precision=0, + )] = (-65, 65) """Range of yaws to include in the transducer fitting search grid, in degrees""" - yaw_step: Annotated[float, OpenLIFUFieldData("Yaw step size (deg)", "Yaw step size when forming the transducer fitting search grid, in degrees")] = 5 + yaw_step: Annotated[float, OpenLIFUFieldData( + name="Yaw step size", + description="Yaw step size when forming the transducer fitting search grid", + units="deg", precision=1, + )] = 5 """Yaw step size when forming the transducer fitting search grid, in degrees""" - planefit_dyaw_extent: Annotated[float, OpenLIFUFieldData("Plane fit yaw extent", "Left and right extents of the point grid to be used for plane fitting along the local yaw axes, in units of `units`")] = 15 + planefit_dyaw_extent: Annotated[float, OpenLIFUFieldData( + name="Plane fit yaw extent", + description="Left and right extents of the point grid to be used for plane fitting along the local yaw axes", + units_field="units", display_units="mm", precision=2, + )] = 15 """Left and right extents of the point grid to be used for plane fitting along the local yaw axes, in units of `units`. The plane fitting point grid will be twice this size, since this is left and right extents. (Note that this has units of length, not angle!)""" - planefit_dyaw_step: Annotated[float, OpenLIFUFieldData("Plane fit yaw step", "Local yaw axis step size to use when constructing plane fitting grids. In spatial units of `units`")] = 3 + planefit_dyaw_step: Annotated[float, OpenLIFUFieldData( + name="Plane fit yaw step", + description="Local yaw axis step size to use when constructing plane fitting grids", + units_field="units", display_units="mm", precision=2, + )] = 3 """Local yaw axis step size to use when constructing plane fitting grids. In spatial units of `units`.""" - planefit_dpitch_extent: Annotated[float, OpenLIFUFieldData("Plane fit pitch extent", "Left and right extents of the point grid to be used for plane fitting along the local pitch axes, in spatial units of `units`")] = 15 + planefit_dpitch_extent: Annotated[float, OpenLIFUFieldData( + name="Plane fit pitch extent", + description="Left and right extents of the point grid to be used for plane fitting along the local pitch axes", + units_field="units", display_units="mm", precision=2, + )] = 15 """Left and right extents of the point grid to be used for plane fitting along the local pitch axes, in spatial units of `units`. The plane fitting point grid will be twice this size, since this is left and right extents.""" - planefit_dpitch_step: Annotated[float, OpenLIFUFieldData("Plane fit pitch step", "Local pitch axis step size to use when constructing plane fitting grids. In spatial units of `units`")] = 3 + planefit_dpitch_step: Annotated[float, OpenLIFUFieldData( + name="Plane fit pitch step", + description="Local pitch axis step size to use when constructing plane fitting grids", + units_field="units", display_units="mm", precision=2, + )] = 3 """Local pitch axis step size to use when constructing plane fitting grids. In spatial units of `units`.""" - top_n_candidates: Annotated[int, OpenLIFUFieldData("No. of candidates returned", "Sets the limit for the number of transducer transform candidates returned by the algorithm.")] = 4 + top_n_candidates: Annotated[int, OpenLIFUFieldData( + name="No. of candidates returned", + description="Sets the limit for the number of transducer transform candidates returned by the algorithm.", + )] = 4 """Sets the limit for the number of transducer transform candidates returned by the algorithm.""" def __post_init__(self): @@ -169,6 +216,15 @@ def from_dict(parameter_dict: Dict[str,Any]) -> VirtualFitOptions: # Override Di parameter_dict["steering_limits"] = tuple(map(tuple,parameter_dict["steering_limits"])) return VirtualFitOptions(**parameter_dict) + def get_summary(self) -> str: + """Return a one-liner summary of the virtual-fit options. + + Returns an empty string: the virtual-fit options are too numerous to + render meaningfully on a collapsible header, so callers should fall + back to showing only the section title. + """ + return "" + def compute_skin_mesh_from_volume( volume_array : np.ndarray, volume_affine_RAS : np.ndarray, diff --git a/src/openlifu/sim/__init__.py b/src/openlifu/sim/__init__.py index 2f7c1b79..f8bb8d49 100644 --- a/src/openlifu/sim/__init__.py +++ b/src/openlifu/sim/__init__.py @@ -2,8 +2,16 @@ from openlifu.sim.kwave_if import run_simulation from openlifu.sim.sim_setup import SimSetup +from openlifu.sim.thermal import ( + compute_heat_source, + generate_pulse_events, + run_thermal_simulation, +) __all__ = [ "SimSetup", - "run_simulation" + "compute_heat_source", + "generate_pulse_events", + "run_simulation", + "run_thermal_simulation", ] diff --git a/src/openlifu/sim/sim_setup.py b/src/openlifu/sim/sim_setup.py index a8cffbd5..d719b919 100644 --- a/src/openlifu/sim/sim_setup.py +++ b/src/openlifu/sim/sim_setup.py @@ -21,31 +21,67 @@ @dataclass class SimSetup(DictMixin): - spacing: Annotated[float, OpenLIFUFieldData("Spacing", "Simulation grid spacing")] = 1.0 + spacing: Annotated[float, OpenLIFUFieldData( + name="Voxel spacing", + description="Simulation grid spacing", + units_field="units", display_units="mm", precision=2, + )] = 1.0 """Simulation grid spacing""" - units: Annotated[str, OpenLIFUFieldData("Spatial units", "Units used for spatial measurements")] = "mm" + units: Annotated[str, OpenLIFUFieldData( + name="Spatial units", + description="Units used for spatial measurements", + unit_options=("mm", "cm", "m"), + )] = "mm" """Units used for spatial measurements""" - x_extent: Annotated[Tuple[float, float], OpenLIFUFieldData("X-extent", "Simulation grid extent along the first dimension")] = (-30., 30.) + x_extent: Annotated[Tuple[float, float], OpenLIFUFieldData( + name="X-extent", + description="Simulation grid extent along the first dimension", + units_field="units", display_units="mm", precision=1, + )] = (-30., 30.) """Simulation grid extent along the first dimension""" - y_extent: Annotated[Tuple[float, float], OpenLIFUFieldData("Y-extent", "Simulation grid extend along the second dimension")] = (-30., 30.) + y_extent: Annotated[Tuple[float, float], OpenLIFUFieldData( + name="Y-extent", + description="Simulation grid extent along the second dimension", + units_field="units", display_units="mm", precision=1, + )] = (-30., 30.) """Simulation grid extend along the second dimension""" - z_extent: Annotated[Tuple[float, float], OpenLIFUFieldData("Z-extent", "Simulation grid extend along the third dimension")] = (-4., 60.) + z_extent: Annotated[Tuple[float, float], OpenLIFUFieldData( + name="Z-extent", + description="Simulation grid extent along the third dimension", + units_field="units", display_units="mm", precision=1, + )] = (-4., 60.) """Simulation grid extend along the third dimension""" - dt: Annotated[float, OpenLIFUFieldData("Time step", "Simulation time step")] = 0. + dt: Annotated[float, OpenLIFUFieldData( + name="Time step", + description="Simulation time step", + units="s", precision=6, + )] = 0. """Simulation time step""" - t_end: Annotated[float, OpenLIFUFieldData("End time", """Simulation end time""")] = 0. + t_end: Annotated[float, OpenLIFUFieldData( + name="End time", + description="Simulation end time", + units="s", precision=6, + )] = 0. """Simulation end time""" - c0: Annotated[float, OpenLIFUFieldData("Speed of Sound (m/s)", "Reference speed of sound for converting distance to time")] = 1500.0 + c0: Annotated[float, OpenLIFUFieldData( + name="Default speed of sound", + description="Reference speed of sound for converting distance to time", + units="m/s", precision=0, + )] = 1500.0 """Reference speed of sound for converting distance to time""" - cfl: Annotated[float, OpenLIFUFieldData("CFL number", "Courant-Friedrichs-Lewy number")] = 0.3 + cfl: Annotated[float, OpenLIFUFieldData( + name="CFL number", + description="Courant-Friedrichs-Lewy number", + precision=2, + )] = 0.3 """Courant-Friedrichs-Lewy number""" options: Annotated[dict[str, str], OpenLIFUFieldData("Simulation options", "Additional simulation options")] = field(default_factory=dict) @@ -205,6 +241,26 @@ def to_table(self) -> pd.DataFrame: ] return pd.DataFrame.from_records(records) + def get_summary(self) -> str: + """Return a one-liner summary of the simulation setup parameters. + + Format: ``"{spacing}mm spacing, [{x0},{x1}]x[{y0},{y1}]x[{z0},{z1}]"``. + Spacing and extents are converted to millimeters regardless of the + configured storage units. + """ + try: + scale = getunitconversion(self.units, "mm") + except ValueError: + scale = 1.0 + spacing_mm = self.spacing * scale + x0, x1 = self.x_extent[0] * scale, self.x_extent[1] * scale + y0, y1 = self.y_extent[0] * scale, self.y_extent[1] * scale + z0, z1 = self.z_extent[0] * scale, self.z_extent[1] * scale + return ( + f"{spacing_mm:g}mm spacing, " + f"[{x0:g},{x1:g}]x[{y0:g},{y1:g}]x[{z0:g},{z1:g}]" + ) + @staticmethod def from_dict(d: dict, on_keyword_mismatch: Literal['warn', 'raise', 'ignore'] = 'warn') -> SimSetup: """Create a SimSetup instance from a dictionary.""" diff --git a/src/openlifu/sim/thermal.py b/src/openlifu/sim/thermal.py new file mode 100644 index 00000000..3254fa28 --- /dev/null +++ b/src/openlifu/sim/thermal.py @@ -0,0 +1,636 @@ +from __future__ import annotations + +import logging +import math +import time +from typing import List, Literal, Sequence, Tuple + +import numpy as np +import xarray as xa + +from openlifu.util.units import getunitconversion + +logger = logging.getLogger(__name__) + + +def db2neper(alpha_db: np.ndarray | float, y: float = 1.0) -> np.ndarray | float: + """Convert attenuation from dB / (MHz^y cm) to Nepers / ((rad/s)^y m). + + Matches the convention used by k-Wave's ``db2neper``. Multiplying the returned + value by ``(2*pi*f)**y`` (with ``f`` in Hz) gives an attenuation coefficient + in Np/m suitable for computing volumetric heat deposition from pressure. + """ + return 100.0 * alpha_db * (1e-6 / (2.0 * math.pi)) ** y / (20.0 * np.log10(np.exp(1))) + + +def compute_heat_source( + params: xa.Dataset, + pressure: np.ndarray, + freq: float, + alpha_power: float = 0.9, +) -> np.ndarray: + """Compute the volumetric rate of heat deposition from an acoustic pressure field. + + Uses the standard absorbed-power relation ``Q = alpha * p^2 / (rho * c)`` + where ``alpha`` is the attenuation coefficient in Np/m, ``p`` is the peak + pressure amplitude in Pa, ``rho`` is density in kg/m^3, and ``c`` is speed + of sound in m/s. The result has units of W/m^3. + + Args: + params: Tissue parameter maps as an xarray Dataset. Must include + ``density`` (kg/m^3), ``sound_speed`` (m/s), and ``attenuation`` + (dB/cm/MHz) variables. + pressure: Peak pressure amplitude (Pa). Shape must be broadcastable to + the parameter maps (typically ``(nx, ny, nz)`` or + ``(num_foci, nx, ny, nz)``). + freq: Frequency in Hz used when converting attenuation to Np/m. + alpha_power: Power law exponent for attenuation. Default 0.9. + + Returns: + Q with the same shape as ``pressure`` and units of W/m^3. + """ + alpha_db = np.asarray(params['attenuation'].data) + density = np.asarray(params['density'].data) + sound_speed = np.asarray(params['sound_speed'].data) + alpha_np = db2neper(alpha_db, y=alpha_power) * (2.0 * np.pi * freq) ** alpha_power + return alpha_np * np.asarray(pressure) ** 2 / (density * sound_speed) + + +def generate_pulse_events( + pulse_count: int, + pulse_interval: float, + pulse_train_count: int, + pulse_train_interval: float, + pulse_duration: float, + num_foci: int, + order: Sequence[int] | None = None, +) -> List[Tuple[float, int]]: + """Generate a list of ``(time, focus_index)`` events for a pulse sequence. + + Positive ``focus_index`` (0-indexed) marks the start of a sonication at that + focus. ``focus_index == -1`` marks the end of a sonication (Q returns to 0). + + Args: + pulse_count: Number of pulses per pulse train. + pulse_interval: Time between the start of successive pulses in a train (s). + pulse_train_count: Number of pulse trains in the sequence. + pulse_train_interval: Time between the start of successive pulse trains (s). + If 0, trains are treated as back-to-back with no gap. + pulse_duration: Duration of a single pulse (s). + num_foci: Number of foci to cycle through. + order: Optional list of 1-indexed focus indices specifying the order in + which foci are sonicated. If None, the foci are cycled through in + order 1, 2, ..., num_foci, matching the MATLAB LOFUgen1 convention. + + Returns: + A list of ``(time_s, focus_index)`` tuples, sorted by time. + """ + if num_foci <= 0: + return [] + if order is None: + order_idx = list(range(num_foci)) + else: + order_idx = [int(i) - 1 for i in order] + if any(i < 0 or i >= num_foci for i in order_idx): + raise ValueError(f"order values must be 1-indexed integers in [1, {num_foci}].") + if len(order_idx) == 0: + return [] + + if pulse_train_interval == 0: + pt_interval = pulse_count * pulse_interval + else: + pt_interval = pulse_train_interval + + events: List[Tuple[float, int]] = [] + for pt in range(pulse_train_count): + pt_start = pt * pt_interval + for p in range(pulse_count): + p_start = pt_start + p * pulse_interval + p_end = p_start + pulse_duration + p_index = order_idx[p % len(order_idx)] + events.append((float(p_start), int(p_index))) + events.append((float(p_end), -1)) + events.sort(key=lambda e: e[0]) + return events + + +def _get_dx(coords: xa.Coordinates) -> Tuple[float, float, float]: + """Return grid spacing along (x, y, z) in meters.""" + dxs = [] + for dim in ('x', 'y', 'z'): + unit = coords[dim].attrs.get('units', 'm') + scl = getunitconversion(unit, 'm') + data = coords[dim].data + if len(data) < 2: + raise ValueError(f"Coordinate '{dim}' must have at least 2 samples.") + dxs.append(float(np.abs(np.diff(data)[0])) * scl) + return tuple(dxs) # type: ignore[return-value] + + +def _diffuse_step( + T: np.ndarray, + dt: float, + dx: float, + dy: float, + dz: float, + k_field: np.ndarray, + inv_rho_cp: np.ndarray, + Q: np.ndarray | float, +) -> np.ndarray: + """One forward-Euler step of the heterogeneous heat equation. + + Solves ``rho*C_p * dT/dt = div(k * grad T) + Q`` with zero-flux (Neumann) + boundary conditions using cell-centered finite differences. The + conductivity at cell faces is the arithmetic mean of neighboring values. + """ + kx = 0.5 * (k_field[1:, :, :] + k_field[:-1, :, :]) + ky = 0.5 * (k_field[:, 1:, :] + k_field[:, :-1, :]) + kz = 0.5 * (k_field[:, :, 1:] + k_field[:, :, :-1]) + + fx = kx * (T[1:, :, :] - T[:-1, :, :]) / dx + fy = ky * (T[:, 1:, :] - T[:, :-1, :]) / dy + fz = kz * (T[:, :, 1:] - T[:, :, :-1]) / dz + + div = np.zeros_like(T) + div[:-1, :, :] += fx / dx + div[1:, :, :] -= fx / dx + div[:, :-1, :] += fy / dy + div[:, 1:, :] -= fy / dy + div[:, :, :-1] += fz / dz + div[:, :, 1:] -= fz / dz + + return T + dt * inv_rho_cp * (div + Q) + + +def _cfl_dt(dx: float, dy: float, dz: float, k_field: np.ndarray, rho_cp: np.ndarray) -> float: + """CFL-stable time step for explicit heat diffusion with heterogeneous coefficients.""" + with np.errstate(divide='ignore', invalid='ignore'): + D = k_field / rho_cp + D_max = float(np.nanmax(D)) + if D_max <= 0 or not np.isfinite(D_max): + return math.inf + return 0.5 / (D_max * (1.0 / dx ** 2 + 1.0 / dy ** 2 + 1.0 / dz ** 2)) + + +def _advance( + T: np.ndarray, + duration: float, + dstep: float, + dx: float, dy: float, dz: float, + k_field: np.ndarray, + inv_rho_cp: np.ndarray, + Q: np.ndarray | float, +) -> np.ndarray: + """Advance ``T`` by ``duration`` using at most ``dstep``-sized forward-Euler steps.""" + if duration <= 0: + return T + n_full = math.floor(duration / dstep + 1e-12) + frac = duration - n_full * dstep + for _ in range(n_full): + T = _diffuse_step(T, dstep, dx, dy, dz, k_field, inv_rho_cp, Q) + if frac > 1e-15: + T = _diffuse_step(T, frac, dx, dy, dz, k_field, inv_rho_cp, Q) + return T + + +def _run_direct( + heat_sources: np.ndarray, + events: List[Tuple[float, int]], + N: int, + dt: float, + dstep: float, + dx: float, dy: float, dz: float, + k_field: np.ndarray, + inv_rho_cp: np.ndarray, + T_init: np.ndarray | None, +) -> np.ndarray: + """Run the thermal simulation via direct time-stepping through events. + + Mirrors the LOFUgen1 MATLAB ``run_thermal_sim`` control flow: at each + output step the simulation is advanced by ``dt`` in ``dstep``-sized + increments, with fractional sub-steps taken to align exactly with event + times where Q changes. + """ + shape = heat_sources.shape[1:] + if T_init is None: + T = np.zeros(shape, dtype=np.float32) + else: + T = T_init.astype(np.float32).copy() + + T_out = np.zeros((N, *shape), dtype=np.float32) + # Sentinel event to terminate the while loop below without special-casing + sorted_events = [*sorted(events, key=lambda e: e[0]), (math.inf, -1)] + ii = 0 + Q: np.ndarray | float = 0.0 + + for i in range(N): + T_out[i] = T + # Compute step boundaries directly from i to avoid the floating-point + # drift that would otherwise accumulate in ``t += dt`` and cause + # events at exact grid points (e.g. t = 0.75 with dt = 0.05) to fall + # into the wrong output step. + t = i * dt + t_target = (i + 1) * dt + if sorted_events[ii][0] > t_target: + T = _advance(T, dt, dstep, dx, dy, dz, k_field, inv_rho_cp, Q) + else: + t1 = t + while sorted_events[ii][0] <= t_target: + e_time = sorted_events[ii][0] + if e_time > t1: + T = _advance(T, e_time - t1, dstep, dx, dy, dz, k_field, inv_rho_cp, Q) + t1 = e_time + e_idx = sorted_events[ii][1] + Q = 0.0 if e_idx < 0 else heat_sources[e_idx] + ii += 1 + if t_target > t1: + T = _advance(T, t_target - t1, dstep, dx, dy, dz, k_field, inv_rho_cp, Q) + return T_out + + +def _run_superpose( + heat_sources: np.ndarray, + events: List[Tuple[float, int]], + N: int, + dt: float, + dstep: float, + dx: float, dy: float, dz: float, + k_field: np.ndarray, + inv_rho_cp: np.ndarray, + T_init: np.ndarray | None, + pulse_duration: float, +) -> np.ndarray: + """Run the thermal simulation by superposition of per-focus single-pulse responses. + + For each focus ``i`` we integrate the diffusion equation with ``Q_i`` on + for ``pulse_duration``, then off, sampling ``T_i(x, k*dt)`` for ``k=0..N-1``. + Because the diffusion equation is linear, the total temperature rise for + the pulse sequence is the sum of appropriately time-shifted single-pulse + responses. + + This is typically much faster than the direct mode when the sequence + contains many pulses (each pulse contributes an O(N * grid) addition, + whereas a direct time-stepped simulation must resolve every pulse edge). + Pulse event times are snapped to the nearest ``dt`` sample; a warning is + logged if the snapping error is significant. + """ + nfoc = heat_sources.shape[0] + shape = heat_sources.shape[1:] + + if T_init is None: + T_out = np.zeros((N, *shape), dtype=np.float32) + else: + T_out = np.broadcast_to(T_init.astype(np.float32), (N, *shape)).copy() + + # Bucket "on" events by focus and snap their times to the output grid. + max_snap_err = 0.0 + per_focus_events: List[List[int]] = [[] for _ in range(nfoc)] + for event_time, focus_idx in events: + if focus_idx < 0: + continue + j_e_float = event_time / dt + j_e = round(j_e_float) + snap_err = abs(j_e - j_e_float) * dt + max_snap_err = max(max_snap_err, snap_err) + if j_e >= N: + continue + per_focus_events[focus_idx].append(j_e) + + if max_snap_err > 0.05 * dt: + logger.warning( + "Event times are not aligned to the output time step (dt=%g); " + "max snap error is %g s. Consider decreasing dt for better precision.", + dt, max_snap_err, + ) + + for i in range(nfoc): + if not per_focus_events[i]: + continue + logger.info("Computing single-pulse response for focus %d/%d ...", i + 1, nfoc) + j_events = np.asarray(sorted(per_focus_events[i]), dtype=np.int64) + T_i = np.zeros(shape, dtype=np.float32) + Q_active = np.asarray(heat_sources[i], dtype=np.float32) + for j in range(N): + # Distribute T_i (which equals T_i(x, j*dt)) into T_out for every + # on-event at this focus: contribution to sample j_e + j. + valid = j_events + j + valid = valid[valid < N] + for j_out in valid: + T_out[j_out] += T_i + # Advance T_i by one output step from tau=j*dt to (j+1)*dt. + tau = j * dt + tau_target = tau + dt + if tau_target <= pulse_duration + 1e-15: + T_i = _advance(T_i, dt, dstep, dx, dy, dz, k_field, inv_rho_cp, Q_active) + elif tau >= pulse_duration - 1e-15: + T_i = _advance(T_i, dt, dstep, dx, dy, dz, k_field, inv_rho_cp, 0.0) + else: + on_dur = pulse_duration - tau + T_i = _advance(T_i, on_dur, dstep, dx, dy, dz, k_field, inv_rho_cp, Q_active) + T_i = _advance(T_i, dt - on_dur, dstep, dx, dy, dz, k_field, inv_rho_cp, 0.0) + return T_out + + +def _run_impulse( + heat_sources: np.ndarray, + events: List[Tuple[float, int]], + N: int, + dt: float, + dstep: float, + dx: float, dy: float, dz: float, + k_field: np.ndarray, + inv_rho_cp: np.ndarray, + T_init: np.ndarray | None, + pulse_duration: float, +) -> np.ndarray: + """Time-marching impulse mode: apply instantaneous deposits at each pulse + event and diffuse between them. + + Marches a single temperature field forward in time (all foci contribute + to the same field). At each output sample time ``t_j = j*dt`` the current + ``T`` is stored *before* any events at that time (matching the + sample-before-events convention used by the other two modes). Between + samples, ``T`` is diffused with ``Q = 0``, and a fractional sub-step is + inserted at every pulse-start event so that the instantaneous deposit + + .. math:: + + E_i(x) = Q_i(x) \\cdot \\mathrm{pulse\\_duration} / (\\rho C_p) + + is applied at exactly the right time. Off-grid events (event times not + aligned to the output-``dt`` grid) are captured *exactly* -- no snapping, + no linear interpolation, no lost energy -- because the sub-step size is + driven by the event stream itself. Multiple events within a single output + step, including events that would fall at fractional positions, all + contribute independently. + + Pulse-end events (``focus_idx == -1``) are ignored: in the impulse limit + the pulse energy has already been deposited instantaneously at the pulse- + start time and no separate turn-off event is needed. + + The internal sub-step size is bounded by ``dstep = dt / oversample``, so + when events are sparse the diffusion is still integrated with a step + small enough to be numerically stable. Total work per focus scales as + ``O((N + n_events) * volume)`` which is dramatically cheaper than the + per-focus superposition scheme when the sequence contains many pulses. + """ + shape = heat_sources.shape[1:] + if T_init is None: + T = np.zeros(shape, dtype=np.float32) + else: + T = T_init.astype(np.float32).copy() + + T_out = np.zeros((N, *shape), dtype=np.float32) + + # Per-focus energy-deposit fields, expressed as an instantaneous + # temperature rise (K). deposit[i](x) = Q_i(x) * pulse_duration / (rho * C_p). + deposits = ( + np.asarray(heat_sources, dtype=np.float32) + * float(pulse_duration) + * inv_rho_cp.astype(np.float32)[None, ...] + ) + + # Sort pulse-start events by time; drop events that would fall outside + # the simulation window. The sentinel event at t=inf lets the inner + # while-loop terminate without a special case. + on_events: List[Tuple[float, int]] = sorted( + ( + (float(t), int(i)) + for t, i in events + if i >= 0 and t >= 0 and t < N * dt + ), + key=lambda e: e[0], + ) + on_events.append((math.inf, -1)) + ii = 0 + for j in range(N): + # Store T at t = j*dt BEFORE any events at this sample time. + T_out[j] = T + # Compute step boundaries directly from j (rather than accumulating + # t += dt) to avoid floating-point drift. With 15 iterations of + # t += 0.05, t drifts to 0.7500000000000001, which then wrongly + # captures events at exactly t=0.75 in the j=14 iteration instead + # of the j=15 iteration. + t = j * dt + t_target = (j + 1) * dt + t1 = t + # Process every event in [t, t_target). Events falling exactly at + # t_target are handled in the next iteration, preserving the + # sample-before-events convention. + while on_events[ii][0] < t_target: + e_time, focus_idx = on_events[ii] + if e_time > t1: + # Diffuse (with Q=0) up to the event time, taking sub-steps + # bounded by dstep. + T = _advance(T, e_time - t1, dstep, dx, dy, dz, k_field, inv_rho_cp, 0.0) + t1 = e_time + # Apply the instantaneous energy deposit for this focus. + T = T + deposits[focus_idx] + ii += 1 + # Diffuse the remainder of the output step to reach t = (j+1)*dt. + if t_target > t1: + T = _advance(T, t_target - t1, dstep, dx, dy, dz, k_field, inv_rho_cp, 0.0) + return T_out + + +def run_thermal_simulation( + heat_sources: xa.DataArray, + params: xa.Dataset, + events: Sequence[Tuple[float, int]], + dt: float = 0.1, + t_end: float | None = None, + oversample: int = 1, + T0: float = 0.0, + T_init: np.ndarray | None = None, + mode: Literal['direct', 'superpose', 'impulse'] = 'impulse', + pulse_duration: float | None = None, +) -> xa.Dataset: + """Run a k-Wave-style thermal diffusion simulation. + + Solves the heterogeneous heat equation + ``rho * C_p * dT/dt = div(k * grad T) + Q(t)`` with zero-flux boundary + conditions using an explicit finite-difference solver. Three evaluation + modes are provided: + + * ``mode='impulse'`` (default): treats each pulse as an instantaneous + energy deposit ``E_i = Q_i * pulse_duration`` (J/m^3) applied at the + pulse-start time. The solver marches a single temperature field + forward, inserting a fractional sub-step at every pulse event so that + the deposit is applied at *exactly* the correct time. Events do NOT + need to be aligned to the output-``dt`` grid: sub-``dt`` timing is + captured exactly, and any number of events per output step -- including + events that would otherwise fall between samples -- deposit their + energy correctly. + * ``mode='superpose'``: precomputes the per-focus temperature response + to a single pulse (with the source held on for ``pulse_duration``, + then off) and superposes shifted copies over the sequence. Use this + when ``pulse_duration`` is comparable to or longer than ``dt`` and the + pulse-duration cannot be treated as instantaneous. + * ``mode='direct'``: mirrors the LOFUgen1 MATLAB ``run_thermal_sim`` and + walks the sequence of pulse-on/pulse-off events, taking fractional + sub-steps to align with event times. + + All three modes use the sample-before-events convention: ``T[j*dt]`` + reports the temperature-rise at time ``j*dt`` *before* any pulse events + that occur at that instant have taken effect. + + Args: + heat_sources: The per-focus volumetric heat deposition rate (W/m^3). + Expected to be an xarray DataArray with a leading + ``focal_point_index`` dimension and spatial dims ``x, y, z``. The + spatial coordinates (with ``units`` attributes) are used to derive + grid spacing. + params: Tissue parameter maps. Must contain ``density`` (kg/m^3), + ``specific_heat`` (J/kg/K), and ``thermal_conductivity`` (W/m/K) + variables on the same grid as ``heat_sources``. + events: A list of ``(time_s, focus_index)`` events. Use + :func:`generate_pulse_events` to build this from an + :class:`openlifu.bf.Sequence`. Off-events (``focus_index == -1``) + are ignored in ``'impulse'`` mode. + dt: Output time step (s) at which the temperature field is stored. + Default 0.1. + t_end: Simulation duration (s). If None, uses the time of the last + event plus one ``dt``. + oversample: Sets the maximum internal computation step: + ``dstep = dt / oversample``. Increase to enforce numerical + stability at fine spatial resolution (see the CFL warning) or to + take smaller diffusion steps between output samples. In + ``'impulse'`` mode the effective internal step is + ``min(dstep, event-to-event gap)``. + T0: Background temperature offset (degC). Only affects the reported + ``Tmax``; the stored ``temperature_rise`` is always + ``T - ambient``. + T_init: Optional initial temperature-rise field ``(nx, ny, nz)``. + mode: ``'impulse'`` (default), ``'superpose'``, or ``'direct'``. + pulse_duration: Duration of a single pulse (s). Required for + ``'impulse'`` and ``'superpose'`` modes. + + Returns: + An xarray Dataset with variables: + + * ``temperature_rise`` (``t, x, y, z``): temperature rise (degC). + * ``Tmax`` (scalar): ``T0 + max(temperature_rise)`` (degC). + * ``sim_time`` (scalar): wall-clock simulation time (s). + + The ``T0`` value used is stored as a Dataset attribute. + """ + if not isinstance(heat_sources, xa.DataArray): + raise TypeError("heat_sources must be an xarray DataArray with a leading 'focal_point_index' dim.") + if 'focal_point_index' not in heat_sources.dims: + raise ValueError("heat_sources must have a 'focal_point_index' dimension.") + heat_sources_t = heat_sources.transpose('focal_point_index', 'x', 'y', 'z') + Q_arr = np.asarray(heat_sources_t.data, dtype=np.float32) + + for name in ('density', 'specific_heat', 'thermal_conductivity'): + if name not in params: + raise ValueError(f"params is missing required variable '{name}'.") + + density = np.asarray(params['density'].data, dtype=np.float32) + specific_heat = np.asarray(params['specific_heat'].data, dtype=np.float32) + thermal_conductivity = np.asarray(params['thermal_conductivity'].data, dtype=np.float32) + rho_cp = density * specific_heat + if np.any(rho_cp <= 0): + raise ValueError("params.density * params.specific_heat must be strictly positive everywhere.") + inv_rho_cp = 1.0 / rho_cp + + # Diagnostic: catch the common "no temperature rise" failure mode where + # the attenuation map is zero (all default openlifu Materials except + # STANDOFF have attenuation=0). Q = alpha * p^2 / (rho * c) is zero, so + # nothing will heat regardless of pulse handling. + q_peak = float(np.max(np.abs(Q_arr))) if Q_arr.size > 0 else 0.0 + if q_peak == 0: + att = params.get('attenuation') + att_is_zero = att is not None and float(np.max(np.abs(np.asarray(att.data)))) == 0 + logger.warning( + "All heat sources are zero%s; no temperature rise will be simulated. " + "Check that params['attenuation'] is nonzero in the region of interest.", + " (params['attenuation'] is zero everywhere)" if att_is_zero else "", + ) + elif pulse_duration is not None: + # Report the impulse-limit peak dT per pulse so the user can sanity + # check whether their setup is expected to produce a visible rise. + dT_peak_per_pulse = q_peak * float(pulse_duration) * float(np.max(inv_rho_cp)) + logger.info( + "Peak Q = %.3g W/m^3; impulse-limit peak dT per pulse = %.3g K " + "(pulse_duration = %g s).", + q_peak, dT_peak_per_pulse, pulse_duration, + ) + + dx, dy, dz = _get_dx(heat_sources_t.coords) + if oversample < 1: + raise ValueError("oversample must be >= 1.") + dstep = dt / oversample + cfl = _cfl_dt(dx, dy, dz, thermal_conductivity, rho_cp) + if dstep > cfl: + logger.warning( + "Requested internal step %g s exceeds explicit-diffusion CFL bound " + "%g s (dx=%g, dy=%g, dz=%g m). The simulation may be unstable; " + "consider increasing 'oversample' to at least %d.", + dstep, cfl, dx, dy, dz, max(1, math.ceil(dt / cfl)), + ) + + events_list = list(events) + if t_end is None: + if events_list: + t_end = max(e[0] for e in events_list) + dt + else: + t_end = dt + if t_end <= 0: + raise ValueError("t_end must be positive.") + N = round(t_end / dt) + 1 + + if mode in ('superpose', 'impulse') and pulse_duration is None: + raise ValueError(f"pulse_duration must be provided when mode='{mode}'.") + + t0_wall = time.perf_counter() + if mode == 'direct': + T_data = _run_direct( + Q_arr, events_list, N, dt, dstep, + dx, dy, dz, thermal_conductivity, inv_rho_cp, T_init, + ) + elif mode == 'superpose': + T_data = _run_superpose( + Q_arr, events_list, N, dt, dstep, + dx, dy, dz, thermal_conductivity, inv_rho_cp, T_init, + float(pulse_duration), + ) + elif mode == 'impulse': + T_data = _run_impulse( + Q_arr, events_list, N, dt, dstep, + dx, dy, dz, thermal_conductivity, inv_rho_cp, T_init, + float(pulse_duration), + ) + else: + raise ValueError( + f"Unknown mode '{mode}'. Expected 'direct', 'superpose', or 'impulse'." + ) + sim_time = time.perf_counter() - t0_wall + + coords_xyz = heat_sources_t.coords + t_coord = xa.DataArray( + np.arange(N) * dt, + dims=['t'], + attrs={'units': 's', 'long_name': 'Time'}, + ) + T_da = xa.DataArray( + T_data, + dims=['t', 'x', 'y', 'z'], + coords={ + 't': t_coord, + 'x': coords_xyz['x'], + 'y': coords_xyz['y'], + 'z': coords_xyz['z'], + }, + attrs={'units': 'degC', 'long_name': 'Temperature rise'}, + ) + Tmax = float(T_data.max()) + float(T0) + ds = xa.Dataset( + { + 'temperature_rise': T_da, + 'Tmax': xa.DataArray(Tmax, attrs={'units': 'degC', 'long_name': 'Peak temperature'}), + 'sim_time': xa.DataArray(sim_time, attrs={'units': 's', 'long_name': 'Wall-clock simulation time'}), + }, + attrs={'T0': float(T0), 'mode': mode}, + ) + logger.info("Thermal simulation complete in %.1f s (mode=%s, N=%d)", sim_time, mode, N) + return ds diff --git a/src/openlifu/util/annotations.py b/src/openlifu/util/annotations.py index 0cd2edda..ab8480d0 100644 --- a/src/openlifu/util/annotations.py +++ b/src/openlifu/util/annotations.py @@ -1,25 +1,66 @@ from __future__ import annotations -from typing import Annotated, NamedTuple +from dataclasses import dataclass, field +from typing import Tuple -class OpenLIFUFieldData(NamedTuple): +@dataclass(frozen=True) +class OpenLIFUFieldData: """ - A lightweight named tuple representing a name and annotation for the fields - of a dataclass. For example, the Graph dataclass may have fields associated - with this type: - - ```python - class Graph: - units: Annotated[str, OpenLIFUFieldData("Units", "The units of the graph")] = "mm" - dim_names: Annotated[ - Tuple[str, str, str], - OpenLIFUFieldData("Dimensions", "The name of the dimensions of the graph."), - ] = ("x", "y", "z") - ``` - - Annotated[] does not interfere with runtime behavior or type compatibility. + Lightweight metadata attached to a dataclass field via :class:`typing.Annotated`, + primarily consumed by GUI editors (e.g. SlicerOpenLIFU) to render fields with + human-friendly labels, units, and tooltips. + + Example:: + + class Pulse: + frequency: Annotated[ + float, + OpenLIFUFieldData( + name="Frequency", + description="Frequency of the pulse", + units="Hz", + display_units="kHz", + precision=1, + ), + ] = 400e3 + + The presence of ``Annotated[]`` does not affect runtime behavior or type + compatibility, and these fields are *not* serialized -- they describe how to + *display* the underlying value. The stored value remains in ``units``. + + Backwards compatibility: callers that historically constructed + ``OpenLIFUFieldData("Name", "Description")`` positionally continue to work + because ``name`` and ``description`` remain the first two fields and all + new fields are optional. + + Attributes: + name: Display label shown next to the field in editors. ``None`` falls + back to the dataclass field's attribute name. + description: Tooltip text shown on hover. ``None`` falls back to a generic + placeholder. + units: Storage units. The unit in which the underlying dataclass value + is stored (e.g. ``"Hz"``, ``"s"``, ``"Pa"``, ``"m"``, ``"deg"``). + ``None`` means no unit semantics are attached. + display_units: Preferred units for human display (e.g. ``"kHz"``, + ``"ms"``, ``"kPa"``, ``"mm"``). When set, editors should convert + from ``units`` to ``display_units`` for display, and back when + saving. ``None`` (default) means display in ``units``. + unit_options: Optional tuple of unit symbols a user may switch between + for display. Reserved for future use (units dropdown). + precision: Number of decimal places to display. ``None`` means use the + editor's default. + units_field: Optional name of a sibling dataclass field on the same + instance whose value provides the storage unit dynamically (e.g. + ``"distance_units"``). When set, editors should NOT auto-convert; + instead they should display the value as-is and label it with the + sibling's unit symbol. Mutually exclusive with ``units``. """ - name: Annotated[str | None, "The name of the dataclass field."] - description: Annotated[str | None, "The description of the dataclass field."] + name: str | None = None + description: str | None = None + units: str | None = None + display_units: str | None = None + unit_options: Tuple[str, ...] = field(default_factory=tuple) + precision: int | None = None + units_field: str | None = None diff --git a/src/openlifu/util/field_display.py b/src/openlifu/util/field_display.py new file mode 100644 index 00000000..b0678fd9 --- /dev/null +++ b/src/openlifu/util/field_display.py @@ -0,0 +1,215 @@ +from __future__ import annotations + +from dataclasses import is_dataclass +from typing import Any, Tuple, get_args, get_origin + +try: + # Python 3.10+ provides typing.get_type_hints with include_extras, but the + # rest of the codebase already uses ``get_type_hints`` from ``typing``. + from typing import get_type_hints as _get_type_hints +except ImportError: # pragma: no cover - defensive + from typing_extensions import get_type_hints as _get_type_hints + +from openlifu.util.annotations import OpenLIFUFieldData +from openlifu.util.units import getunitconversion + + +def get_field_metadata(cls: type, field_name: str) -> OpenLIFUFieldData | None: + """Return the :class:`OpenLIFUFieldData` annotation attached to ``cls.field_name``. + + Returns ``None`` if the field has no such annotation, if it has no + ``Annotated[]`` wrapping, or if the class is not a dataclass. + """ + if not is_dataclass(cls): + return None + try: + hints = _get_type_hints(cls, include_extras=True) + except Exception: + return None + annotated_type = hints.get(field_name) + if annotated_type is None: + return None + args = get_args(annotated_type) + if get_origin(annotated_type) is None or len(args) < 2: + return None + for meta in args[1:]: + if isinstance(meta, OpenLIFUFieldData): + return meta + return None + + +def resolve_units(meta: OpenLIFUFieldData, instance: Any) -> Tuple[str | None, str | None]: + """Return ``(storage_units, display_units)`` for ``meta`` applied to ``instance``. + + * Storage unit: ``instance.`` when ``units_field`` is set, + otherwise ``meta.units``. + * Display unit: ``meta.display_units`` when set, otherwise the storage unit. + + This means ``units_field`` and ``display_units`` may be combined: the value + is *stored* in whatever unit the sibling field reports, but the editor + *displays* it in the fixed ``display_units`` (e.g. always ``"mm"`` for a + distance field, regardless of the protocol's ``distance_units`` setting). + """ + if meta.units_field: + sibling = getattr(instance, meta.units_field, None) if instance is not None else None + else: + sibling = None + storage = sibling if meta.units_field else meta.units + display = meta.display_units or storage + return storage, display + + +def to_display(value: float, meta: OpenLIFUFieldData, instance: Any | None = None) -> float: + """Convert ``value`` from storage units to ``meta``'s display units. + + Falls back to the input value if conversion is not possible (no units + declared, or a unit-conversion failure).""" + if value is None: + return value + storage_unit, display_unit = resolve_units(meta, instance) + if storage_unit is None or display_unit is None or storage_unit == display_unit: + return value + try: + return float(value) * getunitconversion(storage_unit, display_unit) + except Exception: + return value + + +def from_display(value: float, meta: OpenLIFUFieldData, instance: Any | None = None) -> float: + """Inverse of :func:`to_display`: convert from display units to storage units.""" + if value is None: + return value + storage_unit, display_unit = resolve_units(meta, instance) + if storage_unit is None or display_unit is None or storage_unit == display_unit: + return value + try: + return float(value) * getunitconversion(display_unit, storage_unit) + except Exception: + return value + + +def format_value( + value: Any, + meta: OpenLIFUFieldData | None = None, + instance: Any | None = None, +) -> str: + """Format ``value`` using the precision/units in ``meta``. + + Numeric values are converted to ``meta.display_units`` (when known), + rounded to ``meta.precision`` (default 2 decimals for floats, no rounding + for ints), with trailing zeros trimmed, and the unit symbol appended. + Non-numeric values are passed through ``str()``. + """ + if value is None: + return "" + if meta is None: + if isinstance(value, float): + return _format_number(value, None) + return str(value) + + # Tuple/list: format element-wise, share unit suffix at the end + if isinstance(value, tuple | list): + formatted = [format_value(v, meta, instance) for v in value] + # Strip per-element unit so we don't repeat it; we'll add once at the end. + unit_suffix = _display_unit_suffix(meta, instance) + if unit_suffix: + stripped = [s[: -len(unit_suffix)].rstrip() if s.endswith(unit_suffix) else s for s in formatted] + return ", ".join(stripped) + " " + unit_suffix + return ", ".join(formatted) + + if isinstance(value, bool): + return "yes" if value else "no" + + if isinstance(value, int | float): + # Convert to the display unit. We always go through float so that, for + # an int value with display-unit conversion (e.g. an integer count of + # microns displayed in mm), we still get the proper scaled number. + display_value = to_display(float(value), meta, instance) + if isinstance(value, int) and (meta.display_units is None or meta.display_units == meta.units): + # No conversion needed for an int field; keep it integral. + display_value = value + precision = meta.precision + text = _format_number(display_value, precision) + suffix = _display_unit_suffix(meta, instance) + return f"{text} {suffix}" if suffix else text + + return str(value) + + +def _display_unit_suffix(meta: OpenLIFUFieldData, instance: Any | None) -> str: + _, display_unit = ( + resolve_units(meta, instance) if instance is not None else (meta.units, meta.display_units or meta.units) + ) + return display_unit or "" + + +def _strip_trailing_zeros(text: str) -> str: + """Remove trailing zeros (and a trailing dot) from a fixed-precision float string.""" + if "." not in text or "e" in text or "E" in text: + return text + stripped = text.rstrip("0").rstrip(".") + return stripped if stripped not in ("", "-") else "0" + + +def _format_number(value: Any, precision: int | None) -> str: + """Format a number with ``precision`` decimal places, stripping trailing zeros. + + Falls back to ``%g`` formatting whenever a fixed-precision render would + clip a non-zero value to ``"0"`` (for example, ``0.0025`` with + ``precision=1``). This mirrors the intent: don't print + ``300.0000000``, but also don't lose the entire value just because the + declared precision is too coarse for an unusually small number. + """ + if isinstance(value, bool): + return "yes" if value else "no" + if isinstance(value, int): + return str(value) + if not isinstance(value, float): + return str(value) + if precision is None: + return _strip_trailing_zeros(f"{value:.6g}") + fixed = f"{value:.{precision}f}" + fixed = _strip_trailing_zeros(fixed) + if fixed in ("0", "-0") and value != 0.0: + # The declared precision would erase the value entirely; use %g so + # the reader can still see the magnitude. + return _strip_trailing_zeros(f"{value:.6g}") + return fixed + + +def field_summary(instance: Any, field_name: str, label: str | None = None) -> str | None: + """Return ``"