From d7fba581b79e36a967e59e4cee5bd0364f10d8b6 Mon Sep 17 00:00:00 2001 From: Fabian Schneider Date: Tue, 16 Dec 2025 13:54:54 +0100 Subject: [PATCH 1/2] switch from k-Wave (matlab) to k-Wave-python --- pyproject.toml | 1 + .../acoustic_module/k_wave_adapter.py | 182 ++++++++++++++---- .../acoustic_module/simulate_2D.m | 172 ----------------- .../acoustic_module/simulate_3D.m | 173 ----------------- .../reconstruction_module/time_reversal_2D.m | 124 ------------ .../reconstruction_module/time_reversal_3D.m | 118 ------------ .../time_reversal_adapter.py | 124 ++++++++---- simpa/utils/matlab.py | 39 ---- .../automatic_tests/test_additional_flags.py | 22 --- 9 files changed, 240 insertions(+), 715 deletions(-) delete mode 100644 simpa/core/simulation_modules/acoustic_module/simulate_2D.m delete mode 100644 simpa/core/simulation_modules/acoustic_module/simulate_3D.m delete mode 100644 simpa/core/simulation_modules/reconstruction_module/time_reversal_2D.m delete mode 100644 simpa/core/simulation_modules/reconstruction_module/time_reversal_3D.m delete mode 100644 simpa/utils/matlab.py diff --git a/pyproject.toml b/pyproject.toml index 08fe5571..acf5fcf7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -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) ] [project.optional-dependencies] diff --git a/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py b/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py index 2348ec37..839232d4 100644 --- a/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py +++ b/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py @@ -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, @@ -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 @@ -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, @@ -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], }) if isinstance(pa_device, CurvedArrayDetectionGeometry): k_wave_settings[Tags.SENSOR_RADIUS_MM] = pa_device.radius_mm @@ -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: 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') + + # 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))) + + for i in range(elem_pos.shape[1]): + karray.add_rect_element( + tuple(elem_pos[:, i]), + element_width, 0.00001, + 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 + 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))) + + 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, diff --git a/simpa/core/simulation_modules/acoustic_module/simulate_2D.m b/simpa/core/simulation_modules/acoustic_module/simulate_2D.m deleted file mode 100644 index c24e6bde..00000000 --- a/simpa/core/simulation_modules/acoustic_module/simulate_2D.m +++ /dev/null @@ -1,172 +0,0 @@ -%%SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ -%%SPDX-FileCopyrightText: 2021 Janek Groehl -%%SPDX-License-Identifier: MIT - -function [] = simulate_2D(optical_path) - -%% In case of an error, make sure the matlab scripts exits anyway -clean_up = onCleanup(@exit); - -%% Read settings file - -data = load(optical_path); -settings = data.settings; - -%% Read initial pressure - -source.p0 = data.initial_pressure; - -% Choose if the initial pressure should be smoothed before simulation -if isfield(settings, 'initial_pressure_smoothing') == true - p0_smoothing = settings.initial_pressure_smoothing; -else - p0_smoothing = true; -end - -%% Define kWaveGrid - -% add N pixel "gel" to reduce Fourier artifact -GEL_LAYER_HEIGHT = 3; - -source.p0 = padarray(source.p0, [GEL_LAYER_HEIGHT 0], 0, 'pre'); -[Nx, Ny] = size(source.p0); -disp(Nx); -disp(Ny); -if isfield(settings, 'sample') == true - if settings.sample == true - dx = double(settings.voxel_spacing_mm) / (double(settings.upscale_factor) * 1000); - else - dx = double(settings.voxel_spacing_mm) / 1000; % convert from mm to m - end -else - dx = double(settings.voxel_spacing_mm) / 1000; % convert from mm to m -end -kgrid = kWaveGrid(Nx, dx, Ny, dx); - -%% Define medium - -% if a field of the struct "data" is given which describes the sound speed, the array is loaded and is used as medium.sound_speed -if isfield(data, 'sos') == true - medium.sound_speed = data.sos; - medium.sound_speed = padarray(medium.sound_speed, [GEL_LAYER_HEIGHT 0], 'replicate', 'pre'); -else - medium.sound_speed = 1540; -end - -% if a field of the struct "data" is given which describes the attenuation, the array is loaded and is used as medium.alpha_coeff -if isfield(data, 'alpha_coeff') == true - medium.alpha_coeff = data.alpha_coeff; - medium.alpha_coeff = padarray(medium.alpha_coeff, [GEL_LAYER_HEIGHT 0], 'replicate', 'pre'); -else - medium.alpha_coeff = 0.01; -end - -medium.alpha_power = double(settings.medium_alpha_power); % b for a * MHz ^ b -medium.alpha_mode = 'no_dispersion'; - -% if a field of the struct "data" is given which describes the density, the array is loaded and is used as medium.density -if isfield(data, 'density') == true - medium.density = data.density; - medium.density = padarray(medium.density, [GEL_LAYER_HEIGHT 0], 'replicate', 'pre'); -else - medium.density = 1000 * ones(Nx, Ny); -end - -%% Sampling rate - -% load sampling rate from settings -dt = 1.0 / double(settings.sensor_sampling_rate_mhz * 1000000); - -% Simulate as many time steps as a wave takes to traverse diagonally through the entire tissue -Nt = round((sqrt(Ny*Ny+Nx*Nx)*dx / mean(medium.sound_speed, 'all')) / dt); - -estimated_cfl_number = dt / dx * mean(medium.sound_speed, 'all'); -disp(estimated_cfl_number); - -% 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 -if estimated_cfl_number < 0.3 - kgrid.setTime(Nt, dt); - disp("Setting custom time!"); -else - kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3); - disp("Making time!"); -end - -%% Define sensor - -% create empty array -karray = kWaveArray; - -elem_pos = data.sensor_element_positions/1000; - - -elem_pos(1, :) = elem_pos(1, :) - 0.5 * kgrid.x_size + dx * GEL_LAYER_HEIGHT; -elem_pos(2, :) = elem_pos(2, :) - 0.5 * kgrid.y_size; - -num_elements = size(elem_pos, 2); - -element_width = double(settings.detector_element_width_mm)/1000; -angles = data.directivity_angle; - -if isfield(settings, 'sensor_radius_mm') == true - radius_of_curv = double(settings.sensor_radius_mm)/1000; -end - - -% add elements to the array - -for ind = 1:num_elements - x = elem_pos(1, ind); - y = elem_pos(2, ind); - alpha = angles(1, ind); - x = x + 0.5 * (element_width*sin(alpha)); - y = y + 0.5 * (element_width*cos(alpha)); - karray.addRectElement([x, y], element_width, 0.00001, [angles(1, ind)]); -end - -% assign binary mask from karray to the sensor mask -sensor.mask = karray.getArrayBinaryMask(kgrid); - -% model sensor frequency response -if isfield(settings, 'model_sensor_frequency_response') == true - if settings.model_sensor_frequency_response == true - center_freq = double(settings.sensor_center_frequency); % [Hz] - bandwidth = double(settings.sensor_bandwidth); % [%] - sensor.frequency_response = [center_freq, bandwidth]; - end -end - -%% Computation settings - -if settings.gpu == true - datacast = 'gpuArray-single'; -else - datacast = 'single'; -end - -input_args = {'DataCast', datacast, 'PMLInside', settings.pml_inside, ... - 'PMLAlpha', settings.pml_alpha, 'PMLSize', 'auto', ... - 'PlotPML', settings.plot_pml, 'RecordMovie', settings.record_movie, ... - 'MovieName', settings.movie_name, 'PlotScale', [-1, 1], 'LogScale', settings.acoustic_log_scale, ... - 'Smooth', p0_smoothing}; - -if settings.gpu == true - time_series_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}); - time_series_data = gather(time_series_data); -else - time_series_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); -end - -% combine data to give one trace per physical array element -time_series_data = karray.combineSensorData(kgrid, time_series_data); - -%% Write data to mat array -save(optical_path, 'time_series_data')%, '-v7.3') -time_step = kgrid.dt -number_time_steps = kgrid.Nt -save(strcat(optical_path, 'dt.mat'), 'time_step', 'number_time_steps'); - -end diff --git a/simpa/core/simulation_modules/acoustic_module/simulate_3D.m b/simpa/core/simulation_modules/acoustic_module/simulate_3D.m deleted file mode 100644 index 1772c0c8..00000000 --- a/simpa/core/simulation_modules/acoustic_module/simulate_3D.m +++ /dev/null @@ -1,173 +0,0 @@ -%%SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ -%%SPDX-FileCopyrightText: 2021 Janek Groehl -%%SPDX-License-Identifier: MIT - -function [] = simulate_3D(optical_path) - -%% In case of an error, make sure the matlab scripts exits anyway -clean_up = onCleanup(@exit); - -%% Read settings file - -data = load(optical_path); -settings = data.settings; - -%% Read initial pressure - -source.p0 = data.initial_pressure; - -% Choose if the initial pressure should be smoothed before simulation -if isfield(settings, 'initial_pressure_smoothing') == true - p0_smoothing = settings.initial_pressure_smoothing; -else - p0_smoothing = true; -end - -%% Define kWaveGrid - -% add 2 pixel "gel" to reduce Fourier artifact -GEL_LAYER_HEIGHT = 3; - -%source.p0 = padarray(source.p0, [GEL_LAYER_HEIGHT 0], 0, 'pre'); -[Nx, Ny, Nz] = size(source.p0); -if isfield(settings, 'sample') == true - if settings.sample == true - dx = double(settings.voxel_spacing_mm)/(double(settings.upscale_factor) * 1000); - else - dx = double(settings.voxel_spacing_mm)/1000; % convert from mm to m - end -else - dx = double(settings.voxel_spacing_mm)/1000; % convert from mm to m -end -kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx); - -%% Define medium - -% if a field of the struct "data" is given which describes the sound speed, the array is loaded and is used as medium.sound_speed -if isfield(data, 'sos') == true - medium.sound_speed = data.sos; - % add 2 pixel "gel" to reduce Fourier artifact -% medium.sound_speed = padarray(medium.sound_speed, [GEL_LAYER_HEIGHT 0], 'replicate', 'pre'); -else - medium.sound_speed = 1540; -end - -% if a field of the struct "data" is given which describes the attenuation, the array is loaded and is used as medium.alpha_coeff -if isfield(data, 'alpha_coeff') == true - medium.alpha_coeff = data.alpha_coeff; - % add 2 pixel "gel" to reduce Fourier artifact -% medium.alpha_coeff = padarray(medium.alpha_coeff, [GEL_LAYER_HEIGHT 0], 'replicate', 'pre'); -else - medium.alpha_coeff = 0.01; -end - -medium.alpha_power = double(settings.medium_alpha_power); % b for a * MHz ^ b -medium.alpha_mode = 'no_dispersion'; - -% if a field of the struct "data" is given which describes the density, the array is loaded and is used as medium.density -if isfield(data, 'density') == true - medium.density = data.density; - % add 2 pixel "gel" to reduce Fourier artifact -% medium.density = padarray(medium.density, [GEL_LAYER_HEIGHT 0], 'replicate', 'pre'); -else - medium.density = 1000*ones(Nx, Ny, Nz); -end - -%% Sampling rate - -% load sampling rate from settings -dt = 1.0 / double(settings.sensor_sampling_rate_mhz * 1000000); - -% Simulate as many time steps as a wave takes to traverse diagonally through the entire tissue -Nt = round((sqrt(Ny*Ny+Nx*Nx+Nz*Nz)*dx / mean(medium.sound_speed, 'all')) / dt); - -estimated_cfl_number = dt / dx * mean(medium.sound_speed, 'all'); - -% smaller time steps are better for numerical stability in time progressing simulations -% A minimum 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 -if estimated_cfl_number < 0.3 - kgrid.setTime(Nt, dt); -else - kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3); -end - -%% Define sensor - -% create empty array -karray = kWaveArray; - -elem_pos = data.sensor_element_positions/1000; - -elem_pos(1, :) = elem_pos(1, :) - 0.5 * kgrid.x_size + dx * GEL_LAYER_HEIGHT; -elem_pos(2, :) = elem_pos(2, :) - 0.5 * kgrid.y_size; -elem_pos(3, :) = elem_pos(3, :) - 0.5 * kgrid.z_size; -num_elements = size(elem_pos, 2); - -element_width = double(settings.detector_element_width_mm)/1000; -orientation_angles = data.directivity_angle; -euler_angles = data.intrinsic_euler_angle; - -if isfield(settings, 'sensor_radius_mm') == true - radius_of_curv = double(settings.sensor_radius_mm)/1000; -end - -% For addArcElement orient all elements towards the focus -% For the iThera MSOT Acuity Echo, it is [0.008, 0] - -%focus_pos = [0.008, 0]; - -% add elements to the array - -%for ind = 1:num_elements -% karray.addArcElement(elem_pos(:, ind), radius_of_curv, element_width, focus_pos); -%end -for ind = 1:num_elements - elem_pos(:, ind) = elem_pos(:, ind) - 0.5*(element_width*sind(orientation_angles(:, ind))); - karray.addRectElement(elem_pos(:, ind), element_width, 0.0001, euler_angles(ind, :)); -end - -% assign binary mask from karray to the sensor mask -sensor.mask = karray.getArrayBinaryMask(kgrid); - -% model sensor frequency response -if isfield(settings, 'model_sensor_frequency_response') == true - if settings.model_sensor_frequency_response == true - center_freq = double(settings.sensor_center_frequency); % [Hz] - bandwidth = double(settings.sensor_bandwidth); % [%] - sensor.frequency_response = [center_freq, bandwidth]; - end -end - -%% Computation settings - -if settings.gpu == true - datacast = 'gpuArray-single'; -else - datacast = 'single'; -end - -input_args = {'DataCast', datacast, 'PMLInside', settings.pml_inside, ... - 'PMLAlpha', double(settings.pml_alpha), 'PMLSize', 'auto', ... - 'PlotPML', settings.plot_pml, 'RecordMovie', settings.record_movie, ... - 'MovieName', settings.movie_name, 'PlotScale', [-1, 1], 'LogScale', settings.acoustic_log_scale, ... - 'Smooth', p0_smoothing}; - -if settings.gpu == true - time_series_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:}); - time_series_data = gather(time_series_data); -else - time_series_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); -end - -% combine data to give one trace per physical array element -time_series_data = karray.combineSensorData(kgrid, time_series_data); - -%% Write data to mat array -save(optical_path, 'time_series_data')%, '-v7.3') -time_step = kgrid.dt; -number_time_steps = kgrid.Nt; -save(strcat(optical_path, 'dt.mat'), 'time_step', 'number_time_steps'); - -end \ No newline at end of file diff --git a/simpa/core/simulation_modules/reconstruction_module/time_reversal_2D.m b/simpa/core/simulation_modules/reconstruction_module/time_reversal_2D.m deleted file mode 100644 index 990963b0..00000000 --- a/simpa/core/simulation_modules/reconstruction_module/time_reversal_2D.m +++ /dev/null @@ -1,124 +0,0 @@ -%%SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ -%%SPDX-FileCopyrightText: 2021 Janek Groehl -%%SPDX-License-Identifier: MIT - -function [] = time_reversal_2D(acoustic_path) - -%% In case of an error, make sure the matlab scripts exits anyway -clean_up = onCleanup(@exit); - -%% Read settings file -data = load(acoustic_path); - -settings = data.settings; - -%% Read time_series_data -time_series_data = data.time_series_data; - -%% Define kWaveGrid - -[Nx, Ny] = size(data.sensor_mask); -if isfield(settings, 'sample') == true - if settings.sample == true - dx = double(settings.voxel_spacing_mm)/(double(settings.upscale_factor) * 1000); - else - dx = double(settings.voxel_spacing_mm)/1000; % convert from mm to m - end -else - dx = double(settings.voxel_spacing_mm)/1000; % convert from mm to m -end -kgrid = kWaveGrid(Nx, dx, Ny, dx); -source.p0 = 0; - -%% Define medium - -% if a field of the struct "data" is given which describes the sound speed, the array is loaded and is used as medium.sound_speed -if isfield(data, 'sos') == true - medium.sound_speed = double(data.sos); -else - medium.sound_speed = 1540; -end - -% if a field of the struct "data" is given which describes the attenuation, the array is loaded and is used as medium.alpha_coeff -if isfield(data, 'alpha_coeff') == true - medium.alpha_coeff = double(data.alpha_coeff); -else - medium.alpha_coeff = 0.01; -end - -medium.alpha_power = double(settings.medium_alpha_power); % b for a * MHz ^ b - -% if a field of the struct "data" is given which describes the density, the array is loaded and is used as medium.density -if isfield(data, 'density') == true - medium.density = double(data.density); -else - medium.density = 1000*ones(Nx, Ny); -end - -kgrid.setTime(settings.Nt, settings.dt) - -%sound_speed_ref = min(min(medium.sound_speed)); -%kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3); % time array with -% CFL number of 0.3 (advised by manual) -% Using makeTime, dt = CFL*dx/medium.sound_speed and the total -% time is set to the time it would take for an acoustic wave to travel -% across the longest grid diagonal. - -%% Define sensor - -sensor.mask = data.sensor_mask; - - -% if a field of the struct "data" is given which describes the sensor directivity angles, the array is loaded and is used as sensor.directivity_angle -%if isfield(data, 'directivity_angle') == true -% sensor.directivity_angle = data.directivity_angle; -%end -% -%if isfield(data, 'directivity_size') -% sensor.directivity_size = settings.sensor_directivity_size; -%end - -%sensor.directivity_pattern = settings.sensor_directivity_pattern; - -% model sensor frequency response -if isfield(settings, 'model_sensor_frequency_response') == true - if settings.model_sensor_frequency_response == true - center_freq = double(settings.sensor_center_frequency); % [Hz] - bandwidth = double(settings.sensor_bandwidth); % [%] - sensor.frequency_response = [center_freq, bandwidth]; - end -end - -% kwave expects the time series data to be in an order -% "based on the angle that each sensor point makes with the centre of the grid". -[simpa, sort_idx] = reorderSensorData(kgrid, sensor, time_series_data); -len_array = sort(sort_idx); -[simpa, real_idx] = intersect(sort_idx, len_array); -time_series_data = time_series_data(real_idx, :); - -sensor.time_reversal_boundary_data = time_series_data; - -%% Computation settings - -if settings.gpu == true - datacast = 'gpuArray-single'; -else - datacast = 'single'; -end - -input_args = {'DataCast', datacast, 'PMLInside', settings.pml_inside, ... - 'PMLAlpha', settings.pml_alpha, 'PMLSize', 'auto', ... - 'PlotPML', settings.plot_pml, 'RecordMovie', settings.record_movie, ... - 'MovieName', settings.movie_name, 'PlotScale', [0, 1]}; - -if settings.gpu == true - reconstructed_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}); - reconstructed_data = gather(reconstructed_data); -else - reconstructed_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); -end - -%% Write data to mat array -save(strcat(acoustic_path, 'tr.mat'), 'reconstructed_data') - -end \ No newline at end of file diff --git a/simpa/core/simulation_modules/reconstruction_module/time_reversal_3D.m b/simpa/core/simulation_modules/reconstruction_module/time_reversal_3D.m deleted file mode 100644 index edb7371c..00000000 --- a/simpa/core/simulation_modules/reconstruction_module/time_reversal_3D.m +++ /dev/null @@ -1,118 +0,0 @@ -%%SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ -%%SPDX-FileCopyrightText: 2021 Janek Groehl -%%SPDX-License-Identifier: MIT - -function [] = time_reversal_3D(acoustic_path) - -%% In case of an error, make sure the matlab scripts exits anyway -clean_up = onCleanup(@exit); - -%% Read settings file -data = load(acoustic_path); - -settings = data.settings; - -%% Read time_series_data -time_series_data = data.time_series_data; - -%% Define kWaveGrid - -[Nx, Ny, Nz] = size(data.sensor_mask); -if isfield(settings, 'sample') == true - if settings.sample == true - dx = double(settings.voxel_spacing_mm)/(double(settings.upscale_factor) * 1000); - else - dx = double(settings.voxel_spacing_mm)/1000; % convert from mm to m - end -else - dx = double(settings.voxel_spacing_mm)/1000; % convert from mm to m -end -kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx); -source.p0 = 0; - -%% Define medium - -% if a field of the struct "data" is given which describes the sound speed, the array is loaded and is used as medium.sound_speed -if isfield(data, 'sos') == true - medium.sound_speed = double(data.sos); -else - medium.sound_speed = 1540; -end - -% if a field of the struct "data" is given which describes the attenuation, the array is loaded and is used as medium.alpha_coeff -if isfield(data, 'alpha_coeff') == true - medium.alpha_coeff = double(data.alpha_coeff); -else - medium.alpha_coeff = 0.01; -end - -medium.alpha_power = double(settings.medium_alpha_power); % b for a * MHz ^ b - -% if a field of the struct "data" is given which describes the density, the array is loaded and is used as medium.density -if isfield(data, 'density') == true - medium.density = double(data.density); -else - medium.density = 1000*ones(Nx, Ny, Nz); -end - -kgrid.setTime(settings.Nt, settings.dt) - -%sound_speed_ref = min(min(medium.sound_speed)); -%kgrid.dt = 1 / (settings.sensor_sampling_rate_mhz * 10^6); -%kgrid.Nt = ceil((sqrt((Nx*dx)^2+(Ny*dx)^2) / sound_speed_ref) / kgrid.dt); -%kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.3); % time array with -% CFL number of 0.3 (advised by manual) -% Using makeTime, dt = CFL*dx/medium.sound_speed and the total -% time is set to the time it would take for an acoustic wave to travel -% across the longest grid diagonal. - -%% Define sensor - -sensor.mask = data.sensor_mask; - -% if a field of the struct "data" is given which describes the sensor directivity angles, the array is loaded and is used as sensor.directivity_angle -%if isfield(data, 'directivity_angle') == true -% sensor.directivity_angle = data.directivity_angle; -%end -% -%if isfield(data, 'directivity_size') -% sensor.directivity_size = settings.sensor_directivity_size; -%end - -% sensor.directivity_pattern = settings.sensor_directivity_pattern; - -% model sensor frequency response -if isfield(settings, 'model_sensor_frequency_response') == true - if settings.model_sensor_frequency_response == true - center_freq = double(settings.sensor_center_frequency); % [Hz] - bandwidth = double(settings.sensor_bandwidth); % [%] - sensor.frequency_response = [center_freq, bandwidth]; - end -end - -sensor.time_reversal_boundary_data = time_series_data; - -%% Computation settings - -if settings.gpu == true - datacast = 'gpuArray-single'; -else - datacast = 'single'; -end - -input_args = {'DataCast', datacast, 'PMLInside', settings.pml_inside, ... - 'PMLAlpha', settings.pml_alpha, 'PMLSize', 'auto', ... - 'PlotPML', settings.plot_pml, 'RecordMovie', settings.record_movie, ... - 'MovieName', settings.movie_name, 'PlotScale', [0, 1]}; - -if settings.gpu == true - reconstructed_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:}); - reconstructed_data = gather(reconstructed_data); -else - reconstructed_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); -end - -%% Write data to mat array -save(strcat(acoustic_path, 'tr.mat'), 'reconstructed_data') - -end \ No newline at end of file diff --git a/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py b/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py index 5c61ff44..d52897ec 100644 --- a/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py +++ b/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py @@ -2,16 +2,22 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +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 +from kwave.utils.signals import reorder_sensor_data + from simpa.core.simulation_modules.reconstruction_module.reconstruction_utils import compute_image_dimensions from simpa.utils import Tags, round_x5_away_from_zero -from simpa.utils.matlab import generate_matlab_cmd from simpa.utils.settings import Settings from simpa.core.simulation_modules.reconstruction_module import ReconstructionAdapterBase from simpa.core.device_digital_twins import LinearArrayDetectionGeometry import numpy as np -import scipy.io as sio -import subprocess -import os class TimeReversalAdapter(ReconstructionAdapterBase): @@ -129,7 +135,6 @@ def reconstruction_algorithm(self, time_series_sensor_data, detection_geometry): input_data[Tags.DATA_FIELD_TIME_SERIES_DATA] = time_series_sensor_data input_data, spacing_in_mm = self.get_acoustic_properties(input_data, detection_geometry) - acoustic_path = self.global_settings[Tags.SIMPA_OUTPUT_FILE_PATH] + ".mat" possible_k_wave_parameters = [Tags.MODEL_SENSOR_FREQUENCY_RESPONSE, Tags.KWAVE_PROPERTY_ALPHA_POWER, Tags.GPU, Tags.KWAVE_PROPERTY_PMLInside, Tags.KWAVE_PROPERTY_PMLAlpha, Tags.KWAVE_PROPERTY_PlotPML, @@ -152,35 +157,15 @@ def reconstruction_algorithm(self, time_series_sensor_data, detection_geometry): k_wave_settings[parameter] = self.global_settings[parameter] if Tags.K_WAVE_SPECIFIC_DT in self.global_settings and Tags.K_WAVE_SPECIFIC_NT in self.global_settings: - k_wave_settings["dt"] = self.global_settings[Tags.K_WAVE_SPECIFIC_DT] - k_wave_settings["Nt"] = self.global_settings[Tags.K_WAVE_SPECIFIC_NT] + k_wave_settings[Tags.K_WAVE_SPECIFIC_DT] = self.global_settings[Tags.K_WAVE_SPECIFIC_DT] + k_wave_settings[Tags.K_WAVE_SPECIFIC_NT] = self.global_settings[Tags.K_WAVE_SPECIFIC_NT] else: num_samples = time_series_sensor_data.shape[1] time_per_sample_s = 1 / (self.component_settings[Tags.SENSOR_SAMPLING_RATE_MHZ] * 1000000) - k_wave_settings["dt"] = time_per_sample_s - k_wave_settings["Nt"] = num_samples - input_data["settings"] = k_wave_settings - sio.savemat(acoustic_path, input_data, long_field_names=True) - - if Tags.ACOUSTIC_SIMULATION_3D in self.component_settings and \ - self.component_settings[Tags.ACOUSTIC_SIMULATION_3D]: - time_reversal_script = "time_reversal_3D" - axes = (0, 2) - else: - time_reversal_script = "time_reversal_2D" - axes = (0, 1) - - matlab_binary_path = self.component_settings[Tags.ACOUSTIC_MODEL_BINARY_PATH] - cmd = generate_matlab_cmd(matlab_binary_path, time_reversal_script, acoustic_path, self.get_additional_flags()) - - cur_dir = os.getcwd() - os.chdir(self.global_settings[Tags.SIMULATION_PATH]) - self.logger.info(cmd) - subprocess.run(cmd) - - reconstructed_data = sio.loadmat(acoustic_path + "tr.mat")[Tags.DATA_FIELD_RECONSTRUCTED_DATA] + k_wave_settings[Tags.K_WAVE_SPECIFIC_DT] = time_per_sample_s + k_wave_settings[Tags.K_WAVE_SPECIFIC_NT] = num_samples - reconstructed_data = reconstructed_data.T + reconstructed_data = self.run_k_wave_timereversal(input_data, k_wave_settings) field_of_view_mm = detection_geometry.get_field_of_view_mm() _, _, _, xdim_start, xdim_end, ydim_start, ydim_end, zdim_start, zdim_end = compute_image_dimensions( @@ -205,8 +190,81 @@ def reconstruction_algorithm(self, time_series_sensor_data, detection_geometry): self.logger.critical("Unexpected number of dimensions in reconstructed image. " f"Expected 2 or 3 but was {len(np.shape(reconstructed_data))}") - os.chdir(cur_dir) - os.remove(acoustic_path) - os.remove(acoustic_path + "tr.mat") + return reconstructed_data + + def run_k_wave_timereversal(self, data: dict, k_wave_settings: Settings): + """ + Run k-Wave time reversal for the provided medium and sensor configuration. + + Parameters + ---------- + :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) + - Optional: Tags.DATA_FIELD_SPEED_OF_SOUND, Tags.DATA_FIELD_DENSITY, + Tags.DATA_FIELD_ALPHA_COEFF (scalar) + + k_wave_settings : Settings + k-Wave configuration and simulation parameters (spacing, sensor settings, PML, GPU flag, etc.). + + :return: Reconstructed image + """ + + sensor_mask = data[Tags.KWAVE_PROPERTY_SENSOR_MASK] + ndim = sensor_mask.ndim + assert ndim in (2, 3), "Only 2D and 3D simulations are supported." + + # Setup grid + N = np.asarray(sensor_mask.shape, dtype=int) + spacing_mm = k_wave_settings[Tags.SPACING_MM] / 1000.0 + d = spacing_mm * np.ones(ndim) + kgrid = kWaveGrid(N, d) + kgrid.setTime(k_wave_settings[Tags.K_WAVE_SPECIFIC_NT], k_wave_settings[Tags.K_WAVE_SPECIFIC_DT]) + + # Setup source + source = kSource() + # source.p0 = 0 + + # Setup medium properties + medium = kWaveMedium( + sound_speed=data.get(Tags.DATA_FIELD_SPEED_OF_SOUND, 1540), + alpha_coeff=data.get(Tags.DATA_FIELD_ALPHA_COEFF, 0.01), + alpha_power=float(k_wave_settings[Tags.KWAVE_PROPERTY_ALPHA_POWER]), + alpha_mode='no_dispersion', + density=data.get(Tags.DATA_FIELD_DENSITY, 1000), + ) + + # Setup sensor + sensor = kSensor() + sensor.mask = sensor_mask + + 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]) + + time_series_data = data[Tags.DATA_FIELD_TIME_SERIES_DATA] + if ndim == 2: + # kwave expects the time series data to be in an order + # "based on the angle that each sensor point makes with the centre of the grid" + time_series_data = reorder_sensor_data(kgrid, sensor, time_series_data) + + sensor.time_reversal_boundary_data = time_series_data + + simulation_options = SimulationOptions( + data_cast='single', + pml_inside=k_wave_settings[Tags.KWAVE_PROPERTY_PMLInside], + pml_alpha=k_wave_settings[Tags.KWAVE_PROPERTY_PMLAlpha], + 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 + reconstructed_data = kspaceFirstOrdernD(kgrid, source, sensor, medium, simulation_options, execution_options) + reconstructed_data = reconstructed_data['p_final'] return reconstructed_data diff --git a/simpa/utils/matlab.py b/simpa/utils/matlab.py deleted file mode 100644 index 1459ba0b..00000000 --- a/simpa/utils/matlab.py +++ /dev/null @@ -1,39 +0,0 @@ -# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ -# SPDX-FileCopyrightText: 2021 Janek Groehl -# SPDX-License-Identifier: MIT - -import inspect -import os -from typing import List - - -def generate_matlab_cmd(matlab_binary_path: str, simulation_script_path: str, data_path: str, additional_flags: List[str] = []) -> List[str]: - """Generates the MATLAB execution command from the given paths - - :param matlab_binary_path: path to the MATLAB binary file as defined by PathManager - :type matlab_binary_path: str - :param simulation_script_path: path to the MATLAB script that should be run (either simulate_2D.m or simulate_3D.m) - :type simulation_script_path: str - :param data_path: path to the .mat file used for simulating - :type data_path: str - :param additional_flags: list of optional additional flags for MATLAB - :type additional_flags: List[str] - :return: list of command parts - :rtype: List[str] - """ - - # get path of calling script to add to matlab path - base_script_path = os.path.dirname(os.path.abspath(inspect.stack()[1].filename)) - # ensure data path is an absolute path - data_path = os.path.abspath(data_path) - - cmd = list() - cmd.append(matlab_binary_path) - cmd.append("-nodisplay") - cmd.append("-nosplash") - cmd.append("-automation") - cmd.append("-wait") - cmd += additional_flags - cmd.append("-r") - cmd.append(f"addpath('{base_script_path}');{simulation_script_path}('{data_path}');exit;") - return cmd diff --git a/simpa_tests/automatic_tests/test_additional_flags.py b/simpa_tests/automatic_tests/test_additional_flags.py index 48968956..81912024 100644 --- a/simpa_tests/automatic_tests/test_additional_flags.py +++ b/simpa_tests/automatic_tests/test_additional_flags.py @@ -5,7 +5,6 @@ import unittest from simpa import MCXReflectanceAdapter, MCXAdapter, KWaveAdapter, TimeReversalAdapter, Tags, Settings -from simpa.utils.matlab import generate_matlab_cmd class TestAdditionalFlags(unittest.TestCase): @@ -35,27 +34,6 @@ def test_get_cmd_mcx_adapter(self): self.assertIn( flag, cmd, f"{flag} was not in command returned by mcx adapter but was defined as additional flag") - def test_get_cmd_kwave_adapter(self): - self.settings.set_acoustic_settings({ - Tags.ADDITIONAL_FLAGS: self.additional_flags - }) - kwave_adapter = KWaveAdapter(global_settings=self.settings) - cmd = generate_matlab_cmd("./matlab.exe", "simulate_2D.m", "my_hdf5.mat", kwave_adapter.get_additional_flags()) - for flag in self.additional_flags: - self.assertIn( - flag, cmd, f"{flag} was not in command returned by kwave adapter but was defined as additional flag") - - def test_get_cmd_time_reversal_adapter(self): - self.settings.set_reconstruction_settings({ - Tags.ADDITIONAL_FLAGS: self.additional_flags - }) - time_reversal_adapter = TimeReversalAdapter(global_settings=self.settings) - cmd = generate_matlab_cmd("./matlab.exe", "time_reversal_2D.m", "my_hdf5.mat", - time_reversal_adapter.get_additional_flags()) - for flag in self.additional_flags: - self.assertIn( - flag, cmd, f"{flag} was not in command returned by time reversal adapter but was defined as additional flag") - if __name__ == '__main__': unittest.main() From 00614d8e2107224f6c361d90710bb2eed47626fd Mon Sep 17 00:00:00 2001 From: Fabian Schneider Date: Tue, 16 Dec 2025 17:28:46 +0100 Subject: [PATCH 2/2] fix double logs and docstrings --- .../simulation_modules/acoustic_module/k_wave_adapter.py | 4 ++-- .../reconstruction_module/time_reversal_adapter.py | 6 +++--- simpa/log/file_logger.py | 1 + 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py b/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py index 839232d4..ecef8961 100644 --- a/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py +++ b/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py @@ -245,8 +245,8 @@ def run_k_wave_simulation(self, data: dict, k_wave_settings: Settings): - Tags.KWAVE_PROPERTY_DIRECTIVITY_ANGLE (np.ndarray) - 3D only: Tags.KWAVE_PROPERTY_INTRINSIC_EULER_ANGLE (np.ndarray) - :param k_wave_settings: Settings - k-Wave configuration and simulation parameters (spacing, sensor settings, PML, GPU flag, etc.). + :param k_wave_settings: k-Wave configuration and simulation parameters (spacing, sensor + settings, PML, GPU flag, etc.). :return: Combined time series data (one trace per physical array element). :return: Number of simulated time steps. diff --git a/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py b/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py index d52897ec..9626bb72 100644 --- a/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py +++ b/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py @@ -204,8 +204,8 @@ def run_k_wave_timereversal(self, data: dict, k_wave_settings: Settings): - Optional: Tags.DATA_FIELD_SPEED_OF_SOUND, Tags.DATA_FIELD_DENSITY, Tags.DATA_FIELD_ALPHA_COEFF (scalar) - k_wave_settings : Settings - k-Wave configuration and simulation parameters (spacing, sensor settings, PML, GPU flag, etc.). + :input k_wave_settings: k-Wave configuration and simulation parameters (spacing, sensor + settings, PML, GPU flag, etc.). :return: Reconstructed image """ @@ -260,7 +260,7 @@ def run_k_wave_timereversal(self, data: dict, k_wave_settings: Settings): ) execution_options = SimulationExecutionOptions( - is_gpu_simulation=k_wave_settings[Tags.GPU], + is_gpu_simulation=k_wave_settings[Tags.GPU] ) kspaceFirstOrdernD = kspaceFirstOrder3D if ndim == 3 else kspaceFirstOrder2D diff --git a/simpa/log/file_logger.py b/simpa/log/file_logger.py index 5911bdbc..ddaab1b7 100644 --- a/simpa/log/file_logger.py +++ b/simpa/log/file_logger.py @@ -50,6 +50,7 @@ def __new__(cls, path=None, force_new_instance=False, startup_verbose=False): cls._logger.addHandler(console_handler) cls._logger.addHandler(file_handler) + cls._logger.propagate = False if startup_verbose: cls._logger.debug("##############################")