### Inverse Algorithm for Estimating Optical Properties of Two-Layered Samples  
Prior to determining the optical properties of two-layered simulation samples from the spatial-frequency domain reflectance, the effect of varied \(\mu_{a}\) and \(\mu_{s}{{}^{\prime}}\) of each layer on reflectance prediction was investigated, which would be helpful for understanding the two-layered diffusion model. To select the most proper, robust and accurate parameter estimation method for determining the optical properties of two-layered samples, different curve fitting methods that differ in the number of free variables were compared and evaluated (Table 1). First, all five parameters (\(\mu_{a1}\), \(\mu_{s1}{{}^{\prime}}\), \(\mu_{a2}\), \(\mu_{s2}{{}^{\prime}}\) and \(d\)) were estimated simultaneously using the two-layered model in Equation (2) (five-variable fit). Second, the top-layer thickness was assumed to be known and the other four parameters were estimated at once (four-variable fit). Third, \(\mu_{a}\) and \(\mu_{s}{{}^{\prime}}\) of either the top layer or the bottom layer were estimated, depending on which layer had changed variables, while \(\mu_{a}\) and \(\mu_{s}{{}^{\prime}}\) of the other layer and thickness the top layer were treated as known (two-variable fit). Fourth, only \(\mu_{a}\) or \(\mu_{s}{{}^{\prime}}\) was estimated, while all others were treated as known (one-variable fit). Fifth and finally, the one-layered model in Equation (1) was utilized to estimate \(\mu_{a}\) and \(\mu_{s}{{}^{\prime}}\) of the top layer, which was used in our proposed stepwise method [13]. The estimated optical properties were compared with the true values for evaluating the estimation accuracy. A nonlinear least-squares fitting method was used to minimize the sum-of-squares of the difference between the true reflectance and predicted reflectance values with estimated parameters using Equation (1) or Equation (2). In this study, a subspace trust-region method based on the interior-reflective Newton approach was used to achieve the algorithm optimization [26], which is defined by minimizing a quadratic function subject to an ellipsoidal constraint. The approach can generate iterates in the strictly feasible region by using a new affine scaling transformation, and the speed of convergence is accelerated by following a reflective line search technique. The optical property estimation procedure from the spatial-frequency domain reflectance of two-layered simulation samples was implemented using the Toolbox function 'Isqcurvefit' in Matlab 8.4 (The MathWorks, Inc., Natick, MA, USA). Previous studies revealed that there were inherent discrepancies of the reflectance between the diffusion model and MC simulation, which would greatly affect the accuracy of optical property estimation [16; 27]. It was proven that sample-based calibration method could decrease the reflectance discrepancies and improve the estimation accuracy. The method was implemented by first selecting a set of calibration samples depending on the initial estimated values of \(\mu_{s1}\)\({}^{\prime}\). Ratios of diffusion model generated reflectance to MC-generated reflectance for all the calibration samples were obtained, and they were then used for calculating the calibrated reflectance by multiplying the original reflectance for each sample. Since this study is focused on theoretical analysis of a two-layered diffusion model, the sample-based calibration method was thus applied for reflectance correction, which was then used for optical property estimation. More details about the calibration method are referred to [PERSON] et al. (2019) [13].