# Classes for Detector, Absorber, and PA Simulator
class Detector:
def __init__(self,position):
self.position = position
class Absorber:
def __init__(self, position, radius, absorption_coeff):
self.position = position # (z,x)
self.radius = radius
self.absorption_coeff = absorption_coeff
class PA_simulator:
def __init__(self,time,z_length,x_length,dz,laser_energy,pulse_width,sound_speed,Gruneisen,beam_focus_z,beam_width,wavelength,t0):
self.time=time
self.z_length=z_length
self.x_length=x_length
self.dz=dz
self.laser_energy = laser_energy
self.pulse_width = pulse_width
self.sound_speed = sound_speed
self.Gruneisen = Gruneisen
self.beam_focus_z = beam_focus_z
self.beam_width = beam_width
self.t0=t0
self.wavelength=wavelength
self.absorbers = []
self.detectors=[]
def add_absorber(self, absorber):
self.absorbers.append(absorber)
def add_detector(self, detector):
self.detectors.append(detector)
def laser_pulse(self):
sigma = self.pulse_width
time=self.time
t0=self.t0
# Area normalized
return (1/(sigma*np.sqrt(2*np.pi)))*np.exp(-(time - t0)**2/(2 * sigma**2))
def initial_pressure(self,absorbed_energy):
time=self.time
pulse = self.laser_pulse()
dt = np.abs(time[1] - time[0])
dpulse_dt = np.gradient(pulse, dt)
p0= absorbed_energy* self.Gruneisen * dpulse_dt
return p0
def distance(self,p1,p2):
return np.sqrt((p1[0] - p2[0])**2 +(p1[1] - p2[1])**2)
def propagate(self,p1,p2,absorbed_energy):
time=self.time
dt = np.abs(time[1] - time[0])
r=self.distance(p1,p2)
arrival_time = r / self.sound_speed # Time of arrival of the PA signal in microseconds
shift = int(arrival_time / dt)
pressure=self.initial_pressure(absorbed_energy)
propagated = np.zeros_like(pressure)
if shift < len(time):
propagated[shift:] = pressure[:-shift]
propagated /= max(r, 1e-12) # Devided by distance to account for spherical spreading, and avoid division by zero
return propagated
def create_SpaceGrid(self):
x_length=self.x_length
z_length=self.z_length
dx=self.dz
dz=dx
x=np.linspace(0,x_length,int(round(x_length/dx)))
z=np.linspace(0,z_length,int(round(z_length/dz)))
return z,x
# Function to create absorber map
def absorber_map_function(self):
z,x=self.create_SpaceGrid()
X, Z = np.meshgrid(x, z)
absorber_map = np.zeros_like(X)
for absorber in self.absorbers:
x0=absorber.position[1]
z0=absorber.position[0]
mask = ((X-x0)**2 + (Z-z0)**2) <= absorber.radius**2
absorber_map += mask.astype(float)*absorber.absorption_coeff
return absorber_map
def propagate_beam(self):
# Length parameters
z,x=self.create_SpaceGrid()
dz=self.dz
# Initialize beam
beam=np.zeros((len(z),len(x)),dtype=complex)
# Wavevector calculations
wavelength=self.wavelength
wavelength=wavelength*1e-3 # Convert wavelength to mm (from μm)
k=2*np.pi/wavelength # Wave number in radians/mm
kx=fftfreq(len(x),np.diff(x)[0])*2*np.pi # Spatial angular frequency in radians/mm
# Laser paramters
pos_z=self.beam_focus_z
pos_x=self.x_length*0.5 # This make sure that the beam is in the middle of x range
E=np.sqrt(self.laser_energy)
w0=self.beam_width # Beam width
z0=np.pi*w0**2/wavelength # Rayleigh range
zz=-max(pos_z,1e-10) # Distance of z=0 from focused position
W=w0*np.sqrt(1+(zz/z0)**2) # Beam width at z=0, z-z_position
R=zz*(1+(z0/zz)**2)
phase=np.exp(1j*(-k*zz-k*((x-pos_x)**2)*0.5/R+np.atan(zz/z0)))
# Intialize beam for all x at z=0, as a Gaussian beam
beam[0,:]=E/(W*np.sqrt(np.pi))*np.exp(-((x-pos_x)/W)**2)*phase
absorber_map_val=self.absorber_map_function()
# Initialize the heat map which captured energy absorbed * absorption coefficient
heat_map=np.zeros((len(z),len(x)))
for zi in range(len(z)-1):
absorber_x=absorber_map_val[zi,:] # Absorber (coefficient) array along x direction
beam_attenuated=beam[zi,:]*np.exp(-absorber_x*dz*0.5) # 1/2 factor, as the equations are in amplitude
local_energy=(np.abs(beam[zi,:]))**2 # Energy
heat_map[zi+1,:]=absorber_x*(1-np.exp(-local_energy*absorber_x))# Energt absorbed X absorption coefficient
F=np.fft.fft(beam_attenuated) # Beam in spatial frequency
beam[zi+1,:]=np.fft.ifft(F*np.exp(-1j*dz*np.sqrt(k**2-kx**2+0j)))# Propagated by multiplying e^ik_z dz and taking FFT
return (np.abs(beam))**2,heat_map
def detector_signal(self,pos_detector):
time=self.time
total_signal = np.zeros_like(time)
z,x=self.create_SpaceGrid()
beam,heat_map=self.propagate_beam()
z_idx, x_idx = np.where(heat_map > 0)
for indx_z, indx_x in zip(z_idx, x_idx):
pos_absorber=[z[indx_z], x[indx_x]]
# print(pos_absorber)
absorbed_energy = heat_map[indx_z, indx_x] # Energy X absorption_coefficient
# print(pos_absorber)
propagated_signal=self.propagate(pos_absorber,pos_detector,absorbed_energy)
total_signal += propagated_signal
return total_signal
def all_detector_signal(self):
signals=[]
No_detector=0
for detector in self.detectors:
signal = self.detector_signal(detector.position)
signals.append(signal)
No_detector +=1
return No_detector, signals
def plot_2D_function(self, u,plot_title):
z, x = self.create_SpaceGrid()
plt.figure(figsize=(5, 2))
plt.pcolormesh(z, x, u.T, shading='auto')
plt.xlabel('Optic axis, z (mm)',fontsize=10, fontweight='bold')
plt.ylabel('x (mm)',fontsize=10, fontweight='bold')
plt.title(plot_title,fontsize=10, fontweight='bold')
# plt.colorbar(label='Amplitude')
plt.tight_layout()
plt.show()Photoacoustic signal generation and reconstruction algorithm
Forward problem: Photoacoustic signal generation
Process
Laser energy absorption (fluence X absorption coefficient X volume) \(\Delta Q(\mathrm{r})=\left(\int F(\mathrm{r},t) dt\right) \mu(\mathrm{r}) V\)
- Laser function \(F(\mathrm{r},t)=A(\mathrm{r}) I(t)\), \(A(\mathrm{r})\): fluence, I(t): temporal function (area normalized)
\(\to\) Temperature rise, without diffusion for maximum signal generation (1/specific heat X 1/(density X volume)) \(\Delta T=(1/C_p \rho V) \Delta Q\)
- Assuming diffusion is negligible, true for pulse width \(t_{pulse}\ll \frac{d^{2}}{D}\) (D \(\sim\) \(1.4\times 10^{-7}\) \(\text{m}^{2}\)/s),d: dimension of absorption target
\(\to\) Pressure rise, without volume expansion for maximum signal generation \(\Delta p= (\beta / k) \Delta T\)
- \(\beta\): Isothermal compressibility (\(\sim\) \(5\times 10^{-10}\) \(\text{Pa}^{-1}\), for water)
- \(k\): Thermal coefficient of volume expansion (\(\sim\) \(4\times 10^{-4}\) \(\text{K}^{-1}\), for water)
Pressure amplitude at the source
\[ p(\mathbf{r})_0=\frac{\beta}{k\ \rho C_{p}}\mu(\mathbf{r})A(\mathbf{r})=\frac{\beta \ v_{s}^{2}}{C_{P}}\mu(\mathbf{r})A(\mathbf{r}) \] Gruneisen parameter \[\Gamma=\frac{\beta \ v_{s}^{2}}{C_{P}}\]
- \(\Gamma\) \(\sim\) 0.1, for water
The pressure wave equation
\[\left(\nabla^{2}-\frac{1}{v_s^{2}} \frac{ \partial^{2} }{ \partial t^{2} } \right) p(\mathbf{\mathrm{r},t})=-\frac{\Gamma}{v_s^2}A(\mathbf{r})\frac{ \partial I(t) }{ \partial t } \]
Laser propagation function
\[ A(z,x)=IFT[FT[A(z=0,x)]\exp\left(-j2\pi z\sqrt{ \lambda^{-2}-\nu_{x}^{2}}\right)]\]
# Simulation initialization
sim = PA_simulator(
time=np.linspace(-2,10,1000), # Time over which simulation is carried out
z_length=2, # Length of z (which is optics axis)
x_length=0.5, # One of the direction perpendicular to optics axis
dz=0.001,
laser_energy=10.0e-6, # Intial pulse energy in μJ
pulse_width=60e-3, # 60 ns pulse, unit of μs
sound_speed=1.5, # in mm/μs
Gruneisen=0.1, # for water, Dimensionless
beam_focus_z=0, # the spot where the beam is focused
beam_width=40e-3, # in mm
wavelength=1.064, # Wavelength of light in μm
t0=0.0
)
# Absorption coefficient in the unit of mm-1
# sim.add_absorber(
# Absorber(position=(1,0.25), radius=0.01, absorption_coeff=12))
sim.add_absorber(
Absorber(position=(0.5,0.25), radius=0.02, absorption_coeff=8))
x_length=sim.x_length
Detecto_length=2
N_detector=50
spacing_detector=Detecto_length/N_detector
for i in range(N_detector):
sim.add_detector( Detector(position=(0,0.5*(x_length-Detecto_length)+(i-0)*spacing_detector)))u,heat_map=sim.propagate_beam()
sim.plot_2D_function(u,"Beam propagation")
sim.plot_2D_function(heat_map,"Energy absorbed")
# Generate the signals for the detectors and store
z, x = sim.create_SpaceGrid()
t = sim.time
No_detector, signals = sim.all_detector_signal()
signals_hilbert = np.abs(hilbert(signals, axis=1))
# Detector positions
det_pos = np.array([det.position for det in sim.detectors])
sound_speed = sim.sound_speedText(0.5, 1.0, 'PA signal spectra')
Reconstruction algorithm
Delay and sum algorithm
\[I_{DAS}(\mathrm{r})=\sum_{n=0}^N S_m (\lvert r_n-r \rvert/v_s)\] - \(I_{DAS}\): image value at position \(r\), - N: Total number of detectors m: Detector index, - \(S_m\): Signal at detector m - \(r_n\): Position of detector m ### Coherence factor \[ CF=\frac{\lvert \sum_{n=1}^N S_{m} \rvert^2}{N\sum_{n=1}^N\lvert S_{m} \rvert^2} \] - Image=\(I_{DAS}*CF\)
# Reconstruction algorithm multiple methods
@njit(parallel=True)
def reconstruction_all_methods(z, x, t,signals,signals_hilbert,det_pos,sound_speed):
Nz = len(z)
Nx = len(x)
Nd = len(det_pos)
image_DAS = np.zeros((Nz, Nx))
image_DAS_CF = np.zeros((Nz, Nx))
image_envelope = np.zeros((Nz, Nx))
# Parallelize over rows
for iz in prange(Nz):
z0 = z[iz]
for ix in range(Nx):
x0 = x[ix]
das_sum_bp=0.0
das_sum = 0.0
energy_sum = 0.0
energy_sum_bp=0.0
for d in range(Nd):
z_det = det_pos[d, 0]
x_det = det_pos[d, 1]
r = np.sqrt((z_det - z0)**2 +(x_det - x0)**2)
delay = r / sound_speed
# For bipolar signal
sample_bp = np.interp(delay, t, signals[d] )
das_sum_bp += sample_bp
energy_sum_bp += sample_bp * sample_bp
# For hilbert transformed
sample = np.interp(delay, t, signals_hilbert[d] )
das_sum += sample
energy_sum += sample * sample
cf = (np.abs(das_sum **2)) / (Nd * np.abs(energy_sum) + 1e-12)
cf2 = (np.abs(das_sum_bp **2)) / (Nd * np.abs(energy_sum_bp) + 1e-12)
image_DAS[iz, ix] = das_sum_bp
image_DAS_CF[iz, ix] = das_sum_bp*cf2
image_envelope[iz, ix] = das_sum*cf
return image_DAS,image_DAS_CF,image_envelope
image_DAS,image_DAS_CF,image_envelope = reconstruction_all_methods(z, x, t,signals,signals_hilbert,det_pos,sound_speed)Text(0.5, 1.0, 'With envelop')
Effect of pulse duration in image generation
Axial resolution \(R_A\) (along the direction of propagation of acoustis waves: optic axis, z for the current example): \[ R_{A}=0.88\frac{v_s}{\Delta f} \]
- Assuming Gaussian shape, -6 dB bandwidth
- Acoustic bandwidth is, \(\Delta f\sim 0.44/\Delta t\), where \(\Delta t\) is the FWHM pulse width
- Resolution=390 μm, for \(f_0\) = 5 MHz. and fractional bandwidth= 69.48 %, (M309-SU, Olympus)
@njit(parallel=True)
def reconstruction(z, x, t,signals,det_pos,sound_speed):
Nz = len(z)
Nx = len(x)
Nd = len(det_pos)
image = np.zeros((Nz, Nx))
# Parallelize over rows
for iz in prange(Nz):
z0 = z[iz]
for ix in range(Nx):
x0 = x[ix]
das_sum_bp=0.0
das_sum = 0.0
energy_sum = 0.0
for d in range(Nd):
z_det = det_pos[d, 0]
x_det = det_pos[d, 1]
r = np.sqrt((z_det - z0)**2 +(x_det - x0)**2)
delay = r / sound_speed
# For hilbert transformed
sample = np.interp(delay, t, signals[d] )
das_sum += sample
energy_sum += sample * sample
cf = (das_sum * das_sum) / (Nd * energy_sum + 1e-12)
image[iz, ix] = das_sum * cf
return image# Simulation 2 initialization
sim2 = PA_simulator(
time=np.linspace(-2,10,1000), # Time over which simulation is carried out
z_length=2, # Length of z (which is optics axis)
x_length=0.5, # One of the direction perpendicular to optics axis
dz=0.001,
laser_energy=10.0e-6, # Intial pulse energy in μJ
pulse_width=2e-3, # 60 ns pulse, unit of μs
sound_speed=1.5, # in mm/μs
Gruneisen=0.1, # for water, Dimensionless
beam_focus_z=0, # the spot where the beam is focused
beam_width=40e-3, # in mm
wavelength=1.064, # Wavelength of light in μm
t0=0.0
)
# Absorption coefficient in the unit of mm-1
# sim.add_absorber(
# Absorber(position=(1,0.25), radius=0.01, absorption_coeff=12))
sim2.add_absorber(
Absorber(position=(0.5,0.25), radius=0.01, absorption_coeff=8))
x_length=sim2.x_length
Detecto_length=2
N_detector=50
spacing_detector2=Detecto_length/N_detector
for i in range(N_detector):
sim2.add_detector( Detector(position=(0,0.5*(x_length-Detecto_length)+(i-0)*spacing_detector)))# For the second simulation Generate the signals for the detectors and store
No_detector2, signals2 = sim2.all_detector_signal()
signals2 = np.abs(hilbert(signals2, axis=1))
# Detector positions
det_pos2 = np.array([det.position for det in sim2.detectors])
sound_speed = sim.sound_speed
image2=reconstruction(z, x, t,signals2,det_pos,sound_speed)# Plot the reconstructed images
z, x = sim.create_SpaceGrid()
pulse=sim.laser_pulse()
plt.subplots(2,1, figsize=(10,7))
plt.subplot(2,1, 1)
plt.pcolormesh(z, x, ((np.abs(image_envelope))**4).T, shading='auto')
plt.xlabel('Optic axis, z (mm)',fontsize=10, fontweight='bold')
plt.ylabel('x (mm)',fontsize=10, fontweight='bold')
plt.title("Pulse width (rms) = "+str(sim.pulse_width*1e3)+ " ns",fontsize=10, fontweight='bold')
plt.subplot(2,1, 2)
plt.pcolormesh(z, x, ((np.abs(image2))**4).T, shading='auto')
plt.xlabel('Optic axis, z (mm)',fontsize=10, fontweight='bold')
plt.ylabel('x (mm)',fontsize=10, fontweight='bold')
plt.title("Pulse width (rms) = "+str(sim2.pulse_width*1e3)+ " ns",fontsize=10, fontweight='bold')
Text(0.5, 1.0, 'Pulse width (rms) = 2.0 ns')