Coherent thermal wave imaging of subsurface chromophores in biological materials

Thermal wave imaging of discrete subsurface chromophores in biological materials is reported using a phase sensitive coherent detection technique applied to recorded infrared (IR) images. We demonstrate that utilization of a periodically modulated laser source for thermal wave excitation and coherent detection applied to each pixel may be used to compute images of thermal wave amplitude and phase at the laser modulation frequency. In comparison to recorded IR images, the narrow-band detection technique signiﬁcantly improves the quality of thermal wave amplitude images of subsurface chromophores in biological materials. Additionally, the technique provides phase information, which may be used to estimate chromophore depth in tissue. Application of the technique is demonstrated using tissue phantoms and in vivo biological models. We present a theoretical analysis and computer simulations that demonstrate the effect of tissue optical and thermal properties on thermal wave amplitude and phase. In comparison to the pulsed photothermal technique, coherent thermal wave imaging of subsurface chromophores in tissue for diagnostic applications allows reduction of peak incident laser ﬂuence by as much as four orders of magnitude and is safer and more amenable to in vivo imaging.


Introduction
Subsurface chromophores in biological materials affect light penetration and are important in understanding mechanisms of laser-tissue interactions.Identification of light-absorbing discrete structures and homogeneous substances in tissue is important for a variety of medical applications utilizing laser sources.Clinical practice often necessitates determining optical properties non-invasively to avoid non-specific tissue damage and reduce risk of post therapeutic complications.Non-invasive characterization of chromophores in tissue is a challenging problem due to the structural heterogeneity of biological tissues which often consist of multiple layers with different optical properties (e.g.skin) and may contain discrete 0031-9155/02/040657+15$30.00 © 2002 IOP Publishing Ltd Printed in the UK 657 chromophores (e.g., blood vessels, melanin, hair follicles) associated with particular tissue constituents.Various methods have been employed to visualize subsurface chromophores in biological materials (Chen et al 1997, Oraevsky et al 1997, Telenkov et al 1999).The choice of method depends on the problem under investigation and the specific response (photochemical, photomechanical or photothermal) of tissue to laser exposure.In the case of a photothermal response, relaxation of optical excitation results in non-homogeneous heat generation, which is transferred to adjacent structures via a diffusion process.The basic principle of photothermal radiometry (PTR) is that areas of dominant laser absorption in the substrate under study may be detected remotely using an infrared (IR) sensor.This technique has found numerous applications ranging from material testing to determination of optical properties (Leung andTam 1984, Prahl et al 1992).
The inverse heat conduction problem provides the mathematical framework for analysis of the pulsed photothermal response.Solving the inverse problem allows one to reconstruct one-and three-dimensional profiles of subsurface chromophores in a test material (Milner et al 1995, 1996, Telenkov et al 1999).Although this technique provides valuable tomographic information, noise in the recording bandwidth (1 Hz-10 kHz) limits imaging depth of subsurface chromophores in biological materials to 300-400 µm.Imaging of deep chromophores using the pulsed PTR method is difficult when an overlying layer contains a competing chromophore that contributes to the radiometric signal.In human skin, epidermal melanin acts as a competing chromophore to haemoglobin in dermal blood vessels and contributes a large fraction of the measured radiometric signal.To increase the imaging depth using pulsed PTR, various approaches have been tested.Although the simplest approach is to increase laser fluence, risk of chromophore denaturation and thermal damage to surrounding tissue structures limit the maximum fluence.Recently Majaron et al have demonstrated a two-wavelength excitation technique that improves imaging depth and may be applied to distinguish between epidermal melanin heating and dermal blood vessels.
Here, we approach photothermal imaging of subsurface chromophores in biological materials by generation of thermal waves using periodically modulated laser excitation and coherent detection of radiometric temperature increase.Although fundamentals of thermal wave phenomena in condensed matter have been known for decades (Carslaw and Jaeger 1959) and were extensively utilized for nondestructive material testing, we believe application of thermal waves to imaging subsurface chromophores in biological materials is novel.In earlier studies (Busse et al 1992, Busse andWalther 1992), intensity modulated laser radiation was utilized to stimulate thermal waves in opaque industrial materials.Measurement of thermal wave amplitude and phase provided data to visualize discontinuities of thermal properties due to subsurface defects and material inhomogeneities.Murphy et al (1992) analysed excitation of thermal waves by a laser-induced spatially distributed heat source in homogeneous and layered materials.Although an expression for the temperature field at the material surface was determined, photothermal detection and the effect of volumetric emission of IR radiation were not considered.In general, thermal wave studies of inorganic samples are limited to detection of surface temperature increase due to exponentially attenuated laser excitation and thermal wave scattering due to subsurface inhomogeneities.Inasmuch as the objective of our study is to image discrete subsurface chromophores in biological materials, our analysis includes two effects not contemplated in earlier studies.First, laser excitation of the thermal wave is confined to a discrete spatial region corresponding to the position of the subsurface chromophore.Second, our analysis includes effect of volumetric emission of IR radiation from the test material.
In our approach, a time sequence of two-dimensional IR images of a test specimen is recorded during excitation by periodically modulated laser radiation at a wavelength near maximum of chromophore absorption.The time sequence of two-dimensional IR images recorded in response to the excitation by periodically modulated laser radiation is referred to as a 'thermal wave radiometric image'.A coherent signal detection algorithm (Poor 1994) is applied to compute amplitude and phase of the thermal wave radiometric image at the laser modulation frequency.When imaging subsurface chromophores within in vivo specimens, thermal wave imaging provides a number of important advantages as compared to pulsed PTR.First, when the oscillating radiometric signal is processed using a narrow-band coherent detection technique, the signal-to-noise ratio (SNR) can be increased significantly due to suppression of the noise bandwidth.Improved SNR allows decreasing incident laser fluence, an important consideration for biomedical applications involving in-vivo specimens.Second, coherent detection of the oscillating radiometric signal at a single frequency provides phase information, which can be used to estimate chromophore depth.Although chromophore thickness may be estimated by measuring phase at various excitation frequencies using a chirped laser intensity profile, determination of chromophore depth is sufficient in some applications.Third, the coherent thermal wave imaging technique may be used for tissue monitoring over an extended period of time to detect relatively slow variations in optical and thermal properties due to physiological effects (Liu and Xu 1999) or changes in biochemical concentration.
Haemoglobin is a primary chromophore of tissue in the visible spectral range 400-700 nm.We utilize visible laser radiation and coherent thermal wave imaging to visualize subsurface blood vessels in tissue.Specifically, we are interested in imaging blood vessels in the upper dermis of skin, because quantitative information on the vascular structure in skin has important implications for laser treatment of hypervascular lesions (van Gemert et al 1995) such as port wine stain birthmarks.Our experimental studies include imaging of discrete chromophores in epoxy phantoms, blood vessels in the chick chorioallantoic membrane (CAM) model and in vivo dermal blood vessels in a rodent animal model.We present a theoretical analysis of the three-dimensional thermal wave problem and numerical simulations of amplitude and phase of the thermal wave radiometric image.

