# GPU Adding-Doubling Algorithm for Analysis of Optical Spectral Images
## 2 Materials and Methods
### Implementation of the Adding-Doubling (AD) Algorithm
The algorithm is initialized by calling the function **initPR** with the following syntax:  
\[\text{obj}=\text{obj}.\text{initPR(g, ['alpha', alpha], ['SPFtype', SPFtype], ['R0L',  R0L]['nE', nE], ['nS', nS], ['m', m], ['Nsam', Nsam])  
where \(g\) is a matrix of tissue anisotropies at the specified wavelengths _lam_ and samples; the parameters in [] are optional. These optional parameters include:  
1. _alpha_--the Gegenbauer scattering phase function (SPF) exponent at the specified wavelengths _lam_ for each sample. It is empty by default;
2. _SPFtype_--the type of SPF to be used. Currently implemented options are Henyey-Greenstein (HG) SPF and Gegenbauer (GB) SPF. To select HG SPF, the parameter value is 0, and to select GB SPF, the parameter value is 1. If GB SPF is selected, the _alpha_ parameter must also be specified;
3. _R0L_--the reflectance of the bottom substrate, ranging between 0 and 1, if it is different from the upper environmental medium. This parameter is used in cases where a sample is placed on top of another thick substrate, such as a support plate. _R0L_ is  
Figure 3: Flowchart of the proposed GPU-accelerated Adding-Doubling algorithm, which runs exclusively on a graphics processing unit (GPU). The first step, “Initialize object”, a computationally intensive process, is executed only once during the algorithm initialization. The “Set tissue properties” and “Perform Adding Doubling” steps are run whenever there is a change in tissue properties. A detailed description of each algorithm step can be found in Prahl [12].  
a matrix with \(Nl\) rows corresponding to \(lam\) and \(Ns\) columns corresponding to the number of samples;
4. \(nE\)--a column vector of environmental medium refraction indices with \(Nl\) values corresponding to \(lam\). By default, the value is set to 1, representing air;
5. \(nS\)--a column vector of substrate refraction indices with \(Nl\) values corresponding to \(lam\). By default, the value is set to 1.4.+;
6. \(m\)--the number of fluxes (even integer \(>=2\)) used in the algorithm. The default value is 20, as it has been found to provide adequate accuracy for most combinations of optical parameters;
7. \(Nsam\)--the number of different spectra being calculated. By default, it is set to 1, meaning that only one set of different wavelengths is considered.  
An example of initializing the algorithm, providing \(g\) for Gegenbauer SPF and 20 fluxes, can be performed as follows:  
A.initPR(g,'SPFtype',1,'alpha',alpha,'m',20);  
where \(g\) is a column vector of \(Nl\) anisotropy factors, \(alpha\) is a column vector of \(Nl\) GB SPF exponents, and \(m\) is the number of fluxes.  
The variable properties of the tissues are set by calling the **setOP** method. The calling syntax is  
\[\text{obj}=\text{obj}.\text{setOP}(\text{muA},\text{muS},\text{d})\]  
where \(muA\) and \(muS\) represent the matrices of absorption and scattering coefficients in units of [1/cm], and \(d\) is the layer thickness in units of [cm]. The dimensions of \(muA\) and \(muS\) are \(Nl\) rows and \(Ns\) columns, while \(d\) has \(Ns\) values, one for each sample. An example of calling this method is  
A.setOP(muAar,muSar,d);  
Finally, the **solveAD** is called to calculate the reflectance and transmittance values at the specified wavelengths. The calling syntax is  
obj = obj.solveAD()  
Although this method does not directly return any values, it stores the results in the class properties. Currently, only normal incidence of the illumination and complete diffuse reflectance and transmittance (e.g., all radiation in the exit hemispheres, two solid angles) are supported. The results can be retrieved through the following properties:  
* _obj.Rn_--2 diffuse reflectance;
* _obj.Tn_--2 transmittance;
* _obj.Rsp_--specular reflectance. An example of retrieving all these parameters in a MATLAB environment is:  
[A.Rn, A.Tn, A.Rsp]