Imaging laser heated subsurface chromophores in biological materials: determination of lateral physical dimensions.

. We describe a non-contact method using infrared radiometry to determine lateral physical dimensions of laser heated subsurface chromophores in biological materials. An imaging equation is derived that relates measured radiometric temperature change to the reduced two-dimensional temperature increase of laser heated chromophores. From measured images of radiometric temperature change, the lateral physical dimensions of chromophores positioned in an in vitro model of human skin are determined by deconvolution of the derived imaging equation using a non-negative constrained conjugate gradient algorithm. Conditions for optimum spatial resolution are found by analysis of a derived radiometric transfer function and correspond to superﬁcial chromophores and/or weak infrared absorption in a laser irradiated biological material. Analysis indicates that if the infrared attenuation coefﬁcient is sufﬁciently small (i.e., less than 10 mm − 1 ), infrared radiometry in combination with a deconvolution algorithm allows estimation of lateral physical dimensions of laser heated subsurface chromophores in human skin.


Introduction
Successful laser treatment of selected dermatoses requires controlled photothermal destruction of discrete subsurface chromophores in human skin. Because physical dimensions of targeted subsurface chromophores determine thermal relaxation time (τ r ), proper selection of laser pulse duration (t p ) is critical for successful treatment. For example, selective photothermolysis of port wine stain (PWS) blood vessels in human skin (Anderson and Parish 1983) requires matching laser pulse duration to thermal relaxation time of targeted vessels (i.e., t p ∼ τ r ) (Kimel et al 1993). In this case, τ r is proportional to squared diameter, and thus proper selection of t p requires knowledge of the lateral physical dimensions of the targeted blood vessels.
A number of investigators have reported application of pulsed photothermal radiometry using a single-element infrared detector to determine various optical, thermal, and structural properties of semiconductor (Tam 1987) and biological materials (Prahl et al 1992, Vitkin et al 1995, Milner et al 1995. We describe a non-contact method using infrared imaging radiometry to determine lateral physical dimensions of subsurface laser heated 0031-9155/96/010031+14$19.50 c 1996 IOP Publishing Ltd 31 chromophores in an in vitro biological material. In response to a diagnostic laser pulse, the true radiometric temperature increase is measured by a fast infrared focal plane array (IR FPA) camera. Due to inherent limitations of the camera (i.e., size of the collection aperture, lens aberrations, finite number and size of discrete detector elements) and thermal diffusion in the irradiated biological material, lateral physical dimensions of laser heated chromophores can not be directly measured in any one radiometric temperature image. In early images when blurring due to heat diffusion is least, the signal-to-noise ratio (SNR) of the true radiometric temperature increase is reduced when chromophores are imbedded deep within the material. Although SNR may be higher in later images as heat moves to the surface, lateral resolution decreases as thermal diffusion and blurring increase. To estimate lateral physical dimensions of laser heated subsurface chromophores in biological materials, we introduce the reduced two-dimensional temperature increase [ T 2D (x, y, t 0 )]. From a measured image of radiometric temperature change, T 2D (x, y, t 0 ) is determined by deconvolution of the derived imaging equation; lateral physical dimensions of subsurface chromophores are readily deduced from T 2D (x, y, t 0 ).
First, we derive the imaging equation and discuss use of a deconvolution algorithm to determine the reduced two-dimensional temperature increase ( T 2D (x, y, t 0 )) of laser heated subsurface chromophores from a measured image of the radiometric temperature change ( M(x, y, t 0 )). Second, experiments are described that measure the radiometric temperature change following pulsed laser irradiation of chromophores positioned at various depths in an in vitro biological material. Third, the reduced two-dimensional temperature increase of subsurface chromophores is computed by deconvolution of the derived imaging equation using selected measured images of radiometric temperature change ( M(x, y, t 0 )) and a non-negative constrained conjugate gradient algorithm. Following discussion of experimental results, we derive an analytic expression for a radiometric transfer function and examine effects of chromophore depth (z 0 ) and infrared absorption coefficient (µ ir ) on spatial resolution.

Theory
We derive an imaging equation that relates the reduced two-dimensional temperature increase ( T 2D (x, y, t 0 )) of laser heated subsurface chromophores in a biological material to a measured image of radiometric temperature change ( M(x, y, t 0 )). Our analysis includes effects of heat diffusion and inherent limitations of the IR FPA camera used to measure M(x, y, t 0 ). Because the imaging equation cannot be solved uniquely for T 2D (x, y, t 0 ), deconvolution of an image of measured radiometric temperature change ( M(x, y, t 0 )) is necessary. Given the biophysical properties of a laser irradiated material and specifications of the IR FPA camera, application of a numerical algorithm allows deconvolution of the derived imaging equation and thus estimation of lateral physical dimensions of subsurface chromophores.

The imaging equation
For the purpose of computing the thermal component of the imaging equation point spread function, we assume a δ pulse of energy, q[J ], is released within a semi-infinite medium (depth z = ζ 0 , time t = 0) representing a biological material (figure 1) and seek a Green's function solution, T δ , to the bioheat equation (van Gemert et al 1991), In (1) Q (s −1 ) represents the volume fraction of blood perfused tissue per unit time; κ (W m −1 K −1 ) and χ (m 2 s −1 ) are, respectively, thermal conductivity and diffusivity of the laser irradiated biological material which we assume are homogeneous. Additionally, we assume a Robin type boundary condition, in which the heat loss coefficient h (W m −2 K −1 ), represents combined radiative and/or convective thermal energy losses at the air-material interface. Solution of the bioheat equation ( T δ (ξ, η, ζ, t)) represents the increase in temperature above ambient values at position (ξ, η, ζ ) and time (t) in response to a δ pulse of energy released at position (0, 0, ζ 0 ) and time t = 0.
and erfcx(u) = exp(u 2 )erfc(u), where erfc(·) is the complementary error function [1−erf(·)]. Infrared emission from laser heated subsurface chromophores is scattered and absorbed prior to exiting the material. We assume reduced lateral resolution due to scattering of infrared radiation (Gaussian parameter σ s (m)), and an exponential loss (absorption coefficient µ ir (m −1 )), and write the relative increase in infrared emission at position r = (x 2 + y 2 ) 1/2 on the surface due to an infinitesimal volume of excess thermal energy (with unit temperature difference above the background level) located at position (ξ, η, ζ ), The Gaussian parameter, σ s , accounts for lateral blur due to scattering of infrared radiation emitted from the sample. Although the actual functional form of the blur is probably not Gaussian, the approximation in (4) is well founded when one considers scattering phase functions frequently used in light fluence computations in biological tissues (Ishimaru 1978). By integrating the excess thermal energy over the entire material volume and assuming that the blood perfusion relaxation time is much longer than the duration of measurement ((1/Q) t), we find the three-dimensional thermal point spread function (K T ), Substituting the thermal response of tissue ( T δ (ξ, η, ζ, t), (3)) to a δ pulse of energy, q (depth z = ζ 0 ), into equation (5), gives ρ t and C t are, respectively, the mass density (kg m −3 ) and specific heat capacity (J kg −1 K −1 ). Additionally, we have assumed that the infrared absorption coefficient (µ ir ) is independent of depth and constant over relevant temperature changes. Furthermore, the scattering parameter (σ s ) is assumed constant and characteristic of infrared emission from a given depth (∼ 1/µ ir ). Since q/(ρ t C t ) represents the initial temperature increase at position (0, 0, ζ 0 ), K T is written and consists of two physically distinct terms ((9a) and (9b) below) that represent, respectively, heat diffusion along lateral (K r ) and longitudinal (K z ) axes, Determination of K T allows computation of true radiometric temperature increase ( R(x, y, t)) as a convolution integral ((10) below) in terms of the initial three-dimensional temperature increase ( T 3D (x, y, z, t = 0)), The effect of lateral heat diffusion on the true radiometric temperature increase ( R(x, y, t)) can be understood by writing (10) in terms of the reduced two-dimensional subsurface temperature increase ( T 2D (x, y, t 0 ), (12) below) and convolution operator * The lateral distribution of the true radiometric temperature increase (i.e., the x, y dependence of R(x, y, t)) is given by a convolution integral of the thermal point spread function, K r , and T 2D (x, y, t 0 ). T 2D (x, y, t 0 ) (12) represents the reduced two-dimensional radiometric temperature increase of laser heated chromophores without effects of lateral blurring or limitations of the IR FPA camera. The objective of the work presented here is determination of T 2D (x, y, t 0 ) from a measured image of radiometric temperature change ( M(x, y, t 0 )).
Because the true radiometric temperature increase ( R(x, y, t)) is not measured directly, a realistic model should include inherent limitations of the IR FPA camera (i.e., the size of the collection aperture, lens aberrations, and the finite number and size of discrete detector elements). Since an IR FPA camera is a linear imaging system of incoherent light, the measured radiometric temperature change ( M(x, y, t 0 )) of a laser irradiated biological material can be written (Goodman 1968) as a convolution integral ((13) below) of the camera point spread function (K C ) and true radiometric temperature increase, R(x , y , t).
Since R(x , y , t) is expressible as a convolution integral (11), we may write an expression ((14) below) for the imaging equation in terms of the two-dimensional thermal (K r ) and camera (K C ) point spread functions, As written, the imaging equation represents a forward problem in which the recorded time sequence of images of the radiometric temperature change ( M(x, y, t 0 )) may be computed from the reduced two-dimensional temperature ( T 2D (x, y, t 0 )) and known biophysical properties of the material. Our objective, however, is to determine the reduced two-dimensional temperature increase ( T 2D (x, y, t 0 )) immediately following pulsed laser irradiation in order to remove blurring effects due to thermal diffusion or inherent limitations of the IR FPA camera. Inasmuch as only a timed sequence of images of radiometric temperature ( M(x, y, t)) is available as input data, determination of T 2D (x, y, t 0 ) in the imaging equation constitutes a deconvolution problem.

The deconvolution problem
The imaging equation (14) can be viewed as a two-dimensional deconvolution problem in which the unknown ( T 2D (x, y, t 0 ), (12)) is the reduced two-dimensional temperature increase of laser heated subsurface chromophores without blurring effects due to lateral thermal diffusion or the IR FPA camera. The imaging equation can be approximated as a linear matrix problem in which T 2D (x, y, t 0 ) and the point spread function (K = K c * K r ) are represented by, respectively, vector and matrix quantities ((15) below); here, bold symbols are used to represent vector and/or matrix terms.
The discretized imaging equation (15) is ill posed and cannot be solved directly because some singular values of the matrix K are essentially zero (i.e., less than the noise level).
To find a realistic solution to this equation, a regularization term ( T 2D ) is included and the L 2 or Euclidean norm, f , is minimized with respect to T 2D , Here M is a vector with N 2 components representing measured radiometric temperature change at time t 0 by individual detector elements in an N × N IR FPA; K is an N 2 × N 2 matrix representing the point spread function (K = K c * K r ) of the imaging equation; and T 2D is a N 2 -component vector containing the unknown reduced two-dimensional temperature increase of laser heated subsurface chromophores. The ill posedness of the deconvolution problem originates in attenuation of higherorder spatial frequencies due to lateral thermal diffusion ((9a); K r ) and inherent limitations of the IR FPA camera ((14), K c ). A variety of methods exist for solving linear leastsquares imaging problems; however, the large dimension of the deconvolution problem make iterative methods particularly attractive. Inasmuch as iterative spatial domain methods allow treatment of nonlinear problems and minimize 'ringing' effects (i.e., numerical oscillation), they are usually preferred over frequency domain techniques. Spatial domain methods, however, require storage and inversion of a N 2 × N 2 matrix and, therefore, cannot immediately be applied to solve the deconvolution problem. For example, when a 128×128 IR FPA is used, double-precision representation of the matrix K in the discretized imaging equation requires a storage capacity of greater than 2 GB. In contrast, frequency domain methods require less storage space and reduce the number of computations [≈ (N log N) 2 ] by utilizing the fast Fourier transform (FFT). A spatial domain algorithm that implements successive iterations in the frequency domain has been developed (Goodman et al 1992) that significantly reduces computation time. The method is based on conjugate gradients and uses a bending line search (McCormick 1989) and an active set strategy for implementation of the non-negativity constraint (i.e., T 2D (x, y, t 0 ) 0) to give a highly efficient algorithm for image processing problems.

Materials and methods
Experiments are described that record a time sequence of images of radiometric temperature change ( M(x, y, t)) in response to pulsed laser irradiation of an in vitro collagen material that contains subsurface chromophores positioned at various depths. An efficient nonnegative constrained conjugate gradient algorithm is used to deconvolve the imaging equation and thus estimate lateral physical dimensions of the chromophores.
In each experiment, a pulsed laser source was used to expose an in vitro collagen material containing subsurface chromophores. A three-element infrared lens (f/2, 50 mm focal length) imaged true radiometric temperature (R(x, y, t)) onto a 128×128 InSb IR FPA (AE-4128, Amber Engineering, Goleta, CA, USA). The in vitro collagen model was positioned approximately 190 mm from the infrared imaging lens, giving an image to object geometric magnification of 0.4. Efforts to increase spatial resolution by repositioning the lens resulted in unacceptable lens aberrations and large dark-field radiometric temperature variations due to vignetting. The IR FPA camera acquired 217 images of radiometric temperature per second and was externally triggered by a digital delay generator that was optically triggered by a fast silicon photoreceiver. The infrared signal collected by each detector element in the IR FPA was digitized with a 3.5 MHz 12 bit A/D converter, immediately stored in the computer's random access memory, and later downloaded to a magneto-optic disk storage device (figure 2). A bandpass infrared filter (3-5 µm) was positioned near the cold stop of the IR FPA to reduce background fluctuations and, hence, increase SNR, where M is defined as the mean measured radiometric temperature change due to laser heated subsurface chromophores and NE T is the noise equivalent temperature difference (Hopper 1993). Radiometric temperature changes measured in each image ( M(x, y, t 0 )) were calibrated using a resistor-heated aluminium surface coated with highly emissive (ε ≈ 0.967) black paint (TC-303 black, GIE-Tracor, Provo, UT) positioned in object space conjugate to the IR FPA. A surface mount thermistor (model No 8681, Keithley, Cleveland, OH) measured temperature (sensitivity < 0.01 • C) as the surface was slowly heated (1 • C min −1 ). An in vitro model using thin (i.e., 50-150 µm) collagen films (F1310, Collatec, Plainsboro, NJ) containing variable amounts of absorber was used to simulate discrete chromophores buried in multilayered composite human skin.
Chromophores were Figure 2. A schematic diagram of IRT instrumentation to measure the radiometric temperature change ( M(x, y, t 0 )).
constructed by staining a film with triphenylmethane dye (Aldrich Chemical Co., Milwaukee, WI) which absorbs optimally at the wavelength utilized in the experiments. Discrete chromophores were prepared by cutting a stained collagen film (w z = 125 µm) with known optical absorbance (µ a = 40 mm −1 ) into a number of thin strips (w x = 100-300 µm). The linear orientation of the stained collagen strips used in our study was selected to simulate blood vessels present in multilayered human skin. A model skin phantom was constructed by positioning variably spaced (50-700 µm) phantom blood vessels underneath a known thickness of weakly absorbing collagen films. Phantom blood vessels and collagen films were positioned on a 10 mm thick collagen sponge to simulate an infinite half-space as in in vivo skin. Pulsed (t p = 0.45 ms) radiation (λ = 585 nm) released from a flashlamp pumped dye laser (SPTL-1, Candela, Wayland, MA) was incident on the material at 30 • from the surface normal and covered a 30 mm 2 elliptical area. Immediately following exposure to a 1.0-1.5 J laser pulse, a time sequence of 100-200 images of radiometric temperature (M(x, y, t)) was recorded by the IR FPA camera. The relationship between lateral resolution and chromophore depth was investigated by placing the stained collagen strips underneath a variable number of weakly absorbing collagen films of different thicknesses (110, 180, 280, 380 µm). From each recorded time sequence of radiometric temperature (M(x, y, t)), the earliest time (t 0 ) was identified that gave images with greatest contrast. At each selected time, a given number of frames (217 (frames s −1 ) × (t 0 /10)s) centred about t 0 were averaged to give a mean radiometric temperature M(x, y, t 0 ). Representative times and corresponding number of frames for each model system tested (i.e., chromophores at depths of 110, 180, 280, 380 µm) were, respectively, 32, 69, 207, and 295 (ms) and one, one, four, and six frames. The mean radiometric temperature change at time t 0 ( M(x, y, t 0 )), was obtained by subtracting the non-irradiated background temperature and averaging over selected frames.

Results and discussion
The reduced two-dimensional temperature increase T 2D (x, y, t 0 ) was computed (figure 3) for each model by applying the non-negative constrained conjugate gradient algorithm using M(x, y, t 0 ) as input data. Assumed values of thermal diffusivity, conductivity, and heat loss coefficient in all calculations were, respectively, 1.1 × 10 −7 m 2 s −1 , 0.45 W mK −1 , and 15 W m −2 K −1 (Duck 1990). At deeper depths, lateral physical dimensions of subsurface chromophores determined by deconvolution increase because thermal diffusion strongly attenuates higher-order spatial frequencies required for fine spatial resolution.
Lateral physical dimensions are best resolved and closely match prepared values when chromophores are positioned closest to the surface ( figure 3(A)). Individual closely spaced (50-100 µm) chromophores positioned deeper than 280 µm are more difficult to resolve ( figure 3(C,D)). For example, at depths of 280 µm or deeper, resolution of three linear chromophores located in the upper leftmost region of the irradiation zone is not possible (figure 3(C)). A central pair of nearly parallel linear chromophores easily differentiated when irradiated at superficial positions is difficult to resolve when located at a depth of 380 µm ( figure 3(D)).

Discussion of results
Reduction in spatial resolution when imaging deeper chromophores is due to decreased radiometric temperature changes that originate from increased lateral (K r , (9a)) and longitudinal (K z , (9b)) thermal diffusion. Although SNR can increase in later frames as heat moves from chromophores to the material surface, greater lateral thermal diffusion blurs images and limits spatial resolution. Increased SNR may be obtained by use of higher laser irradiance which would produce a larger initial temperature increase. Because protein denaturation can significantly change thermal properties of biological materials (Si et al 1995), however, caution should be exercised when using higher laser irradiance. In deriving (14), for example, we assumed that thermal properties are constant and homogeneous during and following pulsed laser irradiation. When protein denaturation due to excessive laser heating occurs, thermal diffusivity (χ) changes and, consequently, errors and artifacts are introduced into the analysis.
Instrumentation and conditions employed in preliminary experiments reported here were far from ideal. For example, the IR FPA suffered from time-dependent temperature fluctuations. Although IR FPA non-uniformities were corrected prior to experimentation using a standard two-point hot-cold calibration procedure, the time-dependent nature of the fluctuations reduced SNR. In the following analysis, use of an optimized IR FPA camera is assumed and we analyse effects of chromophore depth (z 0 ) and infrared absorption (µ ir ) on lateral spatial resolution.

The radiometric transfer function
A radiometric transfer function ( R(k, t, z 0 )) is derived that relates peak-to-peak amplitude of true radiometric temperature increase in response to uniform heating of a linear array of identical subsurface chromophores in a biological material (figure 1). Specifically, we assume a homogeneous 1 • C initial three-dimensional temperature increase ((18) below) in a uniform array of infinitely long linear chromophores positioned at depth z 0 of rectangular cross section w x × w z with centre-to-centre spacing x 0 , where rect(·) and comb(·) are standard notation for, respectively, rectangle and comb functions (Goodman 1985). From (10), and explicit expressions for lateral (K r , (9a)) and longitudinal (K z , (9b)) point spread functions, the peak-to-peak amplitude of the true radiometric temperature increase at the fundamental spatial frequency (k 0 ) is computed ( R(k 0 = 1/ x 0 , t, z 0 ), (19) below) by integrating T 3D (x, y, z, t = 0) in an arbitrary (y, z) plane, followed by one-dimensional Fourier transformation along the x-axis (i.e., perpendicular to the long axis of each subsurface chromophore), and extraction of the k 0 harmonic.
The magnitude of µ ir is a spectral-and temperature-dependent function (i.e., µ ir = µ ir (λ, T )) of each molecular species present in the laser-irradiated material. Primary constituents in mammalian tissue are collagen and water; for example, in human dermis, water and collagen comprise, respectively, approximately 70% and 30% of in vivo mass (von Zglinicki et al 1993). The infrared absorption coefficient of collagen (Doyle and Blout 1975) lies in a relative local minimum between amide I and amide A-B absorption peaks approximately centred at, respectively, λ = 3.1 and 6.0 µm (Arrondo et al 1993). Similarly, the infrared absorption coefficient of water (µ ir (H 2 O)) also falls to minimum values between two absorption peaks approximately centred at λ = 2.9 and 6.6 µm (Hale andQuerry 1973, Wieliczka et al 1989). At physiological temperatures (∼ 33 • C), µ ir (H 2 O) varies over 10-40 mm −1 in the 4-5.4 µm spectral region and can decrease at least twofold for a moderate increase of 40 • C (Marechal 1991).
If we assume a uniform temperature increase T 3D (x, y, z, t = 0) ≈ 30 • C) in an infinite linear array of square chromophores (w x = w z = 50 µm, k 0 = 10 mm −1 ) positioned at a depth of z 0 = 250 µm with µ ir = 10 mm −1 , then the optimum contrast is R(k 0 , t 0 = 0, z 0 ) = 0.6 • C (figure 6). At a 120 Hz frame rate, use of a commercially available IR FPA camera that is detector noise limited allows NE T < 0.01 • C (Wilson et al 1992) or equivalent SNR ≈ 60. We conclude, therefore, that if µ ir is sufficiently small, use of an optimized IR FPA camera in combination with deconvolution of the derived imaging equation may provide a useful means to determine lateral physical dimensions of subsurface chromophores in human skin. Experiments that measure µ ir of human skin at various temperatures in the infrared spectral region are under way in our laboratory.

