In vivo port-wine stain depth determination with a photoacoustic probe

We have designed a photoacoustic probe for port-wine stain (cid:1) PWS (cid:2) depth measurements consisting of optical ﬁbers for laser light delivery and a piezoelectric element for acoustic detection. We characterized the capabilities and limitations of the probe for proﬁling PWS skin. The probe induced and measured photoacoustic waves in acrylamide tissue phantoms and PWS skin in vivo . The optical properties of the phantoms were chosen to mimic those of PWS skin. We denoised acoustic waves using spline wavelet transforms, then deconvolved with the impulse response of the probe to yield initial subsurface pressure distributions in phantoms and PWS skin. Using the phantoms, we determined that the limit in resolv-ing epidermal and PWS layers was less than 70 (cid:3) m. In addition, we used the phantoms to determine that the maximum epidermal melanin concentration that allowed detection of PWS was between 13 and 20%. In vivo measurements of PWS skin with different epidermal melanin concentrations correlated with the phantoms. Thus the photoacoustic probe can be used to determine PWS depth for most patients receiving laser therapy. © 2003 Optical Society of America


Introduction
Laser treatment of port-wine stain ͑PWS͒ lesions is an effective therapy to restore normal appearance to human skin, although the success rate depends on many factors, including epidermal melanin concentration, lesion depth, and size of the blood vessels.Epidermal melanin, a broadband optical absorber, blocks laser energy and decreases fluence at the lesion.For human skin with a high concentration of melanin ͑skin types III or higher 1,2 ͒, epidermal temperature may rise to such levels that irreversible damage and scarring or dyspigmentation occurs.As a preventive measure, clinicians have utilized skin cooling, such as cryogen spray cooling, 3,4 to prevent epidermal damage while still treating the deeper PWS.With cryogen spray cooling, a cryogen spray is directed onto the skin surface prior to and during laser irradiation, and the epidermis cools nearly instantaneously.The PWS cools at some later time, preferably long after the therapeutic laser pulse has been delivered, causing irreversible thermal damage to the lesion.Knowing the depth profile of the PWS skin, including the spatial relationship between epidermal melanin and blood vessels, the clinician can optimize cryogen spray cooling parameters for treatment on an individual patient basis.Thus we have developed a photoacoustic probe for determining PWS depth quickly and noninvasively.
Photoacoustic depth profiling, as described in this paper, uses pulsed laser irradiation to induce rapid thermoelastic expansion in targeted chromophores.This process is distinct from photoacoustic methods that use modulated continuous-wave irradiation, such as photoacoustic spectroscopy. 5Photoacoustic generation by thermoelastic expansion can be conceptually described as laser energy being quickly absorbed by a small volume such that the resultant heating induces rapid expansion that manifests itself as a transient pulse of acoustic energy.Thermoelastic expansion occurs when the condition of stress confinement is achieved, i.e., where optical energy is deposited before the energy can propagate away acoustically.This condition is expressed as t p Ͻ ␦͞c s , where t p is the laser-pulse duration, ␦ is the absorption depth of laser energy, and c s is the speed of sound in the medium.If the radiant exposure is not excessive, the resulting acoustic waves behave according to the linear wave equation.Furthermore, if the laser spot diameter is much larger than ␦, then a simple plane-wave analysis can be used, allowing acoustic propagation time to be used as an indicator of distance traveled.If stress confinement, linearity, and plane-wave geometry are preserved, depth profiling and imaging of layered tissue can be achieved by simple photoacoustic analysis.
Use of photoacoustic techniques is ideally suited for the application used to determine skin structure because of the robustness and high resolution of the acoustic signal.Hence analysis of photoacoustic waves has been used to investigate layered structures in biological tissue.Oraevsky et al. 6 derived optical properties of absorbing solutions by fitting the photoacoustic signal to Beer's law.Viator et al. 7 derived an iterative scheme to determine the absorption coefficient of layered gels and stained elastin from photoacoustic waves.Viator et al. 8 developed an endoscopic photoacoustic probe for determination of treatment depth after photodynamic therapy for esophageal cancer.Viator et al. 9 developed a probe to determine the epidermal and PWS blood vessel depths in human subjects using a photoacoustic reflectance probe.We present an evaluation of this photoacoustic probe used for PWS depth to determine the limitations of its use in human patients.
The photoacoustic probe, consisting of a small acrylic handpiece that combines two optical fibers and a piezoelectric element, was used to map the depth distribution of epidermal melanin and PWS.We present measurements with tissue phantoms that show the minimum separation with which these two layers can be discriminated.We also use phantoms to show the maximum epidermal melanin content that allows threshold detection of the PWS layer.The raw photoacoustic signals from phantoms and human subjects were denoised by use of wavelet transforms and deconvolved with the probe's impulse response to give the approximate initial pressure distribution after the laser pulse.

