### Monte Carlo Simulations  
Monte Carlo (MC) offers a flexible and accurate approach for modeling light propagation within tissues [18]. In order to investigate the effects of different optical properties on reflectance prediction from the two-layered diffusion model, to evaluate the estimation accuracy of different curve fitting methods, and to analyze the factors influencing optical property extraction of the top layer, a publicly available MC simulation program for multi-layered turbid media was used [19]. In the simulations, a package of five million photons was tracked. The maximum radial distance of the medium was set to 50 mm, which is large enough to be treated as semi-infinite. The spatial resolution for both radial distance and depth was set to 0.1 mm. The average refractive indices of the two-layered media were both chosen to be 1.35, which was typical for most horticultural products (e.g., apple, blueberry, citrus, tomato etc.) [20; 21], while the media above and beneath the tissue were treated as air, with the refractive index of 1.0. The anisotropy factor \(g\) was set to 0.9, and Henyey-Greenstein phase function was used to describe light scattering. Figure 1: Schematic of a two-layered turbid medium under structured illumination. \(\mu_{a1}\), \(\mu_{a2}\), \(\mu_{s1}^{\prime}\) and \(\mu_{s2}^{\prime}\) are absorption and reduced scattering coefficients for the top layer and bottom layer, respectively, \(d\) is the thickness of the top layer, and \(R_{d}\) is the diffuse reflectance at the surface. More details about the phase function can be found in reference by [PERSON] and [PERSON] (1941) [22]. The Monte Carlo multi-layered (MCML) program was first applied to generate the spatially-resolved diffuse reflectance profiles of two-layered samples along the radial distance under the normal incidence of an infinite small light source. Then the 1-D Hankel transform of order zero was used to convert the spatially-resolved reflectance to spatial-frequency domain reflectance [16], which was further used for evaluating the curve fitting methods and investigating the potential influencing factors through optical property estimations. In this study, a total of 20 combinations of \(\mu_{a}\) and \(\mu_{s}{{}^{\prime}}\) for two-layered samples (five samples with varying \(\mu_{a}\) or \(\mu_{s}{{}^{\prime}}\) of each layer with the top-layer thickness of 2 mm, Table S1 in Supplementary Materials) were prepared for investigating the effects of optical properties on reflectance prediction and comparing the estimation accuracy of diverse curve fitting methods. One optical property value was varied, while the other four were held constant. The optical property values for these samples were chosen based on published data [23; 24; 25], covering a large range of horticultural products with 0.001 mm\({}^{-1}\)\(\leq\)\(\mu_{a}\)\(\leq\) 0.1 mm\({}^{-1}\) and 0.5 mm\({}^{-1}\)\(\leq\)\(\mu_{s}{{}^{\prime}}\)\(\leq\) 4 mm\({}^{-1}\). More simulation samples were created (see Tables S2-S5 in Supplementary Materials for more details) for studying the effects of relative values of the optical properties of two layers (\(mfp_{1}{{}^{\prime}}\)/\(mfp_{2}{{}^{\prime}}\), \(\mu_{a1}\)/\(\mu_{a2}\) and \(\mu_{s1}{{}^{\prime}}\)/\(\mu_{s2}{{}^{\prime}}\)) and relative values of \(\mu_{s1}{{}^{\prime}}\)and \(\mu_{a1}\) (\(\mu_{s1}{{}^{\prime}}\)/\(\mu_{a1}\)), for optical property extraction from the top layer. The parameters of \(\mu_{a1}\)/\(\mu_{a2}\) and \(\mu_{s1}{{}^{\prime}}\)/\(\mu_{s2}{{}^{\prime}}\) were selected for studying the contributions of absorption and reduced scattering coefficients of two layers on estimating optical properties of the top layer, respectively, while \(mfp_{1}{{}^{\prime}}\)/\(mfp_{2}{{}^{\prime}}\) was used for studying the combination effect of absorption and reduced scattering coefficients.