Experimental apparatus
Thermal wave generation in biological materials can be achieved using amplitude modulated radiation emitted by a CW laser with a wavelength near the maximum of chromophore absorption.Our experimental apparatus for coherent thermal wave imaging of subsurface blood vessels is shown in figure 1(A).When imaging blood vessels in skin, second harmonic of a Nd:YAG laser (λ = 532 nm) can be used.This excitation wavelength provides absorption contrast of one order of magnitude between blood vessels and the epidermis (absorption coefficients for blood and epidermis α bv = 26.6 mm −1 and α e = 2.1 mm −1 , respectively (van Gemert et al 1995)).Continuous radiation emitted from the Nd:YAG laser is intensity modulated by an electronic shutter, creating a symmetric square waveform for synchronous signal processing (figure 1(B)).Heat generation in subsurface blood vessels and the subsequent heat diffusion process produce a spatial temperature distribution in tissue oscillating with the laser modulation frequency f m = 1/T m (figure 1(C)).The oscillating temperature field can be conveniently described in terms of the Helmholtz wave equation with a complex wave number σ = (iω m /D) 1/2 , where ω m = 2πf m is the angular modulation frequency and D the thermal diffusivity of tissue.Mathematical analysis of thermal wave phenomenon can be found in numerous studies (Carslaw and Jaeger 1959, Busse and Walther 1992, Murphy et al 1992).A harmonic thermal wave in space is characterized by an oscillating temperature with exponentially decreasing amplitude and phase lag proportional to the distance between the or diffusion length (L D ) depends on thermal diffusivity (D) of the media and laser modulation frequency (ω m = 2πf m ).Inasmuch as thermal diffusivity of tissue is relatively low (D = 0.11 mm 2 s −1 ), laser modulation frequencies must also be low to achieve a sufficiently large thermal diffusion length for imaging deep chromophores.For example, a typical depth of blood vessels in dermis of 0.1-0.5 mm gives a range of modulation frequencies f m = 0.14−3.5 Hz.
The phase lag (measured in radians) between temperature oscillations separated by distance |r -r 0 | can be written as Increase in emission of blackbody radiation due to periodic laser excitation is recorded by an IR focal plane array (FPA) camera.In our experiments we use an InSb IR FPA camera (Phoenix model, Indigo Systems, Santa Barbara, CA) sensitive in the 3-5 µm spectral band.The array of 256 × 256 pixels and IR optics allows imaging of a 9 × 9 mm 2 specimen area giving a spatial resolution of 35 µm/pixel.The IR camera frame rate and total number of recorded frames (N ) determine the image acquisition time (t N ).The maximum number of recorded frames is limited only by the quantity of available computer random access memory (RAM).Longer acquisition times, however, create practical difficulties in processing large three-dimensional data arrays.In our experiments, the total recording time was 20.5 s corresponding to 1024 images recorded at a 50 Hz frame rate.We note that in comparison to traditional problems in nondestructive material testing, where samples are often opaque and optical radiation is absorbed in a thin superficial layer, a biological specimen with selectively excited chromophores can be considered as thermally homogeneous media with distributed discrete heat sources.This difference has implications for the detected photothermal response.For traditional problems in nondestructive material testing, contrast in IR images is attributed to thermal waves scattered from a discontinuity of the thermal diffusivity (Nicolaides and Mandelis 1997).For laser excitation of tissue, the temperature increase is due to heat diffusion from discrete chromophores into the surrounding tissue.Because of the evanescent nature of thermal waves, we expect a greater imaging depth compared to traditional nondestructive testing for similar detection conditions.