A. Probe Design
The probe is shown schematically in Fig. 1 and in a photograph in Fig. 2. It consists of two 1500-mdiameter optical fibers terminating into a cylindrical acrylic handpiece that is 22 mm long and 16 mm in diameter.An acoustic transducer element within the probe sent a signal through a 10-⍀ coaxial cable to an oscilloscope ͑Tektronix TDS 3014, Wilsonville, Ore.͒ with a bandwidth of 100 MHz.The input impedance of the oscilloscope was 1 M⍀.The laser spot from the optical fiber was slightly elliptical with a diameter of approximately 2 mm.The radiant exposure at the surface of the targets was 0.70 J͞cm 2 .The optical fibers were steered by set screws within the acrylic handpiece so that both laser spots were coincident.The probe was placed in contact with the targets so that the photoacoustic waves were generated directly below the acoustic detector.A 1.1mm-diameter semirigid coaxial cable ͑UT-43-10, Micro-Coax, Pottstown, Pa.͒ was inserted into the probe with a piezoelectric film ͓polyvinylidene fluoride ͑PVDF͒, K-Tech, Albuquerque, N.M.͔ attached to the end.The active area of the acoustic sensor was 700 m in diameter.Thus the optical fibers irradiated the tissue surface, inducing photoacoustic waves, which were sensed in reflection mode by the PVDF sensor.The fibers and sensor were housed in a water-filled chamber, which allowed for acoustic impedance matching between the target surface and the sensor.The acoustic element was recessed ap-Fig.1. Photoacoustic probe consisting of two 1500-m-diameter optical fibers and a PVDF sensing element.The output ends of the fibers and the PVDF were set within a chamber that was filled with water for acoustic coupling to the target surface.Set screws were used to steer the fiber ends to direct the laser spot directly beneath the PVDF.Fig. 2. Photograph of the probe design.The rigid coaxial cable can be seen extending into the inner chamber of the probe.proximately 3 mm into the probe housing to separate the target surface from the sensor.This separation created an acoustic delay line of approximately 2 s that prevented electrical noise, caused by the laser pulse and occurring at 0 -0.5 s, from contaminating the photoacoustic signal that was transmitted via the coaxial cable to an oscilloscope.