Conclusions
We have derived an imaging equation that relates the measured radiometric temperature change of a laser irradiated biological material to the reduced two-dimensional temperature increase of subsurface chromophores. Estimation of lateral physical dimensions of chromophores requires deconvolution of the derived imaging equation to determine the reduced two-dimensional temperature increase. Measurements of the radiometric temperature change in response to pulsed laser irradiation of an in vitro collagen model of multilayered composite human skin are applied to determine the reduced two-dimensional temperature increase and thus estimate lateral physical dimensions of prepared subsurface chromophores. Deconvolution of the derived imaging equation by an efficient non-negative constrained conjugate gradient algorithm indicates that spatial resolution rapidly degrades with increasing chromophore depth due to longitudinal and lateral heat diffusion. Derivation of a radiometric transfer function predicts that for increasing spatial frequencies a transition from longitudinal-to lateral-dominated diffusive blurring occurs when k 0 = µ ir /2π. Decay of optimum contrast for increasing spatial frequency is most rapid when blurring is dominated by longitudinal thermal diffusion. Sufficient contrast (SNR = 60) of smalldiameter (50 µm) subsurface chromophores (250 µm deep) such as port wine stain blood vessels in human skin heated by a diagnostic laser pulse (i.e., < 2.5 J cm −2 ) may be possible if µ ir is sufficiently small (µ ir ∼ 10 mm −1 ).