1. Introduction

Spectral imaging involves the acquisition of a series of 2D images of reflected light, where each image uses a different wavelength (λ) yielding a “spectral cube”, R(x,y,λ). Each pixel R(x,y) has a spectrum R(λ).

The reflectance spectrum can be described using a variety of light transport computations, and is here called getR(μ_a, μ_s’, n_r), for example diffusion theory [1], Monte Carlo [2], or the P_n method [3], which depends on the two tissue optical properties, the absorption coefficient μ_a [cm⁻¹] and the reduced scattering coefficient μ_s’ [cm⁻¹], and the refractive index mismatch (n_r) at the air/tissue surface. The behavior of getR(μ_a, μ_s’,n_r) is mimicked by the expression:

The method can consider fat/lipid and other absorbers. For this paper, only the blood and water absorption are considered, to simplify the notation. The reduced scattering coefficient is described by an expression that matches experimental data [4,5]:

This paper considers reflectance from skin, so an additional parameter due to epidermal melanin is included, exp(-Mμ_a.melL_epi), where

Also, any variation in the coupling of reflected light from skin into the camera, relative to coupling of reflected light from a reflectance standard into the camera, is considered by a factor K. The expression for R(λ) is:

The factor K is wavelength independent and equals the ratio of reflectance collection from the skin (C_skin) and the reflectance standard (C_std), C_skin/C_std, which may also include a factor cosθ where θ is the angle of observation off the normal to the skin surface (typically θ = 0° if the camera views perpendicular to the skin surface, and K = 1).

The analysis of R(λ) at each x,y pixel can use Eq. (3) in a least-squares fitting routine to find the best values of B,S,W, μ_s’_.500nm, f, and b that specify the μ_a and μ_s’ for getR() to match the data.

In the MATLAB programming environment, we have previously used the function fminsearch() for this least-squares fitting of the spectrum from each pixel. However, the fitting is slow, about 120 ms per pixel, or 8.7 hr to analyze 262,144 (512x512) pixels. This report discusses a faster method that reduces the analysis 100-fold down to 5.2 min.

2. The method

The method considers the optical depth, OD = -ln(R), so that Eq. (4) becomes:

The functional forms of L and L_epi are not discussed in this report, since such a discussion should include a full description of the Monte Carlo simulations.

The algorithm must assume a wavelength dependence for μ_s’, and for skin the literature [4,5] suggests μ_s’_500nm = 43 cm⁻¹, f = 0.64 and b = 0.91. The μ_s’_500nm is allowed to vary so as to scale the strength of this scattering to match the scattering of a particular skin site, but the initial choice is 43 cm⁻¹. Initial values are also assumed for B, S, W, M and K to specify X. Then the average volumetric skin absorption is calculated using the B, S and W values: μ_a = A(:,1:3)*X(1:3), equivalent to Eq. (1). This μ_a spectrum, the μ_s’ using Eq. (3), and the melanin optical depth Mμ_a.mel allow calculation of L(λ) and L_epi(λ), which are inserted into matrix A [Eq. (6)]. Solving for X using Eq. (7) specifies new values of X, which are then used to recalculate L(λ) and L_epi(λ) and update A. These steps are iteratively repeated several times. Typically, the solution for X converges after 5 iterations.

There is one additional step within these iterative steps. The choice of μ_s’_500nm scales the strength of scattering. Since total reflectance R is nonlinearly related to the ratio μ_s’/μ_a, an error in μ_s’ can be compensated by an error in the blood content B and water content W, such that the ratio μ_s’/μ_a yields the correct R. If one chooses a series of μ_s’_500nm values and solves for X, the values B and W will increase linearly versus the chosen value μ_s’_500nm. This linearity can be made useful, if one has a reliable estimate of the true tissue water content. The method is as follows.

At each iteration, the new X specifies the water content, W₁ = X(3). The true water content of the skin is assumed to be 0.65. Therefore, the value of μ_s’_500nm is updated by multiplying the current value of μ_s’_500nm by the ratio 0.65/W₁. After several iterations, the solution for X becomes stable, and μ_s’_500nm equals a value that yields W₁ = 0.65. In other words, the water content of the tissue has been used as an internal standard to allow specification of the strength of scattering, μ_s’_500nm. Consequently, the value of B is also correct. The values of S, M and -ln(K) are relatively insensitive to the choice of μ_s’_500nm.

This algorithm is about 100-fold faster than using fminsearch(), 1.2 ms versus 120 ms, to fit a total reflectance spectrum R(λ) for the parameters B, S, M, K and μ_s’_500nm. The advantage of this method is speed. The disadvantage is that assumed values for W, f and b are required. The algorithm is outlined in Fig. 1 .

 

Fig. 1 The algorithm for spectral analysis. The measured spectrum specifies OD. The array A [see Eq. (5) is initially filled using Eq. (1)] to specify μ_a and using the lookup tables L(μ_a, μ_s’, Mμ_a.mel) and L_epi(μ_a, μ_s’, Mμ_a.mel) based on Monte Carlo simulations to specify L and L_epi for each wavelength. An initial value μ_s’_500nm is assumed. The matrix inversion X = A⁻¹OD yields values for B, S, W that are used to update μ_a. The predicted water W₁ is used to update μ_s’_500nm: μ_s’_500nm = μ_s’_500nm(0.65/W₁). Then μ_a, M and μ_s’_500nm are used to update L and L_epi. The A matrix is updated. The cycle is repeated. After 5 iterations, the values of X converge on stable values.

Download Full Size | PDF

The error in the estimate of tissue water content W is likely to be only ~10%. The values of f and b can be learned for a particular tissue type by slower nonlinear least-squares fitting that allows fitting for f and b. Once these are known for a tissue type, the values of f and b can be used in the faster algorithm. The issue of μ_s’_500nm and B and W interacting via the ratio μ_s’/μ_a is still an issue with slow nonlinear least-squares fitting. The strategy of using water as an internal standard helps both the fast and slow analyses.

3. Experimental example: portwine stain lesion

To illustrate the spectral analysis, a reflectance spectrum from a portwine stain lesion (PWS) was acquired with a hand-held spectral camera (see Fig. 2 ). An optical fiber bundle delivers light through a linear polarizer to illuminate an open port that contacts the skin. A CCD camera (12 bit, Flea, Point Grey Research, Richmond, BC, Canada) views the skin through a second cross-polarized linear polarizer, thereby avoiding surface glare and collecting diffusely backscattered light. The light source is a white light source passing through a liquid crystal tunable filter (LCTF) (CRI Inc., Cambridge, MA, USA). Software control written in Labview^TM scans the LCTF over a set of 21 wavelengths (480:10:650 nm, plus 545, 558, and 585 nm).

 

Fig. 2 A hand-held spectral camera for measuring skin sites. Liquid crystal tunable filter scans to different wavelengths. Camera views reflected skin through cross-polarizer filter to reject surface glare. The open port contacting the skin does not compress the skin site and avoids blanching the vasculature.

Download Full Size | PDF

Spectral images were acquired before and after treatment with a pulsed dye laser, a standard treatment for PWS. The laser treatment elicits thrombi formation in the PWS, causing ischemia, which elicits a wound healing response that removes the PWS. However, if thrombi loosen and blood flow is re-established, then the treatment fails. The goal of the project is to image regions where such reperfusion has occurred and predict treatment failure. In the clinical study, spectral images are typically acquired 15-30 min after treatment. Regions of reperfusion seen in images are later correlated with the success of the PWS clearance at 6 weeks post treatment.

Figure 3 shows the PWS spectra acquired before and 35 min after laser treatment. The laser treatment caused an increase in blood (B) and an decrease in oxygen saturation (S), but no significant effect on the melanin (M) or the constant (K).

 

Fig. 3 Portwine stain lesion (PWS) spectra taken (A) before laser treatment and (B) 35 min after laser treatment (black circles). Red line shows fit by algorithm. Magenta lines show oxygenated blood (S = 0.6-1.0). Cyan lines show deoxygenated blood (S = 0-0.5). The blood content (B) increases and the oxygen saturation (S) decreases as thrombi are formed in the PWS by the laser treatment. The melanin content (M) and the constant (K) did not change significantly. The spectra for scatter only (black), + water (blue), + melanin (green) and + blood (red) are shown.

Download Full Size | PDF

4. Speed improvement for spectral analysis

To illustrate the improvement in speed, a set of spectra were analyzed using the least squares fitting routine fminsearch() in MATLAB^TM, which is based on the Nelder-Mead simplex method. Then the same spectra were analyzed using the algorithm of this report, using 5 iterations. Figure 4 shows histograms of computation times for the spectra using the two methods. The matrix method for spectral analysis of this report required 1.2 ms, while the fminsearch() required 120 ms, which is a 100-fold increase in speed.

 

Fig. 4 Time required for fitting spectra, using (A) the matrix method of this report, and (B) the fminsearch() function in MATLAB^TM. The histograms for computation time (ms) for fitting one spectrum is shown, based on fitting over 300 spectra. The mean time for the matrix method (1.2 ms) is 100-fold faster than the mean time for the fminsearch() method (120 ms).

Download Full Size | PDF

4. Summary

The matrix inversion method provides a rapid spectral analysis for spectral imaging. A 512x512 pixel image would require about 5.2 min to analyze all pixels using a single processor and implemented on MATLAB^TM to yield the blood, oxygen saturation, melanin content and scattering strength of a skin site. This method is being used to analyze spectral images of portwine stain lesions before and after laser treatment to detect reperfusion and predict treatment failure.

This matrix inversion method should be able to be implemented on a graphics programming unit (GPU) to attain another factor of 10 or more in speed. Combining the 100-fold increase in speed of the rapid analysis algorithm with a high speed GPU should enable near real-time spectroscopy of tissues. For spectral imaging, analyzing an image would take ~30 s.