Design, implementation, and analysis of a compressed sensing photoacoustic projection imaging system

Markus Haltmeier; Matthias Ye; Karoline Felbermayer; Florian Hinterleitner; Peter Burgholzer

doi:10.1117/1.JBO.29.S1.S11529

22 April 2024 Design, implementation, and analysis of a compressed sensing photoacoustic projection imaging system

Markus Haltmeier, Matthias Ye, Karoline Felbermayer, Florian Hinterleitner, Peter Burgholzer

Author Affiliations +

Journal of Biomedical Optics, Vol. 29, Issue S1, S11529 (April 2024). https://doi.org/10.1117/1.JBO.29.S1.S11529

Abstract

Significance

Compressed sensing (CS) uses special measurement designs combined with powerful mathematical algorithms to reduce the amount of data to be collected while maintaining image quality. This is relevant to almost any imaging modality, and in this paper we focus on CS in photoacoustic projection imaging (PAPI) with integrating line detectors (ILDs).

Aim

Our previous research involved rather general CS measurements, where each ILD can contribute to any measurement. In the real world, however, the design of CS measurements is subject to practical constraints. In this research, we aim at a CS-PAPI system where each measurement involves only a subset of ILDs, and which can be implemented in a cost-effective manner.

Approach

We extend the existing PAPI with a self-developed CS unit. The system provides structured CS matrices for which the existing recovery theory cannot be applied directly. A random search strategy is applied to select the CS measurement matrix within this class for which we obtain exact sparse recovery.

Results

We implement a CS PAPI system for a compression factor of 4:3, where specific measurements are made on separate groups of 16 ILDs. We algorithmically design optimal CS measurements that have proven sparse CS capabilities. Numerical experiments are used to support our results.

Conclusions

CS with proven sparse recovery capabilities can be integrated into PAPI, and numerical results support this setup. Future work will focus on applying it to experimental data and utilizing data-driven approaches to enhance the compression factor and generalize the signal class.

1. Introduction

Photoacoustic tomography (PAT) is an emerging non-invasive imaging technique that combines the high contrast of optical imaging with the high spatial resolution of ultrasound imaging.¹^–³ It is based on the generation of acoustic waves by illuminating a sample with picosecond or nanosecond optical pulses. The acoustic waves are measured outside the object, and mathematical algorithms are used to reconstruct an image of the inside. While there are many important practical and theoretical aspects along the lines of signal generation, signal detection, system design, image generation and enhancement, in this paper we focus on the measurement and inversion of acoustic waves.⁴^,⁵ Specifically, we focus on PA projection imaging (PAPI) based on integrating line detectors (ILDs).⁶^,⁷ Our goal is to use ideas from compressed sensing (CS) to reduce the number of spatial measurements compared to standard measurements where each ILD is used to record its own time-dependent signal. Specifically, we present our design and development of CS in PAT under physical constraints that naturally arise in the already existing self-developed PAPI system.⁸

1.1.

Photoacoustic Projection Imaging

A PA projection tomograph records the induced acoustic signals with an array of parallel ILDs, with each sensor integrating (averaging) the pressure along the lines of the detectors. The data thus consist of samples of the linear projection of the three-dimensional (3D) acoustic pressure wave in the direction of the ILDs. Reconstruction in two-dimensional (2D) gives a projection of the initial pressure distribution. If a 3D reconstruction is required, the object can be rotated around an axis perpendicular to the fibers, and a 3D reconstruction is computed from the collection of 2D projections by inversion of the 2D Radon transform, which is similar to parallel beam X-ray CT.⁹^,¹⁰ As in X-ray imaging, where in certain situations single projections are sufficient, the same can be said for photoacoustic imaging. We will therefore restrict ourselves to 2D PAPI.

Figure 1 shows a photograph of our self-developed all-optical PAPI system used in this study. The setup is based on fiber optic Mach–Zehnder interferometers (FOMZIs) with graded index polymer optical fibers (GIPOFs). These have a higher bandwidth than glass optical fibers and are more stable for measurements. In the current system, 64 ILDSa are arranged on a circle forming a cylinder. Readout for each sensor requires an analog-to-digital (AD) converter, and four sensors are multiplexed to one AD converter. Thus, to measure all 64 signals, the measurement process must be repeated four times. Our hypothesis is that proper combinations of ILD signals will be advantageous over recording individual signals when used in conjunction with a nonlinear CS recovery algorithm.

Fig. 1

Photographic image of the PA projection tomography with 64 FOMZIs as ILDs forming the basics of the presented research.

1.2.

Compressed Sensing in PAPI

Following the CS paradigm, instead of recording pressure signals $P = {[p_{1}^{T}, \dots, p_{n}^{T}]}^{T}$ where $p_{j}$ is the pressure signal (written as column vector) of the $j$ ’th ILD, we record CS data

Eq. (1)

y_{i} = (A P)_{i} = \sum_{j = 1}^{n} a_{i j} p_{j} for i \in {1,2, \dots, m},

where

A = (a_{i, j})_{i, j} \in R^{m \times n}

denotes the CS measurement matrix. Usually in CS, the measurement matrix is chosen randomly, since this gives an exact recovery of sparse vectors with a high probability for large

n, m

. However, in practice, and specifically in our application, the matrix

A

cannot be chosen completely at random. First, the measurements cannot combine all pressure values if they are not connected to the same controller. Second, the numbers

a_{i, j}

are often restricted to specific values, in our case, for example to 0 and 1. Finally, the dimensionality

n

in our case is small, which limits the applicability of existing asymptotic CS theory that applies to the limit

n, m \to \infty

.

The goal of this work is to design, analyze, and implement a CS strategy that can actually be realized with our PAPI system. Within the considered family of measurements, we investigate the optimal design of matrices. Due to the low dimensionality of CS matrices, even a small compression factor $n / m$ below 2 seems to be a substantial challenge.

1.3.

Outline

In this paper, we present our findings and results in building a CS-PAPI system. This development is based on several steps. First, we provide a rigorous description of the PAPI problem. In this context, we also provide an overview of the most important background knowledge required. Second, we introduce a novel class of CS measurements that are practically feasible and can be realized with the existing self-developed PAPI setup. Third, we present a concept of optimal measurement design that allows researchers and practitioners to strategically select measurements to maximize imaging accuracy for CS in PAPI and other imaging modalities. While these results are developed in the context of sparsity, we present an outlook for the use of more general signal classes potentially enabling data-driven machine learning methods. Finally, we go from theory to practice and show how these results can be translated into the experimental realization of CS-PAPI.