B. Apparatus
The experimental apparatus is shown in Fig. 3 and consisted of a frequency-doubled Nd:YAG laser operating at 532 nm with a 4-ns ͑FWHM͒ pulse duration.Laser light was coupled into two optical fibers by a cube beam splitter.Two fibers were used to increase the total amount of energy delivered, although energy was limited to less than 22 mJ to keep patient discomfort to a minimum.Subsequent photoacoustic analysis was done on a Macintosh computer ͑G4 Powerbook, Apple Computer, Cupertino, California͒.
We calibrated the transducer using the procedure in Viator et al. 9 We calibrated the detector by inducing photoacoustic waves in solutions where the absorption coefficient was known.The PVDF sensor detected the waves in a transmission setup.The free beam of the laser was used, providing a large spot that minimized diffraction and delivered more energy.The detectors were immersed in absorbing solutions and centered directly above the laser spot.
The laser spot diameter was 4.6 mm.The radiant exposure was 0.084 J͞cm 2 , which we calculated by measuring the total energy with a detector ͑Molectron, Beaverton, Ore.͒ and dividing by the spot size.The absorption coefficients of the solutions were 51, 103, 148, 197, and 239 cm Ϫ1 at 532 nm.The equation was used to predict the photoacoustic pressure ͑joules per cubic centimeter͒.⌫ is the Gru ¨neisen coefficient, which models the fraction of optical energy converted to acoustic energy.In this analysis, ⌫ ϭ 0.12 was used. 6H 0 is the radiant exposure ͑joules per square centimeter͒, and a is the absorption coefficient of the solution in inverse centimeters.Finally, we used the conversion of 10 bars ϭ 1 J͞cm 3 to determine a calibration factor of millivolts per bar for the acoustic detector by dividing the amplitude of the acoustic waveform by the calculated pressure.The calibration factor was 1.31 mV͞bar.

C. Layered Acrylamide Phantoms
Tissue phantoms were made with 20% acrylamide in water with added dye for absorption and fat emulsion for scattering.The amount of dye and fat emulsion was chosen to mimic optical properties of human skin.Direct Red 81 ͑Sigma Chemical, St. Louis, Mo.͒ was used to simulate hemoglobin absorption, and Intralipid ͑Abbott Laboratories, North Chicago, Ill.͒ was used to induce light scattering.An Intralipid-1% solution was used to approximate 200 cm Ϫ1 , the approximate scattering coefficient of human skin at 532 nm. 10 -12 The procedure for preparing acrylamide is fully discussed in Viator et al. 9 The two types of phantom prepared consisted of three layers, representing epidermis, bloodless dermis, and PWS.Acrylamide solutions were injected between glass slides with plastic feeler gauge stock of various thicknesses used as spacers ͑feeler gauge stock, McMaster-Carr, Los Angeles, California͒.In the first set of experiments, the minimum distance that the probe could discriminate between the two absorbing layers ͑epidermis and PWS͒ was studied.These epidermal and PWS layers were entirely absorbing, with absorption coefficients of 25 cm Ϫ1 , and the bloodless dermis layer was clear.In the second set of experiments, phantoms were turbid.These experiments were meant to determine the maximum epidermal melanin concentration that would allow detection of a photoacoustic signal from the PWS.Each layer experiment was performed once on each phantom, although the reported measurements were averages over 64 laser pulses.

Layer Discrimination
The epidermal layers were 100 m thick, and the PWS layers were approximately 1 mm thick.The bloodless dermis layers were 70, 140, 350, 450, 770, and 1030 m thick.The three layers were placed on top of each other.Enough moisture was present to ensure acoustic coupling, and any air pockets were removed.After laser irradiation, layers were measured with a micrometer ͑Digital Micrometer, Mitutoyo, Aurora, Ill.͒ with an accuracy of approximately 20 m.The epidermal thickness was subtracted from photoacoustically derived thicknesses to deduce the thicknesses of the bloodless dermis layers.Photoacoustic peak-to-peak times were derived from the epidermal surface to the PWS layer; thus the depths indicated by the peak-to-peak time included the epidermal thickness and needed to be subtracted.

Epidermal Melanin Filtering
To model skin with various epidermal melanin concentrations, turbid phantom layers were used with a scattering coefficient of approximately 200 cm Ϫ1 . 10The optical absorption was varied to model 2, 5, 13, and 20% epidermal melanin volume fractions, 1 corresponding to typical values for skin types I, II, III, and IV.We assumed a melanin absorption coefficient of 400 cm Ϫ1 at 532 nm.The optical depths corresponded to 0.06, 0.15, 0.39, and 0.60, as the layer thickness was 75 m.Optical properties and layer thicknesses of the phantoms are shown in Table 1.

