Skip to content

Switch to k-Wave-python #435

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
faberno wants to merge 2 commits into
base: develop
Choose a base branch
from
Draft

Switch to k-Wave-python #435

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d7fba58
switch from k-Wave (matlab) to k-Wave-python
faberno Dec 16, 2025
00614d8
fix double logs and docstrings
faberno Dec 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ dependencies = [
"pre-commit>=3.2.2", # Uses MIT-License (MIT compatible)
"PyWavelets", # Uses MIT-License (MIT compatible)
"scikit-learn>=1.1.0", # Uses BSD-License (MIT compatible)
"k-wave-python==0.4.0" # Used LGPL-3.0 License (MIT compatible)
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo in comment: "Used LGPL-3.0 License" should be "Uses LGPL-3.0 License" to match the pattern of other license comments in this file.

Suggested change
"k-wave-python==0.4.0" # Used LGPL-3.0 License (MIT compatible)
"k-wave-python==0.4.0" # Uses LGPL-3.0 License (MIT compatible)

Copilot uses AI. Check for mistakes.
]

[project.optional-dependencies]
Expand Down
182 changes: 148 additions & 34 deletions simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,16 @@
# SPDX-FileCopyrightText: 2021 Janek Groehl
# SPDX-License-Identifier: MIT

import gc
import os
import subprocess

from kwave.utils.kwave_array import kWaveArray
from kwave.kgrid import kWaveGrid
from kwave.kmedium import kWaveMedium
from kwave.ksource import kSource
from kwave.ksensor import kSensor
from kwave.options.simulation_execution_options import SimulationExecutionOptions
from kwave.options.simulation_options import SimulationOptions
from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D
from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D
import numpy as np
import scipy.io as sio
from scipy.spatial.transform import Rotation

