Accuracy of subsurface temperature distributions computed from pulsed photothermal radiometry

Pulsed photothermal radiometry (PPTR) is a non-contact method for determining the temperature increase in subsurface chromophore layers immediately following pulsed laser irradiation. In this paper the inherent limitations of PPTR are identified. A time record of infrared emission from a test material due to laser heating of a subsurface chromophore layer is calculated and used as input data for a non-negatively constrained conjugate gradient algorithm. Position and magnitude of temperature increase in a model chromophore layer immediately following pulsed laser irradiation are computed. Differences between simulated and computed temperature increase are reported as a function of thickness, depth and signal-to-noise ratio (SNR). The average depth of the chromophore layer and integral of temperature increase in the test material are accurately predicted by the algorithm. When the thickness/depth ratio is less than 25%, the computed peak temperature increase is always significantly less than the true value. Moreover, the computed thickness of the chromophore layer is much larger than the true value. The accuracy of the computed subsurface temperature distribution is investigated with the singular value decomposition of the kernel matrix. The relatively small number of right singular vectors that may be used (8% of the rank of the kernel matrix) to represent the simulated temperature increase in the test material limits the accuracy of PPTR. We show that relative error between simulated and computed temperature increase is essentially constant for a particular thickness/depth ratio.


Introduction
Pulsed photothermal radiometry (PPTR) is a non-contact method for obtaining information on subsurface chromophores in a test material. A fast infrared detector is used to measure increase in infrared emission at the test material surface following pulsed laser irradiation (Milner et al 1995a, b). From the time record of infrared emission increase, estimates of thickness and depth of subsurface chromophores can be computed. Reported applications of PPTR include the evaluation of surface coating thickness in industrial components (Crostack et al 1989), the identification of subsurface microcracks in aircraft structures (Favro et al 1993), the determination of the optical absorption coefficients in human arteries  and biliary calculi  and the characterization of port wine Also at: Physics and Astronomy Department, University of Canterbury, Christchurch, New Zealand. ¶ Author to whom correspondence should be addressed. 0031-9155/98/092453+11$19.50 c 1998 IOP Publishing Ltd stain birthmarks (Jacques et al 1993). Analysis of the performance of PPTR to accurately compute the size and position of blood vessels in a port wine stain (PWS) is of particular interest. Efficacy of the pulsed laser treatment of PWS can be improved if the size and position of blood vessels is known (Kimel et al 1994).
While several authors have reported on the use of PPTR (Jacques et al 1993, Vitkin et al 1994, Milner et al 1995a, b, 1996a, theoretical limits on the accuracy of the estimated position and magnitude of the temperature increase in discrete laser-heated subsurface chromophores have not been proposed or established. In this paper, we investigate the fundamental limitations of PPTR.
A non-negatively constrained conjugate gradient algorithm is used to compute the position and magnitude of temperature increase within a test material from the time record of infrared emission following pulsed laser irradiation.
Results from the non-negatively constrained conjugate gradient algorithm are verified using a singular value decomposition (SVD) analysis of the PPTR kernel function K(t, z). Since each model temperature increase may be resolved as a superposition of variable number of right singular vectors of the PPTR kernel matrix, limitations on the computed thickness, depth and magnitude of temperature increase are established for various signalto-noise ratios (SNR).

Method
The skin model (figure 1) contains a single chromophore layer. The laser spot diameter at the model surface is assumed to be large relative to the thermal and optical diffusion lengths so that heat transport is only considered along the depth (z) axis. Immediately following pulsed laser irradiation (t = 0), the simulated temperature increase within the skin model is given by T s (z) ( • C). The relationship between increase in PPTR signal, S ( • C), and T s is given by a Fredholm integral of the first kind (equation (1)) (Milner et al 1995a) where time is represented by t (s) and distance into the tissue, measured normal to the surface, is given by z (mm). Including effects of infrared emission and a convective boundary condition, the kernel function K(t, z) (Milner et al 1995a), is where erfcx(u) = exp(u 2 ) × erfc(u) is the exponential complementary error function. Functions µ ± (t, z) and µ 1 (t, z) are defined as Signal noise, n(t), an inherent component of the system is taken as a zero-mean, white, Gaussian distributed function. The signal-to-noise variance n 2 is related to the SNR of the infrared detection system SNR = S(t) / n 2 1/2 (5) where x represents a time average of the quantity x. Physical constants used in the skin model (from Robertson and Williams 1971, Duck 1990, Incropera and DeWitt 1985 are: thermal diffusivity D (1.1 × 10 −2 mm 2 s −1 ), infrared absorption coefficient µ ir (45 × 10 −2 mm −1 ), and free convection heat-transfer coefficient, h (0.03 mm −1 ). A tophat shaped simulated temperature increase distribution (equation (6)) was used to represent the temperature change within tissue. Such a temperature increase distribution is difficult to reconstruct since it contains abrupt changes in value. Consequently this temperature distribution provides an excellent test for the reconstruction algorithm and any deficiencies will become apparent T st (z) = T 0 z z m − r and z z m + r 0 otherwise.
The PPTR signal from the model, ( S(t)), was determined from where and derived from equation (1). Computation of temperature increase in the skin model, given the PPTR signal as input data, is a severely ill-posed inverse problem. The problem is converted to a matrix equation and vector quantities representing S(t) and T (z) are denoted by boldface symbols S and T . A non-negatively constrained conjugate gradient algorithm (see, for example, Goodman et al 1993) was used to compute the temperature increase within skin at t = 0. Briefly, this iterative method computes an ith estimate for the temperature increase within skin ( T cg (i)) (initially a null field for i = 0) and the corresponding PPTR signal. T cg (i) is compared with S, and a revised estimate of ( T cg (i)) is computed based upon conjugate gradient direction. As temperature increase within skin is always positive, the solution is constrained so that negative values are not allowed ( T cg (i) 0).
The iterative process is terminated when the difference between the computed PPTR signal and S is less than the noise level (Groetsch 1984) or the difference ceases to change between successive iterations, or a preset number of iterations have been exceeded (Goodman et al 1993, Frank and Friedman 1993, Stone and Brooks 1990.
In our simulations, 256 uniformly spaced nodes were used to represent temperature increase in a 1.5 mm tissue depth. PPTR signal over a 3 s period immediately following pulsed laser irradiation was stored in a 512-element array.

Results
Computed temperature increase ( T cg ) was calculated (figure 2) after 100 iterations of the non-negatively constrained conjugate gradient algorithm for z m = 0.5 mm and thicknesses (2r) between 0.032 and 0.5 mm. Error in the computed peak temperature increase was less than 15% for chromophore layers of thicknesses 0.25 and 0.50 mm and increased to 77% for the thinnest chromophore layer (2r = 0.032 mm). Note, however, that the depth z m and T cg dz are still predicted accurately, even when the computed peak temperature increase is much less than T 0 .
In each of the following simulations results were obtained after 100 iterations of the nonnegatively constrained conjugate gradient algorithm, and regularized by early termination (Milner et al 1995a). For one particular T ts , repeated solutions of equation (1) yields a slight variation in the computed temperature increase ( T cg ), because S contains a noise component, n. Thus, 20 inversions for each simulated temperature increase were computed and subsequent plots show mean and standard deviation of the displayed values.
The effect of source chromophore layer thickness is shown for z m = 0.2 mm (figure 3). The computed peak temperature increase to T 0 is reported for three values of T 0 (5, 10 and 40 • C), and the chromophore layer thickness 2r is between 0.025 and 0.30 mm. When the chromophore layer thickness is greater than 0.1 mm, the computed peak temperature increase differs from T 0 by less than 20%. As the chromophore layer thickness increases, the difference between the computed peak temperature increase and T 0 becomes less. For values of 2r < 0.1 mm, the difference increases as the thickness of the source chromophore layer decreases. The thickness/depth ratio has a significant effect on computed temperature increase (figure 4). Chromophore layer depths used were 0.1, 0.2, 0.5 and 1.0 mm. For a thickness/depth ratio of 25%, the computed peak temperature increase is always less than T 0 . The difference is significant for a shallow chromophore layer, and especially so when T 0 is less than 25 • C. For a thickness/depth ratio of 50%, the difference between computed peak temperature increase and T 0 is reduced ( figure 4). Only when T 0 is 5 • C or less, and the chromophore layer thickness is 0.025 mm, is the computed peak temperature increase less than 60% of T 0 .
For a chromophore layer thickness of 0.025 mm, centred at a depth of 0.1 mm, influence of SNR is significant. The computed peak temperature increase is greater than 65% of T 0 for SNR ratios of 100, 200 and 1000 (figure 5). When T 0 exceeds 25 • C, SNR becomes irrelevant. For an SNR of 1000, the computed peak temperature increase is essentially unaffected by the source temperature T 0 .
The difference between the computed peak temperature increase and T 0 is greatest when T 0 is less than 25 • C, and SNR is less than 1000. At the extreme, the computed peak temperature increase/T 0 ratio is 12% when T 0 is 2.5 • C and SNR is 100.
In our calculations, the mean (first moment) of the deduced temperature increase distribution and T cg dz were predicted accurately (±10%), except for cases of: (i) thickness/depth ratio less than 10%, (ii) SNR below 10, or (iii) T 0 less than 5 • C. In these cases, the PPTR signal from the model ( S) is very small, and cannot be distinguished by the non-negatively constrained conjugate gradient algorithm from the background radiation level. Figure 4. Ratio of computed peak temperature increase and T 0 , for a chromophore layer at a depth of 0.1 mm (--), 0.2 mm (---), 0.5 mm (----) and 1.0 mm (· · · · · ·). Thickness/depth ratios used: (a) 25% and (b) 50%. In (b), because computed peak temperature increases did not change appreciably for values of T 0 greater than 20 • C, these values are not displayed.

Discussion
It is to be expected that the non-negatively constrained conjugate gradient algorithm will perform poorly in cases where T 0 is low and SNR is high, since the PPTR signal from the model chromophore is obscured by the background radiation. To show why the thickness/depth ratio has a significant influence on the computed result, we use the singular value decomposition (SVD) of the kernel matrix, K. Thus where σ i and τ i are orthonormal left and right singular vectors of K and form a basis for the respective vector spaces containing the PPTR signal ( S) and temperature increase ( T ); δ i is the ith singular value of the PPTR kernel matrix K, where δ i δ j and i < j. As {τ i } forms an orthonormal basis for the vector space spanning the simulated temperature increase ( T s ), a solution estimate ( T svd (k)) may be expanded in terms of a variable number (k) of right singular vectors (equation (10)) Since the singular values, δ i , rapidly decrease with increasing index i, (δ i ≈ A −i ) where A is constant, the PPTR inversion problem is severely ill-posed and it is physically meaningless to include a large number of terms in the solution estimate (equation (10)). The maximum number of singular vectors in the solution estimate is determined by the SNR of the infrared detection system, which is limited by background photon noise. From Milner et al (1995a), when SNR is between 100 and 1000, k is ≈5-10, which is less than 8% of the rank of the kernel matrix. An example of the effect of thickness/depth ratio on the computed temperature increase is given in figure 6, where T svd (k) for various k are indicated. Solution estimates ( T svd ) are computed for superficial and deep chromophore layers of varying thickness. In all cases, greater accuracy of the solution estimate is obtained by increasing the number (k) of right singular vectors (τ i ) included. Accuracy of the solution estimate represented by a thin (0.1 mm) chromophore layer is best at superficial depths; inclusion of seven right singular values is sufficient for a reasonably accurate representation of the simulated temperature increase ( figure 6(a)). In comparison, when the same chromophore layer is positioned 1 mm below the skin model surface and consequently the thickness/depth ratio is reduced, inclusion of 17 right singular vectors results in an estimate with significant error ( figure 6(b)). Alternatively, when the thickness of the same chromophore layer is increased to 0.5 mm, the relative error of the solution estimate is improved significantly (figure 6(c)). Figure 6 illustrates a trend observed: when the thickness/depth ratio is less than 25% the computed peak temperature increase is always significantly less than T 0 . Also, the thickness of the computed temperature increase is always larger than T ts . Figure 6. T svd calculated from a simulated temperature increase, where the chromophore layer has dimensions: (a) 0.1 mm thickness, depth 0.2 mm, (b) 0.1 mm thickness, depth 1.0 mm and (c) 0.5 mm thickness, depth 1.0 mm. Simulated temperature increase (--), 7 vectors (---), 17 vectors (-· · -) and 37 vectors (· · · · · ·).
Next we compute the relative error ( (k), equation (11)) of the solution estimate (T svd ), corresponding to a family of subsurface chromophore layers of varying thickness (2r) and depth (z m ) In figure 7 (k) is depicted in a contour plot for various tophat shaped simulated temperature increase functions ( T st ). Values of 2r and z m which give rise to temperature increase outside of the 0-1.5 mm depth range are not considered. Individual contour lines show a small oscillation that originates from the relatively few singular vectors used in the expansion (10) when calculating the relative error. A similar effect is observed when computing a Fourier expansion of an arbitrary function using relatively few sinusoidal terms. When the thickness/depth ratio is greater than 25%, the relative error ( (k)) is essentially constant. In a PWS, the temperature distribution following pulsed laser irradiation decays approximately as an exponential with depth, commensurate with decreased fluence levels within blood (van Gemert et al 1986a). To simulate this, we describe temperature increase due to laser heated PWS blood vessels by an exponential distribution Here, L 0 is the exponential decay length, d is the minimum depth at which blood is found, T 0 is the peak temperature at z = d. Values for exponential decay length L 0 and depths d which introduced a discontinuity at z = 1.5 mm (i.e. T es (z = 1.5 mm) > 0.1T 0 ) and significantly increased the relative error, (k), were not considered. The threedimensional reconstruction of one PWS (Smithies et al 1997) indicates that the exponential temperature distribution is a reasonable approximation. There was a high concentration of vessels close to the surface, with the blood volume fraction approximately inversely proportional to depth. Following laser irradiation, the average temperature increase at each depth point will be similar to that described by equation (12). In figure 8 (k) is depicted in a contour plot for various exponential shaped simulated temperature increase functions ( T es ). Individual contour lines show a small oscillation that originates from the relatively few singular vectors used in the expansion (10) when calculating the relative error. A similar effect is observed in figure 7. Relative error ( (k)) is essentially constant for a given exponential decay length/depth ratio.

Conclusion
Our results indicate the circumstances when PPTR may be used to accurately predict the thickness, depth and magnitude of temperature increase in the chromophore layer contained in a test material.
The mean (first moment) of the deduced temperature distribution and T cg (z) dz were predicted accurately (±10%), except for cases of (i) thickness/depth ratio less than 10%, (ii) SNR below 10, or (iii) T 0 less than 5 • C.
When thickness/depth ratio is less than 25%, the difference between simulated and computed temperature increase is sufficiently large that the utility of the calculation is minimal. For these cases, the computed peak temperature increase is always significantly less than T 0 , and the computed thickness of the chromophore layer is significantly larger, so that their product remains constant. Conversely, with a thickness/depth ratio greater than 25%, the computed peak temperature increase and chromophore layer thickness is much closer to the true value. Thus, PPTR alone cannot distinguish between a thin deep chromophore layer that reaches a high temperature and a relatively cool, thick and deep chromophore layer.
From the singular value decomposition of the kernel matrix, we verify results computed by application of the non-negatively constrained conjugate gradient algorithm, that the relative error is essentially constant for a particular thickness/depth ratio. The relatively small number of right singular vectors that may be used limits the accuracy of the computed temperature increase.