In vivo Measurements
We tested three healthy human subjects using the photoacoustic probe.All institutional rules governing human subjects were strictly complied with.The subjects had skin phototype I-II, III, and IV.PWS birthmarks were on the face and arms.All subjects indicated little or no sensation during laser irradiation.Individual tests took approximately 10 min and preceded normal laser therapy for PWS.

D. Signal Processing
The raw signal from the photoacoustic probe was a convolution of its impulse response and actual pressure induced by the laser pulse.This assertion is reasonable in that photoacoustic generation and propagation are described by the linear wave equation, 13,14 and hence the probe constitutes a linear, time-invariant system.Analysis of waveforms from phantoms of known dimensions showed a broadening, or smearing, from what would be expected from the pressure signals.This smearing was also evident in the raw signals from human subjects, as the breadth of the epidermal signal exceeded the normal range of skin thicknesses.Smearing was due to the relatively large photoacoustic source, caused by the large laser spot, compared with the small active area of the acoustic detector.Such an extended source is detected by the small detector as a single signal over a time period related to the dimensions of the source ͑laser spot diameter͒.After determining the probe impulse response, we performed a simple deconvolution on all signals to produce the initial pressure distribution of each experiment.

Signal Denoising
Prior to deconvolution, we denoised the photoacoustic signals to obtain optimal deconvolution using our algorithm.We achieved denoising by two means: signal averaging during the experiment and postexperiment using wavelet shrinkage techniques.During the experiment, the signals were averaged over 64 laser pulses to minimize random noise.Longer averages were not taken because of dynamic pro-cesses that could change the photoacoustic signal, such as subject movement.
We accomplished further denoising with wavelet transforms using Wavelet Explorer ͑Wolfram Research, Inc., Urbana, Ill.͒, an add-on of MATHEMATICA.Although research in wavelet transforms is relatively recent, many references exist.An elementary but well-written treatment is by Walker. 15A concise introduction to wavelets in a biomedical context is found in Thakor and Sherman. 16A longer but selfcontained treatment is found in Kaiser. 17Daubechies's volume 18 is extensive and sophisticated, but is not oriented toward applications and requires many notions from functional analysis.Wavelet shrinkage for denoising is treated in Donoho and Johnstone. 19The specifics of wavelet denoising by soft thresholding is presented in Appendix A.
Spline wavelets were chosen because the expected pressure signal was suited to relatively low-order polynomial fits, verified by visual inspection of the noisy signals.The denoising algorithm used fourlevel spline wavelet transforms, and we obtained the threshold level by estimating the noise level on each signal.We selected the threshold by taking a value between the noise level and the smallest signal variation, with the threshold set closer to the noise level.An example of wavelet denoising of a photoacoustic signal is shown in Fig. 4.

Deconvolution of Probe Impulse Response
We determined the initial pressure distribution generated by absorption of laser light by deconvolving the photoacoustic signal with the probe impulse response, which we determined by irradiating a highly absorbing acrylamide gel with the probe.Ideally, an impulse response would require a target of infinite absorption, but, practically, an absorption coefficient higher than the ability of the system's resolution would be sufficient.Probe resolution is limited by the laser-pulse duration p .By use of the speed of sound in tissue c s , approximately 1.5 mm͞s, the resolution limit that is due to the laser-pulse duration is ⑀ p ϭ c s p ϭ 6 m.An absorption depth of 6 m corresponds to an absorption coefficient ͑ a ͒ of approximately 1700 cm Ϫ1 .We used an acrylamide phantom of 1 mm thickness with a a of 1500 cm Ϫ1 because making a phantom of higher absorption was not possible with Direct Red.Radiant exposure to determine the impulse response was approximately 0.5 mJ͞cm 2 .The resulting impulse response is shown in Fig. 5. Diffraction is evident beginning at 2.3 s, although this region was not included in the deconvolution algorithm.
Once the impulse response was determined, we deconvolved the photoacoustic signals by taking the fast-Fourier transforms ͑FFTs͒ of the impulse response and signal: where h͑t͒ is the impulse response and p͑t͒ is the photoacoustic signal.Deconvolution was performed by a smooth division by zero routine 20 : where S͑͒ is the FFT of the true pressure signal and ␦ is an arbitrarily small quantity.For these experiments, ␦ Ϸ 10 Ϫ4 of the transform values.This scheme was used to prevent division by zero from regions of the impulse response's spectrum with no frequency information.The result of this division scheme at such points would result in S͑͒ ϭ 0. The deconvolution scheme was implemented in MATH-EMATICA.