Amplitude
Another important factor in thermal wave imaging of biological materials is absorption coefficient at the IR detection wavelength (α ir ).Inasmuch as the fractional water mass content in tissue can be as high as 70%, attenuation of IR emission is similar to that in water (Hale and Querry 1973).The finite value of α ir (λ) complicates quantitative analysis because an IR sensor records the radiometric signal that originates within a characteristic depth α −1 ir dependent on the emission wavelength.To simplify our analysis, we characterize optical properties of test specimens in the 3-5 µm spectral range (sensitivity of InSb detectors) by a constant absorption coefficient α ir .In a more detailed analysis, total radiometric signal should be modelled by an integral over the IR detection band.
Following image acquisition, the recorded time sequence of two-dimensional frames is processed using a coherent detection algorithm to compute amplitude (A ij ) and phase (θ ij ) at each pixel (i, j) of the thermal wave radiometric image.The block diagram of the coherent detection algorithm applied to the signal s ij (t) from pixel (i, j) is shown in figure 2 and is similar to that utilized in digital lock-in amplifiers (Meade 1983).We implement this algorithm as a software module to compute amplitude and phase at the laser modulation frequency f m for each pixel throughout the sequence of N radiometric image frames.This technique is equivalent to 256 × 256 synchronized lock-in amplifiers detecting radiometric signals from adjacent pixels.A difference between previously reported lock-in methods (Ahmed et al 1989) is that our implementation does not require additional hardware for coherent processing of high-speed video data flow.However, significant computer RAM (up to 500 MB) is required for storage of digital image data recorded by the IR FPA camera.
A sinusoidal reference signal r(t) is generated within the software module in-phase with the square modulation waveform.An additional constant phase shift θ r can be introduced into the reference signal from the software interface to avoid phase wrapping.Following the multiplication operation, the in-phase and quadrature components pass through a lowpass filter with the specified time constant τ 0 which determines the signal detection bandwidth centred at the laser modulation frequency f m .The outputs of the coherent processing algorithm are amplitude [A ij (f m )] and phase [θ ij (f m )] of the thermal wave radiometric image at the laser modulation frequency f m .In our experiments we set τ 0 = t N to maximize SNR and compute an image of both amplitude and phase from the recorded time sequence.An additional feature of the coherent imaging technique arises from the fact that spectrum of the square laser modulation waveform and time-dependent radiometric signal contains odd harmonics of the fundamental frequency (i.e.f m , 3f m , 5f m , etc).This feature of the modulation waveform and time-dependent radiometric signal can be utilized for imaging specimens at higher harmonics using the recorded sequence of IR frames by simply changing frequency of the reference signal.According to equation (1), radiometric signal due to excitation by higher harmonics corresponds to a shorter heat diffusion length and may be used to discriminate contribution of chromophores positioned at different depths without acquiring a new set of image frames.In practice, however, radiometric signal amplitude at frequencies beyond the fifth harmonic (5f m ) is small and difficult to detect.
We demonstrate the application of coherent thermal wave imaging using several test specimens.Initially, light scattering epoxy phantoms with embedded cylindrically shaped discrete chromophores (diameter 50-100 µm) were studied to verify system performance.We used two biological in vivo models for imaging subsurface blood vessels: (1) a CAM and (2) the rodent dorsal skin flap window.Finally, to understand the observed spatial variations in recorded phase images [θ ij (f m )] we solve numerically the three-dimensional thermal wave problem and compute radiometric phase response.

