Image Deconvolution with Hybrid Reweighted Adaptive Total Variation (HRATV) for Optoacoustic Tomography

: Optoacoustic tomography (OAT) is a hybrid biomedical imaging modality that usually employs a transducer array to detect laser-generated ultrasonic signals. The reconstructed image suffers low contrast and degraded resolution due to the limited bandwidth and the spatial directivity of the transducer element. Here, we introduce a modiﬁed image deconvolution method with a hybrid reweighted adaptive total variation tailored to improve the image quality of OAT. The effectiveness and the parameter dependency of the proposed method are veriﬁed on standard test images. The performance of the proposed method in OAT is then characterized on both simulated phantoms and in vivo mice experiments, which demonstrates that the modiﬁed deconvolution algorithm is able to restore the sharp edges and ﬁne details in OAT simultaneously. The signal-to-noise ratios (SNRs) of the target structures in mouse liver and brain were improved by 4.90 and 12.69 dB, respectively. We also investigated the feasibility of using Fourier ring correlation (FRC) as an indicator of the image quality to monitor the deconvolution progress in OAT. Based on the experimental results, a practical guide for image deconvolution in OAT was summarized. We anticipate that the proposed method will be a promising post-processing tool to enhance the visualization of micro-structures in OAT.


Introduction
Optoacoustic tomography (OAT), also known as photoacoustic tomography (PAT), is a fast-evolving biomedical imaging technique that has been widely investigated in fields including vascular morphology [1,2], cancer research [3], cardiology [4], and neuroscience [5,6]. It combines the high optical contrast and ultrasonic resolution by employing the photoacoustic effect, in which ultrasonic waves are generated by irradiating tissue with nano-seconds laser pulse [7,8]. In early OAT systems, the laser induced ultrasonic signals are detected by a scanning single-element transducer around the tissue and then the initial pressure field is reconstructed using analytical or model based algorithms [9]. Later, large detector arrays with an arc or circular shape are employed to achieve fast two-dimensional imaging [10,11]. Recently, a multi spectral optoacoustic tomography (MSOT) system with a 512-element transducer array has been developed for full-view high resolution imaging of small animals [12]. It has been concluded that increasing the number of detectors achieves enhanced detection sensitivity and spatial resolutions, as more detection angles are considered in image reconstruction and smaller elements are utilized without sacrifice of the field of view (FOV). However, parallel acquisition of hundreds of signal channels requires multiple analog-to-digital converters operating at sampling frequency of tens of MHz, and consequently restricts the use of large transducer arrays [13]. Therefore, a trade-off between the image quality and the hardware complexity must be considered in OAT systems. Furthermore, the axial resolution is also constrained by the limited bandwidth of the transducer and the image quality is reduced due to the loss of high frequency components.
As the number and the bandwidth of the transducer element cannot go up forever, an alternative way to improve the quality of OAT is via image deconvolution, which is considered as an energy compact process that accumulates the energy to its original position. In OAT, the lateral blurring can be modeled as the convolution of the spatial directivity of the transducer and the initial pressure field, while the axial blurring can be considered as the convolution of the transducer impulse response and the initial pulsed ultrasonic signal [14]. In total, the spatial blurring can be expressed by the point spread function (PSF) of the OAT system and the blurred image is simplified as the convolution of the PSF and the initial pressure field. Reversely, the observed image can be deblurred via a deconvolution process.
Wang first demonstrated a deconvolution based OAT reconstruction method that directly decomposes the detected signals from the point source signal [15]. This method is computationally efficient, but the deconvolution results are sensitive to the noise and the measurement of the point source signal. Wiener filter deconvolution has also been studied in optoacoustic imaging to restore the initial wideband signals [16][17][18]. More recently, Guo proposed a Wiener filter with empirical mode decomposition (EMD) in linear array based OAT to improve the axial resolution and reduce unexpected artifacts [19]. A deconvolution method with Tikhonov regularization has also been presented to correct the optoacoustic signal based on the impulse response of the transducer [20,21]. The aforementioned works only focus on one-dimensional signals and require either the point source measurement or the impulse response of the transducer as a priori knowledge. Another drawback is that these methods are not robust to high level noise. Therefore, only the axial blurring can be reduced and limited improvement is achieved.
Richardson-Lucy (RL) deconvolution, an iterative image restoration process, has gained increasing popularity in the field of medical imaging [22][23][24], thanks to its implementation of the maximum likelihood estimation and the robustness to the high level noise. When the exact PSF of the system is unknown, the deconvolution is termed as RL blind deconvolution (RLBD) [25]. Because of its success in image restoration, RL deconvolution has also been introduced to the field of optoacoustic imaging to improve the image resolution and SNR. Jetzfellner and Ntziachristos first investigated the feasibility of using RLBD as a post-processing method to increase the resolution of OAT [26]. Later, the group proposed a multi-bandwidth RL deconvolution approach to improve the resolution along the elevation dimension [27]. RL deconvolution has also been employed in optoacoustic microscopy (OAM) to improve both the axial and the lateral resolution [28][29][30]. However, few works have been done on RL deconvolution of OAT images.
The main drawback of RL deconvolution is that it tends to amplify noise after iterations and thus requires regularization. Dey first proposed a regularized RLBD algorithm with total variation (RLBD-TV) filter in confocal microscopy to restore the high resolution image [31], as total variation (TV) is a commonly used regularization filter in image denoising that achieves a good trade-off between edge preserving and noise suppression [32]. A modified RLBD with spatially adaptive TV (RLBD-SATV) that incorporates the spatial information to the regularization term has been proposed in astronomical imaging to achieve automatic parameter selection based on the Hessian matrix of the image [33]. Inspired by their works, a high resolution OAM based on the RLBD-TV method was developed for imaging of retinal microvasculature [34]. However, TV regularization favors a piece-wise constant solution that leads to the so-called "staircase effect" on the restored image. The staircase effect will smooth out the difference between adjacent structures and reduce the resolution. To overcome this problem, a fourth-order partial differential equation (FPDE) filter was proposed in denoising of MRI images [35]. Higher-order PDEs are good at restoring smooth regions and fine textures, while a TV filter performs better at edges and noise suppression. As OAT images typically contain rich blood vessels and micro-structures, the staircase effect caused by the TV filter is not favorable in OAT deconvolution because it will blur the fine textures and degrade the resolution. The FPDE filter can recover the fine details and smooth regions, but does not produce a good solution to edge preservation and noise suppression. Therefore, it is intrinsic to combine both methods in RL deconvolution of OAT. Motivated by the work in [36,37], we first introduced the spatially adaptive fourthorder PDE (SAFPDE) in RLBD to overcome the staircase effect in RLBD-SATV. Then we proposed a modified RLBD with hybrid reweighted adaptive TV (RLBD-HRATV) that automatically combines the SATV and the SAFPDE in the deconvolution process to recover both the edges and fine details. To validate the effectiveness of the proposed method, we tested the deconvolution performance on commonly used standard test images. Then we investigated the effects of the spatial directivity and frequency bandwidth of the transducer on OAT using the k-wave toolbox. Finally, we analyzed the performance of different deconvolution methods on simulated OAT phantoms and real in vivo OAT images. In iterative deconvolution, how to choose appropriate stopping criteria has been a long-standing problem as there is no ground truth in a real scenario. To address this issue, we investigated whether the Fourier ring correlation (FRC) [38] can be used as an image quantifier in OAT to monitor the deconvolution process and terminate the iteration properly. Based on the simulation and experimental results, we discussed the initialization method and parameter selection in OAT deconvolution and provided a practical guide for deconvolution of OAT images in a real scenario.
The main contributions of this study are summarized as follows:

1.
A modified RLBD algorithm with hybrid reweighted adaptive TV (HRATV) is proposed to restore the edges as well as the fine details in OAT.

2.
The feasibility of using Fourier ring correlation (FRC) to monitor the deconvolution of OAT was evaluated.

3.
A practical guide for deconvolution of OAT was summarized based on the simulation and experimental results.
The rest of the paper is organized as follows. Section 2 introduces the typical image restoration framework and the RLBD solution. Subsequently, the proposed RLBD-HRATV method is presented. Section 3 presents the simulation methods and the in vivo OAT experiments on nude mice. Image quality metrics are introduced in this section as well. Section 4 provides the comparative analysis of the three RL deconvolution methods, namely RLBD-SATV, RLBD-SAFPDE and RLBD-HRATV. The initialization method and parameter selection were discussed and a practical guide for OAT deconvolution was provided. Conclusions were drawn in Section 5.

Hybrid Reweighted Adaptive Total Variation Regularization
In this section, we will present the formation of the modified RLBD with hybrid reweighted adaptive total variation (HRATV). Before this, we first introduce the typical image restoration framework and the related RLBD solutions.
Image deblurring is one of the fundamental problems in image processing, which aims to restore the blurred edges, smooth regions and fine textures from the observed image. The observed image is often degraded by a blur kernel and corrupted by noise, which can be modeled as a convolution process under the assumption that the system is linear and spatially invariant: where i and u denote the observed image and the object, h is the PSF of the system, and n is the noise. Here we omit the spatial coordinates for simplicity. Deconvolution is a process to estimate the object u from the degraded image i. As there exists not only one solution that is mathematically correct, the deconvolution is known as an ill-posed inverse problem especially in the presence of noise and requires a priori knowledge to obtain a suitable solution. The Richardson-Lucy (RL) algorithm employs a probabilistic framework that finds the estimation of the object by maximizing the likelihood probability p(i|u) [39,40]. The RL algorithm incorporates a priori that the noise follows a Poisson distribution. Therefore, the likelihood probability can be expressed as: where x is the total set of the spatial coordinates. To maximize Equation (2), one alternative way is to minimize the equivalence of −log p(i|u) : By minimizing Equation (3) using a multiplicative gradient-based iterative algorithm [31], one RL iteration can be derived as: where u k is the estimation of the object at k iteration. From Equation (4), it can be seen that the PSF of the system h must be provided to initialize the iteration. In the blind form RL deconvolution (RLBD) [25], the blur kernel h is updated at each iteration together with Equation (4): where h k denotes the estimation of the PSF at k iteration. The PSF h is constrained to non-negativeness and normalized to 1 as:

Related Works
As RLBD tends to amplify noise after several iterations, it is necessary to regularize the noise during the deconvolution process. TV regularization has been proven to achieve a good trade-off between edge preservation and noise reduction in the field of image denoising, which was defined in [32] as: where |∇u(x)| = (u x ) 2 + u y 2 , x and y represent the spatial coordinates. Dey [31] first proposed the RLBD-TV method by introducing the TV regularization term to the minimization function Equation (3): where λ 1 is the parameter to tune the strength of the regularization. By minimizing Equation (9) in the same way as Equation (3), one iteration of the object u k for RLBD-TV can be obtained as: where div represents the divergence. The performance of the TV regularization is sensitive to parameter λ 1 that has to be manually tuned to obtain a suitable solution, as a large λ 1 would be too strong to preserve the edges while a small λ 1 may not be enough to suppress the noise. Generally, λ 1 should be tuned in the range of [10 −5 , 1]. In Equation (9), it should be noted that a uniform parameter λ 1 is applied on the whole image, which is not ideal since the noise region usually needs stronger regularization than the edges. To solve this problem, Yan proposed the RLBD-SATV [33] method, which introduces a spatially adaptive λ 1 based on the difference eigenvalue edge indicator that automatically applies a unique λ 1 to each pixel of the image. The new minimization function is then: where β is a constant and D(x) represents the difference eigenvalue edge indicator. By minimizing Equation (11) using the multiplicative gradient-based algorithm, one iteration of the object u k for RLBD-SATV can be obtained: It was shown that RLBD-SATV achieves better image quality than RLBD-TV and behaves more robust to the selection of the regularization parameter. However, the staircase effect caused by TV regularization has not been addressed.