A. Layer Discrimination
A montage of photoacoustic waves is shown in Fig. 6 representing laser irradiations of three layered phantoms with clear layer separations of 70, 140, 350, 450, 770, and 1030 m each.The separation of acoustic peaks corresponds to bloodless dermis layer thickness.
The thicknesses of the phantoms, as determined by the micrometer measurements, are compared with the photoacoustically determined thicknesses in Table 2.

B. Epidermal Melanin Content
A montage of the photoacoustic waves from the epidermal filtering experiments is shown in Fig. 7.The second peak, representing the PWS layer, is evident Fig. 5.We determined the impulse response of the photoacoustic probe by irradiating an acrylamide phantom with a ϭ 1500 cm Ϫ1 .The oscillation prior to 1 s is due to electrical noise from the laser.Fig. 6.Separations of the epidermal and PWS peaks.The PWS peak decreases in amplitude as it merges with the negative diffractive signal from the early peak ͑epidermal signal͒.The bloodless dermis thickness is indicated by the value at the top of each graph, and deeper PWSs correspond to greater separation of the two peaks.
in the samples with 2, 5, and 13% melanin concentrations.The PWS layer is obscured in the sample with 20% melanin concentration.

C. In Vivo Measurements
In vivo photoacoustic measurements of PWS in human skin are shown in Fig. 8. Skin types I-II, III, and IV are shown.

Discussion
The photoacoustic probe described in this paper is well suited for determining PWS depth and its relation to the overlying epidermal melanin layer.Photoacoustic propagation is robust in tissue, showing little attenuation and scattering in the depths of a PWS lesion.However, the limits of the probe to investigate PWS depth are dependent on epidermal melanin concentration.Thus we performed experiments on tissue phantoms to determine these limits.

A. Phantom Measurements
Tissue phantoms were used so that the accuracy of depths of absorbing layers could be determined.The typical discrepancy between the actual depth measured by micrometer and the photoacoustic probe was approximately 5%.The greatest discrepancy was less than 13%.The minimum depth discrimination between layers was 70 m because acrylamide layers could not be made any thinner.The actual ability of the probe to discriminate absorbing layers may be 50 m or less, as 70-m phantoms had photoacoustic peaks clearly distinguishable from each other.The theoretical limit of depth discrimination is dependent on the laser-pulse duration.With a pulse duration of 4 ns, the resolution limit is 6 m because the speed of sound in tissue is approximately 1.5 mm͞s.In practice, discriminating the layers also involves signal duration, which depends on absorption depth of a layer and the actual layer thickness.Extrapolating from the phantom experiments, such a limit is approximately 50 m.
Although the phantom experiments showed promising results for depth profiling, it is important to consider the differences between an idealized phantom, in which anatomical structures are represented by planar layers, and the in vivo case.PWS lesions are composed of many individual blood vessels within the dermis and only approximate a layer.In addition, the epidermal-dermal junction is not necessarily planar, but may have numerous papillae, the dimensions of which are comparable to the epidermal thickness.These nonplanar structures would contribute to the differences between the phantom measurements and the in vivo measurements, where diffraction may make substantial contributions to the    Results of the epidermal melanin filtering ͑Fig.7͒ showed that the probe could detect the PWS signal for skin with melanin volume fraction of up to 13%.However, the PWS signal was obscured by a phantom simulating a melanin volume fraction of 20%.The second photoacoustic peak in the 20% signal was not the PWS layer, as the distance between layers indicated by the time separation of the two peaks was approximately 130 m, or approximately half of the true separation.The second peak was merely the PWS signal being obscured by diffraction from the epidermal melanin layer.

