13 October 2017 Noncontact holographic detection for photoacoustic tomography
Author Affiliations +
A holographic method for high-speed, noncontact photoacoustic tomography is introduced and evaluated. Relative changes of the object’s topography, induced by the impact of thermoelastic pressure waves, were determined at nanometer sensitivity without physical contact. The object’s surface was illuminated with nanosecond laser pulses and imaged with a high-speed CMOS camera. From two interferograms measured before and after excitation of the acoustic wave, surface displacement was calculated and then used as the basis for a tomographic reconstruction of the initial pressure caused by optical absorption. The holographic detection scheme enables variable sampling rates of the photoacoustic signal of up to 50 MHz. The total acquisition times for complete volumes with 230 MVoxel is far below 1 s. Measurements of silicone and porcine skin tissue phantoms with embedded artificial absorbers, which served as a model for human subcutaneous vascular networks, were possible. Three-dimensional reconstructions of the absorbing structures show details with a diameter of <inline-formula> <math display="inline" xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mrow> <mn>310</mn> <mtext>  </mtext> <mi>μ</mi> <mi mathvariant="normal">m</mi> </mrow> </math> </inline-formula> up to a depth of 2.5 mm. Theoretical limitations and the experimental sensitivity, as well as the potential for <i>in vivo</i> imaging depending on the detection repetition rate, are analyzed and discussed.



Photoacoustic imaging (PAI) has become very popular in recent years, especially for noninvasive imaging of biological tissue.1 The combination of optics and ultrasound enables deep tissue imaging and overcomes fundamental restrictions of optical imaging in scattering tissues.2 PAI is based on the conversion of pulsed- or intensity-modulated continuous-wave (CW) electromagnetic waves into thermoelastic acoustic waves by absorbing chromophores. The wavelength-dependent optical absorption is the underlying contrast mechanism of PAI, which enables imaging blood vessels,3 vascularization within a tumor, mammography screening,4 and functional imaging of oxygen metabolism5 by utilizing the individual absorption of naturally occurring chromophores.

In contrast to electromagnetic waves, the emitted pressure waves propagate through the homogeneous tissue with much lower refraction and scattering and can be detected at the specimens’ surface. By tomographic imaging, the position and the shape of the absorbing structures can be reconstructed in three dimensions. Different reconstruction techniques were demonstrated. The most commonly used are time reversal algorithms based on the acoustic wave equation,6,7 temporal backward projection using Fourier transform methods,8 and triangulation-based algorithms like delay-and-sum.9 For the detection of the acoustic wave, basically two principles exist.

The first detects the pressure with an adequate ultrasonic transducer. Here, acoustic contact for impedance matching, which defines the application field, is required. Piezoelectric transducers are widely used as a single element in scanning mode with high spatial resolution10 or parallelized in detector arrays like in conventional ultrasound imaging.11 Also, optical techniques have been developed to measure the pressure, e.g., by means of a Fabry–Perot film placed on the specimen12 and making use of the reflection of a sensing laser beam being modulated by the Fabry–Perot film deformation.

The second principle of detection measures the surface deformation caused by the interaction of the arriving pressure waves with the surface without contact. Carp and Venugopalan13 introduced an interferometric noncontact detection approach to measure such surface displacement, with a thin water film on the specimen surface. Rousseau et al.14 proposed an entirely optical noncontact detection method, using a confocal Fabry–Perot interferometer in differential mode. Wang et al.15 introduced a fiber-based low-coherence interferometer as a surface scanning acoustic detector. Berer et al.16 and Hochreiner et al.17 introduced a fiber-based Mach–Zehnder interferometer (MZI) using a low bandwidth CW fiber laser at a wavelength of 1550 nm. With two linear stages, the measuring spot scanned the surface of the specimen. The system was combined with a fiber-based spectral domain optical coherence tomography using wavelength-division multiplexing. Park et al.18 and Eom et al.19 introduced an all-fiber MZI, working in the heterodyne mode, as an acoustic detector. Here, an xy translation stage was used to scan the surface of the specimen. Hajireza et al.20 introduced a noninterferometric, noncontact detection method, where local changes of the refractive index due to the arriving pressure wave are measured. A drawback of these approaches for in vivo imaging is the requirement for mechanical scanning of an optical beam, which consequently leads to long acquisition times.

Our high speed noncontact holographic detection approach for photoacoustic tomography (PAT) also measures the transient surface displacement. The main advantages are as follows: The data acquisition for a full photoacoustic three-dimensional (3-D) dataset can be obtained within a fraction of a second, thus widely suppressing motion artifacts. A physical contact to the imaged object is not necessary. A high versatility exists with respect to scalability of the field-of-view (FoV) and the sampling rate, which can be varied according to the measurement requirements. These features might be an important step toward in vivo imaging of awake animals or humans and, finally, clinically usable devices.

The aim of this work is to enable imaging of subcutaneous vascular networks in tissue in situations that do not allow physical contact, e.g., clinical imaging of wounds or imaging in neurosurgery that does not allow contact approaches (see Fig. 1). The feasibility of noncontact full-field PAT using a high speed camera was shown by Horstmann et al.21 using artificial tissue phantoms containing simple, spherical absorber structures. The calculation of the surface displacement was based on spatial phase shifted electronic speckle pattern interferometry (SPS ESPI). This article focuses on an evaluation of the interference based on off-axis digital holography. In comparison, SPS ESPI has several disadvantages, which are briefly described below and motivate the change: To apply the SPS ESPI algorithm, a minimal speckle size of three pixels is necessary. The larger the speckle size, the worse the lateral resolution. Using algorithms from off-axis digital holography in combination with nonlinear filtering techniques,22 an optimization of the signal usage in the frequency domain that corresponds to a reduction in the speckle size to two pixels is achieved. Due to the Fourier transform used in the evaluation of the digital holograms, the calculation is faster and can easily be further accelerated by parallel GPU processing.

Fig. 1

Schematic illustration of the measurement principle to image subdermal vascular networks in human skin.