Epoxy phantoms with embedded chromophores
A homogeneous epoxy sample with discrete cylindrically shaped chromophores positioned just below the surface is a simple and convenient model of tissue with discrete chromophores.We utilized an epoxy resin with aluminium oxide powder (Master Bond Inc, Hackensack, NJ) to simulate optical and thermal properties of tissue.Similar phantoms have been used by Torres et al (1999) to investigate the thermal response of tissue to cryogen spray cooling.The cured epoxy was highly scattering to visible light with an estimated thermal diffusivity of D = 0.25 mm 2 s −1 (Telenkov et al 2001).Two chromophores were oriented obliquely to the phantom surface so depths (measured from phantom surface to chromophore centre) ranged between 250 and 500 µm.Results of thermal wave imaging at laser modulation frequency f m = 0.5 Hz are shown in figure 3. The image depicted in figure 3(A) demonstrates one of 1024 frames recorded by the IR camera when the shutter was in the open position.Traces of two laser-heated chromophores inside the epoxy phantom are not clearly visible in the recorded image.After coherent processing of the image sequence at f m = 0.5 Hz, amplitude and phase of the thermal wave radiometric image were computed (figures 3(B) and 3(C)).Increase in SNR (estimated maximum improvement factor is 512) due to narrow-band detection allows resolution of individual chromophores in the amplitude image.The phase image shown in figure 3(C) reveals subsurface chromophores as dark linear features, because the grey scale assigns darker colour to smaller phase values.Temperature oscillation of a black paint marker on the surface (M) appears in-phase with laser heating while subsurface features display a nonzero phase lag.Additionally, the chromophore labelled '2' is oriented oblique to the phantom surface, which is indicated by the gradual decrease of signal amplitude and increase in phase towards the top of the image.Estimates of chromophore depth using equation (2) give values, which differ from the true depth (measured from phantom surface to chromophore centre) by as much as 30%.A hypothesis for the observed discrepancy was tested using numerical simulations of the thermal wave response described in section 4.