B. In Vivo Measurements
In vivo measurements were taken from three PWS patients with skin types I-II, III, and IV.Measurements for skin type I-II were previously reported by Viator et al., 9 with the results compared with optical Doppler tomography.The same probe was used for the current in vivo measurements.Actual depth relationships between epidermal melanin layers and PWS lesions were not determined by biopsy, but the acoustic waveforms were consistent with turbid phantom experiments.The temporal delays in the photoacoustic signals from the epidermal melanin and PWS layers corresponded with the phantoms, as well as their relative amplitudes.Ideally, photoacoustic signals, after denoising and deconvolution, should be an accurate representation of PWS anatomy in vivo.

C. Signal Processing
The signal processing in these experiments was performed to reconstruct the initial pressure distributions in tissue and phantoms immediately after laser irradiation.Raw data from the photoacoustic probe clearly did not represent the initial pressure because signals from the layered phantoms indicated layer thicknesses that were greater than the actual phantom materials.In addition, raw signals from in vivo PWS skin showed epidermal melanin layers thicker than the known thickness of the entire epidermis.We performed deconvolution of the raw data with the instrument response function, resulting in a better approximation of the actual pressure distribution, although the deconvolution scheme did not give an entirely accurate reconstruction.We improved the deconvolution scheme by increasing the signal-tonoise ratio using wavelet denoising, although reconstructions were still not perfect.Reconstruction of absorbing tissue phantoms still showed broader signals than expected, although better than the raw signals themselves.
Wavelet thresholding successfully denoised the signals, while still preserving salient features.Threshold was chosen after visual inspection of the signal and estimation of the noise level.We performed denoising using different threshold levels, above and below the estimated noise level.Most thresholds were set at approximately 2-3 times the noise level.A new deconvolution scheme may be implemented in the future to obtain more accurate reconstructions of the photoacoustic pressure distribution.

Conclusion
The photoacoustic probe is suited for depth measurements of PWS subjects with skin types I-III and some subjects with skin type IV, as it has the ability to resolve lesions from epidermal melanin to within 70 m.The probe is noninvasive and causes no patient discomfort because of the low laser fluence.Robust signals from tissue phantoms at a 325-m depth indicate that the probe could be used for most PWS patients, given the limitation of epidermal melanin content.Laser treatment of PWS patients with skin type IV or greater is difficult, so the current probe can perform PWS depth determination for most patients receiving laser therapy.The probe successfully detected PWS lesions in patients with skin types I-IV, giving a suggested depth relationship between the epidermal melanin layer and the PWS.
Finally, better reconstruction algorithms could be implemented to improve on the probe's ability to determine the actual depth profiles of subsurface chromophores, ultimately useful in possible photoacoustic imaging of skin and its pathology.

Appendix A: Denoising Data by Wavelet Soft Thresholding
One of the most useful applications of the wavelet transform is denoising data.Because of the nature of the transform, denoising is often more effective than traditional filtering with Fourier transforms.Denoising with Fourier methods involves cutting out frequency components, thus destroying components of the signal along with noise in those frequency ranges.With wavelet denoising, signal energy is preserved as much as possible, removing only those components of the transform that exist beneath a certain threshold.This method effectively preserves signal structure, while selectively decimating small fluctuations associated with noise.Choosing the threshold is of prime importance, although an effective threshold can be chosen by simple inspection of the noisy signal.
Wavelet theory is described in the references mentioned above, but wavelet denoising with soft thresholding, or wavelet shrinkage, can be succinctly explained by means of the discrete wavelet transform.The material in this appendix is a condensation of important concepts from Donoho and Johnstone, 19 Taswell, 21 and Walker. 15The theory of wavelet denoising is simple and consists of three steps: where X is noisy data, W is the wavelet transform operator, Y is its wavelet transformed data, D is the denoising function that can employ hard or soft thresholding, is the threshold level, Z is the wavelet coefficients after thresholding, and S ˆis the denoised data.The concept can be explained by use of the discrete Haar transform.Although the Haar transform is not usually used for practical purposes, it is ideally suited to explain the fundamentals of wavelet theory.

