Translator Disclaimer
22 May 2014 Deblurring algorithms accounting for the finite detector size in photoacoustic tomography
Author Affiliations +
Abstract
Most reconstruction algorithms for photoacoustic tomography, like back projection or time reversal, work ideally for point-like detectors. For real detectors, which integrate the pressure over their finite size, images reconstructed by these algorithms show some blurring. Iterative reconstruction algorithms using an imaging matrix can take the finite size of real detectors directly into account, but the numerical effort is significantly higher compared to the use of direct algorithms. For spherical or cylindrical detection surfaces, the blurring caused by a finite detector size is proportional to the distance from the rotation center (spin blur) and is equal to the detector size at the detection surface. In this work, we apply deconvolution algorithms to reduce this type of blurring on simulated and on experimental data. Two particular deconvolution methods are compared, which both utilize the fact that a representation of the blurred image in polar coordinates decouples pixels at different radii from the rotation center. Experimental data have been obtained with a flat, rectangular piezoelectric detector measuring signals around a plastisol cylinder containing various small photoacoustic sources with variable distance from the center. Both simulated and experimental results demonstrate a nearly complete elimination of spin blur.

1.

Introduction

In photoacoustic (also known as thermoacoustic or optoacoustic) imaging, a sample is illuminated by a short pulse of electromagnetic radiation, such as a laser or a microwave pulse. The electromagnetic radiation is absorbed to different extents in different regions of the sample, and the resulting weak heating and thermal expansion launch an acoustic wave. This acoustic pressure signal is measured by detectors outside the sample and can be used for reconstruction of the distribution of the absorbed electromagnetic energy inside the sample. Therefore, photoacoustic imaging is a technique that combines the advantages of optical methods (good optical absorption contrast between different types of tissues) and ultrasound (high spatial resolution due to low ultrasound scattering). Utilizing acoustic detectors with large aperture that detect signals simultaneously from all areas of the sample requires tomographic methods for image reconstruction. This is known as photoacoustic tomography (good overviews are provided in Refs. 1, 2, and 3).

Mathematically, photoacoustic wave propagation can be modeled as an initial value problem for the three-dimensional (3-D) wave equation. Within this mathematical framework, image reconstruction algorithms have been developed for detectors so small that they can be modeled as points. Physical implementations of detector arrays use piezoelectric sensors, which have the drawback of decreasing signal-to-noise ratio with decreasing size. Therefore, piezoelectric detector elements are used with a size in the range of several millimeters, leading to some blurring artifacts if reconstruction algorithms for point sensors are applied to measured data. An alternative is to use linear or planar detectors that integrate the acoustic pressure over one or two spatial dimensions. Because of pressure integration over a finite area with known shape, noise is reduced and signal-to-noise ratio is enhanced. Reconstruction algorithms have been adapted to handle such special kinds of detectors. For integrating line detectors, the inversion of the 3-D wave equation is usually split into the following two steps:4,5 first, for a given orientation of the line detectors, a two-dimensional (2-D) wave equation is inverted to recover linear projections of the initial pressure distribution. Second, multiple projections in different directions are used to reconstruct the original 3-D pressure distribution via a 2-D inverse Radon transform.

In this work, we concentrate on the first step of this procedure, i.e., the inversion of the 2-D wave equation. However, also in these modified algorithms, the finite size in the remaining dimensions is usually neglected and artifacts arise in reconstructed images. Our aim is to explain the particular form of these artifacts and to study additional signal processing algorithms that can be applied at the image postprocessing stage to reduce these artifacts. In Ref. 6, two deconvolution methods were proposed for planar measurement geometry: Wiener deconvolution and piecewise polynomial truncated singular value decomposition (PP-TSVD) deconvolution (see Sec. 3). The latter method was actually implemented. For cylindrical geometry, it was mentioned7 that the same methods could be used, but it has not been demonstrated on numerical or experimental data. One of the aims in this paper is to test the assertion in Ref. 7 and to compare the two methods in this context. The implementation as well as the results are different from what was reported for the planar case, which justifies our careful re-examination of these methods for cylindrical geometry.

Matrix-model based reconstruction algorithms8 can build arbitrary shapes and spatial impulse responses of the detectors directly into the model and are well suited to reduce finite detector artifacts, but require high numerical effort for the large matrices, even the application of graphics processing units for parallel computing.9 In our work, we show that these artifacts can be effectively reduced in a matter of a second of additional computation time, at least for some regular geometrical detector shapes, after using standard reconstruction methods like filtered backprojection10 or time reversal11 that require, in the 2-D case, about a minute on a standard PC.

After a short review of photoacoustic tomography, in Sec. 2 we address some issues that arise for image reconstruction if finite-sized detectors are used. Furthermore, we describe the theory behind deconvolution algorithms for compensation of spin-blur artifacts that appear in connection with detectors that are parts of cylindrical surfaces with finite angular aperture, or parts of their approximating tangential planes if the aperture is not too large. The accuracy of this approximation is also discussed in this section.

In Sec. 3, computer-simulation studies are performed to test the two proposed deconvolution algorithms on simulated data from different source profiles located at various distances from imaging center and recorded with detectors of various apertures. Section 4 describes real measurements on a cylindrical plastisol phantom with six thin holes filled with OrangeG as an absorbing liquid. The holes were located on a spiral emanating from the center of the cylinder. The acoustic signals were recorded by a piezoelectric detector [Polyvinylidenfluorid (PVDF) foil] having the form of a narrow long strip. The time reversal reconstruction shows the expected blurring, and we demonstrate that the blurring can be almost completely eliminated by applying the proposed deblurring algorithms. A preliminary version of this work appeared in Ref. 12.

2.

Background

2.1.

Photoacoustic Tomography with Finite-Sized Detectors

In photoacoustic tomography with line detectors, ultrasound wave propagation of the pressure integrated over the z-direction, p(x,y,t), is described by an inhomogeneous 2-D wave equation with a very short initial forcing (see, for example, Ref. 4).

2pt2c2Δp=p0(r)·ddtδ(t).

Here c is the sound speed, p0(r) is the (integrated) initial pressure produced by thermoelastic expansion, and δ(t) is the Dirac delta function. This is equivalent to an initial value problem for the homogeneous 2-D wave equation. 5

2pt2c2Δp=0p(r,0)=p0(r)tp(r,0)=0.

The goal in photoacoustic tomography is to recover the spatial distribution of absorbed energy density inside the sample, which is proportional to p0(r), from acoustic pressure signals p(rS,t) measured outside the sample at the surface S using a detector D (see Fig. 1). Ideal line detectors (in real 3-D space) correspond to ideal point detectors for the 2-D wave equation.

Fig. 1

Illustration of photoacoustic tomography. D represents a detector at position rS scanning across a closed surface S enclosing the initial source p0(r). R is the volume enclosed by S.

JBO_19_5_056011_f001.png

