SPACE DEPENDENT TEMPERATURE INCREASE IN HUMAN SKIN SUBSURFACE CHROMOPHORES IMMEDIATELY FOLLOWING PULSED LASER EXPOSURE

Specifying the distribution of laser energy within a tissue is the first step toward understanding and capitalizing on a variety of laser-tissue interactions. Whether photothermal, photochemical, or photomechanical in nature, laser-tissue interactions begin with the absorption of photon energy. The spatial distribution of photon absorption specifies the required laser exposure to be delivered and the extent of subsequent therapeutic action. Using infrared tomography (IRT), the broad, long term objective of this research is the development of a three-dimensional tomographic reconstruction algorithm (TRA) as a means to determine the: (1) initial space-dependent temperature increase in subsurface chromophores (iTcHR(,1,t,t=O)) immediately following pulsed laser exposure; and (2) depths and physical dimensions of discrete subsurface chromophores. Analysis of the recorded time sequence of infrared emission images by longitudinal inversion and lateral deconvolution algorithms provides a direct means to determine the depths and physical dimensions of subsurface chromophores. Although our research is being shared with workers in a variety of disciplines, and pertinent to many clinical applications involving laser-induced photothermal mechanisms, we are particularly interested in addressing the problems associated with determination of the initial space-dependent temperature increase in subsurface chromophores in human skin in general, and port wine stain (PWS) blood vessels in particular. basis, an extended treatment Histopathological studies of PWS show a plexus of layers of dilated blood vessels located below the skin surface in the upper dermis. PWS blood vessel diameters vary on an individual patient basis, and evenfrom site to site on the same patient, over a range of 10-150 (15) The thermal relaxation time (tr) j5 defined as the time required for the core temperature, produced by the light within the target blood vessel, to cool to one-half of the original value immediately after the laser pulse and may be expressed

The IRT integral equation (Eq. 1) can be written as a multi-dimensional convolution integral that relates the measured time sequence of infrared emission images (LMCHR (x,y,t)) to the initial space-dependent temperature increase (TcHR (i,t= 0) ) in subsurface chromophores immediately following pulsed laser exposure, CHR(x,Y,t) J$KC(x-x'y-y')dx'dy'fflKT(x'-,y'-77,C,t)ATCHR(,77,c,t=O)dd17dc x',v' ,ii4 (1) '"CHR K * KT * ATCHR where KT(x-,y-r,t) and K(x-,y-i) represent, respectively, the thermal and camera point spread functions. Analysis of the recorded time sequence of infrared emission images (M(x,y,t)) by longitudinal inversion and lateral deconvolution algorithms provides a direct means to determine the depths and physical dimensions of subsurface chromophores.

Thermal and Camera Point Spread Functions
The important physical processes represented in the IRT integral equation, the thermal (KT) and camera (K) point spread functions, are schematically illustrated in where the expression for KT (Eq. 2) consists of two distinct terms (Eqs. 3a and 3b) that describe, respectively, heat diffusion along lateral (Kr) and longitudinal (K) axes.
K (x,y,t)= 1 e_22Y(2G4xt) r 2it(2 + 4t) (3a) Because the radiometric temperature at the tissue surface is not measured directly, we account for the inherent limitations of the JR-EPA camera and compute the camera point spread function (Ku) which includes effects due to size of the collection aperture, lens aberrations, finite number and size of discrete detector elements. Since the IR-FPA camera system uses a clear collection aperture and is an incoherent (i.e., thermal radiation) linear imaging system(4), the camera point spread function (Ks) is the Fourier transform (Eq. 4) of the modulation transfer function (MTFc). K = 3(MTFc); MTFC = MTFD .MTFA . MTF0 (4) An analytic expression for MTFc can be formed by a product of transfer functions representing the IR-FPA (MTFD), lens aberrations (MTFA), and camera optics (MTF0).