Chick chorioallantoic membrane (CAM) microvasculature
The chorioallantoic membrane of a fertilized chicken egg is well suited for photothermal imaging of discrete blood vessels in vivo (Kimel et al 1994).A transparent liquid media allows observation of the microvascular network directly through a light microscope.To gain access to the membrane blood vessels, a 25 mm window was cut in the egg shell on the fifth day of embryonic development.The chicken egg was positioned horizontally and FPA camera recorded IR emission images reflected from a dichroic mirror situated just above the egg.Image data recorded during periodic laser excitation at 1.0 Hz is shown in figure 4(A).As with the epoxy phantom, the recorded IR image appears strongly blurred and no discrete blood vessels can be resolved.The phase image depicted in figure 4(C) reveals a microvascular structure which appears as a network of contrasting dark/light lines and contains more detail than that observed in the amplitude image.An interesting feature of the phase response is observed at the blood vessel boundaries which appear as white on the grey-scale image.We hypothesize that these areas appear white due to simultaneous contributions to the radiometric signal from heat generated in blood vessels and background tissue structures.Heat generated in blood vessels diffuses outwards producing an increasing phase at larger lateral distances whereas heating of background tissue structures dominates the signal at large lateral distances from the blood vessel.Our theoretical model for the coherent thermal wave response (section 4) predicts the observed phase behaviour for superficial blood vessels positioned within a thermal diffusion length L D .

Rodent dorsal skin flap window model
To model photothermal response of human skin we utilized a rodent dorsal skin flap (Papenfuss et al 1979) window exposed to periodic laser excitation.This in vivo model allows the simultaneous observation of the epidermal and subdermal sides of the skin while maintaining hydration and function.This model has been used extensively as an in vivo model for optical imaging of blood vessels on the subdermal side of skin while still having access to the epidermal side (Barton et al 1999).
Preparation of the window consists of shaving and epilating the entire dorsal area of a hamster.The skin is pulled away from the body along the dorsal midline and sutured to a vertical c-clamp.A circle of 10 mm in diameter is cut from one single thickness of skin, exposing the subdermal side of the opposed tissue.An aluminium chamber is sutured to both sides of the skin and a glass window is placed over the cutout section to prevent the tissue from dehydrating on the subdermal side.During the experiments, the glass window is removed for improved imaging, however hydration is maintained with application of physiological saline.A subdermal view of the prepared hamster dorsal window is recorded using a CCD camera and is shown in figure 5(A).Blood vessels in the skin of the dorsal window are located in the subdermal fat and connective tissue layers 400-500 µm beneath the epidermal surface and 80-100 µm from the exposed subdermal surface.Arterioles and venules in the dorsal skin range in inner diameter between 50 and 500 µm, with the main arteriole and venule approximately 350-500 µm in diameter (marked 1 and 2 in figure 5(A)).The periodic laser radiation is incident on the epidermal side of the skin providing uniform illumination over a 12 mm diameter spot at a fluence of 0.38 W cm −2 .IR image acquisition is configured for imaging blood vessels through the epidermal side (figure 1).Detected amplitude and phase of thermal wave radiometric image at the laser modulation frequency f m = 1.0 Hz are shown in figures 5(B) and (C) (left-right is inverted).Although an apparent correlation in vascular structure between the optical image and the radiometric response (blood vessels 3 and 4) due to thermal wave generation is observed, small capillaries and relatively deep (greater than 200 µm) blood vessels remain unresolved in the amplitude and phase images.Moreover, the phase image shows an inverted response, i.e. deep blood vessels (1 and 2) appear white on the grey scale whereas superficial capillaries appear dark.Similar to the CAM phase image (figure 4(C)) blood vessel boundaries appear highlighted by an increased phase over the background response.