The aim of this paper is to demonstrate the new algorithm. Experiments on tissue models containing complex absorber structures modeling vascular networks are carried out and discussed. In addition, phase noise due to tissue motion was determined to assess the necessary frame rate for in-vivo measurements. These measurements were carried out on three test persons and a silicone tissue phantom as a reference.


Materials and Methods



For the detection of small surface displacement caused by photoacoustic pressure waves, a modified Mach–Zehnder based interferometer was realized, as schematically illustrated in Fig. 2.

Fig. 2

Modified MZI to detect surface displacement caused by photoacoustic excitation. A scalable surface area is illuminated by the detection laser. Slight variations of the surface topography induced by the pressure wave impact lead to alterations of the phase relation between object and reference waves, which change the captured interferograms.


For photoacoustic detection, the light pulses emitted from the detection laser are expanded and subsequently divided into an object and a reference beam by a polarization-sensitive beam splitter. A λ/2-wave plate is used to modify the polarization angle, which makes the object and reference intensity adjustable. The surface of the specimen is illuminated by the object beam. The backscattered and reflected light is superimposed with the reference beam, which incidents under a small angle slightly different from 0 deg to separate the different signal components (see Sec. 2.2) when imaged with a high speed camera. Due to the coherence of the illumination light and the optical roughness of the surface, the image is superimposed by a speckle pattern. Slight variations of the object’s topography caused by the pressure wave after photoacoustic excitation influence the phase of the speckle field, which is recorded as a digital hologram of the interference between the object and reference wave. To determine the phase shift, which is proportional to the surface displacement, the phase difference between two reconstructed speckle patterns is calculated.


System components

The setup for PAT consists of three main components:

Excitation laser. To generate an emission of acoustic waves under stress confinement, a Nd:YAG laser (Edgewave, BX60-2-G) with a wavelength of 1064 nm, a pulse duration of 10 ns, and a repetition rate of 1 kHz is used. Its beam profile at the specimen’s surface has a top-hat shape with an area size of 16  mm2. The radiant exposure of the excitation HEL is slightly below the maximum permissible exposure (MPE) for human skin of 20  mJ/cm2.23

Detection laser. For the interferometric detection, the surface of the specimen is illuminated with a Nd:YAG laser (Newport Spectra Phyiscs Explorer 532-2Y), with a wavelength of 532 nm, a pulse duration of 15 ns, and a Gaussian lateral beam shape with an expanded diameter of 20 mm FWHM. The radiant exposure of the detection laser HDL was only 0.2  mJ/cm2. Hence, the combination of both lasers, according to Eq. (1), did not violate the MPE for human skin:



High speed camera. A high-speed camera (Photron SA3-120k-M1) is used to record the interferograms. It has a resolution of 768×768  pixel and a pixel pitch of 17  μm.

For synchronization, all electronic devices are triggered by an I/O counter/timer module (NI PCI 6602, National Instruments Corp.). To guarantee a data acquisition time in the ms-range, the repetition rate of the detection laser and the frame rate of the camera are working at 2 kHz. The excitation laser operates with half the repetition rate of the detection laser.

When photons from the object beam incident on the specimen surface, they will be backscattered or, based on the optical properties, penetrate the tissue. Inside the tissue, they will be scattered or absorbed. Without any filtering, the camera would detect all backscattered photons limited by the numerical aperture of the imaging system, which leads to a superposition of scattering from all depths within the penetration depth. For the noncontact detection approach, only the photons backscattered from the tissue–air interface are desired, which are the majority of the incident photons. To distinguish between backscattered light from different penetration depths inside the tissue, the linear polarized object wave changes its polarization by passing through a quarter-wave plate. This results in a circular polarized object beam. The reflection at the surface leads to a reversal of the direction of rotation. Photons that are multiple scattered inside the tissue lose their polarization and will be filtered out by the beam splitter, so their intensity is reduced up to 50%. The lateral resolution of the imaging system to measure the surface deformation is 35  μm.


Phase Determination

The phase determination of the recorded interferograms is essential for the calculation of the relative displacement. We use a digital off-axis holography approach, which is based on the fringe-pattern analysis of Takeda et al.24 The phase determination method used here was first reported by Schnars.25 For optimization, the extended spatial phase filtering technique26,27 is used in the frequency domain to extract relevant information (cross-correlation part). Figure 3 illustrates the computation of the phase maps. In the first step, the recorded interferogram is Fourier transformed. Due to the limited bandwidth of the imaging system and the off-axis configuration, the information carrying cross-correlation part (highlighted by the red circle) can be separated and extracted by multiplication of the Fourier transformed interferogram F{I(m,n)} with a spectral mask W(km,kn). In the last step, the spatial-filtered spectrum will be transformed to the time domain, where the phase information ϕ(m,n) can be extracted for further postprocessing steps:



Fig. 3

Processing to determine the phase based on the captured interferograms. After applying the FFT, the cross-correlation or its conjugated complex (highlighted by the red circle) is filtered out. Subsequently, the phase can be extracted after inverse Fourier transformation (IFFT).


Equation (2) is the mathematical representation of the phase extraction, where m and n are the dimensions in the time domain and km and kn in the frequency domain.


Displacement Determination

To resolve the time-varying surface displacement, a sampling rate in the MHz range is required; however, it is beyond the capabilities of the camera. Therefore, a triple-pulse mode was developed by Horstmann et al.21

Using the triple pulse detection mode, two detection laser pulses are applied to the tissue phantom for measuring the deformation at one point in time. Between these two pulses, the photoacoustic wave is excited by a third pulse. This basic measurement sequence, which is repetitively carried out while increasing the time delay between the excitation laser and the second detection laser pulse, is schematically shown in Fig. 4. The first camera image serves as a reference image and the second image is used to detect the surface, which is deformed by the incoming photoacoustic pressure wave. By subtraction, a difference image, which maps the movements between the two images, is computed. Low-pass filters are applied for further noise reduction. This technique allows the system components to be operated in the kHz regime while recording the photoacoustic signal with a sampling rate in the MHz range.

Fig. 4

Basic sequence of the triple pulse mode for repetitive detection of surface displacement.