A simple and efficient way to solve the 2-D photoacoustic inverse problem (or analogously the 3-D problem) for point detectors and to recover the initial pressure distribution is a time reversal algorithm:11 assume a time T0>0 to be large enough that the pressure satisfies p(r,T0)0 within the whole volume R enclosed by the surface S. As the initial pressure distribution outside the surface S is zero, for the 3-D problem, T0 can be chosen, e.g., as the maximal diameter of R, divided by c. A 2-D wave propagation does not satisfy a strict form of the Huygens principle, which implies that the solution of the wave equation never exactly vanishes after a certain time, but exhibits an infinite, algebraically decaying tail. However, using the Abel transform and its inverse, one can see that the long tail in the 2-D measured signal contains no further information about the initial pressure distribution (see Ref. 12). So it is justified also in 2-D to choose such a time T0 when all pressure transients above a small level of size ε (e.g., 1/1000 of the detected pressure maximum) have passed all detection points. At that time T0, the time evolution of the pressure field is reversed; we start to “rewind the film.” Mathematically, this can be accomplished by setting the pressure values on the surface S equal to p(rS,T0t), the measured pressure values in time-reversed order. A finite difference or Fourier transform scheme for the wave equation is executed on a spatial computational grid with these time-reversed values as boundary conditions. After another time period T0, we obtain an approximation to the initial pressure p0(r). We have recently proposed a pseudospectral variant of time reversal that yields highly accurate results in 2-D and 3-D.13 This combines an exact time discretization with a Fourier (k-space) approximation to the spatial derivatives [pn(k)=p(k,nΔt)].

Eq. (1)

pn+1(k)=[24sin2(ckΔt/2)]·pn(k)pn1(k).

The time needed for the Fourier transforms (enabling high accuracy in the spatial derivatives) is offset by the fact that the Courant-Friedrichs-Lewy condition does not apply here and, thus, larger timesteps can be taken.

In this work, we will assume the more realistic scenario that instead of point-wise pressure data we have pressure values integrated over a spatially extended aperture. In Fig. 2, we illustrate our strategy of accounting for such a finite detector size in the simple case of a planar detector array: according to source-receiver reciprocity for acoustic waves, measuring signals from a point-like source with a planar detector array is equivalent to measuring signals from an extended source with an array of point-like detectors. Hence, we may (almost) recover the imaging result from a point-like detector array by deconvolution of the reconstructed image that was obtained by measuring with an array of finite-size detectors.

Fig. 2

For a planar integrating detector array, by source receiver reciprocity, we can use the well-known reconstruction formulas from point-like detectors followed by a spatial deconvolution to approximately reconstruct the initial pressure distribution, which is a point source in this figure.

JBO_19_5_056011_f002.png

A similar compensation strategy can also be applied for cylindrical or spherical detection geometry: using the translational and rotational invariance of the wave equation, Haltmeier and Zangerl show in Ref. 14 that the reconstructed image becomes blurred with a certain blurring kernel if a finite-sized detector is used instead of a point (3-D) or line (2-D) detector. This has been proven in 2-D for a detector extending to a cylinder of finite radius, which applies to optical detection by a laser beam15 or an elongated optic fiber sensor.16 In our case, a piezoelectric detector is used that extends tangentially to a finite strip (Sec. 4.1). The theory of Ref. 14 can also be adapted to such a situation as we show in the Appendix.

A different kind of finite-size issue is raised by the fact that infinite line detectors in reality can only be approximated by finite lines. In our experiment (Sec. 4.1), we have photoacoustic sources that are located in an xy-plane, detected at a radial distance of rS=16mm from the coordinate origin. The detectors are rectangular strips whose long axis points in the z-direction. To illustrate the influence of the z-extent of the detectors in this setup, Fig. 3 shows simulated signals from a small ball source with 1-mm diameter, measured by an integrating strip detector growing in z-length. For the 3×3-mm detector, one gets an N-shaped 3-D signal similar to a signal observed by a point detector. For the 10-mm-long detector, the signal initially behaves like the 2-D signal, but shows a clear deviation from the 2-D signal after 11 μs. For the 3×25-mm detector, one gets the 2-D pressure signal up to a measurement time of 13 μs; then a small negative peak appears. This peak, which is due to the arrival of the wave at the ends of the strip, could be shifted to longer times by taking an even longer strip; but it is apparent that the signal and, hence, the reconstruction at this z-length, are already sufficiently close to the ideal 2-D situation. A more detailed discussion of the required length of such integrating detectors can be found in Ref. 17. This simulation also demonstrates the gain in signal amplitude by a factor of 1.5 with increasing integration length of the detector.

Fig. 3

Simulated signals from a small ball with 1-mm diameter as a source measured at a distance of 16 mm with an integrating detector of constant x-width (3 mm) and growing in z-length: 3×3mm (solid), 3×10mm (dash), 3×25mm (dash-dot); for comparison infinite strip detector signal (dotted); the gain in positive peak amplitude by a factor of 1.5 mainly occurs between z-lengths of 3 and 6 mm.

JBO_19_5_056011_f003.png

2.2.

Algorithms for Compensation of Spin Blur

As described in the previous subsection, a narrow strip detector of sufficient length measures signals governed by the 2-D wave equation. By rotating the detector along a curve, image reconstruction can, therefore, be performed by inversion of the 2-D wave equation. Due to the width of the strip detectors, application of reconstruction methods designed for ideal line detectors will result in a blurred reconstruction.

Suppose that instead of the pointwise pressure values p(rS,t), integrals pα(rS,t) of the acoustic pressure over an arc centered at points rS on a circular recording curve are measured. Here the superscript indicates that the detector covers an angle 2α on the recording surface (compare Fig. 5). Then applying reconstruction methods, such as time reversal or backprojection, the integrated pressure values pα(rS,t) yield as reconstruction (see the Appendix).

p0α(r,φ)=ααp0(r,φ+α)dα=02πp0(r,φ)kα(φφ)dφ.

Hence, the obtained reconstruction p0α(r,φ) [expressed in polar coordinates r=(r,φ)] is a blurred version of the unblurred initial pressure distribution p0(r,φ), obtained by convolving p0(r,φ) with a purely angular blurring kernel given by the box function

