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

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.


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 [40,38,39].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 [34,36].Specifically, we focus on PA projection imaging (PAPI) based on integrating line detectors (ILDs) [12,33].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 [6].

Photoacoustic projection imaging (PAPI)
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 consists of samples of the linear projection of the 3D acoustic pressure wave in the direction of the ILDs.Reconstruction in 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 [30,10].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, 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 4 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.

Compressed sensing (CS) in PAPI
Following the CS paradigm, instead of recording pressure signals P = [p T 1 , . . ., p T n ] T where p j is the pressure signal (written as column vector) of the j-th ILD, we record CS data with A = (a i,j ) i,j ∈ R m×n denoting the CS measurement matrix.Usually in CS, the measurement matrix chosen randomly, since this gives 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 → ∞.
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.

Outline
In this paper, we present our finding 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 presents 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.

Background
In this section we present the background of our work.This includes PAPI modeling (subsection 2.1), sparse CS theory (subsection 2.2) and the description of the selfdeveloped PAPI system (subsection 2.3).

PA projection imaging (PAPI)
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 [30,10] and in this work we consider the 2D version only.Further, we restrict ourselves to constant sound speed and a circular measurement geometry as illustrated in Figure 2.
Let us denote by u : R 2 → 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 where δ ′ (t) is the first time derivative of the Dirac delta distribution, r ∈ R 2 is the spatial location, t ∈ R the time variable, ∆ r the spatial Laplacian and v s the constant speed of sound.The wave equation ( 2) is augmented with p(r, t) = 0 for t < 0 such that the acoustic pressure is uniquely defined as solution of (2).We rescale time in such a way that v s = 1.
PAPI in circular geometry consist in recovering the function u from measurements of Wu(s, t) = p(s, t) made on C R × (0, ∞).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 [18] u Note the inversion operator in (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 p (s j , t ℓ ) for (j, ℓ) ∈ {1, . . ., n} × {1, . . ., q} , with s j ≜ R(cos (Ω(j − 1)/n) , sin (Ω(j − 1)/n)), t ℓ ≜ 2R(ℓ − 1)/(q − 1), and Ω ≤ 2π denoting the angular covering on the detection circle.The number n of detector positions in ( 4) is directly related to the resolution of the final reconstruction.Namely, n ≥ 2Rλ equally spaced transducers covering the full circle are required to stably recover any PA source u that has maximal essential wavelength λ; see [22].Image reconstruction in this case can be performed by discretizing the inversion formula (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 2Rλ = πN r , where N r × N r is the number of samples for discretizing the object of interest on the square 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.

Compressed sensing (CS)
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 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 [16,14] 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 ∈ N and 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 , and calls x compressible, if σ s (x) decays sufficiently fast with increasing s.

The restricted isometry constant (RIP)
Let s ∈ N and δ ∈ (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 ∈ R m×n is said to satisfy the RIP of order s with constant δ if and write δ s for the smallest constant satisfying (5).Many sparse recovery results have been derived using the RIP.For example, the result derived in [13] states that if A ∈ R m×n satisfies the 2s-RIP with constant δ 2s < 1/2 then for ∥y − Ax∥ 2 ≤ δ any c 2 depending only on δ 2s .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)/ √ 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 ∈ {−1, 1} m×n having independent entries that take the values −1 and 1 with equal probability.A Bernoulli matrix satisfies δ 2s < δ with probability tending to 1 as m → ∞, provided that m ≥ C δ s(log(n/s) + 1) for some constant C δ > 0 as n → ∞.However, such a theory is hardly applicable in our situation due to the small dimension of our measurement matrices.

Binary CS matrices
Another useful type of CS matrices are 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, . . ., n is the set of left vertices, R := {1, . . ., m} the set of right vertices, and E ⊆ L × R is the set of edges.Any element (j, i) ∈ 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 ∈ L is connected to the vertex i ∈ 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 [8,19].However, these results are again of asymptotic nature and are not applicable for PAPI with small CS matrices.

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 [32] or a free-beam as well as fiber-based Fabry-Perot interferometer, [11].To accelerate the data collection process arrays of line detectors have been developed either consisting of a piezoelectric array [31] or an array of FOMZIs introduced in [7,6].Optical and piezoelectric ILDs have been compared in [28].A method where a PA projection image is collected at one shot is the full-field technique [29].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 Figure 3: Schematics of the photoacoustic projection imaging (PAPI) system setup using 64 detector positions.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.1kHz.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 × 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 4 controller cards.A robust analog (bang-bang) controller with digital potentiometers and easy USB control was developed at RECENDT [7].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 2 cards, each with 8 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).

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 build upon an extension of the all-optical PAPI described in Section 2.3 using specific CS measurements that we optimize by introducing the sparse injectivity number (SIN) as a quality measure for CS measurement matrices.

Compressive PAPI
We conduct CS measurements of the pressure P = Wu in the detector domain, ensuring that pressures from different times are not mixed.Thus, instead of collecting m individually sampled signals as in (4), we take CS measurements y i,τ = (AP ) i,τ ≜ n k=1 a i,j p j,ℓ for (i, ℓ) ∈ {1, . . ., m} × {1, . . ., 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 , . . ., (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 where the y i is i-th CS measurement signal.
The aim of CS-PAPI image reconstruction is to recover the unknown u from data in (6).If the matrix A would have Rank n, then (6) would have the solution 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 = AWu 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 [37,23,9].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 (6) does not act in the temporal variable.Thus for any operation Φ acting in temporal variable only, we have the commutation relation A • Φ = Φ •A.This has been the motivation for the two-step image reconstruction approach proposed in [37], 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.

Proposed structured CS measurements
Recall that the PAPI system (see Figure 1) consists of 64 ILDs in total wich naturally come in 16 blocks of 4 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 [ 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 = 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 .
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 where A G ∈ {0, 1} m 0 ×16 has the structure as in (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 ( 7) and ( 8) are experimentally implementable with the CS-PAPI system.
Remark 3.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.

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.Additionally, 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), (8). Figure 4 illustrates the schematic concept of SUM4, consisting of on/off switches, summation over blocks of four, and transmission to the ADC.Additionally, Figure 4 shows a photo of parts of the CS-PAPI system.
SUM4 can be seen as a device for analog signal conditioning and implements the CS Additionally, 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.

Optimal design
The CS-PAPI system with SUM4 allows us to perform any CS measurements of the form (7), (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 × 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 ∈ R m 0 ×16 of the form ( 7) with m 0 < 16 which allow us to recover sparse signals x ∈ R 16×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 ∈ 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 which we think fits better to our aims.

Definition 3.2 (Sparse injectivity number).
For a matrix M ∈ R m 0 ×n 0 and any s we define the s-sparse injectivity number (s-SIN) of M as Alternatively, the s-SIN can be defined as the largest constant Θ ≥ 0 such that 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 . Furthermore, for s ≤ n 0 /2, it is easy to verify that σ s is the smallest singular value among all m 0 × 2s 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.
In Algorithm 1, the function random.sampleselects 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 procedure results in CS matrices with a SIN over 0.1 in a reasonable time.Specifically, we take m 0 = 12 and s = 2 for the results shown below.if SIN > SINopt then 9: LISTopt ← LIST 10: end if 13: end for return SINopt, LISTopt, Mopt ▷ Return optimal CS list, matrix and SIN 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.

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 Wu are first recovered from CS data A[Wu] via iterative methods, and in a second step u is recovered from Wu 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, 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 Wu, potentially after suitable basis transform.On the downside, 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 he forward matrix can hardly be integrated.Hybrid methods such as those proposed in [17] 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 → R q acting in the time domain such that the transformed pressure P Φ T has sparse gradients.Thus an approximation H = [h T 1 , . . ., h T n ] T to P Φ T can be recovered by TV minimization where ∂ 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 formula (3) as we can recover the unknown u from the filtered signals Φ •Wu in the first step.Equation ( 10) and ( 11) constitute the two-step method we use for image reconstruction in this paper Remark 3.3.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 [24,41], 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 [35,21,1,5,27], where various compressive sampling strategies have been used with sparse recovery techniques.Recently, machine learning methods have been used in the context of CSPAT [25,4,20,26,2,15,3].

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 4 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.

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 4 for our PAPI) and the number of blocks within a group (which is also 4 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 nonvanishing s-SIN.So we take 2s = 4 to have at least some generality in the signal class.
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 a 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 100000 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 4.1 (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 approximately 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 2-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.

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 (10).
Reconstruction results from exact and noisy data are shown in Figure 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.The CS measurement data and the added noise are shown in Figure 8.For specific parameter settings, we refer to the Matlab code that will be made publicly available.We consider the FBP recon- struction 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 total variation (TV) minimization.In this way, we avoid a transformation that modifies the signal in that regard, as suggested in [37].
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 (Figure 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.

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 ∈ R m×16 not targeting sparsity but actual real data.This can be done two-fold.First one can train a matrix such that 16 × 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 × 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, Data, and Materials Availability
The code for generating a CS matrix with large SIN as well for producing the shown numerical results are provided upon request.No additional data is required for this study.

Figure 1 :
Figure 1: Photographic image of the PA projection tomography with 64 fiber optic Mach-Zehnder interferometers (FOMZIs) as ILDs forming the basics of the presented research.

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

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

Figure 5 :
Figure 5: Comparison of a measurement matrix for optimized 2-SIN (top left) and a random matrix according to the CS setup (top right), which has been corrected to have non-vanishing SIN.The bottom image shows the computed s-SIN for the sparsity level 2s = 1, 2, 3, 4, 5.

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

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

Figure 8 :
Figure 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.

Figure 9 :
Figure 9: Reconstruction results of sparse object.Top row (from left to right): Data from 64 ILD, reconstruction using an optimized CS matrix and reconstruction using a random matrix.Middle row: corresponding time-transformed data from the pressure.Bottom: 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).