The result is a stack of phase difference images, which can be converted into geometrical displacement Δξ using the following Eq. (3):


with m (columns: m[1,,M]) and n (rows: n[1,,N]) defining the pixel position of the sensor array and tx describing the time of the captured images.


Three-Dimensional Tomographic Reconstruction

For a 3-D tomographic reconstruction, the determined surface displacement data are transferred into pressure data using


which is simplified by neglecting the acoustic impedance Z of the tissue.14 In the second step, an adapted fast Fourier transform (FFT)-based backprojection reconstruction algorithm8 was used to calculate the PAT image. The results are postprocessed (thresholding) and visualized by ParaView (Kitware Inc.). The reconstruction approach overcomes problems regarding accuracy and computational effort of the 3-D delay-and-sum beam forming algorithm,28 which was used for reconstruction by Horstmann et al.21


Sensitivity Analysis


Displacement determination sensitivity

The measurement of the phase shift or surface displacement is fundamentally influenced by the underlying noise level. Using a passively cooled laser and a motion isolating table, detection noise was minimized to shot noise level, basically limited by the amount of available light and the full well capacity (fwc) of the camera pixels. The fwc, which defines the amount of charge a pixel can hold before saturating, was specified by the manufacturer at 42,000e for the camera. To verify this value, the camera chip was illuminated with a super luminescence diode to generate a speckle-free brightness gradient from the minimum to the maximum detectable gray values. The aim was to determine the detector noise, as introduced in Ref. 29, and to define the noise equivalent displacement (NED). The mean and the standard deviation (or STD) of the gray value for each pixel in 1000 captured images were calculated. The total noise function is given by


where n is the number of detected photons and σSHOT describes the STD for the number of detected photons, called “photon shot noise,” and is given by σSHOT=n. σREAD is called “read noise” and includes any noise sources that are signal independent. Usually, σREAD is constant across the whole measuring range (σREAD=cread). σFPN is the fixed pattern noise and describes fluctuations in the charge collection process in the different pixels, which lead to nonrandom pixel-to-pixel sensitivity differences. Its STD is proportional to the number of photon counts σFPN=cFPN·n.

After numerical correction of each pixel, the FPN can be neglected and Eq. (5) can be reduced to



Based on the effective quantum yield ηE of the camera chip, the number of detected photons leads to n=ηEDN¯, where DN represents the digital number as a unit of the generated raw data. Due to the shot noise limitation of the system, the connection between n and DN can be rewritten as follows:



From these observations, the following equation gives the noise as a function of average digital numbers called the photon transfer function (PTF):


with cshot and cread being fit parameters, which are determined using the curve fitting toolbox of Mathworks MATLAB.


Photoacoustic sensitivity

The pressure wave p(x,y,z,t), caused by the optical excitation, propagates through the object as a transient bipolar signal. This can be associated with a particle displacement ξ(x,y,z,t) in the direction of the wavefront normal, which leads to the varying topography at the surface of the object. For a small amplitude of the pressure transient propagating in the z direction perpendicular to the wavefront at time t, pressure and displacement are related by


where E is the appropriate modulus of elasticity for the medium.30 The link to the acoustic impedance Z of the medium leads to Eq. (4). The factor of 2 considers the fact that the particle displacement at the free boundary is twice the particle displacement within the tissue.14 Here, it should be noted that the acoustic impedance depends on the actual acoustic wavelength.

Equation (4) can be used as a transfer function between the estimated NED and the noise equivalent pressure (NEP). For an analysis of the photoacoustic sensitivity, a spherical pressure source is assumed. The spherical symmetry allows the approximation of the transient duration and its signal frequency depending on the size and the speed of sound within the specimen.10 For the simulation of data, a varying displacement function ξ(z,t) with Gaussian shape and a sampling rate of 40 MHz was used; it was adapted with regard to the maximum amplitude corresponding to the NED and the time duration depending on the absorber diameter. Its simulated diameter varied between 60  μm and 3 mm, which corresponds to a signal frequency between 25 and 0.5 MHz at a sound velocity of 1500  m/s and an acoustic impedance Z1.5×106  Pas/m14 for soft tissue. To determine the NEP at the surface of the specimen, the displacement function was transferred into pressure.

The depth of the targeted vessels of the hypodermis (subcutis) and the underlying structures varies depending on the observed region of the body, age, sex, and disease status. As a first estimate, we expect a depth of the photoacoustic sources to range from 1 to 5 mm. To estimate the necessary initial pressure, the frequency depending attenuation has to be taken into account.


describes the pressure attenuation for a single-frequency spherical wave, where z is the depth coordinate, p0 is the initial pressure, and α is the attenuation coefficient (dB/cm). The attenuation coefficient depends on the material properties and the frequency31 and can be approximated by


where α0 is often 0 and y is a power law exponent, which is ranges between 0.9 and 2 for tissue.32 In the case of the skin, a linear relationship y=1 and α1=0.4[dBcmMHz] 3334.35 is assumed. Values for y reported in the literature vary between 1 and 1.17. For acoustical plane waves and a linear relationship, Deán-Ben et al.36 introduced an analytic expression for attenuation.

In the second step, the signal attenuation due to its frequency dependence was simulated for a source depth between 1 and 5 mm to estimate the effect of signal reduction on photoacoustic sensitivity.


Evaluation Experiments

To evaluate our noncontact approach for PAT, two representative tissue phantoms were produced.


Silicone tissue phantom

The first one is made of a two component, highly transparent silicone (Wacker, Elastosil RT 604 A/B). The optical properties of the absorber surrounding the silicone were varied by barium sulfate (BaSO4) to enhance the reduced scattering coefficient close to human skin tissue37 at a wavelength of 1064 nm by mixing silicone and barium sulfate at a ratio of 10 g:1.4 g. For emulation of the blood vessels, silicone tubes (Silastic BioMedical Grade ETR Elastomer Q7-4750) with an inner diameter of 310  μm (outer diameter 640  μm) are used. The silicone tubes contain a blood substitute, which is also made of silicone (Wacker, Elastosil RT 604 A), that stayed fluid due to the omission of the cross-linking agent. The absorption coefficient was increased by adding a black pigment paste. Preliminary investigations yield no quantifiable displacement with an absorption coefficient of hemoglobin (μa=3  cm1) at an excitation wavelength of 1064 nm. For this reason, the absorption was increased to μa=20  cm1, which corresponds to an excitation at 600 nm.37