from simpa.core.device_digital_twins import (CurvedArrayDetectionGeometry,
Expand All @@ -16,7 +20,6 @@
AcousticAdapterBase
from simpa.io_handling.io_hdf5 import load_data_field, save_hdf5
from simpa.utils import Tags
from simpa.utils.matlab import generate_matlab_cmd
from simpa.utils.calculate import rotation_matrix_between_vectors
from simpa.utils.dict_path_manager import generate_dict_path
from simpa.utils.path_manager import PathManager
Expand Down Expand Up @@ -197,9 +200,6 @@ def k_wave_acoustic_forward_model(self, detection_geometry: DetectionGeometryBas
data_dict[Tags.KWAVE_PROPERTY_DIRECTIVITY_ANGLE] = angles
data_dict[Tags.KWAVE_PROPERTY_INTRINSIC_EULER_ANGLE] = intrinsic_euler_angles

optical_path = optical_path + ".mat"
optical_path = os.path.abspath(optical_path)

possible_k_wave_parameters = [Tags.SPACING_MM, Tags.MODEL_SENSOR_FREQUENCY_RESPONSE,
Tags.KWAVE_PROPERTY_ALPHA_POWER, Tags.GPU, Tags.KWAVE_PROPERTY_PMLInside, Tags.KWAVE_PROPERTY_PMLAlpha, Tags.KWAVE_PROPERTY_PlotPML,
Tags.RECORDMOVIE, Tags.MOVIENAME, Tags.ACOUSTIC_LOG_SCALE,
Expand All @@ -210,7 +210,8 @@ def k_wave_acoustic_forward_model(self, detection_geometry: DetectionGeometryBas
Tags.DETECTOR_ELEMENT_WIDTH_MM: pa_device.detector_element_width_mm,
Tags.SENSOR_CENTER_FREQUENCY_HZ: pa_device.center_frequency_Hz,
Tags.SENSOR_BANDWIDTH_PERCENT: pa_device.bandwidth_percent,
Tags.SENSOR_SAMPLING_RATE_MHZ: pa_device.sampling_frequency_MHz
Tags.SENSOR_SAMPLING_RATE_MHZ: pa_device.sampling_frequency_MHz,
Tags.ACOUSTIC_MODEL_BINARY_PATH: self.component_settings[Tags.ACOUSTIC_MODEL_BINARY_PATH],
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tags.ACOUSTIC_MODEL_BINARY_PATH is added to k_wave_settings but doesn't appear to be used in the run_k_wave_simulation method. If the k-wave-python library doesn't use this path directly through the API (it's not passed to SimulationExecutionOptions), this line may be unnecessary. Consider removing it or documenting how the binary path is intended to be used with k-wave-python.

Suggested change
Tags.ACOUSTIC_MODEL_BINARY_PATH: self.component_settings[Tags.ACOUSTIC_MODEL_BINARY_PATH],

Copilot uses AI. Check for mistakes.
})
if isinstance(pa_device, CurvedArrayDetectionGeometry):
k_wave_settings[Tags.SENSOR_RADIUS_MM] = pa_device.radius_mm
Expand All @@ -225,39 +226,152 @@ def k_wave_acoustic_forward_model(self, detection_geometry: DetectionGeometryBas
else:
self.logger.warning(f"Did not find parameter {parameter} in any settings.")

data_dict["settings"] = k_wave_settings
sio.savemat(optical_path, data_dict, long_field_names=True)
raw_time_series_data, Nt, dt = self.run_k_wave_simulation(data_dict, k_wave_settings)

del data_dict, k_wave_settings, detector_positions_mm, pa_device
gc.collect()
self.global_settings[Tags.K_WAVE_SPECIFIC_DT] = dt
self.global_settings[Tags.K_WAVE_SPECIFIC_NT] = Nt

if not simulate_2d:
self.logger.info("Simulating 3D....")
simulation_script_path = "simulate_3D"
else:
self.logger.info("Simulating 2D....")
simulation_script_path = "simulate_2D"
return raw_time_series_data, self.global_settings

def run_k_wave_simulation(self, data: dict, k_wave_settings: Settings):
"""
Run a k-Wave acoustic simulation for the provided medium and sensor configuration.

matlab_binary_path = self.component_settings[Tags.ACOUSTIC_MODEL_BINARY_PATH]
cmd = generate_matlab_cmd(matlab_binary_path, simulation_script_path, optical_path, self.get_additional_flags())
:param data: Dictionary containing simulation input fields. Expected keys include:
- Tags.DATA_FIELD_INITIAL_PRESSURE (np.ndarray)
- Optional: Tags.DATA_FIELD_SPEED_OF_SOUND, Tags.DATA_FIELD_DENSITY,
Tags.DATA_FIELD_ALPHA_COEFF (np.ndarray or scalar)
- Tags.SENSOR_ELEMENT_POSITIONS (np.ndarray)
- Tags.KWAVE_PROPERTY_DIRECTIVITY_ANGLE (np.ndarray)
- 3D only: Tags.KWAVE_PROPERTY_INTRINSIC_EULER_ANGLE (np.ndarray)

cur_dir = os.getcwd()
self.logger.info(cmd)
subprocess.run(cmd)
:param k_wave_settings: k-Wave configuration and simulation parameters (spacing, sensor
settings, PML, GPU flag, etc.).

raw_time_series_data = sio.loadmat(optical_path)[Tags.DATA_FIELD_TIME_SERIES_DATA]
time_grid = sio.loadmat(optical_path + "dt.mat")
num_time_steps = int(np.round(time_grid["number_time_steps"]))
:return: Combined time series data (one trace per physical array element).
:return: Number of simulated time steps.
:return: Time step size [s].

self.global_settings[Tags.K_WAVE_SPECIFIC_DT] = float(time_grid["time_step"])
self.global_settings[Tags.K_WAVE_SPECIFIC_NT] = num_time_steps
"""
GEL_LAYER_HEIGHT = 3

initial_pressure = data[Tags.DATA_FIELD_INITIAL_PRESSURE]
ndim = initial_pressure.ndim
assert ndim in (2, 3), "Only 2D and 3D simulations are supported."

sound_speed = data.get(Tags.DATA_FIELD_SPEED_OF_SOUND, 1540.)
alpha_coeff = data.get(Tags.DATA_FIELD_ALPHA_COEFF, 0.01)
density = data.get(Tags.DATA_FIELD_DENSITY, 1000.)
if ndim == 2:
initial_pressure = np.pad(initial_pressure, ((GEL_LAYER_HEIGHT, 0), (0, 0)))
sound_speed = np.pad(sound_speed, ((GEL_LAYER_HEIGHT, 0), (0, 0)), 'edge')
alpha_coeff = np.pad(alpha_coeff, ((GEL_LAYER_HEIGHT, 0), (0, 0)), 'edge')
density = np.pad(density, ((GEL_LAYER_HEIGHT, 0), (0, 0)), 'edge')
Comment on lines +267 to +269
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The padding operations assume that sound_speed, alpha_coeff, and density are arrays. However, they could be scalars (as indicated by the default values at lines 262-264). If these are scalars, the np.pad operation will fail. Add type checking or ensure these values are converted to arrays before padding.

Copilot uses AI. Check for mistakes.

# Setup source
source = kSource()
source.p0 = initial_pressure

# Setup grid
N = np.asarray(initial_pressure.shape, dtype=int)
spacing_mm = k_wave_settings[Tags.SPACING_MM] / 1000.0
d = spacing_mm * np.ones(ndim)
kgrid = kWaveGrid(N, d)

# Setup medium properties
medium = kWaveMedium(
sound_speed=sound_speed,
alpha_coeff=alpha_coeff,
alpha_power=float(k_wave_settings[Tags.KWAVE_PROPERTY_ALPHA_POWER]),
alpha_mode='no_dispersion',
density=density,
)

# Setup simulation parameters
dt = 1.0 / (k_wave_settings[Tags.SENSOR_SAMPLING_RATE_MHZ] * 1e6)

# Simulate as many time steps as a wave takes to traverse diagonally through the entire tissue
mean_sos = medium.sound_speed.mean()
Nt = int(np.round((np.linalg.norm(N) * spacing_mm / mean_sos) / dt))

# smaller time steps are better for numerical stability in time progressing simulations
# A maximum CFL of 0.3 is advised in the kwave handbook.
# In case we specify something larger, we use a higher sampling rate than anticipated.
# Otherwise we simulate with the target sampling rate
cfl_est = dt / spacing_mm * mean_sos
if cfl_est < 0.3:
kgrid.setTime(Nt, dt)
else:
kgrid.makeTime(sound_speed, cfl=0.3)

os.remove(optical_path)
os.remove(optical_path + "dt.mat")
os.chdir(cur_dir)
# Setup sensor
element_width = float(k_wave_settings[Tags.DETECTOR_ELEMENT_WIDTH_MM]) / 1000.0
angles = data[Tags.KWAVE_PROPERTY_DIRECTIVITY_ANGLE]

return raw_time_series_data, self.global_settings
karray = kWaveArray()
if ndim == 2:
elem_pos = data[Tags.SENSOR_ELEMENT_POSITIONS] / 1000
elem_pos[0] -= 0.5 * kgrid.x_size - spacing_mm * GEL_LAYER_HEIGHT
elem_pos[1] -= 0.5 * kgrid.y_size

elem_pos += 0.5 * element_width * np.vstack((np.sin(angles), np.cos(angles)))
Comment on lines +313 to +317
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In 2D mode, angles is a 1D array in radians (line 183), and it's used directly with np.sin and np.cos without degree conversion. This is correct if the angles are already in radians. However, this differs from 3D mode where angles are in degrees (line 188-190) and require np.deg2rad conversion. Verify that the 2D angles are indeed in radians as expected by this calculation.

Suggested change
elem_pos = data[Tags.SENSOR_ELEMENT_POSITIONS] / 1000
elem_pos[0] -= 0.5 * kgrid.x_size - spacing_mm * GEL_LAYER_HEIGHT
elem_pos[1] -= 0.5 * kgrid.y_size
elem_pos += 0.5 * element_width * np.vstack((np.sin(angles), np.cos(angles)))
# angles are given in degrees; convert to radians for trigonometric calculations
angles_rad = np.deg2rad(angles)
elem_pos = data[Tags.SENSOR_ELEMENT_POSITIONS] / 1000
elem_pos[0] -= 0.5 * kgrid.x_size - spacing_mm * GEL_LAYER_HEIGHT
elem_pos[1] -= 0.5 * kgrid.y_size
elem_pos += 0.5 * element_width * np.vstack((np.sin(angles_rad), np.cos(angles_rad)))

Copilot uses AI. Check for mistakes.
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Inconsistent element position adjustment between 2D and 3D modes. In 2D, the position is adjusted by adding 0.5 * element_width (line 317), while in 3D it's adjusted by subtracting 0.5 * element_width (line 334). This inconsistency may lead to different sensor positioning behavior between 2D and 3D simulations.

Copilot uses AI. Check for mistakes.

for i in range(elem_pos.shape[1]):
karray.add_rect_element(
tuple(elem_pos[:, i]),
element_width, 0.00001,
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The third parameter to add_rect_element (0.00001 in 2D vs 0.0001 in 3D at line 339) differs between 2D and 3D modes. This appears to be the element height/thickness. Consider using a named constant or documenting why these values differ by an order of magnitude.

Copilot uses AI. Check for mistakes.
angles[i]
)

else:
euler_angles = data[Tags.KWAVE_PROPERTY_INTRINSIC_EULER_ANGLE]

elem_pos = data[Tags.SENSOR_ELEMENT_POSITIONS] / 1000
elem_pos[0] -= 0.5 * kgrid.x_size - spacing_mm * GEL_LAYER_HEIGHT
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In 3D mode, the GEL_LAYER_HEIGHT adjustment is applied to elem_pos[0] but the gel layer is not actually added to the initial pressure, sound_speed, alpha_coeff, or density arrays (unlike in 2D mode where padding is applied at lines 266-269). This creates an inconsistency where sensor positions are adjusted for a gel layer that doesn't exist in the simulation domain, which matches the bug mentioned in the PR description (issue #252).

Suggested change
elem_pos[0] -= 0.5 * kgrid.x_size - spacing_mm * GEL_LAYER_HEIGHT
elem_pos[0] -= 0.5 * kgrid.x_size

Copilot uses AI. Check for mistakes.
elem_pos[1] -= 0.5 * kgrid.y_size
elem_pos[2] -= 0.5 * kgrid.z_size

elem_pos -= 0.5 * (element_width * np.sin(np.deg2rad(angles)))
Copy link

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The angles array shape and units are different between 2D and 3D modes. In 2D (line 183), angles is a 1D array in radians. In 3D (line 199-200), angles is a 3xN array in degrees. At line 334, angles is used as if it's a 1D array, but it's actually a 3xN array. This will cause incorrect broadcasting in the element position calculation. The code should likely use only one row of the angles array (e.g., angles[0]) or restructure how angles are stored for 3D.

Suggested change
elem_pos -= 0.5 * (element_width * np.sin(np.deg2rad(angles)))
elem_pos -= 0.5 * (element_width * np.sin(np.deg2rad(angles[0])))

Copilot uses AI. Check for mistakes.

for i in range(elem_pos.shape[1]):
karray.add_rect_element(
tuple(elem_pos[:, i]),
element_width, 0.0001,
tuple(euler_angles[i])
)

# assign binary mask from karray to the sensor mask
sensor = kSensor()
sensor.mask = karray.get_array_binary_mask(kgrid)

if k_wave_settings.get(Tags.MODEL_SENSOR_FREQUENCY_RESPONSE, False):
center_freq = float(k_wave_settings[Tags.SENSOR_CENTER_FREQUENCY_HZ]) # [Hz]
bandwidth = float(k_wave_settings[Tags.SENSOR_BANDWIDTH_PERCENT]) # [%]
sensor.frequency_response([center_freq, bandwidth])

smooth = k_wave_settings.get(Tags.KWAVE_PROPERTY_INITIAL_PRESSURE_SMOOTHING, True)
simulation_options = SimulationOptions(
data_cast='single',
pml_inside=k_wave_settings[Tags.KWAVE_PROPERTY_PMLInside],
pml_alpha=k_wave_settings[Tags.KWAVE_PROPERTY_PMLAlpha],
smooth_p0=smooth,
smooth_c0=smooth,
smooth_rho0=smooth,
pml_auto=True,
save_to_disk=True
)

execution_options = SimulationExecutionOptions(
is_gpu_simulation=k_wave_settings[Tags.GPU]
)

kspaceFirstOrdernD = kspaceFirstOrder3D if ndim == 3 else kspaceFirstOrder2D
sensor_data = kspaceFirstOrdernD(kgrid, source, sensor, medium, simulation_options, execution_options)

# combine data to give one trace per physical array element
time_series_data = karray.combine_sensor_data(kgrid, sensor_data['p'].T)

return time_series_data, int(sensor_data['Nt']), float(sensor_data['dt'])

def perform_k_wave_acoustic_forward_simulation(initial_pressure: np.array,
detection_geometry: DetectionGeometryBase,
Expand Down
Loading
Loading