Examples/at linear array transducer

.gitignore

@@ -17,4 +17,7 @@ src/
 .vscode/
 *.asv
 *.h5
+*.tar.gz
+*.whl
+*.png
 *.pyc

docs/README.md

@@ -46,7 +46,7 @@ This example file steps through the process of:
 This example expects an NVIDIA GPU by default to simulate with k-Wave.
 
 To test the reconstruction on a machine with a GPU,
-set `RUN_SIMULATION` [on line 26 of `bmode_reconstruction_example.py`](https://github.com/waltsims/k-wave-python/blob/master/examples/bmode_reconstruction_example.py#L26)
+set `RUN_SIMULATION` [on line 28 of `bmode_reconstruction_example.py`](https://github.com/waltsims/k-wave-python/blob/master/examples/bmode_reconstruction_example.py#L28)
 to `True`, and the example will run without the pre-computed data.
 
 ## Development

examples/bmode_reconstruction_example.py

@@ -4,14 +4,14 @@
 import numpy as np
 import scipy.io
 
-from kwave.data import Vector
-from kwave.options import SimulationOptions, SimulationExecutionOptions
-
 from example_utils import download_from_gdrive_if_does_not_exist
+from kwave.data import Vector
 from kwave.kgrid import kWaveGrid
 from kwave.kmedium import kWaveMedium
-from kwave.kspaceFirstOrder3D import kspaceFirstOrder3DC
+from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D
 from kwave.ktransducer import NotATransducer, kWaveTransducerSimple
+from kwave.options.simulation_execution_options import SimulationExecutionOptions
+from kwave.options.simulation_options import SimulationOptions
 from kwave.reconstruction.beamform import beamform
 from kwave.reconstruction.converter import build_channel_data
 from kwave.utils.dotdictionary import dotdict
@@ -20,72 +20,38 @@
 if __name__ == '__main__':
     # pathname for the input and output files
     pathname = gettempdir()
+    phantom_data_path = 'phantom_data.mat'
+    PHANTOM_DATA_GDRIVE_ID = '1ZfSdJPe8nufZHz0U9IuwHR4chaOGAWO4'
 
     # simulation settings
     DATA_CAST = 'single'
     RUN_SIMULATION = False
 
-    # =========================================================================
-    # DEFINE THE K-WAVE GRID
-    # =========================================================================
-    print("Setting up the k-wave grid...")
-
-    # set the size of the perfectly matched layer (PML)
     pml_size_points = Vector([20, 10, 10])  # [grid points]
-
-    # set total number of grid points not including the PML
     grid_size_points = Vector([256, 128, 128]) - 2 * pml_size_points  # [grid points]
-
-    # set desired grid size in the x-direction not including the PML
     grid_size_meters = 40e-3  # [m]
-
-    # calculate the spacing between the grid points
     grid_spacing_meters = grid_size_meters / Vector([grid_size_points.x, grid_size_points.x, grid_size_points.x])
 
-    # create the k-space grid
-    kgrid = kWaveGrid(grid_size_points, grid_spacing_meters)
-
-    # =========================================================================
-    # DEFINE THE MEDIUM PARAMETERS
-    # =========================================================================
-    # define the properties of the propagation medium
     c0 = 1540
     rho0 = 1000
-
-    medium = kWaveMedium(
-        sound_speed=None,  # will be set later
-        alpha_coeff=0.75,
-        alpha_power=1.5,
-        BonA=6
-    )
-
-    # create the time array
-    t_end = (grid_size_points.x * grid_spacing_meters.x) * 2.2 / c0  # [s]
-    kgrid.makeTime(c0, t_end=t_end)
-
-    # =========================================================================
-    # DEFINE THE INPUT SIGNAL
-    # =========================================================================
-    print("Defining the input signal...")
-
-    # define properties of the input signal
     source_strength = 1e6  # [Pa]
     tone_burst_freq = 1.5e6  # [Hz]
     tone_burst_cycles = 4
 
-    # create the input signal using tone_burst
-    input_signal = tone_burst(1 / kgrid.dt, tone_burst_freq, tone_burst_cycles)
+    kgrid = kWaveGrid(grid_size_points, grid_spacing_meters)
+    t_end = (grid_size_points.x * grid_spacing_meters.x) * 2.2 / c0  # [s]
+    kgrid.makeTime(c0, t_end=t_end)
 
-    # scale the source magnitude by the source_strength divided by the
-    # impedance (the source is assigned to the particle velocity)
+    input_signal = tone_burst(1 / kgrid.dt, tone_burst_freq, tone_burst_cycles)
     input_signal = (source_strength / (c0 * rho0)) * input_signal
 
-    # =========================================================================
-    # DEFINE THE ULTRASOUND TRANSDUCER
-    # =========================================================================
-    print("Setting up the transducer configuration...")
+    medium = kWaveMedium(
+        sound_speed=None,  # will be set later
+        alpha_coeff=0.75,
+        alpha_power=1.5,
+        BonA=6
+    )
 
-    # physical properties of the transducer
     transducer = dotdict()
     transducer.number_elements = 32  # total number of transducer elements
     transducer.element_width = 2  # width of each element [grid points/voxels]
@@ -95,66 +61,48 @@
 
     # calculate the width of the transducer in grid points
     transducer_width = transducer.number_elements * transducer.element_width + (
-                transducer.number_elements - 1) * transducer.element_spacing
+            transducer.number_elements - 1) * transducer.element_spacing
 
     # use this to position the transducer in the middle of the computational grid
-    transducer.position = np.round([1, grid_size_points.y / 2 - transducer_width / 2, grid_size_points.z / 2 - transducer.element_length / 2])
+    transducer.position = np.round(
+        [1, grid_size_points.y / 2 - transducer_width / 2,
+         grid_size_points.z / 2 - transducer.element_length / 2])
+
+    transducer = kWaveTransducerSimple(kgrid, **transducer)
 
     # properties used to derive the beamforming delays
     not_transducer = dotdict()
     not_transducer.sound_speed = c0  # sound speed [m/s]
     not_transducer.focus_distance = 20e-3  # focus distance [m]
     not_transducer.elevation_focus_distance = 19e-3  # focus distance in the elevation plane [m]
     not_transducer.steering_angle = 0  # steering angle [degrees]
-
-    # apodization
     not_transducer.transmit_apodization = 'Hanning'
     not_transducer.receive_apodization = 'Rectangular'
-
-    # define the transducer elements that are currently active
     not_transducer.active_elements = np.ones((transducer.number_elements, 1))
-
-    # append input signal used to drive the transducer
     not_transducer.input_signal = input_signal
 
-    # create the transducer using the defined settings
-    transducer = kWaveTransducerSimple(kgrid, **transducer)
     not_transducer = NotATransducer(transducer, kgrid, **not_transducer)
 
-    # =========================================================================
-    # DEFINE THE MEDIUM PROPERTIES
-    # =========================================================================
-    # define a large image size to move across
     number_scan_lines = 96
 
     print("Fetching phantom data...")
-    phantom_data_path = 'phantom_data.mat'
-    PHANTOM_DATA_GDRIVE_ID = '1ZfSdJPe8nufZHz0U9IuwHR4chaOGAWO4'
     download_from_gdrive_if_does_not_exist(PHANTOM_DATA_GDRIVE_ID, phantom_data_path)
 
     phantom = scipy.io.loadmat(phantom_data_path)
     sound_speed_map = phantom['sound_speed_map']
     density_map = phantom['density_map']
 
-    # =========================================================================
-    # RUN THE SIMULATION
-    # =========================================================================
     print(f"RUN_SIMULATION set to {RUN_SIMULATION}")
-    # run the simulation if set to true, otherwise, load previous results from disk
 
-    # set medium position
+    # preallocate the storage set medium position
+    scan_lines = np.zeros((number_scan_lines, not_transducer.number_active_elements, kgrid.Nt))
     medium_position = 0
 
-    # preallocate the storage
-    simulation_data = []
-
-    # loop through the scan lines
-    for scan_line_index in range(1, number_scan_lines + 1):
-        # for scan_line_index in range(1, 10):
-        # update the command line status
+    for scan_line_index in range(0, number_scan_lines):
 
         # load the current section of the medium
-        medium.sound_speed = sound_speed_map[:, medium_position:medium_position + grid_size_points.y, :]
+        medium.sound_speed = \
+            sound_speed_map[:, medium_position:medium_position + grid_size_points.y, :]
         medium.density = density_map[:, medium_position:medium_position + grid_size_points.y, :]
 
         # set the input settings
@@ -172,20 +120,22 @@
         )
         # run the simulation
         if RUN_SIMULATION:
-            sensor_data = kspaceFirstOrder3DC(
+            sensor_data = kspaceFirstOrder3D(
                 medium=medium,
                 kgrid=kgrid,
                 source=not_transducer,
                 sensor=not_transducer,
                 simulation_options=simulation_options,
-                execution_options=SimulationExecutionOptions()
+                execution_options=SimulationExecutionOptions(is_gpu_simulation=True)
             )
 
+            scan_lines[scan_line_index, :] = not_transducer.combine_sensor_data(sensor_data['p'].T)
+
         # update medium position
         medium_position = medium_position + transducer.element_width
 
     if RUN_SIMULATION:
-        simulation_data = np.stack(simulation_data, axis=0)
+        simulation_data = scan_lines
         # scipy.io.savemat('sensor_data.mat', {'sensor_data_all_lines': simulation_data})
 
     else:

examples/example_at_linear_array_transducer.py

@@ -0,0 +1,97 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import kwave.data
+from kwave import kWaveGrid, SimulationOptions, kWaveMedium
+from kwave.ksensor import kSensor
+from kwave.ksource import kSource
+from kwave.kspaceFirstOrder3D import kspaceFirstOrder3DC
+from kwave.options.simulation_execution_options import SimulationExecutionOptions
+from kwave.utils.kwave_array import kWaveArray
+from kwave.utils.plot import voxel_plot
+from kwave.utils.signals import tone_burst
+
+# DEFINE LITERALS
+model = 1
+c0 = 1500
+rho0 = 1000
+source_f0 = 1e6
+source_amp = 1e6
+source_cycles = 5
+source_focus = 20e-3
+element_num = 15
+element_width = 1e-3
+element_length = 10e-3
+element_pitch = 2e-3
+translation = kwave.data.Vector([5e-3, 0, 8e-3])
+rotation = kwave.data.Vector([0, 20, 0])
+grid_size_x = 40e-3
+grid_size_y = 20e-3
+grid_size_z = 40e-3
+ppw = 3
+t_end = 35e-6
+cfl = 0.5
+
+# GRID
+dx = c0 / (ppw * source_f0)
+Nx = round(grid_size_x / dx)
+Ny = round(grid_size_y / dx)
+Nz = round(grid_size_z / dx)
+kgrid = kWaveGrid([Nx, Ny, Nz], [dx, dx, dx])
+kgrid.makeTime(c0, cfl, t_end)
+
+# SOURCE
+if element_num % 2 != 0:
+    ids = np.arange(1, element_num + 1) - np.ceil(element_num / 2)
+else:
+    ids = np.arange(1, element_num + 1) - (element_num + 1) / 2
+
+time_delays = -(np.sqrt((ids * element_pitch) ** 2 + source_focus ** 2) - source_focus) / c0
+time_delays = time_delays - min(time_delays)
+
+source_sig = source_amp * tone_burst(1 / kgrid.dt, source_f0, source_cycles,
+                                     signal_offset=np.round(time_delays / kgrid.dt).astype(int))
+karray = kWaveArray(bli_tolerance=0.05, upsampling_rate=10)
+
+for ind in range(1, element_num + 1):
+    x_pos = 0 - (element_num * element_pitch / 2 - element_pitch / 2) + (ind - 1) * element_pitch
+    karray.add_rect_element([x_pos, 0, kgrid.z_vec[0]], element_width, element_length, rotation)
+
+karray.set_array_position(translation, rotation)
+source = kSource()
+source.p_mask = karray.get_array_binary_mask(kgrid)
+voxel_plot(np.single(source.p_mask))
+source.p = karray.get_distributed_source_signal(kgrid, source_sig)
+# MEDIUM
+
+medium = kWaveMedium(sound_speed=c0, density=rho0)
+
+# SENSOR
+sensor_mask = np.zeros((Nx, Ny, Nz))
+sensor_mask[:, Ny // 2, :] = 1
+sensor = kSensor(sensor_mask, record=['p_max'])
+
+# SIMULATION
+simulation_options = SimulationOptions(
+    pml_auto=True,
+    pml_inside=False,
+    save_to_disk=True,
+    data_cast='single',
+)
+
+execution_options = SimulationExecutionOptions(is_gpu_simulation=True)
+
+sensor_data = kspaceFirstOrder3DC(kgrid=kgrid, medium=medium, source=source, sensor=sensor,
+                                  simulation_options=simulation_options, execution_options=execution_options)
+
+p_max = np.reshape(sensor_data['p_max'], (Nx, Nz), order='F')
+
+# VISUALISATION
+plt.figure()
+plt.imshow(1e-6 * p_max, extent=[1e3 * kgrid.x_vec[0][0], 1e3 * kgrid.x_vec[-1][0], 1e3 * kgrid.z_vec[0][0],
+                                 1e3 * kgrid.z_vec[-1][0]], aspect='auto')
+plt.xlabel('z-position [mm]')
+plt.ylabel('x-position [mm]')
+plt.title('Pressure Field')
+plt.colorbar(label='[MPa]')
+plt.show()