Q:how to run a simple photoacoustic image reconstruction example using time reversal?

[vanlinus3068]

hi! This is my first time starting a discussion in this repository! I want to use the kspaceFirstOrder2D function to simulate the time-reversal algorithm, but I’m not sure how to use or define sensor.time_reversal_boundary_data. Could someone please help me with this?

Following is my code:

import cv2
import numpy as np
from scipy.io import savemat
from kwave.data import Vector
from kwave.kgrid import kWaveGrid
from kwave.kmedium import kWaveMedium
from kwave.ksensor import kSensor
from kwave.ksource import kSource
from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D
from kwave.options.simulation_execution_options import SimulationExecutionOptions
from kwave.options.simulation_options import SimulationOptions
from kwave.utils.mapgen import make_cart_circle
import matplotlib.pyplot as plt

if name == 'main':

Nx = 128
Ny = 128
dx = 0.1e-3
dy = 0.1e-3
kgrid = kWaveGrid(N=(Nx, Ny), spacing=(dx, dy))


medium = kWaveMedium(sound_speed=1500,
                     density=1000,
                     alpha_coeff=0.75,
                     alpha_power=1.5, )

img_matrix = cv2.imread('EXAMPLE.png', cv2.IMREAD_GRAYSCALE)

source = kSource()
source.p0 = (1 * img_matrix) / 255

plt.imshow(img_matrix, cmap='gray')
plt.title('Loaded Image')
plt.axis('off')
plt.show()


sensor_radius: float = 4e-3
num_sensor_points: int = 50
sensor_mask = make_cart_circle(radius=sensor_radius, num_points=num_sensor_points, plot_circle=True)
sensor = kSensor(mask=sensor_mask)
kgrid.makeTime(c=medium.sound_speed)

opts = {
    'smooth_p0': True,
    'pml_inside': False
}

simul_opts = SimulationOptions(save_to_disk=True, **opts)
exec_opts = SimulationExecutionOptions(is_gpu_simulation=True)

result = kspaceFirstOrder2D(kgrid=kgrid, source=source, sensor=sensor, medium=medium,
                            simulation_options=simul_opts, execution_options=exec_opts)


print('keys in result:', result.keys())
print('sensor data shape:', result['p'].shape)
print('dt =', result['dt'])

sensor_data = result['p']
np.savez('my_sensor_data.npz', result=sensor_data)
savemat('my_sensor_data.mat', {'sensor_data': sensor_data})

plt.imshow(sensor_data.T, cmap='viridis', aspect='auto')
plt.colorbar()
plt.title("Sensor Data Visualization")
plt.xlabel('time idx')
plt.ylabel('sensor idx')
plt.show()


sensor = kSensor(mask=sensor_mask)
sensor.time_reversal_boundary_data = sensor_data
sensor.p_mask = sensor_mask
source.p0 = np.zeros([128, 128])
source.p = np.zeros([128, 128])

input_args = {'SaveToDisk': False}
reconstructed_data = kspaceFirstOrder2D(kgrid, source, sensor, medium, simulation_options=simul_opts,
                                        execution_options=exec_opts)
initial_pressure = reconstructed_data['p']