Noise Equivalent Temperature Difference (NELT)
Successful development of an IRT instrument system for determining the initial space-dependent temperature increase in subsurface chromophores (TcHR(,1,,t=O)) immediately following pulsed laser exposure requires a realistic determination of the noise equivalent temperature difference (NEST). The accuracy to which our tomographic reconstruction algorithm (TRA) can determine the initial space-dependent temperature increase and depths and physical dimensions of subsurface chromophores depends directly on NEzT.
In an optimized infrared detection system that measures radiometric temperature (i.e., 30-75°C). noise is limited by background emission fluctuations. The NEST (Eq. 5) is defined in terms of the numerical aperture (NA)2 of the infrared imaging optics with geometric magnification factor m; detectivity (D*(X)) of the JR-EPA at infrared wavelengths area (Ar) of each detector element; integration time (t) ofthe camera readout electronics; emissivity [c(A)] of the object under test; and spectral radiant exitance (M(X,T)) of the source as given by Planck's law5. The flashlamp-pumped pulsed dye laser (FLPPDL) has offered a superior approach in therapy due to its ability to destroy selectively port wine stain (PWS) blood vessels68. However, only a small proportion ofpatients obtain 100% fading of their PWS, even after undergoing multiple laser treatments912. The principal reason for poor results or treatment failure is that the pulse duration (0.45 ms) and light dosage ( 10 J/cm2) employed in current FLPPDL are, respectively, too short and too low to reach, and sustain over sufficient time, the critical core intravascular temperature necessary to destroy large PWS blood vessels irreversibly'3'14. Improved therapeutic outcome in many patients is expected from the development of laser systems with longer pulse durations (0.25-15 ms) and higher light dosages ( 25 J/cm2). Now that such laser systems have been developed, how will the clinician select the optiinalpulse duration and light dosage oflaser exposurefor PWS therapy, on an individualpatient basis, throughout an extended treatment protocol?
Histopathological studies of PWS show a plexus of layers of dilated blood vessels located below the skin surface in the upper dermis. PWS blood vessel diameters vary on an individual patient basis, and evenfrom site to site on the same patient, over a range of 10-150 (15) The thermal relaxation time (tr)j5 defined as the time required for the core temperature, produced by the absorbed light energy within the target blood vessel, to cool to one-half of the original value immediately after the laser pulse and may be expressed T -- where dBv is the diameter of the targeted blood vessel and is the thermal diffusivity of skin (1 .1 x 107m2/s). For vessels with diameters of 10, 70 and 150 tm, ;has calculated values of 0.06, 2.8, and 12.8 ms, respectively. The pulse duration of laser exposure (ta) governs the spatial confinement of heat and should, ideally, match the thermal relaxation time ('ç) for the targeted blood vessels'6. If longer pulse durations are employed (t,,>> 're), heat diffuses outside the vessel during laser exposure. The target specificity is reduced resulting in nonspecific thermal damage to adjacent structures. Alternatively, if too short a laser pulse is used (t,<< ;), high-peak intravascular temperature rises can produce explosive vaporization of tissue water, or photoacoustic transients which can result in vessel rupture. In such cases, repair mechanisms may revascularize the PWS. In conclusion, selection of the correct pulse duration of laser exposure is crucial to successful blood vessel destruction. Only when t -r can the critical core intravascular temperature, necessary to destroy large PWS blood vessels irreversibly, be achieved and sustained, for sufficient time.
It is important to be aware that like pulse duration, light dosage is also dependent on blood vessel diameter. For larger vessels, because energy deposition is limited by the optical penetration depth in whole blood (.ç' = 30 pm). Higher light dosages must be administered to reach, and sustain for sufficient time, the necessary critical core intravascular temperature'4"7. Therefore, for the clinician to select the optimal light dosage as a function of pulse duration correctly, knowledge of the initial space-dependent temperature increase in subsuiface PWS blood vessels (4TBv(c17,C,t=O) immediately following pulsed laser exposure is required.
For laser pulses in the millisecond domain, irreversible damage to PWS blood vessels occurs at core temperatures higher than 8O°C'4. Once zTBV has been determined from the IRT record in response to a diagnostic sub-therapeutic laser pulse of duration to and light dosage D0, the therapeutic laser light dosage (D) necessary to destroy the PWS blood vessels can then be computed from equation 7 knowing the critical temperature increase, ET = 50°C (assuming an ambient skin temperature of 30°C).
AT t (1_et0) For example, in response to a diagnostic laser pulse (t0 = 0.25 ms and D0 = 5 J/cm2), the initial space-dependent temperature increase immediately following pulsed laser exposure (zTBV) is determined from our three-dimensional tomographic reconstruction algorithm (TRA) to be 30°C in PWS blood vessels with a mean diameter of 70 im, inside the 5-7 mm circular beam spot of uniform light intensity. The therapeutic pulse duration (tp) iS calculated from equation 6 to be 2.8 ms. The therapeutic laser light dosage (D) required to heat the 70 im diameter PWS blood vessels to the critical core intravascular temperature (iT), necessary for destruction is then 12.6 J/cm2 (Eq. 7).
The rationale for using IRT in the clinical management ofpatients with PWS is that the technique offers a means of documenting the vascular characteristics of PWS on a site to site basis for each patient. Prior to the institution of laser therapy, in response to a sub-therapeutic diagnostic laser pulse, analysis of the measured time sequence of emission images (ziM,/x,y,t)) contained within the IRT record at each PWS site by the TRA provides information on the: (1) initial spacedependent temperature increase (Tij,,t=O)) within the blood vessels immediatelyfollowing pulsed laser exposure; and (2) depths and physical dimensions of the most supeificial absorbing PWS blood vessels inside the 5-7 mm circular beam spot of uniform light intensity. Such information can then be used by the clinician to help select the optimal parameters (pulse duration and light dosage)for theJkff laser treatment. Once the mean diameter ofPWS blood vessels (dBv) inside the 5-7 mm spot is determined from the deconvolved image, tp can be calculated using equation 6. With 4TB,/m,t=O) determined, the therapeutic laser light dosage (Dp5) to destroy the PWS blood vessels can then be computedfrom equation 7. At the next patient visit, a second IRT record is made. Because the majority of the most superficial PWS blood vessels were obliterated by the first laser treatment, the diagnostic laser pulse now images a deeper absorbing layer ofPWS blood vessels. The optimal parameters (pulse duration and light dosage) for the second laser exposure are selected based on analysis of the IRT record and a second treatment is peiformed. Proceeding in this manner, the IRT record probes deeper and deeper absorbing layers of the PWS after the most superficial layers are removed by each successive treatment at the same site.