Spatially Adaptive Fourth-Order Partial Differential Equation (SAFPDE)
In OAT, it is preferable to remove the staircase effect, as it produces piece-wise constant solutions on blood vessels and micro-structures. To overcome the intrinsic drawback of RLBD-SATV and restore the fine details, we introduce the fourth-order partial differential equation (FPDE) filter to constrain the RLBD process. Unlike Equation (8), which computes the total variation norm of u, the FPDE filter is actually the total variation norm of ∇u, which was defined as [35]: where λ 2 is the regularization parameter. Inspired by the idea of the spatially adaptive parameter in SATV, we introduce the same edge indicator to FPDE, then the minimization function can be expressed as: By minimizing Equation (14) using the multiplicative gradient-based algorithm, one iteration of the object u k for RLBD-SAFPDE can be derived as: where F(u(x)) =

Hybrid Reweighted Adaptive Total Variation (HRATV)
Since SATV and SAFPDE regularizations have complementary merits in image restoration, it is appealing to combine the two filters together. Motivated by the combination strategy proposed in [36], we develop an RLBD algorithm with hybrid reweighted adap-tive total variation (RLBD-HRATV) that uses a weighting function (defined in Equations (18) and (19)) to update the weight of SATV and SAFPDE in a pixel-wise manner at each iteration. To avoid ambiguity, we denote the solution u(x) in Equations (12) and (15) as f (x) and g(x), respectively. Here, we rewrite Equations (12) and (15) as: where F(g(x)) = g xx Equations (16) and (17) is defined as: where u k is the combined solution at k iteration and θ k is the weighting matrix defined as: where c is a constant in the interval [0, 1]. The motivation of Equation (19) is to assign a large weight θ k to SATV on the edges and background regions, while assigning a large weight (1 − θ k ) to SAFPDE in the smooth regions.

Implementation Details
To discretize Equations (16) and (17), we adopt the finite difference scheme that is used to solve partial differential equations by approximating the space and time derivatives with finite differences [35,37]. The implementation details are provided in Supplementary Note S1. The implementation of the proposed RLBD-HRATV method is summarized as follows on Algorithm 1.
3: update object estimation f k+1 and g k+1 : 5: do combination: Output: restored image u k+1 By fixing the weighting matrix θ to 1, the RLBD-HRATV method degrades to RLBD-SATV and by fixing the weighting matrix θ to 0, the RLBD-HRATV method degrades to RLBD-SAFPDE. For standard test images, we use the normalized mean square error (NMSE) as the stopping criteria. While for OAT images, we use the Fourier ring correlation (FRC) as the stopping criteria since there is no ground truth in the real scenario. We will introduce the NMSE and the FRC in the following section.

Materials and Methods
The three RL deconvolution methods were tested on three datasets including the standard test images, simulated OAT phantoms, and in vivo OAT images. The algorithm was implemented using MATLAB R2018b on a 64-bit Windows 7 desktop with an Intel CPU at 3.10 GHz and 16 GB RAM.

Standard Test Images
whereû and u denote the estimation and the reference image, respectively. The SSIM was calculated using the MATLAB function ssim(). During iteration, the NMSE was used to monitor the image quality and tune the regularization parameters. All results are shown at the iteration when NMSE reaches the minimum. The initial guess of the PSF h 0 was set to 1 with the same size of the blur kernel, which was implemented as ones (size(blur kernel)) in MATLAB. The edge indicator constant β was set to 2 for all images.

OAT Simulations
Simulations of OAT were implemented using the k-Wave MATLAB toolbox, which is designed for the simulation and reconstruction of optoacoustic wave-fields using a k-space pseudo-spectral solution. Simulations were done in 2D to save computation time. The computational grid was set to 300 × 300 pixels with the total grid size of 10 × 10 mm 2 , as illustrated in Figure 1b. A perfectly matching layer (PML) with the thickness of 20 grid points was defined as the absorbing boundary. The transducer array used for signal detection contains 128 elements evenly distributed on a 270-degree-arc that centers at the center of the computational grid. The radius of the arc was 4.5 mm. The spatial directivity of the transducer element was simulated by setting its directivity angle pointing to the arc center. The directivity size of the transducer element was 10 pixels. The sound velocity of the medium was set to 1500m/s. The total computation time duration was automatically set to allow waves traveling from one corner of the computational grid to the opposite corner. The data collected by the transducer elements were contaminated by 1% noise to simulate a 40 dB SNR. A low-pass filter with cut-off frequency of 6.75 MHz was applied to the noisy sensor data to simulate the limited frequency bandwidth of the transducer. The OAT images were reconstructed using the time reversal method in k-Wave.
shown at the iteration when NMSE reaches the minimum. The initial guess of the PSF ℎ was set to 1 with the same size of the blur kernel, which was implemented as ones (size(blur kernel)) in MATLAB. The edge indicator constant was set to 2 for all images.