kα(φφ)={1if|φφ|α0if|φφ|>α.

We emphasize that from a mathematical point of view, there is no restriction of our proposed deconvolution methods to a box function as blurring kernel. The algorithms can be applied with minor modifications to a more general kernel, i.e., a general signal impulse response function describing a variable sensitivity within the detector foil.

Since the spatial extent of a circular arc with aperture 2α is proportional to its radius, the blurring effect increases linearly for increasing distance to the center of the reconstruction domain and is equal to the detector width close to detection circle. We refer to the purely angular blurring proportional to the distance of the center of rotation as spin blurring because the same type of blurring arises in photographic imaging when a camera is spinning around the lens axis during exposure. Figure 4 shows an initial pressure distribution consisting of four points (smoothed by box function convolution) located on a spiral emanating from the detection center (left), and the time reversal reconstruction using signals from an arc-shaped detector with 2α=45deg in Cartesian (center) and polar (right) coordinates.

Fig. 4

Initial pressure (a) and two-dimensional (2-D) time reversal reconstruction in Cartesian (b) and polar (c) coordinates for a circular integrating detector of aperture 2α=45deg (bold white arc), which is rotated in 80 equidistant positions around the sample (white dots). The best resolution is near the rotation center and blurring increases toward the detector circle.

JBO_19_5_056011_f004.png

While in Cartesian coordinates, pixels are coupled in a very intricate way by the spin, in polar coordinates, the blurring is constant in φ and independent of r (Fig. 4, right) and the problem simplifies considerably, with pixels at a given radius decoupled from all other pixels. We convert the n×n image, which is a discrete representation of p0α(r), given in Cartesian pixels, to a polar coordinate image Ypol, using spline interpolation to obtain pixel values at [n/2] values of r and 2n values of φ. Each column of Ypol corresponds to a particular radius, and therefore, the blur of pixels in that column is independent of pixels in every other column. Thus, we have [n/2] independent blurring problems, each of size 2n and all involving the same blurring kernel.

Let us, for the moment, omit the radial variable and write y(φ) for the blurred initial pressure p0α(r,φ) and x0(φ) for the desired unblurred version p0(r,φ). We then face the problem of estimating the signal φx0(φ) from the data

y(φ)=(x0*kα)(φ)+n(φ)forφ[0,2π].

Here n(·) is the deterministic or stochastic noise and kα(·) is the blurring kernel. The blurred initial pressure is itself the result of a reconstruction algorithm, such as time reversal, and therefore, the noise will be nonstationary and correlated even if the noise in the original data is white and Gaussian.

Deconvolution is a typical inverse problem, where the exact solution may either not exist or the noise is severely amplified. In order to make the reconstruction process well-posed, regularization techniques have to be applied, where approximate solutions are constructed that are stable with respect to data perturbations. The most common technique for solving a linear inverse problem is Tikhonov regularization. There an estimate for the unknown is defined as the unique minimizer of the Tikhonov functional

xx*kαy2+λx2.

The minimizer of the strictly convex Tikhonov functional can be found by setting its gradient to zero. Recalling that the convolution operation is diagonalized by the Fourier transform, the unique minimizer of the Tikhonov functional is given by x=y*gα, where the filter gα has the Fourier representation

F[gα](ω)=F[kα](ω)¯|F[kα](ω)|2+λ.

Here F[h](ω)=02πh(φ)exp(iωφ)dφ is the Fourier transform of some function h and the bar indicates complex conjugation. For signal and image deconvolution with additive random noise, this estimator is also known as Wiener deconvolution. The angular extent of the blurring kernel is given by 2α=2arctan(w/2rS), where w is the detector width and rS is the radius of the detection circle (Fig. 5). In the context of angular spin deblurring, Wiener deconvolution works best when combined with reconstruction algorithms that already use polar coordinates, such as the Fourier domain algorithms derived in Refs. 18 and 19. However, it may also be applied separately on images reconstructed by any other algorithm.

Fig. 5

Schematic of curved detector versus flat detector.

JBO_19_5_056011_f005.png

For our numerical implementation, we use a slightly different formulation:20,21 the [n/2] blurring problems are written in discretized form as Kxpolypol, where ypol is one column of Ypol. The blur matrix K is a 2n×2n circulant band matrix with band size proportional to the blurring angle and has an SVD K=UΣVT, where U and V are orthogonal matrices and Σ is diagonal. We again define regularized solutions as minimizers on R2n of the Tikhonov function xKxypol2+λx2. Using the SVD of the matrix K, the minimizer of the Tikhonov function is given by

xpol=(KTK+λI)1KTypol=V(Σ2+λI)1ΣUTypol.

Since any circulant N×N matrix has complex eigenvectors

vm=[1,exp(2πimN),,exp(2πim(N1)N)]form=0,,N1,
the SVD in this case essentially amounts to a discrete Fourier transform. A small difference to Wiener deconvolution as described above is that here discretization to a finite vector is done before Tikhonov regularization. In the context of general least squares problems, the first approach is often referred to as “solving the normal equations by eigen-decomposition” as opposed to “solving via SVD.”22 The matrix formulation is also convenient to apply the standard generalized cross-validation (GCV) algorithm23 to determine an optimal regularization parameter λ for each image and to diminish the influence of noise. The columns xpol give us a matrix Xpol, which we transform back to Cartesian pixels X, again using spline interpolation. Pixels in X outside the circular domain are set to zero. The image X is then a discrete approximation of the unblurred initial pressure p0(x,y) in Cartesian coordinates. We note that with minor modifications, the squared 2-norm x2 in the Tikhonov functional can be replaced by the squared 2-norm of a solution derivative. However, replacing x2 by a nonquadratic penalty makes the minimization of the resulting Tikhonov functional more difficult.

In Ref. 24, Hansen described such a variant of TSVD that seeks to regularize with the 1-norm of a solution derivative (in our case, the first derivative) instead of regularizing with the 2-norm of the solution. Recall that the 1-norm of a vector is defined by the sum of the absolute values of all components. The 1-norm regularization produces a vector of delta functions (sparse spike train) for the first derivative, and therefore, the solutions will be piecewise constant functions with at most k discontinuities, where k is the truncation index in the TSVD. If higher solution derivatives are used for regularization, piecewise polynomial solutions are obtained, hence the name PP-TSVD. Obviously, this deblurring method will work best for images that are themselves better described by piecewise constant functions (having sharp edges) than by smooth functions. The method was reported in the context of photoacoustic imaging with planar detection geometry in Ref. 6, where it was demonstrated on a phantom of two cylinders of constant optical absorbance. In Ref. 7, it is mentioned that it could also be transferred to a cylindrical detection geometry.

In Ref. 25, Hansen discusses PP-TSVD for deconvolution of 2-D images in the case of a separable 2-D convolution kernel. In our case, the kernel (when expressed in polar coordinates) is not only separable, but also independent of r, so that the radial part after discretization involves only the identity matrix. The extension to 2-D as described by Hansen (using the Kronecker product of the two kernel factor matrices), in this case, reads as follows: apply the deconvolution method using the same blurring matrix K as for Wiener deblurring; for each column representing a fixed radius, execute a one-dimensional PP-TSVD, where the algorithm is taken from Ref. 26. The optimal choice of the regularization parameter k was not considered in Ref. 25. In the section on reconstructions from computer-simulated data, we will determine an optimal k empirically by comparing reconstruction quality for a sample image and a series of k-values.

2.3.

Approximation of a Curved Detector by a Flat Detector

For the measurements described in Sec. 4.1, a 3×25-mm PVDF-foil detector was used, which is a good approximation of an extended integrating line detector (Fig. 3). However, this detector is not curved but flat. In this section, we investigate the difference between curved and flat aperture analytically.

Omitting the z-coordinate, the pressure data from the experiment are integrated over a line of length w tangent to a detection circle of radius rS (Fig. 5). Denoting by l=l(α)=rStan(α) the distance of a point on the line segment from its center, using a first-order Taylor approximation of p(r,φ,t) in the radial variable (around the center r=rS), and applying the approximation 1+tan(α)2=1+[(α2)/2]+O(α4) for small α, we obtain

w/2w/2p(r,φ,t)dl=ααp[rS1+tan(α)2,φ0+α,t][1+tan(α)2]rSdα=ααp(rS,φ0+α,t)rSdα+ααp(rS,φ0+α,t)α2rSdα+ααpr(rS,φ0+α,t)α22rS2dα+O(α4).

This implies that the difference

w/2w/2p(r,φ,t)dlααp(rS,φ0+α,t)rSdα,
between the pressure integrals over curved detector and flat detector is of order O(α3). As a consequence, the measured signals for curved and flat detector shapes agree up to second order in α. Hence, for small apertures 2α (in our experiment, 2α=10.7deg), the measured signals will differ little for the two detector types and, therefore, also the blurred reconstructions (which follows, for example, from the stability of the time reversal algorithm).

However, if the aperture gets larger (>20deg), a blur that was of constant width in the radial coordinate will split up at both ends for the flat detector. In fact, the ends look like a pair of tongs grasping two balls with radii proportional to the largest deviation of the flat detector from the detection circle [Fig. 6(b)].

Fig. 6

Reconstruction of a ball of 1-mm diameter located 12 mm from detection center (rS=16mm) from a flat strip detector with aperture 20 deg (a) and 40 deg (b). The bold white bar indicates a flat detector that is tangent to the detection circle (dashed white). In both cases, we applied a reconstruction algorithm that is theoretically exact for point-like data.

JBO_19_5_056011_f006.png

3.

Computer-Simulation Studies

For testing the two deconvolution methods described in the previous section, we use a phantom similar to the four-point phantom shown in Fig. 4. The points are refined in their structure and described by spatial profiles that are either hollow cylinders (HC) or Gaussians (G). The HC profile type was slightly smoothed by convolving with a box function to make it look more realistic and to make it more interesting for PP-TSVD. We want to examine reconstruction quality in dependence of source distance from center and of detector aperture. Quality will be assessed by visual inspection of 2-D profiles; after visual inspection, it will become apparent which quantitative measures are additionally useful.

In all numerical studies, the solution of the 2-D wave equation has been computed by numerically solving the wave equation. The image reconstruction (prior to deconvolution) has been performed by time reversal. The computational domain for the simulations is 2×2mm, the detection circle radius rS=0.8mm, and the spatial grid size Δx is 2/500mm for the forward computation and 2/441mm for time-reversed computation. Also, the timestep size Δt as well as the ending and starting times T0 in the forward and time-reversed computations are chosen somewhat different to avoid an inverse crime. Both the forward computations (using p(k,nΔt)=p(k,0)·cos(cknΔt), Δt=1.33ns, sound speed c=1500m/s) and the time-reversed computations [using Eq. (1) and Δt=1.68ns] are performed in k-space. Pressure signals are recorded at 720 equidistant points around the detection circle. A curved integrating detector of aperture 2α degrees is simulated by summing the recorded signals of consecutive measurement points lying within the designated aperture and dividing by the number of these points (in order to make signal sizes comparable to the point detector case, although physically a summation and not an averaging is done). The detector aperture 2α (in degrees) is chosen in a way that 2α/360·720=4α is an integer, which is the case for selecting 2α=10deg or 2α=20deg. After image reconstruction, we present profiles through the center of a point object of the phantom. These profiles are drawn over abscissas representing a pixel sequence in either the radial or tangential direction. Note that the profiles in the radial and tangential direction are equal for the initial pressure, but will be different for the reconstructions. For demonstrating the effect of source distance from center, we display the profiles for the first and the fourth point on the spiral (counted from center). The effect of detector aperture is shown by displaying the results for 2α=10deg and 2α=20deg. In summary, our tests with simulated data include 32 simulations (two choices for methods, points, profiles, apertures, and directions).

At this point, we want to determine an optimal parameter k for PP-TSVD deconvolution. The effect of k on deconvolution is examined using k=80, k=200 k=n=442 (the reconstructed image is n×npixles), and k=850. In Fig. 7, a close-up of the third and fourth point reconstructions for the HC profile is provided. For the annular xy-profiles of the HC sources, the reconstruction quality is distinctly better along a cross-section in radial direction than along a cross-section in tangential direction. This feature can be observed in all deblurred reconstructions, including those using Wiener deconvolution. However, with PP-TSVD the shape of the reconstruction is more sharp-edged for k=80, and is smoother for the higher-k reconstructions. On the other hand, computing time increases about linearly with k (ranging from 30 min for k=80 to >3h for k=850); moreover, the background noise level increases somewhat with k, becoming intolerable when k approaches 2n. Inspecting Fig. 7, we decided to use k=n=442 in our simulations. In any case, choosing k=n, i.e., truncating at half matrix size, seems to be reasonable for the other images as well in our scenario if no additional noise is added.

Fig. 7

Comparison of reconstructions for third and fourth point [hollow cylinders (HC) profile] for varying piecewise polynomial truncated singular value decomposition (PP-TSVD) regularization parameter: (a) k=80, (b) k=200, (c) k=442, and (d) k=850.

JBO_19_5_056011_f007.png

For the visual inspection, we display the results on four 4-subplot figures. Each four-subplot figure shows four curves corresponding to initial pressure (solid), reconstruction without deconvolution (dash-dotted), reconstruction after Wiener deconvolution (dashed), and reconstruction deconvolved with PP-TSVD (dotted). The four subplots in each figure show the profiles through the two chosen points in the two directions. The amplitude range for display is always [0.1 1.1], assuming that the initial pressure is normalized to [0,1]. The following general observations pertain to reconstructed amplitudes: with point detectors, our algorithm returns amplitude values in the range [0.12, 0.90] for our phantoms. The values returned with the integrating (here averaging) detectors lie in the range [0.08, 0.78] because they are smeared over a larger area. Deconvolution provides partial amplitude compensation: values lie in the range [0.13, 0.83]. Note that these numbers apply only to the point closest to the center that is reconstructed best; for the other three points, the amplitude degradation is even higher.

The inspection of Figs. 8 and 9 yields two surprising insights: both the deconvolution method and the detector aperture hardly have any influence on the deblurred results. There is some amplitude decrease between the first and fourth point as mentioned before, and the difference of reconstruction quality in radial and tangential directions is very pronounced and affects both amplitude and full width at half maximum (FWHM) of the profiles. A distinctive feature is also the much stronger smoothing of the hollow interior of the HC profile for the fourth point, compared to the first point. In Figs. 10 and 11, the same simulations were performed for Gaussian (G) initial profiles. The two insights gained with the HC profiles are confirmed once more: there is virtually no difference among (deblurred) results for the two deconvolution methods and two detector apertures. The amplitude decrease from first point to fourth point is as expected and the difference between radial and tangential profiles is again pronounced, but here it mainly affects the width of the Gaussian profile.

Fig. 8

Detector aperture 2α=10deg: four subplots, each showing profiles of initial pressure HC (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point r-direction (a), first point φ-direction (b), fourth point r-direction (c), and fourth point φ-direction (d).

JBO_19_5_056011_f008.png

Fig. 9

Detector aperture 2α=20deg: four subplots, each showing profiles of initial pressure HC (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point r-direction (a), first point φ-direction (b), fourth point r-direction (c), and fourth point φ-direction (d).

JBO_19_5_056011_f009.png

Fig. 10

Detector aperture 2α=10deg: four subplots, each showing profiles of Gaussian initial pressure (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point r-direction (a), first point φ-direction (b), fourth point r-direction (c), and fourth point φ-direction (d).

JBO_19_5_056011_f010.png

Fig. 11

Detector aperture 2α=20deg: four subplots, each showing profiles of Gaussian initial pressure (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point r-direction (a), first point φ-direction (b), fourth point r-direction (c), and fourth point φ-direction (d).

JBO_19_5_056011_f011.png

As a quantitative measure we list in Table 1 the FWHM of all the Gaussian-like profiles in Figs. 10 and 11. The measures are given in pixel units Δx (spatial step size of the time-reversed grid). Linear interpolation was done between the integer pixel grid, so that accuracy is ±0.05Δx. As can already be seen from visual inspection, the deconvolved reconstruction significantly improves resolution compared to the blurred image in terms of reducing the FWHM in φ-direction of isolated objects. However, there are other aspects of degraded image quality, such as the filling up of holes, which can be seen in Figs. 8(c), 8(d), 9(c), and 9(d), that cannot be avoided even with deconvolution and that are strongly dependent on source distance from center.

Table 1

Full width at half maximum (FWHM) in pixel units of the profiles displayed in Figs. 10 and 11.

FWHM in pixelsFig. (a) P1 - rFig. (b) P1 - φFig. (c) P4 - rFig. (d) P4 - φ
2α=10deg, initial11.4511.4511.4511.45
2α=10deg, blurred10.6014.2011.7519.10
2α=10deg, Wiener10.5013.2011.8015.50
2α=10deg, PP-TSVD10.5013.2011.4015.50
2α=20deg., initial11.4511.4511.4511.45
2α=20deg., blurred11.0518.3011.60>40
2α=20deg, Wiener10.5512.5012.0016.25
2α=20deg, PP-TSVD10.5512.5011.6515.75
Note: PP-TSVD, piecewise polynomial truncated singular value decomposition.

An explanation of the similarity between PP-TSVD and Wiener deconvolution is that the high number (442) of possible discontinuities produces quasi-smooth solutions in the PP-TSVD case also (note that our plots show only 27 or 33 pixel values) so that the different regularizers produce very similar results. The constancy in results for varying aperture could be confirmed for much higher apertures also. The noticeable amplitude degradation for an increased distance of the source from the center may result from the fact that a uniform grid spacing in polar coordinates transforms to a grid spacing in Cartesian coordinates that decreases in density with distance from center. Since our algorithms work separately on each radius, using a denser angular grid for large radii might reduce this effect. Such an investigation is beyond the scope of our present work.

Computation time for Wiener deblurring is 1 s; for PP-TSVD deblurring with an appropriate k, it is at least an hour if the code pptsvd.m from Ref. 26 is used. The documentation for this algorithm (see Ref. 25, Sec. 5) reveals that it was written especially for large-sized problems. For our small-sized problems, the following measures would make the algorithm faster: (1) instead of MATLAB® function svds, use svd to perform SVD; (2) optimize the function l1c to solve a discrete 1 linear optimization. We estimate that with these measures, the algorithm for 442×442pixel reconstruction in this section and for 251×251pixel reconstruction in Sec. 4.2 could be made as fast as a few seconds. Still, the Wiener deconvolution is faster and provides the same image quality.

So we can currently recommend the Wiener deblurring for problems of the type presented in this work. It must be stressed, however, that image-specific regularization parameter optimization via GCV or a similar algorithm is essential for this excellent quality of Wiener deblurring. With these provisions, we were not able to observe any instabilities or oscillations as was reported for planar geometry in Ref. 6. Further numerical studies shown below demonstrate that these conclusions are also valid in the case of noisy data.

It is well known that the performance of deblurring algorithms often strongly depends on the noise statistics in the data. For our algorithms, we investigated this issue numerically by adding Gaussian white noise to the numerically computed solution of the 2-D wave equation and applying time reversal to the noisy data. The results of Wiener deblurring and the PP-TSVD deblurring algorithm on the intermediate images are shown in Fig. 12. For Wiener deconvolution, the regularization parameter has again been computed in a data-driven way by applying the GCV algorithm. Note that the regularization parameter computed by GCV increases with increasing noise level. The result is λ=9.2×104 for the simulated data without noise, λ=2.5×103 for 3% noise, and λ=2.7×103 for 10% noise. Such a behavior is expected for an ill-conditioned problem since the regularization parameter acts as trade-off between accuracy for exact data and stability with respect to noise. For the PP-TSVD algorithm, the parameter 1/k acts as regularization parameter, and therefore, k should be chosen smaller for higher noise levels. Optimal values for k depending on the noise level have been found empirically by performing similar experiments as in the case of noise-free data. For the result shown in Fig. 12, we found k=300 to be optimal for 3% noise and k=120 for 10% noise.

Fig. 12

Reconstruction results using detector aperture 2α=10deg after adding 3 and 10% Gaussian white noise, respectively, to the solution of the 2-D wave equation. (Here by p% noise we mean that the standard deviation of the added random variable is p% of the maximal data value.) The regularization parameters for the Wiener deconvolution (red, dashed) have been determined by the generalized cross-validation. For PP-TSVD deblurring (black, dotted), we have taken k=300 in the case of 3% noise and k=120 in the case of 10% noise, compared to k=442 (half dimension of the matrix) without noise in Figs. 10(c) and 10(d).

JBO_19_5_056011_f012.png

It can be seen that both deblurring algorithms still provide a comparable image quality. Since time reversal is a linear algorithm, the noise in the blurred intermediate images is still Gaussian but neither uncorrelated nor stationary. However, because of the rotational invariance of the image reconstruction problem, the noise is stationary in the angular direction (i.e., long the angular direction, the variance of the noise is constant and the correlation depends only on angular spacing). This may be exploited to further improve the performance of the Wiener deconvolution algorithm. A precise analysis of the noise statistics, however, is beyond the scope of this paper.

4.

Experimental Results

4.1.

Experimental Setup

Top and front views of the measurement setup are shown in Fig. 13. A homemade piezoelectric ultrasound detector was used to record the acoustic signals. The detector was made of a strip of PVDF-foil with a thickness of 28 μm, with electrodes deposited on either side attached on a polymer substrate (backing material). The width w and the length L of the stripe (sensitive area) were 3 and 25 mm, respectively. The signals were recorded by connecting the detector with an active probe to the digital storage oscilloscope. To improve the SNR, the signals were recorded with bandwidth limitation (dc up to 20 MHz) by the oscilloscope. This is the only limitation to the bandwidth of the ultrasonic transducer.

Fig. 13

Top and front view of the experimental setup.

JBO_19_5_056011_f013.png

The cylindrical phantom sample (D=20mm) was made of plastisol with six channels oriented along z-direction, distributed along a spiral curve in the cross-section area (see top view in Fig. 13). The sound speed of plastisol is 1400m/s; further acoustic properties can be found in Ref. 27. The channels, each with a diameter of 1 mm, were filled with a light-absorbing dye solution (OrangeG, absorption coefficient 100cm1).

For photoacoustic excitation, the sample was illuminated by 10-ns laser pulses from a frequency-doubled, Q-switched Nd:YAG laser with a repetition rate of 10 Hz. The output beam had a wavelength of 532 nm and was split into two beams to illuminate the sample from two opposite sides parallel to the y-direction. Furthermore, the beams were focused with cylindrical lenses to narrow the extension of the illumination along z-direction (1.5mm), forming a line-shaped illumination pattern on the surface of the sample (see front view in Fig. 13). The rotation axis of the sample was positioned at a normal distance rS=16mm to the detector and was oriented parallel to the z-direction. The z-position of the line-shaped illumination was fixed to the height of the center of the z-extension of the detector. Signals were acquired while rotating the sample over five full rotations with an angular increment of 0.9 deg. The five times averaged temporal signals are shown for the 400 orientations over 360 deg in Fig. 14. The dye solution in the holes is strongly absorbing, so the absorbed energy density decreases rapidly inside the holes. This can already be seen from the sinogram in Fig. 14 (signals are generated only near the boundaries of the absorbers), but more clearly in the reconstructed images in Sec. 4.2.

Fig. 14

Measured pressure signals (sinogram) as a function of the rotation angle (x axis) and the time; on the y axis, the time has been multiplied by the sound speed in water (c=1500m/s). Already in the sinogram it can be seen that the laser light is strongly absorbed by the dye solution in the holes.

JBO_19_5_056011_f014.png

4.2.

Reconstruction and Deblurring Results

The sinogram in Fig. 14 shows data measured over 400 equidistant angular positions using a sampling interval of Δt=10ns. Note that the y axis displays r=ct (time multiplied by sound speed in water). Beginning at r0=7mm, 1210 samples were recorded, ending at r1=25mm. The computational domain was 33.6×33.6mm, and the spatial grid size was 33.6/400mm. A 2-D time-reversal reconstruction was performed using Eq. (1). The acoustic properties of the plastisol27 were considered to be similar enough to those of water so that the two media were not modeled separately.

For closer inspection of the reconstructed sources, which are located in the inner half of the detection area (rS=16mm), only the inner 251×251pixel square was cut out of the original 401×401pixel reconstructed image; then the deblurring algorithms were applied to this inner square. Computation time is 0.5 s for Wiener deconvolution and 10min for PP-TSVD deconvolution with pptsvd.m. The choice for the regularization parameter k=n is also a reasonable one here.

In Figs. 15 and 16, we see cross-sections through the centers of the six holes in the plastisol phantom in the radial and tangential directions, respectively. Again, the similarity between the results from the two deblurring methods is remarkable, although there are more differences from the simulated data (see fifth hole in Fig. 15, for instance). Amplitude compensation by deblurring is somewhat better in the radial direction than in the tangential direction, whereas FWHM compensation and, thus, resolution is similar for the two directions. In the reconstructions from our experimental data, we observe some unusually large negative pixel values. Reasons that this occurs may be as follows: modeling the plastisol medium the same as water, modeling the temporal impulse response of the piezoelectric detector as a delta function, disturbances or feedback during electronic acquisition. However, looking at the 2-D reconstructed images (Fig. 17, where negative pixel values were set to zero), we conclude that these imperfections do not seriously affect the overall quality of the information about the optical absorbance of the dye solution in the holes. The reconstructions show that most of the absorption occurs near the boundaries of the holes, which is in agreement with the absorption properties of Orange G given in Sec. 4.1 and with the images of the original phantom drawn in Fig. 13.

Fig. 15

Slices through the six reconstructed holes in (r,z)-direction with hole distance from center increasing from left to right: before deblurring (blue, solid), after Wiener deblurring (red, dashed), and after PP-TSVD deblurring (black, dotted).

JBO_19_5_056011_f015.png

Fig. 16

Slices through the six reconstructed holes in (φ,z)-direction with hole distance from center increasing from left to right: before deblurring (blue, solid), after Wiener deblurring (red, dashed), and after PP-TSVD deblurring (black, dotted).

JBO_19_5_056011_f016.png

Fig. 17

Reconstructed xy-profiles of the entire six-hole phantom: before deblurring (left), after Wiener deblurring (center), and after PP-TSVD deblurring (right).

JBO_19_5_056011_f017.png

5.

Conclusion

The mathematical model of photoacoustic wave propagation is the 3-D linear wave equation. If line detectors are used that extend to infinity in one direction (or are at least sufficiently long), the reconstruction process typically involves the inversion of a 2-D wave equation for the pressure projections onto the remaining dimensions, followed by a 2-D inverse Radon transform. Considering only the first problem, we are confronted with a cylindrical geometry. Then the ideal 2-D model still is an abstraction from reality if (piezoelectric) detectors have sizeable extent in angular coordinates.

If the detector size is increased from an (idealized) line detector to a finite width, the resolution starts to deteriorate in parts of the image. Xu and Wang have first quantified this effect of the detector size on the spatial resolution in Ref. 28. Haltmeier and Zangerl have explained these effects as a consequence of the rotational and translational invariance of the wave equation.14 They have described the effect of the finite detector size in terms of a blurring kernel that has to be convolved with the image obtained from point or line detectors. This suggests the use of deconvolution methods to reverse the effect.

We have demonstrated with computer simulations and with an example on experimental data that both the Wiener and the PP-TSVD deconvolution methods yield an excellent compensation of the blurring. The computational cost of Wiener deblurring is very small (<1s); the cost for PP-TSVD deblurring is significantly larger, but still in the range of a few seconds (assuming image sizes of a few hundred pixels squared). So there is every reason to apply Wiener deblurring, together with appropriate optimization of the regularization parameter, for any blurring kernel that is explicitly known. We did not observe any of the problems with Wiener deblurring that were reported in Ref. 6. If a 3-D object is to be imaged, this deconvolution can be applied at the first step that involves the inversion of the 2-D wave equation. After correction of the 2-D projections, the final step using the inverse Radon transform can be applied as usual for obtaining a full 3-D image.

Matrix-model based reconstruction algorithms8,9 (or Refs. 29 and 30, using a different model) have the capability to incorporate arbitrary detector shapes and spatial impulse responses into the imaging model and, thus, into the reconstruction algorithm itself. However, such approaches require solving optimization algorithms working iteratively with very large-sized matrices and entail high computational effort. The conclusion of Ref. 9 states that “current implementations of (3-D) iterative image reconstruction algorithms still require several days to process a densely sampled data set.” In Ref. 29, we read “the run time for constructing and inverting the (2-D) model matrix for a 81×81grid was approximately one hour.” While this cannot be avoided for more complicated ultrasonic transducer shapes, we demonstrated that effective reduction of imaging artifacts in 2-D caused by the finite size of some regular-shaped detectors is possible in the framework of traditional reconstruction algorithms (like time-reversal or filtered backprojection) and takes only a fraction of the time needed with matrix-model algorithms.

Appendices

Appendix:

Derivation of the Blurring Kernel

Let p(r,φ,t) denote acoustic pressure expressed in polar coordinates, which is the unique solution of the two-dimensional wave equation 2/t2p(r,φ,t)=c2Δp(r,φ,t), with initial conditions p(r,φ,0)=p0(r,φ) and /tp(r,φ,0)=0. If we know all pointwise values p(rS,φ,t) on a detection circle of radius rS, then we can uniquely recover the initial pressure p0(r,φ) by either time reversal, Fourier domain algorithms, or filtered backprojection algorithms. Suppose now that we observe integrated values pα(rS,φ,t)=ααp(rS,φ+α,t)kα(α)dα of the pressure over a finite angular region, with kα(α) denoting some angular function. Any reconstruction method R0 that is exact for data p(rS,φ,t) will result in blurred reconstruction when applied to the data pα(rS,φ,t). In this appendix, we show that R0[pα](r,α)=p0α(r,α) with

p0α(r,α)ααp0(r,φ+α)kα(α)dα.

In fact, this identity is an easy consequence of the rotational invariance of the wave equation. To see this, one considers the function (r,φ,t)pα(r,φ,t)=ααp(r,φ+α,t)kα(α)dα. Then the rotational invariance and the linearity of the wave equation imply that pα(r,φ,t) satisfies the wave equation 2/t2pα(r,φ,t)=c2Δpα(r,φ,t). Further, we have pα(r,φ,0)=p0α(r,α) and /tpα(r,φ,0)=0. So the data pα(rS,φ,t) is the same as the signal that would be measured from the initial pressure p0α(r,α). Hence, as claimed, application of any reconstruction method R0 to the data pα(rS,φ,t) yields R0[pα](r,α)=p0α(r,α).

Acknowledgments

This work has been supported by the Austrian Science Fund (FWF), project numbers TRP102-N20, S10502-N20, S10503-N20, and S10508-N20, the European Regional Development Fund (EFRE) in the framework of the EU-program Regio 13, and the federal state Upper Austria. The work of O’Leary was supported by the USA National Science Foundation under grant DMS 1016266.

References

1. 

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

2. 

L. V. WangS. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 (6075), 1458 –1462 (2012). http://dx.doi.org/10.1126/science.1216210 SCIEAS 0036-8075 Google Scholar

3. 

C. LutzweilerD. Razansky, “Optoacoustic imaging and tomography: reconstruction approaches and outstanding challenges in image performance and quantification,” Sensor, 13 (6), 7345 –7384 (2013). http://dx.doi.org/10.3390/s130607345 SNSRES 0746-9462 Google Scholar

4. 

P. Burgholzeret al., “Thermoacoustic tomography with integrating area and line detectors,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 52 (9), 1577 –1583 (2005). http://dx.doi.org/10.1109/TUFFC.2005.1516030 ITUCER 0885-3010 Google Scholar

5. 

P. Burgholzeret al., “Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Probl., 23 (6), S65 –S80 (2007). http://dx.doi.org/10.1088/0266-5611/23/6/S06 INPEEY 0266-5611 Google Scholar

6. 

Y. XuD. FengL. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic—I: Planar geometry,” IEEE Trans. Med. Imaging, 21 (7), 823 –828 (2002). http://dx.doi.org/10.1109/TMI.2002.801172 ITMID4 0278-0062 Google Scholar

7. 

Y. XuM. XuL. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic—II: Cylindrical geometry,” IEEE Trans. Med. Imaging, 21 (7), 829 –833 (2002). http://dx.doi.org/10.1109/TMI.2002.801171 ITMID4 0278-0062 Google Scholar

8. 

K. Wanget al., “An imaging model incorporating ultrasonic transducer properties for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 30 (2), 203 –214 (2011). http://dx.doi.org/10.1109/TMI.2010.2072514 ITMID4 0278-0062 Google Scholar

9. 

K. Wanget al., “Accelerating image reconstruction in three-dimensional optoacoustic tomography on graphics processing units,” Med. Phys., 40 (2), 023301 (2013). http://dx.doi.org/10.1118/1.4774361 MPHYA6 0094-2405 Google Scholar

10. 

M. XuL. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 (1), 016706 (2005). http://dx.doi.org/10.1103/PhysRevE.71.016706 PLEEE8 1063-651X Google Scholar

11. 

P. Burgholzeret al., “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E, 75 (4 Pt 2), 046706 (2007). http://dx.doi.org/10.1103/PhysRevE.75.046706 PLEEE8 1063-651X Google Scholar

12. 

P. Burgholzeret al., “Deconvolution algorithms for photoacoustic tomography to reduce blurring caused by finite sized detectors,” Proc. SPIE, 8581 858137 (2013). http://dx.doi.org/10.1117/12.2003889 PSISDG 0277-786X Google Scholar

13. 

T. Bereret al., “Remote photoacoustic imaging on non-flat surfaces and appropriate reconstruction algorithms,” Proc. SPIE, 8581 858134 (2013). http://dx.doi.org/10.1117/12.2003566 PSISDG 0277-786X Google Scholar

14. 

M. HaltmeierG. Zangerl, “Spatial resolution in photoacoustic tomography: effects of detector size and detector bandwidth,” Inverse Probl., 26 (12), 125002 (2010). http://dx.doi.org/10.1088/0266-5611/26/12/125002 INPEEY 0266-5611 Google Scholar

15. 

G. Paltaufet al., “Photoacoustic tomography using a Mach-Zehnder interferometer as an acoustic line detector,” Appl. Opt., 46 (16), 3352 –3358 (2007). http://dx.doi.org/10.1364/AO.46.003352 APOPAI 0003-6935 Google Scholar

16. 

T. Bereret al., “Characterization of broadband fiber optic line detectors for photoacoustic tomography,” J. Biophotonics, 5 (7), 518 –528 (2012). http://dx.doi.org/10.1002/jbio.201100110 JBOIBX 1864-063X Google Scholar

17. 

G. PaltaufR. NusterP. Burgholzer, “Characterization of integrating ultrasound detectors for photoacoustic tomography,” J. Appl. Phys., 105 102026 (2009). http://dx.doi.org/10.1063/1.3116133 JAPIAU 0021-8979 Google Scholar

18. 

L. Kunyansky, “Fast reconstruction algorithms for the thermoacoustic tomography in certain domains with cylindrical or spherical symmetries,” Inverse Probl. Imaging, 6 (1), 111 –131 (2012). http://dx.doi.org/10.3934/ipi 1930-8337 Google Scholar

19. 

M. Haltmeieret al., “Thermoacoustic tomography and the circular Radon transform: exact inversion formula,” Math. Models Methods Appl. Sci., 17 (4), 635 –655 (2007). http://dx.doi.org/10.1142/S0218202507002054 MMMSEU 0218-2025 Google Scholar

20. 

J. M. ChungG. EasleyD. P. O’Leary, “Inverting radial and angular blur using polar pixels,” Google Scholar

21. 

S. RibaricM. MilaniZ. Kalafatic, “Restoration of images blurred by circular motion,” in First Int. Workshop on Image and Signal Processing and Analysis, 53 –60 (2000). Google Scholar

22. 

J. NocedalS. Wright, Numerical Optimization, 2nd ed.Springer, Wien, New York (2006). Google Scholar

23. 

G. H. GolubM. HeathG. Wahba, “Generalized cross-validation as a method for choosing a good ridge parameter,” Technometrics, 21 (2), 215 –223 (1979). http://dx.doi.org/10.1080/00401706.1979.10489751 TCMTA2 0040-1706 Google Scholar

24. 

P. C. Hansen, “Piecewise polynomial solutions without a priori break points,” Numer. Linear Algebra Appl., 3 (6), 513 –524 (1996). http://dx.doi.org/10.1002/(ISSN)1099-1506 1070-5325 Google Scholar

25. 

P. C. Hansen, “The PP-TSVD algorithm for image restoration problems,” Methods and Applications of Inversion, 171 –186 (2013) Google Scholar

26. 

P. C. Hansen, “Regularisation tools: a Matlab package for analysis and solution of discrete ill-posed problems,” Numer. Algorithms, 6 (1), 1 –35 (1994). http://dx.doi.org/10.1007/BF02149761 NUALEG 1017-1398 Google Scholar

27. 

G. M. Spironet al., “Optical and acoustic properties at 1064 nm of polyvinyl chloride-plastisol for use as a tissue phantom in biomedical optoacoustics,” Phys. Med. Biol., 50 (14), N141 –N153 (2005). http://dx.doi.org/10.1088/0031-9155/50/14/N01 PHMBA7 0031-9155 Google Scholar

28. 

M. XuL.V. Wang, “Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction,” Phys. Rev. E, 67 (1), 056605 (2003). PLEEE8 1063-651X Google Scholar

29. 

A. RosenthalD. RazanskyV. Ntziachristos, “Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography,” IEEE Trans. Med. Imaging, 29 (6), 1275 –1285 (2010). http://dx.doi.org/10.1109/TMI.2010.2044584 ITMID4 0278-0062 Google Scholar

30. 

A. RosenthalV. NtziachristosD. Ranzansky, “Model-based optoacoustic inversion with arbitrary-shape detectors,” Med. Phys., 38 (7), 4285 –4295 (2011). http://dx.doi.org/10.1118/1.3589141 MPHYA6 0094-2405 Google Scholar

Biography

Heinz Roitner received his PhD in applied mathematics from the University of Arizona in 1991. From 1998 to 2009, he was a software engineer with Siemens AG in Vienna, before he joined the photoacoustic research group at RECENDT. His research interests include partial differential equations and inverse problems applied to tomography as well as signal and image processing.

Markus Haltmeier received his PhD in mathematics from the University of Innsbruck, Austria, in 2007. Since 2012, he has been a full professor for mathematics at the University of Innsbruck.

Robert Nuster received his PhD in experimental physics from Karl-Franzens-University Graz, Austria, in 2007 with a thesis on development and application of optical sensors for laser induced ultrasound detection. From 2008 to 2011, he was a postdoctoral research fellow at the Karl-Franzens-University Graz. Since 2011, he has been the leader of the project: “Optical parallel detection for photoacoustic tomography.” His current research interests include photoacoustic imaging, characterization of materials by laser ultrasound and ultrasound sensor development.

Dianne P. O’Leary is a professor of computer science at the University of Maryland and also holds an appointment in the university’s Institute for Advanced Computer Studies (UMIACS). She earned a BS from Purdue University and a PhD from Stanford University. Her research is in computational linear algebra and optimization, with applications including solution of ill-posed problems, image deblurring, information retrieval, and quantum computing. She is a fellow of SIAM and ACM.

Thomas Berer received his PhD degree in technical science from the Department of Semiconductor Physics of the Johannes Kepler University of Linz, Austria, in 2007. Since 2007, he has been working for the Sensor Department of the Upper Austrian Research GmbH, which became the RECENDT GmbH in 2009. Since 2010, he has been the head of the photoacoustic imaging group. His research interests include photoacoustic imaging and laser ultrasound.

Guenther Paltauf received his PhD in physics from Karl- Franzens-University Graz, Austria, in 1993. After postdoctoral studies in Graz, the University of Bern, Switzerland, and the Oregon Medical Laser Center in Portland, Oregon, he is now associate professor at the Department of Physics at the University of Graz. His research interests are photoacoustic imaging, laser ultrasonics, and photomechanical laser-tissue interactions.

Hubert Grün has been working for the Sensor Department of the Upper Austrian Research GmbH, which became the Research Center for Non-Destructive Testing (RECENDT) in 2009, where he is now the head of the strategic department. His research interests include photoacoustic imaging and laser ultrasonics.

Peter Burgholzer received his PhD in physics from Johannes Kepler University Linz, Austria, in 1993. After postdoctoral studies in Linz, Austria, and Ispra, Italy, at the Joint Research Center (JRC) of the European Commission he has been head of the Sensor Department of the Upper Austrian Research from 2000 to 2009. Now, he is director of a Christian Doppler Laboratory for Photoacoustic Imaging and Laser Ultrasonics and CEO of RECENDT.

© 2014 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2014/$25.00 © 2014 SPIE
Heinz Roitner, Markus Haltmeier, Robert Nuster, Dianne P. O’Leary, Thomas Berer, Guenther Paltauf, Hubert Grün, and Peter Burgholzer "Deblurring algorithms accounting for the finite detector size in photoacoustic tomography," Journal of Biomedical Optics 19(5), 056011 (22 May 2014). https://doi.org/10.1117/1.JBO.19.5.056011
Published: 22 May 2014
JOURNAL ARTICLE
13 PAGES


SHARE
Lens.org Logo
CITATIONS
Cited by 18 scholarly publications.
Advertisement
Advertisement
Back to Top