Barium sulfate was chosen as a scattering component because it generates an X-ray contrast. This enables a coevaluation based on the shape and the position of the reconstructed structures in comparison to high resolution computed tomography (μCT) recordings (GE imagination at work, phoenix nanotom m, nanoCT) with a voxel resolution of 25  μm. Figure 5 shows the maximum intensity projection of the reconstructed μCT data to illustrate the position arrangement of the silicone tubes.

Fig. 5

Maximum intensity projection (MIP) of the μCT data to illustrate the position of the silicone tubes. The green and red dashed squares show the detection FoV (7×7  mm) and the excitation area (5×5  mm).


To increase the visibility in the PAT images, a threshold filter is applied to the data. The choice of an appropriate threshold is a compromise between the aims to reduce the artifacts while preserving the reconstructed structures. The reconstructed PAT and μCT data are registered to each other using MeVisLab (MeVis Medical Solutions AG) to illustrate the quality of conformity.

In consideration of the MPE limitations, an area of 5×5  mm2 was irradiated. This area was centrally arranged with respect to the detection FoV, which has a size of 7.7×7.7  mm2.


Porcine skin phantom

To create a more realistic tissue phantom in terms of mechanical and optical properties, the second phantom was made of porcine skin with embedded silicone tubes with an inner diameter of 310  μm to imitate blood vessels. A thin layer (about 2 mm) was dissected from a porcine ear and folded to embed the silicone tube in between. After the dissection, the hairs were removed. Ultrasound gel was used to decrease the impedance mismatch between the silicone tube and the surrounding skin. In this experiment, the same absorber material was used as described above. The arrangement of the silicone tubing in the phantom is schematically illustrated in Fig. 6. Here, in consideration of the MPE for skin, an area of 4×4  mm2 was irradiated.

Fig. 6

Schematic illustration of the tissue phantom made of porcine skin with an embedded silicone tube that contains liquid black silicone absorber with an absorption coefficient of 20  cm1. The arrangement of the tube is similar to a loop. The region at the apex of the loop was measured and analyzed.



Potential for In Vivo Application of the Detection System

Phase noise during in vivo imaging, which limits sensitivity of our PAT system, was investigated on three volunteers. The volunteers were informed about the risks and agreed to the experiments. The MPE was not exceeded during the experiments. The duration of the time interval between the reference image and the image of the deformed surface, which significantly determines the presence of motion artifacts, was varied, and the influence of the detection rate on detection noise (STD) was determined. For the volunteers, a skin region at the inner side of the forearm close to the wrist was selected for comparison to a silicone tissue phantom. For the measurements, the arm was placed on a table below the detection system and kept as still as possible. Ten measurements were carried out for all repetitions rates, including 2, 4, 8, 16, and 32 kHz. Increase of imaging speed significantly reduces the FoV on the camera, e.g., 128×128  pixel at 32 kHz. The measured radiant exposure of the detection laser varied between 0.2  mJ/cm2 (HMZB=1.8  mJ/cm2) and 0.04  mJ/cm2 (HMZB=0.1  mJ/cm2) due to the maximum possible pulse energy at the used repetition rate. There were no conflicts with the MPE. For the analysis of each measurement, the STD was calculated for all pixels over time. In the second step, the detection noise was calculated as the STD of the elevation of the tissue in the absence of an optoacoustic excitation. The averaged values of all subjects and the silicone phantom were analyzed with respect to the repetition rate.



To give an impression of how surface displacements look, the results for a silicone phantom with an embedded spherical absorber are presented in Fig. 7. The absorber has a diameter of 1 mm, and its center was placed 1 mm below the surface.

Fig. 7

Selection of surface displacement images at different time delays Δt: (a) Δt=0  μs, (b) Δt=1  μs, (c) Δt=1.5  μs, (d) Δt=2  μs, (e) Δt=2.5  μs, and (f) Δt=3.5  μs. White arrows indicate RoIs that were further analyzed (see Fig. 8).


A closer look at the surface displacement images reveals that, in the center of the incoming pressure wave, a localized deformation remains after the transition of the acoustic wave. Figure 7(f) shows the displacement 3.5  μs after the excitation pulse. Within three regions of interest (RoIs), the displacement is shown over time in Fig. 8. Inside the central RoI, the displacement increases until a time of 1.4  μs when the maximum of 19.5 nm is reached followed by a decrease, which does not return to the zero level within the observation period. Depending on the size of the absorbing structure, the excitation strength, and the mechanical properties of the phantom, the decrease slows down and reaches the zero level within milliseconds. In the outer region of the FoV, this effect also exists but is less pronounced.

Fig. 8

Analysis of the transient surface displacement at the three different RoIs, which are marked in Figs. 7(e) and 7(f).



Silicone Tissue Phantom

Figure 9 presents the reconstruction of the measurements of a silicone tissue phantom with vessel mimicking absorber tubes. The first three pictures show maximum intensity projections from different sides of view. The last is a thresholded and 3-D rendered illustration. In particular, in the yz- and xy-projection, the interlaced configuration of the absorbing silicone tubing is easily visible.

Fig. 9

MIPs from (a) xy-view, (b) yz-view, (c) xy-view, and (d) one 3-D rendered and thresholded view of the reconstruction results of the silicone tissue phantom.


Due to the limited aperture of the detector geometry and the smaller excitation area compared with the detection FoV, the reconstruction in the peripheral regions of the FoV image quality deteriorates.38 In addition, a decrease of the signal amplitude with increasing depth and corresponding reduction of contrast is detectable. Nevertheless, absorber structures, which are clearly distinguishable from the background, can be reconstructed up to a depth of 2.5 mm.

