## 2 Materials and Methods
### Image Reconstruction Algorithm  
The spatial distribution of absorption and scattering coefficients in finger tissue was reconstructed with the reduced Hessian sequential quadratic programming (rSQP) algorithm of PDE-constrained optimization [33] tailored to use the FD-SP\({}_{3}\) model as the forward model of light propagation in tissue.  
Following the framework of PDE-constrained optimization, the FD-SP\({}_{3}\) image reconstruction finds the composite moments vector \(u=(\phi_{1},\phi_{2})\) and the optical property vector \(\mu=(\mu_{a},\mu_{\rm s}^{\prime})\) in an all-at-once manner such that:  
\[\begin{array}{c}\text{min}f(u,\mu)=\frac{1}{2}(M-Qu)^{T}\overline{(M-Qu)}\\ s.t.Au=b\end{array} \tag{4}\]  
where \(Qu=P\) is the prediction of measurement \(M\) and the operator (\({}^{-}\)) denotes the complex conjugate of the complex value.  
Given the current estimate of forward and inverse variables \((u,\mu)\), the rSQP scheme gives the new iterate for both forward and inverse variables:  
\[(u,\mu)^{k+1}=(u,\mu)^{k}+\alpha^{k}\Delta(u,\mu) \tag{5}\]  
where a step length \(\alpha^{k}\) provides a sufficient decrease in the \(l_{1}\) merit function, and a search direction \(\Delta p=(\Delta u,\Delta\mu)^{T}\) can be obtained by solving the following quadratic programming problem:  
\[\begin{array}{c}\text{min}\Delta p^{T}\
abla f+\frac{1}{2}\Delta p^{T}W \Delta p\\ s.t.\
abla c^{T}\Delta p+(Au-b)=0\end{array} \tag{6}\]  
where \(\
abla f\) denotes the gradient of \(f(u,\mu)\), \(W\) denotes the reduced Hessian of the Lagrangian functional \(L(u,\mu,\lambda)=f(u,\mu)+\lambda^{T}(Au-b)\), and \(\
abla c\) represents the gradient of the constraint \(c=Au-b\). Details on the theory and implementation of the rSQP method can be found in [33].  
Note that the existing clinical data has no reference measurement taken--therefore, standard calibration is not possible. Instead, we normalized the data by dividing the measurement at each detector location by the average intensity taken over all source-detector measurements [34]. This procedure eliminates the effects by various unknown factors such as unknown power source and unknown detector response function. As a result, the normalized \(M\) and \(P\) are given as  
\[\hat{M}=\big{(}\overline{M}\big{)}^{-1}M_{i,j}\;;\;\hat{P}=\big{(}\overline{P }\big{)}^{-1}P_{i,j} \tag{7}\]  
The indices \(i\;\in\;S\) and \(j\;\in\;D\) denote all possible sources and detectors, respectively. Here, \(\overline{M}\) and \(\overline{P}\) represent the average measured intensity and the averaged predicted intensity, respectively, and are given by  
\[\overline{M}=\frac{1}{S\cdot D}\sum_{i=1}^{S}\sum_{j=1}^{D}M_{i,j}\;;\; \overline{P}=\frac{1}{S\cdot D}\sum_{i=1}^{S}\sum_{j=1}^{D}P_{i,j} \tag{8}\]  
To improve convergence, each element of the product \(\big{(}\hat{M}-\hat{P}\big{)}^{T}\overline{(\hat{M}-\hat{P})}\) in Equation (4) is further normalized by each element of the measurement norm \(\hat{M}\overline{\hat{M}}\) as  
\[f(u,\mu)=\frac{1}{2}\frac{\big{(}\hat{M}-\hat{P}\big{)}^{T}\overline{(\hat{M} -\hat{P})}}{\hat{M}\overline{\hat{M}}} \tag{9}\]  
which is implemented into Equations (4) and (6) of the rSQP algorithm with SP\({}_{3}\) model.