OAT Simulations
Simulations of OAT were implemented using the k-Wave MATLAB toolbox, which is designed for the simulation and reconstruction of optoacoustic wave-fields using a kspace pseudo-spectral solution. Simulations were done in 2D to save computation time The computational grid was set to 300 × 300 pixels with the total grid size of 10 × 10 mm 2 , as illustrated in Figure 1b. A perfectly matching layer (PML) with the thickness of 20 grid points was defined as the absorbing boundary. The transducer array used for signal detection contains 128 elements evenly distributed on a 270-degree-arc that centers at the center of the computational grid. The radius of the arc was 4.5 mm. The spatial directivity of the transducer element was simulated by setting its directivity angle pointing to the arc center. The directivity size of the transducer element was 10 pixels. The sound velocity of the medium was set to 1500m/s. The total computation time duration was automatically set to allow waves traveling from one corner of the computational grid to the opposite corner. The data collected by the transducer elements were contaminated by 1% noise to simulate a 40 dB SNR. A low-pass filter with cut-off frequency of 6.75 MHz was applied to the noisy sensor data to simulate the limited frequency bandwidth of the transducer. The OAT images were reconstructed using the time reversal method in k-Wave. To study the effects of the frequency bandwidth and spatial directivity on OAT, we simulated the PSF functions of the imaging system under different setups. The calculation of the PSF was done by placing a point source in the imaging field and then reconstructing the image of this point. Here, four points were placed around the grid center with the distance to the center as 80, 60, 40, and 20 pixels. The initial pressure magnitude of the source was set to 1.
To characterize the deconvolution performance in OAT, a blood vessel phantom was simulated in k-Wave. The reconstructed image was then processed using RLBD-SATV RLBD-SAFPDE and RLBD-HRATV, respectively. As there is no reference in real applications, the aforementioned NMSE is not applicable to observe the deconvolution process and the SSIM and PSNR cannot be used to quantify the image quality as well. Fourier ring correlation (FRC) is a quantifier that can be calculated from a single image without any reference and has been demonstrated to be a powerful metric to monitor the progress of iterative deconvolution in fluorescence microscopy [38]. Here, we investigated whether FRC could be used in OAT to observe the deconvolution process. The calculation of the FRC is shown in Figure 2. An OAT image is first split into two sub-images using the check- To study the effects of the frequency bandwidth and spatial directivity on OAT, we simulated the PSF functions of the imaging system under different setups. The calculation of the PSF was done by placing a point source in the imaging field and then reconstructing the image of this point. Here, four points were placed around the grid center with the distance to the center as 80, 60, 40, and 20 pixels. The initial pressure magnitude of the source was set to 1.
To characterize the deconvolution performance in OAT, a blood vessel phantom was simulated in k-Wave. The reconstructed image was then processed using RLBD-SATV, RLBD-SAFPDE and RLBD-HRATV, respectively. As there is no reference in real applications, the aforementioned NMSE is not applicable to observe the deconvolution process and the SSIM and PSNR cannot be used to quantify the image quality as well. Fourier ring correlation (FRC) is a quantifier that can be calculated from a single image without any reference and has been demonstrated to be a powerful metric to monitor the progress of iterative deconvolution in fluorescence microscopy [38]. Here, we investigated whether FRC could be used in OAT to observe the deconvolution process. The calculation of the FRC is shown in Figure 2. An OAT image is first split into two sub-images using the checkerboard pattern: one sub-image consists of the even columns and even rows, another consists of the odd columns and odd rows. Afterward, the 2D Fourier spectrum of each sub-image is calculated. Then the correlation of two spectrums is computed for each frequency bin of the spectrum using Equation (22): where F 1 and F 2 are the 2D Fourier spectrums of two sub-images, r i is the ith frequency bin, and r denotes the coordinates of the points on r i . The FRC curve reveals the change of the signal correlation along the spatial frequency axis. The high frequency components correspond to the details of the image. When the correlation drops to a pre-defined threshold, it assumes that there are no resolvable details and the noise dominates. The spatial frequency r cutoff at which the FRC reaches the threshold is considered as the resolution (1/r cutoff ), indicated by the red marker on the FRC curve in Figure 2. The 1/7 threshold is used in this work as it has been demonstrated to give a relatively reasonable estimation of the resolution [41]. The FRC was used to estimate the resolution of the image at each iteration of the deconvolution process. When the estimated resolution stops increasing, the iteration is terminated. erboard pattern: one sub-image consists of the even columns and even rows, another consists of the odd columns and odd rows. Afterward, the 2D Fourier spectrum of each subimage is calculated. Then the correlation of two spectrums is computed for each frequency bin of the spectrum using Equation (22): where and are the 2D Fourier spectrums of two sub-images, is the ith frequency bin, and denotes the coordinates of the points on . The FRC curve reveals the change of the signal correlation along the spatial frequency axis. The high frequency components correspond to the details of the image. When the correlation drops to a pre-defined threshold, it assumes that there are no resolvable details and the noise dominates. The spatial frequency at which the FRC reaches the threshold is considered as the resolution (1/ ), indicated by the red marker on the FRC curve in Figure 2. The 1/7 threshold is used in this work as it has been demonstrated to give a relatively reasonable estimation of the resolution [41]. The FRC was used to estimate the resolution of the image at each iteration of the deconvolution process. When the estimated resolution stops increasing, the iteration is terminated.

In vivo OAT experiments
To test the performance of the RL deconvolution methods in real OAT applications, nude mice were imaged using a commercial OAT system (MSOT inVision 256-TF, iTher-aMedical, Munich, Germany). All mice experiments were approved by the Ethics Committee of Suzhou Institute of Biomedical Engineering and Technology, Chinese Academy of Sciences. Specific pathogen (SPF) grade BALB/c nude mice (6-8 weeks old) were provided by the laboratory of the Animal Center of Soochow University. To acquire the OAT images, the mice were anesthetized with 2% isoflurane and kept on a custom-design station for animal preparation. A layer of medical ultrasound gel was applied on the skin of the mice for acoustic wave coupling. Afterward, the mice were covered by a thin transparent polyethylene film for waterproofing and transferred to the imaging chamber of the MSOT system. The system utilizes a 5 MHz transducer array with 256 elements arranged on a 270-degree arc for optoacoustic signal detection. The laser beam generated by a tunable optical parametric oscillator with pulse width of ~10 ns at 10 Hz was guided through fiber bundles to achieve 360° light illumination. The excitation wavelength employed for imaging was set to 850 nm. Each image was obtained by averaging 4 successive measurements to improve the image quality. The size of the image is 500 × 500 pixels with the pixel size of 40 µm, corresponding to a field of view (FOV) of 20 × 20 mm. The OAT images exported from the MSOT system were normalized to [0,1] and then processed using RLBD-SATV, RLBD-SAFPDE and RLBD-HRATV, respectively. FRC was used to monitor the deconvolution progress. Two regions of interest (ROI) were selected as the signal and the background noise to calculate the SNR.