Figure 10 shows a selection of different depth layers of reconstructed PAT and μCT data as an overlay. The white rotated square represents the registered FoV of the PAT data.

Fig. 10

Selection of slices from varying depth of the overlaid PAT and μCT data to illustrate the conformity of the reconstructed structures. The white square emphasizes the registered FoV of the PAT data. (a) Depth=0.9  mm, (b) depth=1.5  mm, (c) depth=1.7  mm, (d) depth=2.0  mm, (e) depth=2.1  mm, and (f) depth=2.3  mm.



Porcine Skin Phantom

Figure 11(a) illustrates the measured displacement 1.5  μs after photoacoustic excitation of the porcine skin tissue phantom. The displacement data were temporally smoothed, and a histogram spread was performed in the presented images to increase the contrast.

Fig. 11

Analysis of the measured displacement data of the porcine skin tissue phantom. (a) Surface displacement 1.5  μs after photoacoustic excitation. The image was temporally filtered for smoothing the signal, and a histogram spread was performed to increase the contrast. (b) Displacement data extracted from the RoI, which is located in a region where the displacement reaches a maximum value [see (a) red box].


The red square marks a RoI, which is located in a region where the displacement reaches its maximum value. The time course of the displacement within this RoI is shown in Fig. 11(b). The ratio of the maximum value (green dashed line) to the new zero level (blue dashed line) after the decay of the surface deformation shows the influence of the mechanical properties of the specimen. Here, the maximum measured displacement was 7.1 nm and the remaining background level after 30  μs was 1.6 nm, which corresponds to 24%. For comparison, in the silicone-based phantom introduced at the beginning of this section (see Fig. 8), the maximum measured displacement decreased from 19.5 to 10.7 nm, which corresponds to a background of 55%.

Figure 12 shows the results of the reconstruction of the PAT measurements. The first three figures show the maximum intensity projection from different directions. In addition, a 3-D visualization of binary segmented data (light gray) is shown. Above the absorber close to the surface of the skin, a second structure is visible. A third structure is located at the bottom of the volume.

Fig. 12

MIPs of (a) xz-view, (b) yz-view, (c) xy-view, and (d) the 3-D reconstruction of the blood vessel imitating silicone tube that was arranged as a loop in between the folded porcine skin layer. The silicone tube has an inner diameter of 310  μm and contains the liquid absorber material, which has an absorption coefficient of μa=20  cm1.



Displacement Determination Sensitivity

Figure 13 shows signal mean and STDs for each pixel (black dots) of the camera, when exposed for a certain period of time to a light gradient. The red solid line shows a PTF curve fitting, which yields the parameters, according to Eq. (8): cshot=6 and cread=44. From the fit parameter, the full well capacity is calculated as



Fig. 13

Logarithmic representation of the mean and STD values (black dots) of the individual pixels to which the PTF (red line) is fitted.


The value is smaller by a factor of 1.7 than the value specified by the manufacturer.

The shot noise at maximum exposure is fwc=158e, which corresponds to 0.6% of the dynamic range. The intensities of a captured interferogram are represented as a 12-bit (4096 counts) number. This results in a statistical fluctuation of 24.6 counts. Due to error propagation, this fluctuation of intensity values can induce a deviation in the displacement measurement sensitivity in the range of


σ(ξ)=λ2πfwc=532  nm2π250000.5  nm.

A measurement STD of ±1  nm was determined by analyzing measurements without optical excitation. For further noise reduction, a temporal low-pass filter, which sets the frequencies above a certain threshold frequency to zero, is applied. This reduces the noise close to 0.5 nm.

The ratio of the measured displacement to the noise floor is defined here as the signal-to-noise ratio (SNR). For a 1-mm spherical absorber, which is introduced at the beginning of Sec. 3, the maximum geometrical out-of-plane deformation was around 20 nm (see Fig. 8). Assuming a noise level of ±1  nm, the resulting SNR is about 20.


Photoacoustic Sensitivity

Figure 14 shows the simulation to estimate the photoacoustic sensitivity based on the relation between pressure and surface displacement that is given by Eq. (4). The frequency-dependent maximum pressure in kPa at the surface of the specimen to create a 0.5-nm displacement ξ (NED) is shown in Fig. 14(a). Due to the reciprocal relation between absorber diameter and pressure, the course follows a 1d-decrease. The two solid lines in Fig. 14(b) show the simulation of the frequency-dependent transmission of the acoustic wave (attenuation) for a spherical photoacoustic source located in 1 and 5 mm depths. Attenuation increases with both signal frequency and depth. The combination of the frequency-dependent attenuation and the absorber diameter-dependent NEP is shown for five absorber depth positions in Fig. 14(c). For a better visibility, the diagram was limited to a pressure range of up to a maximum 100 kPa. The simulation demonstrates the difficulty of detecting structures that decrease in size and increase in depth. In other words, small objects lying deeper are difficult to detect.

Fig. 14

NEP as a function of the signal frequency. (a) Simulation based on the relation between pressure and surface displacement [Eq. (4)] and the determined surface displacement noise (0.5 nm) introduced in the previous section in a frequency range of 0.5 to 25 MHz. (b) Simulation of the frequency depending attenuation of the pressure signal presented for 1 and 5 mm depth positions of the absorber. (c) NEP as a combination of (a) and (b) for five different absorber depth positions.



Potential for In Vivo Application of the Detection System

Figure 15 shows the results of the detection noise depending on the repetition rate. The silicone phantom (blue dots and blue dotted line) as a reference shows a mean noise of 0.4 nm with a STD of ±0.5  nm, which increases slightly with increasing repetition rates. The results of the in vivo measurements show a significantly higher noise using a lower repetition rate (see Fig. 15, green dots and green dotted line), which decreases for each measurement subject to less than 1 nm at 32 kHz. In addition, the STD decreases, which indicates an enhancement of the measurement stability. A slight increase of noise between 16 and 32 kHz is also visible in the in vivo measurements.

Fig. 15

