This repository contains the code for performing large-scale full-wave simulations of several scattering-based imaging methods:
"Numerical experiments of tomographic optical imaging inside scattering media"
- Scattering matrix tomography (SMT) in reflection mode
- Reflectance confocal microscopy (RCM)
- Optical coherence tomography (OCT)
- Optical coherence microscopy (OCM)
- Interferometric synthetic aperture microscopy (ISAM)
"Deep imaging inside scattering media through virtual spatiotemporal wavefront shaping"
- Synthetic OCT
- Synthetic OCM
We use the scattering matrix to perform numerical modeling. The scattering matrix encapsulates the sample’s complete linear response. Given the scattering matrix of the sample, we can obtain the scattering field from arbitrary input and model any scattering-based imaging method. More importantly, the computation of scattering matrices is accelerated by several orders-of-magnitude thanks to the introduction of a new technique called augmented partial factorizaion, making it possible to perform the full-wave simulation on a large system. Such numerical modeling can be helpful for developing new imaging methods by providing the ground truth, the flexibility to tailor the system and the imaging scheme, and the ease of comparing different methods.
The code is written in MATLAB so no compilation is required. But you need the following dependencies to run different components of the code:
- System setup: MESTI, METIS (optional)
- Reflection matrix computation: MESTI, MUMPS (optional but highly recommended)
- Image reconstruction: MESTI, FINUFFT (optional)
- Field profile computattion: MESTI, MUMPS
Note that if you plan to use one component of the code, there is no need to install dependencies of other components.
MESTI is required for all components of the code. The installation of MESTI is straightforward. Simply download the source code from here and add it to the MATLAB search path.
METIS produces high-quality orderings for matrix factorization, which accelerates the reflection matrix computation. In our large system, the computation time is shorten by about 20%.
Download the source code from here. No need to use the OpenMP version. You can compile METIS using the default option:
$ make config
$ make
265s
Once the compilation is finished, you can find the METIS library under metis-x.x.x/build/*/libmetis. The library path and header file path (metis-x.x.x/include/) will be used later for the compilation of MUMPS. There is no need to add METIS to the MATLAB search path.
In the reflection matrix computation, we need MUMPS to perform the partial factorization in the APF method. It is possible to run the code without MUMPS but it will use another method and the computation time will be much longer.
The detailed instructions for compiling MUMPS can be found here.
Most of our image reconstruction code uses non-uniform fast Fourier transforms (NUFFTs). FINUFFT is a library to perform the NUFFTs efficiently. MATLAB also implements its own NUFFT function in R2020a and improve the performance in R2022a and R2023b. The code will use the MATLAB NUFFT if FINUFFT is not installed. We found that on R2023b the reconstruction time is comparable with FINUFFT or MATLAB NUFFT. Thus, you can simply use the MATLAB NUFFT if your MATLAB version is R2023b or later.
To install FINUFFT, download the source code from here. There are two routes to compile: the new CMake-based route and the old GNU makefile-based route. Note that the CMake-based route is still under development, which does not support many building options. You can take the makefile-based route for now.
There is a caveat in compiling the MATLAB interface on Mac. You may get a warning of license has not been accepted from Xcode and a following error of no supported compiler was found. The same error can happen to the MATLAB interface of MUMPS. The simple solution is typing
$ /usr/libexec/PlistBuddy -c 'Add :IDEXcodeVersionForAgreedToGMLicense string 10.0' ~/Library/Preferences/com.apple.dt.Xcode.plist
in the command line.
Remember to add finufft-x.x.x/matlab to the MATLAB search path after the compilation.
After installing the required dependencies above, you can simply download the code and add all folders to the MATLAB search path. Make sure that the working directory is /path/to/Imaging-simulations when you run the code.
To get started, we suggest running the image reconstruction code for the system described in our paper. After downloading the precomputed reflection matrix, system_data.mat file and moving them to Image-simulations/data/large_system, you can use recon_all.m to reconstruct all images. Before running the script, make sure that the working directory is /path/to/Imaging-simulations and all scripts are in the MATLAB search path. It takes less than one minute to reconstruct each image on a MacBook Air with Apple M1 chip.
After obtaining the reconstructed images, you may download the z-dependent weights and use the scripts/plotting/plot_all_images.m to plot all images.
In addition, we upload a small system example where all computation can be done on a local machine. The total computation time for the reflection matrix is about 15 mins on a MacBook Air with Apple M1 chip. The image reconstruction takes a few seconds.
Here, we give an overview of the functionality and usage of all components. For more detailed information, please refer to comments at the beginning of each script.
This component (Imaging-simulations/scripts/system_setup) builds the system and outputs a system_data.mat file, which contains all the necessary data for reflection matrix computation and image reconstruction. In addition, it also outputs METIS ordering files if METIS is installed and the option of producing METIS orderings is set to true.
To build the system, you should first create a directory
data_dir and a script data_dir/set_simulation_parameters.m. All the simulation parameters shall be defined in that script. You can use our sample systems as the template. The data_dir will be used later for saving all simulation data of this system, including the reflection matrices and reconstructed images. In this repository, we use data_dir=Imaging-simulations/data/[system_name]. After defining the parameters, you can specify your data_dir in build_system.m and run the script. The script builds the system and saves all relevant data in data_dir/system_data.mat. It also saves optional METIS ordering files under data_dir/orderings.
This component (Imaging-simulations/reflection_matrix_computation) computes hyperspectral angular/spatial reflection matrices. There are three scripts in this componenet. The compute_angular_R.m computes angular reflection matrices at multiple frequencies while compute_spatial_R.m computes a subset of the spatial reflection metrics at multiple frequencies. The planckwin.m is used for generating a Planck-taper window profile for smoothing the line source profile.
The R computation function loads necessary data from system_data.mat and the METIS orderings if exist. After the computation is done, the reflection matrix will be saved under data_dir/hyperspectral_reflection_matrices. Note that reflection matrices are converted to single-precision for reducing the storage space and accelerating the image reconstrcution speed.
This component (Imaging-simulations/image_reconstruction) consists of functions for reconstructing images from hyperspectral reflection matrices. Each function implements one or two imaging methods menitioned in the introduction section.
To use those functions, you need to specify the data directory data_dir and/or the imaging method in the arguments. The function will load necessay data from data_dir/system_data.mat and reflection matrices from data_dir/hyperspectral_reflection_matrices, reconstruct the image, and save the reconstructed image data under data_dir/reconstructed_images. Note that the image data is saved in single-precision, which is sufficient for most imaging applications.
This component (Imaging-simulations/depth_dependent_weights) generates depth-dependent weights for reducing the intensity variation across different image depths. The weight is an exponential function of a piecewise polynomial. The polynomial coefficients are determined by manual tuning to fit the maximum image intensity at each depth.
The weight parameters for all images are defined in data_dir/z_z_dependent_weights/get_weight_parameters. You specify the imaging method in the script. The script will generate the weight and plot the weight and the maximum image intensity at each depth. As a rule of thumb, you should tune the weight parameters such that trend of both curves matches. You can specify the option of saveing the weight to true after the tunning. The weight will be saved under data_dir/z_dependent_weights.
This component (Imaging-simulations/field_profile_computation) contains the script for computing field profiles with the plane wave or focused Gaussiab beam input.
This component (Imaging-simulations/plotting) consists of some plotting scripts. the ground truth, reconstructed images and field profiles.
You should specify the data directory data_dir and the list of imaging methods. If you have generated depth-dependent weights for those methods, you can set the option of applying weights to true. The script will load the image data from data_dir/reconstructed_images, divide the image intensity by the depth-dependent weight if applicable, and plot the images. It also plots the ground truth if you includes "ground_truth" in the list of imaging methods. By setting the option of saving figures to true, you can save the JPEG image under Imaging-simulations/figs/[system_name].
In this section, we show some image reconstruction results.