IRT Integral Equation
The IRT integral equation (Eq. 1) forms the basis for our tomographic reconstruction algorithm (TRA) that estimates the initial space-dependent temperature increase in subsurface chromophores (TcHR(ii,,t=O)) immediately following pulsed laser exposure. As written, the IRT integral equation represents aforward problem in which the time sequence of measured infrared emission images (M(x,y,t)) may be computed from the initial space-dependent temperature increase (TcHR(1,,t= 0)) and the known biophysical properties of human skin. Our objective, however, is to determine the: (1) initial space-dependent temperature increase (i.e., TcHR(i,t=O)) immediately following pulsed laser exposure; and (2) depths and physical dimensions of subsurface chromophores. Inasmuch as only a time sequence of measured infrared emission images (M(x,y,t)) is available as input data, determination of TcHR(,i,t=O)) in the IRT integral equation constitutes an inverse problem (Eq. 8). The size of the IRT problem restricts the selection of candidate algorithms; a typical IRT record of 100 infrared emission images, collected with a 128 x 128 IR-FPA camera, represents a vector containing N = 1.6 million components. We seek a method that incorporates the physical constraints of the problem to find a solution in a reasonable time period utilizing affordable computing technology. We have developed computationally efficient one-dimensional 'ongitudinal inversion and two-dimensional lateral deconvolution algorithms. Both algorithms are based on conjugate gradients and use a nonnegativity constraint (i.e., TcHR(,i,,t=O) 0) for improved solution estimates.