Noise depending on the repetition rate for three in vivo measurements (green dot and green dotted line), in which the noise could be reduced to a level of less than 1 nm by increasing the repetition rate. Blue dots and blue dotted line: silicone phantom measurement as reference with an almost constant noise of 0.4 nm.




Algorithms from off-axis digital holography to calculate surface displacements, which are caused by the interaction of the arriving pressure transients with the tissue surface, were successfully evaluated.


Evaluation Experiments

Maximum intensity projections of reconstructed volumes measured with the silicone tissue phantom are presented in Fig. 9. The three vertically straight as well as the three horizontal twisted tubes can be identified well in the x/y-view representation. To improve the visibility of the structures, segmentation is used in the postprocessing step. The decrease in the signal amplitude with increasing depth, which is expected from the simulations shown in Fig. 14, complicates the choice of the threshold value. Nevertheless, structures can clearly be distinguished from the background up to a depth of 2.5 mm [see Figs. 9(a) and 9(b)].

Comparison with μCT data shows the geometrical consistency of the reconstructed PAT structures. With increasing imaging depth, the signal leaks within the PAT data increase due to the detection sensitivity and signal shadowing.

As a more realistic phantom for vascular structures in skin, an artificial absorber was embedded in porcine skin. In the central part of the reconstructed volume, the loop-like run of the absorber structure is very clearly delineated from the background. It matches the tubes in shape and dimensions. Due to the natural pigmentation (cutaneous melanin) of the porcine skin, which depends on the volume fraction and the wavelength depending absorption coefficient of cutaneous melanosome,37 a reconstructible photoacoustic signal is expected from these structures. The thin laterally spread structure above the constructed vessel structure can be assigned to the mentioned melanin layer. The flat structure at the bottom of the volume could be a reflection artifact of the ultrasonic wave of the lower part of the folded skin flap. The reconstruction results of both phantom types visualized as xz-view and yz-view MIP show periodic artifactual structures. The analysis to determine the source revealed an influence of the cutoff frequency of the low pass filter on the shape and the periodicity of the artifacts but without suppressing them. At the moment, we are evaluating new filter techniques to reduce the artifacts.

The influence of the different viscoelastic properties between silicone and skin becomes visible by comparing Figs. 8 and 11(b). After the surface displacement reached its maximum, the deformation does not decrease to zero but stays until the end of the measurement at a level of 25% to 50% of the maximum signal. In previous studies on silicone phantoms, the decrease to the considered zero level took 1 ms. It is assumed that the elevated background is caused by a viscoelastic surface deformation superimposed on the photoacoustic signal. It is expected that the viscoelastic deformation depends on the mechanical material properties since different levels were observed between silicone and porcine skin. To exclude a thermal expansion due to heating by the repetitive excitation, measurements were carried out in the reversed scanning mode. Here, the measurement started with a delay, corresponding to a maximum depth, and continued to a smaller delay. If thermal expansion drives the background, an increase signal above the surface is expected. However, these measurements showed the same characteristics as in the forward scanning mode. Heating of the sample during measurement has no effect on the PAT signal. Nevertheless, the viscoelastic deformation has no influence on the results of the reconstruction since the temporal derivative of the deformation is used to calculate the initial pressure field.


Motion Artifacts

The triple pulse measurement scheme uses a temporal separation of the two camera frames to determine the phase difference before and after excitation. This unnecessarily long time contributes significantly to motion artifacts. The results of the noise analysis of the in vivo measurements demonstrate the necessity of an enhancement of the repetition rate in the current detection scheme to reduce the noise to carry out in vivo imaging. The enhancement leads to a decrease of the used camera resolution to 128×128  pixel at 32 kHz. This small detection area is not sufficient for imaging. Either faster cameras with repetition rates 32  kHz or different detection techniques are needed to overcome this problem. Imaging both holograms, before and after excitation, in one frame by angular multiplexing techniques for phase calculation39,40 reduces the camera speed to the excitation frequency.


Photoacoustic Sensitivity

The simulation results in Fig. 14 show the NEP at the surface of the specimen using the relation between pressure and surface displacement that is given by Eq. (4) as well as the signal attenuation, both depending on the signal frequency. The combination of varying efficiencies of the sound generation, detection, and attenuation in the tissue defines the necessary initial pressure to detect an absorbing volume. For example, a spherical absorber 1 mm in diameter results in a frequency of the pressure transient of about 1.5 MHz and a NEP of 2.4 kPa at the surface of the specimen. At a depth position of, for example, 2 mm, an initial pressure of 16  kPa would be necessary. Knowing the optical properties of the specimen and the absorber and the excitation wavelength allows an estimation of the minimal excitation to generate a detectable signal. Eventually, sensitivity of the detection and the restrictions of the MPE limit the imaging depth range and the minimum absorber size that can be detected. Further investigation should determine experimentally the minimum absorber size that can be detected and compare this with direct acoustic detection of the PAT signals.



In conclusion, algorithms from digital holography have been successfully used to calculate surface deformation in a noncontact PAT detection approach. The general feasibility of noncontact area detection of the photoacoustic wave was shown by Horstmann et al.21 using tissue phantoms with an embedded spherical absorber. Here, we successfully imaged tissue phantoms with an embedded complex absorber structure, which imitates vascular networks. The 3-D reconstruction was significantly improved by a FFT-based backprojection technique. Comparison of μCT recordings as a ground truth showed a high morphological matching without any visible distortions. Furthermore, absorbing structures with a diameter of 320  μm could be reconstructed up to a depth of 2.5 mm. The noise analysis to verify the feasibility for in vivo imaging shows the necessity to either increase the repetition rate to 32 kHz or to introduce a multiplexing scheme for the detection of the holograms. For clinical use, angular multiplexing seems to be the more promising approach since it will work with frame rates of 100 Hz, which is used in machine vision and even smart phone cameras. Currently, experiments are under way to prove this concept.


J. H. and R. B. are named as inventors of a patent related to the method of full-field noncontact photoacoustic imaging, which is owned by Medizinisches Laserzentrum Lübeck GmbH. R. B. is the CEO of Medizinisches Laserzentrum Lübeck GmbH.


