1.IntroductionPhotoacoustic 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.1–3 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.4,5 Specifically, we focus on PA projection imaging (PAPI) based on integrating line detectors (ILDs).6,7 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.8 1.1.Photoacoustic Projection ImagingA 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.9,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 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. 1.2.Compressed Sensing in PAPIFollowing the CS paradigm, instead of recording pressure signals where is the pressure signal (written as column vector) of the ’th ILD, we record CS data where 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 . However, in practice, and specifically in our application, the matrix 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 are often restricted to specific values, in our case, for example to 0 and 1. Finally, the dimensionality in our case is small, which limits the applicability of existing asymptotic CS theory that applies to the limit .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 below 2 seems to be a substantial challenge. 1.3.OutlineIn 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.BackgroundIn 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 ImagingPA 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 problem9,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 shown in Fig. 2. Let us denote by the 2D PA source distribution which is our image of interest and supposed to be enclosed by a circle of radius . The 2D projected pressure satisfies the 2D wave equation where is the first time derivative of the Dirac delta distribution, is the spatial location, is the time variable, is the spatial Laplacian, and is the constant speed of sound. The wave Eq. (2) is augmented with for such that the acoustic pressure is uniquely defined as solution of Eq. (2). We rescale time in such a way that .PAPI in circular geometry consists of recovering the function from measurements of made on . In the case of full data, exact and stable PA image reconstruction is possible and several efficient methods for recovering are available. We will use the FBP formula derived in Ref. 11: Note the inversion operator in Eq. (3) is also the adjoint of the forward operator . This in particular implies that inverting 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 with , , and denoting the angular covering on the detection circle. The number of detector positions in Eq. (4) is directly related to the resolution of the final reconstruction. Namely, equally spaced transducers covering the full circle are required to stably recover any PA source 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 , where is the number of samples for discretizing the object of interest on the square . In this case, we get 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 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 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 SensingThe 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 works13,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 and . The vector is called -sparse, if where stands for the number of elements in a set . Signals of practical interest are often not sparse in the strict sense, but can be well approximated by sparse vectors. One calls the best -term approximation error of and calls compressible, if decays sufficiently fast with increasing . 2.2.1.Restricted isometry constantLet and . 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 is said to satisfy the RIP of order with constant if and write 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 satisfies the -RIP with constant then for any satisfies for constants depending only on . This implies stable and robust recovery for measurement matrices satisfying the RIP. The error estimate consists of two terms: The term is due to the data noise and accounts for the fact that the unknown may not be strictly -sparse.No deterministic construction is known providing large measurement matrices satisfying the RIP with near-optimal . 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 having independent entries that take the values and 1 with equal probability. A Bernoulli matrix satisfies with probability tending to 1 as , provided that for some constant as . However, such a theory is hardly applicable in our situation due to the small dimension of our measurement matrices. 2.2.2.Binary CS matricesAnother 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 where is the set of left vertices, the set of right vertices, and is the set of edges. Any element can be interpreted as an edge joining vertices and . The left vertices represent the sensors, and the right vertices model each measurement. The vertex is connected to the vertex if sensor contributes to measurement . 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.16,17 However, these results are again asymptotic and are not applicable for PAPI with small CS matrices. 2.3.All-Optical PA Projection TomographIn 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 interferometer9 or a free-beam as well as fiber-based Fabry–Perot interferometer.18 To accelerate the data collection process arrays of line detectors have been developed either consisting of a piezoelectric array19 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.22 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 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.20 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 , 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 AnalysisIn 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 PAPIWe conduct CS measurements of the pressure in the detector domain, ensuring that pressures from different times are not mixed. Thus, instead of collecting individually sampled signals as in Eq. (4), we take CS measurements for with . Recall that is the number of sensors, the number of measurements, and the number of temporal samples. If we write as a block column vector where the ’th row is the signal of the ’th ILD, the CS-PAPI data can be written as where is ’th CS measurement signal.The aim of CS-PAPI image reconstruction is to recover the unknown from data in Eq. (6). If the matrix would have rank , then Eq. (6) would have the solution , where is a numerical realization of the inversion formula of the wave equation and is the least square inverse of . In the case of compressive measurements, however, we have and the matrix is singular. Thus, solving 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 have been suggested for PAT.23–25 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 in Eq. (6) does not act in the temporal variable. Thus, for any operation acting in temporal variable only, we have the commutation relation . 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 MeasurementsRecall 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 , the second in group , the third in group , and the fourth in group . In every measurement, there is at most one ILD active within one block and every other sensor is inactive. Making such measurements, results for each group in a binary matrix where each block 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 measurements isAccording 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 . The overall CS matrix acting on the 64 sensors arranged in four groups takes the block diagonal form where 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 RealizationIn 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. 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 DesignThe 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 sub-matrix of the form (7). Second, we focus on optimal design in the context of sparse recovery. Thus, we aim for binary matrices of the form (7) with which allow us to recover sparse signals from data . 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 is the injectivity of over the set of -sparse elements. However, injectivity alone is not sufficient in the sense that and can get close to each other for sparse signals very different from each other. Thus, we actually need to bound the difference in order to sufficiently separate and . 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 and any we define the -sparse injectivity number (-SIN) of as Alternatively, the -SIN can be defined as the largest constant such that for all -sparse signals . The -SIN is strictly positive if and only if the matrix is injective on the set of all -sparse elements. Unlike the usual RIP, it only asks for the one-sided estimate . Furthermore, for , it is easy to verify that is the smallest singular value among all sub-matrices of . A good CS matrix is a CS matrix with large relative to . Values of greater than 0.1 have been empirically observed to result in stable and robust signal reconstruction. Randomly selecting from our class of matrices turned out to very often yield (almost) vanishing -SIN. On the other hand, computing the -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 from our CS matrix class and update the matrix whenever the -SIN is increased. This procedure is summarized in Algorithm 1, where for PAPI we have . Algorithm 1Optimized detector selection for CS matrix with large s-SIM
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 -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 and 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 ReconstructionDue 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 are first recovered from CS data via iterative methods, and in a second step is recovered from 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 . 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 and its adjoint. Moreover, CS properties of the matrix can be exploited together with the sparsity of , potentially after a suitable basis transform. On the downside, the image structure of 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 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 acting in the time domain such that the transformed pressure has sparse gradients. Thus, an approximation to can be recovered by TV minimization: where is the derivative in the sensor direction. Problem (10) can be solved by a series of 1D TV minimization problems for the 1D signals and is numerically efficient. Further, by writing the FBP Eq. (3) as we can recover the unknown from the filtered signals in the first step. Equations (10) and (11) constitute the two-step method we use for image reconstruction in this paper.Remark 2Let 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,27,28 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. 2930.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.34–40 4.Numerical ExperimentsDue to the restricted CS matrices, it is challenging to achieve even a small compression factor . 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 . Numerically, it turns out that we need 12 measurements to obtain a non-singular SIN with the algorithm outlined above. 4.1.Measurement Design ResultsWe 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 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 -SIN. The larger , the more general the signal class, but the less likely it is to get a non-vanishing -SIN. So we take to have at least some generality in the signal class. Running Algorithm 1 we found that a non-singular SIN could be found for 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 we could find such matrices after a longer search. However, we could not increase the compression factor further in the sense that for , 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 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 measurements results in a nonsingular s-SIN of , demonstrating an increased compression factor of . 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 ResultsFor 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. 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. Table 1Relative ℓ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).
5.Conclusion and OutlookIn 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 not targeting sparsity but actual real data. This can be done two-fold. First one can train a matrix such that 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 , the patches are actually correlated since they originate from the same initial source. Deep learning and neural networks are natural candidates unveiling such hidden correlation. Code and Data AvailabilityThe 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. AcknowledgmentsThis work has been supported by the Austrian Science Fund (FWF), projects P 30747-N32 and P 33019-N. ReferencesM. 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
K. Wang, M. A. Anastasio,
“Photoacoustic and thermoacoustic tomography: image formation principles,”
Handbook of Mathematical Methods in Imaging, Springer, New York
(2011). Google Scholar
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
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
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
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
G. Paltauf et al.,
“Photoacoustic tomography with integrating area and line detectors,”
Photoacoustic Imaging and Spectroscopy, 251
–264 CRC Press(
(2017). Google Scholar
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
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
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
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
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
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
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
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
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
S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing, Springer(
(2013). Google Scholar
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
|
Matrices
Design
Sensors
Image restoration
Compressed sensing
Biological research
Imaging systems