Skip to content

Commit 3b76420

Browse files
peterhollenderpeterhollender
authored
enhance analysis for #334
fix pulse train interval type and make text consistent
1 parent 7f1505b commit 3b76420

File tree

2 files changed

+98
-50
lines changed

2 files changed

+98
-50
lines changed

‎src/openlifu/plan/solution.py‎

Lines changed: 31 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
from openlifu.plan.solution_analysis import (
1818
SolutionAnalysis,
1919
SolutionAnalysisOptions,
20+
find_centroid,
2021
get_beamwidth,
2122
get_mask,
2223
)
@@ -172,8 +173,12 @@ def analyze(self,
172173
# xyz = np.stack(np.meshgrid(*coords, indexing="xy"), axis=-1) #TODO: if fus.Axis is defined, coords.ndgrid(dim="ax")
173174
# z_mask = xyz[..., -1] >= options.sidelobe_zmin #TODO: probably wrong here, should be z{1}>=options.sidelobe_zmin;
174175

175-
pulsetrain_dutycycle = self.get_pulsetrain_dutycycle()
176-
treatment_dutycycle = self.get_treatment_dutycycle()
176+
solution_analysis.duty_cycle_pulse_train_pct = self.get_pulsetrain_dutycycle()*100
177+
solution_analysis.duty_cycle_sequence_pct = self.get_sequence_dutycycle()*100
178+
if self.sequence.pulse_train_interval == 0:
179+
solution_analysis.sequence_duration_s = float(self.sequence.pulse_interval * self.sequence.pulse_count * self.sequence.pulse_train_count)
180+
else:
181+
solution_analysis.sequence_duration_s = float(self.sequence.pulse_train_interval * self.sequence.pulse_train_count)
177182
ita_mWcm2 = rescale_coords(self.get_ita(units="mW/cm^2"), options.distance_units)
178183

179184
power_W = np.zeros(self.num_foci())
@@ -182,6 +187,10 @@ def analyze(self,
182187
pnp_MPa = pnp_MPa_all.isel(focal_point_index=focus_index)
183188
ipa_Wcm2 = ipa_Wcm2_all.isel(focal_point_index=focus_index)
184189
focus = self.foci[focus_index].get_position(units=options.distance_units)
190+
focus_mm = self.foci[focus_index].get_position(units="mm")
191+
solution_analysis.target_position_lat_mm += [focus_mm[0]]
192+
solution_analysis.target_position_ele_mm += [focus_mm[1]]
193+
solution_analysis.target_position_ax_mm += [focus_mm[2]] # TODO: this should be a list, not a single value
185194
apodization = self.apodizations[focus_index]
186195
origin = transducer.get_effective_origin(apodizations=apodization, units=options.distance_units)
187196

@@ -193,16 +202,6 @@ def analyze(self,
193202

194203
p0_Pa = np.max(output_signal_Pa, axis=1)
195204

196-
# get focus region masks (for mainlobe, sidelobe and beamwidth)
197-
# mainlobe_mask = mask_focus(
198-
# self.simulation_result,
199-
# foc,
200-
# options.mainlobe_radius,
201-
# mask_op=MaskOp.LESS_EQUAL,
202-
# units=options.distance_units,
203-
# aspect_ratio=options.mainlobe_aspect_ratio
204-
# )
205-
206205
mainlobe_mask = get_mask(
207206
pnp_MPa,
208207
focus = focus,
@@ -223,7 +222,13 @@ def analyze(self,
223222
z_mask = pnp_MPa.coords['ax'] > options.sidelobe_zmin
224223
sidelobe_mask = sidelobe_mask.where(z_mask, False)
225224

226-
pk = float(pnp_MPa.where(mainlobe_mask).max())
225+
pnp_mainlobe = pnp_MPa.where(mainlobe_mask)
226+
pk = float(pnp_mainlobe.max())
227+
mainlobe_focus = find_centroid(pnp_mainlobe, pk*10**(-3/20), "mm")
228+
solution_analysis.focal_centroid_lat_mm += [mainlobe_focus[0]]
229+
solution_analysis.focal_centroid_ele_mm += [mainlobe_focus[1]]
230+
solution_analysis.focal_centroid_ax_mm += [mainlobe_focus[2]]
231+
227232
solution_analysis.mainlobe_pnp_MPa += [pk]
228233

229234
for dim, scale in zip(pnp_MPa.dims, options.mainlobe_aspect_ratio):
@@ -250,15 +255,16 @@ def analyze(self,
250255
solution_analysis.global_pnp_MPa += [float(pnp_MPa.where(z_mask).max())]
251256
solution_analysis.global_isppa_Wcm2 += [float(ipa_Wcm2.where(z_mask).max())]
252257
i0_Wcm2 = (p0_Pa**2 / (2*standoff_Z)) * 1e-4
253-
i0ta_Wcm2 = i0_Wcm2 * pulsetrain_dutycycle * treatment_dutycycle
258+
i0ta_Wcm2 = i0_Wcm2 * solution_analysis.duty_cycle_sequence_pct/100
254259
power_W[focus_index] = np.mean(np.sum(i0ta_Wcm2 * ele_sizes_cm2 * self.apodizations[focus_index, :]))
255260
TIC[focus_index] = power_W[focus_index] / (d_eq_cm * c_tic)
256261
solution_analysis.p0_MPa += [1e-6*np.max(p0_Pa)]
262+
solution_analysis.global_ispta_mWcm2 = float((ita_mWcm2*z_mask).max())
263+
solution_analysis.MI = (np.max(solution_analysis.mainlobe_pnp_MPa)/np.sqrt(self.pulse.frequency*1e-6))
257264
solution_analysis.TIC = np.mean(TIC)
258-
solution_analysis.power_W = np.mean(power_W)
259265
solution_analysis.voltage_V = self.voltage
260-
solution_analysis.MI = (np.max(solution_analysis.mainlobe_pnp_MPa)/np.sqrt(self.pulse.frequency*1e-6))
261-
solution_analysis.global_ispta_mWcm2 = float((ita_mWcm2*z_mask).max())
266+
solution_analysis.power_W = np.mean(power_W)
267+
262268
solution_analysis.param_constraints = param_constraints
263269
return solution_analysis
264270

@@ -323,7 +329,7 @@ def scale(
323329

324330
def get_pulsetrain_dutycycle(self) -> float:
325331
"""
326-
Compute the pulsetrain dutycycle given a sequence.
332+
Compute the pulse train dutycycle given a sequence.
327333
328334
Returns:
329335
A float.
@@ -332,19 +338,19 @@ def get_pulsetrain_dutycycle(self) -> float:
332338

333339
return pulsetrain_dutycycle
334340

335-
def get_treatment_dutycycle(self) -> float:
341+
def get_sequence_dutycycle(self) -> float:
336342
"""
337-
Compute the treatment dutycycle given a sequence.
343+
Compute the overall sequence dutycycle given a sequence.
338344
339345
Returns:
340346
A float.
341347
"""
342348
if self.sequence.pulse_train_interval == 0:
343-
treatment_dutycycle = 1
349+
between_pulsetrain_duty_cycle = 1
344350
else:
345-
treatment_dutycycle = (self.sequence.pulse_count * self.sequence.pulse_interval) / self.sequence.pulse_train_interval
346-
347-
return treatment_dutycycle
351+
between_pulsetrain_duty_cycle = (self.sequence.pulse_count * self.sequence.pulse_interval) / self.sequence.pulse_train_interval
352+
sequence_duty_cycle = self.get_pulsetrain_dutycycle() * between_pulsetrain_duty_cycle
353+
return sequence_duty_cycle
348354

349355
def get_ita(self, units: str = "mW/cm^2") -> xa.DataArray:
350356
"""
@@ -360,7 +366,7 @@ def get_ita(self, units: str = "mW/cm^2") -> xa.DataArray:
360366
"""
361367
intensity_scaled = rescale_data_arr(self.simulation_result['intensity'], units)
362368
pulsetrain_dutycycle = self.get_pulsetrain_dutycycle()
363-
treatment_dutycycle = self.get_treatment_dutycycle()
369+
treatment_dutycycle = self.get_sequence_dutycycle()
364370
pulse_seq = (np.arange(self.sequence.pulse_count) - 1) % self.num_foci() + 1
365371
counts = np.zeros((1, 1, 1, self.num_foci()))
366372
for i in range(self.num_foci()):

‎src/openlifu/plan/solution_analysis.py‎

Lines changed: 67 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -11,30 +11,39 @@
1111
from openlifu.plan.param_constraint import PARAM_STATUS_SYMBOLS, ParameterConstraint
1212
from openlifu.util.annotations import OpenLIFUFieldData
1313
from openlifu.util.dict_conversion import DictMixin
14-
from openlifu.util.units import getunittype
14+
from openlifu.util.units import getunitconversion, getunittype
1515

1616
DEFAULT_ORIGIN = np.zeros(3)
1717

1818
PARAM_FORMATS = {
1919
"mainlobe_pnp_MPa": ["max", "0.3f", "MPa", "Mainlobe Peak Negative Pressure"],
20-
"mainlobe_isppa_Wcm2": ["max", "0.3f", "W/cm^2", "Mainlobe I_SPPA"],
21-
"mainlobe_ispta_mWcm2": ["mean", "0.3f", "mW/cm^2", "Mainlobe I_SPTA"],
22-
"beamwidth_lat_3dB_mm": ["mean", "0.3f", "mm", "3dB Lateral Beamwidth"],
23-
"beamwidth_ele_3dB_mm": ["mean", "0.3f", "mm", "3dB Elevational Beamwidth"],
24-
"beamwidth_ax_3dB_mm": ["mean", "0.3f", "mm", "3dB Axial Beamwidth"],
25-
"beamwidth_lat_6dB_mm": ["mean", "0.3f", "mm", "6dB Lateral Beamwidth"],
26-
"beamwidth_ele_6dB_mm": ["mean", "0.3f", "mm", "6dB Elevational Beamwidth"],
27-
"beamwidth_ax_6dB_mm": ["mean", "0.3f", "mm", "6dB Axial Beamwidth"],
20+
"mainlobe_isppa_Wcm2": ["max", "0.1f", "W/cm^2", "Mainlobe I_SPPA"],
21+
"mainlobe_ispta_mWcm2": ["mean", "0.1f", "mW/cm^2", "Mainlobe I_SPTA"],
22+
"target_position_lat_mm": ["mean", "0.1f", "mm", "Target Position (Lateral)"],
23+
"target_position_ele_mm": ["mean", "0.1f", "mm", "Target Position (Elevation)"],
24+
"target_position_ax_mm": ["mean", "0.1f", "mm", "Target Position (Axial)"],
25+
"focal_centroid_lat_mm": ["mean", "0.1f", "mm", "Focal Centroid (Lateral)"],
26+
"focal_centroid_ele_mm": ["mean", "0.1f", "mm", "Focal Centroid (Elevation)"],
27+
"focal_centroid_ax_mm": ["mean", "0.1f", "mm", "Focal Centroid (Axial)"],
28+
"beamwidth_lat_3dB_mm": ["mean", "0.2f", "mm", "3dB Beamwidth (Lateral)"],
29+
"beamwidth_ele_3dB_mm": ["mean", "0.2f", "mm", "3dB Beamwidth (Elevational)"],
30+
"beamwidth_ax_3dB_mm": ["mean", "0.2f", "mm", "3dB Beamwidth (Axial)"],
31+
"beamwidth_lat_6dB_mm": ["mean", "0.2f", "mm", "6dB Beamwidth (Lateral)"],
32+
"beamwidth_ele_6dB_mm": ["mean", "0.2f", "mm", "6dB Beamwidth (Elevational)"],
33+
"beamwidth_ax_6dB_mm": ["mean", "0.2f", "mm", "6dB Beamwidth (Axial)"],
2834
"sidelobe_pnp_MPa": ["max", "0.3f", "MPa", "Sidelobe Peak Negative Pressure"],
29-
"sidelobe_isppa_Wcm2": ["max", "0.3f", "W/cm^2", "Sidelobe I_SPPA"],
35+
"sidelobe_isppa_Wcm2": ["max", "0.1f", "W/cm^2", "Sidelobe I_SPPA"],
3036
"global_pnp_MPa": ["max", "0.3f", "MPa", "Global Peak Negative Pressure"],
31-
"global_isppa_Wcm2": ["max", "0.3f", "W/cm^2", "Global I_SPPA"],
32-
"global_ispta_mWcm2": [None, "0.3f", "mW/cm^2", "Global I_SPTA"],
37+
"global_isppa_Wcm2": ["max", "0.1f", "W/cm^2", "Global I_SPPA"],
38+
"global_ispta_mWcm2": [None, "0.1f", "mW/cm^2", "Global I_SPTA"],
39+
"MI": [None, "0.2f", "", "MI"],
40+
"TIC": [None, "0.2f", "", "TIC"],
41+
"voltage_V": [None, "0.1f", "V", "Voltage"],
3342
"p0_MPa": ["max", "0.3f", "MPa", "Emitted Pressure"],
34-
"power_W": [None, "0.3f", "W", "Emitted Power"],
35-
"voltage_V": [None, "0.3f", "V", "Voltage"],
36-
"TIC": [None, "0.3f", "", "TIC"],
37-
"MI": [None, "0.3f", "", "MI"]}
43+
"power_W": [None, "0.2f", "W", "Emitted Power"],
44+
"duty_cycle_pulse_train_pct": [None, "0.1f", "%", "Pulse Train Duty Cycle"],
45+
"duty_cycle_sequence_pct": [None, "0.1f", "%", "Sequence Duty Cycle"],
46+
"sequence_duration_s": [None, "0.0f", "s", "Sequence Duration"],}
3847

3948
@dataclass
4049
class SolutionAnalysis(DictMixin):
@@ -47,6 +56,24 @@ class SolutionAnalysis(DictMixin):
4756
mainlobe_ispta_mWcm2: Annotated[list[float], OpenLIFUFieldData("Mainlobe ISPTA", "Spatial peak time average intensity in the mainlobe, in mW/cm²")] = field(default_factory=list)
4857
"""Spatial peak time average intensity in the mainlobe, in mW/cm²"""
4958

59+
target_position_lat_mm: Annotated[list[float], OpenLIFUFieldData("Target position lateral (mm)", "Lateral position of the target, in mm")] = field(default_factory=list)
60+
"""Lateral position of the target, in mm"""
61+
62+
target_position_ele_mm: Annotated[list[float], OpenLIFUFieldData("Target position elevation (mm)", "Elevation position of the target, in mm")] = field(default_factory=list)
63+
"""Elevation position of the target, in mm"""
64+
65+
target_position_ax_mm: Annotated[list[float], OpenLIFUFieldData("Target position axial (mm)", "Axial position of the target, in mm")] = field(default_factory=list)
66+
"""Axial position of the target, in mm"""
67+
68+
focal_centroid_lat_mm: Annotated[list[float], OpenLIFUFieldData("Focal centroid lateral (mm)", "Lateral centroid of the focus, in mm")] = field(default_factory=list)
69+
"""Lateral centroid of the focus, in mm"""
70+
71+
focal_centroid_ele_mm: Annotated[list[float], OpenLIFUFieldData("Focal centroid elevation (mm)", "Elevation centroid of the focus, in mm")] = field(default_factory=list)
72+
"""Elevation centroid of the focus, in mm"""
73+
74+
focal_centroid_ax_mm: Annotated[list[float], OpenLIFUFieldData("Focal centroid axial (mm)", "Axial centroid of the focus, in mm")] = field(default_factory=list)
75+
"""Axial centroid of the focus, in mm"""
76+
5077
beamwidth_lat_3dB_mm: Annotated[list[float], OpenLIFUFieldData("3dB lateral beamwidth", "Lateral beamwidth at -3 dB, in mm")] = field(default_factory=list)
5178
"""Lateral beamwidth at -3 dB, in mm"""
5279

@@ -77,23 +104,33 @@ class SolutionAnalysis(DictMixin):
77104
global_isppa_Wcm2: Annotated[list[float], OpenLIFUFieldData("Global ISPPA", "Maximum spatial peak pulse average intensity in the entire field, in W/cm²")] = field(default_factory=list)
78105
"""Maximum spatial peak pulse average intensity in the entire field, in W/cm²"""
79106

80-
p0_MPa: Annotated[list[float], OpenLIFUFieldData("Emitted pressure (MPa)", "Initial pressure values in the field, in MPa")] = field(default_factory=list)
81-
"""Initial pressure values in the field (MPa)"""
107+
global_ispta_mWcm2: Annotated[float | None, OpenLIFUFieldData("Global ISPTA (mW/cm²)", "Global Intensity at Spatial-Peak, Time-Average (I_SPTA) (mW/cm²)")] = None
108+
"""Global Intensity at Spatial-Peak, Time-Average (I_SPTA) (mW/cm²)"""
109+
110+
MI: Annotated[float | None, OpenLIFUFieldData("Mechanical index (MI)", "Mechanical index (MI)")] = None
111+
"""Mechanical index (MI)"""
82112

83113
TIC: Annotated[float | None, OpenLIFUFieldData("Thermal index (TIC)", "Thermal index in cranium (TIC)")] = None
84114
"""Thermal index in cranium (TIC)"""
85115

116+
voltage_V: Annotated[float | None, OpenLIFUFieldData("Voltage (V)", "Voltage applied to the transducer (V)")] = None
117+
"""Voltage applied to the transducer (V)"""
118+
119+
p0_MPa: Annotated[list[float], OpenLIFUFieldData("Emitted pressure (MPa)", "Initial pressure values in the field, in MPa")] = field(default_factory=list)
120+
"""Initial pressure values in the field (MPa)"""
121+
86122
power_W: Annotated[float | None, OpenLIFUFieldData("Emitted Power (W)", "Emitted power from the transducer face (W)")] = None
87123
"""Emitted power from the transducer face (W)"""
88124

89-
voltage_V: Annotated[float | None, OpenLIFUFieldData("Voltage (V)", "Voltage applied to the transducer (V)")] = None
90-
"""Voltage applied to the transducer (V)"""
125+
duty_cycle_pulse_train_pct: Annotated[float | None, OpenLIFUFieldData("Duty cycle pulse train", "Duty cycle within a pulse train (0-1)")] = None
126+
"""Duty cycle within a pulse train (0-1)"""
91127

92-
MI: Annotated[float | None, OpenLIFUFieldData("Mechanical index (MI)", "Mechanical index (MI)")] = None
93-
"""Mechanical index (MI)"""
128+
duty_cycle_sequence_pct: Annotated[float | None, OpenLIFUFieldData("Duty cycle sequence", "Duty cycle of the overall sequence (0-1)")] = None
129+
"""Duty cycle of the overall sequence (0-1)"""
130+
131+
sequence_duration_s: Annotated[float | None, OpenLIFUFieldData("Sequence duration (s)", "Total duration of the sequence (s)")] = None
132+
"""Total duration of the sequence (s)"""
94133

95-
global_ispta_mWcm2: Annotated[float | None, OpenLIFUFieldData("Global ISPTA (mW/cm²)", "Global Intensity at Spatial-Peak, Time-Average (I_SPTA) (mW/cm²)")] = None
96-
"""Global Intensity at Spatial-Peak, Time-Average (I_SPTA) (mW/cm²)"""
97134

98135
param_constraints: Annotated[Dict[str, ParameterConstraint], OpenLIFUFieldData("Parameter constraints", None)] = field(default_factory=dict)
99136
"""TODO: Add description"""
@@ -258,12 +295,17 @@ def from_dict(cls: Type[SolutionAnalysisOptions], parameter_dict: Dict[str, Any]
258295

259296
return cls(**parameter_dict)
260297

261-
def find_centroid(da: xa.DataArray, cutoff:float) -> np.ndarray:
298+
def find_centroid(da: xa.DataArray, cutoff:float, units:None) -> np.ndarray:
262299
"""Find the centroid of a thresholded region of a DataArray"""
300+
if units is not None and getunittype(units) != 'distance':
301+
raise ValueError(f"Units must be a length unit, got {units}")
263302
da = da.where(da > cutoff, 0)
264303
dims = list(da.dims)
265304
coords = np.meshgrid(*[da.coords[coord] for coord in dims], indexing='ij')
266305
centroid = np.array([np.sum(da*coords[dims.index(dim)])/da.sum() for dim in da.dims])
306+
if units is not None:
307+
da_units = [da.coords[dim].attrs.get('units', None) for dim in dims]
308+
centroid = np.array([getunitconversion(coord_units, units) * c for coord_units, c in zip(da_units, centroid)])
267309
return centroid
268310

269311
def get_focus_matrix(focus, origin=[0,0,0]) -> np.ndarray:

0 commit comments

Comments
 (0)