UNIVERSITY OF LJUBLJANA FACULTY OF MATHEMATICS AND PHYSICS DEPARTMENT OF PHYSICS Matija Milanič DEVELOPMENT AND EVALUATION OF PULSED PHOTOTHERMAL RADIOMETRY FOR TEMPERATURE PROFILING IN BIOLOGICAL TISSUES Doctoral thesis ADVISER: Assist. Prof. Boris Majaron Ljubljana, 2008 UNIVERZA V LJUBLJANI FAKULTETA ZA MATEMATIKO IN FIZIKO ODDELEK ZA FIZIKO Matija Milanič RAZVOJ IN EVALUACIJA METODE SUNKOVNE FOTOTERMALNE RADIOMETRIJE ZA TEMPERATURNO PROFILOMETRIJO BIOLOŠKIH TKIV Doktorska disertacija MENTOR: doc. dr. Boris Majaron Ljubljana, 2008 Abstract In this dissertation we present the development and evaluation of pulsed photothermal radiometric (PPTR) temperature depth profiling in biological tissues. Motivation for this work is incomplete success in laser therapy of port-wine stain birthmarks (PWS). PPTR technique, which utilizes infrared (IR) emission from materials following pulsed laser exposure, can provide information about PWS depth and epidermal thickness, both required to optimize the therapy. We develop three original reconstruction codes, which are based on truncated singular value decomposition (TSVD), conjugated gradient (CG) and the ?-method, respectively. All codes involve a non-negativity constraint to the sought temperature vector and automatic regularization. When applied to different test objects, all three codes produce reconstruction results, which are much more accurate than the results published in earlier PPTR studies. Calibration of PPTR signals and the error due to linearization of PPTR signal expression is analyzed. We find that the linearization error depends on temperature amplitude, absorber depth, acquisition time, and the spectral acquisition band. PPTR measurements commonly employ broad-band signal acquisition to increase signal-to-noise ratio (SNR), but all reported studies use a fixed effective IR absorption coefficient (µeff). We show that in samples with large spectral variation of µ(?) in mid-IR, which includes most biological tissues, selection of µeff strongly affects the accuracy of the results. A novel analytical approach to determination of optimal µeff from spectral properties of the sample and radiation detector is presented. In extensive numerical simulation of PPTR temperature profiling in human skin using different IR detectors and spectral bands, we demonstrate that our approach predicts viable values of µeff. The influence of spectral filtering on the accuracy of temperature profiles is studied by a systematic experimental comparison of PPTR temperature profiling on agar tissue phantoms utilizing the customary spectral band of the InSb detector (? = 3.0-5.6 µm) and a narrowed acquisition spectral band (? = 4.5-5.6 µm). To support our experimental observations, we present also a detailed numerical simulation of the experimental procedure utilizing spectral acquisition bands with the lower wavelength limit varied between ?l = 3.0 µm and 5.0 µm and the upper wavelength limit fixed at 5.6 µm. The experimental and numerical simulation results indicate that spectral filtering reduces reconstruction error and broadening of temperature profiles, especially for shallower and more complex absorbing structures. Analogously, we performed the experiment and numerical simulation involving gelatin tissue phantoms, which more closely resemble human skin. Again, we find that a suitable spectral filtering (?l = 4.0–4.5 µm) is beneficial, despite the associated reduction of SNR. We determine experimentally the accuracy of PPTR temperature depth profiling in custom tissue phantoms composed of agar gel layers separated by single very thin absorbing layers. The laser-induced temperature depth profiles, reconstructed from measured PPTR signals, correlate very well with absorber depths determined by magnetic resonance imaging and optical microscopy. We observe significant broadening and attenuation of reconstructed profiles with increasing depth of absorbing layer. Corresponding numerical analysis indicates that the broadening equals to ~13% of the absorber depth. Using a numerical simulation we also analyzed, how the accuracy of reconstructed temperature profiles depends on sampling frequency. PACS: 02.30.Zz, 42.30.Wb, 42.62.Be, 87.57.Gg, 87.63.it, 87.64.Cc, 87.50.W-, 87.61.Ff Keywords: pulsed photothermal radiometry (PPTR), temperature depth profiling, image reconstruction, infrared absorption Povzetek V disertaciji je predstavljen razvoj in evalvacija sunkovne fototermalne radiometrije (SFTR) za temperaturno profilometrijo v bioloških tkivih. Motivacija za to delo je nepopoln uspeh laserske terapije ognjenih znamenj (OZ). SFTR tehnika, ki izrablja infrardeče (IR) sevanje snovi pogrete zaradi laserskega sunka, omogoča določitev globine OZ in debeline epidermisa, podatkov pomembnih za optimizacijo zdravljenja OZ. Zato smo razvili rekonstrukcijske algoritme, ki omogočajo reševanje inverznega problema SFTR. Algoritmi temeljijo na metodah dekompozicije po singularnih vrednostih (TSVD), konjugiranih gradientov (CG) in ?-metodi. V vse algoritme smo vključili pogoj nenegativnosti in avtomatsko regularizacijo. Ko smo jih preizkusili na različnih testnih objektih, vsi trije algoritmi vrnejo rekonstruirane rešitve, ki so precej natančnejše od rezultatov objavljenih v predhodnih študijah SFTR. Nato smo analizirali kalibracijo SFTR signalov in napako linearizacije izraza SFTR signala. Ugotovili smo, da je napaka odvisna od amplitude in globine temperaturnih profilov, časa zajemanja in spektralnega pasu zajemanja. SFTR meritve se običajno opravlja v širokem spektralnem pasu, saj se s tem izboljša razmerje med signalom in šumom. Vendar pa se hkrati uporablja konstantna vrednost IR absorpcijskega koeficienta µeff. Pokazali smo, da pri vzorcih z veliko variacijo IR spektra µ(?), kar vključuje večino bioloških vzorcev, izbira µeff bistveno vpliva na kvaliteto rekonstrukcije. Predstavili smo nov pristop analitičnega določanja optimalne vrednosti µeff, ki vključuje spekter µ(?) vzorca in spektralno občutljivost IR detektorja. Z obsežno numerično simulacijo smo pokazali, da naš analitični pristop da optimalne vrednosti µeff. Vpliv spektralnega filtriranja na točnost rekonstruiranih temperaturnih profilov smo proučili s sistematično eksperimentalno primerjavo temperaturne profilometrije SFTR v tkivnih fantomih iz agarja, pri čemer smo uporabili celoten (? = 3,0–5,6 µm) in zožen spektralni pas (? = 4,5–5,6 µm) detektorja. V podporo eksperimentalnim rezultatom so izvedli tudi numerično simulacijo eksperimenta za različne spektralne pasove med ? = 3,0 in 5,6 µm. Tako eksperimentalni kot numerični rezultati nakazujejo, da zoženje spektralnega pasu zmanjša napako rekonstrukcije in razširitev rekonstruiranih temperaturnih profilov, še posebej za plitve in kompleksnejše absorbirajoče strukture. Izvedli smo tudi eksperiment in numerično simulacijo za kolagenske tkivne fantome, ki posnemajo lastnosti človeške kože bolje kot agar. Tudi v tem primeru smo ugotovili, da spektralno filtriranje izboljša rezultate rekonstrukcije, navkljub zmanjšanem razmerju signal-šum. Eksperimentalno smo določili natančnost našega sistema SFTR na fantomih iz agarja, ki so vsebovali eno tanko plast absorberja na različnih globinah. Globine rekonstruiranih temperaturnih profilov se zelo dobro ujemajo z globinami, ki smo jih določili z optično mikroskopijo ali mikro magnetno resonančnim slikanjem. Očitna pa je razširitev in atenuacija temperaturnih profilov z naraščajočo globino absorbirajoče plasti. Izvedli smo tudi numerično simulacijo in na podlagi eksperimentalnih in numeričnig rezultatov ugotovili, da efekt razširitve znaša do 15% globine absorberja. PACS: 02.30.Zz, 42.30.Wb, 42.62.Be, 87.57.Gg, 87.63.it, 87.64.Cc, 87.50.W-, 87.61.Ff Ključne besede: sunkovna fototermalna radiometrija (SFTR), temperaturna globinska profilometrija, rekonstrukcija slike, infrardeča absorpcija I would like to sincerely thank my advisor assist. prof. Boris Majaron for his guiding, discussions and support. I also thank prof. dr. J. Stuart Nelson from Department of Biomedical Engineering, University of California, Irvine, for inviting me to his group where I spent a fruitful month experimenting on PPTR temperature profilometry and OCT imaging. I am very grateful to assist. prof. Igor Serša for help with micro-MRI imaging, and to Nadja Derviševič who helped me with tissue phantoms and PPTR measurements. I also thank all my colleagues for all the fun I had during my time at the Complex Matter Department Finally, I owe special gratitude to my family for continuous and unconditional support and to Nina for her enduring patience, understanding, and love. Contents 1 Introduction ....................................................................................................................................... 13 2 Theoretical background .................................................................................................................... 15 3 Experimental setup ........................................................................................................................... 17 3.1 Pulsed laser source .............................................................................................................. 17 3.2 IR radiation detector ............................................................................................................ 18 3.3 Infrared optics ..................................................................................................................... 20 3.4 Radiometric signal ............................................................................................................... 21 3.5 Experimental characterization of our PPTR system ............................................................ 21 4 Reconstruction algorithms ............................................................................................................... 25 4.1 Discrete ill-posed problems ................................................................................................. 25 4.2 Reconstruction algorithms ................................................................................................... 27 4.3 Regularization ..................................................................................................................... 29 4.4 Non-negativity constraint .................................................................................................... 31 4.5 Performance of constrained reconstruction algorithms ....................................................... 33 4.6 Conclusions ......................................................................................................................... 37 5 Calibration of radiometric signal ..................................................................................................... 39 5.1 Linearization error in monochromatic approximation ........................................................ 39 5.2 Broad-band signal acquisition ............................................................................................. 43 5.3 Conclusions ......................................................................................................................... 45 6 Effective infrared absorption coefficient ......................................................................................... 47 6.1 InSb detector ....................................................................................................................... 47 6.2 HgCdTe detectors ............................................................................................................... 52 6.3 Discussion ........................................................................................................................... 53 6.4 Conclusions ......................................................................................................................... 54 7 Sampling rate ..................................................................................................................................... 55 7.1 Methods ................................................................................................................................ 55 7.2 Results .................................................................................................................................. 56 7.3 Discussion ............................................................................................................................ 60 7.4 Conclusions .......................................................................................................................... 60 8 Tissue phantoms ................................................................................................................................ 61 8.1 Hydrogel layer ..................................................................................................................... 61 8.2 Absorbing layer ................................................................................................................... 63 9 Spectral filtering ................................................................................................................................ 65 9.1 Experiments in agar tissue phantoms .................................................................................. 65 9.2 Numerical simulation ......................................................................................................... 70 9.3 Discussion ........................................................................................................................... 75 9.4 Conclusions ........................................................................................................................ 76 10 Spectral filtering in collagen samples ........................................................................................... 77 10.1 Numerical simulation ....................................................................................................... 77 10.2 Experiment ....................................................................................................................... 85 10.3 Discussion ......................................................................................................................... 88 10.4 Conclusions ...................................................................................................................... 88 11 Accuracy of temperature profiling ............................................................................................... 89 11.1 Experiments on agar tissue phantoms .............................................................................. 89 11.2 Numerical simulation ....................................................................................................... 93 11.3 Discussion ......................................................................................................................... 94 11.4 Conclusions ...................................................................................................................... 94 12 Summary and conclusions ............................................................................................................. 94 A Filter factors .................................................................................................................................... 94 B Chebyshev approximation for erfcx(x) ......................................................................................... 94 C Povzetek disertacije v slovenščini .................................................................................................. 94 C.1 Teoretično ozadje ............................................................................................................... 94 C.2 Postavitev eksperimenta ..................................................................................................... 94 C.3 Rekonstrukcijski algoritmi ................................................................................................. 94 C.4 Kalibracija .......................................................................................................................... 94 C.5 Efektivni IR absorpcijski koeficient .................................................................................. 94 C.6 Frekvenca zajemanja .......................................................................................................... 94 C.7 Tkivni fantomi ................................................................................................................... 94 C.8 Spektralno filtriranje .......................................................................................................... 94 C.9 Spektralno filtriranje v kolagenskih vzorcih ...................................................................... 94 C.10 Natančnost temperaturne profilometrije SFTR .................................................................. 94 Bibliography ......................................................................................................................................... 94 Chapter 1 Introduction Nowadays, cooperation between physics and medicine has become tighter than ever before. Physics provides medicine with new explanations, tools for diagnostics and treatment, while medicine presents physics with challenges from a completely new world. There are countless areas of cooperation; just to name some of them: x-rays, MRI, biomedical optics.1 The advent of lasers had a great impact on modern surgery, diagnostics of eye diseases, cancer, etc. In the dissertation we are presenting the development and evaluation of a radiometric technique known as pulsed photothermal radiometry (PPTR). The motivation for this work is incomplete success in laser therapy of port-wine stain (PWS) birthmarks. PWS are permanent hypervascular lesions in human skin, which consist of an excess of ectacic blood vessels. These are usually fully contained within the most superficial millimeter of the skin. The exact depth varies from patient to patient, but on average, the highest fractional blood content is found 200–400 µm below the epidermal-dermal junction.2 PWS are currently treated by selective photocoagulation of the ectatic vasculature using pulsed green or yellow/orange laser. In order to optimize laser therapy parameters (pulse duration, radiant exposure, and light wavelength) on an individual patient basis, determination of PWS depth, epidermal thickness, and epidermal heating is required.3,4 PPTR is a noncontact technique, which utilizes infrared (IR) emission from materials following pulsed laser exposure. Selective absorption of laser radiation in subsurface cromophores results in localized heating and may be detected as transient increase in IR emission from tissue surface. When thermal properties of the sample are known, the laser-induced temperature profile can be reconstructed from acquired radiometric signals. Such PPTR temperature profiling was recognized as a promising technique for non-invasive determination of structure and cromophore distribution in strongly scattering biological tissues and was extensively investigated.5–12 In contrast to alternative diagnostic techniques, which utilize detection of scattered or frequency converted light from the irradiated tissue (e.g., optical coherence tomography,13 diffuse optical tomography,14 ultrasound-modulated optical tomography15), PPTR provides a signal amplitude that is directly related to the initial space-dependent temperature increase in targeted cromophores. An important advantage of PPTR over the photoacustic technique16 is that the measurement is made without touching the sample. For example, resting a photoacustic probe on the skin surface for several seconds may change the hydratation of stratum corneum, while slightly heavier contact will change the hemodynamics in the area under the probe. In the dissertation we present development and evaluation of a PPTR system for temperature profiling in biological tissues performed at the Complex Matter Department at the Jožef Stefan Institute in Ljubljana, Slovenia. The work involved development of novel reconstruction algorithms, numerical analysis of PPTR profiling, construction and optimization of a laboratory PPTR system and experiments on tissue phantoms. The dissertation is structured as follows. The PPTR temperature profiling inverse problem is derived in Chapter 2. Chapter 3 describes the actual PPTR temperature profiling setup. Expressions for radiometric signal and measurement noise are derived, which are then used in experimental characterization of our PPTR temperature profiling system. In Chapter 4 we develop original reconstruction codes for solving the PPTR temperature profiling inverse problem. The calibration procedure and error estimation due to simplifications of PPTR signal expression are presented in Chapter 5. Chapter 6 presents the effect of the monochromatic approximation deficiency on 13 14 Chapter 1 Introduction reconstructed temperature profile for different spectral acquisition bands, and an analytical approach for determination of the optimal effective absorption coefficient. Chapter 7 analyzes the influence of the sampling rate. Second part of the dissertation presents experimental evaluation of the PPTR system performance. We begin with construction of tissue phantoms, which serve as test samples in our experiments (Chapter 8). The experiments and numerical simulations on agar (Chapter 9) and collagen gel tissue phantoms (Chapter 10) demonstrate the effect of spectral filtering on reconstructed temperature profiles. In Chapter 11 we analyze the accuracy of reconstructed temperature profiles, specifically the depth and width of reconstructed temperature profiles. In Chapter 12 we summarize the main ideas and conclusions. Chapter 2 Theoretical background The mechanism of PPTR signal generation was given in early work by Leung and Tam.17 Basic relations for PPTR temperature depth profiling were first derived in one dimension and monochromatic approximation,6–9 and later on extended to account for spectral variation of the sample IR absorption coefficient µ(X) in the mid-IR detection range.18 The one-dimensional theory was also extended to three-dimensions.19 We derive an expression for PPTR signal amplitude AS(t) in terms of the initial temperature depth profile AT(z, 0) in a tissue immediately following pulsed laser irradiation. For the purpose of our analysis, we assume that the tissue occupies a semi-infinite half-space. We use the one-dimensional heat equation for the temperature increase AT(z, t)20 ^\2 A T ^ A T D----- +---------QAT = 0 (2.1) dz 2 dt with a mixed boundary condition at the air-tissue interface dAT haT -A\ =0 (2.2) z=0 z=0 dz where D denotes the thermal diffusivity of tissue, which we assume is homogenous, h represents the reduced heat transfer coefficient at the air-tissue interface. In addition to heat diffusion, heat can be removed from the irradiated tissue also by the blood flow. The effect of blood flow is included in (2.2) as the blood perfusion rate Q. We do not include heat loss due to radiation in (2.1), because it is significantly smaller than heat conduction for typical experimental conditions. The Green’s function solution to (2.1)-(2.2)21 represents temperature increase in the tissue at depth z and time t in response to instantaneous release of a planar impulse heat source at depth z ’ and time t = 0 G t (z' z , t ) = 4^Dt e (W) /4Dt+e^z) /4Dt [1-h/4^erfcx(u)] (2.3) where z + z ' + hDt and erfcx(u) = exp(u2) erfc(u), where erfc(u) is the complementary error function. Because the relaxation time due to blood perfusion (1/Q ) in the microvasculature is much longer than the time of measurement, the exponent function involving Q can be neglected from (2.3). Using the Green's function solution, the temperature depth profile AT(z, t) at an arbitrary depth z and time t is written as (•CO AT(z,t)=\ AT(z,0)GT(z',z, t)dz (2.4) Jz'=0 The regular expression for radiometric signal of blackbody involving constant temperature Tb, can be generalized to correspond to objects with non-homogenous temperature distribution by including z-integral with factor µ(X) exp(– µ(X) z). Thus, the measured radiometric signal S(t ) is then given by 15 16 Chapter 2 Theoretical background integrating Planck's expression for radiative emission BK(Tb + AT(z, t )) over all depths \ oo S(t ) = C h r(X) ju(A) [ BÄ[Tb + AT(z, t)] e'KX)z dzdX (2.5) \ z=0 where h and K are the lower and upper limit of the detection spectral range, respectively, while R(X) describes spectral sensitivity of the radiation detector. The constant C accounts for sample emissivity and other experimental specifics (e.g., losses of collection optics, radiation detector field of view, etc.). However, we obtain the regular expression for radiometric signal of blackbody by setting AT(z, t ) = 0 and integrating the z-integral. When the induced temperature rise AT(z, t) is significantly smaller than Tb, we can expand Bk(Tb) in Taylor series, which leads to the linearized expression for the transient part of the spectrally composite radiometric signal \ CO AS(t ) = C h R(Ä) B'x (Tb) ju(A) [ AT(z, t) e^(X)zdzdX (2.6) \ z=0 where Bx'(Tb) represents the temperature derivative of Bk(Tb). From (2.4) the PPTR signal is related to the laser-induced temperature profile AT(z,0) by a simple convolution CO AS(t ) = [ K(z,t)AT(z,0)dz (2.7) z=0 with the kernel function K(z,t) defined by (2.3)-(2.6): \ CO K(z,t) = C h R(X)B^(Tb) /u(Ä) JGT(z',z,t ) e^Wz'dz' dÄ (2.8) When using a spectrally invariant value µ is justified (or assumed), K(z,t) can be simplified through factorization of the double integral in (2.8): jj CO K(z,t) = C h R(Ä)B^(Tb)dA ju\GT(z',z,t)e^dz' (2.9) where, for given experimental conditions, the first integral yields a constant. This relation allows us to establish a direct correspondence with earlier reports on PPTR depth profiling, disregarding the spectral variation µ^). By inserting GT(z’, z, t), the second integral results in6 K(z,t) = -exp[-z 2 /(4D t )] ]erfcx(u ) + erfcx(u+) +—[erfcx( u )-erfcx( u 1)] (2.10) 2 [ /u-h with erfcx(u) = [1 - erf(u)] exp(u2), u± = HDt + z/(2 Dt), u1 = hDt + z/(2 Dt). In experimental practice, PPTR signals are represented by vectors, and relation (2.7) becomes multiplication of the initial temperature profile vector T (Tj = AT(zj,0)) with kernel matrix K (Ki,j= S = KT (2.11) 0 0 Chapter 3 Experimental setup Figure 3.1 shows a schematic of an experimental PPTR system. The test object is irradiated with a single pulse from a pulsed laser. The irradiation spot size must be large (typical diameter ~ 5 mm) as compared to the surface of the studied volume of the sample (typical diameter ~ 1 mm), because one-dimensional analysis is used. IR emission is monitored at normal incidence to the sample surface by a radiation detector. IR radiation is collected on the detector by an IR collection lens. The electrical response of the detector, after preamplification, is monitored by a fast analog-to-digital converter. A silicon photodiode is used to detect the laser pulse onset. Digitized signals are stored and processed in a personal computer. ruiocu idoci Signal acquisition and processing Figure 3.1: Experimental setup for PPTR temperature profiling in biological tissues. 3.1 Pulsed laser source Any pulsed light source, which is selectively absorbed in the sample, can be used to induce the initial temperature profile. For PPTR temperature profiling in human skin, we most often apply the same light sources as used for laser therapy: pulsed green (KTP/YAG, ? = 532 nm) or yellow/orange lasers (flashlamp-pumped dye, 577 nm and 585-600 nm). Since the above wavelengths are preferentially absorbed in melanin and hemoglobin, which are located in epidermis and blood vessels, we obtain information about epidermal thickness and blood vessel distribution. Furthermore, temperature profiles similar to those induced during laser therapy are obtained when the above therapeutic lasers are used. The required pulse energies depend on the irradiation spot size, detector sensitivity and cromophore absorption. Clearly, it is desirable to avoid any thermal damage. Pulse length ? should ideally be about 1 ms or less to prevent excessive heat diffusion from induced temperature profiles. Table 3.1 lists some combinations of excitation wavelength ?, radiant exposure H and pulse length ? used for PPTR temperature profiling in human skin. In our experiments we used a pulsed-dye laser (PDL) at wavelength 585 nm with pulse length 1.5 ms (ScleroPlus, Candela, Wayland, MA, USA) and a KTP laser generating 1 ms long 532 nm laser pulses (DualisVP, Fotona, Ljubljana). 17 18 Chapter 3 Experimental setup Table 3.1: Some combinations of excitation wavelength ?, radiant exposure H and pulse length ? used for PPTR temperature profiling in human skin. PDL – pulsed dye laser, KTP – second harmonic Nd:YAG laser. Laser type A(nm) H (J/cm2) r(ms) 0.45 References PDL 585 7.0 6 PDL 577 0.5 0.001 9 KTP 532 3.4 2-10 19 PDL 585-600 5.0-6.0 1.5 18 3.2 IR radiation detector Commercially available IR radiation detectors are divided into two groups: thermal and quantum detectors. The former are not applicable to PPTR temperature profiling, because they have a low signal-to-noise (SNR) ratios and they are not convenient for measurement of transient phenomena. Common IR detectors involved in PPTR temperature profiling are liquid nitrogen cooled photovoltaic InSb and photoconductive HgCdTe (MCT) quantum detectors. The main disadvantage of the photoconductive detector is that it requires a mechanical chopper to modulate IR radiation. The chopper limits the sampling rate to a few 1000 s-1, introduces additional measurement noise to the radiometric signal, requires a lock-in amplifier and generates air current, which increases the heat loss at the sample surface (h). Furthermore, InSb detector offers a more predictable responsivity,22 which is especially favorable for numerical simulations. In contrast, MCT detector yields larger PPTR signals, since wider spectral acquisition bands can be used, and the temperature derivative of Planck’s formula B’(Tb) is larger for these spectral bands (see Fig. 3.3, dashed line). Considering the advantages and weaknesses of both detectors, we have decided to use the InSb detector in our experimental setup. Figure 3.2: Current mode preamplifier circuit for InSb detector. ? represents IR radiation, is signal current, Rf feedback resistor, and Vs signal voltage. Photovoltaic detectors are diodes made from semiconductor materials. InSb material typically has a band gap of 0.22 eV at temperature T = 77 K, thus only photons with ? ? 5.5 µm are detected. This maximal detectable wavelength is called the cut-off wavelength (?c). When the diode is exposed to IR radiation, it generates a current proportional to the photon arrival rate. Photovoltaic detectors are commonly used in a current preamplifier circuit (Fig. 3.2), where voltage drop across the detector is practically zero. When monochromatic radiation of wavelength ? and power P irradiates the photovoltaic detector, the signal current is equals22 hc P (3.1) 3.2 IR radiation detector 19 where ? denotes quantum efficiency (about 0.64 for InSb), e0 is the electron charge, h represents the Planck’s constant and c the speed of light. The spectral responsivity of a radiation detector R(?) is defined as the signal output is divided by the radiant input P, thus the theoretical spectral responsivity of InSb detector is easily deduced from (3.1) R(Ä) = ?je0 X hc (3.2) Figure 3.3 shows relative spectral responsivity R(?)/Rp for a specific InSb detector (P5968-100, Hamamatsu) with the peak responsivity Rp = 2.5 A/W at ?p = 5.3 µm and a specific HgCdTe detector (P3257, Hamamatsu) with the peak responsivity at ?p = 10 µm. 0.5 ?p (µm) Figure 3.3: Relative spectral responsivity R(?) of InSb and HgCdTe (MCT) detectors (solid lines) with peak responsivities at ?p = 5.3 µm and 10.0 µm, respectively. Temperature derivative of Planck’s radiation formula B?’(Tb) at Tb = 303 K (dashed line) has the maximum at 8 µm. 3.2.1 Measurement noise Noise present in an IR detection system originates from a number of sources.2223 Shot noise originates from the discrete nature of photodetection process. Its amplitude is (3.3) where ?f represents frequency bandwidth. Frequency bandwidth is calculated as 1/(2tint), where tint is the integration time. Johnson noise originates from random motion of free charges. Its amplitude is 4 kBTdAf R (3.4) f where kB is the Boltzmann’s constant, Td represents temperature of the detector, and Rf is the feedback resistance (see Fig. 3.2). Other contributions to measurement noise include amplifier noise, digitalization noise and external noises. Inasmuch as multiple noise sources are present, we can use the central-limit theorem to estimate square of the total noise amplitude nt as the sum of squares of amplitudes of all noise contributions. While most noise contributions are spectrally invariant (“white”), so-called 1/f noise is often present in radiometric signals. In general, the presence of 1/f noise is characterized by the corner frequency fc and exponent ?. The total noise spectral density n%t is given by22 1.0 0.0 nJ = 20 Chapter 3 Experimental setup n = n white f\a 1 yf) (3.5) where n%white denotes the spectral density of white noise and f represents frequency. For discrete radiometric signals, the total noise amplitude nt and the noise spectral density n%t are correlated using the discrete Parseval’s theorem24 N i 1 S n t2(fi) N2 (3.6) where fi represent the discrete frequencies and N denotes the number of radiometric signal and spectral density data points. 3.3 Infrared optics Figure 3.4 presents the transmittance T(?) of three materials commonly utilized for mid-IR optics (? = 3–5 µm). 10 0 90 80 70 60 50 /^ ^ CaF2 ' Si (BBAR) /\ Sapphire ^v 2 3 4 5 /l (um) Figure 3.4: Transmittance of broadband antireflex coated silicon (MellesGriot), and CaF2 (ThorLabs). 9 10 (Tydex, St. Petersburg, Russia), sapphire In our setup we use the broadband antireflex silicon IR optics because of large transmittance in the k = 3-6 µm spectral band and good mechanical characteristics. The collection optics consists of two planoconvex silicon antireflex coated lenses (Galvoptics, Essex, UK) with transmittance T > 98% at the 3-5 µm spectral band. The lenses are positioned in such a manner that magnification is M = 1 (Fig. 3.5). n t = Figure 3.5: Schematics of the 2-lens collection optics with magnification M = 1. The source area As equals the detector active area A, and the collection angle at the source 0s matches the collection angle 0 at the IR detector. Arrows indicate direction of IR radiation. 3.4 Radiometric signal 21 3.4 Radiometric signal Figure 3.5 presents a schematic of radiometric signal collection using the 2-lens collection optics with M = 1. The source area As equals the IR detector active area A and the collection angle at the source 0s matches the collection angle 0 of the IR detector. Expression for 0 follows from Fig. 3.5 tan9 = — (3.7) f where r is the lens radius (i.e., aperture) and f represents its focal length. The amount of radiant power emitted from the source area As at temperature Tb at wavelength k and which falls within the solid angle dQ around the direction specified by 0s is dPÄ = —BÄ(Tb)cos9s As dQ (3.8) where e denotes the sample emissivity. The radiation power varies as the cosine of 0s (Lambert’s law). The total radiation power Pk emitted into the solid angle 0s is found by integrating (3.8) PÄ = ssin2ßsAs BÄ(Tb) (3.9) In absences of transmission lenses, the emitted power equals the collected power on the detector. However, the power must be reduced for the collection optic losses. Hence, the analogous expression to (3.8) for the collected radiation power is obtained PÄ = sT(Ä) sin29 A BÄ(Tb) (3.10) Accordingly, the radiometric signal is is detected in a broad spectral band between k and Ah K is = esin29A h T(X) R(X) BÄ(Tb) dÄ (3.11) \ 3.5 Experimental characterization of our PPTR system In our system, IR radiation is detected by a single-element InSb detector (P5968-100, Hamamatsu). The detector properties, as specified by the manufacturer, are active area A = 0.79 mm2 (d = 1 mm), half-angle of the field of view 0 = 22.5°, and spectral responsivity R(X) with peak value Rp = 2.5 A/W at Ap = 5.3 µm (Fig. 3.3), while the cut-off wavelength is Xc = 5.6 µm. In addition, the manufacturer has specified measurement noise under specific conditions, where small radiometric signal from a 500 K blackbody is detected in a background-limited regime (BLIP) at modulation frequency of 1200 Hz, and with a small frequency bandwidth Af (1 Hz). The peak detectivity Dp*, as a standard figure of merit22 can be calculated from the reported noise standard deviation and the above detector parameters, resulting in D* = 2.8 × 1011 cm Hz1/2/W. This value does not include 1/f noise, which may dominate at lower frequencies, and other noise sources (i.e., amplifier noise, external noise sources), which may be present in an actual PPTR system. Therefore, we have determined experimentally the noise parameters for our PPTR system. We measured the response of the IR detector using a blackbody (BB701, Omega Engineering, Stamford, CT, USA) set to different temperatures, TBB = 288-232 K. The entire detector's field of view was filled with radiation from the blackbody. The detector response at each TBB was acquired using the entire spectral band of the IR detector (X = 3.0-5.6 µm). The radiometric values were acquired at a sampling rate of 50,000 s1. We have numerically reduced the sampling rate by calculating the average value of 50 consecutive values (t int = 1 ms). This yielded radiometric signals 22 Chapter 3 Experimental setup s(TBB) with the sampling rate f = 1000 s-1 and frequency bandwidth ?f = 500 s-1. We determined theoretical values of s(TBB) by inserting the blackbody temperature TBB and the system parameters into (3.11) and completing the ? integral. The blackbody emissivity was ? = 0.94 and transmittance T = 0.98 (both specified by the manufacturers). Figure 3.6 presents the average of experimental radiometric signals (circles), which matches perfectly the theoretically predicted signal values (line); thus (3.11) adequately predicts the PPTR system response when the above R(?) is used. 7.0x10-6 6.0x10 5.0x10 4.0x10 3.0x10" 2.0x10 290 300 310 T (K) 320 Figure 3.6: Experimental (circles) and theoretical (line) radiometric signal values is(TBB) determined for TBB = 288–323 K and for the entire spectral band. To obtain the noise spectrum each determined signal s(TBB) was Fourier transformed using the fft function implemented in Matlab 14 (Mathworks, Natick, MA, USA). We fit (3.5) to the noise spectrum to determine the noise parameters fc and ?. The total noise amplitude nt was determined as the standard deviation of s(TBB), albeit it could be found also from the noise spectrum using the discrete Parseval’s theorem (Eq. 3.6). An example of a noise spectrum and the corresponding fit of (3.7) is presented in Fig. 3.7. The 1/f noise evidently dominates at frequencies below ~20 s-1. 10 10 10 f (S ) Figure 3.7: An example of a noise spectrum measured at TBB = 303 K using the entire spectral band (black line) and the best fit of function n2t (Eq. 3.5) (gray line). The determined noise parameters are total noise amplitude nt = 3.2×1010 A, corner frequency fc = 15 s-1, and exponent a = 1.5. In order to accurately discern between detector-specific noise nd and shot noise nsh, the detector response at each TBB was acquired also using a long-pass IR filter (Barr Associates, Westford, MA; see Fig. 3.8) placed between the blackbody and the detector. 3.5 Experimental characterization of our PPTR system 23 1.0 0.5 0.0 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Figure 3.8: Transmittance of a long-pass IR filter with cut-on wavelength of ? = 5 µm (Barr Associates, Westford, MA). Then, we calculated shot noise amplitudes nsh for each temperature TBB and both spectral acquisition bands by substituting the measured radiometric signals s(TBB) into (3.3). Figure 3.9 presents the experimental noise amplitudes nt and shot noise amplitudes nsh as a function of TBB. For the entire spectral band (Fig. 3.9a), values nt are scattered around the average value 3.1 × 10-10 A with a standard deviation of ~10-11 A. Shot noise amplitudes nsh increase monotonically with increasing TBB, but are approximately 10 times smaller than nt. For the reduced spectral band (Fig. 3.9b), nt is (2.9 ± 0.1) × 10-10 A. Since detector-specific noise obviously dominates in our PPTR system and the average nt is almost identical for both spectral bands, a constant value nt = 3 × 10-10 A suitably characterizes the total noise in our PPTR signals at ?f = 500 s-1. 3.2x10" 3.0x10 2.8x10 2.DX10 4.0x10 - I 2.0x10 0.0 a O O O o O 0 O O O nt A n sh A A A A A A A A 3.2x10" 3.0x10 2.8x10 2.DX10 - I 4.0x10 2.0x10 0.0 A A A A A 290 320 300 310 320 290 300 310 TBB (K) TBB (K) Figure 3.9: Determined total noise amplitude nt (circles) and shot noise amplitude nsh (triangles) as functions of TBB for (a) the entire spectral band and (b) the reduced spectral band. Radiometric signals measured by our PPTR system are dominated by the 1/f-noise at low frequencies (Fig. 3.7). The average corner frequency is fc = 15 ± 10 s-1 and 10 ± 8 s-1 for the entire and reduced spectral acquisition band, respectively. The average exponent is ? = 1.3 ± 0.4 for both spectral acquisition bands. D o n a n Sh Chapter 4 Reconstruction algorithms Ill-posed inverse problems arise quite often when one is interested in determining the internal structure of a physical system from the system's measured behavior. It is well known that a small perturbation of signal (i.e. noise) can substantially change the solution. Ordinary methods for solving algebraic equations (i.e., inversion) do not succeed in solving ill-posed problems therefore special methods must be used. A good survey of numerical methods applicable to general classes of algebraic equations is provided by Björck,25 while reconstruction algorithms suitable for the ill-posed problems are described by Hansen26 In this chapter, we analyze the PPTR inverse problem using the singular value decomposition (SVD). Then, we describe the development original reconstruction codes dedicated to solving the PPTR inverse problem. Finally, we show that our reconstruction codes perform significantly better than the general purpose algorithms available on the market and earlier dedicated codes for PPTR profilometry. 4.1 Discrete ill-posed problems Basic relations of PPTR temperature profiling were derived in Chapter 2. The expression for PPTR signal (2.7) represents a Fredholm integral equation of the first kind, which results in a discrete forward problem (2.11). The calculation of kernel matrix components Ki,j (see Eq. 2.10) is a delicate task, since it involves mathematical operations on very large (i.e., exponential function) and very small numbers (i.e., error function). Hence, one must consider the order of mathematical operations and include adequate approximations to avoid unwanted numerical errors (see the Appendix B). Equation (2.11) represents the PPTR inverse problem for T, which is commonly solved by iterative minimization of the residual norm || S - K T ||2, yielding the best approximate solution T. A statistical approach to the minimization problem can be found in Calvetti.28 The process of finding the approximate solution T is also referred to as reconstruction. 4.1.1 Singular Value Decomposition Singular value decomposition (SVD) is a useful tool for analysis of discrete ill-posed problems.2526 The SVD of the matrix K is a unique decomposition of the form n K = ULV = ^ui 1, and the convergence is fastest if ||K||2 is slightly smaller than 1. Therefore we assured the convergence by substituting K and S with a modified kernel matrix and signal vector, i.e. A and b (see 1st and 2nd line of Algorithm 4.1). ), where e is a small real constant (e.g., e = 10-3). The u–method features the semiconvergence, where the regularization parameter is number of iterations /. Algorithm 4.1: The v-method. Vector Treg minimizes the expression ||S - K Tt||2. Vector S is provided as input. A = K / (||K||2 +e); b = S / (||K||2 +e); i = 0; T (0) = 0; r (0) = b; repeat: Hi = 1 + (/' - 1) (2/ - 3) (2/ + 2v - 1) / [(/ + 2v -1) (2/ + 4w -1) (2/ + 2v -3)]; mi = 4 (2/ - 2w - 1) (/ -1> - 1) / [(/ - 2w - 1) (2/ + 4w -1)]; T ( ) =ni T (!-;) + (1 -fii) T ( - 2) + (Oi AT r ( - 1); r (,) = b - A T (!; i = i + 1; until(stopping criterion) Treg = T (i); 4.2.3 Conjugate Gradients (CG) The basic CG iteration scheme for non-selfadjoint matrices K is well known,32 but several implementation variants exist. A comprehensive survey of different CG variants was provided by Hanke,33 while a simple and understandable text about CG method was published by Shewchuk.30 Based on preliminary testing, we have selected the CG least-square (CGLS) algorithm presented in Algorithm 4.2.26 An essential property of the CG iterates T (i) with residual vectors r (i) = S – K T (i) is that the corresponding vectors KT r (i) are orthogonal. An important consequence is that if the starting vector T (0) is zero, then the solution norm ||T (i)||2 increases monotonically with i and the residual norm ||r(i)||2 decreases monotonically with i.33 The CG method produces iteration vectors in which the spectral components associated with the large singular values converge faster than the remaining components, hence the CG features an inherent semiconvergence. 4.3 Regularization 29 Algorithm 4.2: The CG least-square (CGLS) algorithm. Vector Treg minimizes expression ||S – K Ti||2. Vector S is provided as input. i = 0; T (0) = 0; r (0) = S – K T (0); d (0) = KT r (0); repeat: ?i = ||KT r (i – 1)||22 / ||K d (i – 1)||22; T (i) = T (i – 1) + ?i KTr (i – 1); r (i) = r (i – 1) – ?i K d (i – 1); ßi = ||KT r (i)||22 / ||KT r (i – 1)||22; d (i) = KT r (i) + ßi d (i – 1) ; i = i + 1; until(stopping criterion); Treg = T (i); 4.3 Regularization An ill-posed problem can be regularized by adding a regularization term to the problem (e.g. Tikhonov regularization),34 which penalizes large solution estimates. This regularization approach suffers from low computational efficiency, because the reconstruction is repeated several times for different values of regularization parameter. But the early termination approach significantly reduces computational costs since each successive reconstruction algorithm step represents a new regularized solution. Regarding the accuracy of reconstruction results of PPTR inverse problem, both approaches are equal,6 thus we use early termination. The regularization parameter in TSVD is the number p of right singular vectors vi included in the regularized solution T (p), while the iteration number i serves as the regularization parameter in CG and the ?-method. Since the experiments can often not be repeated, we would like to extract as much information as possible from the given data. If we choose a regularization parameter smaller than the optimal parameter, we leave out too much information and the regularization error dominates. In this case the solution is over-regularized. On the other hand, if we choose a regularization parameter larger than optimal, then the perturbation error dominates and the solution is under-regularized. So far we have presented three different reconstruction algorithms, but no method for selection of the regularization parameter. Although many parameter-selection methods exist,26 we discuss below the discrepancy principle, the generalized cross validation and the L-curve method, which can be successfully applied to our inverse problem. 4.3.1 Discrepancy Principle In practice various types of errors are present in S and K. Common sources of error are measurement noise, approximations, errors due to the discretization process, and also the rounding errors. We denote the vector including all errors as e. If knowledge or an estimate of the norm of signal perturbation ||e||2 is available, the discrepancy principle suggests stopping the iteration when the corresponding residual norm is approximately equal to ||e||2 30 Chapter 4 Reconstruction algorithms S -KT (i ) = Se (4.5) where ?e ? ||e||2. In general, the discrepancy principle tends to produce over-regularized solutions. 4.3.2 Generalized Cross-Validation Generalized cross-validation (GCV) was first presented by Craven35 and Golub.36 GCV is an ||e||2-free method for choosing the regularization parameter based on statistical considerations. The optimal regularization parameter iopt is found at the minimum of the following GCV function G(i) (i) G(i) Ls-Kt [n-/>(i)] (4.6) where ?(i) is the sum of all filter factors in the regularized solution T(i) (see Appendix A). Unfortunately, the minimum can be very flat, which leads to numerical difficulties in computing the minimum of G(i). For TSVD, ?(i) = i and thus the GCV function is simplified to G(i) S-KT (i ) (n-i) (4.7) A different approach to computing the denominator in the GCV function, which is particularly attractive for iterative methods, is to approximate the denominator by a statistical estimate. This approach is called Monte Carlo GCV.37 The idea is to simultaneously run the iterative method on the given left-hand side S and on a random left-hand side S whose elements have a zero mean and standard deviation ?0. If T (i) denotes the additional iteration vector corresponding to S, then er; [n-/>(i)] = ST(S-KT(i)) (4.8) This estimate can be used in the GCV function. However, the Monte Carlo CGV doubles the amount of work in an iterative regularization method. 10-3 °--0-O -O-o --0--0-0 -O -0-O--0 --0--0-0-o 10 Figure 4.3: G(i) as a function of added singular vectors i, when the unconstrained TSVD is applied to a simulated PPTR signal with SNR = 100 (Fig. 4.6). Arrow indicates the minimum. Figure 4.3 shows the GCV function G(i) (4.7) when the unconstrained TSVD is applied to a simulated PPTR signal with SNR = 100 (see Fig. 4.6). The minimum of the curve is indicated by the arrow (i = 8). 2 2 2 2 2 -2 10 -4 10 0 5 15 20 4.4 Non-negativity constraint 31 4.3.3 L-Curve The L-curve criterion was first introduced by Lawson.38 It is based on a parametric plot of logarithm of the norm ||T(i)||2 versus logarithm of the corresponding residual norm ||S – K T(i)||2, with the iteration count i as the parameter. The L-curve theoretically consists of an almost horizontal and a steeply ascending part. The horizontal part corresponds to over-regularized solutions and the ascending part corresponds to under-regularized solutions. The L-shaped corner of the L-curve appears for regularization parameter close to the optimal value, which balances the regularization and perturbation errors in the regularized solution T(see Hansen,26 chapters 4 and 7). Figure 4.4 shows the L-curve for the same example as in Fig. 4.3. It features a characteristic L-shape with a distinct corner at i = 8. y\- 20 8 t) 6 / / - 8 / = 1 «t>.............................o..........o-............................y 0.0 0.5 1.0 log(||S - KT.||2) Figure 4.4: L-curve for the unconstrained TSVD applied to a PPTR signal with SNR = 100 (Fig. 4.6). Arrow indicates the corner of L-curve (i = 8). 4.3.3 Regularization and reconstruction algorithms: summary To find which combinations of regularization selection method and reconstruction algorithm are optimal, we performed an extensive numerical simulation. We have combined the unconstrained TSVD, ?-method and CGLS with the above regularization parameter selection methods, and applied these combined algorithms to PPTR signals for different test objects presented in the last section of this chapter. Since the study was extensive, and the study details are not essential for this thesis, we only summarize the results. Optimal results were obtained when the ?-method was combined with the Monte Carlo GCV, and CGLS was combined with the L-curve method, while TSVD yielded comparable results when combined with either GCV or L-curve. In contrast, the discrepancy principle was only useful if accurate noise information was available. We have found that GCV and L-curve occasionally yield under-regularized solutions, which was also reported by Hansen.26 Hence, we have included in our reconstruction codes also the discrepancy principle to improve their robustness. 4.4 Non-negativity constraint One of the shortcomings of the discussed reconstruction algorithms is presence of unrealistic negative temperatures in the solutions. Because we know that the temperature changes in PPTR profilometry are strictly non-negative, it is very desirable that the reconstructed temperature profile T has only nonnegative values. Thus we search for solution of a constrained minimization problem min S - KT 2 , subjec to Ti ? 0, foralli (4.9) 8 6 4 2 32 Chapter 4 Reconstruction algorithms Milner et al6 reported a successful implementation of the non-negativity constraint into the CG method. They observed significantly improved accuracy and increased robustness of reconstruction results, when the non-negatively constrained reconstruction algorithm was applied to the PPTR inverse problem. Unfortunately, the authors did not reveal the implementation details. Different techniques for computing regularized non-negative solutions (images) to linear discrete ill-posed problems have been recently proposed in the literature. The projected Landweber method was discussed by Bertero and Boccacci.39 Hanke et al suggested the use of a nonlinear transformation of the unknown variable, which eliminates the need to explicitly impose the constraint.40 Nagy et al suggested a modified steepest descent method, referred to as MRNSD.41 Calvetti et al proposed two iterative schemes, based on the generalized minimum residual method (GMRES) or the CGLS method.42 Verkruysse et al demonstrated implementation of the non-negativity constraint into TSVD.46 The efficiency of various iterative reconstruction algorithms was studied by Favati et al.43 In the following we present two approaches for implementation of the non-negativity constraint into the presented reconstruction algorithms. In addition to these, we have also tested three algorithms that do not require any special non-negativity strategy. The expectation-maximization algorithm (EM)44,41 and the image space reconstruction algorithm (ISRA)45, both popular in astronomy and medical imaging, have yielded poor reconstruction results. The MRNSD41 algorithm resulted in somewhat better solutions, but still significantly less accurate as compared to the reconstruction codes discussed below. 4.4.1 Projection When projection is implemented into the reconstruction algorithm, all the negative components of iterative solutions T(i) are set to zero in each step of the algorithm (i.e., T(i) is cropped). This strategy can be successfully implemented into the ?-method, where semiconvergence is preserved.26 The non-negatively constrained reconstruction algorithm based on the ?-method equals Algorithm 4.1, except that T (i) is cropped before the new residual vector r (i) is calculated (line 10). In contrast, a direct implementation of non-negativity constraint in CGLS conflicts with underlying assumptions of the method and severely degrades its convergence and stability. 4.4.2 Projected-Restarted approach In the projected-restarted (PR) approach, the unconstrained reconstruction algorithm runs in the inner loop, while the non-negativity constraint is applied in the outer iteration loop.42 We now describe two new algorithms, where the PR approach is applied to the regularized CGLS and TSVD. We begin with an unbiased starting approximation (i.e., T (0) = 0), which makes the initial residual vector (r (0)) equal to the signal vector S. One full run of the inner loop (one of the unconstrained regularization methods) is then performed, yielding a temporary unconstrained solution t. This vector is truncated to become the first non-negative approximate solution of the outer loop (T (1)). Only the new residual vector (r (1) = S – KT (1)) is then passed to the inner loop for another run. In this way, each successive run of the inner loop result t provides an unconstrained correction to the current outer solution. After the sum of the two is cropped (i.e., (T (1) + t) ? 0), we obtain an improved non-negative solution (T(2)). This sequence is repeated until a predefined convergence criterion is reached, or until the number of outer iteration steps exceeds a preset maximal number (imax). The last T (i) represents the regularized solution of the minimization problem (2.11). Algorithm 4.3 presents a practical realization of the PR approach for the regularized CGLS and TSVD reconstruction algorithms. 4.5 Performance of constrained reconstruction algorithms 33 Algorithm 4.3: Projected-restarted (PR) approach combined with the CGLS or TSVD reconstruction algorithm. The non-negativity constraint is implemented in the outer loop, while the inner loop performs the unconstrained regularization using either unconstrained regularization algorithm. (...)+ indicates the cropping operator. T (0) = 0; r (0) = S; i = 0; repeat: t = CGLS(r (i)) or TSVD(r (i)); T (i + 1) = (T (i) + t)+; r (i + 1) = S – K T (i + 1); i = i + 1; until (||T (i + 1) – T (i)||/|| T (i)|| < ? or i < imax); Treg = T (i); 4.5 Performance of constrained reconstruction algorithms In this section we present reconstruction results of different test objects obtained by our non-negatively constrained reconstruction algorithms with automatic regularization and results of other reported PPTR studies. 4.5.1 Hyper-Gaussian test profile The performance of the dedicated CGLS reconstruction code is demonstrated here by re-evaluation of a test example from an earlier PPTR study,18 where a general-purpose commercial optimization software was used (Solver, part of Microsoft ExcelTM). In order to provide an objective evaluation of the algorithm performance, we present the basic outline and simulation parameters from the earlier study. The initial temperature profiles have a hyper-Gaussian form: AT(z, 0) = AT0 exp[-2(z-z0)4/w4] with AT0 = 10 K, w = 0.1 mm and z0 = 0.3 mm. PPTR signal is computed according to (2.11) and is represented as vector S with 250 equidistant components (At = 2 ms), augmented by normally distributed white noise (SNR = 300). Mismatch between the image and test object is evaluated by calculating normalized quadratic norm of the difference between the image and object vectors {iT—T II represents a time average of either PPTR signal ?S(t) or noise e2. From each signal vector S, the temperature profile T is reconstructed using the kernel matrix K and three reconstruction algorithms. We applied the PR-TSVD algorithm with the GCV regularization parameter selection method, the PR-CGLS algorithm with L-curve, and the ?-method with the Monte Carlo GCV, since numerical simulation results show that these combinations are optimal (section 4.3). The maximum number of iterations is imax = 20,000 and the convergence criterion is set to ? = 10-5. a 20 10 0.0 0.2 0.4 0.6 z (mm) 1.0 0.0 0.5 1.0 t (s) 1.5 2.0 Figure 4.6: (a) Temperature profiles T0 and (b) corresponding theoretical PPTR signals S0 for the test object featuring two absorbing layers Figure 4.7 presents a statistical analysis of 20 regularized solutions. Solid black lines connect the average solution values, while the gray bars indicate the standard deviation. At a given SNR, all three reconstruction algorithms result in comparable solutions. Standard deviation of the solution T decreases with increasing SNR for the PR-TSVD and PR-CGLS, but increases for the projected ?-method – likely due to more right singular vectors vi included in the solution T. SNR = 100 SNR = 500 SNR = 1000 40 20 0.0 0.5 0.0 0.5 0.0 0.5 1.0 z (mm) Figure 4.7: Statistical analysis of 20 regularized solutions T for the second test object obtained using the PR TSVD (top row), PR CGLS (middle row), and projected ?-method (bottom row). SNR levels are 100 (left), 500 (center), and 1000 (right). The actual object is presented for comparison (dark-gray line). Gray bars represent standard deviations. The corresponding relative image errors are presented in Fig. 4.8. The ?-method results in the smallest image error (dark-gray bars; ? = 0.27, 0.26, and 0.24 at SNR = 100, 500, and 1000, respectively). With this method, we also obtain minimal standard deviations of ? and minimal ?max, which is the maximal obtained image error value for the same example and SNR. However, the differences between ? obtained by different methods are insignificant at SNR = 500 and 1000. 1/2 0 60 40 20 0 0 40 20 0 36 Chapter 4 Reconstruction algorithms 0.4 0.3 0.2 0.1 0.0 100 500 SNR 1000 Figure 4.8: Relative image errors ? (Eq. 4.10) of images of the second test object obtained using the PR TSVD (light-gray bar), PR CGLS (gray bars), and ?-method (dark-gray bars). Corresponding maximal obtained image error ? is represented by dashed-line bar. The temperature profiles reconstructed by our algorithms resemble the actual test object better than temperature profiles presented by Milner et al,6 especially the superficial temperature peak. Furthermore, the average image error error for the ?-method and SNR = 100 (? = 0.27) is markedly smaller than the relative image error (? = 0.35) reported in the previous study, albeit the previous study was performed on the more accurate computer platform (64-bit workstation). Therefore we conclude that our dedicated algorithms perform better than the algorithm used in their study. 4.5.2 Simulated port-wine stain lesion The third test object is the temperature profile obtained from a Monte Carlo simulation of optical transport in digitalized histology of port wine stain lesion.12 This should provide a more realistic estimation of the reconstruction algorithm performance for real port-wine stain PPTR signals. The temperature profile is represented as a vector T0 with 750 equidistant components (Az = 2 µm) over a maximum depth of 1.5 mm (Fig. 4.9a). There are 2000 equidistant components (At = 1 ms) in PPTR signal vector S0 (Fig. 4.9b) over a 2 s time interval. The assumed tissue constants are D = 0.11 mm2 s-1, fi = 26 mm1, and h = 0.02 mm1. 10 5 0.0 0.2 0.4 0.6 0. z (mm) 1.0 t (s ) Figure 4.9: (a) Temperature profile T0 and (b) corresponding theoretical PPTR signal S0 for test object MCPWS. Figure 4.10 shows a statistical analysis of 20 reconstruction results for the MCPWS test object and SNR = 100, 500 and 1000. Similarly to the above example, the reconstruction robustness increases with SNR, most markedly for the ?-method (Fig. 4.10, bottom row). At SNR = 100, PRCGLS and PRTSVD produce solutions with more detail than the average solution obtained by the ?-method – likely because of stronger regularization. At SNR = 500 and SNR = 1000, all reconstruction algorithms yield very similar results. As seen in Figure 4.11, the PR CGLS algorithm results in minimal image errors (? = 0.27 ± 0.03 and 0.19 ± 0.02) at SNR = 100 and 1000, while the ?-method is optimal at SNR = 500 (? = 0.22 ± 0.01). CGLS and the ?-method also present acceptably small ?max (dashed line). In contrast, the PR TSVD produces reconstruction results with largest standard deviation of ? (?? = 0.06), and largest ?max. 4.6 Conclusions 37 SNR = 100 SNR = 500 SNR = 1000 0.5 z (mm) Figure 4.10: Statistical analysis of 20 regularized solutions T for MCPWS test object obtained using the PR TSVD (top row), PR CGLS (middle row), and projected ?-method (bottom row). SNR levels are 100 (rigt), 500 (center), and 1000 (left). The exact solution is presented for comparison (dark-gray line). Light-gray bars represent standard deviations. PR TSVD 0.5 ! |____ PR CGLS ! j ¦ projected u 0.4 0.3 I 71 I _L 0.1 100 500 SNR 1000 Figure 4.11: Relative image errors ? (Eq. 4.12) of images of MCPWS test object obtained using the PR TSVD (light-gray bar), PR CGLS (gray bars), and ?-method (dark-gray bars). Corresponding maximal obtained image error ? is represented by dashed-line bar. 4.6 Conclusions Our dedicated reconstruction algorithms with implemented non-negativity constraint and regularization parameter selection method produce markedly more accurate solutions of the PPTR inverse problem than the general-purpose optimization software. Furthermore, our algorithms result in temperature profiles with smaller image errors than temperature profiles presented in an earlier PPTR study and obtained also by a dedicated non-negatively constrained reconstruction algorithm. All three presented algorithms generate similar reconstruction results, but the PR CGLS and projected ?-method yield somewhat smaller image errors and more stable reconstructions as compared to PR TSVD. 10 5 5 5 L' Chapter 5 Calibration of radiometric signal In Chapter 2 we derived the expression for radiometric signal s(t) corresponding to the sample temperature Tb + AT(z, t) (Eq. 2.5). The absolute temperature rise AT(z, t) can not be determined from s(t), unless the exact value of constant C is known. Moreover, (2.5) must be linearized (Eq. 2.6) to obtain the PPTR inverse problem (2.11) which can be solved using the reconstruction algorithms described in Chapter 4. Yet, radiometric signal s(t) does not depend linearly on AT(z, t). In order to overcome these limitations, it is customary to transform radiometric signals s(t) into radiometric temperature using a calibration process. The calibration eliminates the constant C and reduces the nonlinear dependence of AT(z, t) on s(t), but some error due to linearization of (2.6) may still exist. Therefore, in all reported PPTR studies authors justied utilization of (2.11) only for relatively small temperature rises.6918 Because moderate temperature changes as compared to Tb can also be observed in temperature profiles reconstructed from PPTR signals measured on PWS birthmarks, we analyze the effect of linearization on the accuracy of reconstructed temperature profiles. 5.1 Linearization error in monochromatic approximation The calibration process consists of fitting measured PPTR signals to a PPTR setup response to IR radiation emitted by a blackbody set to different temperatures TBB. Here we assume that an IR detector detects emitted IR radiation in a narrow spectral band, thus constant values for wavelength A, responsivity R(X), IR absorption coefficient µ, and constant C can be used. Such monochromatic approximation is used in almost all reported PPTR studies. From (2.5), the expression for radiometric signal sBB(TBB) due to blackbody radiation at TBB is sBB(TBB)=C R(A) BÄ(TBB)AÄ + D (5.1) where D accounts for dark current. 0.005 0.004 0.003 0.002 0.001 0.000 280 300 320 340 360 TBB (K) Figure 5.1: Radiometric signal sBB(TBB) as a function of blackbody temperature TBB for three acquisition wavelengths ? (see labels). The detector responsivity is modeled as R(?) = Rp ?/?p with ?p = 5.0 µm and D = 0. Figure 5.1 presents blackbody radiometric signals sBB(TBB) as a function of TBB for three wavelengths ? = 3.0, 4.0, and 5.0 µm. The detector responsivity R(?) is modeled as R(?) = Rp ?/?p with ?p = 5.0 µm. 39 40 Chapter 5 Calibration Linearization of (5.1) at a fixed TBB results in an error, which increases with increasing temperature difference AT. So, when temperature is increased by AT = 20 K, the expression linearized at TBB = 300 K yields 13% smaller values for X = 5.0 µm, and even 24% smaller values for X = 3.0 µm (see Fig. 5.1) as compared to the actual value of SBB. The expression (5.1) is fitted to the measured PPTR setup response to the IR radiation emited by the blackbody, and the constants C and D are determined, thus yielding a calibration curve sBB(T). When calibrating radiometric signal s(t), the corresponding radiometric temperature S(t ) is then found by matching each s(ti) to the calibration curve sBB(T) sBB(S(ti)) = s(ti) (5.2) Thus the calibration procedure eliminates the constants C and D, and establishes a relation between s(t) and radiometric temperature S(t). An expression for Planck’s law of radiation is 2nhc2 B* ( T ) = r----------------------1 (5.3a) A5 (exp(hc/AkBT)-) which can be approximated with hc Bji(T)=5 c 2 e W (5.3b) A since the exponential in B?(T) is much larger than 1. Using (5.3b) we find an expression for a monochromatic PPTR signal 2nCAAjLiR(A)hc2 hc - /Hz s( t ) = e AkB ( Tb+aT ( z , t )) dz + D (5 4a) I5 z=0 and an expression for the calibration curve 2 hc sBB(T) = 27rCAÄR 5 hc e *kbt + D (5 4a) A We substitute (5.4a) and (5.4b) into (5.2) to get the expression for calibrated PPTR signal s(t) in radiometric temperature units AS(t) =-----------------------hc---------------------T (5.5) AkBlnj Me "B(T b4T z , t)e z d z The calibration process eliminates the constant factor 2?C??hc2R(?) and constant D. We expand the expression for ?S(t) into Taylor series and neglect terms of the order ?T(z, t)2 and higher, which results in a linearized expression for the PPTR signal ASL (t) = ju\ &T(z,t) e~Mz dz (5.6) All factors involving k and Tb are canceled in (5.6), so that the linearized calibrated signal ASL(t) depends on acquisition wavelength X only indirectly, through the value µ. z=0 z=0 We can estimate the linearization error by considering the quadratic term in Taylor series 5.1 Linearization error in monochromatic approximation 41 N2(xT(z, t))=1 hc 2XkT V CO j *T(z,t) b J 2 e^z dz y z=0 CO j AT(z,t ) e^zdz (5.7) Assuming that N2(AT(z, t )) (5.7) is the largest contribution to the linearization error, we find an expression for the relative linearization error. First, we write the temperature rise as AT(z, t) = Tav f(z, t), where Tav is average temperature rise and f(z, t) is a normalized form function in the sense CO fi J f(z,t) e-"z dz = 1. Second, we divide (5.7) by (5.6) to obtain N2 (AT ( z,t )) = Tav ASL(t ) Tb hc 2AkBTb M CO Jf2( z , t )e dz -1 (5.8) Evidently, the linearization error increases proportionally with Tav. Furthermore, the error increases with decreasing Tb and L The above expression shows that for a sample with constant temperature T (i.e., f(z, t) = 1) the linearization error is zero; which indicates that the linearized expression ASL(t) is exact for homogenously heated samples. However, for non-homogenously heated samples (i.e., f(z, t) + 0) the linearization error persist in spite of optimal calibration. For example, we evaluate (5.8) for a Gaussian initial temperature profile AT(z) = T0 exp(-(z - z0)2/2w2) with T0 = 20 K, z0 = 100 µm, and w = 30 µm for X = 3 and 5 µm. Absorption coefficient µ = 26.5 µm. Thus, we obtain the linearization error N2/ASL = 20% and 11% for X = 3 and 5 µm, respectively. These errors are smaller than the corresponding errors for non-calibrated linearized radiometric signal (see Fig. 5.1). Since PPTR signal involves transient temperature rise AT(z, t), we perform a numerical simulation to consider evolution of the linearization error with time. 5.1.1 Numerical simulation The initial temperature profiles have a hyper-Gaussian form: ?T(z,0) = Tp exp[-2(z–z0)4/w4] with Tp = 1 – 100 K, w = 30 µm, and central depths z0 = 100, 200, and 300 µm. The temperature profile ?T(z’, t) at depth z’ and time t is calculated as a convolution of ?T(z, 0) and the Green’s function GT(z, z’, t) (Eq. 2.3). Then, ?S(t) and ?SL(t) are calculated by inserting ?T(z, t) into (5.5) and (5.6) and completing the z integrals. We use numerical integration routines implemented in Mathematica 4.1 (Wolfram research, Champaign, IL) to compute signal vectors with 1000 components representing signal values acquired at a sampling rate 1000 s-1 for a total acquisition time of 1s. PPTR signals are simulated for ? = 4.5 µm, µ = 26.5 mm-1, and Tb = 303 K. Thermal constants in GT(z, z’, t) are D = 0.11 mm2/s and h = 0.02 mm-1. Figure 5.2 presents PPTR signals calculated using the exact (solid line) and the linearized (dashed line) expressions for test objects with z0 = 100 µm and Tp = 1, 10, and 100 K (see the labels). For Tp = 100 K, the exact signal features an initial spike (see the arrow), which is not present in the linearized signal. This difference between the simulated signals appears only at small times (t < 20 ms); while both signals are equal at later times. For the lower Tp (1 and 10 K), both signals are equal at al t. z=0 42 Chapter 5 Calibration 20 16 12 8 2 0.0 0.2 0.4 0.6 0.8 1.0 t (s) Figure 5.2: Theoretical calibrated PPTR signals for hyper-Gaussian test objects with z0 = 100 µm, and Tp = 1, 10, and 100 K (see labels). The exact ?S(t) (Eq. 5.5; solid line) and linearized PPTR signals ?SL(t) (Eq. 5.6; dashed line) are presented. In order to find how the linearization error affects temperature profiles, we simulate the exact PPTR signals (Eq. 5.5). Then we reconstruct initial temperature profiles ?T(z, 0) from these signals using the linearized expression for PPTR signal (Eq. 2.11) and the monochromatic approximation (Eq. 2.9). The reconstruction algorithm used is the ?-method (chapter 4). We simulate PPTR signals for test objects located at z0 = 100 µm with Tp = 1–100 K, and test objects with Tp = 50 K located at z0 = 100-300 µm. Figure 5.3 presents the reconstructed temperature profiles (solid line) and the actual test objects (dashed line) for comparison. 0.0 0.8- A T p = 1 k p a 0.4- ] 0.0- , i T = 10 k 8- 4- i\ 0: 80- / / ¦'P T p = 100 k 40- t 0.1 0.2 0.3 0.4 0.5 A b 40 i \ 20 V ;1 z0 = 100 um 0 40 A 20 / \ / I \ z0 = 200 um 0 40 A 20 i '\ z0 = 300 um 0.0 0.1 z (mm) 0.2 0.3 z (mm) 0.4 0.5 Figure 5.3: Reconstructed temperature profiles ?T obtained from the exact PPTR signals ?S(t) for (a) test objects located at z0 = 100 µm with Tp = 1–100 K (see labels), and (b) test objects located at z0 = 100–300 µm with Tp = 50 K (see labels). The actual objects (dashed line) are plotted for comparison. Arrows indicate surface artifacts due to the linearization error. When reconstructed test object is located at z0 = 100 µm (Fig. 5.3a), the reconstruction results are not significantly deteriorated for Tp = 1 and 10 K, but we can observe a small artifact near the object surface (see the arrows). In contrast, reconstruction of test object with Tp = 100 K (bottom) results in a severely deformed temperature profile, and a significant surface artifact. Moreover, the temperature profile is shifted deeper as compared to original depth (z0 = 100 µm). 0 5.2 Broad-band signal acquisition 43 The effect of linearization error is reduced as z0 increases (Fig. 5.3b). The distinct surface artifact in reconstructed temperature profile of test object with z0 = 100 µm (top) and Tp = 50 K, diminishes as z0 is increased to 200 µm (center), and disappears for z0 = 300 µm (bottom). The shift of temperature peak (~3 µm) is also present only for the temperature profile located at z0 = 100 µm. Occasionally, we observe such initial spikes in measured PPTR signals. Because the temperature profile reconstruction is performed using the linearized expression (5.6), superficial artifacts due to the linearization error can appear in reconstructed temperature profiles, especially when large temperature amplitudes are involved. In addition to the above reconstruction results, we present a systematic analysis of the linearization error. Figure 5.4 shows the normalized linearization error defined as (?S(t) – ?SL(t))/ for all simulated hyper-Gaussian test objects. In general, the normalized linearization error increases with increasing Tp. However, the error practically vanishes in 10–100 ms for all depths z0 and temperatures Tp. Specifically, for Tp = 100 K the linearization error is reduced below 1% in 20 ms and 46 ms for z0 = 100 and 200 µm, respectively. The linearization errors for the test objects with z0 = 300 µm feature peak values at t = 40–50 ms, but are significantly smaller than for z0 = 100 and 200 µm. 0.04 0.03 0.02 0.01 0.00 3.0 jim 3.5 jim 4.0 jim 4.5 jim 4.9 jim 0.000 0.005 0.010 t (s) 0.015 0.020 0.00 0.00 0.01 0.02 0.03 t (s) 0.00 0.01 0.02 0.03 Figure 5.4: The normalized linearization error for all simulated hyper-Gaussian test objects (see labels) calculated as (AS(t) - ASL(t))/, where is the average value of AS(f). The peak temperatures are Tp = 1 - 100 K (see legend), the wavelengths = 4.5 µm and the IR absorption coefficient ii = 26.5 mm-1. t (s) 5.2 Broad-band signal acquisition Since the monochromatic calibration is accurate only for very narrow spectral bands,18,47 we analyze next the calibration process involving broadband signal acquisition. From (2.5), radiometric signal due to blackbody radiation at TBB is (7bb) = c J R(fy BÄ(TBB) d/t + D (5.9) Tp = 1 K Tp = 10 K Tp = 100 K s 44 Chapter 5 Calibration where Xl and Xh are the lower and the upper limit wavelengths, respectively. Analogous to the monochromatic calibration, the expression (5.9) is fitted to the measured broad-band PPTR setup response to the IR radiation emited by the blackbody, and the constants C and D are determined, thus yielding a broad-band calibration curve sBB(T). The broad-band radiometric temperature S(t ) is then found by matching each s(ti) to the calibration curve sBB(T) according to (5.2). By using the Taylor series expansion, the linearized expression for the calibrated spectrally composite PPTR signal is \ oo I \ ASL(t) = h R(X) Bl (Tb) KX) J AT(z,t) e~"wz dzdÄ / h R(X) BÄ' (Tb) dÄ (5.10) With broad-band signal acquisition, the factors involving specific spectral properties (i.e., R(X), Bk’(Tb), µ(X)) and Tb do not cancel out in the linearized calibrated PPTR signal ASL(t) (5.12). Because we can not derive the linearization in the spectrally composite PPTR signals AS(t) analytically, we perform a numerical simulation. T = 1 K 0.04 q p 3.0 )jm ....... 3.5 )jm i ....... 4.0 )jm 0.03 , 4.5 )jm ------- 4.9 )jm 0.02 I 0.01 V^iw—— 0.000 0.005 0.010 0.015 0.020 t (s) t (s) Tp = 100 K 0.00 0.01 0.02 0.03 t (s) Figure 5.5: The normalized PPTR signal linearization error for all hyper-Gaussian test objects with peak temperatures are Tp = 1-100 K (see the labels) and the spectral bands with^l = 3.0–4.9 µm (see the legend). Figure 5.5 present the obtained relative linearization errors calculated as in subsection 5.1.1. Similarly to the monochromatic case, the linearization error practically vanishes in 20-30 ms, and increases with the temperature amplitude Tp for all spectral bands. In addition, the linearization error progressively decreases, when the spectral acquisition band is reduced from k = 3.0-5.0 µm to k = 4.9-5.0 µm. The initial temperature profiles AT(z, 0) have a hyper-Gaussian form, as in subsection 5.1.1. The profile central depth is fixed at z0 =100 µm, where the monochromatic linearization error was found to be largest. Detector spectral responsivity is modeled as R(X) = Rp l/lp, and the IR absorption coefficient µ^) is modeled as 75% of the values for water48 sampled with spectral interval of 0.1 µm. The exact spectrally composite PPTR signals AS(t) are computed in three steps. First, the radiometric signal s(t) is found by completion of the z and X integrals (Eq. 2.5). Second, the blackbody radiometric signal AsBB(TBB) is determined for TBB = 273-373 K by completion of the X integral (Eq. 5.9). Finally, AS(t ) is calculated by fitting s(ti) to the calibration curve sBB(T) (Eq. 5.2) at each time t i = i At. Tp = 10 K 5.3 Conclusions 45 Meanwhile, the corresponding linearized spectrally composite PPTR signal ?SL(t) is computed by completion of the z and ? integrals in (5.10). All computations were performed using Mathematica 4.1. Table 5.1 lists the maximum linearization error and the time t1%, when the linearization error drops below 1% of the maximum value for the entire-spectral band of the InSb detector (? = 3.0–5.0 µm) and a reduced spectral band (? = 4.5–5.0 µm).18,49 Maximal linearization error for ?l = 3.0 µm is on average 50% larger than the corresponding error for ?l = 4.5 µm. t1% increases with Tp, and is smaller for ?l = 4.5 µm than for ?l = 3.0 µm. Table 5.1: Maximal normalized linearization error (?S(t)– ?SL(t))/< ?S(t)> and time t1%, when the linearization error drops below 1% of the maximum. Results are presented for spectral bands with ?l = 3.0 µm and 4.5 µm, and for the monochromatic PPTR signals with ? = 4.5 µm and z0 = 100 µm (Fig. 5.5). X = 3.0-5.0 µm X = 4.5-5.0 µm T0 (K) Error t1% (ms) Error t1% (ms) 1 0.04 1 0.01 0 2 0.05 2 0.03 1 5 0.06 3 0.03 2 10 0.11 8 0.07 6 20 0.15 13 0.13 9 50 0.52 16 0.35 14 100 1.10 20 0.75 16 5.3 Conclusions All nonlinearity between the radiometric signal and the radiometric temperature is eliminated, when a calibration in monochromatic approximation is performed. In linearized calibrated PPTR signals the factors involving the acquisition wavelength (i.e., B?(Tb) and R(?)) cancel out from the resulting expression. The linearized expression is exact for samples with homogenous temperature distribution, but for samples with non-homogenous temperature distribution the linearization produces a linearization error. In general, the linearization error is significant for absorbing structures located at shallow depths and short after laser pulse. The error increases with increasing peak temperature, which is predicted by the analytical expression (5.8). This expression also indicates that linearization error decreases with increasing acquisition wavelength ? and baseline temperature Tb. Due to the linearization error we can observe in reconstructed temperature profiles superficial artifacts, which can be easily mistaken for superficial absorbing structures. Same trends are obtained for the spectrally composite calibrated PPTR signals. In addition, we observe also that the linearization error increases with spectral acquisition band broadening. Since temperature amplitudes are in general moderate, and most of absorbing structures (i.e., blood vessels) are deeper in the sample, the linearization error in typical PPTR temperature profiling in human skin is small. The error can be additionally reduced by appropriate spectral filtering. Chapter 6 Effective infrared absorption coefficient Reconstruction of temperature profiles from PPTR signals involves tissue absorption coefficient at the detected IR wavelength, µ(X). In general, soft biological tissues feature a pronounced spectral variation µ(X) in mid-IR. Yet, all reported PPTR studies utilized a fixed effective value (µeff) to reduce computational complexity of the reconstruction process, albeit broad-band signal acquisition was used almost invariably to increase signal-to-noise ratio (SNR). In earlier PPTR studies two approaches of analytical estimation of µeff were presented,69 but no study has been done to find which approach is better. However, implications of monochromatic approximation for PPTR temperature profiling were analyzed numerically for specific case of two spectral bands of an InSb detector.18 In this chapter, we focus on determination of the optimal effective value µeff for samples that exhibit a significant spectral variation µ(X) within the IR detection band. Taking into account spectral characteristics of three common IR radiation detectors, we simulate realistic PPTR signals for four hyper-Gaussian temperature profiles centered at different depths below the sample surface. From each PPTR signal, the initial temperature profile is then reconstructed using the approximation with a constant µeff. Analysis of the mismatch between the actual and reconstructed temperature profiles at different values µeff enables determination of the optimal effective value, µopt. Because determination of µ opt from numerical simulations is very tedious, we propose a novel analytical approach to the same task. Based on an implicit equation derived earlier1118, this approach enables a direct estimation of µopt from spectral properties of the sample and radiation detector. In a systematic analysis, involving multiple acquisition spectral bands for each IR radiation detector, our analytically assessed values µ opt match the numerical results very well, unlike the previously used analytical estimates.69 6.1 InSb detector 6.1.1 Numerical simulation All initial temperature profiles (“objects”) in our simulation examples have a hyper-Gaussian form: AT(z, 0) = AT0 exp[-2 (z-z0)4/w4] with AT0 = 10 K and w = 0.1 mm. With z0 varied from 0.1 to 0.4 mm, these profiles represent laser-induced temperature rise in subsurface vascular plexus (such as PWS) at different depths in skin. The corresponding vectors T0 represent temperature values at 100 equidistant positions over a depth of 1 mm. Each vector S0 has 1000 components, representing signal values acquired at a sampling rate of 1000 s1 (for a total acquisition time of 1 s). It is obtained by multiplying T0 with a spectrally composite matrix K Based on (2.7) and (5.10), each matrix element of the latter is computed as a sum of N = 100 monochromatic kernel functions k(z, t) (Eq. 2.10) evaluated at equidistant wavelengths h within the spectral acquisition band h–K N ^R(Ak)B^(Tb) M(Ak) ^i , tj^(A,)] Az K(zi ,tj) = k 1 N--------------------------------------- (6.1) k=1 47 48 Chapter 6 Effective infrared absorption coefficient where Az denotes depth discretization. The absorption spectrum of human skin is modeled by adding 75% of µ(X) for water48 and 25% of µ(X) for collagen51 (Fig. 6.1a). The thermal constants in v(z,t) are set to52 D = 1.1 x 10-7 m2s-1 and h = 0.02 mm1. InSb quantum radiation detector has peak sensitivity wavelength of ^ = 5.3 µm and cut-off at lc = 5.6 µm. The corresponding spectral responsivity R(A) is modeled as increasing linearly with X up to lp, then linearly decreasing to 10% of the maximal value at lc (Fig. 6.1b). The lower spectral limit (A) is varied between 3.0 um and V while Ah is fixed at the respective lc. This simulates application of different cut-on filters to narrow the detection spectral band1850. No noise is added to the simulated signals, because it would induce inconsistencies in the analysis. From each signal vector, the initial temperature profile is reconstructed using a range of simplified kernel matrices based on (2.9) and employing different effective values µeff. We apply the PR-CG method (chapter 4) to solve the inverse problem (2.11). The minimum in dependence of relative image error S (4.10) on µeff indicates the optimal effective value (µopt) for a given combination of object, detector and acquisition spectral band. 3 4 5 6 7 8 9 10 11 12 13 14 Wavelength [um] Figure 6.1: (a) Model absorption coefficient of skin in mid-IR spectral region, /u (solid line), computed as 75% of µ for water48 (short dash) increased by 25% of µ for collagen51. (b) Temperature derivative of Planck’s radiation formula at Tb = 300 K (solid line) and spectral sensitivity R(X) of a typical InSb detector (dashed line), and two HgCdTe detectors with peak sensitivity at 10 µm and 12 µm (dash-dot and short-dash line, respectively). Figure 6.2 presents three results, obtained from a simulated PPTR signal for the hyper-Gaussian object centered at z0 = 0.3 mm (dashed line) and employing the entire acquisition spectral band of the InSb detector (3.0-5.6 µm). Reconstructions were performed using the simplified, quasi-monochromatic kernel matrices based on (2.9) with varying µeff. At a sub-optimal value (µeff = 20.0 mm1, top graph), the tails of the profile are squeezed and peak temperature is markedly overshot. The best match is obtained with µeff = 22.3 mm-1 (center), while an overestimated value (µeff = 25.0 mm1, bottom) results in a narrowed peak with blurred tails. Figure 6.3 presents the dependence of relative image error S (Eq. 4.10) on µeff, as obtained for four hyper-Gaussian profiles centered at different depths z0 (see labels in the graphs). The optimal µeff for each object is indicated by the minimum of the respective dependence (arrows). The results in Fig. 6.3a represent the case where the entire spectral band of the InSb detector is used. With a narrowed acquisition spectral band (A{= 4.5 µm; Fig. 6.3b), the four optimal values nopt(z0) lie closer together and the respective image errors are significantly smaller. 6.1 InSb detector 49 0.5 eff = 20.0 mm-1 " ; / \ \ 22.3 mm-1 ¦ ¦ / V 25.0 mm-1 " 0.4 0.6 Depth [mm] Figure 6.2: Temperature depth profiles (solid lines) reconstructed from simulated PPTR signal employing the entire spectral band of the InSb detector (3.0–5.6 µm). The values of µeff used in reconstruction were, respectively, too low (top), near-optimal (center), and too high (bottom; values indicated in the graphs). The initial temperature profile is plotted for comparison (dashed line). 0.12 A A (a) 0.10 A aA ^ 0.4 mm A^fp 0.08 eff (z i )e --( zi zi i In computing both sums, only the value µ eff(zi) from the solution branch that is closer to µ opt(0) is considered at each depth zi. Using the obtained value, µopt(1), as improved approximation, the procedure is then repeated in the same manner. The successive approximations µopt(m) converge toward the optimal effective absorption coefficient, µopt. The exponential weight function, exp[- µeff(z) z], was selected because the contribution of emitters from various depths z to the radiometric signal (Eq. 2.6) at any time t is governed by Beer’s law. After figuring out how to deal with the obtained relation (Eq. 6.2) involving the two-valued function µ eff(z), everything fell into place. Figure 6.5 presents solutions µeff(z) of the implicit Eq. (6.2) for four acquisition spectral bands of the InSb detector (see the caption). The variation of µeff with depth (z) is most prominent for the widest acquisition band (Fig. 6.5a), and almost absent in the nearly monochromatic case (Fig. 6.5d). The dashed lines indicate the corresponding optimal values µopt assessed using our analytical approach (Eq. 6.2). For the four presented examples they amount to 23.0 mm1, 23.2 mm1, 24.5 mm1, and 20.7 mm-1, respectively. It may seem odd that the optimal µeff depends also on the object depth, z0. But this is just a reminder that no constant value µeff used in Eq. (2.9) can entirely replicate the effect of the spectrally composite K(z,t) (2.8). 6.1 InSb detector 51 40 30 20 10 40 30 20 10 0 40 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 Depth [mm] Depth [mm] 40 30 20 10 30 20- — 10 0.0 0.1 0.2 Depth [mm] 0.3 0.4 0.0 0.1 0.2 Depth [mm] 0.3 0.4 Figure 6.5: Effective absorption coefficient /ze m z) (Eq. 6.2) computed for InSb detector and detection band of: (a) 3.0-5.6 µm, (b) 3.9-5.6 µm, (c) 4.5-5.6 ur, and (d) 5.1-5.6 µm. Dashed lines indicate the corresponding values /uopt determined using Eq. (6.3). Figure 6.6 presents the analytically assessed optimal values, /opt (Eq. 6.3), as a function of the lower limit wavelength, k (solid line). For all included values 4 the analytical estimate falls within the range of values /opt(z0) obtained for the four test objects in numerical simulations (grey vertical stripes). Difference between the arithmetic mean of the latter (open circles) and/opt (Eq. 6.3) does not exceed 0.6 mm1 for any spectral band under test. For the most important case, employing the entire spectral band of the InSb detector (k = 3.0 |am), the difference is only 0.25 mm1 or 1.1% of the target value. For objective PPTR profiling of arbitrary samples, one would namely wish to obtain equally good performance at all subsurface depths. In reality, however, images of deeper objects are increasingly blurred due to higher susceptibility to experimental noise (chapters 7 and 11). To account for this, we consider also a weighted average which favors shallower test objects (black squares). The weight z0-1 was selected because it approximates scaling of the PPTR signal amplitude from a thin subsurface layer. 9 By considering a weighted average of /opt(z0) at each 4 /opt (Eq. 6.3) falls within 0.3 mm1 off the mark for all detection bands included in the analysis. Two formerly applied analytical estimates perform significantly worse in such a comparison. The weighted average of /u(X) with the weight set to Rß)53,9 amounts to 62.3 mm1 in the most relevant case of k= 3.0 |am, overshooting the average of numerical values /opt(z0) by 174% (Fig. 6.6, short-dash line). The more elaborate expression used by Milner et al6,7 ßeff = j R(Ä) BÄTb) M/0 dÄ JR(Ä)B^(Tb)dÄ (6.4) 0 provides a better match to the numerical data (dash-dot line). Nevertheless, it exceeds the optimal 52 Chapter 6 Effective infrared absorption coefficient value for unfiltered InSb detector by a significant margin of 2.4 mm-1 (or 13%). 22 ^ /^ f Simple mean Weighted average Analytical prediction 3.0 3.5 4.0 A[l-im] 4.5 5.0 Figure 6.6: Analytically assessed values µopt as a function of the lower limit wavelength (t) (solid line). Grey stripes indicate the range of µopt(z0) values obtained for the four test objects in numerical simulations; simple means and weighted averages of these values are marked by open circles and black squares, respectively. Two formerly used analytical estimates are plotted for comparison (dash-dot and short-dash lines). 6.2 HgCdTe detectors We present similar analysis to the above for two HgCdTe detectors, also commonly used for PPTR temperature profiling. First HgCdTe detector has peak sensitivity at ?p = 10 µm, the upper spectral limit (?h) is fixed at the cut-off wavelength of the detector, ?c = 12 µm, while the lower limit is varied between 3.0 µm and ?p. 70 60 50 Simple mean Weighted average Analytical prediction 3456789 10 A[ m] Figure 6.7: Analytically assessed µopt (Eq. 6.2) as a function of Xl (solid line) for the HgCdTe detector with peak sensitivity at 10 um. Grey stripes indicate the range of µopt(z0) values obtained for the four test objects in numerical simulations; simple means and weighted averages of these values are marked by open circles and black squares, respectively. Two formerly used analytical estimates are plotted for comparison (dash-dot and short-dash lines; see text for details). Figure 6.7 presents the analytically assessed optimal absorption coefficient value, µopt (Eq. 6.3), as a function of the lower limit wavelength, k (solid line). For all included spectral bands, the analytical estimate falls within (or just marginally deviates from) the range of values µ opt(z0) obtained for the four test objects in numerical simulations (grey vertical stripes). For the maximal spectral band (k = 3.0 |am), the mismatch between µ opt (Eq. 6.3) and simple mean of the latter (open circles) is 3.5 mm1. If 6.3 Discussion 53 the weighted average is taken as the target value (black squares), the mismatch is only 2.0 mm1 (3.8%). Clearly, our analytical estimate /uopt (Eq. 6.3) provides a much better fit to the numerical data than the earlier used fieff (Eq. 6.4) or convolution of /i(X) with R(X) (dash-dot and short-dash lines, respectively) for all 36 spectral bands included in this part of the study. For the HgCdTe detector with peak sensitivity at lp = 12 µm, our analytical estimate /opt (Eq. 6.3) (Fig. 6.8, solid line) falls within the range of numerically determined values /opt(z0) (grey stripes) for each of the included spectral bands. The match is particularly good if the weighted average of the four /opt(z0) is considered at each Xl (black squares). In contrast, the earlier used analytical estimates fieff (Eq. 6.4) and weighted average of fi(X) with weight R(X) (dash-dot and short-dash lines, respectively) deviate significantly from the numerical results even for moderately wide detection bands (e.g., Xl = 10 |am). With further decrease of 4 the mismatch between fieff (Eq. 6.4) (dash-dot line) and simple mean of the numerical /opt(z0) (open circles) stays large and amounts to 25 mm1 (or 43% of the target value) at Al = 3.0|am. 200 150 100 50 Simple mean Weighted average Analytical prediction 3 4 5 6 7 8 9 10 11 12 ?l [µm] Figure 6.8: Same as Fig. 9.7, but for HgCdTe detector with peak sensitivity at ?p = 12 µm. 6.3 Discussion Figure 6.2 illustrates how a ~10% deviation from the optimal value /I& deteriorates the reconstructed temperature profile (top, bottom). In soft biological tissue, such as human skin with /j{X) in mid-IR varying by two orders of magnitude (Fig. 6.1), determining the optimal value ^ff to such a narrow margin is certainly a nontrivial task. In a separate numerical study,47 white noise in the simulated signal (SNR = 300) suspended the high dependence of relative image error (S) on #eff only within a very narrow interval around ^pt. At deviations larger than 2-3% of the optimum value, the influence on S was clearly present, rather than buried in the effect of noise, as is believed quite often. Being aware of the spectral variation /u{K) in their gel samples, Prahl et al53,9 estimated the effective value fieff by averaging /u(X) of water over the 8-12 |am detection band, using spectral sensitivity R(X) of the HgCdTe detector as a weight. The closest match to such experimental conditions is found in Fig. 6.8 (short dash line, h = 8.0 ^m), which suggests that their estimate of fieff may have been too high. Of course, this comparison does not provide any indication of the possible implications of such deviation for their results. In fact, at such high values of ^ we would expect only a minimal effect on profiling, likely limited to the most superficial 10 |am of the sample. On the other hand, Fig. 6.6 leaves no doubt that applying the same approach to the commonly used InSb detector with Xl = 3.0 |am would be disastrous. Milner et at,7 applied a more elaborate expression ßeff (Eq. 6.4) to HgCdTe detector used at 7-11 |^m, and 10-14 |am, respectively. As indicated by Figs. 6.7 and 6.8 (dash-dot line), the 54 Chapter 6 Effective infrared absorption coefficient improvement over the previous approach is only marginal. This can be attributed to the weak spectral variation of B* (Tb) in this part of the spectrum (Fig. 6.1b). On the other hand, a dramatic improvement is obtained for InSb detector at 3-5 um (see Fig. 6.6, dash-dot vs. short-dash line) where BÄ’(Tb) exhibits a large spectral variation. Nevertheless, the remaining discrepancy between /2eff (Fig. 6.6, dash-dot line) and numerically determined µ opt(z0) at h = 3.0 um (grey stripe) suggests that the value used for temperature profiling in human skin by Milner et al6 may have been overestimated. The possible effect of such a discrepancy can be guessed from Fig. 6.2 and extrapolation of data in Fig. 6.3(a) - while keeping in mind the possible differences between the model functions rfX) and R(X) in the two studies. The values µ opt (Eq. 6.3) obtained with our novel analytical approach in general provide a much better fit to the numerical results than the previously used estimates (Figs. 6.6-6.8; solid, short-dash, and dash-dot lines, respectively). (The results converge towards the right end of Figs. 6.6 and 6.8). Moreover, µopt (Eq. 6.3) fall within the range of values µopt(z0) obtained in the numerical simulations (grey stripes) for most of the 68 detector/spectral band combinations tested in this study. The same is rarely the case for either earlier approach. Spectral variation of /li{X) in water is most pronounced in the lower part of detection band for all three detectors in this study (Fig. 6.1). Narrowing of the acquisition band with a cut-on filter significantly improves the validity of quasi-monochromatic approximation (Eq. 2.9). This is evidenced by the decrease of relative image error (S) for all test objects (Figs. 6.3b, 6.4b). The four optimal values µopt(z0) lie closer together (Figs. 6.3b, 6.4a), which should further improve the accuracy of depth profiling when absorbers are distributed at varying (or multiple) subsurface depths. On the other hand, narrower spectral bands will in general result in smaller PPTR signal amplitude and thus decrease SNR. The optimal detection band for given experimental conditions (i.e., µ (ß), Rß)) can thus be determined only by consideration of the related experimental noise. Such analysis is presented in Chapters 9 and 10. Finally, broad-band IR detection is commonly employed also in other PTR techniques, such as measurements of thermal diffusivity, depth profiling of sample thermal conductivity, or nondestructive evaluation using measurements in frequency domain. To our knowledge, most of these techniques also assume an effective absorption coefficient value µeff in the involved signal analysis. Selection of µeff and its influence on the results deserve some thought when these techniques are applied to biological tissues or other materials with significant variation of µ(X) in the detection band. Because rather general relations describing formation of the radiometric signal were involved in our determination of µ opt (Eq. 6.3), the same or slightly adapted approach might provide viable values might be applicable also to some PTR techniques beyond PPTR temperature profiling. 6.4 Conclusions The presented approach enables analytical determination of viable effective IR absorption coefficients to be used in reconstruction of temperature profiles from PPTR signals at any combination of sample spectral properties, radiation detector, and acquisition spectral band. This is particularly important in samples with large spectral variation µ(X) in mid-IR, such as most biological tissues, where previously used analytical estimates are not sufficiently accurate. Chapter 7 Sampling rate Signal sampling rate is one of the experimental parameters involved in PPTR temperature profiling. Reported PPTR studies in biological tissues67181949 utilized a wide range of sampling rates (f = 125-1700 s1), but the reasons for selecting the actual sampling rates were not discussed. Basic relation, which illustrates how the axial resolution Az depends on sampling time At, is the well known expression for the heat diffusion length Az = J4DAt. Evidently, the axial resolution increases with sampling frequency. However, this relation does not include the effects of µ, h and noise. In this chapter we study the effect of sampling rate on the accuracy of reconstructed PPTR temperature profiles. We simulate PPTR signals with different sampling rates for different temperature profiles, and then reconstruct initial temperature profiles. 7.1 Methods We consider a temperature profile resulting from radiative absorption in a layer of width w located at depth z1 fAT 0exp(-k\z-z 1), z :-- — o— 200 p 0.6 5^ — •— 400 p 0.4 "°--"n ^'%;si Ss==e^; ====:»=E===;e 0.2 """"°-------° 100 1000 f (s-1) 130 b 50 fj.m — n— 100 fj.m --¦— 200 fj.m 110 ¦ "\ "¦¦ \ -- o— 400 fj.m 90 "¦-.V '"¦ 70 ^*"*--.^*-- "^- \ 50 ^ Ä= =:?"*=--------i —- JL_::====G 10000 100 1000 f (s-1) 10000 Figure 7.2: (a) Relative image error ? (Eq. 4.12) and (b) full-width at half-maximum of the temperature profiles W. Results for samplig rates f = 100–10 000 s-1 are presented (see legend). Gray line in (b) indicates the actual width (50 µm). 7.2.2 Signals with noise Figure 7.3 presents simulated PPTR signals with realistic noise for f = 100 (left), 1000 (center) and 10,000 s-1 (right). Four test objects located at z1 = 50–400 µm (see labels) are considered. Evidently, PPTR signals become more noisy as f increases (Fig. 7.4). 3 2 0.0 f = 100 s f = 10,000 s \ 50 ^m / ^\^"*\ 100)JT1 A \ / ^ /^ 0.2 0.4 0.6 0.8 0.0 z (mm) 0.2 0.4 0.6 0.8 1.0 Figure 7.3: Simulated PPTR signals for test objects with z1 = 50, 100, 200 and 400 µm (see labels) and sampling rates of 100 (left), 1000 (center) and 10,000 s-1 (right) augmented by realistic noise. S;-- — ¦— 50 fj.m •-... ---*- "¦-.. — d— 100 fj.m — •— 200 fj.m 1000 0.. '-*.. '" ---o— 400 fj.m "~'°-.. "xa. 100 ^"'O 100 1000 f (s-1 ) 10000 4 0 Figure 7.4: SNR of the simulated PPTR signals as a function of sampling frequency f for test objects with z1 50 (solid squares), 100 (open squares), 200 (solid circles) and 400 µm (open circles). 7.2 Results 59 Figure 7.5 presents statistical analysis of the reconstruction results for three sampling rates. While images for f = 100 s-1 (left column) are similar to images reconstructed from noiseless signals with f = 100 s-1 (Fig. 7.1, left column) for all depths, the images of objects located at z = 200 and 400 µm for f = 1000 and 10,000 s-1 are evidently broadened and attenuated as compared to corresponding images reconstructed from noiseless signals (Fig. 7.1, center and right). Both effects are more expressed for f = 10,000 s-1 due to smaller SNR. 15 10 5 0 15 10 5 f = 100 s — -inn/-» _-1 f - 1000 S f- 10,000 S z = 50 um i r 0 15 10 5 0 15 10 5 z = 100 um z, = 200 um 1 ' z = 400 um \ 1 ' 0.0 0.2 0.4 0.6 0.0 0.2 0.4 z (mm) 0.6 0.0 0.2 0.4 0.6 Figure 7.5: Average temperature profiles (black lines) and standard deviations (light-gray spots) reconstructed from 10 simulated PPTR signals with different noise realizations. Images of test objects with z1 = 50–400 µm (see labels) are reconstructed from PPTR signals with f = 100 (left), 1000 s-1 (center) and 10,000 s-1 (right). The actual test objects (dark-gray line) are plotted for comparison. In general, the accuracy of determined peak temperature depths Z is similar for f = 1000 and 10,000 s-1. However, for the deepest test object (z1 less accurate than for f = 1000 s-1. 400 µm) Z obtained for f = 10,000 s-1 is markedly Relative image error ? and temperature peak width W as function of sampling rate for PPTR signals with noise are presented in Fig. 7.6. For test objects with z1 = 50 and 200 µm, we obtain the smallest image errors (? = 0.29 and 0.31 for z1 = 50 and 200 µm, respectively) at f = 1000 s-1 (Fig. 7.6a), while the largest sampling rate f = 10,000 s-1 is optimal for object with z1 = 100 µm (? = 0.28). On average, with f increased above the optimal f image error ? increases due to decreased SNR, and with f decreased below the optimal f image error ? decreases due to insufficient information in PPTR signal. Similar trends are observed for width W (Fig. 7.6b). In contrast to results obtained from noiseless signals, low SNR significantly deteriorates images reconstructed from PPTR signals with large sampling rates. While shallow test objects (z1 = 100 µm) benefit from large f, deeper test objects (z1 = 200 µm) are reconstructed more accurately at moderate sampling rates (f = 1000–2000 s-1). Low sampling rates (f = 100–200 s-1) yield robust and comparable reconstruction results for all depths z1, but markedly less accurate as compared to these obtained at moderate sampling rates. o 60 Chapter 7 Sampling rate ,.-« 0.8- a 1". ,''*"'' --*'" -"""" 0.6- '"^i:— ?-, ^ —H —K 0.4- — d— 50 jim / l 0.2- — ¦—100 jim — o—200 jim — •—400 jim 100 1000 f (s-1) 10000 b — d— 50 jim ,1» 200 --¦— 100 jim -- o— 200 jim — •— 400 jim 150 ^» 100 fc.. / _..#'' 50 '-L:- 100 1000 f (s-1) 10000 Figure 7.6: (a) Relative image error ?, (b) full-width at half-maximum W, and (c) peak temperature Tp as a function of sampling rate f for test objects located at z1 = 50–400 µm (see legend) and PPTR signals with noise. Gray line indicates (b) the actual width (50 µm) and (c) the actual peak temperature. Arrows point to optimal values. 7.3 Discussion Temperature profiles reconstructed from PPTR signals augmented with noise are markedly less accurate when large f is used. The results in Fig. 6.5 indicate that high sampling rate (f = 10,000 s-1) is preferred only for shallow object (z1 = 100 µm), while deeper temperature profiles are reconstructed more accurately with moderate sampling rates (f = 1000–2000 s-1). Equal trends are observed also for parameters of temperature profiles (Fig. 6.6). The image error ? as well the excessive width W are more prominent at high sampling rates (f = 5000 and 10,000 s-1) as compared to moderate sampling rates (f = 1000 and 2000 s-1). Both, ? and W, increase with increasing depth Z for these sampling rates, while they are almost independent of Z for small sampling rates (f = 100 and 200 s-1). Thus, low sampling rates can be used to reconstruct deep objects (z1 > 500 µm), where signals with larger f are significantly deteriorated by noise. In general, the optimal sampling frequency depends also on noise amplitude and depth discretization. Considering the results of this study, a good alternative to fixed sampling rate used in PPTR temperature profiling would be a non-uniform spacing in time (i.e., binning) as suggested by Sathyam and Prahl9. This would also reduce the signal vector lengths and kernel matrix size, thus reducing computational costs. 7.4 Conclusions When temperature profiles are reconstructed from PPTR signals with large SNR, higher sampling rates always yield better reconstruction results. But in presence of realistic noise, high sampling rates (f ? 10,000 s-1) are optimal for shallowest temperature profiles (z ? 100 µm) only, while moderate sampling rates (f ? 1000 s-1) are optimal for other temperature profiles. However, for our PPTR system f = 1000 s-1 offers acceptable reconstruction results for all depths. In general, optimal sampling rate depends on specifics of experimental system and properties of the studied samples. Chapter 8 Tissue phantoms It is important for evaluation of PPTR temperature profilometry of human skin to develop reliable phantoms with well defined geometry that suitably mimic optical and thermal properties as well infrared absorption of human skin. The development of tissue phantom for PPTR involves the choice of matrix composition and absorber. In addition, one may also add scattering particles. An excellent review of tissue phantoms for medical optics is provided by Pouge and Patterson54, while Pifferi et al55 suggested directions for designing and characterization of tissue phantoms for biomedical optics. Here we show preparation of agar and collagen tissue phantoms, including gel layer and absorption layer preparation. In addition, we present measured IR spectrum of both types of gel and write corresponding heat diffusion constants. 8.1 Hydrogel layer Since the main component of human skin and both gels is water, matrix composition of tissue phantoms used for PPTR temperature profiling should have similar thermal and IR absorption properties as human skin. Tissue phantoms composed of thin collagen films are a perfect match for real skin.7 Due to good optical and thermal properties collagen and polcrylamide gels were also used for a PPTR temperature profiling.9 In our studies we use agar and collagen gels. 8.1.1 Agar gel Agar gel was prepared by dissolving 0.15 mg of agar powder in 6 ml of distilled water, thus obtaining the agar solution with weight percent of agar 2.5 wt.%. When the agar powder was completely dissolved, the polymerization was initiated by heating the mixture to the boiling point in a microwave oven. Individual gel layers were produced by injecting the agar solution onto a wetted microscope slide with two identical spacers positioned near the ends of the slide (Fig. 8.1). To prevent the entrance of air bubbles, which were present mostly on the surface of the agar mixture, into the agar layer and to precisely control the quantity of the agar solution, a syringe was used. A second slide was placed on top of the agar mixture and gently pressed against the spacers. When polymerization was complete, the top slide was carefully removed, exposing the gel layer of uniform thickness (Fig. 8.2). Figure 8.1: Schematic of agar layer preparation. Figure 8.2: Prepared agar layer with scatterrers (TiO2) on a microscope glass slide. 61 62 Chapter 8 Tissue phantoms Agar layers without scatterrer are transparent to visible light, because they are mostly composed of water. In contrast, agar gel is not transparent to infrared radiation. The sample infrared spectrum is required for PPTR temperature profiling, therefore we measured the IR spectrum of the prepared agar gel using an IR spectroscope. The measured agar gel absorption coefficient 11 as a function of wavelength X is presented in Figure 8.3. The agar spectrum (black line) agrees well with the spectrum of water (gray line) for X = 3-6 µm, while at larger wavelengths protein peaks are present. Another important quantity is the thermal diffusivity constant D of agar gel. Reported values of diffusivities are D = 0.140 ± 0.007,56 0.1427,57 and 0.138 mm2/s 58 for 2.5% agar gel. All reported values are a little smaller when compared to the thermal diffusivity constant of water (D = 0.145 mm2/s at T = 25° C) 59. In our studies, we have found that D = 0.143 mm2/s for agar gel yields optimal reconstruction results. 1000- -------agar gel — pure water 100- r^y^^~ \\ / *J^ "v^L^_^-^/ ^^, 67 /t(ixm) 10 Figure 8.3: Measured IR absorption coefficient ii of the agar gel (black line). Absorption spectrum of water (gray line)48 is plotted for comparison. 7.1.2 Collagen gel To obtain the gelatin solution with 25 weight percent of gelatin, we dissolved 1.25 g of gelatin powder (bovine skin, Sigma-Aldrich) in 3.75 ml of water with 0.2% of formaldehyde. Inclusion of formaldehyde increases the melting temperature of the gelatin matrix by increasing the crosslinking of the fibers while preserving the thermal and spectral properties.54 This allows the collagen phantoms to be used at room temperature. Figure 8.4: Schematic of gelatin layer preparation. We prepared the collagen layer by dissolving the gelatin powder in water at 70° C in water bath.60 When the gelatin powder was completely dissolved under vigorous stirring for 2 min, the homogenous viscous solution was carefully poured onto a microscope slide with two identical spacers positioned near the ends of the slide (Fig. 8.4). The slide was covered by a thin Teflon stripe, since collagen gel sticks on glass substrate very strongly. A second microscope slide covered by a Teflon stripe was placed on top of the gelatin solution and gently pressed against the spacers. When the gelatin was cooled down to room temperature, we first carefully removed the top slide and then the top Teflon stripe to expose the gel layer. Because of the large viscosity of gelatin, we could not use a syringe. 8.2 Absorbing layer 63 Figure 8.5 presents the measured IR spectrum of the gelatin gel. The gelatin spectrum (black line) is about 75% of water spectrum (gray line) at X = 3-6 µm, while at larger wavelengths protein absorption peaks dominate. The thermal diffusivity D of the gelatin gel is assumed to be that of human skin (D = 0.11 mm2/s), based on similar water content to that in human skin. 1000 100 10 1 3456789 10 /t(ixm) Figure 8.5: Measured IR absorption coefficient fi of the gelatin gel (black line). Absorption spectrum of water (gray line)48 is plotted for comparison. 8.1.3 Scatterrer Human skin is a highly turbid medium with the reduced scattering coefficient of dermis61 #s’ = 4 mm1 and of epidermis /*s’ = 12 mm1. To better mimic skin properties, scattering particles were added to our tissue phantoms. Since gelatin and agar phantoms were used for a short period only and then discarded, we used inexpensive titanium dioxide (TiO2) scattering particles. To obtain the above scattering coefficient we mixed 4 mg and 12 mg of TiO2 (Sigma-Aldrich) into 1 ml of gel solution. These concentrations of TiO2 were found experimentally by measuring spectrum of visible light transmitted through thin gel layers with different concentrations of scatterers. Figs. 8.2, 8.7 and 8.9 presents agar layers with TiO2 particles. 8.2 Absorbing layer When we first started preparing tissue phantoms, we prepared the absorbing layers by powdering the surface of the gel layer with small amounts of fine carbon black powder, which was then cowered by another gel layer (Fig. 8.6). Carbon black powder was selected because it is hydrophobic and, therefore, does not diffuse into the gel, enabling preparation of stable thin absorbing layers.49 Figure 8.6: Fine carbon black powder absorbing layer on a agar substrate covered by a thin (~100 µm) agar gel layer with scatterrers. Very thin absorbing layers can be prepared by placing a thin absorbing foil over the hydrogel layer. We used polyethylene foil of thickness ~10 µm. Since foil is crushed during preparation of tissue phantom (see Fig. 8.7), the effective thickness ((~20 µm) is larger than the actual thickness. 64 Chapter 8 Tissue phantoms Figure 8.7: Absorbing foil covered by a thin (~100 µm) agar gel layer with scatterrers. Thicker absorbing layers (d > 30 µm) were prepared as a hydrogel layer with addition of black India ink as the absorber. Both techniques of absorbing layer preparation resulted in absorbing layers of well defined geometry and of homogenous absorber distribution. Figure 8.8: Schematic of a tissue phantom with a single absorbing layer. Water bath prevents the formation of air bubbles between adjacent layers. We put the prepared absorbing and hydrogel layers in water bath, where we constructed the tissue phantom (Fig. 8.8). Water bath effectively prevents formation of air bubbles between the adjacent layers. Finally, the composed tissue phantom was carefully removed from the water bath. Figure 8.9 presents a finished tissue phantom with a single absorbing foil as the absorbing layer. Figure 8.9: Finished tissue phantom with a single absorbing foil as the absorbing layer. If tissue phantoms are not immediately used in an experiment, they must be kept sealed in airtight enclosures such as plastic bags or containers to prevent drying. Keeping the phantoms in vegetable oil has also been reported as an excellent way to preserve the water content.62 Chapter 9 Spectral filtering Numerical simulations in Chapter 6 show how narrowing of the acquisition band improves the validity of quasi-monochromatic approximation (Eq. 2.9) and therefore improves the reconstruction results. However, this causes a reduction in signal, therefore lower SNR, which adversly affects the reconstruction results. Hence, the optimal detection band for given experimental conditions can be determined only by considering the related experimental noise. 9.1 Experiments in agar tissue phantoms 9.1.1 Materials and methods Agar gel tissue phantoms Three tissue phantoms evaluated in this study (samples A, B, and C) consisted of a 1–2 mm thick gel substrate, thin absorbing layer, and one superficial gel layer of varying thickness. Absorbing layers were prepared by powdering the agar gel layer surface with fine carbon black powder. Details of the agar layer preparation were presented in Chapter 8. The subsurface depths of the absorbing layers were approximately 130, 280, and 450 µm, respectively, which corresponded to the location of a vascular network in shallow, medium, and deep port-wine stain birthmarks, respectively. One tissue phantom (sample D) included two absorbing layers at approximate depths of 240 and 440 µm. In this sample, TiO2 powder was homogeneously dispersed in the substrate to enhance optical scattering. The increased light fluence in the deeper absorbing layer resulted in two temperature peaks with comparable amplitudes, despite strong attenuation of the incident laser pulse by the upper absorbing layer.63 Pulsed photothermal profiling For each PPTR measurement, the sample was irradiated with a single 1.5 ms long 585 nm pulse from a pulsed dye laser. Radiant exposure near the center of a 10 mm diameter laser spot was ~3 J/cm2. Radiation emitted from the center of the irradiated area was collected on the focal-plane array of a fast IR camera (Phoenix, Indigo, Santa Barbara, CA, USA) using a macro IR objective with magnification M = 1. By limiting the data read-out to a 128 × 64 pixel sub-window and setting the integration time to tint = 0.5 ms, the acquisition rate was 1083 frames per second. The acquisition time was set to 1 s after the laser pulse. Radiometric signals were obtained from 3 different sites on each sample, separated by a few millimeters to prevent thermal interference between successive measurements. On each test site, up to three radiometric signals were acquired using the entire spectral band of the IR camera (? = 3.0–5.6 µm) and also with a custom long-pass IR filter (cut-on at 4.5 µm, Barr Associates, Westford, MA) fitted to the collection optics.18 The response of each array element was calibrated using a computer-controlled black body (BB701, Omega Engineering, Stamford, CT, USA). Finally, the PPTR signals S (in absolute radiometric temperature units) were obtained by averaging data from 40 × 40 detector elements (active area A = 1.2 × 1.2 mm2) and subtracting the baseline value. Initial temperature profiles T (images) were reconstructed using the monochromatic approximation (Eq. 2.9) and the PR-CG reconstruction algorithm (chapter 4). Elements of the monochromatic kernel matrix K were calculated using the thermal parameter values D = 0.134 mm2/s and h = 0.02 mm-1. The 65 66 Chapter 9 Spectral filtering effective absorption coefficient was determined as /ueff = 28.0 mm1 for broad-band signal acquisition (k = 3.0-5.6 µm) and fieff = 30.2 mm1 for the narrowed spectral band (4.5-5.6 µm). These values were determined from IR spectral properties of the sample (Fig. 8.3) and radiation detector (Fig. 9.1) following the approach presented in Chapter 6. Each reconstruction result consisted of 140 temperature values over a depth range of 0.7 mm (discretization step Az = 5 |am). 1.0 0.8 0^ 0.6 5 0.4 0.2 0.0 3.0 3.5 4.0 4.5 5.0 5.5 A (p,m) Figure 9.1: Relative responsivity of an InSb detector (J10, Judson Technologies, Montgomeryville, PA) in mid-IR spectral region. Optical coherence tomography and histology Several cross-sectional images per sample were acquired using an OCT system64 with a central source wavelength of 1.3 µm. The axial and lateral scanning lengths were set to 800 µm and 2 mm, respectively. The images were saved in JPEG format (400×634 pixels) for further analysis. For purposes of presentation and analysis, the axial dimensions within the sample were corrected using an estimated index of refraction (1.32). Distance from the sample surface to the center of the absorbing layer was determined at six equidistant positions in each image. Finally, thin vertical sections were cut from the sample using a pair of blades and placed onto a clean microscope slide. The sections were inspected under a microscope at magnifications 4, 20, and 40, and photographed using a charge-coupled device (CCD) camera (resolution 552×744). Depth of the absorbing layer was determined from the microphotographs at ten locations per sample. 9.1.2 Experimental results Pulsed photothermal profiling 6 5 4 3 2 1 0 0.0 0.2 0.4 0.6 0.8 1.0 t (s) Figure 10.2: PPTR signals acquired from agar gel samples A and C using the entire spectral band (? = 3.0–5.6 µm; dashed lines) and the reduced spectral band (4.5–5.6 µm; solid lines). Sample A Sample C 9.1 Experiments in agar tissue phantoms 67 Figure 9.2 presents PPTR signals acquired from samples A and C. PPTR signals acquired using the entire spectral band (? = 3.0–5.6 µm; dashed lines) differ in shape from those obtained using the reduced spectral band (? = 4.5–5.6 µm; solid lines). Noise amplitude in PPTR signals is NE?T = 6.5 and 10.9 mK for the full and reduced spectral band, respectively. SNR values are summarized in Table 9.1 for all measured PPTR signals. Table 9.1: Signal-to-noise ratios (SNR) for PPTR signals obtained from four agar gel samples (A–D) using the full (left column) and reduced spectral bands (right column). SNR Sample 1 = 3.0-5.6 µm 1 = 4.5-5.6 µm A 496 245 B 456 285 C 295 155 D 482 277 30- Sample A 20- A 10- I \ 0-30- i * A 20- / \ 10- \ 0-30- // \ A 20- A 10- /\ 0.0 0.1 0.2 0.3 0.4 z (mm) 0.5 0.6 0.7 30- Sample B 20- A 10-0-30- / Y 20- r\ 10- // v-, 0-30- /^ 20- A 10- / \ 0.0 0.1 0.2 0.3 0.4 z (mm) 0.5 0.6 0.7 30- Sample C 20- 10- //~~~ x. 0-30- // N^.. 20- 10- \ 20- 10- /T~ ^v^ X;... 0.0 0.1 0.2 0.3 0.4 z (mm) 0.5 0.6 0.7 40 30 20 10 30 20 10 0 0.0 0.1 0.2 0.3 0.4 z (mm) 0.5 0.6 0.7 Figure 9.3: Reconstructed temperature profiles from three sites on samples A–D. Corresponding PPTR signals were obtained using the full-spectrum acquisition (? = 3.0–5.6 µm, dashed line) and reduced-spectrum acquisition (? = 4.5–5.6 µm, solid line), respectively. 68 Chapter 9 Spectral filtering Figure 9.3 presents temperature profiles reconstructed from PPTR measurements on three different sites on samples A–D, using the entire (dashed line) and reduced spectral band (solid line). Peak temperature depths Z determined on the same site with the two approaches do not differ significantly. The observed variation between different sites on the same sample is due to non-uniform thickness of the superficial gel layer. Yet, temperature profiles obtained with the reduced-spectrum acquisition appear narrower as compared to the full-spectrum acquisition. This effect is particularly evident in samples A and D. 200 150 100 50 - full reduced C Sample D1 D2 Figure 9.4: Average widths of absorbing layers W and standard deviations (error bars) for samples A–C and both absorbing layers in sample D (D1 and D2). Temperature profiles were reconstructed from PPTR signals acquired using the full-spectrum (light gray) and the reduced-spectrum acquisition (dark gray). Table 9.2: Average peak temperature depths (Z) and full-widths at half-maximum (W) of the absorbing layers as determined for the full- (left column) and reduced-spectrum acquisition (right column). Sample X = 3.0-5.6 µm X = 4.5-5.6 µm Z (µm) W (µm) Z (µm) W (µm) A 123 45 126 40 146 55 143 45 135 60 138 50 B 276 95 273 70 283 115 278 90 288 85 288 85 C 448 165 458 140 453 140 458 150 453 185 468 160 D 1st 243 45 223 36 253 53 258 35 208 46 213 43 2nd 443 90 431 77 443 75 453 48 388 67 396 61 Average peak temperature depths Z and full-widths at half-maximums W as determined from all test sites using both experimental approaches are summarized in Table 9.2. Depths Z determined using the full-spectrum reconstruction results do not differ significantly from depths determined from the reduced-spectrum acquisition. In samples A and B, the difference between the corresponding depths is smaller than the temperature profile discretization (?z = 5 µm), and slightly larger in samples C and D (~10 µm). But, lobe widths W are significantly larger when the full-spectrum acquisition is used (Table 9.2, left column), as compared to reduced-spectrum acquisition (right column). 0 9.1 Experiments in agar tissue phantoms 69 This difference in the broadening effect is illustrated in Figure 9.4, where average widths W and standard deviations are presented for all samples and both acquisition approaches. Clearly, widths W for the reduced-spectrum acquisition (dark gray) are smaller as compared to widths for the full-spectrum acquisition (light gray). Optical coherence tomography and histology An OCT cross-sectional image of sample A (Fig. 9.5a) shows clearly the sample surface (upper arrow) and the absorbing layer (lower arrow) due to strong scattering of incident laser light at these two boundaries. Blurring of both boundary lines, which amounts to ~20 µm precludes, accurate and reliable determination of the top layer’s thickness. If the center of each line is selected to represent the boundary location, the average depth of the absorbing layer is determined as 120 µm with a standard deviation of 15 µm. Sample B (Fig. 9.5b) displays the largest variation of the top layer thickness, which is determined as 290 ± 32 µm. Characteristic “ringing” artifacts are present around the surface line in the image of sample D (Fig. 9.5c; top arrow). Nevertheless, both absorbing layers are easily discernible (mid- and bottom arrow, respectively), and their average depths are determined at 240 and 412 µm. Pronounced optical scattering due to TiO2 particles in the gel substrate underneath the deeper absorbing layer is clearly visible in the image. The OCT results from all samples are presented in Table 9.3. Figure 9.5: OCT images of tissue phantoms: (a) sample A and (b) sample B. The upper arrow indicates the sample surface, and the lower one indicates the absorbing layer. (c) Sample D: The middle arrow indicates the first absorbing layer, and the bottom one indicates the deeper absorbing layer. The substrate (underneath the latter) shows scattering due to TiO2 particles. Figure 9.6a presents histology of sample C under an optical microscope. Sample surface is indicated by the top arrow. Carbon powder granules are confined to the boundary between the top agar layer and the thicker substrate layer (note the 500 µm scale bar). The absorbing layer depth varies with location, resulting in an average value of 450 µm with a standard deviation of 30 µm. Histology of sample D is presented in Fig. 9.6b. The upper arrow indicates the sample surface and the middle arrow the first absorbing layer. The optically scattering substrate layer underneath the second absorbing layer (bottom arrow) appears darker in this transillumination microphotograph. The absorbing layer depths as determined from histology in all samples are presented in Table 9.3. 70 Chapter 9 Spectral filtering Figure 9.6: Histology of (a) sample C: The upper arrow indicates the sample surface, and the lower arrow indicates the absorbing layer. (b) Sample D: The upper arrow indicates the sample surface, the middle arrow the upper absorbing layer, and the bottom arrow the deeper absorbing layer. The substrate layer in sample D appears darker because of pronounced optical scattering (microscope objective, 4×). Table 9.3: Average depths and widths of the absorbing layers as determined with three measurement techniques. PPTR,3-5 µm PPTR,3-5 µm OCT W (µm) Histology Z (µm) W (µm) Z (µm) W (µm) Z (µm) Z(µm) 135 ± 12 53 ± 8 136 ± 9 45 ± 5 120 ± 15 20 125 ± 15 282 ± 6 98 ± 15 280 ± 8 82 ± 10 290 ± 32 20 287 ± 15 451 ± 3 163 ± 23 461 ± 6 150 ± 10 472 ± 15 30 450 ± 30 235 ± 24 48 ± 4 231 ± 24 38 ± 4 240 ± 3 40 200 ± 4 425 ± 32 77 ± 12 427 ± 29 62 ± 15 412 ± 7 40 340 ± 4 9.2 Numerical simulation 9.2.1 Methods Initial temperature profiles in our numerical simulation have a hyper-Gaussian form: AT(z, 0) = AT0 exp[-2 (z-z0)4/w04]. To simulate the experimental data, we select the parameter values AT0 = 30 K and w0 = 30 µm and set z0 to 133 µm (object A), 282 µm (object B), and 468 µm (object C). An additional test profile is composed of two hyper-Gaussian lobes of width w0 centered at z1 = 223 µm and z2 = 403 µm (object D). The corresponding vectors T0 consist of 140 values evaluated at equidistant depths within a depth of 0.7 mm. Theoretical signals vectors S0 are calculated from T0 using Eq. (2.11). These have 1083 components which represent PPTR signal values acquired at a sampling rate of 1083 s–1. We simulated different spectral acquisition bands, with the lower wavelength limit k varied from 3.0 to 5.0 µm and the upper wavelength limit fixed at the InSb radiation detector cut-off wavelength (Fig. 9.1), 4 = 5.6 µm. The corresponding kernel matrices K are calculated by dividing each spectral acquisition band into N intervals of width 0.02 µm and adding up their contributions using (6.1), where µ(X) is IR absorption spectrum of agar gel. Each theoretical PPTR signal S0 is augmented by realistic noise. To calculate noise equivalent temperature rise NEAT for each simulated spectral band, we must first determine the total noise amplitude nt. In the following, we apply the parameters of the experimental system: active area of the detector A = 1.44 × 10–2 cm2, collection half angle 6 = 11.3°, frequency bandwidth Af = 1000 Hz, and baseline temperature Tb = 298 K. Peak responsivity of the InSb radiation detector is estimated to Rp = 3.0 A/W. For simplicity, we set both e and C to 1 (Eq. 3.11), because their influence on the simulation results is minimal. Using the relation, which follows from (2.6) and (3.11) 9.2 Numerical simulation 71 NEAT h CsAsin20 [ R(Ä)B;(Tb)dÄ (9.1) we can determine nt from the experimentally determined NEAT values (6.5 and 10.9 mK for spectral bands of X = 3.0-5.6 µm and 4.5-5.6 µm, respectively). The result is almost identical for both spectral bands, nt = 3.0 × 10–10 A, since the shot noise amplitudes nsh computed by using (3.4) and (3.11) are markedly smaller than nt (nsh = 2.0 × 10–11 A and 1.7 × 10–11 A for Xl = 3.0 and 4.5 µm, respectively). Hence, we calculate NEAT using (9.1) with nt = 3.0 × 1010 A for all simulated spectral bands. In addition, the presence of 1/ f noise is characterized by a corner frequency fc = 330 Hz and an exponent a = 1 (Eq. 3.5), so we simulate noise that contains appropriate contributions of zero-mean white noise and 1/ f noise. Initial temperature profiles T are reconstructed from simulated PPTR signals using the monochromatic approximation. Elements of the monochromatic kernel matrix K (2.9) are calculated using the effective absorption coefficient values µeff, determined separately for each spectral band as presented in Chapter 6.1.2. Because reconstruction results are very sensitive to specific noise realizations, each theoretical signal S0 is augmented with 10 different realizations of noise and the results are analyzed statistically. 9.2.2 Simulation results Simulated PPTR signals for test objects A–C are presented in Fig. 10.6 for acquisition spectral bands with ? = 3.0–5.6 µm, 4.5–5.6 µm and 5.0–5.6 µm. The narrowest spectral band (right) presents a significantly larger NE?T as compared to the entire spectral band (left). 2 p A ä= 3.0 - 5.6 jim r\ A ä= 4.5 - 5.6 jim / ^v A ä= 5.0 - 5.6 jim 1 / B ^-^___^ f / B ^^^Z^ZT"";——*-—--^__ B ¦^^^ / / C ^~^——" -~~ C j——-~*~" / C j/*" "^ *^^ 0 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0 t (s) Figure 9.6: Simulated PPTR signals for test objects A-C and spectral bands X = 3.0-5.6 µm (left), 4.5-5.6 µm (center) and 5.0-5.6 µm (right). All signals are augmented by realistic noise. 600 400 200 A...............A~~....A... O...............O------o- O..............-O—-O..............-O- 3.0 3.5 4.0 4.5 2 (um) 20 10 5.0 n 4 B 5 0 0 Figure 9.7: NE?T for the simulated spectral acquisition bands (solid circles) as a function of ?l. SNR of the simulated PPTR signals (open symbols) decreases with ?l for all test objects (A–D). 72 Chapter 9 Spectral filtering Figure 9.7 presents the NEAT values (9.1) as a function of A (solid circles). NEAT increases monotonically from 6.9 mK at Al = 3.0 µm to 21.1 mK at Al = 5.0 µm. Consequently, SNR of simulated PPTR signals decreases with Al for all test objects (open symbols). Effective IR absorption coefficients µeff (6.3) for the simulated spectral acquisition bands are presented in Fig. 9.8. The values µeff vary between 25.0 mm1 and 30.2 mm1, as dictated by spectral dependences of µ(X), R(X) and Bx(Tb) in the corresponding spectral bands. 31 30 »- • 29 28 27 •..._ y* • 26 25 • 3.0 3.5 4.0 4.5 5.0 \ (u.m) Figure 9.8: Effective IR absorption coefficient µeff as a function of Xl for the simulated spectral bands. Figure 9.9 presents statistical analysis of reconstruction results for three spectral acquisition bands: k = 3.0-5.6 µm (left column), 4.5-5.6 µm (center), and 5.0-5.6 µm (right). In each panel, black lines connect the average temperature values and light-gray bars indicate standard deviations. The actual test objects are depicted for comparison (dark-gray lines). 40 20 0 20 Ä = 3.0 - 5.6 |am Ä = 4.5 - 5.6 |am X= 5.0 - 5.6 \xm 0 20 0 20 0.4 z (mm) Figure 9.9: Average temperature profiles (black lines) and standard deviations (light-gray bars) reconstructed from 10 simulated PPTR signals with different noise realizations. Images of four test objects (see the labels) are reconstructed from PPTR signals with spectral bands: ? = 3.0–5.6 µm (left), 4.0–5.6 µm (center) and 5.0–5.6 µm (right). The actual test objects are plotted for comparison (dashed lines). For the reduced spectral bands (center and right), reconstructed temperature profiles of object A (top row) are narrower and higher as compared to the full spectral band (left). Similar trends are observed for objects B and C (2nd and 3rd row, respectively), while object D (bottom row) is reconstructed optimally, when the spectral band with ? = 4.5–5.6 µm is used. But spectral band reduction compromises the stability of the reconstruction results (i.e., increased standard deviation), due to reduced SNR. Figure 9.10a presents relative image errors ? as a function of ?l. For object A (circles), the minimal average error is obtained at ?l = 3.8 µm (? = 0.088; note the arrow). With ?l increased above ?l = 4.0 µm, ? increases progressively due to decreasing SNR (? = 0.19 at ?l = 5.0 µm). With ?l decreased 9.2 Numerical simulation 73 below 3.8 µm, S increases monotonically due to increasing deficiency of the monochromatic approximation, reaching S = 0.14 at Xl = 3.0 µm. The same trend is present also for standard deviation of S. Standard deviation presents a minimum at kl = 3.8 µm (as = 0.02) and significantly increases when spectral band is ether broadened (o* = 0.05 at Xl = 3.0 µm) or narrowed (o* = 0.10 at Xl = 5.0 µm). 1.0 0.8 0.6 0.4 0.2 0.0 0.4 0.2 0.0 3.0 3.5 4.0 4.5 5.0 150 100 50 100 50 D2 D1 //V^ 4 ±= -I—I -0-----0—------o-----o—------§¦-----6--------A-— 3.0 3.5 4.0 1 (um) 4.5 5.0 Figure 9.10: (a) Relative image error ? as a function of ?l for test objects A (circle), B (square), C (diamond) and D (triangle). (b) Analogous results for the full-width at half maximum of reconstructed temperature peaks (W). The widths of both peaks are analyzed for object D (D1, D2; bottom). Standard deviations are presented as error bars. Dotted lines in (b) indicate the actual object width, W = 47 µm. Arrows indicate minimums of ? and W. Similar trend is observed for object B (Fig. 9.10a, squares), where the minimum S = 0.36 at h = 3.8 µm is markedly larger as compared to test object A due to the broadening. Standard deviation 0-8 tends to increase with Xl (