Real-Time 3D PET Image with Pseudoinverse Reconstruction

Real-time positron emission tomography (PET) may provide information from first-shot images, enable PET-guided biopsies, and allow awake animal studies. Fully-3D iterative reconstructions yield the best images in PET, but they are too slow for real-time imaging. Analytical methods such as Fourier back projection (FBP) are very fast, but yield images of poor quality with artifacts due to noise or data incompleteness. In this work, an image reconstruction based on the pseudoinverse of the system response matrix (SRM) is presented. w. To implement the pseudoinverse method, the reconstruction problem is separated into two stages. First, the axial part of the SRM is pseudo-inverted (PINV) to rebin the 3D data into 2D datasets. Then, the resulting 2D slices can be reconstructed with analytical methods or by applying the pseudoinverse algorithm again. The proposed two-step PINV reconstruction yielded good-quality images at a rate of several frames per second, compatible with real time applications. Furthermore, extremely fast direct PINV reconstruction of projections of the 3D image collapsed along specific directions can be implemented.


Motivation
Positron emission tomography (PET) is widely used in medical imaging. In many applications of PET such as first-shot tests for chemicals [1] or PET image-guided surgery or biopsy [2], it is crucial to obtain images as quickly as possible. The frame rate achievable will depend on the sensitivity of the acquisition system, the hardware used for data acquisition, and the software employed for data processing and image reconstruction [3]. State-of-the-art PET systems [4][5][6] can process up to 10 million coincidences per second, providing useful images in very short time frames. The challenge to be considered here is to provide useful quasi real-time images, that is, at a rate of a few frames per second with appropriate image quality. To derive the axial SRM, we considered an image space discretized in axial (z) and length coordinates (w) along the LOR (Figure 1b). The probability for each voxel to be connected to a given LOR was obtained as follows: ( , , , ) = · ( , , , ) , (2) where ( , , , ) is the minimum distance between the voxel at ( , ) and the LOR given by ( , ) (Figure 1b).
In a similar way, each element of the transaxial SRM can be modeled as ( Figure 1c where ( , , , ) is the distance from the voxel ( , ) to the LOR ( , ) [39]. In order to ease operations with these matrices, we combined some variables in the same bin index. We used ( , ) and ( , , ) instead of ( , , , ) and , ( , , , ) . In what follows, we will use this notation.

Pseudoinverse of the SRM
The solution of the least-squares problem [40] for a system with linear equations ( = • ) can be obtained using the pseudoinverse of the matrix [41][42][43]: || − · || ≥ || − · · || → (|| · || ) = · , To derive the axial SRM, we considered an image space discretized in axial (z) and length coordinates (w) along the LOR (Figure 1b). The probability for each voxel to be connected to a given LOR was obtained as follows: where d(z, w, z 1 , z 2 ) is the minimum distance between the voxel at (z, w) and the LOR given by (z 1 , z 2 ) (Figure 1b). In a similar way, each element of the transaxial SRM can be modeled as (Figure 1c): where d(z, y, ρ, θ) is the distance from the voxel (x, y) to the LOR (ρ, θ) [39]. In order to ease operations with these matrices, we combined some variables in the same bin index. We used A z (zw, z 1 z 2 ) and A xy (x, y, ρθ) instead of A z (z, w, z 1 , z 2 ) and A x,y (x, y, ρ, θ). In what follows, we will use this notation.

Pseudoinverse of the SRM
The solution of the least-squares problem [40] for a system with linear equations (Y = A·X) can be obtained using the pseudoinverse A † of the matrix A [41][42][43]: Appl. Sci. 2020, 10, 2829 4 of 18 The singular value decomposition (SVD) [44][45][46][47][48] was used here to compute the pseudoinverse of the matrix. In this method, matrix A is first decomposed into the product of two orthonormal matrices U and V and a diagonal matrix S containing the singular values (s ≥ 0) of A: With that decomposition, the pseudoinverse A † of the matrix results: where the reciprocal singular values of the matrix S are given by S † .

Regularization of the Pseudoinverse
As can be noted, when singular values are very small or null, their reciprocals are not well defined. Small singular values correspond to high frequency elements, which contain the finer details of the image (interesting to keep) as well as most of the noise of the data [15,16,49] (which amplifies the noise in the reconstructed image). To deal with this problem, we used regularization methods [50]. A common regularization approach is truncated singular value decomposition (TSVD), which is obtained by setting to 0 the reciprocal singular values larger than a chosen threshold [51].
A less extreme approach is Tikhonov regularization [50], which replaces the reciprocal singular values 1/s by: where k is a regularization parameter [52,53]. Very often in iterative reconstruction methods, the number of iterations serves as a regularization parameter [54]. A small number of iterations corresponds to large regularizations. In fact, it was found that a PINV regularization scheme such as the resulting images coincided with the ones of the iterative Landweber algorithm [55]. The Landweber regularization would replace the reciprocal values by (see Appendix A) [56]: where n is the number of Landweber iterations. This implies that through the PINV, an iterative Landweber reconstruction can be performed in one single step. Furthermore, the number of iterations in a Landweber regularization ( Figure 2) can be related with the regularization parameter k in Tikhonov regularization as: We present our results in terms of equivalent Landweber iterations, allowing us to study resolution recovery and noise in the image in terms of the number of equivalent iterations, as is usually done for iterative methods.

Axial Rebinning and Transaxial Reconstruction with the Pseudoinverse
Following the way the SRM was modeled into two independent components, the reconstruction process can be split into two smaller problems: axial rebinning (PINV z ) and 2D reconstruction of the corresponding rebinned slices (2D-PINV).
A † z (z, z 1 z 2 ) was obtained as the pseudoinverse of A z (zw, z 1 z 2 ), obtained from Equation (2) and then collapsing (summing over) in wt. On the other hand, where A † xy (xy, ρθ) is the pseudoinverse of A xy (xy, ρθ), obtained from Equation (3).

Axial Rebinning and Transaxial Reconstruction with the Pseudoinverse
Following the way the SRM was modeled into two independent components, the reconstruction process can be split into two smaller problems: axial rebinning (PINVz) and 2D reconstruction of the corresponding rebinned slices (2D-PINV).
Since both processes are independent and only involve matrix products, the reconstruction can be performed in a single step in which the 3D data are multiplied by each matrix from the left and right: The workflow of the reconstruction process is represented schematically in Figure 3. Since both processes are independent and only involve matrix products, the reconstruction can be performed in a single step in which the 3D data are multiplied by each matrix from the left and right: The workflow of the reconstruction process is represented schematically in Figure 3.

Axial Rebinning and Transaxial Reconstruction with the Pseudoinverse
Following the way the SRM was modeled into two independent components, the reconstruction process can be split into two smaller problems: axial rebinning (PINVz) and 2D reconstruction of the corresponding rebinned slices (2D-PINV).
Since both processes are independent and only involve matrix products, the reconstruction can be performed in a single step in which the 3D data are multiplied by each matrix from the left and right: The workflow of the reconstruction process is represented schematically in Figure 3.

List-Mode Data and Event by Event Reconstruction
As it has been shown in previous works [3], the pseudoinverse of the SRM can be used for list-mode data reconstruction. In a 3D problem, each column of the PINV of the SRM, A † (:, i), contains the region of the image that would produce a coincidence in the given LOR i. Taking this into account, the image can be refreshed after every event. A similar approach is also feasible in the proposed scheme. Each coincidence with coordinates (ρ, θ, z 1 , z 2 ) can be separated into transaxial and axial coordinates. Using A † z , for every coincidence with z 1 z 2 axial coordinates, each z slice (z) of the image is weighted by the corresponding element of the pseudoinverse of the axial SRM A † z (z 1 z 2 , z). Then, for each axial slice, coordinates (x, y) of the image are updated event by event using each row of the transaxial SRM pseudoinverse A † xy (xy, ρθ). Although this is fast and feasible for a 2D case as presented by Selivanov and Lecompte [3], it takes longer than reconstructing with the whole sinogram, except for when the number of coincidences is very small. In our case, reading the list-mode coincidences and histograming them in a sonogram is faster than event by event reconstruction, even for frames as small as 0.1 s for the preclinical SuperArgus PET-CT scanner introduced in Section 2.2.1.

Projection over Planes
In many cases, real-time images are used to follow radiotracer distribution along the body of a patient for which the sum of the image collapsed in specific directions instead of a 3D presentation of the full image may be enough. This way of presenting real-time PET images is also convenient as short frame images would typically contain a very low number of counts, and 2D projections will be less noisy in this case.
In general, to obtain images collapsed over planes, a 3D image is needed, which is later projected. However, using the pseudoinverse, this process can be reduced to one single step, if the pseudoinverse, instead of the image, is collapsed once and for all. This leads to much smaller PINV matrices, and thus much faster matrix products.
In this work, we used 2D-PINV matrices collapsed over the X and Y directions. These matrices, combined with PINV z , produces projections in the XZ and YZ planes, respectively. Note that if we sum over all the slices obtained from the rebinning step, and then reconstruct the resulting slice, we obtain the equivalent of the 3D image collapsed into the XY plane.

Simulated Data
PeneloPET [57] was used for the simulations presented in this work. It is a Monte-Carlo simulator, based on PENELOPE [58]. It considers the main physical effects during PET acquisitions such as positron range [19], non-collinearity, nuclear decay of the isotope, etc. 18 F point sources (radius 0.5 mm) were simulated for a geometry similar to the Biograph TPTV (Siemens) scanner [58] to evaluate the axial resolution of the different rebinning strategies discussed. The image noise was evaluated with another simulation of two homogeneous cylinders. The axial resolution was obtained from the full width at half maximum (FWHM) of the 1D-Gaussian fit to the axial profile obtained in the center of each source. The noise level was obtained as the ratio of the standard deviation to the average value in spherical Regions Of Interest (ROI's) with a radius of 5 cm located inside the simulated cylinders. We used the same sinogram dimensions as the actual clinical scanner data described below. The rebinned 2D slices were reconstructed using FBP with a Hamming filter (cutoff 0.5).
A thin slice Defrise phantom [59] was simulated to study axial rebinning with different methods. In this case, the simulated geometry was the one corresponding to the preclinical SuperArgus PET-CT scanner (SEDECAL) [6]. The simulated phantom was composed of several disks of 8 cm diameter and 3 mm thickness. The gap between disks was 3 mm. We used the same sinogram dimensions as the actual preclinical scanner data, as described in Section 2.2.1.

Preclinical PET/CT Scanner
We used the six ring version of the SuperArgus PET/CT scanner (Figure 4. (left)). This configuration consists of 24 arrays of 13 × 13 crystals with 1.55 mm pitch size. The number of sinogram bins was 175 (radial) × 128 (angular). Radial bins were 0.5 mm in size, configuring a 8.75 cm transaxial FOV. A total of 195 slices of 0.775 mm each were used in the axial direction. The maximum ring difference [28] was set to 97 and a span factor of 19 was applied, resulting in a total of 1185 sinograms. The axial pseudoinverse matrices A † z (z 1 z 2 , z) had a size of 1185 × 195. For the 2D matrices, 175 voxels were chosen for the image FOV. The size of 2D matrices for the SuperArgus PET/CT was around 2.6 Gb. Both σ z and σ xy of the Gaussian distributions of the TOR were 0.6 mm.
We used the six ring version of the SuperArgus PET/CT scanner (Figure 4. (left)). This configuration consists of 24 arrays of 13 × 13 crystals with 1.55 mm pitch size. The number of sinogram bins was 175 (radial) × 128 (angular). Radial bins were 0.5 mm in size, configuring a 8.75 cm transaxial FOV. A total of 195 slices of 0.775 mm each were used in the axial direction. The maximum ring difference [28] was set to 97 and a span factor of 19 was applied, resulting in a total of 1185 sinograms. The axial pseudoinverse matrices ( , ) had a size of 1185 × 195. For the 2D matrices, 175 voxels were chosen for the image FOV. The size of 2D matrices for the SuperArgus PET/CT was around 2.6 Gb. Both of the Gaussian distributions of the TOR were 0.6 mm. For the evaluation of the image quality we used the IQ NEMA and Derenzo phantoms [60]. The rod diameters (Derenzo phantom) were 1.2, 1.6, 2.4, 3.2, 4, and 4.8 mm.
Two animal acquisitions injected with fluorodeoxyglucose (FDG) were also used to evaluate the performance of the PINV (PINVz + 2D-PINV) methodology: a 200 g rat and a cardiac acquisition of a 15 g mouse. Over the cardiac mouse acquisition, a quantitative comparison of 2D-PINV and FBP in real preclinical acquisitions was made. A profile over the heart was fit to two Gaussians (one to each side of the heart). The average of the standard deviation of both Gaussians is shown and the relative improvement of 2D-PINV over FBP was calculated as | − |/ . Moreover, we used a first pass acquisition of a mouse injected with FDG to demonstrate the real-time capabilities of PINV reconstructions.
Acquisitions were performed at the Instituto de Investigación Sanitaria Gregorio Marañón   For the evaluation of the image quality we used the IQ NEMA and Derenzo phantoms [60]. The rod diameters (Derenzo phantom) were 1.2, 1.6, 2.4, 3.2, 4, and 4.8 mm.
Two animal acquisitions injected with fluorodeoxyglucose (FDG) were also used to evaluate the performance of the PINV (PINV z + 2D-PINV) methodology: a 200 g rat and a cardiac acquisition of a 15 g mouse. Over the cardiac mouse acquisition, a quantitative comparison of 2D-PINV and FBP in real preclinical acquisitions was made. A profile over the heart was fit to two Gaussians (one to each side of the heart). The average of the standard deviation of both Gaussians is shown and the relative improvement of 2D-PINV over FBP was calculated as |σ 2D−PINV − σ FBP |/σ FBP . Moreover, we used a first pass acquisition of a mouse injected with FDG to demonstrate the real-time capabilities of PINV reconstructions.
Acquisitions were performed at the Instituto de Investigación Sanitaria Gregorio Marañón (Madrid, Spain). All experimental procedures were done in compliance with the European Communities Council

Clinical PET/CT Scanner
For the clinical data, we used 18 F acquisitions obtained with a Biograph TPTV [58,61] (Figure 4. (right)) scanner at the Medical University of Vienna (Austria) with a 34 cm radial field of view, four rings of 48 arrays with 13 × 13 scintillator crystals, with an axial length of 4 mm of each crystal. The total number of radial and angular sinogram bins was 336 × 336 [61]. Radial bins were 2.024 mm in size. A total of 109 slices of 2 mm each were used in the axial direction. The maximum ring difference was set to 38 and a span factor of 11 was applied, yielded into 559 sinograms. The axial matrices A † z (z 1 z 2 , z) had a size of 559 × 109 (238 KB). For the 2D matrices, 225 voxels were chosen, providing a FOV of 45.54 cm. The size of the 2D matrices was around 21 GB. The σ xy and σ z for the TORs in the matrices were 2 mm.

Image Quality Evaluation
We evaluated the reconstructed images of the phantoms following the National Electrical Manufacturers Association (NEMA) NU 2-2007 and NU 4-2008 protocols for clinical and preclinical phantoms, respectively [62,63].

•
Noise level in the image (%): evaluated as the standard deviation of the activity in a uniform region divided by the average activity in that region. • Resolution (%), following the NEMA standards [62,63].
For the preclinical scanner, the mouse-size IQ phantom was used. Sinograms were corrected by attenuation and random coincidences. In the case of the clinical NEMA, scatter corrections were also included.
We present a comparison of SSRB + FBP or FORE + FBP with the proposed method (PINV z + 2D-PINV). Unless otherwise noted, execution times correspond to one thread of a CPU E5-2640 v4 @ 2.40 GHz in a Linux computer with Centos 7 OS.  In Figure 6, we can see a coronal view of the reconstruction of a Defrise disk phantom where an improvement in the pseudoinverse rebinning can be appreciated when compared to SSRB. We can see that FORE and PINV z also provided similar detail. Figure 5. Noise vs. resolution recovery for different rebinning methods (SSRB, FORE, and PINVz). Point sources were simulated in a scanner with a geometry similar to Biograph TPTV with MRD 38 and SPAN 11. Numbers in brackets at the FORE and SSRB points indicate the distance (R,Z) to the center of the scanner in cm, for example, (10,0) represents a source at a radial distance of 10 cm from the center of the scanner at the axis. Several points for PINVz show the axial resolution vs. recovery for a different Landweber iterations, from right to left 5, 8, 10, 15, and 20 iterations. Outside the axis of the scanner, SSRB was not accurate, while FORE provided the best results overall, although with some sensitivity to distance to the scanner axis. PINVz keeps a more uniform resolution across the whole FOV.

Preclinical Data
Results of a PINV reconstruction (PINV z + 2D-PINV) were compared with the FORE + FBP images. We also studied the isolated 2D reconstruction with FBP and PINV from the same FORE rebinned sinogram for the IQ-NEMA phantom. These results are presented in Figure 7 and Table 1. In Table 1, we compared the recovery coefficients for the IQ-NEMA phantom at approximately the same level of uniformity for every rebinning + 2D reconstruction combination. The 2D-PINV slice reconstruction systematically outperformed (around 5-10% better for the smallest capillary and up to 15% in the 2 mm capillary). As expected, FORE rebinning provided the best recovery, improving marginally on the pseudoinverse rebinning results.
Furthermore, in Figure 8, we present the results for a Derenzo filled with 250 µCi of 18 FDG measured with a well-counter with an accuracy of ±5% [11]. It can be seen as an improvement in resolution for 2D-PINV reconstruction when compared to FBP.  In Table 1, we compared the recovery coefficients for the IQ-NEMA phantom at approximately the same level of uniformity for every rebinning + 2D reconstruction combination. The 2D-PINV slice reconstruction systematically outperformed (around 5-10% better for the smallest capillary and up to 15% in the 2 mm capillary). As expected, FORE rebinning provided the best recovery, improving marginally on the pseudoinverse rebinning results.
Furthermore, in Figure 8, we present the results for a Derenzo filled with 250 μCi of 18 FDG measured with a well-counter with an accuracy of ±5% [11]. It can be seen as an improvement in resolution for 2D-PINV reconstruction when compared to FBP.  Figure 9 shows the cardiac region of a 15 g mouse injected with 170 μCi of 18 FDG, with slices reconstructed with FBP and 2D-PINV from a PINV rebinned sinogram. Improvement in spatial resolution (a relative 15%) for 2D-PINV reconstructions can be clearly seen. Figure 10 shows the PINVz + 2D-PINV versus SSRB + FBP reconstructions of a 200 g rat injected with 300 μCi of 18 FDG where the improvement using PINV was evident. Images from both 30 min (bottom) and 5 s (top) acquisitions are shown, the second demonstrating the capabilities for near real time imaging.  Figure 9 shows the cardiac region of a 15 g mouse injected with 170 µCi of 18 FDG, with slices reconstructed with FBP and 2D-PINV from a PINV rebinned sinogram. Improvement in spatial resolution (a relative 15%) for 2D-PINV reconstructions can be clearly seen.
reconstructed with FBP and 2D-PINV from a PINV rebinned sinogram. Improvement in spatial resolution (a relative 15%) for 2D-PINV reconstructions can be clearly seen. Figure 10 shows the PINVz + 2D-PINV versus SSRB + FBP reconstructions of a 200 g rat injected with 300 μCi of 18 FDG where the improvement using PINV was evident. Images from both 30 min (bottom) and 5 s (top) acquisitions are shown, the second demonstrating the capabilities for near real time imaging.

Clinical Scanner
We studied the recovery coefficients for the clinical IQ NEMA phantom from images rebinned with SSRB, FORE, and PINVz and reconstructed with the same parameters as in Section 3.1 (i.e., FBP with a Hamming filter). Additionally, we studied a reconstruction PINVz + 2D-PINV. Table 2 shows the recovery coefficients and uniformity over the six spheres following NEMA 2-2007. These results are summarized in Figure 11 and Table 2.

Clinical Scanner
We studied the recovery coefficients for the clinical IQ NEMA phantom from images rebinned with SSRB, FORE, and PINV z and reconstructed with the same parameters as in Section 3.1 (i.e., FBP with a Hamming filter). Additionally, we studied a reconstruction PINV z + 2D-PINV. Table 2 shows the recovery coefficients and uniformity over the six spheres following NEMA 2-2007. These results are summarized in Figure 11 and Table 2.  Figure 11. Transverse FBP and 2D-PINV images from data rebinned with PINV, FORE, and SSRB for the clinical IQ NEMA phantom.

Reconstruction Time
Reconstruction times are shown in Tables 3 and 4. Time to generate the pseudoinverse matrix from the SRM and pre-load it to RAM memory was not included as part of the reconstruction, as both processes are performed before it.

Real-Time PET Image Reconstruction
As described in the previous sections, pseudoinverse methods allow for extremely fast, direct reconstruction of projected or collapsed images into any plane. Figure 12 shows the results of a reconstruction collapsed onto the YZ plane during a first shot acquisition for a 13 g mouse injected with 120 μCi of 18 FDG, on the SuperArgus scanner. Six frames of one second each are presented, clearly showing the activity coming in along the tail vein and further into the lungs, vascularity, and remaining circulatory system. The reconstruction of the accumulated 6 s is also shown for reference.

Reconstruction Time
Reconstruction times are shown in Tables 3 and 4. Time to generate the pseudoinverse matrix from the SRM and pre-load it to RAM memory was not included as part of the reconstruction, as both processes are performed before it.

Real-Time PET Image Reconstruction
As described in the previous sections, pseudoinverse methods allow for extremely fast, direct reconstruction of projected or collapsed images into any plane. Figure 12 shows the results of a reconstruction collapsed onto the YZ plane during a first shot acquisition for a 13 g mouse injected with 120 µCi of 18 FDG, on the SuperArgus scanner. Six frames of one second each are presented, clearly showing the activity coming in along the tail vein and further into the lungs, vascularity, and remaining circulatory system. The reconstruction of the accumulated 6 s is also shown for reference.

Real-Time PET Image Reconstruction
As described in the previous sections, pseudoinverse methods allow for extremely fast, direct reconstruction of projected or collapsed images into any plane. Figure 12 shows the results of a reconstruction collapsed onto the YZ plane during a first shot acquisition for a 13 g mouse injected with 120 μCi of 18 FDG, on the SuperArgus scanner. Six frames of one second each are presented, clearly showing the activity coming in along the tail vein and further into the lungs, vascularity, and remaining circulatory system. The reconstruction of the accumulated 6 s is also shown for reference.   Table 5 presents the total time needed to obtain these reconstructions using a multithreaded version of the code.

Discussion
In this work, we showed that pseudoinverse reconstruction of 3D PET data using PINV z + 2D-PINV provided better (Tables 1 and 2) and faster reconstructions than SSRB + FBP (Tables 3 and 4). It is worth noting that in FBP, it is further needed to estimate the missing data (filling in the gaps in between detectors) [22], and this time was left out from Table 4. On the other hand, pseudoinverse reconstruction only requires knowing the position of the missing data in order to skip those, further speeding the process. Pseudoinverse axial rebinning was much faster than FORE, with just a moderately worse performance in terms of noise/resolution. As has been noted in previous works [3], there is a 2D-PINV matrix that would be fully equivalent to FBP. However, in the construction of the 2D-PINV, it is also possible to add geometrical information and physics effects, which introduce blurring to the TORs. This is the main reason for the improved resolution of in-slice PINV images compared to FBP ones. The matrices for PINV rebinning and slice reconstruction can be precomputed and stored once and for all for a given scanner, making it easy for the user to produce their own reconstructions.
Although it would be desirable to apply PINV methods in 3D altogether instead of the division into axial and transaxial components employed in this work, building a 3D PINV matrix currently seems impractical as 3D PINV arrays for modern complex scanners will not fit in the memory of common computers. Furthermore, as we have presented here, the noise exhibited by PINV rebinning methods is slightly larger than that of FORE. Better regularization techniques [50] or improvements in the SRM could remedy this.
New regularization techniques could involve filtering of the reconstructed image or even the data before reconstruction. Any linear filter (F) can be applied during PINV reconstruction from its matrix representation: X A + ·Y → X = F·X = F·A + ·Y = F·A + ·Y, Therefore, applying the filter in the pseudoinverse before the reconstruction would result in the same image as using the filter over the reconstructed image. The same applies to any linear filter applied in the data domain: Matrices to represent the SRM in this work were obtained analytically for convenience, but they can also be obtained from Monte Carlo simulations, with likely improvement in image quality without affecting the speed of the method.
Reconstruction times presented in Section 3.4 corresponded to only one CPU core. Multi-CPU or Graphics Processing Unit (GPU) strategies can be used to further accelerate the reconstructions. For instance, we have verified a very straightforward method to use multicore 2D by assigning each slice to a CPU-thread, reducing the total time almost linearly within the number of threads employed. In Tables 3 and 4, reconstruction times for smaller sinograms were also presented. Although some resolution was lost with these smaller sinograms, their very short reconstruction times make them suitable for real time imaging.
Moreover, for real time applications, most often, only projection images into predefined planes will be needed. This projection operation can be incorporated into the pseudoinverse array before reconstruction, drastically reducing the size of the matrices needed as well as the reconstruction time as shown in Figure 12 and Table 5.

Conclusions
We showed that PINV z outperformed SSRB in terms of resolution for off-axis sources and provided a more homogeneous resolution across the FOV than SSRB or FORE. Furthermore, pseudoinverse methods are fast and can be implemented in real time applications, while other methods such as FORE are not currently fast enough.
Additionally, 2D-PINV reconstruction can take into account the physical and geometrical information in the SRM, which may result in improved resolution compared to analytical reconstructions, without any impact in the reconstruction time. Funding: Part of the calculations in this work were performed in the "Clúster de Cálculo para Técnicas Físicas" funded in part by UCM and in part by UE Regional Funds. We acknowledge the support from the Spanish Government (FPA2015-65035-P, RTC-2015-3772-2, and RTI2018-095800-A-I00), Comunidad de Madrid (S2013/MIT-3024 TOPUS-CM, B2017/BMD-3888 PRONTO-CM), and European Regional Funds. This work was also supported by the EU's H2020 under MediNet, a Networking Activity of ENSAR-2 (grant agreement 654002), and by a NIH R01 CA215700-2 grant. The CNIC is supported by the Ministerio de Ciencia, Innovación y Universidades and the Pro CNIC Foundation, and is a Severo Ochoa Center of Excellence (SEV-2015-0505).