In Vivo OAT Experiments
To test the performance of the RL deconvolution methods in real OAT applications, nude mice were imaged using a commercial OAT system (MSOT inVision 256-TF, iTher-aMedical, Munich, Germany). All mice experiments were approved by the Ethics Committee of Suzhou Institute of Biomedical Engineering and Technology, Chinese Academy of Sciences. Specific pathogen (SPF) grade BALB/c nude mice (6-8 weeks old) were provided by the laboratory of the Animal Center of Soochow University. To acquire the OAT images, the mice were anesthetized with 2% isoflurane and kept on a custom-design station for animal preparation. A layer of medical ultrasound gel was applied on the skin of the mice for acoustic wave coupling. Afterward, the mice were covered by a thin transparent polyethylene film for waterproofing and transferred to the imaging chamber of the MSOT system. The system utilizes a 5 MHz transducer array with 256 elements arranged on a 270-degree arc for optoacoustic signal detection. The laser beam generated by a tunable optical parametric oscillator with pulse width of~10 ns at 10 Hz was guided through fiber bundles to achieve 360 • light illumination. The excitation wavelength employed for imaging was set to 850 nm. Each image was obtained by averaging 4 successive measurements to improve the image quality. The size of the image is 500 × 500 pixels with the pixel size of 40 µm, corresponding to a field of view (FOV) of 20 × 20 mm. The OAT images exported from the MSOT system were normalized to [0, 1] and then processed using RLBD-SATV, RLBD-SAFPDE and RLBD-HRATV, respectively. FRC was used to monitor the deconvolution progress. Two regions of interest (ROI) were selected as the signal and the background noise to calculate the SNR.

Results and Discussion
In this section, we first show the performance of the three RL deconvolution methods on standard test images and then present the deconvolution results on both simulated and experimental OAT data. A comparative analysis regarding the edge preservation, details restoration and noise suppression was done among the three deconvolution methods. Last, the initialization procedure and parameters selection were discussed and a practical guide for deconvolution of OAT was summarized.

Simulated Test Images
As presented in Section 3.1, simulated test images were obtained by applying Gaussian blur and Poisson noise to the Barbara, Lena and Circuit image, respectively. The NMSE was used to monitor the deconvolution progress and tune the regularization parameters. The degraded image i was used as the initial guess of the image u 0 . The initial PSF h 0 was defined as a matrix of the same size with the Gaussian blur kernel and the value of one. Table 1 summarizes the simulation methods, the SSIM and PSNR of the degraded images and the decomposed images processed with RLBD-SATV, RLBD-SAFPDE and RLBD-HRATV, respectively. All results are shown when the NMSE error reaches the minimum. It can be seen that the proposed hybrid method achieves the highest SSIM and PSNR among all three methods. The details of the Barbara image are presented in Figure 3 as an example. The degraded image was obtained by blurring the original image with a Gaussian kernel of (7,3) and adding Poisson noise with the maximum mean intensity of 1e4, as shown in Figure 3a. For simplicity, we denote three RLBD methods using the regularization term only as SATV, SAFPDE, and HRATV, respectively. It can be seen that the edges were restored and the noise was reduced after deconvolutions. Perceptually, the image restored by SATV seems sharper than that restored by SAFPDE and HRATV, as TV regularizer favors discontinuity and often leads to false edges that do not exist in the original image. As shown in the close-ups of the face of Barbara, the degraded image was severely blurred such that the contours of the nose and the eyes cannot be resolved, but were restored after deconvolution. However, the SATV method generates false edges on the nose and the cheek, as TV tends to find piece-wise constant solutions in the smooth region. Here, the smooth region means the region in which intensities change gradually. For comparison, the SAFPDE method recovers the transition from dark to bright on the cheek, but the edge between hair and scarf is not as clearly resolved as the SATV method. The hybrid method finds a good combination of the preservation of the edges and the removal of the staircase effect Photonics 2021, 8, x FOR PEER REVIEW 11 of 21 To further show the details, the intensity profile along the nose of Barbara was extracted for analysis. From Figure 3b, it is seen that the profile restored by SATV shows a stair-shape on the rising slope, as indicated by the arrow. Furthermore, the sharp peaks are blocked due to the staircase effect. In contrast, SAFPDE restores smooth slopes and sharp peaks. HRATV also removes the staircase effect, and recovers the peaks better than SATV. To compare the edge preservation ability of the deconvolution methods, the closeups of a desk leg in the Barbara image are shown in Figure 4. From line profiles crossing the desk leg (Figure 4b-d), it can be seen that SATV and HRATV restore the sharp edges that match well with the original profiles. In comparison, SAFPDE tends to smooth the corners, as indicated by the arrow in Figure 4c. To further show the details, the intensity profile along the nose of Barbara was extracted for analysis. From Figure 3b, it is seen that the profile restored by SATV shows a stair-shape on the rising slope, as indicated by the arrow. Furthermore, the sharp peaks are blocked due to the staircase effect. In contrast, SAFPDE restores smooth slopes and sharp peaks. HRATV also removes the staircase effect, and recovers the peaks better than SATV. To compare the edge preservation ability of the deconvolution methods, the close-ups of a desk leg in the Barbara image are shown in Figure 4. From line profiles crossing the desk leg (Figure 4b-d), it can be seen that SATV and HRATV restore the sharp edges that match well with the original profiles. In comparison, SAFPDE tends to smooth the corners, as indicated by the arrow in Figure 4c.
The convergence and the dependence on regularization parameters were investigated, as the choice of the parameters is important to achieve the optimal results. In the case of image Barbara, it takes 220, 209 and 306 steps for SATV, SAFPDE and HRATV to converge, respectively. The HRATV method achieves the lowest NMSE error. It can be seen from Figure 4e that all methods need to be stopped properly to obtain the best performance. The HRATV method is more robust to the variations of the regularization parameter λ 1 or λ 2 compared with the SATV or SAFPDE method, as shown in Figure 4f,g.
From the above results, we conclude that the SATV method preserves sharp edges and avoids amplification of noise, while the SAFPDE method retains smooth regions and peaks. The HRATV method inherits both merits in image restoration and simultaneously recovers edges and fine details. From Table 1, it can also be seen that the HRATV achieves the best performance among all three methods. In OAT, images usually contain not only edges, but also fine textures. Thus, it is desirable to use the hybrid regularization in OAT to improve both the resolution and the SNR. The convergence and the dependence on regularization parameters were investigated, as the choice of the parameters is important to achieve the optimal results. In the case of image Barbara, it takes 220, 209 and 306 steps for SATV, SAFPDE and HRATV to converge, respectively. The HRATV method achieves the lowest NMSE error. It can be seen from Figure 4e that all methods need to be stopped properly to obtain the best performance. The HRATV method is more robust to the variations of the regularization parameter or compared with the SATV or SAFPDE method, as shown in Figure 4f,g. From the above results, we conclude that the SATV method preserves sharp edges and avoids amplification of noise, while the SAFPDE method retains smooth regions and peaks. The HRATV method inherits both merits in image restoration and simultaneously recovers edges and fine details. From Table 1, it can also be seen that the HRATV achieves the best performance among all three methods. In OAT, images usually contain not only edges, but also fine textures. Thus, it is desirable to use the hybrid regularization in OAT to improve both the resolution and the SNR.