Windowed Fourier Transform versus Wavelet Transform
The wavelet transform is similar in concept to the windowed Fourier transform ͑WFT͒.With the notation of Daubechies, 18 the WFT of a square integrable function f can be written as where g is a windowing function, which localizes the analysis of the Fourier transform in time.Analogously, the continuous wavelet transform can be written as where is the wavelet and a and b are dilation and translation parameters, respectively.Thus a allows the wavelet transform to change the frequency range, whereas b translates the wavelet in time.The similarity between the WFT and wavelet transform is obvious; analysis can be performed at different neighborhoods in time.The wavelet transform, however, not only translates in time, but each time window is analyzed with different frequency ranges.The WFT has a constant-shaped window, whereas the wavelet transform has ". . .time-widths adapted to their frequency." 18This basic fact allows wavelet methods to achieve superior results in applications such as the denoising of signals.

Discrete Haar Transform
Wavelet transforms can be accomplished with a number of different wavelets, but the Haar wavelet is the simplest and most useful for building an understanding of the subject.The notion of the Haar wavelet is built on the simple geometry of the Box function ͓Box͑t͒ ϭ 1 t a Ͼ t Ͼ t b , 0 otherwise͔.Discrete wavelet transforms in one dimension take a signal and divide it into two other signals, each half of the length of the original signal.The first signal is the trend and the second signal is the fluctuation.The trend is a running average of the original signal, whereas the fluctuation is a running difference.These two operations correspond to Haar scaling functions and Haar wavelets.Successive application of the scaling function and wavelet along the length of the signal give the trend and fluctuation, respectively, which can be written as where f is the original signal and a 1 is the trend and d 1 is the fluctuation of the first-level Haar transform.The symbol ͉ denotes concatenation.Subsequent levels are obtained when the Haar transformation is run on successive trend signals, such as so that a two-level Haar transform is represented by As an example, consider the Heavisine function of Donoho and Johnstone. 19The function is given by As can be seen, the first 1024 elements constitute the signal trend, which closely resembles the original signal, and the second 1024 elements constitute the signal fluctuation.To perform denoising on a signal, these transform elements are applied to a thresholding function, described in Section 3 of this appendix.

Hard Thresholding
We performed hard thresholding on wavelet transformed data to remove noise by cutting out transform data smaller than a predetermined threshold value.The hard thresholding function is given by where c i is the ith wavelet coefficient and is the threshold level.This thresholding function is discontinuous and usually results in inferior denoising results.After thresholding is applied, the inverse wavelet transform is performed, yielding a denoised signal.The hard threshold denoised version of Heavisine is shown in Fig. 12͑a͒.

Soft Thresholding
In contrast to hard thresholding, soft thresholding not only eliminates wavelet coefficients below threshold, but it provides a smooth transition in the coefficient values near threshold.Therefore, well above threshold, the coefficient values are unchanged, well below threshold they are set to zero, and in between the values are smoothly forced toward zero.This difference in the hard and soft thresholding can be seen in Fig. 13.A typical shrinkage function and the one used in this paper is The result of soft thresholding the noisy Heavisine can be seen in Fig. 12͑b͒.Using the same threshold value, we reduced the noise level in both cases, although soft thresholding shows a signal much closer to the noiseless Heavisine function.This result is the case for most continuous functions.A counterexample for the superiority of soft thresholding would be a composite function of rectangular functions.In such a case, although the function is continuous almost everywhere, the more prevalent the discontinuities, the better the denoising with hard thresholding would be owing to the discontinuous nature of hard thresholding.However, in most cases of experimental interest, including photoacoustic signals in tissue, the functions are continuous and hence are more suitable for soft thresholding.

