3D-FWI-MultiRowRingArrayUST

3D Full-Waveform Inversion for a Multi-Row Ring-Array UST System Simulated in k-Wave

Ultrasound tomography (UST) is a medical imaging system that uses the transmission of ultrasound through tissue to create images of the speed of sound. Currently, UST is based on a single-row elevation-focused ring-array transducer. This ring-array is translated vertically to reconstruct a series of 2D slices that are assembled into a 3D volume. The main algorithm used to reconstruct these 2D slices is 2D frequency-domain full-waveform inversion (FWI): see our previous work on 2D slice-wise FWI (rehmanali1994/WaveformInversionUST).

In this work, we extend the ring-array geometry to multiple rows of elements. We simulate a 22-cm diameter ring array with 32 rows of 256 elements placed circumferentially around the ring (2.4 mm between each row). By emitting cylindrical-wave transmits from this multi-row ring array, we can keep the number of transmisssions to a minimum while still capturing the full 3D insonification needed to reconstruct the full volume. With this imaging geometry, 2D FWI can be applied to each row of receive elements to reconstruct an image slice for each row. Here, we compare 3D FWI to 2D slicewise FWI in numerical breast phantoms.

We provide sample data and algorithms presented in the open-access IEEE OJ-UFFC letters paper with BibTeX reference below

@article{ali20253d,
  title={3D Frequency-Domain Full Waveform Inversion for Whole-Breast Imaging with a Multi-Row Ring Array},
  author={Ali, Rehman and Jin, Gaofei and Singh, Melanie and Mitcham, Trevor and Duric, Nebojsa},
  journal={IEEE Open Journal of Ultrasonics, Ferroelectrics, and Frequency Control},
  year={2025},
  publisher={IEEE}
}

You can also reference a static version of this code by its DOI number:

In this work, our FWI algorithm relies on the convergent Born series method to efficiently solve the 2D and 3D Helmholtz equations for 2D slicewise and 3D FWI, respectively. In addition to the reference above, please cite the following paper on the convergent Born series method from which our solveHelmholtzBornSeries.m and solveHelmholtzBornSeries3D.m functions were adapted:

@article{osnabrugge2016convergent,
  title={A convergent Born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media},
  author={Osnabrugge, Gerwin and Leedumrongwatthanakun, Saroch and Vellekoop, Ivo M},
  journal={Journal of computational physics},
  volume={322},
  pages={113--124},
  year={2016},
  publisher={Elsevier}
}

Running the k-Wave Simulations

The code used to run the k-Wave simulation involves a 2-step process:

1. GenKWaveSimInfoFromMRI.m creates a MAT file that is stored in the sim_info folder. This MAT file (either LeftBreastMRI.mat or RightBreastMRI.mat) contains all the information (medium, multi-row ring-array geometry, and pulse excitation) needed to simulate the UST system. GenKWaveSimInfoFromMRI.m uses sampled_cylinder.m to create the multi-row ring-array transducer and BreastPhantomVolumetricMRI.m to produce the material properties used in the simulation. BreastPhantomVolumetricMRI.m relies on the BreastPhantomFromMRI.mat found under the releases tab. Please place BreastPhantomFromMRI.mat into the phantoms folder. Note that GenKWaveSimInfoFromMRI.m can run two different simulations: one corresponding to the left breast (leftBreast = true) and the other corresponding to the right breast (leftBreast = false). This option is passed into BreastPhantomVolumetricMRI.m function where it selects one of the breasts for the phantom.
2. After generating the MAT file in sim_info, we run the actual k-Wave simulation using kWaveSimLauncher.m, which loops through each single-element transmit (function calls to GenRFDataSingleTxKWave.m). The simulated data for each transmit is then stored in MAT files in the sim_data folder. To conserve space, the traces simulated for each transmit are Fourier transformed, and the scripts in the visualization folder can be used to visualize the time-domain data for each transmit. At the end of kWaveSimLauncher.m, the function call to AssembleSimData.m assembles the simulated data from each indivdual transmit/MAT-file in the sim_data folder into a single MAT file (either LeftBreastMRI.mat or RightBreastMRI.mat) containing the 3D UST dataset needed for FWI (stored in the datasets folder). Note that the complete time-series data is too large to store in a single file for the full 3D dataset. Instead, we extract only the frequencies listed in kWaveSimLauncher.m to keep the file size for the full 3D dataset small.

