In regions of multiple scattering, the direction of photon propagation becomes stochastic, resulting in a loss of imaging capability. Consequently, the majority of optical techniques used for imaging deep within tissues rely on single-scattering photon imaging^[1–4]. To overcome the challenges of multiple scattering, optical coherence tomography (OCT) techniques employing coherent gating have been developed to isolate and extract single-scattering photons more effectively^[5–9]. Compared to traditional microscopy, OCT significantly enhances imaging depth. However, as the depth increases, the number of single-scattering photons decreases exponentially, limiting its imaging capability to within a few transport mean free paths in scattering tissues^[10,11]. The recently introduced Smart-OCT has enhanced the imaging depth of OCT using a reflection matrix (RM)^[12]. Current methods for measuring the RM in scattering media involve techniques such as phase-shifting^[13–15] or heterodyne detection^[16–18]. However, these methods require detectors to measure all sampling points to ensure accuracy and high resolution, which increases the burden of transmission and processing, thereby limiting measurement speed. Furthermore, the theoretical basis for the RM in scattering media assumes that photon propagation in static media is a linear and time-invariant process^[19–21]. As the measurement time increases, imaging systems become more susceptible to environmental changes^[22]. This is especially true in dynamic environments like blood flow and heartbeat, where the properties of imaging targets or scattering media can change rapidly^[17]. This leads to image capture delays and an inability to accurately reflect instantaneous states. In large-scale imaging, the system must balance high resolution with wide coverage while also capturing fine details, such as detecting small lesions in medical imaging. Therefore, reducing RM measurement time is crucial for improving imaging efficiency and accuracy.

Compressed sensing (CS) technology is a powerful tool for image reconstruction. By leveraging the sparsity of signals, it accurately reconstructs data from minimal random sampling, thereby reducing the volume of data required for collection, transmission, and storage. Initially proposed as a mathematical theory^[23,24], CS has been applied to medical imaging fields such as magnetic resonance imaging and photoacoustic imaging^[25,26]. Subsequently, CS was introduced into the field of OCT^[27] and experimentally validated for image data reconstruction^[28]. This includes improvements in CS methods and sparse camera sampling approaches^[29,30]. To accelerate the CS reconstruction process, computer systems equipped with graphics processing units (GPUs) have been utilized, enabling real-time, rapid CS reconstruction^[31,32].

In this paper, a novel technique is introduced for measuring the RM of scattering media using CS and applied to reflection matrix optical coherence tomography (RM-OCT) imaging. The sparsity of RM in scattering media was first validated through an electric field Monte Carlo model, confirming the feasibility of using CS for measurement. Experimental measurements and imaging of the RM in scattering media were then conducted. By exploiting the RM’s sparsity, the matrix can be reconstructed from a subset of sampling points. Results indicate that this method requires less than 50% of the data used in traditional sampling methods to recover imaging target information in scattering media, significantly reducing data acquisition time. This approach has substantial potential to enhance imaging efficiency in scattering media.

CS is a technique based on the sparsity of signals^[23,24]. By utilizing signal sparsity in specific domains, CS can reconstruct the original signal with fewer measurements than the Nyquist limit. In scattering media, multiple scattering leads to highly random photon trajectories, resulting in a relatively small proportion of single scattering components within the RM^{[12,14,33,34]}. However, single-scattering events are closely related to specific light paths and the local structure of the medium. Theoretically, information related to single scattering is primarily concentrated along the diagonal and nearby regions of the RM, while the effects of multiple scattering are represented by the off-diagonal elements^{[14,18,35,36]}. This characteristic suggests that the RM has a sparse structure, providing a theoretical basis for applying CS measurement techniques to the RM.

Figure 1(a) illustrates the schematic of the RM-OCT system. The system uses a laser centered at 790 nm as its light source. A spatial light modulator (SLM) performs biaxial scanning by projecting a set of phase patterns onto the area of interest. Subsequently, a CCD is used to collect backscattered photons from the illuminated area on the imaging plane. The detailed description of the RM-OCT experimental setup is shown in Refs. [16–18].