Fig. 7 .
Fig. 7. Epidermal melanin layers were modeled by means of absorbing superficial layers in the tissue phantoms.Although PWS signals are evident for the ͑a͒ 2, ͑b͒ 5, and ͑c͒ 13% melanin concentration, at ͑d͒ 20% melanin concentration the PWS is not observed.

Fig. 8 .
Fig. 8. Four PWS patient measurements.Epidermal and PWS signals are shown for skin types ͑a͒ I-II, ͑b͒ III, and ͑c͒ IV, although one measurement of skin type ͑d͒ IV is difficult to interpret.
photoacoustic signal.A comparison of the waveforms in Figs.7 and 8illustrates this difference, where the in vivo measurements show considerable diffraction, particularly in the latter part of the signal.Phantoms were constructed with optical properties closely mimicking human skin.The photoacoustic signal from a three-layer phantom is shown next to a signal from in vivo human PWS skin ͑Fig.9͒.The phantom was made to simulate 5% melanin volume fraction, corresponding to skin type II.The signal from the PWS skin was taken from a human subject with skin type I-II.Although the actual pressures are different between the phantom and the PWS skin, the two signals are similar.The pressures are different because the signals were taken at different times and the probe was not necessarily optimized during the phantom measurement.Nevertheless, the similarity justifies use of turbid acrylamide phantoms to simulate PWS skin.Results of the epidermal melanin filtering ͑Fig.7͒ showed that the probe could detect the PWS signal for skin with melanin volume fraction of up to 13%.However, the PWS signal was obscured by a phantom simulating a melanin volume fraction of 20%.The second photoacoustic peak in the 20% signal was not the PWS layer, as the distance between layers indicated by the time separation of the two peaks was approximately 130 m, or approximately half of the true separation.The second peak was merely the PWS signal being obscured by diffraction from the epidermal melanin layer.

Fig. 9 .
Fig. 9. Photoacoustic signals from ͑a͒ the turbid acrylamide phantom with 5% epidermal melanin content and ͑b͒ PWS type I-II human skin.Although the depths and pressures are different, the signals show similar relationships in the epidermal and PWS profiles.
)where sgn is the signum function.The function is shown is Fig.10͑a͒.Although Haar wavelets are not the best choice for a function that is continuous al-most everywhere, Heavisine was chosen to illustrate the advantage of soft thresholding over hard thresholding.The function was composed of 2048 data points.Random noise was added to the function with an average amplitude of 0.3 ͓Fig.10͑b͔͒.The one-level Haar transform uses the coefficients ͓͑1͞ ͌ 2͒, ͑1͞ ͌ 2͔͒ for the running average and difference, so that successive elements are added and subtracted, respectively, to yield the trend and fluctuation signals, each half of the duration of the original signal.The ͌ 2 is introduced into the coefficients to preserve signal energy.The one-level Haar transform of Heavisine is shown in Fig.11.

Fig. 10 .
Fig. 10.͑a͒ Heavisine function used to illustrate the Haar transform.͑b͒ The noisy Heavisine function created when random noise is added to Heavisine.

Fig. 11 .
Fig. 11.One-level Haar transform of noiseless Heavisine shows the trend signal, appearing as a compressed version of Heavisine, followed by the fluctuation signal.The fluctuation signal is nearly featureless because the Heavisine is continuous almost everywhere.The fluctuation signal at approximately the element number of 1400 is due to the discontinuity of Heavisine at element number 800.

Fig. 12 .Fig. 13 .
Fig. 12.The denoised signal was performed on the noisy Heavisine signal ͑a͒ with hard thresholding and ͑b͒ with soft thresholding.Soft thresholding gives better denoising, even though the same threshold level was used in both denoising functions.