Numerical results and discussion
In order to gain a better understanding of the coherent thermal wave response, we evaluate amplitude and phase of thermal wave radiometric images by solving the three-dimensional heat conduction equation for a periodic volumetric source, and compare theoretical results with the experimental data.Since our method detects amplitude and phase of the thermal wave radiometric image at the laser excitation frequency, analysis of the imaging problem in the frequency domain is convenient.The temporal Fourier transform is used to establish a relationship between the time sequence of recorded IR frames and the spectral content of the radiometric response.The objective of our theoretical analysis is to compute an analytic expression for the IR radiometric image R(x, y, f m ) at the laser modulation frequency, f m , for a given configuration of the volumetric source.Since an analytical expression for R(x, y, f m ) is given in complex form, amplitude and phase of the thermal wave radiometric image can be computed after extracting real and imaginary parts: The IR absorption coefficient (α ir ) over the detection wavelength range relates radiometric images and the space-dependent thermal wave temperature field in tissue T (x, y, z, f m ) at the excitation frequency via an integral: R(x, y, f m ) ∝ α ir ∞ 0 T (x, y, z, f m ) e −α ir z dz.
(4) Equation ( 4) assumes that the temperature increase is small enough to use a linear approximation for the radiometric response, and takes into account exponentially attenuated radiometric signal from all depths in tissue.The space-dependent temperature field T (x, y, z, f m ) can be determined from the heat conduction equation, which in the frequency domain can be written in a form similar to the Helmholtz wave equation: where σ = (i2πf m /D) 1/2 is a complex wave vector, κ thermal conductivity and Q(x, y, z, f m ) spectral component of the laser-induced volumetric heat source at the laser modulation frequency f m .Additionally, we assume a homogeneous Neumann boundary condition at the surface (z = 0), i.e.
Equation ( 5) with the boundary condition (6) can be solved using a Green's function in the frequency domain (Mandelis 1995).Since the Green's function represents the temperature response due to a point heat source positioned at r 0 = {x 0 , y 0 , z 0 }, the temperature field T (x, y, z, f m ) due to an arbitrary source is given by a convolution integral of the Green's function and the spatial distribution of heat sources.However, this straightforward approach cannot be applied to compute the radiometric image from equation (4) because of the singularity in the three-dimensional Green's function at the source position r = r 0 .To overcome the singularity problem,we analyse transversal heat diffusion in the spatial frequency domain and replace σ by an effective wave vector where k 1 and k 2 are components of vector k in the two-dimensional spatial frequency space (Murphy et al 1992).The effective wave vector σ e describes combined heat diffusion in a lateral plane and in the normal direction.By calculating the two-dimensional Fourier transform of the temperature field, equation ( 5) can be rewritten in a form similar to the one-dimensional equation in z-space: where are the two-dimensional Fourier transforms of the temperature field and the source function in a lateral plane: Solution to the differential equation ( 7) for transverse spectrum of temperature T (k 1 , k 2 , z, f m ) is given in terms of a regular one-dimensional Green's function with the effective wave vector z+z ) .(9) In this approach, the Green's function (9) represents the transverse spectrum of temperature at the coordinate depth z due to an oscillating point heat source positioned at z .Substitution of equation ( 8) for the transverse spectrum of temperature as a function of z, into equation (4) allows us to compute the spatial frequency spectrum of the radiometric temperature image: where is the radiometric response in the spatial frequency domain due to an infinitesimally thin planar heat source at the coordinate z .Integral (10) can be evaluated using result (11) and the spatial distribution of the heat source.In the case of an arbitrary source transverse spectrum Q(k 1 , k 2 , z, f m ), equation ( 10) can be computed numerically.Finally, the complex radiometric image R(x, y, f m ) is determined using the two-dimensional inverse Fourier transform applied to the spatial spectrum (10), and amplitude and phase of the thermal wave radiometric image are computed from equation ( 3).This mathematical framework can be used to compute the spatial frequency spectrum of the radiometric image, as well as images of amplitude and phase for a given three-dimensional heat source distribution [ Q(x, y, z, f m )] avoiding the difficulties associated with the singularity of the three-dimensional Green's function.
We compute radiometric images for the geometry shown in figure 6(A).The simulated specimen occupies a half-space and contains a single chromophore with dimensions 0.2 × 3.0 × 0.1 mm 3 positioned at the depth d = 300 µm (from the specimen surface to the chromophore centre).We assume that heat is generated due to absorption of harmonically modulated laser radiation both in the chromophore (absorption coefficient α ch = 10 mm −1 ) and surrounding tissue structures (absorption coefficient α 0 = 0.5 mm −1 ).The rate of heat generation is given by Q(x, y, z, t) = α 0 e −α 0 z + α ch e −α 0 z ch −α ch z g(x, y, z) I 0 e i2πf m t (12 where I 0 is the intensity of incident laser radiation, z ch depth of the chromophore top surface and g(x, y, z) the chromophore indicator function, i.e. g = 1 for {x, y, z} within chromophore and g = 0, for outside.We neglect non-uniform laser fluence in a lateral plane because in the experiments we use a flat-top laser beam with the diameter exceeding the camera field-of-view.Plot of the heat generation rate along the z-axis passing through the subsurface chromophore is shown in figure 6(B).Output of the numerical routine implementing equations ( 3)-( 11) for laser modulation frequency f m = 1.0 Hz and IR absorption coefficient α ir = 10 mm −1 is presented in figure 7. Amplitude image (figure 7(A)) shows the subsurface chromophore, which appears blurred due to heat diffusion in the transverse plane.Phase image (figure 7(B)) indicates position of the chromophore as a dark band, which corresponds to smaller values of phase shift.A characteristic feature of the computed phase image is the transition area from within the chromophore to background areas, which appear as a white halo on the grey-scale image.The halo prediction of the theoretical model is consistent with the observed phase response in CAM and rodent skin flap model (figures 4(C) and 5(C)).The phase profile along a horizontal line passing through the image centre (figure 8(A)) indicates a phase minimum at the chromophore and two maxima at both edges.However, for the ideal case of no background absorption  (α 0 = 0), our model predicts a linear increase of phase in both directions (figure 8(B)), which follows from the theoretical description of thermal waves (2).These numerical results suggest that the nonspecific deposition of laser energy in tissue structures surrounding the target choromophore, produce a radiometric signal due to background heating that can eventually dominate the thermal wave response.Competition between the two signals results in the unusual phase response observed in our experiments and predicted by the numerical simulations.
We further investigate the thermal wave phase response using our numerical model by varying the chromophore depth.The thermal wave phase image was computed for three chromophore depths d (175, 275 and 375 µm); phase profiles are shown in figure 9(A).These data demonstrate that light absorption in surrounding media (α 0 ) and thermal diffusion length  chromophore centre and two maxima near the edges (figure 9(A)).In the phase response of deeper chromophores (d > L D ), no central minimum is observed and phase is greater than that of the background response (figure 9(A), curve 3).Our numerical model suggests this effect should invert the appearance of chromophores in grey-scale phase images with respect to a uniform background and was indeed observed in phase images of blood vessels in the rodent skin model (figure 5(C)).
Since we are interested in determining chromophore depth from the phase measurements, dependence of phase as a function of chromophore position along the z-axis is important.Computed phase at a position centred directly above the chromophore versus depth for three laser modulation frequencies (0.5, 1.0 and 2.0 Hz) are shown in figure 9(B).The numerical simulations indicate that radiometric phase can be approximated by a linear function of depth with a slope proportional to the square root of the laser modulation frequency f 1/2 m only for superficial chromophores.The thermal wave phase becomes increasingly non-linear for deeper chromophores and saturates at z ≈ L D .Extrapolation of the linear dependence to z L D , will result in underestimation of the chromophore depth.We believe this effect is responsible for the observed discrepancy of depth measurements in our experiments with epoxy phantoms.
An important consideration in diagnostic imaging using photothermal techniques for imaging in vivo blood vessels in human patients is the peak incident laser fluence.Large peak laser fluences can result in irreversible thermal denaturation of proteins that comprise important structural components in the tissue.Although pulsed photothermal techniques reportedly utilize a sub-therapeutic light dose, typical peak laser fluences are 10 kW cm −2 .For comparison, in the rodent dorsal skin flap window model experiments described above, a peak laser fluence of 0.38 W cm −2 was applied for thermal wave imaging.The substantial reduction in peak laser fluence associated with the thermal wave approach is an important advantage when considering diagnostic imaging of in vivo blood vessels in human patients.

Conclusions
We have demonstrated that coherent thermal wave imaging provides a significant improvement in the quality of the recorded IR images of subsurface chromophores in biological materials.Our technique utilizes a sequence of two-dimensional IR images of a test specimen and a software-implemented coherent signal processing algorithm to determine amplitude and phase of the thermal wave radiometric image which may be used to determine physical dimensions and depth of subsurface chromophores in tissue.A numerical model was developed to simulate radiometric response due to periodic laser excitation and known distribution of a volumetric heat source.Our numerical simulations indicate that optical and thermal properties of the test specimen can strongly influence the detected phase response and should be taken into account for a correct interpretation and analysis of the recorded data.The substantial reduction in the peak laser fluence obtained using the thermal wave methodology is an important advantage for diagnostic imaging in vivo blood vessels in human patients.

Figure 1 .
Figure 1.(A) Schematic diagram of experimental apparatus for thermal wave imaging of blood vessels in biological materials.Laser excitation waveform (B) and radiometric response (C) of a single pixel from the recorded sequence of IR images.

Figure 2 .
Figure 2. Block-diagram of a software module implementing coherent signal processing algorithm applied to each pixel (i, j) in recorded IR frames.

Figure 3 .Figure 4 .
Figure 3. Thermal wave imaging of the epoxy phantom with embedded chromophores.(A) Original image from a sequence of 1024 recorded IR frames, M is a black paint marker on the phantom surface, (B) reconstructed images of amplitude (A ij ) and (C) phase (θ ij ) at laser modulation frequency f m = 0.5 Hz.Cylindrically shaped subsurface chromophores 1 and 2 are oriented obliquely to the phantom surface and positioned at depths between 250 and 500 µm.Scale bar is 1 mm.
After coherent processing of the image sequence at f m = 1.0 Hz, structure of the CAM vascular network (figures 4(B) and (C)) is visible.Although the amplitude image (figure 4(B)) shows capillaries as small as 150-200 µm, finite thermal diffusion length at 1.0 Hz (L D = 187 µm) gives blood vessels a blurred appearance.Although use of higher laser excitation frequencies renders sharper images of superficial capillaries, signal from deeper blood vessels is strongly attenuated and they disappear from the image.

Figure 5 .
Figure 5. Subdermal view (A) of the rodent dorsal skin flap model (window diameter is 12 mm, square indicates IR imaging area).Reconstructed amplitude (A ij ) (B) and phase (θ ij ) (C) of thermal wave radiometric image at laser modulation frequency f m = 1.0 Hz.Thermal wave images are inverted (left-right) so blood vessels appear in the same orientation as the visible subdermal view.

Figure 6 .
Figure 6.Numerical simulation of the thermal wave response due to periodic heating of subsurface chromophore in a weakly absorbing tissue.Rectangular parallelepiped-shaped chromophore (A) with dimensions 0.2 × 3.0 × 0.1 mm 3 positioned d = 300 µm deep below air tissue interface.Distribution of heat generation rate in tissue along z-axis passing through the chromophore (B).

Figure 7 .
Figure 7. Simulated images of amplitude (A ij ) (A) and phase (θ ij ) (B) computed for the subsurface chromophore shown in figure 6(A) and laser modulation frequency f m = 1.0 Hz.Scale bar is 1 mm.

Figure 8 .
Figure 8. Phase profile (A) across the centre of a simulated image shown in figure 7(B) and phase profile (B) without background absorption (α 0 = 0).
L D define thermal wave phase response.For laser modulation frequency f m = 1.0 Hz and thermal diffusivity D = 0.11 mm 2 s −1 , the thermal diffusion length L D = 187 µm.Phase response of superficial chromophores (d < L D ) in the transversal plane has a minimum at the