FWI Reconstruction Code

Two different scripts are used to run FWI (either 2D slicewise or 3D FWI) and their results are compared (see Compare2Dvs3D.m):

1. MultiFrequencyWaveformInvFromData.m is used to run 2D slicewise FWI and outputs the result into results/results2D.
2. MultiFrequencyWaveformInvFromData3D.m is used to run 3D FWI and outputs the result into results/results3D.

The key functions used in the FWI scripts above are:

1. solveHelmholtzBornSeries.m - uses the convergent Born series method to solve the 2D Helmholtz equation for multiple sources. This function is used in MultiFrequencyWaveformInvFromData.m to run 2D slicewise FWI. In the presence of a GPU, this code will use MATLAB's built gpuArray type to accelerate the FFTs involved.
2. ringingRemovalFilt.m - Helper function called during the 2D slicewise FWI script (MultiFrequencyWaveformInvFromData.m) to remove ringing artifacts in the images.
3. solveHelmholtzBornSeries3D.m - uses the convergent Born series method to solve the 3D Helmholtz equation for multiple sources. In the presence of a GPU, this code will use MATLAB's built gpuArray type to accelerate the FFTs involved.
4. solveHelmholtzBornSeries3D_Batches.m - Note that the GPU may not be able to compute the solution for all the sources at the same time, so we break up the sources into multiple batches: solveHelmholtzBornSeries3D_Batches.m has an extra input argument for the number of batches. Inside solveHelmholtzBornSeries3D_Batches.m, we loop solveHelmholtzBornSeries3D.m over the batches to solve the 3D Helmholtz equation for all the sources. This solveHelmholtzBornSeries3D_Batches.m function is used in MultiFrequencyWaveformInvFromData3D.m to run 3D FWI. See the NUM_BATCHES variable inside the MultiFrequencyWaveformInvFromData3D.m. Because each batch is computed serially, using fewest possible batches should compute the solutions in the least amount of time. However, if using too few batches, the number of sources in each batch may be too large for the GPU to actually compute the solution (because FFTs need space). Please adjust NUM_BATCHES as needed. Note: NUM_BATCHES should be a power of 2.
5. ringingRemovalFilt3D.m - Helper function called during the 3D FWI script (MultiFrequencyWaveformInvFromData3D.m) to remove ringing artifacts in the images.
6. subplus.m - Helper function for PML in 2D and 3D Helmholtz equation. This function already comes with MATLAB's Curve Fitting Toolbox; however, we include this here in case you lack access to the Curve Fitting Toolbox.

We recommend running this code with a beefy GPU that has lots of GPU RAM like the NVIDIA RTX A6000 (48 GB GDDR6) or the NVIDIA GeForce RTX 4090 (24 GB GDDR6) so that you can keep the NUM_BATCHES in MultiFrequencyWaveformInvFromData3D.m as small as possible. Even more importantly, your computer will need lots of CPU RAM (minimum 320 GB) just to store intermediate variables for 3D FWI.

Sample Results

For either the LeftBreastMRI.mat or RightBreastMRI.mat (fileID = 1 or fileID = 2 respectively), Compare2Dvs3D.m loads the results of 2D slicewise and 3D FWI from the results/results2D and results/results3D folders. Compare2Dvs3D.m will then display the ground truth sound speed (left), 2D slicewise FWI (middle), and 3D FWI (right) images for a given orthographic slice view: transverse (viewID = 1), sagittal (viewID = 2), and coronal (viewID = 3). Compare2Dvs3D.m also calculates the root-mean-square error (RMSE) and Pearson's correlation coefficient (PCC) with the ground truth sound speed image for both the 2D slicewise and 3D FWI reconstructions.

See the results below compiled across the k-Wave simulations of the left and right breast: