Various implementations of a simple reconstruction procedure for acoustic resolution optoacoustic microscopy based on the virtual detector concept. Including a small literature review.

## Introduction
Acoustic resolution optoacoustic microscopy relies on scanning a spherically focused ultrasound transducer on a two dimensional plane over the volume of interest. At each individual position, the volume of interest is illuminated by a nanosecond pulsed laser resulting in the emission of an ultrasound wave by absorbing structures. To recover spatial resolution in out of focus planes, different reconstruction techniques are applied. The most commonly applied technique is called synthetic aperture focusing technique which was originally introduced as the virtual detector concept. Thereby the spherically focused transducer is reduced to a simple point detector positioned within the focus of the transducer and reconstruction is performed according to the delay and sum approach around this virtual detector point.
For the procedure, we go to each individual output voxel of the volume, calculate the delay time to neighboring virtual detector points, and sum up the corresponding signals. Here we can already see that the procedure is highly parallizable since the individual output voxels are independent from each other.
The procedure is known to work best in the region below the focal plane. Around the focus, the arcs observed in measurements deviate strongly from the ones calculated through the simplification of the virtual detector concept. Similarly, the region between transducer and focal place is more affected by the finite size of the detector. Alternative approaches to the technique include hybrid domain or model-based reconstruction.
Papers about SAFT in acoustic resolution optoacoustic microscopy:
* M. Li et al.: Improved in vivo photoacoustic microscopy based on a virtual-detector concept in Optics Letters 31 (4), 2006, DOI: [10.1364/ol.31.000474](https://doi.org/10.1364/ol.31.000474)
* Z. Deng et al.: Two-dimensional synthetic-aperture focusing technique in photoacoustic microscopy in Journal of Applied Physics 109, 2011, DOI: [10.1063/1.3585828]( https://doi.org/10.1063/1.3585828)
* J. Turner et al.: Improved optoacoustic microscopy through three-dimensional spatial impulse response synthetic aperture focusing technique in Optics Letters 39 (12), 2014, DOI: [10.1364/OL.39.003390](https://doi.org/10.1364/OL.39.003390)
## Expansions to the technique
### Coherence factor weighting
Coherence factor weighting is an extension to the SAFT technique allowing incresaed signal-to-noise ratio and improved lateral and axial resolution by taking the ratio between the coherent and incoherent sum into account. In simple terms: it looks how much the values summed up during the procedure have the same sign. If we are reconstructing along an arc of an absorber the sign should stay the same resulting in a coherence factor of 1.
* M. Li et al.: Improved in vivo photoacoustic microscopy based on a virtual-detector concept in Optics Letters 31 (4), 2006, DOI: [10.1364/ol.31.000474](https://doi.org/10.1364/ol.31.000474)
### Sensitivity field weighting
Turner et al. introduced a SAFT procedure which is supposed to reduce the bias introduced over depth. According to their publication, there is a depth bias during reconstruction. By calculating a base-line SAFT, they take out the error introduced during the reconstruction procedure.
* J. Turner et al.: Universal weighted synthetic aperture focusing technique (W-SAFT) for scanning acoustic microscopy in Optica 4 (7), 2017, DOI: [10.1364/OPTICA.4.000770](https://doi.org/10.1364/OPTICA.4.000770)
### Deconvolution after SAFT
There were several papers published about a further increase in lateral and axial resolution if a deconvolution is applied after SAFT. The most straight forward kernel to use for this would be a reconstructed microsphere.
* D. Cai et al.: In vivo deconvolution acoustic-resolution photoacoustic microscopy in three dimensions in Biomedical Optics Express 7 (2), 2016, DOI: [10.1364/BOE.7.000369](https://doi.org/10.1364/BOE.7.000369)
### Combination of SAFT and fluence compensation
Some papers covered compensation for fluence decreasing over depth during the SAFT procedure
* V. V. Perekatova et al.: Combination of virtual point detector concept and fluence compensation in acoustic resolution photoacoustic microscopy in Journal of Biomedical Optics 23 (6), 2018, DOI: [10.1117/1.JBO.23.9.091414](http://doi.org/10.1117/1.JBO.23.9.091414)
### Simulation based delay times
Another approach to overcome the limitations of the simple virtual detector concept is to simulate the finite size of the transducer surface and then simply use the peaks estimated from the simulation for SAFT instead of the simple method through the point detector.
* K. Peng et al.: Three-dimensional synthetic aperture focusing photoacoustic microscopy based on the acoustic simulation generated delay time and weighted factor in Applied Optics 59 (32), 2020, [10.1364/AO.396272](https://doi.org/10.1364/AO.396272)
### Missing expansions or literature
If there is some important paper missing in the overview please feel free to leave a comment and I will include it into the list.
## Code examples
Below are examples for
* 2D SAFT without coherence factor weighting
* 2D SAFT with coherence factor weighting
* 3D SAFT without coherence factor weighting
* 3D SAFT with coherence factor weighting
* 3D SAFT implementation for `mexcuda`
Please note that the procedures described below take a large range of neighbouring A scans into account. Consider limiting this range to a radial range depending on the opening angle of your transducer and the current depth.
Special care should be taken while defining the time offset `dr(1)` and the resolution of the dataset `dr`.
It might be of advantage to preprocess the datasets to
* Compensate for laser energy fluctuations
* Bandpass filter waveforms to the frequency range of interest
* Correct for laser jitter (fluctuations of pulsing event along time)
### Code preamble
“`MATLAB
% define file and load to workspace
% load your file to the workspace
% vol – signal matrix with [it, ix, iy]
% dr – 3 element vector for resolution as [dt, dx, dy]
% origin – 3 element vector for origin of dataset [t0, x0, y0]
% get dimensions and generate vectors
[nT, nX, nY] = size(vol); % get dimensions of volume
vecX = (0:(nX – 1)) * dr(2); % generate x vector [m]
vecY = (0:(nY – 1)) * dr(3); % generate y vector [m]
vecT = (0:(nT – 1)) * dr(1) + origin(1); % generate time vector [s]
% calculate delay times for each timepoint
SOS = 1495; % speed of sound [m/s]
fd = 7e-3; % focal distance of transducer [m]
rAp = 3.2e-3; % aperture radius of transducer [m] <-- unused
vecZ = vecT * SOS;
```
### 2D SAFT without coherence factor weighting
```MATLAB
% perform 2d saft
sliceData = squeeze(vol(:, :, iY)); % extract slice for 2D saft
sliceSAFT = zeros(size(sliceData)); % allocate memory for output slice
for idxZ = 1:nT % for each z position in the output volume
zDepth = vecZ(idxZ); % get depth of currently reconstructed layer
deltaZ = zDepth - fd; % distance between focal plane and z layer
for idxX = 1:nX % for each x position in the output volume
xRel = (1:nX) - idxX; % relative shift along x
xRel = xRel * dr(2);
for (idxRel = 1:nX)
% calculate distance between
distTot = sqrt(deltaZ * deltaZ + xRel(idxRel) * xRel(idxRel));
deltaT = round(distTot / SOS / dr(1));
if (zDepth > fd)
deltatIdx = idxFoc + deltaT;
else
deltatIdx = idxFoc – deltaT;
end
if (deltatIdx >= 1) && (deltatIdx <= nT)
sliceSAFT(idxZ, idxX) = sliceSAFT(idxZ, idxX) + sliceData(deltatIdx, idxRel);
end
end
end
end
```
### 2D SAFT with coherence factor weighting
```MATLAB
% perform 2d saft with coherence factor weighting
sliceSAFTCoh = zeros(size(sliceData));
cfTemp = 0;
saftTemp = 0;
for idxZ = 1:nT % for each z position in the output volume
zDepth = vecZ(idxZ);
deltaZ = zDepth - fd;
for idxX = 1:nX % for each x position in the output volume
xRel = (1:nX) - idxX;
xRel = xRel * dr(2);
cfTemp = 0;
saftTemp = 0;
for (idxRel = 1:nX)
% calculate distance between
distTot = sqrt(deltaZ * deltaZ + xRel(idxRel) * xRel(idxRel));
deltaT = round(distTot / SOS / dr(1));
if (zDepth > fd)
deltatIdx = idxFoc + deltaT;
else
deltatIdx = idxFoc – deltaT;
end
if (deltatIdx >= 1) && (deltatIdx <= nT)
saftTemp = saftTemp + sliceData(deltatIdx, idxRel);
cfTemp = cfTemp + abs(sliceData(deltatIdx, idxRel));
end
end
cf = saftTemp^2 / cfTemp^2;
sliceSAFTCoh(idxZ, idxX) = saftTemp * cf;
end
end
```
### 3D SAFT without coherence factor weighting
```MATLAB
% perform 3d saft without coherence factor weighting
sliceSAFT3 = zeros(size(vol));
yRel = (1:nY) - iY;
yRel = yRel * dr(3);
for idxZ = 1:nT % for each z position in the output volume
zDepth = vecZ(idxZ);
deltaZ = zDepth - fd;
for idxX = 1:nX % for each x position in the output volume
xRel = (1:nX) - idxX;
xRel = xRel * dr(2);
for (xi = 1:nX)
for (yi = 1:nY)
% calculate distance between
distTot = sqrt(deltaZ^2 + xRel(xi)^2 + yRel(yi)^2);
deltaT = round(distTot / SOS / dr(1));
if (zDepth > fd)
deltatIdx = idxFoc + deltaT;
else
deltatIdx = idxFoc – deltaT;
end
if (deltatIdx >= 1) && (deltatIdx <= nT)
sliceSAFT3(idxZ, idxX) = sliceSAFT3(idxZ, idxX) + vol(deltatIdx, xi, yi);
end
end
end
end
end
```
### 3D SAFT with coherence factor weighting
```MATLAB
sliceSAFT3Cf = zeros(size(vol));
yRel = (1:nY) - iY; % yRel is constant because we only recon a slice
yRel = yRel * dr(3);
for idxZ = 1:nT % for each z position in the output volume
zDepth = vecZ(idxZ);
deltaZ = zDepth - fd;
for idxX = 1:nX % for each x position in the output volume
xRel = (1:nX) - idxX;
xRel = xRel * dr(2);
cfTemp = 0;
saftTemp = 0;
for (xi = 1:nX) % for each relative x position
for (yi = 1:nY) % for each relative y position
% calculate distance between
distTot = sqrt(deltaZ^2 + xRel(xi)^2 + yRel(yi)^2);
deltaT = round(distTot / SOS / dr(1));
if (zDepth > fd)
deltatIdx = idxFoc + deltaT;
else
deltatIdx = idxFoc – deltaT;
end
if (deltatIdx >= 1) && (deltatIdx <= nT)
saftTemp = saftTemp + vol(deltatIdx, xi, yi);
cfTemp = cfTemp + abs(vol(deltatIdx, xi, yi));
end
end
end
cf = saftTemp^2 / cfTemp^2;
sliceSAFT3Cf(idxZ, idxX) = saftTemp * cf;
end
end
```
### GPU version of 3D SAFT with coherence factor weighting
```MATLAB
% perform 3D saft on GPU
% perform 3d saft with coherence factor weighting
tic
recon = simpleSaft(...
single(vol), ... % volume matrix sorted as [it, ix, iy]
single(dr), ... % resolution sorted as [dt, dx, dy]
single(fd), ... % focal distance of trasnducer
single(origin(1)), ... % zero timepoint
single(1495)); % speed of sound in medium
toc
```
mexcuda file:
```CUDA
#include