This work was funded by the German Federal Ministry of Education and Research in the project OptoAk and LUMEN (FKZ 13EZ1140A/B, FKZ 13N12533). LUMEN is a joint research project of Lübeck University of Applied Sciences and Universität zu Lübeck and represents a branch of the Graduate School for Computing in Medicine and Life Sciences of Universität zu Lübeck. We also want to thank Prof. Dr. Botterweck (Lübeck University of Applied Sciences) for the use of the μCT system.


1. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus 1, 602–631 (2011). http://dx.doi.org/10.1098/rsfs.2011.0028 Google Scholar

2. L. V. Wang and L. Gao, “Photoacoustic microscopy and computed tomography: from bench to bedside,” Ann. Rev. Biomed. Eng. 16, 155–185 (2014).ARBEF71523-9829 http://dx.doi.org/10.1146/annurev-bioeng-071813-104553 Google Scholar

3. S. Hu and L. V. Wang, “Photoacoustic imaging and characterization of the microvasculature,” J. Biomed. Opt. 15(1), 011101 (2010).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.3281673 Google Scholar

4. M. Heijblom et al., “Imaging tumor vascularization for detection and diagnosis of breast cancer,” Technol. Cancer Res. Treat. 10(6), 607–623 (2011). http://dx.doi.org/10.7785/tcrt.2012.500227 Google Scholar

5. J. Yao et al., “Label-free oxygen-metabolic photoacoustic microscopy in vivo,” J. Biomed. Opt. 16(7), 076003 (2011).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.3594786 Google Scholar

6. Y. Xu and L. V. Wang, “Time reversal and its application to tomography with diffracting sources,” Phys. Rev. Lett. 92, 033902 (2004).PRLTAO0031-9007 http://dx.doi.org/10.1103/PhysRevLett.92.033902 Google Scholar

7. B. E. Treeby, E. Z. Zhang and B. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Probl. 26(11), 115003 (2010).INPEEY0266-5611 http://dx.doi.org/10.1088/0266-5611/26/11/115003 Google Scholar

8. K. P. Köstli et al., “Temporal backward projection of optoacoustic pressure transients using Fourier transform methods,” Phys. Med. Biol. 46(7), 1863–1872 (2001).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/46/7/309 Google Scholar

9. C. G. Hoelen and F. F. de Mul, “Image reconstruction for photoacoustic scanning of tissue structures,” Appl. Opt. 39, 5872–5883 (2000).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.39.005872 Google Scholar

10. J. Xia, J. Yao and L. V. Wang, “Photoacoustic tomography: principles and advances,” Prog. Electromagn. Res. 147, 1–22 (2014). http://dx.doi.org/10.2528/PIER14032303 Google Scholar

11. K. Daoudi et al., “Handheld probe integrating laser diode and ultrasound transducer array for ultrasound/photoacoustic dual modality imaging,” Opt. Express 22, 26365–26374 (2014).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.22.026365 Google Scholar

12. J. Laufer et al., “Three-dimensional noninvasive imaging of the vasculature in the mouse brain using a high resolution photoacoustic scanner,” Appl. Opt. 48(10), D299–D306 (2009).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.48.00D299 Google Scholar

13. S. A. Carp and V. Venugopalan, “Optoacoustic imaging based on the interferometric measurement of surface displacement,” J. Biomed. Opt. 12(6), 064001 (2007).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2812665 Google Scholar

14. G. Rousseau et al., “Non-contact biomedical photoacoustic and ultrasound imaging,” J. Biomed. Opt. 17(6), 061217 (2012).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.17.6.061217 Google Scholar

15. Y. Wang, C. Li and R. K. Wang, “Noncontact photoacoustic imaging achieved by using a low-coherence interferometer as the acoustic detector,” Opt. Lett. 36, 3975–3977 (2011).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.36.003975 Google Scholar

16. T. Berer et al., “Multimodal noncontact photoacoustic and optical coherence tomography imaging using wavelength-division multiplexing,” J. Biomed. Opt. 20(4), 046013 (2015).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.20.4.046013 Google Scholar

17. A. Hochreiner et al., “Non-contact photoacoustic imaging using a fiber based interferometer with optical amplification,” Biomed. Opt. Express 4, 2322–2331 (2013).BOEICL2156-7085 http://dx.doi.org/10.1364/BOE.4.002322 Google Scholar

18. S. J. Park et al., “Noncontact photoacoustic imaging based on all-fiber heterodyne interferometer,” Opt. Lett. 39, 4903–4906 (2014).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.39.004903 Google Scholar

19. J. Eom, S. J. Park and B. H. Lee, “Noncontact photoacoustic tomography of in vivo chicken chorioallantoic membrane based on all-fiber heterodyne interferometry,” J. Biomed. Opt. 20(10), 106007 (2015).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.20.10.106007 Google Scholar

20. P. Hajireza et al., “Non-interferometric photoacoustic remote sensing microscopy,” Light Sci. Appl. 6(6), e16278 (2017). http://dx.doi.org/10.1038/lsa.2016.278 Google Scholar

21. J. Horstmann et al., “Full-field speckle interferometry for non-contact photoacoustic tomography,” Phys. Med. Biol. 60(10), 4045–4058 (2015).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/60/10/4045 Google Scholar

22. C. S. Seelamantula et al., “Zero-order-free image reconstruction in digital holographic microscopy,” in IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro, pp. 201–204, IEEE (2009). http://dx.doi.org/10.1109/ISBI.2009.5193018 Google Scholar

23. A. N. Standard and Laser Institute of America, ANSI Z136. 1 Safe Use of Lasers (2014), Laser Institute of America (2014). Google Scholar

24. M. Takeda, H. Ina and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72, 156–160 (1982).JOSAAH0030-3941 http://dx.doi.org/10.1364/JOSA.72.000156 Google Scholar

25. U. Schnars, “Direct phase determination in hologram interferometry with use of digitally recorded holograms,” J. Opt. Soc. Am. A 11, 2011–2015 (1994).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.11.002011 Google Scholar