Longitudinal Inversion Problem
If the laser beam diameter is larger than both the thermal diffusion length and size of the infrared object field and a large number of closely spaced subsurface chromophores occupy the region of laser exposure, then heat propagation in the lateral direction can be averaged over a small central area (GA) and the IRT integral equation approximates a one-dimensional longitudinal integral equation (Eq. 9) One-dimensional quantities M(t) and AT4(z, t -0) represent, respectively, the infrared emission image (M(x,y,t)) and the initial temperature increase (TcHR(,1,,t=O)) averaged over area öA positioned at the center of laser exposure (see computed simulation in Figure 3). The temperature change due to all subsurface chromophores at depth z (Eq. lOa), which is the solution to the longitudinal inverse problem [Eq. 9; A7(z,t = 0)], is then equal to the area-weighted average of the temperature increase in discrete subsurface chromophores immediately following pulsed laser exposure (Eq. lOb) where zT(z) and AThA represent, respectively, the temperature increase and fractional area in the (,i) plane at depth z of the ith chromophore.
The one-dimensional longitudinal inversion equation (Eq. 9) can be approximated as a linear matrix problem in which the initial one-dimensional temperature increase (AT(z,t = 0)) in subsurface chromophores and the longitudinal point spread function [Kz; Eq. 3b I are represented by, respectively, discrete vector and matrix quantities (Eq. 1 1); here, boldface symbols are used to represent vector and/or matrix terms.
'CHRKZATCHR (11) The difficulty in computing a solution estimate for the initial one-dimensional temperature increase (Ai;(z, t = 0)) in Eq. We have developed a very fast method that estimates a solution to Eq. 12 by applying our algorithm to the unaugmented problem (A=O) and achieve regularization by early termination (20) In this method, successive estimates of A7j(A=O) are computed using a non-negative constrained conjugate gradient algorithm and the optimum iterate23 is identified as the corner of the corresponding L-curve (i.e., plot of II T,(A)II vs. II M;:-K A7,(A)II for each iteration).

Lateral Deconvolution Problem
The IRT integral equation can be reduced to a two-dimensional inverse problem (Eq. 13a) in which the unknown [AT (ta); Eq. 13b] represents the lateral space dependent temperature increase at time to without blurring effects due to lateral thermal diffusion and the infrared optics.
where the lateral point spread function (KL) (Eq. 13c) accounts for blurring of subsurface chromophore images in human skin caused by lateral heat diffusion and limitations of the IR-FPA camera system.

KL K * KR (13c)
As in the case of the one-dimensional longitudinal inverse problem, a regularized matrix equation (Eq. 14) is written, Here "CHR is a vector with N2 components representing the infrared emission image recorded at time to by individual detector elements in an N x N IR-FPA; KL is an N2 x N2 matrix representing the lateral point spread function (KL, Eq. 13c); and AT' is a N2-component vector containing the unknown initial space dependent temperature increase in each discrete subsurface chromophore following pulsed laser exposure. The ill-posedness of the lateral inverse problem is due to attenuation of higher order spatial frequencies by inherent limitations of the IR-FPA camera system (Eq. 4; K) and lateral thermal diffusion (Eq. 3a; Kr). A variety of methods exist for solving linear least-squares imaging problems; however, the large dimension of the lateral inverse problem make our constrained conjugate gradient iterative method particularly attractive. We have developed an efficient spatial domain algorithm that implements successive iterations in the frequency domain, thus significantly reducing the number ofcomputations20). The algorithm is based on conjugate gradients and uses a non-negativity constraint so as to eliminate 'ringing" effects that often arise when using competing methods.

Layer-Stripping Method
Solutions to the longitudinal inverse problem [ AT(z,t=0)] can be applied to improve estimates of the lateral deconvolution problem when a uniform surface absorbing layer (e.g., epidermal melanin) overlies a distribution of deeper discrete chromophores (e.g., PWS blood vessels). In such cases, the temperature increase immediately following pulsed laser exposure in the surface absorbing layer [A7°(z,t=0)] is computed by solving the one-dimensional longitudinal inverse problem (Eq. 12) using our non-negative constrained conjugate gradient algorithm regularized by early termination. Knowledge of AT'°(z,t=0) allows computation (Eq. 15 of the infrared emission image [zM(x,y,t0)J exclusively due to the surface absorbing layer.
AM(x,y,t0) = $ K(t0, z)AT"(z,t = O)dz Infrared emission exclusively from the deeper discrete chromophores is then computed by subtracting 1Ms(X,Y,to) from a recorded infrared emission image [MHR(X,y,t)]. Solution of the lateral deconvolution problem [7 ;; Eq. 14] using the subtracted image as input data allows determination of the physical dimensions (6A) of discrete laser-heated subsurface chromophores. The mean temperature increase (4T0) of discrete subsurface chromophores that consist of a single absorber (e.g., hemoglobin in blood) is computed (Eq. 16) from the fractional area [f=(AThA)] and onedimensional temperature increase [ iT(z,t=O)] as determined from, respectively, the solutions of the longitudinal inversion and lateral deconvolution problems.

Model longitudinal inversion and lateral deconvolution algorithms
To test our longitudinal inversion and lateral deconvolution algorithms, we have formulated a general model system consisting of an arbitrary user-specified distribution of discrete subsurface chromophores in a composite multilayered structure. A model IRT problem has been formulated that allows numerical testing of our non-negative constrained conjugate gradient longitudinal inversion and lateral deconvolution algorithms. In our model, we assume that discrete subsurface chromophores of variable physical dimensions (2ai,2bi,2ci=50-350 j.tm) are positioned beneath the skin surface (zCHR=300±5O m) and heated (TcHR= 4_5°C) in response to pulsed laser exposure. To simulate light absorption by epidermal melanin, we assume a uniform temperature increase (55°C) in an absorbing surface layer (z=1O-75 tim). The temperature increase, depths, and physical dimensions of each chromophore are allowed to vary about mean values according to a uniform probability distribution. We have examined simulated infrared emission images (zM(x,y,t)) collected by the camera system using f/1.5 infrared imaging optics with 2X magnification interfaced to a 128 x 128 IR-FPA with 50 x 50 mm2 detector elements sensitive in the 3-5 j.tm spectral region.
From the time sequence of recorded infrared emission images (iMCHR(X,y,t)), we compute in Figure 3A the one-dimensional quantity AM;11 (t) by averaging over an area (i.e., A = 2 x 2 mm2) positioned at the center of laser exposure. Using our non-negative constrained conjugate gradient algorithm, we solve in Figure 3B the longitudinal inverse problem (Eq. 9) for the initial one-dimensional temperature increase A7,(z,t=O). Temperature increases in the surface absorbing layer and deeper chromophores are clearly observed and differentiated in the computed one-dimensional solution LZJJJ. Depths of the absorbing surface layer (z=5O jtm) and deeper chromophores (ZCHR=300±5O pm) determined by our longitudinal inversion algorithm ( Figure 3B) closely approximate the specified computer simulation depths. For the same distribution of chromophores used above in the one-dimensional longitudinal inverse problem, we compute (Eq. 1 3a) in Figure 4A j"cHR(X, y, t0 75 ms) . Infrared emission due to the surface absorbing layer [iM5(x,y,t)] is computed (Eq. 15) from the one-dimensional initial temperature increase ATD(z,t=O) as determined by our longitudinal inversion algorithm. Infrared emission exclusively from the deeper discrete chromophores is computed by subtracting M5(x,y,t=75 26 I SPIE Vol. 2623 40 N 1 ms) from a recorded infrared emission image [MHR(x,y,t=75 ms)] in Figure 4B. Solution of the lateral deconvolution problem (A;:; Eq. 14) using the subtracted image [McHR(x,y,t)-Ms(x,y,to)] as input data allows determination of the physical dimensions (6A1) of discrete laser-heated subsurface chromophores in Figure 4C. Physical dimensions of discrete subsuiface chromophores (dcHR=5O-35O /lm) determined by our lateral deconvolution algorithm match the specified computer simulation dimensions (d=50-350 jim). The mean temperature increase (T0 = 35°C) of laser-heated discrete subsuiface chromophores was computed (Eq. 16) from the fractional area [f=O.32] and one-dimensional temperature increase [AT,' (z,t=O)] as determined from, respectively, the solutions of the longitudinal inversion and lateral deconvolution problems.

EXPERIMENTAL RESULTS
Using in vitro and in vivo model chromophore systems, we have conducted a number of preliminary experiments that test the performance of our non-negative constrained conjugate gradient longitudinal inversion and lateral deconvolution algorithms. In each experiment, a 0.45 ms pulsed laser source (X=585 nm) was used to expose a test material containing subsurface absorbing structures that simulate PWS blood vessels. A compound infrared lens (f/5, 50 mm diameter) imaged the surface increase in infrared emission intensity from the test material onto a 128 x 128 InSb IR-FPA. The camera system acquired 217 infrared emission images per second and was externally triggered by a digital delay generator that was optically triggered by a fast silicon photo receiver. The infrared signal collected by each detector element was digitized with a 3.5 MHz 12-bit (0-4,095) A/D converter and immediately stored in the computer's random access memory ( Figure 5).  In vitro PWS Model A model phantom was constructed consisting of a number of closely spaced 175-225 jim diameter stained collagen strips (11a20 mm') positioned underneath two non-absorbing collagen films (total thickness; ZCHR=l5O±lO tim) and overlying a thick collagen sponge. Immediately following exposure to a 0.4 J laser pulse, a timed sequence of 100 infrared emission images (M(x,y,t)) was recorded over 460 ms. From this sequence, we determine in Figure 6A the one-dimensional quantity AM (t) by averaging over a circular area (A =3.0 mm2) positioned at the center of laser exposure. Using our nonnegative constrained conjugate gradient algorithm (33 iterations), we solve in Figure 6B the longitudinal inverse problem (Eq. 9) for the initial one-dimensional temperature increase, A7'(z,t=O), due to light absorption in the buried collagen strips. Depths of the collagen strips (zCHR=l5O±2O im) determined by our longitudinal inversion algorithm ( Figure 6B) match the prepared model phantom (zCHR=l5O±lO jtm). For the same phantom, infrared emission images (iiviCHR (x, y, t0 )) at 4.6 and 1 66 ms following pulsed laser exposure are shown in Figures 7A and 7B, respectively. Use of our non-negative constrained conjugate gradient algorithm (25 iterations) to solve the lateral deconvolution problem ( ; Eq. 14) allows determination of the physical dimensions (A) of discrete laser-heated buried collagen strips in Figure 7C by deconvolving the measured infrared emission image shown in Figure 7B to remove the effect of lateral blurring. Physical dimensions ofdiscrete collagen strips (dcHR=l75-225 jun) determined by our lateral deconvolution algorithm match the prepared model phantom dimensions (dcHR=l75-225 tm). The mean temperature increase (T0 = 27.5 oC) in discrete laser-heated collagen strips was computed (Eq. 16) from the onedimensional temperature increase I AT(Z,t=O)] and the fractional area [f=O.40] as determined from, respectively, the solutions ofthe longitudinal inversion and lateral deconvolution problems. For the same plexus of deep venules, we measured in Figure 9A iMBV(x, y, t0 =32 ms) . Infrared emission exclusively from the plexus of deep (zBv= 300 rim) lying venules is computed by subtracting M5(x,y,t0=32 ms) from a recorded infrared emission image [MvI5v(x,y,to=32 ms)] in Figure 9B. Solution of the lateral deconvolution problem (A7; Eq. 14) using the subtracted image [MBV(x,y,t)-MS(x,y,t)] as input data allows determination of the physical dimensions (A) of discrete laser-heated deep venules in Figure 9C. Physical dimensions ofdiscrete deep venules (dBv=l5O-200 jim) determined by our lateral deconvolution algorithm match the dimensions measured directly using light microscopy. The mean temperature increase (4T0 = 33°C) of laser-heated discrete venules was computed (Eq. 16) from thefractional area [f=O.32] and onedimensional temperature increase I tT,(,t=O)] as determinedfrom, respectively, the solutions ofthe longitudinal inversion and lateral deconvolution problems.