Simulated OAT Images
In this section, we first investigated the effects of frequency bandwidth and spatial directivity of the transducer on OAT. Then the performance of the RL deconvolution methods was analyzed via simulated phantom data sets. Figure 5b presents the reconstructed PSFs of the imaging system with different bandwidth and directivities. The simulation setup is illustrated in Figure 5a. It can be seen that

Simulated OAT Images
In this section, we first investigated the effects of frequency bandwidth and spatial directivity of the transducer on OAT. Then the performance of the RL deconvolution methods was analyzed via simulated phantom data sets. Figure 5b presents the reconstructed PSFs of the imaging system with different bandwidth and directivities. The simulation setup is illustrated in Figure 5a. It can be seen that the ideal setup (point detector and non-limited bandwidth) achieves the highest resolution and SNR, while the limited directivity and bandwidth tends to blur the PSFs. The maximum intensity projections (MIP) along y-axis and x-axis are shown in Figure 5c,d, respectively. The limited bandwidth leads to broadening of the pulse width in the time domain and loss of signal in the frequency domain, thus the edges are smoothed and the SNR decreases. Point detectors are ideal for image formation. However, in practice, a tradeoff between the resolution and the detection sensitivity must be considered. Generally, one wants to increase the detection sensitivity by increasing the element size, which will unfortunately degrade the lateral resolution. As shown in Figure 5b, a transducer array with limited bandwidth and spatial directivity achieves the worst resolution and SNR. It has been shown that capacitive micro-machined ultrasound transducers (CMUTs) are more suitable for optoacoustic applications [42]. However, they are not widely available and involve complex fabrication procedures in contrast to the well-established PZT transducers. It is also possible to approximate the point detector by increasing the element density. However, it will increase the burden to the hardware. Within the hardware limits, image deconvolution can be a promising method to improve the image quality. maximum intensity projections (MIP) along y-axis and x-axis are shown in Figure 5c,d, respectively. The limited bandwidth leads to broadening of the pulse width in the time domain and loss of signal in the frequency domain, thus the edges are smoothed and the SNR decreases. Point detectors are ideal for image formation. However, in practice, a trade-off between the resolution and the detection sensitivity must be considered. Generally, one wants to increase the detection sensitivity by increasing the element size, which will unfortunately degrade the lateral resolution. As shown in Figure 5b, a transducer array with limited bandwidth and spatial directivity achieves the worst resolution and SNR. It has been shown that capacitive micro-machined ultrasound transducers (CMUTs) are more suitable for optoacoustic applications [42]. However, they are not widely available and involve complex fabrication procedures in contrast to the well-established PZT transducers. It is also possible to approximate the point detector by increasing the element density. However, it will increase the burden to the hardware. Within the hardware limits, image deconvolution can be a promising method to improve the image quality. To show the performance of the deconvolution methods in OAT, a vessel phantom was simulated in k-Wave. The resolution estimated by the FRC measurement is used to observe the iteration progress. The deconvolution is terminated when the resolution stops increasing. The original reconstructed image is used as the initial guess of the image . A Gaussian kernel with the size of 25 × 25 is used as the initial guess of the PSF ℎ . To avoid over deconvolution, the maximum iteration step is set to 50. The SNR was defined as: To show the performance of the deconvolution methods in OAT, a vessel phantom was simulated in k-Wave. The resolution estimated by the FRC measurement is used to observe the iteration progress. The deconvolution is terminated when the resolution stops increasing. The original reconstructed image i is used as the initial guess of the image u 0 . A Gaussian kernel with the size of 25 × 25 is used as the initial guess of the PSF h 0 . To avoid over deconvolution, the maximum iteration step is set to 50. The SNR was defined as: where µ s is the mean value of the signals and σ n is the standard deviation of the background noise. The deconvolution results of the three methods are shown in Figure 6. The vessel edges are clearly resolved by the SATV method, but the pixel intensities of the vessel appear to be piece-wise constant due to the staircase effect. SAFPDE recovers the peaks, but starts to fit the background noise after few iterations. Therefore, it achieves lower SNR (45.53 dB) than SATV (46.14 dB) and HRATV (46.37 dB). HRATV achieves comparable SNR with SATV, but higher resolution, because it restores the sharp peaks better. Since this image does not contain large smooth regions, the staircase effect caused by SATV is not as strong as that in the simulated test images. However, it can still be found on the vessels where the peaks are blocked. The FRC curves after deconvolutions are presented in Figure  6b with the resolution indicated by the intersection point with the 1/7 threshold (dashed line). The resolution obtained by the SAFPDE method nearly hits the limit of the spatial frequency range, since SAFPDE starts to fit the background noise after iterations. Figure 6c shows the line profiles along the vessel indicated by the arrow in Figure 6a. It can be seen that SATV tends to block the peaks, but produces sharp edges. HRATV finds a solution combining both the peak and the edge restoration. Figure 6d demonstrates the feasibility of using the FRC estimated resolution to monitor the deconvolution progress. In this case, SAFPDE and HRATV take seven steps to converge, while SATV takes six steps to converge.

