|
| 1 | +# --- |
| 2 | +# jupyter: |
| 3 | +# jupytext: |
| 4 | +# text_representation: |
| 5 | +# extension: .py |
| 6 | +# format_name: percent |
| 7 | +# format_version: '1.3' |
| 8 | +# jupytext_version: 1.19.1 |
| 9 | +# kernelspec: |
| 10 | +# display_name: env (3.11.4) |
| 11 | +# language: python |
| 12 | +# name: python3 |
| 13 | +# --- |
| 14 | + |
| 15 | +# %% [markdown] |
| 16 | +# # 03: Solution Generation and Analysis |
| 17 | +# |
| 18 | +# This notebook demonstrates how to generate an acoustic `Solution` using a `Protocol`, `Target`, and `Transducer`. It also covers running acoustic simulations and analyzing the results, including checking against safety parameters. |
| 19 | +# |
| 20 | +# We will build upon the concepts from Notebook 01 (Object Creation) and Notebook 02 (Database Interaction/Transducer Loading). |
| 21 | + |
| 22 | +# %% [markdown] |
| 23 | +# ## 1. Setup and Imports |
| 24 | +# |
| 25 | +# First, let's import the necessary classes. |
| 26 | + |
| 27 | +# %% |
| 28 | +from __future__ import annotations |
| 29 | + |
| 30 | +import numpy as np |
| 31 | + |
| 32 | +# For displaying tables |
| 33 | +# Core OpenLIFU objects |
| 34 | +from openlifu.bf import Pulse, Sequence, apod_methods, focal_patterns |
| 35 | +from openlifu.geo import Point |
| 36 | +from openlifu.plan import Protocol |
| 37 | +from openlifu.plan.param_constraint import ParameterConstraint |
| 38 | +from openlifu.sim import SimSetup |
| 39 | +from openlifu.xdc import Transducer |
| 40 | + |
| 41 | + |
| 42 | +# %% [markdown] |
| 43 | +# ## 2. Defining Components for Solution Calculation |
| 44 | +# |
| 45 | +# To calculate a solution, we need: |
| 46 | +# 1. A `Target`: Where we want to focus the ultrasound. |
| 47 | +# 2. A `Transducer`: The physical device that will generate the ultrasound. |
| 48 | +# 3. A `Protocol`: The "recipe" defining pulse characteristics, sequence, focal pattern, etc. |
| 49 | + |
| 50 | +# %% [markdown] |
| 51 | +# ### 2.1. Define a Target |
| 52 | +# This is a `Point` object representing the desired focal location. |
| 53 | + |
| 54 | +# %% |
| 55 | +target = Point(position=np.array([0, 0, 50]), units="mm", radius=0.5) # 50mm depth, small radius |
| 56 | +print(f"Target: {target}") |
| 57 | + |
| 58 | +# %% [markdown] |
| 59 | +# ### 2.2. Load or Define a Transducer |
| 60 | +# For this example, we'll first try to load a transducer from the database, as shown in Notebook 02. |
| 61 | +# If that fails (e.g., database not found), we'll fall back to a programmatically generated generic transducer for demonstration purposes. |
| 62 | + |
| 63 | +# %% |
| 64 | +transducer = Transducer.gen_matrix_array( |
| 65 | + nx=16, ny=16, pitch=3, kerf=0.1, |
| 66 | + id="generic_16x16", name="Generic 16x16 Array", |
| 67 | + sensitivity=2000, # Pa/V |
| 68 | +) |
| 69 | + |
| 70 | +print(f"Using Transducer: {transducer.id}, Number of elements: {transducer.numelements()}") |
| 71 | + |
| 72 | +# %% [markdown] |
| 73 | +# ### 2.3. Define a Protocol |
| 74 | +# The protocol specifies *how* the sonication should be performed. |
| 75 | + |
| 76 | +# %% |
| 77 | +# Pulse definition |
| 78 | +f0 = 400e3 # Use transducer f0 if available |
| 79 | +pulse = Pulse(frequency=f0, duration=10e-3) # 10 cycles duration |
| 80 | + |
| 81 | +# Sequence definition |
| 82 | +sequence = Sequence(pulse_interval=100e-3, pulse_count=9, pulse_train_interval=0, pulse_train_count=1) |
| 83 | + |
| 84 | +# Focal Pattern: Let's use a SinglePoint focus for this example. |
| 85 | +# The actual target point is provided during calc_solution. |
| 86 | +# target_pressure is an optional parameter for scaling. |
| 87 | +focal_pattern = focal_patterns.SinglePoint(target_pressure=1.0e6) # Target 1 MPa |
| 88 | + |
| 89 | +# Apodization Method |
| 90 | +apod_method = apod_methods.MaxAngle(max_angle=30) # Limit elements to a 30-degree cone |
| 91 | + |
| 92 | +# Simulation Setup: Defines the grid for acoustic simulation |
| 93 | +sim_setup = SimSetup( |
| 94 | + x_extent=(-25, 25), y_extent=(-25, 25), z_extent=(-5, 70), # in mm |
| 95 | + spacing=1.0 # 1 mm resolution |
| 96 | +) |
| 97 | + |
| 98 | +# Create the Protocol object |
| 99 | +protocol1 = Protocol( |
| 100 | + id='example_protocol_prog', |
| 101 | + name='Example Protocol (Programmatic)', |
| 102 | + pulse=pulse, |
| 103 | + sequence=sequence, |
| 104 | + focal_pattern=focal_pattern, # Store the type |
| 105 | + apod_method=apod_method, |
| 106 | + sim_setup=sim_setup |
| 107 | +) |
| 108 | +print(f"Defined Protocol: {protocol1.name}") |
| 109 | + |
| 110 | +# %% [markdown] |
| 111 | +# ## 3. Calculating the Solution |
| 112 | +# |
| 113 | +# With the `target`, `transducer`, and `protocol` defined, we can now calculate the `Solution`. |
| 114 | +# The `calc_solution` method returns: |
| 115 | +# * `solution`: The `Solution` object containing delays, apodizations, voltage, etc. |
| 116 | +# * `sim_res`: An `xa.Dataset` simulation result if `simulate=True`, aggregated over the focal points. |
| 117 | +# * `analysis`: A `SolutionAnalysis` object if `simulate=True`. |
| 118 | +# |
| 119 | +# Setting `simulate=True` will run an acoustic simulation. |
| 120 | +# Setting `scale=True` will attempt to scale the output pressure to match `target_pressure` defined in the focal pattern or protocol, and returns a `scaled_analysis`. |
| 121 | + |
| 122 | +# %% |
| 123 | +print(f"\nCalculating solution for protocol '{protocol1.name}' and target '{target.name}'...") |
| 124 | +solution1, sim_res1, analysis1 = protocol1.calc_solution( |
| 125 | + target=target, |
| 126 | + transducer=transducer, |
| 127 | + simulate=True, |
| 128 | + scale=True # Try to scale to target_pressure |
| 129 | +) |
| 130 | + |
| 131 | +print(f"\nSolution calculated: {solution1.id}") |
| 132 | +print(f" Calculated Voltage: {solution1.voltage:.2f} V (this is a relative/normalized value before hardware calibration)") |
| 133 | +# print(f" Delays (first 5 elements): {solution1.delays[0, :5]}") |
| 134 | +# print(f" Apodizations (first 5 elements): {solution1.apodizations[0, :5]}") |
| 135 | + |
| 136 | +print("\nSimulation Result object created.") |
| 137 | +# sim_res1 contains the raw simulation grid and pressure data. |
| 138 | +# For example, to get the peak pressure and its location: |
| 139 | +# peak_pressure_Pa, peak_loc_mm = sim_res1.get_peak_pressure() |
| 140 | +# print(f" Peak pressure in simulation: {peak_pressure_Pa/1e6:.2f} MPa at {peak_loc_mm} mm") |
| 141 | +# sim_res1.plot_slices() # This would require matplotlib and a GUI backend |
| 142 | + |
| 143 | +print("\nSolutionAnalysis object created (scaled):") |
| 144 | +# The SolutionAnalysis object provides various calculated acoustic parameters. |
| 145 | +# We can display it as a table: |
| 146 | +analysis_table = analysis1.to_table() |
| 147 | +analysis_table.set_index('Param')[['Value', 'Units', 'Status']] |
| 148 | + |
| 149 | + |
| 150 | +# %% |
| 151 | +solution1.simulation_result['p_min'].sel(focal_point_index=0).sel(y=0).plot.imshow() |
| 152 | + |
| 153 | +# %% [markdown] |
| 154 | +# We can also run the simulation separately if needed, by calling the `simulate` method on the `Solution` object. This also allows for more control over simulation parameters, such as modifying apodizations, delays, or the acoustic parameters being used for the simulation. |
| 155 | + |
| 156 | +# %% |
| 157 | +solution2, sim_res2, analysis2 = protocol1.calc_solution( |
| 158 | + target=target, |
| 159 | + transducer=transducer, |
| 160 | + simulate=False, |
| 161 | + scale=False # Try to scale to target_pressure |
| 162 | +) |
| 163 | +solution2.apodizations[0,:128] = 0 |
| 164 | +params = protocol1.sim_setup.setup_sim_scene(protocol1.seg_method) |
| 165 | +simulation_result = solution2.simulate(params=params) |
| 166 | + |
| 167 | +# %% |
| 168 | +simulation_result['p_min'].sel(focal_point_index=0).sel(y=0).plot.imshow() |
| 169 | +solution2.analyze(simulation_result).to_table().set_index('Param')[['Value', 'Units', 'Status']] |
| 170 | + |
| 171 | +# %% [markdown] |
| 172 | +# ## 4. Using Parameter Constraints in Analysis |
| 173 | +# |
| 174 | +# We can define constraints for various parameters (like MI, TIC, Isppa) and see if the solution meets them. |
| 175 | + |
| 176 | +# %% |
| 177 | +# Define some example parameter constraints |
| 178 | +constraints = { |
| 179 | + "MI": ParameterConstraint('<', 1.8, 1.85), # Mechanical Index should be < 1.8 (error if > 1.85) |
| 180 | + "TIC": ParameterConstraint('<', 2.0), # Thermal Index (cranial) should be < 2.0 |
| 181 | + "global_isppa_Wcm2": ParameterConstraint('within', error_value=(50, 200)) # Isppa between 50-200 W/cm^2 |
| 182 | +} |
| 183 | + |
| 184 | +print("\nAnalysis table with constraints:") |
| 185 | +# The to_table method can accept these constraints directly |
| 186 | +constrained_table = analysis1.to_table(constraints=constraints) |
| 187 | +constrained_table.set_index('Param')[['Value', 'Units', 'Status']] |
| 188 | + |
| 189 | + |
| 190 | +# %% [markdown] |
| 191 | +# ## Summary and Next Steps |
| 192 | +# |
| 193 | +# This notebook showed how to: |
| 194 | +# * Define or load the necessary components (`Target`, `Transducer`, `Protocol`). |
| 195 | +# * Calculate a `Solution` using `protocol.calc_solution()`. |
| 196 | +# * Enable acoustic simulation and obtain `SimResult` and `SolutionAnalysis` objects. |
| 197 | +# * Use `ParameterConstraint` to evaluate the solution against safety or performance criteria. |
| 198 | +# |
| 199 | +# The `Solution` object is key for hardware interaction. The next notebook, `04_Connecting_to_Hardware.py`, will introduce how to establish communication with OpenLIFU hardware. Following that, `05_Solution_to_Hardware_Basic.py` will demonstrate sending a calculated solution to the device. |
| 200 | + |
| 201 | +# %% [markdown] |
| 202 | +# End of Notebook 03 |
0 commit comments