26. H. O. Saldner, N.-E. Molin and K. A. Stetson, “Fourier-transform evaluation of phase data in spatially phase-biased TV holograms,” Appl. Opt. 35, 332–336 (1996).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.35.000332 Google Scholar

27. E. Cuche, P. Marquet and C. Depeursinge, “Spatial filtering for zero-order and twin-image elimination in digital off-axis holography,” Appl. Opt. 39, 4070–4075 (2000).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.39.004070 Google Scholar

28. T. Berer et al., “Reconstruction algorithms for remote photoacoustic imaging,” in IEEE Int. Ultrasonics Symp. (IUS), pp. 1–4, IEEE (2012). http://dx.doi.org/10.1109/ULTSYM.2012.0579 Google Scholar

29. J. R. Janesick, Photon Transfer DN to [lambda], Vol. 170, SPIE Press, Bellingham, Washington (2007). Google Scholar

30. R. W. Speirs and A. I. Bishop, “Photoacoustic tomography using a Michelson interferometer with quadrature phase detection,” Appl. Phys. Lett. 103(5), 053501 (2013).APPLAB0003-6951 http://dx.doi.org/10.1063/1.4816427 Google Scholar

31. T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out, 2nd ed., Academic Press, New York (2013). Google Scholar

32. F. A. Duck, Physical Properties of Tissues: A Comprehensive Reference Book, Academic Press, London (2013). Google Scholar

33. B. I. Raju and M. A. Srinivasan, “High-frequency ultrasonic attenuation and backscatter coefficients of in vivo normal human dermis and subcutaneous fat,” Ultrasound Med. Biol. 27(11), 1543–1556 (2001).USMBA30301-5629 http://dx.doi.org/10.1016/S0301-5629(01)00456-2 Google Scholar

34. C. Guittet et al., “High-frequency estimation of the ultrasonic attenuation coefficient slope obtained in human skin: simulation and in vivo results,” Ultrasound Med. Biol. 25(3), 421–429 (1999).USMBA30301-5629 http://dx.doi.org/10.1016/S0301-5629(98)00176-8 Google Scholar

35. C. M. Moran, N. L. Bush and J. C. Bamber, “Ultrasonic propagation properties of excised human skin,” Ultrasound Med. Biol. 21(9), 1177–1190 (1995).USMBA30301-5629 http://dx.doi.org/10.1016/0301-5629(95)00049-6 Google Scholar

36. X. L. Deán-Ben, D. Razansky and V. Ntziachristos, “The effects of acoustic attenuation in optoacoustic signals,” Phys. Med. Biol. 56(18), 6129–6148 (2011).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/56/18/021 Google Scholar

37. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/58/11/R37 Google Scholar

38. B. Cox, S. Arridge and P. Beard, “Photoacoustic tomography with a limited-aperture planare sensor and a reverberant cavity,” Inverse Probl. 23(6), S95–S112 (2007).INPEEY0266-5611 http://dx.doi.org/10.1088/0266-5611/23/6/S08 Google Scholar

39. P. Girshovitz and N. T. Shaked, “Fast phase processing in off-axis holography using multiplexing with complex encoding and live-cell fluctuation map calculation in real-time,” Opt. Express 23, 8773–8787 (2015).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.23.008773 Google Scholar

40. P. Wang, “Improved angular multiplexing method with polarization interference for enhancing the resolution of digital holographic image,” Optik 131, 312–316 (2017).OTIKAJ0030-4026 http://dx.doi.org/10.1016/j.ijleo.2016.11.100 Google Scholar


Christian Buj is a PhD student at the Universität zu Lübeck. He received his diploma degree in computer science from Carl von Ossietzky University of Oldenburg in 2009 and his MSc degree in medical engineering science from the Universität zu Lübeck in 2013. His current research interests include biomedical imaging, photoacoustic tomography, interferometry, optical coherence tomography, holography, and tissue optics.

Michael Münter is a PhD candidate at the Medical Laser Center in Lübeck. He received his MSc degree in medical engineering science from the Universität zu Lübeck, Lübeck, Germany, in 2015. His current research interests include biomedical imaging, noncontact photoacoustic imaging, optical coherence tomography, image reconstruction, data processing, and visualization.

Benedikt Schmarbeck is an engineer at the Medical Laser Center in Lübeck. He holds a BSc degree in medical engineering science from the Universität zu Lübeck since 2013. In 2015, he received his MSc degree in medical engineering science from the Universität zu Lübeck. His current fields of work are biomedical imaging, interferometry, photoacoustic measurement techniques, retinal photocoagulation, and tissue optics.

Jens Horstmann is a postdoc in the Ocular Surface Group, Department of Ophthalmology, University of Cologne, Germany. His research is concentrated on experimental imaging and clinical translation. He received his PhD from the University of Lübeck, Germany.

Gereon Hüttmann graduated in physics and physical chemistry from the University of Göttingen, Germany. From 1992 to 2005, he worked at the Medical Laser Center Lübeck. Since 2005, he leads a research group at the Institute of Biomedical Optics, University of Lübeck. He is a member of the German Center for Lung Research and vice director of the Institute of Biomedical Optics in Lübeck. His main research interests are optical coherence tomography and intravital microscopy.

Ralf Brinkmann studied physics at the University of Hannover, Germany, with a focus on quantum optics and lasers. After a five year industrial interim period, he joined the Medical Laser Center in Lübeck, Germany, in 1993 and received his PhD at the University of Lübeck. Since 2005, he has held a permanent position as a faculty member at the University Institute of Biomedical Optics and leads the Medical Laser Center Lübeck as CEO.

© 2017 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2017/$25.00 © 2017 SPIE
Christian Buj, Michael Münter, Benedikt Schmarbeck, Jens Horstmann, Gereon Hüttmann, and Ralf Brinkmann "Noncontact holographic detection for photoacoustic tomography," Journal of Biomedical Optics 22(10), 106007 (13 October 2017). https://doi.org/10.1117/1.JBO.22.10.106007
Received: 31 May 2017; Accepted: 22 September 2017; Published: 13 October 2017

Back to Top