Spectral variation of the infrared absorption coefficient in pulsed photothermal profiling of biological samples

Pulsed photothermal radiometry can be used for non-invasive depth profiling of optically scattering samples, including biological tissues such as human skin. Computational reconstruction of the laser-induced temperature profile from recorded radiometric signals is sensitive to the value of the tissue absorption coefficient in the infrared detection band (μIR). While assumed constant in reported reconstruction algorithms, μIR of human skin varies by two orders of magnitude in the commonly used 3−5 μm detection band. We analyse the problem of selecting the effective absorption coefficient value to be used with such algorithms. In a numerical simulation of photothermal profiling we demonstrate that results can be markedly impaired, unless the reconstruction algorithm is augmented by accounting for spectral variation μIR(λ). Alternatively, narrowing the detection band to 4.5–5 μm reduces the spectral variation μIR(λ) to a level that permits the use of the simpler, un-augmented algorithm. Implementation of the latter approach for depth profiling of port wine stain birthmarks in vivo is presented and discussed.


Introduction
Pulsed photothermal radiometry (PPTR) is based on the time-resolved acquisition of infrared (IR) emission from a sample after pulsed laser exposure. When optical and thermal properties of the sample are known, the laser-induced temperature profile can be reconstructed from measured radiometric signals, providing information on chromophore distribution in the sampled volume (Tam andSullivan 1983, Imhof et al 1984). Such a PPTR technique was applied to depth profiling of strongly scattering biological tissues and tissue phantoms (Long et al 1987, Bindra et al 1994, Vitkin et al 1995, Prahl 1997, including port wine stain (PWS) birthmarks in human skin (Jacques et al 1993, Milner et al 1995, 1996. Furthermore, three-dimensional imaging of embedded structures by recording radiative emission from the surface with high spatial resolution has been demonstrated (Telenkov et al 1998). PWS are subsurface hypervascular lesions, located within the most superficial millimetre of human skin. Depth distribution of PWS blood vessels varies from patient to patient, with the highest fractional blood content on average 0.2−0.4 mm below the epidermal-dermal junction (Barsky et al 1980). PWS therapy currently utilizes selective photocoagulation of the ectatic vasculature using pulsed green or yellow lasers. Recent studies have indicated that such therapy could benefit substantially from the optimization on an individual patient basis (Kimel et al 1994), in particular when involving cryogen spray cooling , Verkruysse et al 2000.
For the purpose of PWS characterization, existing non-invasive techniques exhibit various disadvantages, including insufficient spatial resolution (ultrasound), poor contrast (optical coherence tomography) and prohibitive cost (magnetic resonance imaging). The recently developed optical Doppler tomography technique (Zhao et al 2000) compares favourably with the above techniques, enabling three-dimensional imaging of in vivo PWS (Zhao et al 2001). However, while the above approaches may provide reliable structural information, PPTR also measures epidermal temperature rise following a diagnostic laser pulse. This enables assessment of the highest permissible radiant exposure (varying strongly with individual patient's skin type), which is a key to the success and safety of PWS therapy. The same information could be also obtained by photoacoustic profiling, which has been recently applied to the same problem (Viator et al 2002).
Analysis of the recorded radiometric signals involves the tissue absorption coefficient at detected wavelength, µ IR . It is therefore important to note that µ IR of water-the primary absorber in the mid-IR spectral region in many biological tissues, including human skin-varies by two orders of magnitude, between 1070 and 11.8 mm −1 at wavelengths of 3.00 and 3.82 µm, respectively (figure 1) (Palik 1997).
The influence of the spectrally varying value µ IR on PPTR signals, as well as other wavelength-related effects, has been investigated in detail in the past (e.g., Prahl et al 1992. Milner et al (1995) used the singular-value decomposition approach to the inverse problem of PPTR depth profiling, and showed that it is slightly less ill-posed at high µ IR values. Interference filters centred at water absorption peaks (6.05 and 13 µm) were applied by Bindra et al (1994) in profiling of in vivo human skin.
Almost generally, however, broadband detection is utilized to enhance the signal-to-noise ratio (SNR) in PPTR signals. When studying water-based tissue phantoms, Prahl et al (1992) calculated the effective absorption coefficient value by averaging µ IR (λ), weighted by detector sensitivity, D(λ). Milner et al (1995) used both D(λ) and Planck's law of radiation as weights in the calculation of the effective µ IR in similar systems.
Recently, Zuerlein et al (1999) accounted explicitly for the (approximately known) spectral dependence µ IR (λ) in PPTR determination of optical properties in dental enamel. All reported PPTR profiling and imaging experiments on biological systems, however, use the approximation of a fixed µ IR . The main purpose of the present work is, therefore, to investigate the influence of strongly varying µ IR (λ) within a broad detection band on the inverse problem of PPTR profiling.
In particular, we analyse the effects of spectral variation µ IR (λ) on PPTR profiling of human skin using a 3-5 µm radiation detector (InSb). First, we address the problem of selecting the effective absorption coefficient (µ eff ) to be used with image reconstruction algorithms that require a fixed value. Next, we introduce an augmented reconstruction algorithm, which accounts for the spectral variation of µ IR as well as detector response, and compare the performance of the two algorithms in a numerical simulation (in section 3).
Alternatively, our theoretical analysis suggests that narrowing the detection range to 4.5-5.0 µm permits the use of the usual, un-augmented reconstruction algorithm. An experimental implementation of such approach for depth profiling of PWS in vivo is presented and discussed in section 4. The experiment involves the recently introduced dual-wavelength excitation (DWE) technique, which explores spectral differences between melanin and haemoglobin to separate the epidermal and vascular contributions to the radiometric signals (Majaron et al 2000a(Majaron et al , 2000b. In this way, PWS depth can be assessed even when it lies in close proximity to the epidermal-dermal junction, despite the inherently limited spatial resolution of the PPTR technique (Milner et al 1995, Sathyam and Prahl 1997.

PPTR signal
Detailed discussion of PPTR profiling in a semi-infinite medium can be found elsewhere (Milner et al 1995, Prahl 1997. In short, one-dimensional analysis is justified because the laser-irradiated spot in our experiments is much larger (>5 mm in diameter) than the involved optical penetration depths and thermal diffusion lengths. Owing to the finite penetration depth of emitted IR radiation µ IR −1 , a radiometric signal from a non-uniformly heated sample is not representative of its surface temperature, but rather is composed of correspondingly attenuated contributions from all depths z. Finally, Planck's law of radiation, which describes the locally emitted power density at wavelength λ (in W m −2 per unit of wavelength) is linearized around the initial skin temperature T 0 : B λ (T ) ≈ B λ (T 0 ) + B λ (T 0 ) T . As a result, the radiometric signal S λ (t) is expressed as where constant C accounts for optical properties of the sample surface, numerical aperture and losses of the collection optics, and detector specifics (such as size, sensitivity, integration time and spectral bandwidth). For the purpose of PPTR, the temperature field evolution in the sample, T (z, t), is determined uniquely by the laser-induced temperature rise, T (z, t = 0), and the thermal point-spread (Green's) function, K T (z , z, t): In the following, we apply K T (z , z, t) as derived by Milner et al (1995) for a semi-infinite medium with a convective/radiative heat exchange at the surface Dt. D marks the thermal diffusivity of skin (D = 1.1 × 10 −7 m 2 s −1 ; Duck 1990) and h = 0.02 mm −1 is the reduced heat transfer coefficient, defined as the ratio between the customary heat transfer coefficient at the interface (in units of W m −2 K −1 ) (e.g. Carslaw and Jaeger 1959), and thermal conductivity of skin.
By inserting (2) into (1) and re-ordering, the dynamic part of the radiometric signal is Since K T (z , z, t) (3) is known explicitly, integration over z (which marks a general tissue depth from which the signal is collected, while z is the depth coordinate in the initial temperature profile) can be performed, yielding The PPTR signal S λ (t) is related to the initial temperature profile T (z, 0) in a linear manner, with definition of the corresponding kernel function K(z, t) obvious from equations (4) and (5).

Spectral variation of µ IR
Although the radiometric signal expression (4) is valid rigorously only for the case of monochromatic radiation detection, it is adopted quite generally in PPTR studies to reduce the mathematical complexity. To analyse possible deficiencies of such an approximation in the presence of spectral variation µ IR (λ), we calculate the radiometric signal by summing spectral components within the detection band (λ l to λ h ): At the same time, we take into account the detector spectral sensitivity, D(λ). (Note the re-defined constant, C .) As an example, let us consider PPTR profiling of human skin using an InSb detector. For µ IR (λ), we take 70% of the values for water (figure 1, solid line), and disregard the small protein contribution (Yannas 1972). As is common to quantum radiation detectors, the spectral sensitivity of InSb is proportional to λ up to the long-wavelength cut-off around 5 µm. At the short-wavelength end, the detection band is usually limited by a long-pass filter  at 3 µm. Although detailed spectral characteristics of individual detectors may vary, we will demonstrate below that the overall influence of D(λ) on PPTR signals is much smaller than the effect of µ IR (λ). For the purpose of our analysis, we thus set D(λ) = λ µm −1 between λ l = 3 µm and λ h = 5 µm, and D(λ) = 0 elsewhere (figure 1, dashed line). (Evaluation of the radiometric signal amplitude in absolute physical units will not be required in this study.) The spectral dependence of B λ (T 0 ) at T 0 = 310 K (36.8 • C) is plotted in the same figure for completeness (dotted line).
Interchanging the integration order in (6) yields In (7b), we have introduced an attenuation function, R(z), which represents the amplitude of the spectrally composite radiometric signal, originating from a unity temperature increase in a delta layer at varying depths z. Figure 2 presents such an attenuation function R(z), calculated for the assumed spectral dependencies µ IR (λ) and D(λ) (solid line). Its functional dependence differs from the simple exponential decrease, which is predicted when neglecting the spectral variation of µ IR : by introducing a spectrally constant value µ eff into (7a), integration over λ can be taken out of the integral over z, yielding an expression essentially equivalent to the monochromatic approximation (1) Two such exponential attenuation functions, corresponding to µ eff values of 22 and 25 mm −1 , are plotted in figure 2 for comparison (dashed and dotted lines, respectively).
To assess which value µ eff is optimal to use in PPTR algorithms that disregard the spectral variation µ IR (λ), we equate the radiometric signal expression (7a) with the approximate  formula (8), which results in an implicit equation A numerical solution of (9) for the discussed spectral dependencies is presented in figure 3 (solid lines). The optimal absorption coefficient value µ eff varies strongly with the depth of radiation origin (z). This indicates that, under discussed conditions, using (8) with a constant value µ eff will result in inaccurate PPTR signal predictions. Figure 3 also presents the solution of (9) for a hypothetical case of spectrally flat detector response, D(λ) = constant (dotted line). The similarity with the previous solution demonstrates that the deviation of attenuation function R(z) from simple exponential dependence in figure 2 originates mainly from the spectral variation of µ IR (λ), not detector sensitivity, D(λ). In contrast, a practically constant value (µ eff ≈ 26.5 mm −1 ) is obtained with D(λ) as described above, if the detection range is reduced to 4.5-5 µm. Such narrowing of the detection band thus reduces the spectral variation µ IR (λ) to a level which permits use of the approximate expression (8) with a constant value µ eff .
The above analysis suggests that special care must be taken with regard to spectral variation of µ IR (λ) when applying PPTR techniques to biological samples. In the following section, we analyse the effects of the customary, approximate expression (8) in a numerical simulation of PPTR profiling.

Numerical simulation
After defining the initial temperature distribution T (z, 0) ('object'), a spectrally composite PPTR signal S(t) (7a) is calculated using the assumed spectral variations µ IR (λ) and D(λ) within the 3-5 µm detection band. To simulate experimental conditions, the signal is presented as vector S with 250 equidistant components (S i = S(t i ), t = 2 ms), to which a certain amount of normally distributed white noise was added (SNR = 300).
The linear relation (5) is thereby represented by a straightforward multiplication with kernel matrix K: Spectrally composite PPTR signals were calculated using (7a) for a hyper-Gaussian initial temperature profile centred at depth z = 300 µm (dashed line).
Reconstruction of the initial temperature profile T (T j = T (z j , 0)) from a given PPTR signal S, on the other hand, is a severely ill-posed inverse problem (Milner et al 1995). We solve this problem by applying iterative minimization of the quadratic norm of the residue, ε = S − KT (n) , i.e. the difference between the signal vector S and that predicted from the nth iterative solution T (n) ('image'). The minimization is performed using a non-negatively constrained conjugate-gradient algorithm, starting from T (0) = 0, which has been previously identified as an efficient and reliable approach to this problem (Milner et al 1995).
In the presented examples, the depth discretization is z = 10 µm, requiring the optimization of 100 components T j to cover a 1 mm thick superficial layer of skin. Since the kernel matrix elements K i,j (10) do not vary during the iteration run, and do not depend on the object specifics, they must be calculated only once for each kernel function K(z, t) and values t and z. The computational load is thereby significantly reduced, which enables us to perform image reconstructions on a personal computer (AMD Thunderbird 900). We set up a computation of the residual norm (ε) in Microsoft Excel, and utilize the built-in conjugate-gradient algorithm ('Solver' routine) to minimize its value. The applied algorithm settings are non-negativity constraints, tangent estimates and forward derivatives.

Test of the improved kernel function
Possible effects of disregarded spectral variation µ IR (λ) are analysed in a direct comparison of two image reconstructions for each composite PPTR signal S. In the first reconstruction, the kernel matrix elements K i,j are calculated using the approximate expression (8) with a constant µ IR = 22 mm −1 . In the second, the augmented formula (7a), accounting for the spectral variation µ IR (λ) is used to construct the matrix.
Results from one such comparison are presented in figure 4. The object is defined by a hyper-Gaussian function, which resembles the temperature profile of a laser-heated subsurface PWS layer: . Simulated multi-spectral PPTR signal (SNR = 300), used as input in image reconstructions in figure 4 (solid line), and the signal predicted after 100 steps of the approximate reconstruction algorithm using µ IR = 22 mm −1 (bold solid line). The inset shows expanded view of the residue (i.e. difference between these two signals). Radiometric signal calculated from the actual object using the approximate expression (8) with µ IR = 22 mm −1 is plotted for comparison (short dash). w = 0.1 mm (dashed curves). When using the approximate kernel matrix (figure 4(a)), the solution after 50 iteration steps (top graph; solid line) is a blurred image of the actual temperature profile, with a substantial deep artefact. By increasing the number of iteration steps, the image converges slowly towards the actual object. At n = 100 (figure 4(a), middle), when image improvement begins to stall, the peak temperature reaches only 73% of the actual value (and 77% at n = 200; bottom), and a weak artefactual temperature rise starts to appear near the skin surface (see the arrow). The shallower edge of the temperature profile, as assessed at half its maximum height, is indicated as being ∼40 µm too shallow.
When using the improved kernel matrix, the solution converges towards a much more representative image of the actual temperature profile. After 50 iteration steps (figure 4(b), top), the image is already equivalent, if not superior, to the final result of the former algorithm. At n = 100, the depth of PWS top boundary is reproduced within 25 µm (and 15 µm at n = 200), and the peak temperature is within 10% (3%) of the actual value (figure 4(b), middle and bottom). But, the difference between these two approaches is not limited to the convergence rate. Figure 5 presents the radiometric signal as predicted after 100 iteration steps of the former, approximate reconstruction algorithm, i.e. KT (100) (bold solid line). Considering the ill-posedness of the image reconstruction problem, this signal is a rather poor match to the simulated input signal S (thinner, noisy line). The residue, presented in the inset, displays features in the order of 2% of the overall signal amplitude-well above the noise level (∼0.3%).
To better understand the origin of such a discrepancy, consider the PPTR signal calculated from the true object using the approximate expression (8) ( figure 5, dashed line). This is the signal shape towards which the predicted signal would ideally converge when the approximate kernel function is used in the reconstruction. Since the composite signal provided as input is significantly different (thinner solid line), in particular at early times (t < 0.1 ms), even a prolonged iteration with this kernel matrix cannot improve the image.
In contrast, when the augmented kernel function is used in the reconstruction, the residual norm after 100 iteration steps (ε = 4.1 × 10 −6 ) is 3.4 times lower than with the approximate kernel function. By n = 200, it drops to ε = 1.6 × 10 −7 and the residual becomes indistinguishable from signal noise. Figure 6 presents a similar comparison of the two image reconstruction approaches, with an object temperature profile resembling a laser-heated epidermis (∼0.1 mm thick superficial layer of skin): T (z, 0) = T 0 exp −2(z − z 0 ) 4 w 4 , T 0 = 10 • C, z 0 = 0, w = 0.1 mm (dashed lines). When using the approximate kernel function, based on (8), the image is a severely distorted, steeply falling and almost linear temperature profile. It overshoots the surface temperature by 8.5% (figure 6(a), n = 100) and makes assessment of the epidermal thickness unreliable. The solution obtained by using the improved kernel matrix is a closer match to the actual temperature profile ( figure 6(b)). At 100 iteration steps, it overshoots the highest temperature by no more than 5.5% (3.5% at the surface). The residual norm reaches ε = 6.6 × 10 −6 and continues to decrease.
Under experimental conditions involving higher noise levels, such long iterations may be prevented by computational instabilities due to ill-posedness of the image reconstruction problem. Consequently, the benefit of the augmented approach may not be as evident as in the presented examples. Nevertheless, we are confident that disregarding the spectral variation µ IR (λ) can adversely affect PPTR profiling in several realistic scenarios. One such example is presented in the following section.

Depth profiling of PWS birthmarks in vivo
A PWS lesion on a volunteer patient was irradiated with a 1.5 ms pulse from a flashlamppumped dye laser (ScleroPLUS, Candela, Wayland, MA), which allowed selection of wavelength from 585 to 600 nm in increments of 5 nm. Incident radiant exposure at the centre of a 6 mm diameter laser spot was kept at a sub-therapeutic level (∼5 J cm −2 ).
The resulting transient increase in IR radiant emission from the central 1.9 × 1.9 mm 2 area was recorded with an InSb focal plane array camera (Galileo, Raytheon, Dallas, TX) at a rate of 500 frames per second. Due to the 3 µm cut-on filter at the cold stop, the camera was sensitive to the wavelength range of ∼3-5 µm. PPTR signals S(t) were prepared by calibrating the system response with a computer-controlled blackbody (BB701, Omega Engineering, Stamford, CT), averaging over the used subset of 64 × 64 array elements, and subtracting the background signal level.
The initial temperature profiles T (z, 0) were reconstructed using a dedicated computer code running a non-negatively constrained conjugate-gradient iterative algorithm and regularized by early termination (Milner et al 1995). Calculating the kernel matrix elements K i,j anew in each reconstruction run, the presented examples involving a signal vector S with 400 elements and a 64-component image vector T , required 5-10 s CPU time on a Sun Sparc II workstation.

Narrowing the detection band
As discussed in section 2.2, an alternative approach to using the improved,more complex kernel function in the image reconstruction algorithm is a suitable narrowing of the detection band, permitting use of the un-augmented signal expression (4). Based on calculations of µ eff (z) from (9) for several wavelength ranges, and the spectral dependence of B λ (T 0 ) (figure 1), we have decided to narrow the detection band of our system to 4.5-5 µm (see figure 3). A custom long-pass filter with a cut-on wavelength of 4.5 µm (Barr Associates, Westford, MA) was fitted to the collection lens of the IR camera, as the customary installation at the camera cold stop would be impractical for direct comparisons of PPTR profiling using filtered and unfiltered signal acquisitions. Basic tests of the modified system, including irradiating the filter with the excitation laser pulse, revealed no adverse influence of the filter on camera performance.
PPTR signals from the same spot on the patient's PWS lesion were acquired using the full 3-5 µm detection bandwidth of the IR camera and, subsequently, using the cut-on filter to narrow the detection band. The two signals, obtained following 585 nm irradiation with the same radiant exposure, are presented in insets in figures 7(a) and (b) respectively (the amplitude scale is in radiometric temperature degrees). The camera integration time was increased from 0.5 to 0.7 ms for the filtered acquisition, to compensate for the reduced amount of collected IR radiation. In accordance with the theoretical estimations, the raw signal reduction was around 50%.
The presented signals were used as input for the image reconstruction algorithm based on the customary formula with a constant µ IR (4). In accordance with section 2.2, the value µ IR = 22 mm −1 was used to reconstruct the initial temperature profile from the full-bandwidth signal (figure 7(a)), and µ IR = 26 mm −1 for the filtered signal ( figure 7(b)). For each experimental approach, we present a family of results, from blurred images of the temperature profile in under-iterated solutions (n = 2, dashed lines), to near-optimal (solid), and over-iterated solutions (dotted lines). The latter can be recognized by characteristic oscillations or other nonphysical artefacts, developed as the reconstruction algorithm starts amplifying experimental noise and other mismatches between the actual and theoretically predicted signals (Milner et al 1995). Such oscillations are more pronounced in figure 7(b), which can be attributed to the lower SNR in the filtered PPTR signal (see the inset).
Despite the poor convergence of the reconstruction algorithm, it is evident that the narrowing of the detection band has modified the obtained temperature profile. In comparison with the latter, the image obtained from the full-bandwidth signals (figure 7(a)) displays a reduced peak temperature of the deeper (PWS) and increased amplitude of the superficial (epidermal) part of the temperature profile-qualitatively the same effects as observed when disregarding µ IR (λ) in the numerical simulation (section 3.1).

Dual-wavelength excitation
As seen in figure 7, the PWS temperature profile cannot be resolved from the epidermal temperature rise when they are in close proximity, or even overlapping. we can determine the PWS depth and epidermal thickness by implementing the dualwavelength excitation (DWE) technique (Majaron et al 2000a). In this approach, we combine irradiation at 585 and 600 nm wavelengths, and explore the spectral differences between the haemoglobin in PWS blood and epidermal melanin.
Due to the linearity of relation (5), the PPTR signal following 585 nm irradiation is a sum of two contributions, originating from the PWS layer [x(t)] and epidermis [y(t)]: S 585 (t) = x(t) + y(t). In comparison with 585 nm, the haemoglobin absorption at 600 nm is several times weaker (van Kampen and Zijlstra 1965), while the absorption and scattering properties of the epidermis and dermis are essentially equivalent at either wavelength (Wan et al 1981). As a result, the PPTR signal following 600 nm irradiation features a smaller PWS contribution, while the epidermal signal component is nearly the same as at 585 nm, which can be approximated to the first order by where α < 1 and β ≈ 1. The PWS and epidermal contributions to S 585 (t) are thereby given by According to (12), the epidermal contribution to S 585 (t) is removed by subtracting a suitably normalized signal S 600 (t). Based on the anatomical fact that no blood is present in the epidermis, we seek the lowest value β −1 that just eliminates the epidermal temperature rise in the image reconstructed from signal S 585 (t) − β −1 S 600 (t). The image reconstructed from such a signal is proportional to the temperature profile in the superficial part of the PWS, and thus reveals the location of its top boundary (see insets in figures 8(a) and (b)). (The described assessment of optimal β is an improvement upon our earlier reports on the DWE technique, and will be discussed in more detail elsewhere.) Analogous to β, the value of α is determined from prior anatomical knowledge that no melanin is present below the epidermal-dermal junction. Images reconstructed from signals y(t) (13) with increasing values of α differ little in the epidermal region, while the temperature rise deeper in skin decreases (Majaron et al 2000a). The lowest α that yields a zero temperature rise beneath the epidermal-dermal junction is selected, and is used to adjust the amplitudes of the PWS profiles according to (12).
Final DWE results obtained from the customary (3-5 µm) and spectrally narrowed PPTR signals are presented in figures 8(a) and (b) respectively. From the filtered signals, the PWS depth, as assessed at half maximum of the reconstructed temperature profile, is ∼140 µm ( figure 8(b), solid line). When utilizing the full camera bandwidth, in comparison, the PWS depth is determined to be ∼40 µm shallower, the temperature profile is broader, and has a 32% lower peak temperature (figure 8(a), solid line). Most importantly, the narrowing of the detection band resulted in a 30% decrease of the residual norm (from 16.8 to 11.5), despite the higher noise level in the filtered PPTR signals. This demonstrates an improved match between the measured and theoretically predicted signals with the spectrally narrowed signal acquisition.
The epidermal temperature profile reconstructed from unfiltered signals (figure 8(a), dashed line) decreases steeply from the skin surface inward. Since the overall epidermal absorption appears to be moderate (∼10% of the deposited energy), and the melanin concentration is known to be rather constant, or even increasing towards the epidermaldermal junction, such a temperature profile is difficult to reconcile with prior knowledge on light propagation in scattering tissues. A much more plausible result is obtained from the spectrally narrowed signals ( figure 8(b), dashed line), featuring a rather uniform epidermal temperature rise with a sub-surface peak at a depth of 50-60 µm, and a drop of ∼3 • C towards the skin surface.
Note that the images in figure 8  y(t) image (dashed line) begins to appear (Majaron et al 2000a). The latter is not related to melanin absorption, but results from a difference between the temperature profiles induced in PWS by the 585 and less-absorbed 600 nm irradiations. At greater depths, the proportionality between the two PWS temperature profiles, assumed in (11), deteriorates and an increasing share of the blood-related signal is present in y(t). Based on relation S 585 (t) = x(t) + y(t), and due to the linearity of (5), the actual PWS temperature profile (following 585 nm excitation) in this part essentially corresponds to the sum of both solutions.

Discussion
In reported PPTR image reconstruction algorithms, the computational complexity is reduced by using a fixed IR absorption coefficient value µ IR . Our more rigorous calculation of the attenuation function R(z) (7) indicates that such an approximation could be problematic in the presence of strong spectral variation µ IR (λ), characteristic of human skin in the 3-5 µm wavelength range (figures 1 and 2). Even determination of the optimal effective value µ eff is a non-trivial task, and any selection is somewhat arbitrary (figures 2 and 3). Calculating the average of µ IR (λ) weighted by detector sensitivity D(λ) (Prahl et al 1992), or by both D(λ) and B λ (T 0 ) (Milner et al 1995), for example, are but two approximations of the composite signal expression (7a). In our analysis, the attenuation function R(z) and resulting µ eff (z) also depend on D(λ). In the example discussed, however, this effect is dwarfed by the influence of µ IR (λ) (figure 3). Our numerical simulations (section 3) show that the spectral variation µ IR (λ) should be considered in the PPTR profiling of human skin using the 3-5 µm detection band. Using realistic, multi-spectral simulated signals, considerable distortions and artefacts appear in the reconstructed images of simple, single-lobed objects when a constant value µ IR is assumed in the reconstruction algorithm (figures 4(a) and 6(a)). In the example simulating a PWS layer (figure 4(a)), even a prolonged iteration yields a blurred (i.e. broader and lower) replica of the actual temperature profile, leading to an under-estimation of its depth by ∼30%. The simulated epidermal temperature profile is distorted to nearly a straight line when using such an approximate approach to PPTR profiling (figure 6(a)). (Note that a less rigorous approach to the numerical simulations was used in our earlier analysis, yielding less realistic and reliable results in Majaron et al (2000b).) We can attribute the above-mentioned image distortions to the mismatch between the spectrally composite signals (7a) and those predicted by the approximate expression (8) based on the following evidence: first, the features in the residual, much larger than the noise level, persist even after a prolonged iteration (figure 5, inset). In earlier PPTR experiments using the 3-5 µm detection band, we have regularly observed similar features in the residual. Secondly, the composite signal, used as an input to the reconstruction algorithm in the simulation, differs considerably from that predicted for the same object by the approximate expression (8) (figure 5; bold solid and dashed lines, respectively). And most importantly, faster convergence, lower residual norm and an improved match between the reconstructed images and actual objects resulted when the spectral variation µ IR (λ) was accounted for in the reconstruction algorithm (figures 4(b) and 6(b)).
In experimental situations involving higher noise levels than those used in our simulations, or other discrepancies between the measured and predicted signals, one may have to use a smaller number of iteration steps (or regularize the image reconstruction problem in another way). As a result, the benefit of the improved forward-problem formulation may be less pronounced as in the presented examples. Nevertheless, since the amplitude of the residual features in figure 5 amounts to 2% of the signal amplitude, we can expect the adverse effects of neglecting the µ IR (λ) under the conditions discussed to be present at SNR levels as low as 50 (see, e.g. figure 7). Moreover, since all sources of error combine in increasing the uncertainty of the initial temperature estimate, we believe that the improved forward-problem formulation as presented in this work could improve the stability in solving this inverse problem even at relatively high noise levels.
While the simulations demonstrate the deficiency of the constant µ IR approximation in PPTR profiling of human skin, using the more rigorous expression (7a) may not always be the optimal solution. In addition to increased computational complexity, this approach assumes detailed knowledge of the IR spectral characteristics of the tissue, detector and optics. Without further studies, it is difficult to predict the effects of possible inaccuracies in these specifications.
Instead, a suitable narrowing of the IR detection band may often be a more practical and also reliable approach. Specifically, the advantages of using the 4.5-5 µm detection band for PPTR profiling of human skin are two-fold. First, a constant value (µ eff = 26.5 mm −1 ) adequately describes the radiant emission under such conditions, which permits use of simpler reconstruction algorithms, based on the approximate forward problem formulation (4). This applies, in particular, to three-dimensional PPTR imaging, which is computationally much more intensive than depth profiling discussed in this work . Secondly, details of spectral characteristics of the optical system and detector thus have a negligible influence on PPTR signals. We found, for example, that even setting D(λ) = constant did not yield any observable change in µ eff (z) for the 4.5-5 µm detection band (figure 3, dashed line).
A downside of the latter approach is the associated reduction of the collected IR emission. In the discussed example, a ∼50% decrease was indicated by a theoretical estimate and verified experimentally. In a background-limited detection system, this would result in a ∼30% deterioration of the SNR. The singular value decomposition analysis, however, suggests that the related loss in PPTR profiling performance may not be critical (Milner et al 1995). Moreover, whenever the detector integration time has to be limited to prevent signal saturation, the reduction in the collected IR power due to spectral filtering can be compensated, at least in part, by increasing the integration time. In general, the advantages and disadvantages of the two suggested approaches to improved PPTR profiling performance should be weighed for a specific experimental set-up and application.
Our experimental test (section 4) indicates improved PPTR performance with a narrowed detection band. Admittedly, since the exact geometric and spectroscopic properties of the patient's PWS are unknown, direct assessment of image accuracy is not possible. However, if the proposed adverse effect of µ IR (λ) was negligible, equivalent images should be obtained with either detection band. Instead, several characteristics of the reconstructed temperature profiles have changed considerably with the implementation of spectral filtering. Most importantly, however, the residual norm was substantially reduced, despite the impaired SNR in the filtered signals, which convincingly demonstrates the improved match between the measured and theoretically predicted PPTR signals with the narrowed detection band.
Moreover, assuming that the temperature profiles reconstructed from the spectrally narrowed signals (figure 8(b)) are similar to the actual signals, the distortions observed in the full-bandwidth example (figure 8(a)) closely resemble those predicted by numerical simulations. The amplitude of the PWS profile is decreased by ∼30%, and the depth of PWS top boundary is determined as being ∼40 µm too shallow (compare with figure 4(a)). In addition, the epidermal temperature rise is distorted into a non-physical, almost linear profile, with an over-estimated surface temperature and overall epidermal thickness-exactly as seen in figure 6(a). The latter effect can also be observed in prior reports on PPTR profiling of human skin by Milner et al (1995Milner et al ( , 1996, and Majaron et al (2000a)-all utilizing the 3-5 µm detection range.
The dual-wavelength approach introduces an additional uncertainty in the values of parameters α and β, potentially inducing errors in the amplitudes of the epidermal and PWS components of the image, in accordance with (12) and (13). It has been reported that the area under the PPTR-determined temperature profiles correlates accurately with deposited laser energy (Sathyam andPrahl 1997, Smithies et al 1998). In our example, the sum of the energies indicated by the images reconstructed from x(t) and y(t) (figure 8) is within 4% of that indicated by the image from S 585 (t) alone ( figure 7).
Moreover, small variations of α and β have a minimal influence on the shapes of reconstructed profiles (Majaron et al 2000a). Determination of epidermal thickness and PWS depth, which are most important for the intended clinical application, are thereby only marginally affected. The presented example displays a good match between the sum of the images from x(t) and y(t) and the temperature profile reconstructed from S 585 (t) alone-within the error margin, with which the latter could be assessed (figures 7 and 8). This indicates that the possible inaccuracies in α and β are indeed non-critical. At the same time, the DWE approach exhibits improved reconstruction stability and accuracy of the PWS profile over the single-wavelength PPTR profiling, in particular when PWS blood vessels lie in close proximity to the epidermal-dermal junction ( figure 7).
The overlap of the epidermal and PWS temperature profiles indicated in figure 8(b), is consistent with prior anatomical knowledge of an uneven boundary between the epidermis and dermis. However, it may, in part, also result from the broadening of the two profiles due to heat diffusion. From the laser pulse duration of 1.5 ms, and the jitter (i.e. uncontrolled delay between the laser exposure and first acquired signal point, limited by the inverse sampling rate of 2 ms), the diffusion length amounts to approximately √ Dt ∼ 13-20 µm. In addition, the ultimate resolution of PPTR profiling (which is governed by heat diffusion between consecutive signal points, and thus, no less than √ Dt ∼ 15 µm in our example) is known to deteriorate quickly with increasing object depth (Milner et al 1995, Sathyam and Prahl 1997. The deeper part of the image reconstructed from y(t) (dashed lines in figure 8) does not represent melanin absorption, but results directly from different haemoglobin absorptions at 585 and 600 nm. As discussed in more detail elsewhere (Majaron et al 2000a), the proportionality between the two PWS temperature profiles, which is assumed in (11), thus deteriorates with depth and an increasing share of the blood-related signal is present in y(t).
Owing to linearity of (5), however, the actual PWS temperature profile in this part corresponds to the sum of both solutions (in accordance with relation S 585 (t) = x(t) + y(t)). But most importantly, inasmuch as the two PWS profiles can be linearized near its top boundary, equation (11) is valid in a shallow skin layer, including the superficial part of the PWS. The images in figure 8 are thus representative of the PWS and epidermal temperature profiles induced by 585 nm irradiation up to a depth of ∼0.3 mm, defined by the emergence of the deeper part of the y(t) image (dashed line). Based on the above considerations and experimental evidence (e.g. Majaron et al 2000a), the DWE-PPTR technique enables determination of PWS depth, epidermal thickness and epidermal temperature rise with sufficient accuracy to guide the laser therapy. A systematic study of DWE-PPTR profiling under realistic conditions is under way in our group (Choi et al 2002).
In conclusion, spectral variation of µ IR (λ) should be considered when applying PPTR to biological tissues. In particular, depth profiling of human skin using the 3-5 µm detection band benefits from accounting for the water absorption spectrum in the image reconstruction algorithm. Alternatively, narrowing the detection band to 4.5-5 µm allows the use of the customary approximation (i.e. fixed µ IR ), with improved performance. Under these provisions, PPTR has a potential for non-invasive characterization of PWS to optimize laser therapy on an individual patient basis. and SLO-US-2001/19). TM acknowledges support from the National Center for Research Resources at the NIH (RR14069). Institutional support from the Office of Naval Research, Department of Energy, NIH, and the Beckman Laser Institute and Medical Clinic Endowment is gratefully acknowledged.