Figure 1(b) illustrates the principle of RM reconstruction based on CS. The traditional method for constructing an RM first obtains the light field information at each focal point on the imaging location through uniform scanning. Next, the acquired light field information is subjected to a Fourier transform, and finally, it is reordered according to the scanning sequence to generate a two-dimensional reflectance matrix. To reduce sampling time, a reflectance matrix measurement technique is proposed based on CS. A small amount of light field information obtained through random sampling is used to reconstruct the complete reflectance matrix. The measurement process for the reflectance matrix based on CS can be modeled as^[37]Y=A×RM,(1)where Y is an (m×n, m<n) observation matrix obtained from sampling, A is an (m×n) measurement matrix, and RM is an (n×n) reflectance matrix to be solved.

Based on prior information, the RM primarily reflects the information related to single scattering events along the diagonal and its adjacent regions, while the elements away from the diagonal are mainly associated with multiple scattering effects^[35]. The following optimization problem is formulated: minRM∥Y−A×RM∥F2+λ∑i,jWi,j|RMi,j|.(2)

Here, the first term ∥Y−A×RM∥F2 represents the data fidelity term, ensuring the accuracy of the reconstructed RM, while the second term is a position-weighted L1 regularization term, promoting sparsity in the matrix RM, particularly for elements away from the diagonal. The notation ∥·∥F represents the Frobenius norm. The Frobenius norm is a matrix norm used to measure the square root of the sum of the squares of the matrix elements. The element RMi,j represents the value at the ith row and jth column of the RM, where i,j=1,2,…,n. The summation ∑i,j applies to all elements of the matrix, ensuring that the regularization term covers all possible index pairs. The position weights Wi,j are designed as Wi,j=exp(α|i−j|),(3)where α is a parameter controlling the rate of decay, ensuring that elements near the diagonal receive less penalization, while elements further from the diagonal are penalized more heavily. This careful weighting helps balance the risk of overfitting by emphasizing the importance of the significant elements close to the diagonal while discouraging the influence of less relevant elements further away, thus enhancing the robustness of the reconstruction.

To address this sparse constraint problem, we employ the alternating direction method of multipliers (ADMM). By introducing an auxiliary variable Z, the optimization problem can be reformulated as ^[38–40]minRM,Z∥Y−A×RM∥F2+λ∑i,jWi,j|Zi,j|,(4)subject toZ=RM.(5)

The ADMM-solving process consists of the following steps:

1. 1)Updating RM: RMk+1=argminRM∥Y−A×RM∥F2+ρ2∥RM−Zk+Uk∥F2.(6)This step is a quadratic optimization problem and can be solved using gradient descent methods^[41,42].
2. 2)Updating Z: Zk+1=argminZλ∑i,jWi,j|Zi,j|+ρ2∥RMk+1−Z+Uk∥F2.(7)This step is solved via the soft-thresholding operation:^[43,44]Zi,jk+1=sign(RMi,jk+1+Ui,jk)·max(|RMi,jk+1+Ui,jk|−λWi,jρ,0).(8)The soft-thresholding operation encourages sparsity in Z, particularly in elements away from the diagonal.
3. 3)Updating the Lagrange multiplier U:^[45,46]Uk+1=Uk+(RMk+1−Zk+1).(9)This step ensures that the constraint Z=RM is gradually enforced during the iterative process.

The above steps are iterated until convergence.

After reconstructing the RM, it is subjected to singular value decomposition (SVD), as shown in Fig. 2(a). By performing SVD on the RM, the statistical characteristics of the scattering pattern can be obtained. The SVD of the RM can be expressed as^[16]RM=UΣV*(10)where U and V are the orthogonal matrices, V* represents the conjugate transpose of matrix V, and Σ is a diagonal matrix consisting of the real positive singular values in decreasing order, σ1>σ2>…>σn, where n is the size of Σ. Larger singular values correspond to single-scattering photons, while smaller singular values correspond to multiple-scattering photons. By selecting the appropriate singular vectors for reconstruction^[12], it is possible to separately reconstruct images of the single-scattering light field (Iss) and the multiple-scattering light field (Ims), as shown in Fig. 2(b)^[12,14], Iss=∑i=1tσi|U¯i°V¯i|,(11)Ims=∑i=tnσi|U¯i°V¯i|,(12)where U¯i and V¯i are the wavefronts reshaped from the vectors Ui and Vi, respectively, and the reshaping order matches the scanning trajectory of the light field sampling window. The symbol ° denotes the Hadamard product. t is the threshold for distinguishing between single- and multiple-scattering light fields.

To evaluate the feasibility of using CS for RM measurement in complex scattering environments, we first conducted a simulation experiment using the electric field Monte Carlo method, which was developed by Xu and his colleagues^[47–49]. Unlike traditional intensity-based Monte Carlo simulations, the electromagnetic compatibility (EMC) method tracks the changes in the electric field during each scattering event, preserving phase information. These changes are analyzed based on Mie scattering theory, allowing for a more comprehensive analysis of the interactions of light within complex media. For simplicity, a “sandwich” structure sample was designed. Under linearly polarized coherent illumination with a central wavelength of λ=790nm, the top layer contains scatterers with absorption characterized by a mean free path lt=1.05μm, a refractive index of n=1.60+0.25i, an anisotropic factor of g=0.210, and a thickness of 3lt. The middle layer contains two different types of scatterers: (1) those identical to the top layer and (2) those with reduced absorption arranged in a USAF target pattern. The scatterers in the USAF target pattern have a mean free path of lt′=19.00μm, a refractive index of n′=1.60+0.002i, an anisotropic factor of g′=0.050, and a thickness of 2lt. The bottom layer has the same parameters as the top layer. To compare the effects of different sampling methods on RM measurement, both uniform and random sampling were performed on the sample, as shown in Fig. 3(a). The USAF target in the middle layer is illustrated in Fig. 3(b), and the image obtained under coherent illumination is presented in Fig. 3(c). Figures 3(d) and 3(e), respectively, show the RM obtained from uniform sampling (with 51×51 sampling points) and 50% random sampling. The RM was then subjected to SVD, with the singular value distribution displayed in Fig. 3(f). The top 20 singular values were selected for the light field image reconstruction, as shown in Figs. 3(g) and 3(h), which present the reconstructed light field images from uniform and sparse sampling, respectively. Figure 3(i) shows the intensity distribution along the white dashed lines in Figs. 3(g) and 3(h), highlighting the details of the reconstruction quality. To further evaluate the impact of different sampling rates on RM measurement, the mean squared error (MSE)^[50], structural similarity index (SSIM)^[51], and correlation coefficient (CC)^[52] were used for quantitative analysis. The MSE measures the difference between the reconstructed and reference images (uniform sampling), the SSIM evaluates the similarity in brightness, contrast, and structure, and the CC measures the linear correlation between the two images. From the quantitative analysis of Figs. 3(g) and 3(h), the MSE, SSIM, and CC values are 0.0038, 0.9528, and 0.975, respectively, indicating that the RM reconstructed using CS retains consistent information with uniform sampling. This supports the feasibility of using CS for RM measurement in scattering media.

To further validate the effectiveness of CS-based RM measurement in scattering media, experiments were designed and conducted at different sampling rates. In the experiment, a silicone rubber phantom material with a thickness of approximately 2 mm was used, with 10 g of titanium dioxide added to every 100 g of silicone rubber to enhance its scattering characteristics. The light beam was focused at a depth of approximately 0.7 mm within the sample to ensure a more accurate measurement of scattering properties deep within the sample. According to a predefined scanning strategy, both uniform sampling and random sampling were performed on the sample to compare the effects of different sampling methods on RM measurement. Uniform sampling utilized a total of 961 scanning points (31×31), covering an area of approximately 70μm×70μm in the region of interest. The acquisition time for each scanning point was 80 ms. Figures 4(a)–4(c) show the RM distributions under 30% and 50% random sampling rates, as well as uniform sampling. Figure 4(d) presents the singular value distribution of the RM under different sampling rates. Although there are some differences in the singular value distribution between lower sampling rates (e.g.,  50%) and uniform sampling, the major high singular values remain relatively consistent. This indicates that, even at lower sampling rates, the main features of the signal can still be effectively recovered, as shown by the light field distributions corresponding to the dashed lines in Figs. 4(a) and 4(b). Figure 4(e) illustrates the intensity distribution along the blue dashed lines in Figs. 4(a)–4(c). Figure 4(f) shows the deviation in intensity between the 30% and 50% sampling rates and uniform sampling at the blue dashed line positions. Through a quantitative analysis of Figs. 4(a)–4(c), the MSE, SSIM, and CC for the 30% random sampling rate compared to uniform sampling are 0.0024, 0.635, and 0.6985, respectively. For the 50% random sampling rate, these values are 0.0008, 0.8549, and 0.9153, respectively. Quantitative experimental results indicate that CS can effectively recover RM features with reduced sampling points. These results validate the effectiveness of the CS-based RM measurement method in scattering media.

To verify the imaging effectiveness of CS RM-OCT in scattering media, a housefly wing was selected as the imaging target. In the experiment, scattering was introduced by covering the surface of the housefly wing with a layer of subcutaneous chicken mucosa approximately 20μm thick, as shown in Fig. 5(a). Figure 5(b) illustrates the intensity distribution of the light field illuminating the sample directly, where the target information is obscured by the scattered light field. The intensity distribution along the black dashed line is depicted in Fig. 5(c). To mitigate the effects of multiple scattered photons, imaging based on the RM was performed. Figures 5(d) and 5(e) present the reflection matrices obtained using 50% compressed sampling and uniform sampling, respectively. The singular value distribution, derived through the SVD of the RM, is shown in Fig. 5(f). Figures 5(g) and 5(h) display the image reconstruction results using the top 1000 singular values from the reflection matrices for compressed sampling and uniform sampling, respectively. By reconstructing the light field image from the RM, the influence of multiple scattered photons is reduced, resulting in an image contrast that is doubled compared to the light field in Fig. 5(b). Figure 5(i) shows the intensity distribution along the white dashed line, highlighting the reconstruction quality. A quantitative analysis of Figs. 5(g) and 5(h) reveals that the MSE, SSIM, and CC for the 50% random sampling rate compared to uniform sampling are 0.0074, 0.807, and 0.9063, respectively. These results indicate that even with a 50% sampling rate, the imaging target’s information can be accurately recovered. This confirms the feasibility of the CS-based RM imaging method.

This study presents a method that utilizes CS technology to accelerate the measurement of the RM while reducing the required number of sampling points. Traditional high sampling rates provide more information, improving reconstruction accuracy, but this also increases data acquisition time and computational burden. In contrast, lower sampling rates can reduce acquisition time and processing requirements but may lead to degraded image quality due to insufficient data. In our experiments, a 50% sampling rate achieved a good balance between image quality and computational efficiency, while a 30% sampling rate offered faster processing but resulted in a decline in image quality. Lower sampling rates may lead to incomplete reconstruction of the RM, as they fail to capture sufficient information to accurately restore the original light field, causing loss of image details and increased noise. To mitigate these issues, noise filtering techniques, such as Gaussian aperture filtering, can be applied during data processing to reduce the impact of noise. Additionally, during the RM reconstruction process, we can utilize prior information about the RM to constrain the procedure, thereby reducing overfitting and computational errors. Although lower sampling rates may introduce errors and decrease imaging quality, these errors can be significantly minimized through effective strategies and methods, ultimately providing reliable high-quality imaging results.

Although significant progress has been made in applying CS to RM-OCT, further exploration is needed regarding its use in complex and dynamic scenarios where rapid data acquisition is crucial for effective information retrieval. According to the scanning strategy (see Fig. 4), uniform sampling for the RM requires 961 scanning points (31×31), with each point taking approximately 80 ms to acquire, resulting in a total sampling time of about 76 s. In contrast, using a 50% sampling rate reduces the number of sampling points to approximately 480, shortening the sampling time to about 38 s. Despite the complex computational operations involved and the recovery process being proportional to the data volume, CS significantly enhances sampling efficiency. By leveraging the sparsity of the RM and applying CS techniques to reduce computational complexity, the CS program currently operates on a 2.80 GHz Intel Core i9 10900 CPU, with the recovery process taking about 100 s. The total time for RM measurement based on CS is approximately 138 s, including both data acquisition and recovery. While the overall processing time with CS is longer, it improves sampling efficiency by reducing the number of sampling points. Optimizing processing and utilizing sparsity not only reduces the required sampling points but also enhances imaging efficiency. Compared to central processing units (CPUs), GPUs offer significant advantages in matrix operations. Therefore, future work will focus on implementing GPU-based CS to further reduce processing time and enhance imaging speed, ultimately optimizing the overall imaging process.

In summary, we propose a method for measuring the RM of scattering media using CS technology to reduce measurement time. Experimental results demonstrate that the RM of scattering media is sparse, enabling effective CS sampling. This approach reduces the number of sampling points to less than half compared to traditional uniform sampling methods. CS technology provides a more efficient method for RM measurement, opening new possibilities for the application of CS-related techniques and showing potential for enhancing imaging efficiency in scattering media.