2. Background

In this section, we present the background of our work. This includes PAPI modeling (Sec. 2.1), sparse CS theory (Sec. 2.2), and the description of the self-developed PAPI system (Sec. 2.3).

2.1.

PA Projection Imaging

PA tomography is based on generating an acoustic wave inside some investigated object using short optical pulses. When measuring the pressure with ILDs, the imaging problem reduces to a 2D version of the standard problem⁹^,¹⁰ and in this work we consider the 2D version only. Further, we restrict ourselves to constant sound speed and a circular measurement geometry as shown in Fig. 2.

Fig. 2

(a) An object is illuminated with a short optical pulse; (b) the absorbed light distribution causes an acoustic pressure; and (c) the acoustic pressure is measured with ILD arranged on a circle.

Let us denote by $u : R^{2} \to R$ the 2D PA source distribution which is our image of interest and supposed to be enclosed by a circle $C_{R}$ of radius $R$ . The 2D projected pressure satisfies the 2D wave equation

Eq. (2)

\partial_{t}^{2} p (r, t) - v_{s}^{2} Δ_{r} p (r, t) = δ^{'} (t) u (r) for (r, t) \in R^{2} \times R_{+},

where

δ^{'} (t)

is the first time derivative of the Dirac delta distribution,

r \in R^{2}

is the spatial location,

t \in R

is the time variable,

Δ_{r}

is the spatial Laplacian, and

v_{s}

is the constant speed of sound. The wave Eq. (2) is augmented with

p (r, t) = 0

for

t < 0

such that the acoustic pressure is uniquely defined as solution of Eq. (2). We rescale time in such a way that

v_{s} = 1

.

PAPI in circular geometry consists of recovering the function $u$ from measurements of $W u (s, t) = p (s, t)$ made on $C_{R} \times (0, \infty)$ . In the case of full data, exact and stable PA image reconstruction is possible and several efficient methods for recovering $u$ are available. We will use the FBP formula derived in Ref. 11:

Eq. (3)

u (r) = - \frac{1}{π R} \int_{C_{R}} \int_{| r - s |}^{\infty} \frac{(\partial_{t} t W u) (s, t)}{\sqrt{t^{2} - {| r - s |}^{2}}} d t d C (s) .

Note the inversion operator in Eq. (3) is also the adjoint of the forward operator $W$ . This in particular implies that inverting $W$ is stable.

In practical applications, the acoustic pressure can only be measured with a finite number of acoustic detectors. The standard sampling scheme in a circular geometry assumes uniformly sampled values

Eq. (4)

p (s_{j}, t_{ℓ}) for (j, ℓ) \in {1, \dots, n} \times {1, \dots, q},

with

s_{j} ≜ R (\cos (Ω (j - 1) / n), \sin (Ω (j - 1) / n))

,

t_{ℓ} ≜ 2 R (ℓ - 1) / (q - 1)

, and

Ω \leq 2 π

denoting the angular covering on the detection circle. The number

n

of detector positions in Eq. (4) is directly related to the resolution of the final reconstruction. Namely,

n \geq 2 R λ

equally spaced transducers covering the full circle are required to stably recover any PA source

u

that has maximal essential wavelength

λ

; see Ref. 12. Image reconstruction in this case can be performed by discretizing the inversion Eq. (3). The sampling condition requires a very high sampling rate, especially when the PA source contains narrow features, such as blood vessels or sharp interfaces. Commonly,

λ

will be determined by the spatial sampling via the Nyquist condition, such that

2 R λ = π N_{r}

, where

N_{r} \times N_{r}

is the number of samples for discretizing the object of interest on the square

[- R, R] \times [- R, R]

. In this case, we get

n = round (π N_{r} / 2)

for correct sampling according to Shannon sampling theory.

Note that temporal samples can easily be collected at a high sampling rate compared to the spatial sampling, where each sample requires a separate sensor. It is therefore beneficial to keep $n$ as small as possible by using tools that overcome the limitations of classical Shannon sampling theory. Consequently, full sampling is costly and time-consuming, and strategies for reducing the number of detector locations are desirable. In this study, we use $n = 64$ samples, which does not satisfy the Nyqvist criteria for the targeted discretization. However, the image quality in this case is still reasonable. To further reduce the number of measurements while preserving image quality, we use CS techniques.

2.2.

Compressed Sensing

The traditional approach to signal and image processing is to first collect a large number of point-like samples, which are then compressed and transmitted with minimal information loss. The basic idea of CS is to combine signal acquisition and compression by using specific indirect measurements together with mathematical algorithms that exploit the inherent structure of the image. In this way, a high-quality image can be recovered from a smaller number of measurements than required for point sampling at the same resolution. In particular, the seminal works¹³^,¹⁴ invented a theory of CS based on the sparsity of the signal to be recovered and the randomness of the measurements. Subsequent research has identified properties of the measurement matrix, such as the restricted isometry property (RIP), as key elements for stable and robust recovery.

The first basic ingredient of CS is sparsity, that is defined as follows. Let $s \in N$ and $x \in R^{n}$ . The vector $x$ is called $s$ -sparse, if ${‖ x ‖}_{0} ≔ # ({i \in {1, \dots, n} | x [i] \neq 0}) \leq s$ where $# (S)$ stands for the number of elements in a set $S$ . Signals of practical interest are often not sparse in the strict sense, but can be well approximated by sparse vectors. One calls $σ_{s} (x) ≔ \inf {{‖ x - x_{s} ‖}_{1} | x_{s} \in R^{n} is s - sparse}$ the best $s$ -term approximation error of $x \in R^{n}$ and calls $x$ compressible, if $σ_{s} (x)$ decays sufficiently fast with increasing $s$ .

2.2.1.

Restricted isometry constant

Let $s \in N$ and $δ \in (0,1)$ . Stable and robust recovery of sparse vectors requires the measurement matrix to well separate sparse vectors. The RIP guarantees such a separation. We recall that the measurement matrix $A \in R^{m \times n}$ is said to satisfy the RIP of order $s$ with constant $δ$ if

Eq. (5)

(1 - δ) {‖ x ‖}_{2}^{2} \leq {‖ Ax ‖}_{2}^{2} \leq (1 + δ) {‖ x ‖}_{2}^{2} for all s - sparse x \in R^{n},

and write

δ_{s}

for the smallest constant satisfying Eq. (5). Many sparse recovery results have been derived using the RIP. For example, the result derived in Ref. 15 states that if

A \in R^{m \times n}

satisfies the

2 s

-RIP with constant

δ_{2 s} < 1 / 2

then for

{‖ y - Ax ‖}_{2} \leq δ

any

x_{⋆} \in \arg \min {{‖ z ‖}_{1} | {‖ Az - y ‖}_{2}}

satisfies

{‖ x - x_{⋆} ‖}_{2} \leq c_{1} σ_{s} (x) / \sqrt{s} + c_{2} δ

for constants

c_{1}, c_{2}

depending only on

δ_{2 s}

. This implies stable and robust recovery for measurement matrices satisfying the RIP. The error estimate consists of two terms: The term

c_{2} ϵ

is due to the data noise and

c_{1} σ_{s} (x) / \sqrt{s}

accounts for the fact that the unknown may not be strictly

s

-sparse.

No deterministic construction is known providing large measurement matrices satisfying the RIP with near-optimal $s$ . However, several types of random matrices are known to satisfy the RIP with high probability. An important example of a random matrix that satisfies the RIP is the Bernoulli matrix, which is a random matrix $B \in {- 1,1}^{m \times n}$ having independent entries that take the values $- 1$ and 1 with equal probability. A Bernoulli matrix satisfies $δ_{2 s} < δ$ with probability tending to 1 as $m \to \infty$ , provided that $m \geq C_{δ} s (\log (n / s) + 1)$ for some constant $C_{δ} > 0$ as $n \to \infty$ . However, such a theory is hardly applicable in our situation due to the small dimension of our measurement matrices.

2.2.2.

Binary CS matrices

Another useful type of CS matrices is binary matrices having entries 0 or 1. Such measurement matrices can be interpreted as the adjacency matrix of a bipartite graph $(L, R, E)$ where $L ≔ {1, \dots, n}$ is the set of left vertices, $R ≔ {1, \dots, m}$ the set of right vertices, and $E \subseteq L \times R$ is the set of edges. Any element $(j, i) \in E$ can be interpreted as an edge joining vertices $j$ and $i$ . The left vertices $L$ represent the sensors, and the right vertices $R$ model each measurement. The vertex $j \in L$ is connected to the vertex $i \in R$ if sensor $j$ contributes to measurement $i$ . For our application, we have this type of binary measurement matrices.

Specific binary measurements are lossless expanders for which a stable and robust recovery theory exists.¹⁶^,¹⁷ However, these results are again asymptotic and are not applicable for PAPI with small CS matrices.

2.3.

All-Optical PA Projection Tomograph

In order to realize photoacoustic projection tomography, one needs one or several ILDs that integrate the pressure along one dimension. Initial setups used a single line detector that is moved around the object either using a free-beam Mach–Zehnder interferometer⁹ or a free-beam as well as fiber-based Fabry–Perot interferometer.¹⁸ To accelerate the data collection process arrays of line detectors have been developed either consisting of a piezoelectric array¹⁹ or an array of FOMZIs introduced in Refs. 8 and 20. Optical and piezoelectric ILDs have been compared in Ref. 21. A method where a PA projection image is collected at one shot is the full-field technique.²² In this paper, we use the FOMZI array reviewed below.

The PAPI setup consists of 18 individually designed (CAD) parts, for a total of 750 mechanical components. The fiber cage of the system is built with 64 GIPOFs, and each GIPOF has two end faces/ferrules, and five glue points, making a total of 128 end faces and 320 glue points. The fiber laser used is an NKT Koheras AdjustiK E15 with a maximum power of 200 mW and a line width of 0.1 kHz. A 1:2 fiber splitter directly after the fiber laser splits the optical path into a reference arm with 20% laser power and a measurement arm with 80% laser power. The 80/20 splitting is used because the measurement arm is split into 64 beams using a 1:64 fiber splitter whereas the reference arm is only split into 16 beams. Thus, each of the 80 fibers receives 1.25% of the overall laser power. The measurement arm consists of 64 GIPOFs arranged in a circular configuration and multiplexed with sixteen $4 \times 1$ fast fiber optic switches from Sercalo. The 16 fiber optic switches are controlled by the measurement software.

For working point stabilization of the FOMZIs, 16-fiber phase switches are integrated on four controller cards. A robust analog (bang-bang) controller with digital potentiometers and easy USB control was developed at RECENDT.²⁰ The reference and measurement arms are connected by sixteen 2:2 50/50 fiber couplers and the 16 self-developed balanced photodetectors detect the optical signal and provide two electrical signals. A low-frequency (LF) signal is employed for working point stabilization, while a high-frequency (HF) signal represents the actual data. The 16 PA signals are sampled by a National Instrument (NI) device with two cards, each with eight channels resulting in 16 channels in total. Each card has a maximum sampling rate of $60 MS s^{- 1}$ , 12 bit depth, and 128 MB on-board memory. The whole system is controlled by a PC with our own control and measurement software (NI LabWindows).

3. System Design, Implementation, and Analysis

In this section, we present details on the design, implementation, and analysis of our self-developed CS-PAPI device. It is built upon an extension of the all-optical PAPI described in Sec. 2.3 using specific CS measurements that we optimize by introducing the sparse injectivity number (SIN) as a quality measure for CS measurement matrices.

3.1.

Compressive PAPI

We conduct CS measurements of the pressure $P = W u$ in the detector domain, ensuring that pressures from different times are not mixed. Thus, instead of collecting $m$ individually sampled signals as in Eq. (4), we take CS measurements $y_{i, τ} = (A P)_{i, τ} ≜ \sum_{k = 1}^{n} a_{i, j} p_{j, ℓ}$ for $(i, ℓ) \in {1, \dots, m} \times {1, \dots, q}$ with $m < n$ . Recall that $n$ is the number of sensors, $m$ the number of measurements, and $q$ the number of temporal samples. If we write $Wx = [(W_{1} x)^{T}, \dots, (W_{n} x)^{T}]^{T}$ as a block column vector where the $j$ ’th row is the signal of the $j$ ’th ILD, the CS-PAPI data can be written as

Eq. (6)

Y = [\begin{matrix} y_{1} \\ ⋮ \\ y_{m} \end{matrix}] = A [\begin{matrix} W_{1} u \\ ⋮ \\ W_{n} u \end{matrix}] = AW u,

where

y_{i}

is

i

’th CS measurement signal.

The aim of CS-PAPI image reconstruction is to recover the unknown $u$ from data in Eq. (6). If the matrix $A$ would have rank $n$ , then Eq. (6) would have the solution $u = W^{#} [(A^{T} A)^{- 1} A^{T} Y]$ , where $W^{#}$ is a numerical realization of the inversion formula of the wave equation and $(A^{T} A)^{- 1} A^{T}$ is the least square inverse of $A$ . In the case of compressive measurements, however, we have $m < n$ and the matrix $A^{T} A$ is singular. Thus, solving $Y = AW u$ becomes underdetermined and reconstruction algorithms using specific prior information are required. Following the prime CS strategy, we use sparsity for that purpose.

Several choices for the CS measurement matrix $A$ have been suggested for PAT.²³^–²⁵ Specifically, for PAPI with ILDs binary CS matrices are often most easily realized in practice. In this case sparsifying transformations in the detector domain may negatively affect stable recovery results. Note that the CS measurement matrix $A$ in Eq. (6) does not act in the temporal variable. Thus, for any operation $Φ$ acting in temporal variable only, we have the commutation relation $A \circ Φ = Φ \circ A$ . This has been the motivation for the two-step image reconstruction approach proposed in Ref. 23, based on sparsifying temporal transforms, which we essentially follow here. However, in contrast to that paper, we use a structured CS measurement matrix where only certain sensor combinations are allowed to be guided by the experimental design.

3.2.

Proposed Structured CS Measurements

Recall that the PAPI system (see Fig. 1) consists of 64 ILDs in total, which naturally come in 16 blocks of four sensors each, where each of these blocks is characterized by sensors being connected to the same switch. We form CS measurements by selecting at most one sensor of each block and summing the signals over four neighboring blocks. In that way we make four CS measurements in parallel where the first measurement uses detectors in group $[1] = {1, \dots, 16}$ , the second in group $[2] = {16, \dots, 32}$ , the third in group $[3] = {33, \dots, 48}$ , and the fourth in group $[4] = {49, \dots, 64}$ . In every measurement, there is at most one ILD active within one block and every other sensor is inactive. Making $m_{0}$ such measurements, results for each group in a binary $m_{0} \times 16$ matrix

Eq. (7)

A_{G} = [A_{G, 1} | A_{G, 2} | A_{G, 3} | A_{G, 4}] \in {0,1}^{1 \times 16} for G = [1], [2], [3], [4],

where each block

A_{G, b}

has at most one non-vanishing entry. Entry 1 means the corresponding sensor is active and 0 means that the sensor is inactive. An example for such a matrix with

m_{0} = 2

measurements is

A_{G} = [\begin{matrix} \begin{matrix} 1 & 0 & 0 & 0 \end{matrix} & \begin{matrix} 0 & 0 & 0 & 0 \end{matrix} & \begin{matrix} 1 & 0 & 0 & 0 \end{matrix} & \begin{matrix} 0 & 1 & 0 & 0 \end{matrix} \\ \begin{matrix} 0 & 0 & 0 & 1 \end{matrix} & \begin{matrix} 1 & 0 & 0 & 0 \end{matrix} & \begin{matrix} 0 & 0 & 0 & 1 \end{matrix} & \begin{matrix} 0 & 0 & 0 & 0 \end{matrix} \end{matrix}] .

According to the general construction, each row is characterized to have at most one non-vanishing entry in each of the four blocks and the number of rows corresponds to the number of measurements for any group $G = [1], [2], [3], [4]$ .

The overall CS matrix acting on the 64 sensors arranged in four groups takes the block diagonal form

Eq. (8)

A = [\begin{matrix} A_{1} & 0 & 0 & 0 \\ 0 & A_{2} & 0 & 0 \\ 0 & 0 & A_{3} & 0 \\ 0 & 0 & 0 & A_{4} \end{matrix}] \in {0,1}^{4 m_{0} \times 64},

where

A_{G} \in {0,1}^{m_{0} \times 16}

has the structure as in Eq. (7). For these types of CS measurements combined with the sparsity paradigm, we address both the unique recovery question and the optimal design question. All matrices of the form Eqs. (7) and (8) are experimentally implementable with the CS-PAPI system.

Remark 1

(Selection of block size and group size) The parameters guiding the types of CS measurements are the block size (sensors having the same switch) and the number of blocks per group. The product of these numbers gives the group size. The specific choices are determined by the current PAPI setting (block size four and four groups per block); however, they can be adjusted according to different experimental designs. For example, by fixing the group size to 16, another choice is a block size of two and eight groups per block. Such measurements are found to improve CS capabilities. However, on the downside, this doubles the number of fiber phase switches. Our framework is completely flexible in terms of group number and block size. The concrete choice should be determined by practical considerations.

3.3.

Experimental Realization

In order to technically implement CS on PAPI, a plug-and-play concept was developed by designing and implementing a CS module named SUM4 (for summing over 4) that can be integrated into the PAPI system. Recall that before AD conversion, PAPI has 16 acoustic signals, where each signal corresponds to the ILD selected in the 16 blocks by the switch. As a first step, we extend PAPI by enabling the arbitrary selection of ILDs within each group. In addition, we construct SUM4, where signals from four neighboring blocks are summed, resulting in four electrical signals that are sampled by the NI card. Before summation, each signal can be potentially be switched off, resulting in CS measurements of the form (7) and (8). Figure 3 shows the schematic concept of SUM4, consisting of on/off switches, summation over blocks of four, and transmission to the ADC. In addition, Fig. 4 shows a photo of parts of the CS-PAPI system.

Fig. 3

Schematics of the PAPI system setup using 64 detector positions.

Fig. 4

(a) Illustration of the CS detection principle and its technical realization by SUM4.(b) Photo of the CS-PAPI system. Note that the CS module is located in the lower right part of Fig. 3, before the A/D converter. In addition, the 4:1 switches are modified to allow for variable ILD selection.

SUM4 can be seen as a device for analog signal conditioning and implements the CS aspects in the analog electrical domain. It allows the arbitrary superposition of up to four analog signals by switchable addition of the input signals. In addition, the design permits the compensation of system-related losses in signal amplitudes, such as those caused by impedance matching. The low-noise design of the analog signal paths results in a signal-to-noise ratio of 80 dBV, corresponding to a resolution of at least 13 bits. The selection of electrical signals to be superimposed is done via the USB port. This involves implementing a virtual COM port with a custom control protocol. This intuitive control protocol facilitates easy integration of the device into a larger network of instruments via USB. To ensure optimal integration with PAPI, the quad summers combine four separate summing groups in one device, allowing 16 input signals to be routed to four independent outputs.

3.4.

Optimal Design

The CS-PAPI system with SUM4 allows us to perform any CS measurements of the form (7) and (8). The aim in this section is to present a strategy for selecting optimal measurements within this class based on exact reconstruction. For that purpose, we first note that the measurements between the subgroups are independent and thus we aim for optimal design of each $m_{0} \times 16$ sub-matrix of the form (7). Second, we focus on optimal design in the context of sparse recovery. Thus, we aim for binary matrices $M \in R^{m_{0} \times 16}$ of the form (7) with $m_{0} < 16$ which allow us to recover sparse signals $x \in R^{16 \times 1}$ from data $Mx$ . Because the signal size is small, selecting these matrices at random (as in standard CS) resulted in matrices not enabling sparse recovery. We therefore designed a quality measure and a strategy to construct matrices enabling sparse recovery.

A minimal requirement for the identifiability of sparse elements $x \in R^{16}$ is the injectivity of $M$ over the set of $s$ -sparse elements. However, injectivity alone is not sufficient in the sense that ${Mx}_{1}$ and ${Mx}_{2}$ can get close to each other for sparse signals $x_{1}, x_{2}$ very different from each other. Thus, we actually need to bound the difference ${‖ {Mx}_{1} - {Mx}_{2} ‖}_{2}$ in order to sufficiently separate $x_{1}$ and $x_{2}$ . While this is essentially also included in the RIP constant, in this paper, we introduce a different concept that we think fits our aims better.

[SIN] For a matrix $M \in R^{m_{0} \times n_{0}}$ and any $s$ we define the $s$ -sparse injectivity number ( $s$ -SIN) of $M$ as

Eq. (9)

Θ_{s} (M) ≔ \inf {\frac{‖ {Mx}_{1} - {Mx}_{2} ‖_{2}}{‖ x_{1} - x_{2} ‖_{2}} | x_{1} \neq x_{2} \in R^{n_{0}} are s - sparse} .

Alternatively, the $s$ -SIN can be defined as the largest constant $Θ \geq 0$ such that $‖ {Mx}_{1} - {Mx}_{2} ‖_{2} \geq Θ ‖ x_{1} - x_{2} ‖_{2}$ for all $s$ -sparse signals $x_{1}, x_{2} \in R^{n_{0}}$ .

The $s$ -SIN is strictly positive if and only if the matrix $M$ is injective on the set of all $s$ -sparse elements. Unlike the usual RIP, it only asks for the one-sided estimate $‖ {Mx}_{1} - {Mx}_{2} ‖_{2} \geq Θ ‖ x_{1} - x_{2} ‖_{2}$ . Furthermore, for $s \leq n_{0} / 2$ , it is easy to verify that $σ_{s}$ is the smallest singular value among all $m_{0} \times 2 s$ sub-matrices of $M_{0}$ .

A good CS matrix is a CS matrix with $Θ_{s} (M)$ large relative to $‖ M ‖$ . Values of $Θ_{s} (M)$ greater than 0.1 have been empirically observed to result in stable and robust signal reconstruction. Randomly selecting $M$ from our class of matrices turned out to very often yield (almost) vanishing $s$ -SIN. On the other hand, computing the $s$ -SIN for all admissible matrices to make an optimal selection is computationally infeasible. Therefore, to determine a suitable CS matrix, we use a simple algorithm where we repeatedly randomly select $M$ from our CS matrix class and update the matrix whenever the $s$ -SIN is increased. This procedure is summarized in Algorithm 1, where for PAPI we have $n_{0} = 16$ .

Algorithm 1

Optimized detector selection for CS matrix with large s-SIM

1: SINopt ← 0

2: LISTopt ← zeros (1,4)

3: Mopt ← zeros

(m_{0}, n_{0})

4: fori in [1, Niter] do

5: LIST ← random.sample

(m_{0}, n_{0})

▹ Draw

m_{0}

lists of active sensors

6: M ← makeCSMatrix(LIST) ▹ Build the CS matrix

7: SIN ← getSIN(M,s) ▹ Compute the s-SIN of M

8: ifSIN > SINoptthen

9: LISTopt ← LIST

10: Mopt ← M

11: SINopt ← SIN

12: end if

13: end for

ReturnSINopt, LISTopt, Mopt ▹ Return optimal CS list, matrix and SIN

In Algorithm 1, the function random.sample selects a feasible list of sensors and the function makeCSMatrix forms the corresponding CS matrix. Furthermore, getSIN computes the $s$ -SIN. We have found empirically that the procedure results in CS matrices with an SIN over 0.1 in a reasonable time. Specifically, we take $m_{0} = 12$ and $s = 2$ for the results shown below.

Algorithm 1 can be extended to use block sizes other than four and numbers of blocks other than four. The only limiting factor is the increasing numerical complexity with increasing dimensions.

3.5.

Two-Step CS Image Reconstruction

Due to the separable nature of the image reconstruction problem (6), there are naturally two types of reconstruction methods, namely one-step image reconstruction and two-step image reconstruction. In the two-step methods, the complete data $W u$ are first recovered from CS data $A [W u]$ via iterative methods, and in a second step $u$ is recovered from $W u$ via wave inversion such as the FBP inversion formula. In the one-step approach, the initial pressure is directly recovered from CS data using iterative methods applied with the full forward operator $AW$ . Both classes of methods come with certain strengths and limitations. The two-step approach is fast as iterative signal reconstruction and is separated from the computationally costly evaluation of $W$ and its adjoint. Moreover, CS properties of the matrix $A$ can be exploited together with the sparsity of $W u$ , potentially after a suitable basis transform. On the downside, the image structure of $u$ cannot be directly integrated in the image reconstruction. The one-step approach, on the other hand, allows for easy integration of prior information about the image to be generated. However, CS reconstruction theory based on sparsity and specific properties of the forward matrix can hardly be integrated. Hybrid methods such as those proposed in Ref. 26 might overcome such issues. Another drawback of one-step approaches is that they necessitate the repeated use of the time-consuming evaluation of $W$ and its adjoint.

Due to its clear interpretability and computational efficiency in this study, we work with the two-step approach. Specifically, we utilize temporal transforms in combination with 1D total variation (TV) minimization. For that purpose, we apply a transform $Φ : R^{q} \to R^{q}$ acting in the time domain such that the transformed pressure $P Φ^{T}$ has sparse gradients. Thus, an approximation $H = [h_{1}^{T}, \dots, h_{n}^{T}]^{T}$ to $P Φ^{T}$ can be recovered by TV minimization:

Eq. (10)

{‖ A H - Y Φ^{T} ‖}^{2} + {‖ \partial_{1} H ‖}_{1} = \sum_{ℓ = 1}^{q} {‖ A h_{ℓ} - (Y Φ^{T})_{ℓ} ‖}^{2} + {‖ \partial_{1} h_{ℓ} ‖}_{1} \to \min_{H},

where

\partial_{1}

is the derivative in the sensor direction. Problem (10) can be solved by a series of 1D TV minimization problems for the 1D signals

h_{ℓ}

and is numerically efficient. Further, by writing the FBP Eq. (3) as

Eq. (11)

u (r) = - \frac{1}{π R} \int_{\partial B_{R}} \int_{| r - z |}^{\infty} \frac{(\partial_{t} t [Φ^{- 1} \circ Φ \circ W u]) (s, t)}{\sqrt{t^{2} - | r - s |^{2}}} d t d S (s),

we can recover the unknown

u

from the filtered signals

Φ \circ W u

in the first step. Equations (10) and (11) constitute the two-step method we use for image reconstruction in this paper.

Remark 2

Let us mention some further work on image reconstruction in CSPAT. Using intertwining relations between spatial and temporal operations for the wave equation, we extended the sparsifying transform approach to the image domain,²⁷^,²⁸ enabling one-step inversion. This and the two-step method can also be applied to CSPAT with standard point-like measurements. Other early work on CSPAT has been done in Refs. 29 30.31.32.–33, where various compressive sampling strategies have been used with sparse recovery techniques. Recently, machine learning methods have been used in the context of CSPAT.³⁴^–⁴⁰

4. Numerical Experiments

Due to the restricted CS matrices, it is challenging to achieve even a small compression factor $n / m$ . Note that for our structured CS matrices, we require sparsity within the four groups of 16 sensors each. For the following numerical investigation, we use a sparsity level of $s = 2$ . Numerically, it turns out that we need 12 measurements to obtain a non-singular SIN with the algorithm outlined above.

4.1.

Measurement Design Results

We use the parameters of our PAPI system, where measurements on a group of detectors have the form (7), whose structure is determined by the size of the blocks (which is four for our PAPI) and the number of blocks within a group (which is also four for our PAPI). The goal of CS is to keep the number $m_{0}$ of measurements small, while allowing the unique recovery of certain elements. Following the sparsity paradigm, our approach is to use Algorithm 1 to find a matrix with non-singular $s$ -SIN. The larger $s$ , the more general the signal class, but the less likely it is to get a non-vanishing $s$ -SIN. So we take $2 s = 4$ to have at least some generality in the signal class.

Fig. 5

Reconstruction results from exact data. (a) (from left to right): Data from 64 ILD, reconstruction using an optimized CS matrix and reconstruction using a random matrix. (b) Corresponding time-transformed data from the pressure. (c) Corresponding FBB reconstructions. In the first two rows, the horizontal direction represents the spatial dimension, while the vertical direction represents time.

Running Algorithm 1 we found that a non-singular SIN could be found for $m_{0} = 12$ measurements. In particular, in almost every test run with 100 iterations, we could find a matrix with an SIN of about 0.14, which we then selected. Even for $m = 11$ we could find such matrices after a longer search. However, we could not increase the compression factor further in the sense that for $m = 10$ , even after 100,000 iterations, no SIN larger than machine precision could be found. Roughly speaking, our work demonstrates a compression factor of at least $4 / 3$ for block size 4 and group size 16.

Remark 3

(Variable block size and group size) In order to put our work into a broader perspective, it is worth investigating whether different block sizes and numbers of blocks result in a larger compression factor. Testing our algorithm with the same group size but a block size of two, we found that indeed, using $m = 10$ measurements results in a nonsingular s-SIN of $\sim 0.21$ , demonstrating an increased compression factor of $8 / 5$ . A similar effect has been observed when keeping the block size constant while increasing the group size.

Having a non-singular 2-SIR allows for theoretical exact recovery of two-sparse signals from exact data. In reality, robustness regarding noise and stability concerning the sparsity level using specific reconstruction algorithms are central. While this is not part of our theory, we expect similar results to the (unfortunately asymptotic) theory of CS. Our numerical results below support this.

4.2.

Image Reconstruction Results

For image reconstruction, we use the two-step sparse recovery method described above. The key there is to apply a temporal transform to obtain sparsity. Here, we use a phantom such that the spherical means are piecewise constant. Thus, in the first step, we use the Abel transform as the time transform and recover the spherical means using TV minimization Eq. (10).

Reconstruction results from exact and noisy data are shown in Figs. 5 and 6. We use two different measurement matrices, the first one is found by our algorithm and the second one is a randomly selected matrix from the CSPAT family that we corrected by educated guess to get non-vanishing 2-SIN (see Fig. 7). The CS measurement data and the added noise are shown in Fig. 8. For specific parameter settings, we refer to the MATLAB code that will be made publicly available. We consider the FBP reconstruction as our ground truth because our aim is to approximate the image quality achieved with the full sensor array (64 sensors). Our ground truth phantom consists of circles, but they are not homogeneous. The profile has been chosen such that the spherical means of the circular regions are piecewise constant, making it well-suited for TV minimization. In this way, we avoid a transformation that modifies the signal in that regard, as suggested in Ref. 23.

Fig. 6

Reconstruction results from noisy data. (a) (from left to right): Data from 64 ILD, reconstruction using an optimized CS matrix and reconstruction using a random matrix. (b) Corresponding time-transformed data from the pressure. (c) Corresponding FBB reconstructions. In the first two rows, the horizontal direction represents the spatial dimension, while the vertical direction represents time.

Fig. 7

(a) Comparison of a measurement matrix for optimized 2-SIN and (b) a random matrix according to the CS setup, which has been corrected to have non-vanishing SIN. (c) The computed $s$ -SIN for the sparsity level $2 s = 1,2, 3,4, 5$ .

Fig. 8

Exact and noisy data for the optimized and the random matrix. The horizontal direction represents the spatial dimension, while the vertical direction represents time.

We find that the reconstruction procedure is indeed very stable and robust. In particular, the noise had a small negative impact on the results. The reconstruction artifacts are due to the failure of the strict 2-sparsity assumption. To support such a claim, we also show results (Fig. 9) for a simple phantom where 2-sparsity on the 16-groups almost holds. In this case, the CS reconstruction hardly differs from the ground truth. For precise relative error values, see Table 1. All reconstruction results demonstrate stability and robustness.

Fig. 9

Reconstruction results of sparse object. (a) (from left to right): Data from 64 ILD, reconstruction using an optimized CS matrix and reconstruction using a random matrix. (b) Corresponding time-transformed data from the pressure. (c) Corresponding FBB reconstructions. In the first two rows, the horizontal direction represents the spatial dimension, while the vertical direction represents time.

Table 1

Relative ℓ2 error in the data (row 1), the CS reconstruction (row 2), and the final FBP reconstruction error using the optimized and random matrix (row 3).

	Optimized A			Random A
	Data	CS	FBP	Data	CS	FBP
Non-sparse phantom (noisy)	0.0901	0.3159	0.6686	0.0766	0.3414	0.6810
Non-sparse phantom (exact)	x	0.2594	0.6049	x	0.2849	0.6571
Sparse phantom (exact)	x	0.0003	0.0010	x	0.0107	0.0509

5. Conclusion and Outlook

In this paper, we presented the experimental realization of a CS-PAPI system extending the existing tomograph. We demonstrated that the specific setup allows perfect recovery of sparse signals. However, for that purpose, we could not select an admissible matrix uniformly at random, but a systematic strategy exploiting the SIM.

One future task is to go beyond the sparsity model. Thus our aim is to find CS matrices $A \in R^{m \times 16}$ not targeting sparsity but actual real data. This can be done two-fold. First one can train a matrix such that $16 \times 1$ pieces in data domain are optimally separated. Second optimization can be improved by optimizing over the image space. This allows us to consider that, due to the forward map $W$ , the $16 \times 1$ patches are actually correlated since they originate from the same initial source. Deep learning and neural networks are natural candidates unveiling such hidden correlation.

Disclosures

None of the authors have any potential conflicts of interest to disclose.

Code and Data Availability

The code for generating a CS matrix with large SIN as well for producing the shown numerical results are available upon request. No additional data are required for this study.

Acknowledgments

This work has been supported by the Austrian Science Fund (FWF), projects P 30747-N32 and P 33019-N.

References

1.

M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 https://doi.org/10.1063/1.2195024 RSINAK 0034-6748 (2006). Google Scholar

2.

K. Wang, M. A. Anastasio, “Photoacoustic and thermoacoustic tomography: image formation principles,” Handbook of Mathematical Methods in Imaging, Springer, New York (2011). Google Scholar

3.

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

4.

J. Poudel, Y. Lou and M. A. Anastasio, “A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography,” Phys. Med. Biol., 64 (14), 14TR01 https://doi.org/10.1088/1361-6560/ab2017 PHMBA7 0031-9155 (2019). Google Scholar

5.

A. Rosenthal, V. Ntziachristos and D. Razansky, “Acoustic inversion in optoacoustic tomography: a review,” Curr. Med. Imaging, 9 (4), 318 –336 https://doi.org/10.2174/15734056113096660006 (2013). Google Scholar

6.

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

7.

G. Paltauf et al., “Photoacoustic tomography with integrating area and line detectors,” Photoacoustic Imaging and Spectroscopy, 251 –264 CRC Press( (2017). Google Scholar

8.

J. Bauer-Marschallinger, K. Felbermayer and T. Berer, “All-optical photoacoustic projection imaging,” Biomed. Opt. Express, 8 (9), 3938 –3951 https://doi.org/10.1364/BOE.8.003938 BOEICL 2156-7085 (2017). Google Scholar

9.

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

10.

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

11.

D. Finch, M. Haltmeier and Rakesh, “Inversion of spherical means and the wave equation in even dimensions,” SIAM J. Appl. Math., 68 (2), 392 –412 https://doi.org/10.1137/070682137 SMJMAP 0036-1399 (2007). Google Scholar

12.

M. Haltmeier, “Sampling conditions for the circular radon transform,” IEEE Trans. Image Process., 25 (6), 2910 –2919 https://doi.org/10.1109/TIP.2016.2551364 IIPRE4 1057-7149 (2016). Google Scholar

13.

D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, 52 (4), 1289 –1306 https://doi.org/10.1109/TIT.2006.871582 IETTAW 0018-9448 (2006). Google Scholar

14.

E. J. Candès, J. Romberg and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, 52 (2), 489 –509 https://doi.org/10.1109/TIT.2005.862083 IETTAW 0018-9448 (2006). Google Scholar

15.

T. T. Cai and A. Zhang, “Sharp RIP bound for sparse signal and low-rank matrix recovery,” Appl. Comput. Harmon. Anal., 35 (1), 74 –93 https://doi.org/10.1016/j.acha.2012.07.010 ACOHE9 1063-5203 (2013). Google Scholar

16.

R. Berinde et al., “Combining geometry and combinatorics: a unified approach to sparse signal recovery,” in 46th Annu. Allerton Conf. on Commun., Control, and Comput., 798 –805 (2008). https://doi.org/10.1109/ALLERTON.2008.4797639 Google Scholar

17.

S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing, Springer( (2013). Google Scholar

18.

P. Burgholzer et al., “Thermoacoustic tomography using a fiber-based Fabry-Perot interferometer as an integrating line detector,” Proc. SPIE, 6086 60861N https://doi.org/10.1117/12.644047 PSISDG 0277-786X (2006). Google Scholar

19.

G. Paltauf et al., “Piezoelectric line detector array for photoacoustic tomography,” Photoacoustics, 8 28 –36 https://doi.org/10.1016/j.pacs.2017.09.002 (2017). Google Scholar

20.

J. Bauer-Marschallinger et al., “Photoacoustic projection imaging using a 64-channel fiber optic detector array,” Proc. SPIE, 9323 93233U https://doi.org/10.1117/12.2077206 PSISDG 0277-786X (2015). Google Scholar

21.

R. Nuster and G. Paltauf, “Comparison of piezoelectric and optical projection imaging for three-dimensional in vivo photoacoustic tomography,” J. Imaging, 5 (1), 15 https://doi.org/10.3390/jimaging5010015 (2019). Google Scholar

22.

R. Nuster et al., “Full field detection in photoacoustic tomography,” Opt. Express, 18 (6), 6288 –6299 https://doi.org/10.1364/OE.18.006288 OPEXFF 1094-4087 (2010). Google Scholar

23.

M. Sandbichler et al., “A novel compressed sensing scheme for photoacoustic tomography,” SIAM J. Appl. Math., 75 (6), 2475 –2494 https://doi.org/10.1137/141001408 SMJMAP 0036-1399 (2015). Google Scholar

24.

M. Haltmeier et al., “Compressed sensing and sparsity in photoacoustic tomography,” J. Opt., 18 (11), 114004 https://doi.org/10.1088/2040-8978/18/11/114004 (2016). Google Scholar

25.

M. M. Betcke et al., “Acoustic wave field reconstruction from compressed measurements with application in photoacoustic tomography,” IEEE Trans. Comput. Imaging, 3 710 –721 https://doi.org/10.1109/TCI.2017.2706029 (2017). Google Scholar

26.

A. Ebner and M. Haltmeier, “Convergence rates for the joint solution of inverse problems with compressed sensing data,” Inverse Probl., 39 (1), 015011 https://doi.org/10.1088/1361-6420/aca5ae INPEEY 0266-5611 (2022). Google Scholar

27.

M. Haltmeier et al., “A sparsification and reconstruction strategy for compressed sensing photoacoustic tomography,” J. Acoust. Soc. Am., 143 (6), 3838 –3848 https://doi.org/10.1121/1.5042230 JASMAN 0001-4966 (2018). Google Scholar

28.

G. Zangerl and M. Haltmeier, “Multiscale factorization of the wave equation with application to compressed sensing photoacoustic tomography,” SIAM J. Imaging Sci., 14 (2), 558 –579 https://doi.org/10.1137/20M1356154 (2021). Google Scholar

29.

J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imaging, 28 (4), 585 –594 https://doi.org/10.1109/TMI.2008.2007825 ITMID4 0278-0062 (2009). Google Scholar

30.

Z. Guo et al., “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt., 15 (2), 021311 https://doi.org/10.1117/1.3381187 JBOPFO 1083-3668 (2010). Google Scholar

31.

G. S. Alberti, P. Campodonico and M. Santacesaria, “Compressed sensing photoacoustic tomography reduces to compressed sensing for undersampled Fourier measurements,” SIAM J. Imaging Sci., 14 (3), 1039 –1077 https://doi.org/10.1137/20M1375152 (2021). Google Scholar

32.

S. Arridge et al., “Accelerated high-resolution photoacoustic tomography via compressed sensing,” Phys. Med. Biol., 61 (24), 8908 https://doi.org/10.1088/1361-6560/61/24/8908 PHMBA7 0031-9155 (2016). Google Scholar

33.

N. Huynh et al., “Single-pixel camera photoacoustic tomography,” J. Biomed. Opt., 24 (12), 121907 https://doi.org/10.1117/1.JBO.24.12.121907 JBOPFO 1083-3668 (2019). Google Scholar

34.

A. Hauptmann et al., “Model-based learning for accelerated, limited-view 3-D photoacoustic tomography,” IEEE Trans. Med. Imaging, 37 (6), 1382 –1393 https://doi.org/10.1109/TMI.2018.2820382 ITMID4 0278-0062 (2018). Google Scholar

35.

S. Antholzer et al., “NETT regularization for compressed sensing photoacoustic tomography,” Proc. SPIE, 10878 108783B https://doi.org/10.1117/12.2508486 PSISDG 0277-786X (2019). Google Scholar

36.

J. Gröhl et al., “Deep learning for biomedical photoacoustic imaging: a review,” Photoacoustics, 22 100241 https://doi.org/10.1016/j.pacs.2021.100241 (2021). Google Scholar

37.

A. Hauptmann and B. Cox, “Deep learning in photoacoustic tomography: current approaches and future directions,” J. Biomed. Opt., 25 (11), 112903 https://doi.org/10.1117/1.JBO.25.11.112903 JBOPFO 1083-3668 (2020). Google Scholar

38.

M. A. Anastasio, “Deep learning and photoacoustic image formation: promises and challenges,” Proc. SPIE, PC12379 PC1237901 https://doi.org/10.1117/12.2674089 PSISDG 0277-786X (2023). Google Scholar

39.

N. Davoudi, X. L. Deán-Ben and D. Razansky, “Deep learning optoacoustic tomography with sparse data,” Nat. Mach. Intell., 1 (10), 453 –460 https://doi.org/10.1038/s42256-019-0095-3 (2019). Google Scholar

40.

S. Antholzer, M. Haltmeier and J. Schwab, “Deep learning for photoacoustic tomography from sparse data,” Inverse Probl. Sci. Eng., 27 (7), 987 –1005 https://doi.org/10.1080/17415977.2018.1518444 (2019). Google Scholar

Biographies of the authors are not available.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Markus Haltmeier, Matthias Ye, Karoline Felbermayer, Florian Hinterleitner, and Peter Burgholzer "Design, implementation, and analysis of a compressed sensing photoacoustic projection imaging system," Journal of Biomedical Optics 29(S1), S11529 (22 April 2024). https://doi.org/10.1117/1.JBO.29.S1.S11529

Received: 24 October 2023; Accepted: 28 February 2024; Published: 22 April 2024

Access the abstract

JOURNAL ARTICLE
16 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

KEYWORDS

Matrices

Design

Sensors

Image restoration

Compressed sensing

Biological research

Imaging systems

Significance

Aim

Approach

Results

Conclusions

1.

Introduction

1.1.

Photoacoustic Projection Imaging

Fig. 1

1.2.

Compressed Sensing in PAPI

Eq. (1)

1.3.

Outline

2.

Background

2.1.

PA Projection Imaging

Fig. 2

Eq. (2)

Eq. (3)

Eq. (4)

2.2.

Compressed Sensing

2.2.1.

Restricted isometry constant

Eq. (5)

2.2.2.

Binary CS matrices

2.3.

All-Optical PA Projection Tomograph

3.

System Design, Implementation, and Analysis

3.1.

Compressive PAPI

Eq. (6)

3.2.

Proposed Structured CS Measurements

Eq. (7)

Eq. (8)

Remark 1

3.3.

Experimental Realization

Fig. 3

Fig. 4

3.4.

Optimal Design

Eq. (9)

Algorithm 1

3.5.

Two-Step CS Image Reconstruction

Eq. (10)

Eq. (11)

Remark 2

4.

Numerical Experiments

4.1.

Measurement Design Results

Fig. 5

Remark 3

4.2.

Image Reconstruction Results

Fig. 6

Fig. 7

Fig. 8

Fig. 9

Table 1

5.

Conclusion and Outlook

Disclosures

Code and Data Availability

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years