In vivo OAT Experiments
In vivo OAT imaging of nude mice was done to test the performance of three RL deconvolution methods on real data. The deconvolution results of the liver are shown in Figure 7. The vessel indicated by the red arrow is chosen to calculate the SNR. The results are shown beside the image titles. From Figure 7a, it was clearly seen that the blood vessels and the multi-layer skin structures are enhanced after deconvolutions. However, as shown in the close-ups, the SATV method generates piece-wise constant effect on the ves-

In Vivo OAT Experiments
In vivo OAT imaging of nude mice was done to test the performance of three RL deconvolution methods on real data. The deconvolution results of the liver are shown in Figure 7. The vessel indicated by the red arrow is chosen to calculate the SNR. The results are shown beside the image titles. From Figure 7a, it was clearly seen that the blood vessels and the multi-layer skin structures are enhanced after deconvolutions. However, as shown in the close-ups, the SATV method generates piece-wise constant effect on the vessels and layered structures. This effect produces perceptually sharp edges, but actually obscures the fine details and deteriorates the resolution. In contrast, both SAFPDE and HRATV methods restore the liver structures naturally instead of generating false edges. The FRC measurements (Figure 7b) further verify the observed results, as the estimated image resolution after deconvolutions exceeds the original resolution by a large margin. In the meantime, SAFPDE and HRATV achieve higher resolution than SATV. The FRC curve of SAFPDE already saturates as it reaches the maximum detectable frequency (1/(2 × pixel size)). As a result, the SAFPDE method starts to fit the background noise after a few iterations. It can be seen that the noise level of the background in SAFPDE is higher than that in SATV and HRATV. Therefore, even SAFPDE recovers the peak signals better, it achieves lower SNR than the other two methods. The line profile of the vessels indicated by the yellow arrow was extracted from the three deconvolution results for analysis ( Figure 7c). The SAFPDE method is able to recover the sharp peaks, but slightly smooth the edges. On the contrary, SATV flattens the peaks, but generates sharp edges. HRATV keeps the sharp edges and restores the peaks by combing both SATV and SAFPDE regularization. Figure 7b shows the progress of the estimated resolution during iterations. The three deconvolution methods converge at 11 (SATV), 9 (SAFPDE) and 10 (HRATV) steps, respectively.
The deconvolution results of the mouse brain are presented in Figure 8. Different from the liver image that contains dense blood vessels and smooth regions, the brain image has mainly micro-structures. The target vessel (indicated by the red arrow) was selected to calculate the SNR. After deconvolution, the SNR is significantly improved by 9.42 dB (SATV), 6.53 dB (SAFPDE), and 12.69 dB (HRATV), respectively. From the close-ups, it can be found that the structures indicated by the yellow arrow cannot be discriminated from each other in the original image, while they were clearly separated after deconvolutions, especially in the images obtained by SAFPDE and HRATV. This observation is further verified by inspecting the line profiles, as indicated by the arrow in Figure 8c. The FRC measurements show that the SAFPDE and HRATV methods achieve higher resolution than the SATV method, which agrees with the deconvolution results of the simulated vessel phantom and the mouse liver. The brain image is similar to the simulated vessel phantom as both of them are sparsely represented and barely contain smooth regions. Thus, the staircase effect is not as strong as that in the liver image, but still can be observed on the structures shown in the close-ups. The deconvolution results demonstrate that the HRATV method is able to combine both texture restoration and noise suppression, and achieves the highest SNR among all three methods. Figure 8d shows that the iterations stop at 13 (SATV), 11 (SAFPDE) and 16 (HRATV) steps, respectively. by the yellow arrow was extracted from the three deconvolution results for analysis (Figure 7c). The SAFPDE method is able to recover the sharp peaks, but slightly smooth the edges. On the contrary, SATV flattens the peaks, but generates sharp edges. HRATV keeps the sharp edges and restores the peaks by combing both SATV and SAFPDE regularization. Figure 7b shows the progress of the estimated resolution during iterations. The three deconvolution methods converge at 11 (SATV), 9 (SAFPDE) and 10 (HRATV) steps, respectively. The deconvolution results of the mouse brain are presented in Figure 8. Different from the liver image that contains dense blood vessels and smooth regions, the brain im- The deconvolution results of the mouse kidney are shown in Supplementary Figure S1 to further showcase the performance of the proposed method in post-processing of OAT images. The improvement after deconvolution in OAT is not as significant as that in simulated test images. The reasons are two-fold. First, a uniform Gaussian kernel with known size and standard deviation is applied on the test image to simulate the blur. However, the practical OAT systems are more complex, as the PSF of the system is spatially variant and the size of the PSF is unknown. Furthermore, the blur kernel is not uniform Gaussian as the axial and the lateral resolution come from the transducer's frequency bandwidth and spatial directivity, respectively. Therefore, assuming a spatial-invariant, Gaussian blur kernel does not match well with the real scenario, which also inspires us to introduce the spatial dependency to the estimation of the blur kernel in future work. Second, the noise presented in the OAT system often comes from multiple sources, instead of single Poisson or Gaussian noise. The acoustic medium, the transducer and the receiving circuits all contribute to the noise imposed on the detected signals [43]. Thus, it can be beneficial to incorporate the mixed noise removal technique in the OAT deconvolution method. Figure 8c. The FRC measurements show that the SAFPDE and HRATV methods achieve higher resolution than the SATV method, which agrees with the deconvolution results of the simulated vessel phantom and the mouse liver. The brain image is similar to the simulated vessel phantom as both of them are sparsely represented and barely contain smooth regions. Thus, the staircase effect is not as strong as that in the liver image, but still can be observed on the structures shown in the close-ups. The deconvolution results demonstrate that the HRATV method is able to combine both texture restoration and noise suppression, and achieves the highest SNR among all three methods. Figure 8d shows that the iterations stop at 13 (SATV), 11 (SAFPDE) and 16 (HRATV) steps, respectively. One of the long-standing questions in iterative image deconvolution is when to stop the iteration and how to monitor the image quality in the absence of ground truth. The most used indicator is the root mean square (rms) error between two successive iterations and the iteration is terminated when the rms error drops to a pre-defined threshold [26,29,31,33,34]. The threshold is usually found by trial and cannot quantitatively measure the image quality. This stopping criterion is also widely adopted in deconvolution of optoacoustic images. It was shown that the deconvolution has to be stopped early to avoid excessive iteration because the rms error may have oscillations [28]. In this work, we showcased the feasibility of using FRC measurements to observe the iteration progress in OAT deconvolution. The iteration can be stopped when the estimated resolution stops increasing. Furthermore, the resolution estimated by FRC also provides a priori to the initial guess of the image PSF h 0 . In the deconvolution of the simulated and experimental OAT images, the relationship between the standard deviation σ of the initial Gaussian kernel and the estimated resolution R by FRC is expressed as: In addition to iteration steps, regularization parameters are key to the optimal deconvolution performance as well. In the proposed RLBD-HRATV method, there are three parameters to tune: λ 1 , λ 2 and c. λ 1 and λ 2 are parameters to adjust the strength of the SATV and the SAFPDE, separately. Generally, λ 1 and λ 2 can be found in the range of [10 −5 , 1] by trial and error, tuning the parameter value and observing the image quality. It is laborious to search for a suitable value in such a wide range. In this work, we chose the parameter that provides the highest SNR. The drawback of this method is that the signal and noise regions need to be pre-defined in order to calculate the SNR. We found that λ 1 , λ 2 ∈ 10e −4 , 0.1 is generally appropriate for OAT deconvolution. The larger value of the parameter indicates stronger regularization, which means the noise will be well suppressed at the cost of undermining the edges and textures. Therefore, a trade-off must be tuned to achieve optimal SNR. As demonstrated on the Barbara image, SATV regularization favors discontinuities and generates blocky solutions on slopes, while the SAFPDE filter produces smooth images and restores details. Typically, OAT images contain the mixture of edges, regions with smooth change in intensity, and textures, which encourages the use of hybrid regularization on the whole image. However, if only patches containing mainly edges or textures like rich blood vessels are of particular interest, the regularization of SATV or SAFPDE can be reduced correspondingly. The parameter c is used to tune the combination of the SATV and SAFPDE regularizations. It is image dependent and should be chosen such that the weighting matrix θ approaches zero in the smooth region. Choosing zero value of c means that the hybrid method RLBD-HRATV shrinks to RLBD-SATV and the edges of the image will be largely preserved, while pushing c to one indicates that the image is dominated by smooth regions. In this work, we found a proper c can be found in the range of [0.05, 0.4] for OAT when |∇u k | is normalized to [0, 1]. The parameter c is important to the deconvolution performance, as it affects the results of edge preservation and staircase removal.

