Coverage for kwave/kspaceFirstOrder2D.py: 20%
93 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-03-24 12:06 -0700
1from typing import Union
3import numpy as np
4from beartype import beartype as typechecker
6from kwave.executor import Executor
7from kwave.kgrid import kWaveGrid
8from kwave.kmedium import kWaveMedium
9from kwave.ksensor import kSensor
10from kwave.ksource import kSource
11from kwave.ktransducer import NotATransducer
12from kwave.kWaveSimulation import kWaveSimulation
13from kwave.kWaveSimulation_helper import retract_transducer_grid_size, save_to_disk_func
14from kwave.options.simulation_execution_options import SimulationExecutionOptions
15from kwave.options.simulation_options import SimulationOptions
16from kwave.utils.dotdictionary import dotdict
17from kwave.utils.interp import interpolate2d
18from kwave.utils.pml import get_pml
19from kwave.utils.tictoc import TicToc
22@typechecker
23def kspace_first_order_2d_gpu(
24 kgrid: kWaveGrid,
25 source: kSource,
26 sensor: Union[NotATransducer, kSensor, None],
27 medium: kWaveMedium,
28 simulation_options: SimulationOptions,
29 execution_options: SimulationExecutionOptions,
30) -> np.ndarray:
31 """
32 2D ime-domain simulation of wave propagation on a GPU using C++ CUDA code.
34 kspaceFirstOrder2DG provides a blind interface to the C++/CUDA
35 version of kspaceFirstOrder2D (called kspaceFirstOrder-CUDA) in the
36 same way as kspaceFirstOrder3DC. Note, the C++ code does not support
37 all input options, and all display options are ignored (only command
38 line outputs are given). See the k-Wave user manual for more
39 information.
40 The function works by appending the optional input 'save_to_disk' to
41 the user inputs and then calling kspaceFirstOrder2D to save the input
42 files to disk. The contents of sensor.record (if set) are parsed as
43 input flags, and the C++ code is run using the system command. The
44 output files are then automatically loaded from disk and returned in
45 the same fashion as kspaceFirstOrder2D. The input and output files
46 are saved to the temporary directory native to the operating system,
47 and are deleted after the function runs.
48 This function requires the C++ binary/executable of
49 kspaceFirstOrder-CUDA to be downloaded from
50 http://www.k-wave.org/download.php and placed in the "binaries"
51 directory of the k-Wave toolbox (the 2D and 3D code use the same
52 binary). Alternatively, the name and location of the binary can be
53 specified using the optional input parameters 'BinaryName' and
54 'BinariesPath'.
56 This function is essentially a wrapper and directly uses the capabilities
57 of kspaceFirstOrder3DC by replacing the binary name with the name of the
58 GPU binary.
59 """
60 execution_options.is_gpu_simulation = True # force to GPU
61 assert isinstance(kgrid, kWaveGrid), "kgrid must be a kWaveGrid object"
62 assert isinstance(medium, kWaveMedium), "medium must be a kWaveMedium object"
63 assert isinstance(simulation_options, SimulationOptions), "simulation_options must be a SimulationOptions object"
64 assert isinstance(execution_options, SimulationExecutionOptions), "execution_options must be a SimulationExecutionOptions object"
66 sensor_data = kspaceFirstOrder2D(
67 kgrid=kgrid, source=source, sensor=sensor, medium=medium, simulation_options=simulation_options, execution_options=execution_options
68 ) # pass inputs to CPU version
69 return sensor_data
72def kspaceFirstOrder2DC(
73 kgrid: kWaveGrid,
74 source: kSource,
75 sensor: Union[NotATransducer, kSensor, None],
76 medium: kWaveMedium,
77 simulation_options: SimulationOptions,
78 execution_options: SimulationExecutionOptions,
79):
80 """
81 2D time-domain simulation of wave propagation using C++ code.
83 kspaceFirstOrder2DC provides a blind interface to the C++ version of
84 kspaceFirstOrder2D (called kspaceFirstOrder-OMP) in the same way as
85 kspaceFirstOrder3DC. Note, the C++ code does not support all input
86 options, and all display options are ignored (only command line
87 outputs are given). See the k-Wave user manual for more information.
88 The function works by appending the optional input 'save_to_disk' to
89 the user inputs and then calling kspaceFirstOrder2D to save the input
90 files to disk. The contents of sensor.record (if set) are parsed as
91 input flags, and the C++ code is run using the system command. The
92 output files are then automatically loaded from disk and returned in
93 the same fashion as kspaceFirstOrder2D. The input and output files
94 are saved to the temporary directory native to the operating system,
95 and are deleted after the function runs.
96 For small simulations, running the simulation on a smaller number of
97 cores can improve performance as the matrices are often small enough
98 to fit within cache. It is recommended to adjust the value of
99 'NumThreads' to optimise performance for a given simulation size and
100 computer hardware. By default, simulations smaller than 128^2 are
101 set to run using a single thread (this behaviour can be over-ridden
102 using the 'NumThreads' option). In some circumstances, for very small
103 simulations, the C++ code can be slower than the MATLAB code.
104 This function requires the C++ binary/executable of
105 kspaceFirstOrder-OMP to be downloaded from
106 http://www.k-wave.org/download.php and placed in the "binaries"
107 directory of the k-Wave toolbox (the same binary is used for
108 simulations in 2D, 3D, and axisymmetric coordinates). Alternatively,
109 the name and location of the binary can be specified using the
110 optional input parameters 'BinaryName' and 'BinariesPath'.
112 This function is essentially a wrapper and directly uses the capabilities
113 of kspaceFirstOrder3DC by replacing the binary name with the name of the
114 GPU binary.
116 Args:
117 kgrid: kWaveGrid instance
118 source: kWaveSource instance
119 sensor: NotATransducer or kSensor instance or None
120 medium: kWaveMedium instance
121 simulation_options: SimulationOptions instance
122 execution_options: SimulationExecutionOptions instance
124 Returns:
125 Sensor data as a numpy array
126 """
127 execution_options.is_gpu_simulation = False # force to CPU
128 # generate the input file and save to disk
129 sensor_data = kspaceFirstOrder2D(
130 kgrid=kgrid, source=source, sensor=sensor, medium=medium, simulation_options=simulation_options, execution_options=execution_options
131 )
132 return sensor_data
135@typechecker
136def kspaceFirstOrder2D(
137 kgrid: kWaveGrid,
138 source: kSource,
139 sensor: Union[NotATransducer, kSensor, None],
140 medium: kWaveMedium,
141 simulation_options: SimulationOptions,
142 execution_options: SimulationExecutionOptions,
143):
144 """
145 2D time-domain simulation of wave propagation.
147 kspaceFirstOrder2D simulates the time-domain propagation of
148 compressional waves through a two-dimensional homogeneous or
149 heterogeneous acoustic medium given four input structures: kgrid,
150 medium, source, and sensor. The computation is based on a first-order
151 k-space model which accounts for power law absorption and a
152 heterogeneous sound speed and density. If medium.BonA is specified,
153 cumulative nonlinear effects are also modelled. At each time-step
154 (defined by kgrid.dt and kgrid.Nt or kgrid.t_array), the acoustic
155 field parameters at the positions defined by sensor.mask are recorded
156 and stored. If kgrid.t_array is set to 'auto', this array is
157 automatically generated using the makeTime method of the kWaveGrid
158 class. An anisotropic absorbing boundary layer called a perfectly
159 matched layer (PML) is implemented to prevent waves that leave one
160 side of the domain being reintroduced from the opposite side (a
161 consequence of using the FFT to compute the spatial derivatives in
162 the wave equation). This allows infinite domain simulations to be
163 computed using small computational grids.
165 For a homogeneous medium the formulation is exact and the time-steps
166 are only limited by the effectiveness of the perfectly matched layer.
167 For a heterogeneous medium, the solution represents a leap-frog
168 pseudospectral method with a k-space correction that improves the
169 accuracy of computing the temporal derivatives. This allows larger
170 time-steps to be taken for the same level of accuracy compared to
171 conventional pseudospectral time-domain methods. The computational
172 grids are staggered both spatially and temporally.
174 An initial pressure distribution can be specified by assigning a
175 matrix (the same size as the computational grid) of arbitrary numeric
176 values to source.p0. A time varying pressure source can similarly be
177 specified by assigning a binary matrix (i.e., a matrix of 1's and 0's
178 with the same dimensions as the computational grid) to source.p_mask
179 where the 1's represent the grid points that form part of the source.
180 The time varying input signals are then assigned to source.p. This
181 can be a single time series (in which case it is applied to all
182 source elements), or a matrix of time series following the source
183 elements using MATLAB's standard column-wise linear matrix index
184 ordering. A time varying velocity source can be specified in an
185 analogous fashion, where the source location is specified by
186 source.u_mask, and the time varying input velocity is assigned to
187 source.ux and source.uy.
189 The field values are returned as arrays of time series at the sensor
190 locations defined by sensor.mask. This can be defined in three
191 different ways. (1) As a binary matrix (i.e., a matrix of 1's and 0's
192 with the same dimensions as the computational grid) representing the
193 grid points within the computational grid that will collect the data.
194 (2) As the grid coordinates of two opposing corners of a rectangle in
195 the form [x1; y1; x2; y2]. This is equivalent to using a binary
196 sensor mask covering the same region, however, the output is indexed
197 differently as discussed below. (3) As a series of Cartesian
198 coordinates within the grid which specify the location of the
199 pressure values stored at each time step. If the Cartesian
200 coordinates don't exactly match the coordinates of a grid point, the
201 output values are calculated via interpolation. The Cartesian points
202 must be given as a 2 by N matrix corresponding to the x and y
203 positions, respectively, where the Cartesian origin is assumed to be
204 in the center of the grid. If no output is required, the sensor input
205 can be replaced with an empty array [].
207 If sensor.mask is given as a set of Cartesian coordinates, the
208 computed sensor_data is returned in the same order. If sensor.mask is
209 given as a binary matrix, sensor_data is returned using MATLAB's
210 standard column-wise linear matrix index ordering. In both cases, the
211 recorded data is indexed as sensor_data(sensor_point_index,
212 time_index). For a binary sensor mask, the field values at a
213 particular time can be restored to the sensor positions within the
214 computation grid using unmaskSensorData. If sensor.mask is given as a
215 list of opposing corners of a rectangle, the recorded data is indexed
216 as sensor_data(rect_index).p(x_index, y_index, time_index), where
217 x_index and y_index correspond to the grid index within the
218 rectangle, and rect_index corresponds to the number of rectangles if
219 more than one is specified.
221 By default, the recorded acoustic pressure field is passed directly
222 to the output sensor_data. However, other acoustic parameters can
223 also be recorded by setting sensor.record to a cell array of the form
224 {'p', 'u', 'p_max', ...}. For example, both the particle velocity and
225 the acoustic pressure can be returned by setting sensor.record =
226 {'p', 'u'}. If sensor.record is given, the output sensor_data is
227 returned as a structure with the different outputs appended as
228 structure fields. For example, if sensor.record = {'p', 'p_final',
229 'p_max', 'u'}, the output would contain fields sensor_data.p,
230 sensor_data.p_final, sensor_data.p_max, sensor_data.ux, and
231 sensor_data.uy. Most of the output parameters are recorded at the
232 given sensor positions and are indexed as
233 sensor_data.field(sensor_point_index, time_index) or
234 sensor_data(rect_index).field(x_index, y_index, time_index) if using
235 a sensor mask defined as opposing rectangular corners. The exceptions
236 are the averaged quantities ('p_max', 'p_rms', 'u_max', 'p_rms',
237 'I_avg'), the 'all' quantities ('p_max_all', 'p_min_all',
238 'u_max_all', 'u_min_all'), and the final quantities ('p_final',
239 'u_final'). The averaged quantities are indexed as
240 sensor_data.p_max(sensor_point_index) or
241 sensor_data(rect_index).p_max(x_index, y_index) if using rectangular
242 corners, while the final and 'all' quantities are returned over the
243 entire grid and are always indexed as sensor_data.p_final(nx, ny),
244 regardless of the type of sensor mask.
246 kspaceFirstOrder2D may also be used for time reversal image
247 reconstruction by assigning the time varying pressure recorded over
248 an arbitrary sensor surface to the input field
249 sensor.time_reversal_boundary_data. This data is then enforced in
250 time reversed order as a time varying Dirichlet boundary condition
251 over the sensor surface given by sensor.mask. The boundary data must
252 be indexed as sensor.time_reversal_boundary_data(sensor_point_index,
253 time_index). If sensor.mask is given as a set of Cartesian
254 coordinates, the boundary data must be given in the same order. An
255 equivalent binary sensor mask (computed using nearest neighbour
256 interpolation) is then used to place the pressure values into the
257 computational grid at each time step. If sensor.mask is given as a
258 binary matrix of sensor points, the boundary data must be ordered
259 using MATLAB's standard column-wise linear matrix indexing. If no
260 additional inputs are required, the source input can be replaced with
261 an empty array [].
263 Acoustic attenuation compensation can also be included during time
264 reversal image reconstruction by assigning the absorption parameters
265 medium.alpha_coeff and medium.alpha_power and reversing the sign of
266 the absorption term by setting medium.alpha_sign = [-1, 1]. This
267 forces the propagating waves to grow according to the absorption
268 parameters instead of decay. The reconstruction should then be
269 regularised by assigning a filter to medium.alpha_filter (this can be
270 created using getAlphaFilter).
272 Note: To run a simple photoacoustic image reconstruction example
273 using time reversal (that commits the 'inverse crime' of using the
274 same numerical parameters and model for data simulation and image
275 reconstruction), the sensor_data returned from a k-Wave simulation
276 can be passed directly to sensor.time_reversal_boundary_data with the
277 input fields source.p0 and source.p removed or set to zero.
279 Args:
280 kgrid: kWaveGrid instance
281 medium: kWaveMedium instance
282 source: kWaveSource instance
283 sensor: kWaveSensor instance or None
285 Returns:
287 """
288 # start the timer and store the start time
289 TicToc.tic()
291 # Currently we only support binary execution, meaning all simulations must be saved to disk.
292 if not simulation_options.save_to_disk:
293 if execution_options.is_gpu_simulation:
294 raise ValueError("GPU simulation requires saving to disk. Please set SimulationOptions.save_to_disk=True")
295 else:
296 raise ValueError("CPU simulation requires saving to disk. Please set SimulationOptions.save_to_disk=True")
298 k_sim = kWaveSimulation(kgrid=kgrid, source=source, sensor=sensor, medium=medium, simulation_options=simulation_options)
300 k_sim.input_checking("kspaceFirstOrder2D")
302 # =========================================================================
303 # CALCULATE MEDIUM PROPERTIES ON STAGGERED GRID
304 # =========================================================================
305 options = k_sim.options
307 # interpolate the values of the density at the staggered grid locations
308 # where sgx = (x + dx/2, y, z), sgy = (x, y + dy/2, z), sgz = (x, y, z + dz/2)
309 k_sim.rho0 = np.atleast_1d(k_sim.rho0)
310 if k_sim.rho0.ndim == 2 and options.use_sg:
311 # rho0 is heterogeneous and staggered grids are used
312 grid_points = [k_sim.kgrid.x, k_sim.kgrid.y]
313 k_sim.rho0_sgx = interpolate2d(grid_points, k_sim.rho0, [k_sim.kgrid.x + k_sim.kgrid.dx / 2, k_sim.kgrid.y])
314 k_sim.rho0_sgy = interpolate2d(grid_points, k_sim.rho0, [k_sim.kgrid.x, k_sim.kgrid.y + k_sim.kgrid.dy / 2])
315 else:
316 # rho0 is homogeneous or staggered grids are not used
317 k_sim.rho0_sgx = k_sim.rho0
318 k_sim.rho0_sgy = k_sim.rho0
319 k_sim.rho0_sgz = None
321 # invert rho0 so it doesn't have to be done each time step
322 k_sim.rho0_sgx_inv = 1 / k_sim.rho0_sgx
323 k_sim.rho0_sgy_inv = 1 / k_sim.rho0_sgy
325 # clear unused variables if not using them in _saveToDisk
326 if not options.save_to_disk:
327 del k_sim.rho0_sgx
328 del k_sim.rho0_sgy
330 # =========================================================================
331 # PREPARE DERIVATIVE AND PML OPERATORS
332 # =========================================================================
334 # get the PML operators based on the reference sound speed and PML settings
335 Nx, Ny = k_sim.kgrid.Nx, k_sim.kgrid.Ny
336 dx, dy = k_sim.kgrid.dx, k_sim.kgrid.dy
337 dt = k_sim.kgrid.dt
338 pml_x_alpha, pml_y_alpha = options.pml_x_alpha, options.pml_y_alpha
339 pml_x_size, pml_y_size = options.pml_x_size, options.pml_y_size
340 c_ref = k_sim.c_ref
342 k_sim.pml_x = get_pml(Nx, dx, dt, c_ref, pml_x_size, pml_x_alpha, False, 1)
343 k_sim.pml_x_sgx = get_pml(Nx, dx, dt, c_ref, pml_x_size, pml_x_alpha, True and options.use_sg, 1)
344 k_sim.pml_y = get_pml(Ny, dy, dt, c_ref, pml_y_size, pml_y_alpha, False, 2)
345 k_sim.pml_y_sgy = get_pml(Ny, dy, dt, c_ref, pml_y_size, pml_y_alpha, True and options.use_sg, 2)
347 # define the k-space derivative operators, multiply by the staggered
348 # grid shift operators, and then re-order using ifftshift (the option
349 # flgs.use_sg exists for debugging)
350 kx_vec, ky_vec = k_sim.kgrid.k_vec
351 kx_vec, ky_vec = np.array(kx_vec), np.array(ky_vec)
352 if options.use_sg:
353 k_sim.ddx_k_shift_pos = np.fft.ifftshift(1j * kx_vec * np.exp(1j * kx_vec * dx / 2))[None, :]
354 k_sim.ddx_k_shift_neg = np.fft.ifftshift(1j * kx_vec * np.exp(-1j * kx_vec * dx / 2))[None, :]
355 k_sim.ddy_k_shift_pos = np.fft.ifftshift(1j * ky_vec * np.exp(1j * ky_vec * dy / 2))[None, :]
356 k_sim.ddy_k_shift_neg = np.fft.ifftshift(1j * ky_vec * np.exp(-1j * ky_vec * dy / 2))[None, :]
357 else:
358 k_sim.ddx_k_shift_pos = np.fft.ifftshift(1j * kx_vec)[None, :]
359 k_sim.ddx_k_shift_neg = np.fft.ifftshift(1j * kx_vec)[None, :]
360 k_sim.ddy_k_shift_pos = np.fft.ifftshift(1j * ky_vec)[None, :]
361 k_sim.ddy_k_shift_neg = np.fft.ifftshift(1j * ky_vec)[None, :]
363 # force the derivative and shift operators to be in the correct direction for use with BSXFUN
364 k_sim.ddy_k_shift_pos = k_sim.ddy_k_shift_pos.T
365 k_sim.ddy_k_shift_neg = k_sim.ddy_k_shift_neg.T
367 # create k-space operators (the option flgs.use_kspace exists for debugging)
368 if options.use_kspace:
369 k = k_sim.kgrid.k
370 k_sim.kappa = np.fft.ifftshift(np.sinc(c_ref * k * dt / 2))
371 if (k_sim.source_p and k_sim.source.p_mode == "additive") or (
372 (k_sim.source_ux or k_sim.source_uy or k_sim.source_uz) and k_sim.source.u_mode == "additive"
373 ):
374 k_sim.source_kappa = np.fft.ifftshift(np.cos(c_ref * k * dt / 2))
375 else:
376 k_sim.kappa = 1
377 k_sim.source_kappa = 1
379 # =========================================================================
380 # SAVE DATA TO DISK FOR RUNNING SIMULATION EXTERNAL TO MATLAB
381 # =========================================================================
383 # save to disk option for saving the input matrices to disk for running
384 # simulations using k-Wave++
385 if options.save_to_disk:
386 # store the pml size for resizing transducer object below
387 retract_size = [[options.pml_x_size, options.pml_y_size, options.pml_z_size]]
389 # run subscript to save files to disk
390 save_to_disk_func(
391 k_sim.kgrid,
392 k_sim.medium,
393 k_sim.source,
394 k_sim.options,
395 execution_options.auto_chunking,
396 dotdict(
397 {
398 "ddx_k_shift_pos": k_sim.ddx_k_shift_pos,
399 "ddx_k_shift_neg": k_sim.ddx_k_shift_neg,
400 "dt": k_sim.dt,
401 "c0": k_sim.c0,
402 "c_ref": k_sim.c_ref,
403 "rho0": k_sim.rho0,
404 "rho0_sgx": k_sim.rho0_sgx,
405 "rho0_sgy": k_sim.rho0_sgy,
406 "rho0_sgz": k_sim.rho0_sgz,
407 "p_source_pos_index": k_sim.p_source_pos_index,
408 "u_source_pos_index": k_sim.u_source_pos_index,
409 "s_source_pos_index": k_sim.s_source_pos_index,
410 "transducer_input_signal": k_sim.transducer_input_signal,
411 "delay_mask": k_sim.delay_mask,
412 "sensor_mask_index": k_sim.sensor_mask_index,
413 "record": k_sim.record,
414 }
415 ),
416 dotdict(
417 {
418 "source_p": k_sim.source_p,
419 "source_p0": k_sim.source_p0,
420 "source_ux": k_sim.source_ux,
421 "source_uy": k_sim.source_uy,
422 "source_uz": k_sim.source_uz,
423 "source_sxx": k_sim.source_sxx,
424 "source_syy": k_sim.source_syy,
425 "source_szz": k_sim.source_szz,
426 "source_sxy": k_sim.source_sxy,
427 "source_sxz": k_sim.source_sxz,
428 "source_syz": k_sim.source_syz,
429 "transducer_source": k_sim.transducer_source,
430 "nonuniform_grid": k_sim.nonuniform_grid,
431 "elastic_code": k_sim.options.simulation_type.is_elastic_simulation(),
432 "axisymmetric": k_sim.options.simulation_type.is_axisymmetric(),
433 "cuboid_corners": k_sim.cuboid_corners,
434 }
435 ),
436 )
438 # run subscript to resize the transducer object if the grid has been expanded
439 retract_transducer_grid_size(k_sim.source, k_sim.sensor, retract_size, k_sim.options.pml_inside)
441 # exit matlab computation if required
442 if options.save_to_disk_exit:
443 return
445 executor = Executor(simulation_options=simulation_options, execution_options=execution_options)
446 executor_options = execution_options.as_list(sensor=k_sim.sensor)
447 sensor_data = executor.run_simulation(k_sim.options.input_filename, k_sim.options.output_filename, options=executor_options)
448 return sensor_data