Simple MATLAB implementations of SAFT

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 // used to interact with matlab #include // used to run code on GPU #include #include #include // used to export variables into files to lateron debug them __global__ void SAFT ( float* outputVol, // reconstructed output volume [it, ix, iy] const float* inputVol, // input signal matrix [it, ix, iy] const int nT, // number of elemetns in time const int nX, // number of elements along x const int nY, // number of elements along y const float dt, // resolution in time domain const float dx, // resolution along x [m] const float dy, // resolution along y [m] const float t0, // time of first index const float c0, // speed of sound in medium [m/s] const float fd // focal distance of transducer [m] ) { // get index in output volume const int zIm = blockIdx.x * blockDim.x + threadIdx.x; const int xIm = blockIdx.y * blockDim.y + threadIdx.y; const int yIm = blockIdx.z * blockDim.z + threadIdx.z; if ((xIm < nX) && (yIm < nY) && (zIm < nT)) { const float zDepth = (t0 + dt * (float) zIm) * c0; const float deltaZ = zDepth - fd; float rfsaft = 0; // coherent saft sum float rfsaftabs = 0; // incoherent saft sum const int idxFoc = round((fd / c0 - t0) / dt); for (int iY = 0; iY < nY; iY++) // run over all other a scans in y { const float yRel = dy * (float) (iY - yIm); for (int iX = 0; iX < nX; iX++) // run over all other a scans in x { const float xRel = dx * (float) (iX - xIm); // calculate distance, delay time, and index const float distTot = sqrt(xRel * xRel + yRel * yRel + deltaZ * deltaZ); const int deltaT = round(distTot / c0 / dt); const int deltatIdx = (zDepth > fd) ? idxFoc + deltaT : idxFoc – deltaT; if ((deltatIdx >= 0) && (deltatIdx < nT)) { const int dataIdx = deltatIdx + nT * (iX + nX * iY); rfsaft += inputVol[dataIdx]; rfsaftabs += abs(inputVol[dataIdx]); } } } if (rfsaftabs > 0) { const float cf = rfsaft * rfsaft / rfsaftabs / rfsaftabs; outputVol[zIm + nT * (xIm + nX * yIm)] = rfsaft * cf; } else { outputVol[zIm + nT * (xIm + nX * yIm)] = 0; } } return; } // input arguments // [0] input signal matrix sorted as [t, x, y] // [1] resolution vector sorted as [dt, dx, dy] // [2] focal distance of transducer // [3] t0 time offset // [4] speed of sound void mexFunction( int nlhs, mxArray *plhs[], // variables returned to MATLAB int nrhs, const mxArray *prhs[]){ // variables passed from MATLAB to C // get input data and dimension const float* inputVol = (float *) mxGetData(prhs[0]); const int nt = (int) mxGetDimensions(prhs[0])[0]; const int nx = (int) mxGetDimensions(prhs[0])[1]; const int ny = (int) mxGetDimensions(prhs[0])[2]; const int nElements = nt * nx * ny; // get resolution of dataset const float* res = (float*) mxGetData(prhs[1]); // get focal distance of transducer const float fd = (float) mxGetScalar(prhs[2]); // get time offset const float t0 = (float) mxGetScalar(prhs[3]); // get speed of sound const float c0 = (float) mxGetScalar(prhs[4]); // allocate memory on GPU for input and output volume float* inputVol_dev; float* outputVol_dev; bool m1 = (cudaSuccess != cudaMalloc( (void**)&inputVol_dev, nElements * sizeof(float) )); m1 |= (cudaSuccess != cudaMalloc( (void**)&outputVol_dev, nElements * sizeof(float) )); if (m1) { mexErrMsgTxt(“Could not allocate memory on GPU”); return; } // copy signal matrix over to GPU and check if successful bool cpy1 = (cudaSuccess != cudaMemcpy(inputVol_dev, inputVol, nElements * sizeof(float), cudaMemcpyHostToDevice)); if (cpy1) { mexErrMsgTxt(“Could not copy array to GPU”); return; } // define kernel size dim3 blockSize(4, 4, 4); dim3 gridSize( (nt + blockSize.x – 1) / blockSize.x, (nx + blockSize.y – 1) / blockSize.y, (ny + blockSize.z – 1) / blockSize.z); // execute actual kernel SAFT<<< gridSize, blockSize >>> ( outputVol_dev, inputVol_dev, nt, nx, ny, res[0], res[1], res[2], t0, c0, fd ); cudaDeviceSynchronize(); // check if kernel execution was successful cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { printf(“error during kernel execution:\n”); mexErrMsgTxt(cudaGetErrorString(err)); } // allocate memory for matlab and copy data backwards const mwSize dims[3] = {nt, nx, ny}; mxArray* volArray = mxCreateNumericArray(3, dims, mxSINGLE_CLASS, mxREAL); plhs[0] = volArray; float* outData = (float*) mxGetData(volArray); cudaMemcpy(outData, outputVol_dev, nElements * sizeof(float), cudaMemcpyDeviceToHost ); // free gpu memory again cudaFree(outputVol_dev); cudaFree(inputVol_dev); return; } “` ## Conclusion SAFT is (sadly) still the standard tool to reconstruct images in acoustic resolution optoacoustic microscopy. With the code above I show how simple its implementation is both for CPU and GPU through CUDA. One of the major weaknesses and unsolved problems is the focal plane in which the transducer shows its highest sensitivity but the reconstruction through SAFT fails here. Promising alternative approaches include in my opinion model based reconstruction and the hybrid domain approach.

Leave a Reply

Your email address will not be published. Required fields are marked *