Conclusions
In this work, we proposed a modified blind RL method tailored for deconvolution of OAT. The modified method incorporates a hybrid regularization term (HRATV) combining the spatially adaptive total variation (SATV) and the fourth-order partial differential equation (SAFPDE) to the RL iteration function. The motivation is to recover both the fine details and sharp edges in OAT. The proposed method was validated on a wide range of images from standard test images to simulated OAT data and in vivo experimental data. It was shown that the RLBD-HRATV method not only produces the best deconvolution results, but performs more robustly to the variation of regularization parameters. The blood vessels and micro-structures are clearly enhanced after deconvolutions. The hybrid method is able to restore the fine textures of OAT images by reducing the staircase effect and improve the SNR by achieving a trade-off between edge preservation and noise suppression. In addition, we applied the FRC measurements to evaluate the deconvolution progress in OAT and provide quantitative evaluation of the image quality. Based on the results, we provided a practical guide on the initialization procedure and selection of regularization parameters in deconvolution of OAT using the proposed RLBD-HRATV method. Future work will focus on incorporating the spatial variant PSF of the OAT system and mixed noise removal technique to deconvolutions.
Author Contributions: Conceptualization, C.Y.; methodology, C.Y. and Y.J.; validation, C.Y.; formal analysis, C.Y. and Y.J.; investigation, C.Y., and Y.J; resources, Y.C. and X.J.; writing-original draft preparation, C.Y.; writing-review and editing, X.J. and Y.C.; supervision, Y.C.; project administration, Y.C.; funding acquisition, X.J. and Y.C. All authors have read and agreed to the published version of the manuscript. Institutional Review Board Statement: The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of Suzhou Institute of Biomedical Engineering and Technology, Chinese Academy of Sciences (protocol code: 2020-A23 and date of approval: 15 October 2020).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.