In vivo Chick Chorioallantoic Membrane (CAM) Model
An aluminum mirror positioned directly above and oriented at 45° to the irradiated CAM surface imaged the infrared emission intensity from the microvasculature into the IR-FPA camera system. Immediately following exposure to a 1.0 J laser pulse, a timed sequence of 140 infrared emission images (M(x,y,t)) was recorded over 660 ms. From this sequence, we determine in Figure 8A the one-dimensional quantity M1(t) by averaging over an area (A =1.5 mm2) positioned at the center of laser exposure. Using our non-negative constrained conjugate gradient algorithm (20 iterations), we solve in Figure 8B the longitudinal inverse problem (Eq. 9) for the initial one-dimensional temperature increases, AT(z,t=O), in the surface absorbing layer of capillaries (zBv 25 j.m) and deep lying venules (zBv = 300±75jtm). For the same PWS site, we measured in Figure 1 1A AMBV(x, y, t0 = 77 ms) . Infrared emission exclusively from the deep Pws blood vessels is computed by subtracting zM5(x,y,t=77 ms) from a recorded infrared emission image [MvI(x,y,t=77 ms)] in Figure 1 lB. Solution of the lateral deconvolution problem ( A7' ; Eq. 14) using the subtracted image [M(x,y,t)-M(x,y,t)J as input data allows determination of the physical dimensions @A) of discrete laser-heated deep PWS blood vessels in Figure 1 1C. Physical dimensions ofdiscrete PWS blood vessels were determined by our lateral deconvolution algorithm to be 150-200 um. The mean temperature increase (T0 = 4&C) oflaser-heated discrete PWS blood vessels was computed (Eq. 16) from the fractional area [f=O.42] and one-dimensional temperature increase [AT(z,t=O)] as determinedfrom, respectively, the solutions of the longitudinal inversion and lateral deconvolution problems.

PWS Patient IRT Measurements
We have verified the use of the high-speed IR-FPA camera system to image individual laser-heated PWS blood vessels. Immediately following exposure to a 1 .0 J laser pulse (5.0 J/cm2), the IR-FPA camera system recorded a timed sequence of 126 infrared emission images (M4,(x,y,t)). From this sequence, we determine in Figure 1OA the one-dimensional quantity &w(t) by averaging over an area (A = 1.5 mm2) positioned at the center of laser exposure. Using our non-negative constrained conjugate gradient algorithm, we solve in Figure lOB the longitudinal inverse problem (Eq. 9) for the initial onedimensional temperature increases, AT,(z,t=O), in the epidermal melanin layer and deeper PWS blood vessels. Depths of the epidermal melanin layer (z=5O±25 tim) and deeper PWS blood vessels (zBv=275±75 tim) determined by our longitudinal inversion algorithm (Figure lOB) match the mean depths measured directly using histopathology.