Switch to k-Wave-python #435

faberno · 2025-12-16T13:36:39Z

This PR replaces all previous functionality to call k-Wave in matlab by its python counterpart k-Wave-python.
Additionally we could unify the 2d and 3d functionality (#397).

Please check the following before creating the pull request (PR):

Did you run automatic tests?
Did you run manual tests?
Is the code provided in the PR still backwards compatible to previous SIMPA versions?

List any specific code review questions
@jgroehl @kdreher Right now the changes should reproduce the old behaviour exactly. But we should think about fixing some bugs/differences between 2D and 3D when we are already refactoring. For example:

2D uses a GEL_LAYER of 3 pixels, 3D does not
3D doesn't pad the GEL_LAYER, but uses it anyways to calculate its element positions (removed gel_layer_height in simulate_3D.m #252)

Other things to clarify:

we could refactor Tags.ACOUSTIC_MODEL_BINARY_PATH from the matlab path to the k-wave binaries path. The k-wave binaries that are currently shipped are not yet compiled for 50xx Nvidia cards ([BUG] Can not use kspaceFirstOrder-CUDA waltsims/k-wave-python#622). So it could be useful to be able to link to your self compiled binaries.
k-wave-python does not support record_movie. If this was important we could readd support for the normal k-wave. Otherwise we can delete all the Tags relating to it (RECORDMOVIE, MOVIENAME, ACOUSTIC_LOG_SCALE, ...)

List any special testing requirements
To test you can for example run three_vs_two_dimensional_simulation_example.py

Provide issue / feature request fixed by this PR

Fixes #402

jgroehl · 2025-12-18T14:38:16Z

Fantastic, thanks for the pull request!
Will take a look in the second week of January and provide feedback.

Copilot

Pull request overview

This PR replaces the MATLAB-based k-Wave acoustic simulation implementation with the Python k-Wave library (k-wave-python v0.4.0), unifying the previously separate 2D and 3D simulation code paths. The migration eliminates the MATLAB dependency and associated subprocess calls, replacing them with direct Python API calls to k-Wave.

Key changes:

Replaced MATLAB script execution with native Python k-Wave API calls in both acoustic forward modeling and time reversal reconstruction
Unified 2D and 3D simulation logic into single methods that branch based on dimensionality
Removed MATLAB utility functions and associated test cases

Reviewed changes

Copilot reviewed 10 out of 10 changed files in this pull request and generated 10 comments.

Show a summary per file

File	Description
simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py	Replaced MATLAB subprocess execution with direct k-Wave Python API calls; unified 2D/3D simulation in run_k_wave_simulation method
simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py	Replaced MATLAB-based time reversal with k-Wave Python implementation in run_k_wave_timereversal method
simpa/core/simulation_modules/acoustic_module/simulate_2D.m	Removed obsolete MATLAB script for 2D acoustic simulation
simpa/core/simulation_modules/acoustic_module/simulate_3D.m	Removed obsolete MATLAB script for 3D acoustic simulation
simpa/core/simulation_modules/reconstruction_module/time_reversal_2D.m	Removed obsolete MATLAB script for 2D time reversal reconstruction
simpa/core/simulation_modules/reconstruction_module/time_reversal_3D.m	Removed obsolete MATLAB script for 3D time reversal reconstruction
simpa/utils/matlab.py	Removed entire MATLAB utility module including command generation function
simpa_tests/automatic_tests/test_additional_flags.py	Removed tests for KWaveAdapter and TimeReversalAdapter MATLAB command generation
simpa/log/file_logger.py	Added propagate=False to prevent log message propagation to root logger
pyproject.toml	Added k-wave-python v0.4.0 dependency

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2025-12-18T14:48:08Z

simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py

+        ----------
+        :param data: Dictionary containing simulation input fields. Expected keys include:
+            - Tags.DATA_FIELD_TIME_SERIES_DATA (np.ndarray)
+            - Tags.Tags.KWAVE_PROPERTY_SENSOR_MASK (np.ndarray)


Duplicate "Tags." prefix in the documentation. Should be "Tags.KWAVE_PROPERTY_SENSOR_MASK" instead of "Tags.Tags.KWAVE_PROPERTY_SENSOR_MASK".

Suggested change

- Tags.Tags.KWAVE_PROPERTY_SENSOR_MASK (np.ndarray)

- Tags.KWAVE_PROPERTY_SENSOR_MASK (np.ndarray)

Copilot · 2025-12-18T14:48:09Z

simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py

+
+        # Setup source
+        source = kSource()
+        # source.p0 = 0


Commented-out code should be removed. If this line (source.p0 = 0) is not needed, it should be deleted rather than left as a comment.

Suggested change

# source.p0 = 0

Copilot · 2025-12-18T14:48:09Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

+            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)))


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 · 2025-12-18T14:48:09Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

+            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)))


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 · 2025-12-18T14:48:10Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

+            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


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 · 2025-12-18T14:48:10Z

pyproject.toml

    "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)


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 · 2025-12-18T14:48:11Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

+            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)))


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 · 2025-12-18T14:48:11Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

            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],


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 · 2025-12-18T14:48:11Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

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


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 · 2025-12-18T14:48:12Z

simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py

+            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')


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.

faberno added 2 commits December 16, 2025 13:54

switch from k-Wave (matlab) to k-Wave-python

d7fba58

fix double logs and docstrings

00614d8

jgroehl requested a review from Copilot December 18, 2025 14:38

Copilot started reviewing on behalf of jgroehl December 18, 2025 14:38 View session

Copilot AI reviewed Dec 18, 2025

View reviewed changes

	- Tags.Tags.KWAVE_PROPERTY_SENSOR_MASK (np.ndarray)
	- Tags.KWAVE_PROPERTY_SENSOR_MASK (np.ndarray)

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

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

	"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)

Switch to k-Wave-python #435

Are you sure you want to change the base?

Switch to k-Wave-python #435

Uh oh!

Conversation

faberno commented Dec 16, 2025

Uh oh!

jgroehl commented Dec 18, 2025

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 18, 2025

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants