Depth profiling of laser-heated chromophores in biological tissues by pulsed photothermal radiometry.

A solution method is proposed to the inverse problem of determining the unknown initial temperature distribution in a laser-exposed test material from measurements provided by infrared radiometry. A Fredholm integral equation of the first kind is derived that relates the temporal evolution of the infrared signal amplitude to the unknown initial temperature distribution in the exposed test material. The singular-value decomposition is used to demonstrate the severely ill-posed nature of the derived inverse problem. Three inversion methods are used to estimate solutions for the initial temperature distribution. A nonnegatively constrained conjugate-gradient algorithm using early termination is found superior to unconstrained inversion methods and is applied to image the depth of laser-heated chromophores in human skin.


INTRODUCTION
Pulsed photothermal radiometry (PPTR) is a noncontact technique that utilizes an infrared detector to measure temperature changes induced in a test material exposed to pulsed radiation. Heat generated as a result of light absorption by subsurface chromophores in the test material diffuses to the surface and results in increased infrared radiant emission levels. By collecting and concentrating the emitted radiation onto an infrared detector, one obtains a PPTR signal that represents the time evolution of temperature near the test material's surface. Useful information regarding the test material may be deduced from analysis of the PPTR signal; reported applications include (1) evaluation of surface coating thickness in industrial components, 1 (2) identification of subsurface microcracks in aircraft structures, 2 (3) determination of the optical absorption coefficients in human artery 3 and biliary calculi, 4 and (4) characterization of port-wine stain birthmarks. 5 Because coupled radiative-thermal effects are present during and after pulsed irradiation, determination of the optical and/or the thermal properties of the test material from a PPTR measurement can be complex. In the following analysis we assume that the thermal diffusivity in a one-dimensional layered test material is constant following pulsed laser irradiation and derive an integral equation relating measured PPTR signal amplitude to the initial temperature distribution. Because we seek only to find the initial temperature distribution immediately following laser irradiation, analysis of the pulsed radiative transfer in the test material is not included. We deduce the initial temperature distribution from the measured PPTR signal amplitude by finding a solution estimate to the derived integral equation. Computed singular values of the kernel or point-spread function in the derived integral equation demonstrate the severely ill-posed nature of the PPTR inverse problem.
Since the derived integral equation is ill posed, estimates of the initial temperature distribution can change substantially as a result of small perturbations (i.e., noise) in the measured PPTR signal amplitude. We present three inversion methods that are used to compute solution estimates of the initial temperature distribution. The singular-value decomposition is used to analyze the characteristic properties of the kernel function. First, we use a least-squares truncated singularvalue decomposition (LSTSVD) to compute a solution estimate for the unknown initial temperature distribution in the exposed test material. Although inferior to those obtained from the conjugate-gradient methods discussed below, results from LSTSVD demonstrate the effect of the PPTR signal-to-noise ratio on (1) the degree of regularization, (2) the accuracy of the estimated solutions, and (3) the instabilities that arise from small singular values of the kernel function. Second, we use the method of unconstrained conjugate gradients to find solution estimates to the PPTR inverse problem. A shortcoming of unconstrained methods is the presence of unrealistic negative temperature changes in computed solution estimates. Third, we apply a nonnegatively constrained conjugate-gradient algorithm 6 to determine the initial temperature distribution from a measured PPTR signal. Because temperature changes in the test material resulting from radiative absorption are strictly greater than or equal to zero, we constrain the inversion method to search for nonnegative solutions of the integral equation and obtain estimates of greater accuracy.

A. PPTR Signal
We derive an expression for PPTR signal amplitude ͓DS͑t͔͒ in terms of the initial temperature distribution ͓DT ͑z, t 0͔͒ in a test material immediately following pulsed laser irradiation. For the purpose of our analysis we assume that the test material occupies semi-infinite half-space and that the plane positioned at z 0 is coincident with the boundary surface ( Fig. 1). Radiative and convective thermal losses at the air -material interface are modeled with a Robin boundary condition [Eq. (1)]: ( At the test material surface ͑z 0͒, heat flux Q 2kdDT͞dz is proportional to the surface temperature increase Q 2khDT ; proportionality constants include thermal conductivity k and heat-loss coefficient h. Boundary conditions that represent a perfectly insulating ͑Q 0͒ and constant-temperature ͓DT ͑0, t͒ 0͔ surface correspond, respectively, to limiting values of h 0 and h `. From the Green's function solution to the one-dimensional heat-transfer problem 7 the temperature distribution DT ͑z, t͒ at an arbitrary depth z and time t is written as where erfc͑u͒ is the complementary error function and u is defined [Eq. (3)] in terms of the boundary loss coefficient h and the thermal diffusivity D, which we assume is homogeneous: In a physical sense, the three terms in the integrand (from left to right) of Eq. (2) represent (1) a planar thermal source with amplitude DT ͑z 0 , 0͒ located at position z 0 , (2) an equal-amplitude image source located at position 2z 0 , and (3) a continuous distribution of planar image sinks extending from z 2z 0 to z 2`of strength 2DTh exp͓h͑z 1 z 0 ͔͒dz that account for thermal energy losses from the test material surface. The ideal PPTR signal amplitude DS͑t͒ at any time t is proportional to a superposition [Eq. (4)] of the depth-dependent temperature changes DT ͑z, t͒ weighted by the transmission loss in the exposed test material that is due to infrared absorption: Here m a is the infrared absorption coefficient in the exposed test material at the detected wavelength(s) and C d is a proportionality constant determined by the infrared detection system. In addition, we have assumed in Eq. (4) that observed temperature changes are much less than the initial background temperature ͑DT ,, T 0 ͒, so increases in the infrared radiant exitance, DE͑ T 0 1 DT͒ 4 2 T 0 4 ഠ 4T 0 3 DT, are linear with temperature change. An expression for the ideal PPTR signal amplitude is found [Eq. (5)] by substitution of the temperature distribution DT ͑z, t͒ defined in Eq. (2) into Eq. (4) and completion of the z integral: Here, erfcx͑u͒ exp͑u 2 ͒erfc͑u͒ is the exponential complementary error function and u 6,1 are functions [Eq. (6)] of space and time and the relevant material coefficients D, h, and m a : We view the measured PPTR signal amplitude as a blurred image of the initial spatial temperature distribution; our objective is to determine DT ͑z 0 ͒ from the blurred image ͓DS͑t͔͒.

B. Inverse Problem
Equation (5) is a Fredholm integral equation of the first kind for the unknown initial temperature distribution, Fig. 1. Geometry assumed for analysis of the inversion problem; measured PPTR signal ͓DS͑t͔͒ and initial temperature distribution ͓DT ͑z, t 0͔͒ solution.
DT ͑z 0 , 0͒. Determination of a solution for DT ͑z 0 , 0͒ requires evaluation of the kernel or point-spread function [Eq. (7)] and measurement of the PPTR signal amplitude DS͑t͒: Inasmuch as Fredholm integral equations of the first kind arise in numerous remote-sensing problems of practical interest, many inversion methods have been formulated for their solution. 8 Most problems are ill posed in the sense that small perturbations (i.e., noise) in the measured data can result in large variations in a computed solution estimate. In the following, we seek a solution for the initial temperature distribution DT ͑z 0 , 0͒, using a known kernel function K͑z 0 , t͒ and a measured PPTR signal DS͑t͒, which is corrupted with noise. Knowledge of the initial temperature distribution is the first step toward determining the spatial distribution of laser-heated subsurface chromophores in the exposed test material. Because all observed physical systems (i.e., source, transmission channel, and receiver) exhibit intrinsic random fluctuations that contribute to measurement noise, the ideal PPTR signal DS͑t͒ obtained by numerical evaluation of Eq. (5) is a mathematical construct that does not exist in practice. Noise present in an infrared detection system originates from a number of sources 9 : (1) fluctuations in the arrival rate of background photons, (2) shot noise that is due to the discrete nature of the photodetection process, (3) thermal-induced fluctuations or Johnson noise present in any resistive elements utilized in the readout electronics, (4) generation -recombination noise that arises when photoinduced carriers recombine at impurity sites within the detector element, (5) phonon noise that is due to temperature fluctuations in the detector element, and (6) 1͞f noise that is due to processes not completely understood. Inasmuch as multiple noise sources are present, we invoke the central-limit theorem 10 and add a zero-mean, Gaussian distributed noise value e to the ideal PPTR signal amplitude DS͑t͒ in Eq. (5). The variance of the noise distribution is determined by the signal-to-noise ratio [SNR; Eq. (8)] of the infrared detection system: where ‫͘ء͗‬ represents a time average of the quantity *. We begin seeking a solution estimate by approximating the PPTR signal, the initial temperature change in the exposed test material, and the kernel function with respective discrete vector and matrix quantities, DS DS͑t i ͒, DT DT ͑z j 0 , 0͒, and K K͑t i , z j 0 ͒Dz 0 ; additionally, we assume that t i and z j are uniformly spaced, so (1) t i 2 t i21 t net ͑͞N t 2 1͒, with t 1 . 0, (2) Dz z i 2 z i21 ͑͞N z 2 1͒, with z 1 . 0, and (3) N t $ N z so the problem is overdetermined. Boldface terms are used to represent vector and/or matrix quantities. The measured PPTR signal amplitude DS differs from the ideal DS by zeromean, white Gaussian noise vector e. The linear matrix equation that approximates the continuous integral The simplest means to estimate a discrete approximation of the initial temperature distribution is to select the DT that minimizes the Euclidean norm of the residual vector ͑jjrjj jjDS 2 KDT jj͒, which is the minimum-variance unbiased least-squares estimate. Unfortunately, because Eq. (9) is ill posed, the high variance of this estimate makes it essentially useless. One obtains meaningful estimates by augmenting the unbiased least-squares problem with various constraints. We consider three increasingly sophisticated methods that permit computation of initial temperature distribution estimates DT .

INVERSION METHODS
To determine the degree to which a problem is ill posed, a first step 11 is examination of the singular values of the kernel matrix K. The fundamental theorem of linear algebra 12 provides the mathematical basis for the existence and characteristic properties of the singular values and vectors for an arbitrary matrix over a finitedimensional real vector space. The singular-value decomposition 13 allows for an arbitrary kernel matrix K of rank r to be expanded [Eq. (10)] in terms of bases of orthonormal left ͑s i ͒ and right ͑t i ͒ singular vectors: The basic difficulty in computing a solution estimate for the initial temperature distribution DT in Eq. (9) stems from the relatively large number of small singular values d i of the kernel matrix K. Any selected inversion method can be viewed as forming a linear superposition of right singular vectors of the kernel matrix K for the unknown initial temperature distribution DT . In the superposition, each right singular vector is weighted by the inverse of an associated singular value. Large oscillations in such a solution estimate can occur if expansion terms are not selected judiciously. Because smaller singular values correspond to higher-order right singular vectors of increasing spatial frequency (see, for example, Fig. 5 below), efforts to extract edges and sharp features of the unknown initial temperature distribution DT can be difficult and lead to instabilities. Thus any realistic inversion method must deal with the instability problem that originates from small singular values of the kernel matrix and their associated high-spatial-frequency right singular vectors.

A. Least-Squares Truncated Singular-Value Decomposition
If there exists some integer m , r so that the singular values d l ͑l . m͒ are small [Eq. (11)], expansion of the kernel matrix in Eq. (10) for values greater than m is not physically meaningful because the energy in successive terms is less than the noise energy N t ͗e i 2 ͘: We use Eq. (11) to truncate the singular-value expansion [Eq. (10)] after the first m terms and seek a solution estimate that is a linear combination of the first m orthonormal right singular vectors t i , with coefficients c i selected to minimize the Euclidean norm of the residual, f ͑c 1 , c 2 , . . . , c m ͒, in Eq. (12): The concept of the LSTSVD method is to truncate the singular-value decomposition expansion [Eq. (10)] before small singular values [Eq. (11)] begin to dominate the estimate of the initial temperature distribution. Later we apply the LSTSVD method to compute solution estimates of the unknown initial temperature distribution in a model problem.
The truncation integer m in Eq. (11) can be viewed as analogous to a regularization parameter. The optimum value of m cannot be computed directly from Eq. (11) because knowledge of the unknown initial temperature distribution DT is required. We determine the truncation integer by using an L-curve analysis. 14 In this method we find the optimum truncation integer by plotting the Euclidean 2-norms of the computed solution estimate jjDT ͑ p͒jj versus that of the residual vector jjDS 2 KDT ͑ p͒jj for a number of truncation integers p. The value of p that corresponds to the corner of the resulting curve -which often resembles the letter L -is selected as the optimum truncation integer m.

B. Method of Conjugate Gradients
Originally developed by Hestenes and Stiefel, 15 the method of conjugate gradients is an iterative technique that can be successfully applied to minimize f ͑DT , l͒: or, equivalently, The additional regularization term ljjLDT jj 2 is added to the squared norm of the residual, jjKDT 2 DSjj 2 , to penalize large solution estimates. First, we consider the unconstrained case in which the unknown solution can represent positive or negative temperature changes. Second, we find a better approach in a constrained conjugate-gradient algorithm 6 that requires that the computed estimate for the initial temperature distribution be nonnegative (i.e., DT $ 0). The nonnegativity constraint permits computation of physically realistic (i.e., representative of radiative absorption) and more-accurate estimates of the initial temperature distribution.
When constrained and unconstrained conjugategradient methods are used the quality of the estimate strongly depends on the degree of regularization as specified by the real parameter l. A number of techniques and philosophies exist for specifying the degree of regularization. The most common are (1) the discrepancy principle, 16 (2) generalized cross validation, 17 and (3) the L curve. 15 The L curve has intuitive appeal because it graphically represents the variance-versus-bias trade-off that is involved in specifying the degree of regularization. For small l the estimate for DT overfits the measured PPTR signal and models measurement noise that is inherent in the infrared detection process. The result is a large, noisy, high-variance estimate for DT . For large l the result is a small, smooth, high-bias estimate that does not adequately model the initial temperature distribution; consequently the residual norm is very large. Plotting the estimate norm jjDT ͑l͒jj against the residual norm jjKDT ͑l͒ 2 DSjj was first suggested by Lawson and Hanson 14 and yields an L-shaped plot, where the optimum regularization parameter is at the corner of the L. Use of the L-curve technique has an advantage over use of the discrepancy principle in that prior knowledge of the noise variance is not required. The advantages over generalized cross validation are that it is easily adapted to the nonnegatively constrained case and appears not to require the same degree of orthogonality between the signal and noise. Various simulations 18 suggest the advantages of the L-curve technique over generalized cross validation.
A method that is much faster than either LSTSVD or solving Eq. (13b) is to apply the conjugate gradientalgorithm to the unaugmented problem ͑l 0͒ of Eq. (13a), starting with DT 0, and to regularize by early termination. In the chemometrics literature the method is called partial least squares and has been given a sound theoretical basis. 19,20 In these references the best iterate is selected by use of ordinary or generalized cross validation. 21 Furthermore, both the discrepancy principle and L curve can be used to determine the best iterate. Although the theory behind this approach is too complex to discuss here, we have observed in numerous computations that this method is fast and works well.

METHODS AND MATERIALS
We formulate a model PPTR problem that allows for numerical testing of each of the three inversion methods. We assume a temperature increase DT i ͑z, t 0͒ in the test material as a result of radiative absorption in an ith subsurface layer [Eq. (14)] and derive an analytic expression [Eq. (15)] for the corresponding ideal PPTR signal ͓DS i ͑t͔͒: 2m a k͑m a 1 h͒erfcx͑u 2 ͒ ͑m a 1 k͒͑m a 2 k͒͑h 2 k͒ We derive Eq. (15) by solving the Green's function integral [Eq. (2)] for the time-dependent temperature distribution DT ͑z, t͒ and computing the ideal PPTR signal [Eq. (4)]; u 6,1 are defined in Eqs. (6), u 2 in Eq. (16): Although noise in any infrared detection system is due to a number of sources, best system performance is attained when the dominant source originates from intrinsic fluctuations in the photon emission rate (i.e., radiant power) from the background field; in this case the detection system is said to be photon-noise limited or a background-limited photodetector (BLIP). In practice, BLIP systems require more stringent design criteria, which frequently include cooling the detector element to liquid-nitrogen or -helium temperatures and use of an ultralow-noise electronic preamplifier.
For the purpose of determining the optimum SNR we examine a BLIP infrared detection system and write an expression 9 for SNR blip : SNR blip is expressed [Eq. (17)] in terms of the background-limited normalized spectral detectivity D ‫ء‬ ͑l͒, the temperature derivative of the Planck expression for the spectral radiant exitance M l ͑l, T ͒ of a unit emissivity blackbody, the transmission coefficient t͑l͒ from the test material surface to the detector element (i.e., through the intervening air and optical components), the infrared region of spectral sensitivity ͑l 1 , l 2 ͒, the mean signal increase above background levels at the test material surface in response to pulsed radiative exposure ͑͗DS͒͘, the collection solid angle V 0 in object space, the detector surface area A d , the frequency bandwidth of the detection system Df , and geometric magnification m of the optical system. We compute SNR blip for two infrared detectors with a peak responsivity wavelength (l p 11 mm) that are sensitive in the 7 -11mm spectral range. For liquid-nitrogen-(77 K) and thermoelectric-͑275 ± C͒ cooled detectors, the respective normalized spectral detectivities are D ‫ء‬ ͑l p ͒ 3 3 10 10 and 2 3 10 8 ͑cm Hz 1/2 W 21 ͒. We assume a measured temperature increase ͗DS͘ 5 ± C and a unit magnification optical system with a 25-mmdiameter, f͞1 lens stopped down to a 6-mm-diameter entrance pupil that is antireflection coated so that ͗t͑l͒͘ ഠ 0.8. For the purpose of model inversion calculations, allowance for the presence of multiple noise sources is made for a liquid-nitrogen-cooled infrared detector that reduces the SNR up to three times below the BLIP limit so that SNR 100, corresponding to a Df 1 kHz system bandwidth (Fig. 2). The SNR is compromised by 2 orders of magnitude when a thermoelectric-cooled detector is used; for a 200-Hz system bandwidth the optimum SNR is 10. A numerical routine 22 that computes random numbers following a zero-mean Gaussian white-noise distribution is used to generate noise values e i . The variance s 2 of the Gaussian distribution is fixed to require that s ͗e i 2 ͘ 1/2 ͗DS i ͘͞SNR. Accurate evaluation of Eq. (15) requires realistic estimates of the diffusion constant D, the infrared absorption coefficient m a , and the heat-loss coefficient h. Inasmuch as our primary interest concerns the response in biological tissue immediately following pulsed laser irradiation, we assume a diffusion constant ͑D 0.11 mm 2 s 21 ͒ representative of human skin. 23 The infrared attenuation coefficient in human skin displays a strong dependence on water and protein. Except in the superficial stratum corneum (upper 5 -10 mm), water is the dominant biochemical species in human skin and is present at a concentration 24 of 70%. We assume a 70% waterdominated absorption spectrum, 25,26 which is weighted by the detector responsivity and blackbody emission and averaged over a spectral bandwidth (7 -11 mm) to yield m a ഠ 50 mm 21 . In the absence of surface cooling agents, heat loss from the skin surface is due primarily to radiation and free-air convection. We assume skin emissivities close to unity and typical free convection heat-transfer coefficients 27 to find h ഠ 0.03 mm 21 , representative of 20-Wm 22 K 21 surface losses.
A simulated PPTR signal is computed to test each inversion method. A temperature increase DT 1 20 ± C in a thin ͑z 2 2 z 1 150 mm͒, superficial ͑z 1 300 mm͒, and absorptive ͑k 3 mm 21 ͒ layer (Fig. 3A) is assumed to give a simulated signal (Fig. 3B). The size of the kernel matrix K used in each inversion computation is fixed at 256 3 128, corresponding to a digitized signal DS sampled at 8-ms intervals over 2 s and an initial temperature distribution vector DT with 8-mm spatial resolution extending to a depth of 1 mm.
Error ͑e͒ of a solution estimate DT is computed [Eq. (18)] with respect to the true solution DT ± by use of the L 1 vector norm: e jjDT 2 DT ± jj 1 jjDT ± jj 1 .
Error computations allow for objective comparison of each inversion method. The nonnegatively constrained conjugate-gradient algorithm is applied to solve a biomedical problem frequently encountered in laser treatment of port-wine stain (PWS) birthmarks. 28 In this approach, pulsed laser light partially transmits through the epidermis and is preferentially absorbed by hemoglobin (the major chromophore in blood) contained in ectactic blood vessels located in the upper dermis. There the incoming laser energy is converted to heat, causing thermal injury and thrombosis in the targeted blood vessels. The normal overlying epidermis is not totally spared because of the partial absorption of pulsed laser energy therein by melanin. For improved Fig. 2. SNR blip of liquid-nitrogen-͑· · · · · ·͒ and thermoelectric-(----) cooled infrared detection systems; see text for the specific configuration. laser treatments we seek to know the position and the relative concentration of melanin in the epidermis and of hemoglobin with PWS blood vessels in the upper dermis. PPTR has been used to measure the level of infrared radiant emission from PWS skin in response to pulsed laser exposure. 5 One motivation for the present study is determination of the initial temperature distribution in PWS skin immediately following pulsed laser exposure -given only a diagnostic PPTR measurement. Knowledge of the initial temperature distribution gives diagnostic information that can be used for designing improved laser treatments.
To aid in interpreting clinical results, we first consider a two-layer PWS model that consists of a thin ͑z 2 2 z 1 40 mm͒ superficial layer ͑z 1 10 mm deep) of melanin in the epidermis overlying a thicker ͑z 2 2 z 1 650 mm͒ and deeper ͑z 2 350 mm͒ plexus of laser-heated PWS blood vessels in the upper dermis. Because incoming laser radiation is strongly scattered, absorption of incident pulsed laser light in the melanin layer is augmented by backscattered light 29 and is assumed to give a nearly uniform ͑k ഠ 0͒ temperature increase ͑DT 1 40 ± C͒. We assume an initial temperature increase DT 2 20 ± C in the PWS layer that decreases with depth ͑k 5 mm 21 ͒ as a result of light absorption and scattering (Fig. 9 below, solid curve). Then we apply the inversion problem to a PPTR signal recorded following pulsed laser exposure of in vivo PWS skin. The recorded signal represents the infrared response collected over a 1-mm-diameter circular area on the skin surface following pulsed (0.45-ms, 1.4-J) laser ͑l 585 nm͒ irradiation ͑7 J cm 22 ͒ of a PWS; a detailed description of the apparatus and irradiation procedure was given previously. 5

RESULTS AND DISCUSSION
We present computational results for the PPTR inversion problems and discuss their interpretation. First we compute the singular values and right singular vectors of the kernel matrix [Eq. (7)]. We compute solution estimates of the PPTR inverse problem, (1) using LSTSVD, (2) unconstrained, and (3) using nonnegatively constrained conjugate gradients. For each inversion method tested no prior knowledge of the initial temperature distribution is assumed; the initial estimate is always DT 0. For each inversion method we use an L-curve analysis to determine the optimum degree of regularization. The error of the computed solution estimates [Eq. (18)] permits an objective comparison of each inversion method.

A. Singular-Value Decomposition
The degree to which the PPTR inverse problem is ill posed is determined by examination of the singular values of the kernel matrix [Eq. (7)]. Solution error depends not only on the degree to which the problem itself (i.e., the kernel function) is ill posed but also on the SNR of the detection system. We plot the magnitude of the first 15 singular values of the kernel matrix corresponding to the previously defined model problem and examine the dependence on the infrared absorption ͑m a ͒ and heat-loss (h) coefficients.
a marginal increase in the number of included singular values and right singular vectors is obtained for an order-of-magnitude increase in SNR. This feature is particularly important when one is considering SNR limitations of infrared detection systems. The best SNR is achieved for a BLIP detection system and is subject to limitations in the frequency bandwidth Df . For example, the SNR of a BLIP liquid-nitrogen-cooled detector system with a 10-kHz bandwidth, using a 1-mm 2 detector that is sensitive in the 7 -11-mm range, does not exceed 100 (Fig. 2); in this case the number of useful singular values is limited to eight or fewer, approximately 6% of the total. Higher SNR (500) can be attained by reduction of the system bandwidth (Df 500 Hz); however, such an increase permits inclusion of only one additional singular value. Third, if the infrared absorption coefficient m a in the test material increases, singular values of the kernel matrix also increase (Fig. 4A), and the inversion problem is slightly less ill posed. We understand this by considering a relatively small infrared absorption coefficient; in this case the penetration depth 1͞m a into the exposed test material is deeper, and the infrared measurement samples the temperature distribution DT ͑z, t͒ over a greater depth. Hence a function [e.g., exp͑2m a z͒ in Eq. (4)] with a larger spread (i.e., smaller m a ) blurs the sampled temperature distribution, which accentuates the ill posedness of the inversion problem. Numerical experiments test the dependence of the singular values on the magnitude of the heat-loss coefficient h. Reduction of the largest singular values is observed when the heat-loss coefficient is increased to values ͓h ഠ 20 mm 21 , or 670 ͑W m 22 K 21 ͔͒ corresponding to forced liquid convective cooling (Fig. 4B). In such circumstances, large heat fluxes force a large temperature gradient near the test material surface so that the measured infrared signal amplitude DS͑t͒ is less indicative of the temperature distribution at deeper positions [Eq. (4)].
We plot (Fig. 5) the first five right singular vectors ͑t i ͒ of the kernel matrix corresponding to the model problem. As noted above, smaller singular values are associated with higher-spatial-frequency right singular vectors. When the initial temperature distribution contains large thermal gradients representative of discrete chromophores, high-spatial-frequency right singular vectors with associated small singular values must be included in the solution for accurate edge reconstruction. Because the SNR limits the number of high-frequency right singular vectors that can be included in a solution estimate, accurate reconstruction of corners and edges that may be inherent in the initial temperature distribution is impossible.

B. LSTSVD Solutions
Solution estimates for the initial temperature distribution DT ͑z, 0͒ by LSTSVD are computed, given a PPTR signal that is representative of the defined model problem (Fig. 3B). We examine the effects of noise in the infrared detection system by considering SNR 100 and SNR1 1000. Proper regularization or selection of the appropriate truncation integer (m) in Eq. (11) is necessary for best results. To determine the optimum truncation integer, we plot the estimate ͑jjDT jj͒ versus the residual ͑jjKDT 2 Sjj͒ Euclidean norms for each SNR (Fig. 6A). Points clustered near the corner of the L curve give a truncation integer; for both SNR 100 and SNR 1000, we find that m 10. Minimum solution error [Eq. (18)] in the computed estimates is 0.80 and 0.62 for respective SNR's (100 and 1000). We plot (Fig. 6B) both solution estimates and the true initial temperature distribution.
When a solution estimate is formed with few right singular vectors, the residual error is large, and the corresponding point lies on the lower right-hand edge of the L curve. As more right singular vectors are included in the solution estimate, solution error e decreases, and points move toward the corner of the L curve where residual and estimate norms are relatively balanced in magnitude. As successive right singular vectors are included, the solution estimate becomes unstable and the norm jjDT jj increases rapidly, with little reduction of the residual norm. In this regime, signal noise magnifies the  influence of the small singular values associated with high-spatial-frequency right singular vectors (Fig. 5). The effect is evident in the L-curve analysis (Fig. 6A); inclusion of 11 or more right singular vectors produces unstable solution estimates with large norms (i.e., jjDT jj . 20). Because the inversion problem is severely ill posed, an order-of-magnitude increase in SNR of the infrared detection system gives a small decrease in computed solution error (i.e., from e 0.80 to e 0.62) despite an order-of-magnitude decrease in the residual norm. Furthermore, reconstruction of corners and edges that are inherent in the initial temperature distribution (Fig. 6B) is difficult because a sufficient number of right singular vectors cannot be included in a solution estimate, given the SNR and bandwidth limitations of the detection system.

C. Conjugate-Gradient Solutions
We determine estimates for the inverse problem [Eq. (5)], using constrained and unconstrained conjugate gradients, by minimizing the regularized norm [Eqs. (13)]. Solution estimates are computed by conventional and early-termination regularization methods. For the conventional case we estimate a solution DT by solving the augmented problem [Eq. (13b)] with L I for various regularization parameters ͑l 0.01 -3 3 10 28 ͒. When using early termination with conjugate gradients, we solve the unaugmented problem [Eq. (13a), l 0] and regularize by selecting the terminating iteration number.
We find characteristic L curves by plotting (Fig. 7A) the estimate versus residual norms corresponding to regu-larization by both conventional and early termination. For each method we find that respective parameters (l 1 3 10 26 , 15th iterate) represent a reasonable balance between signal error and excessive regularization. Conventional and early-termination regularization give the respective solution errors e 0.74 and e 0.75. Both methods of regularization are nearly equivalent and are observed to give solution estimates that are similar to those found when LSTSVD is used. The 15 iterations used in the unconstrained conjugategradient method is approximately equal to the number of right singular vectors (i.e., 10) included in the corresponding LSTSVD solution estimate. When conventional regularization is used, however, each value of l used to compute a single point on the L curve requires a large number of iterations, greatly increasing the total computation time; for example, generation of the L curve for conventional regularization requires in excess of 7500 iterations. Use of early termination simplifies the computation because each successive iteration represents a new and decreased degree of regularization and the optimum point is found with substantially fewer iterations (i.e., 300).
As with the LSTSVD method, however, the unconstrained conjugate-gradient solution estimates exhibit unrealistic negative temperature changes that are not representative of the true solution. Examination of each solution estimate (Figs. 6B and 7B) points to an underlying limitation of both methods. When the initial temperature distribution contains corners and sharp edges representative of large thermal gradients, higher-order right singular vectors must be included in a solution estimate for accurate reconstruction. Inasmuch as an arbitrarily large number of right singular vectors cannot be included because of fundamental SNR and bandwidth limitations, broadened and physically unrealistic negative temperature changes are inevitably present and contribute to increased solution error e.
The problem with negative temperature changes in unconstrained solutions is resolved by introduction of a nonnegativity constraint. We consider the model problem for two infrared detection systems with SNR 10 and SNR 100 and regularize by early termination. The residual and solution estimate norms follow an L curve that we had observed earlier, using the unconstrained methods. The constrained conjugate-gradient solution estimates (Fig. 8) substantially reduce the solution error;   for SNR 10 and SNR 100 the respective values are e 0.65 and e 0.27. Solution error when the nonnegativity constraint is used is reduced by 2.5 times compared with that for the unconstrained method for equivalent SNR 100. When SNR 10, the solution error of the constrained algorithm is comparable with that obtained with the unconstrained method and SNR 100. Furthermore, unrealistic negative temperature changes are not present in solution estimates computed with the nonnegativity constraint, and the initial temperature gradient is reproduced with greater accuracy.
Next we consider a medical diagnostic problem. A PPTR signal ͑DS͒ representative of a laser-irradiated simulated PWS is computed by application of Eq. (15). The nonnegatively constrained conjugate-gradient algorithm is applied to compute a solution estimate of the initial temperature distribution. The terminating iteration number (200th iterate) is identified by use of an L-curve analysis. The position and the amplitude of both the epidermal melanin and PWS blood vessel temperature increases are reproduced (Fig. 9) with a small solution error e 0.35.
Finally, we consider a PPTR signal recorded following pulsed laser irradiation of in vivo PWS skin. Important features in the recorded signal are evident (Fig. 10A): a rapid initial temperature increase that is due to light absorption in the epidermis by melanin and the presence of a delayed thermal wave that is due to heat generated within subsurface PWS blood vessels diffusing to the skin surface. The terminating iteration number (32nd iterate) is determined with an L-curve analysis. Heating of both epidermal melanin and PWS blood vessels is clearly distinguished in the initial temperature distribution computed (Fig. 10B) by the nonnegatively constrained conjugate gradient algorithm. The computed temperature increase in the epidermis ͑35 ± C͒ provides a means to estimate a threshold laser dose (16 J cm 22 ) that corresponds to the temperature of melanosome explosion ͑110 ± C͒. 31 Laser doses above the threshold increase the risk of epidermal injury, long-term changes in pigmentation, and scarring; such a regime is to be avoided in clinical practice.
The depth (200 mm) and the extent (400 mm) of blood vessel heating are also deduced from the computed solution estimate. Inasmuch as blood vessels are discrete absorbers of limited spatial extent, the computed solution represents a mean areal temperature and does not provide sufficient information to permit estimation of the clinically important temperature increase of individual blood vessels. Noninvasive optical methods that allow for accurate estimates of the fractional blood vessel area at a given depth on a site-to-site basis are under investigation.

CONCLUSIONS
The PPTR inversion problem is severely ill posed and inherently difficult to solve. A large infrared absorption m a and/or a small heat-loss coefficient h increases the magnitude of the singular values and thus alleviates the severity of the inverse problem, albeit slightly. The SNR of the background-limited infrared detection system and corresponding bandwidth constraints fundamentally restrict the number of singular values and vectors that can be used to reconstruct the initial temperature distribution. The difficulty is compounded when large thermal gradients representative of discrete chromophores are present. Unconstrained inversion methods are not well adapted to reconstruct edges because -of necessity -SNR and bandwidth limitations restrict the number of included right singular vectors, resulting in physically unrealistic negative temperature changes and increased error in computed solution estimates. Solution error is substantially reduced ͑2.53͒ by inclusion of a nonnegativity constraint. Regularization by early termination permits faster ͑303͒ computation of solution estimates over conventional methods. Inversion of PPTR signals recorded following pulsed laser irradiation of PWS skin provides for direct estimates of the magnitude of melanin heating and the limiting laser dose as well as the depth and the spatial extent of laser-heated blood vessels.