Convolutional sparse coding for compressed sensing photoacoustic CT reconstruction with partially known support

In photoacoustic tomography (PAT), imaging speed is an essential metric that is restricted by the pulse laser repetition rate and the number of channels on the data acquisition card (DAQ). Reconstructing the initial sound pressure distribution with fewer elements can significantly reduce hardware costs and back-end acquisition pressure. However, undersampling will result in artefacts in the photoacoustic image, degrading its quality. Dictionary learning (DL) has been utilised for various image reconstruction techniques, but they disregard the uniformity of pixels in overlapping blocks. Therefore, we propose a compressive sensing (CS) reconstruction algorithm for circular array PAT based on gradient domain convolutional sparse coding (CSCGR). A small number of non-zero signal positions in the sparsely encoded feature map are used as partially known support (PKS) in the reconstruction procedure. The CS-CSCGR-PKS-based reconstruction algorithm can use fewer ultrasound transducers for signal acquisition while maintaining image fidelity. We demonstrated the effectiveness of this algorithm in sparse imaging through imaging experiments on the mouse torso, brain, and human fingers. Reducing the number of array elements while ensuring imaging quality effectively reduces equipment hardware costs and improves imaging speed.


Introduction
Photoacoustic imaging (PAI) technology has flourished in recent years.[1][2][3][4][5][6] PAT stimulates the photoacoustic signal in the imaging region with wide-field illumination and obtains an image of the target area through inverse reconstruction.It can quantify blood oxygen saturation by measuring the concentration of deoxyhemoglobin (HbR) and oxyhemoglobin (HBO2).PAT combines high ultrasonic resolution and strong optical contrast in a single modality, providing a safe, efficient and inexpensive method for functional imaging with imaging depths of up to several centimetres [7].PAT has been effectively utilized to the imaging of in vivo vascular structures, functional hemodynamic imaging, and the visualization of breast [8,9] and brain tumors [10][11][12], among others.It holds tremendous promise for clinical diagnosis [13].
PAT's implementation modes include linear array, circular array, and convex array, among others.Each mode has its own application-specific scenario [14].This article addresses the issue of high-quality imaging using sparse circular array.That is, a set of sparsely distributed ultrasonic transducers are positioned around the imaging target to detect the photoacoustic signal in the circular area and reconstruct the image using the inverse algorithm.In practical system applications, a comprehensive information collection within the imaging area requires a dense distribution of array elements.Single element scanning imaging systems require hundreds or even thousands of scans.In addition, it places greater demands on the data collection system, which substantially increases the costs of the equipment.Sparse imaging has therefore become a research priority in the subject of array photoacoustic imaging, which necessitates the use of a small number of array elements to collect photoacoustic signals in the imaging area and perform high-quality image.Under the sparse sampling condition, traditional methods which include filtered back-projection (FBP) reconstruction algorithm will generate under-sampling artifacts, such as streaking artifacts or grating lobes [15].
To recover the original signal or image from under-sampled measurements, numerous sparse reconstruction algorithms have been developed.There are numerous classic iterative reconstruction algorithms that can tackle sparse reconstruction problems, such as the algebraic reconstruction algorithm (ART), [16] the joint algebraic reconstruction algorithm (SART), etc [17].In recent years, the theory of CS has swiftly evolved to restore the original signal from a sparse signal effectively and precisely [18].Behind its complex mathematical expressions, CS theory contains very subtle notions.CS has generated an enchanted outcome by violating the Nyquist sampling law on the basis of a creative concept supported by rigorous mathematical proof [19,20].During the course of signal sampling, using a small number of sampling points achieves the same effect as complete sampling.
The DL-based image reconstruction algorithm is a significant branch of computer science that has produced numerous results [21][22][23].However, the majority of DL-based methods typically learn offset versions of features that contain the same features, which will produce related artefacts in the reconstructed image, and learning dictionaries has a substantial effect on the imaging results [24].Liu et al. proposed an advanced reconstruction model using the K-VSD dictionary learning technique and adapt the model into the 3D PACT system [22].Brendt Wohlberg proposed a new solution to the problem of high computational complexity in Convolutional Spatial Computing (CSC) in the Fourier domain with a novel method [25].The concept of CSC was initially introduced by Zeiler et al [26].Peng Bao et al. introduced CSC to CT reconstruction in 2019 [27].Considering that the image will be sparser in the gradient domain, we add the gradient domain constraint to the model, which we refer to as CS-CSCGR.The theory of partial known support (PKS), which alludes to the position of non-zero pixels in the sparse transform domain, was recently introduced in CS.Recently, N. Vaswani et al. devised a modified-CS approach where an extra prior beyond sparsity was used: partially known support (PKS), that is, some known nonzero positions in the transformed sparse domain [28].PKS can further reduce the amount of measurements required for high-fidelity reconstruction, a considerable benefit that PACT is also very interested in, according to theoretical studies and applications in MRI.
Our work introduced PKS into the CS-CSCGR reconstruction framework and enhanced the reconstruction effect.Simulation and experimental data indicate that CS-CSCGR-PKS is superior to other image reconstruction algorithms due to its ability to achieve high-quality photoacoustic image reconstruction under sparse under-sampling and converge quickly after a small number of iterations.These benefits of CS-CSCGR-PKS can enhance imaging speed, reduce the cost of devices, and be utilized in a variety of biomedical applications.

Penalized least squares for circular-array PA
The following equation describes the formation and transmission of photoacoustic signals during the photoacoustic imaging process [29].
where p(r, t) refers to the photoacoustic pressure signal at the position r and time t, δ(t) is the pulsed laser excitation, p 0 (r) refers to the initial acoustic pressure, and c refers to the sound velocity.We can obtain Eq. ( 2) according to Eq. (1).
where r ′ represents the location of the area of interest.
The measurement matrix K of the ring array photoacoustic tomography system is shown below.
where r i,j indicates coordinate values of the region of interest, r m and p indicates the position and the total number of the ultrasonic transducer, and q s indicates the sampling point of the photoacoustic signal.The basic CS reconstruction algorithm model is as follows [30,31].
where y indicates the photoacoustic signal measurements, u indicates the photoacoustic image, K refers to the measurement matrix, β refers to a hyperparameter, and R(u) indicates the regularization term.

Convolutional sparse coding with partially known support
The CSC-based method considers the reconstructed image to be the sum of convolutions of the among a group of filters and the corresponding feature map [27].
Where u indicates the photoacoustic image to be reconstructed Let M i ∈ R N be the feature maps corresponding to filters f i , and the support is where T i0 ⊂ {1, 2, . . ., N} is the known support of T i and ∆ = {1, 2, . . ., N} is the unknown support of T i .The introduction of PKS can constrain feasible solutions to a smaller range, thus producing accurate reconstructions with fewer samples.Therefore, the CS optimization process of PKS is shown below.
where M i∆ indicates the feature map based on the known support and ε is the allowable noise level.

Method
As demonstrated below, CSC and CSC-PKS were devised for circular array PAT based on the CS and PKS theories described previously.
where λ is the regularization parameter.
It is preferable to impose gradient constraints on the feature maps rather than the image domain constraints.Next, we enhance the optimisation problem by imposing additional constraints on CS-CSC-PKS via gradient regularisation on the characteristic graph.
where g 0 and g 1 are the gradient-calculating filters across the axes of x and y, respectively, and τ acts as a hyperparameter.In this work, the filters {f i } were trained by the full-sampled photoaccoustic dataset.So Eq. ( 9) becomes the new problem: We can use alternating-direction multiplier method (ADMM) to solve the Eq. ( 1)0.The first step is to use a group of initialized feature maps {M i } to obtain the intermediate reconstructed image u.Equation (1)0 becomes: [27] arg min The second step is to use a fixed filter {f i } to represent the reconstructed image u obtained in the previous calculation, which means calculating {M i }.Then we update the optimization equation.

arg min
The definition of M i∆ in Eq. (1)0 is given by: To solve Eq. ( 10), we introduce a matrix W i whose elements are set 1 when the corresponding gradient domain feature image pixels belong to the unknown support ∆ in the (k − 1) th iteration, and set 0 otherwise.Equation ( 12) is rewritten as.arg min We refer to W i as a partial support matrix by T i0 , which can be obtained by setting a threshold.This paper sets the position greater than a certain value on the gradient domain feature map as a known support.
where x (k) j is the j th element of the feature map M i , k refers to the total number of repetitions and ξ refers to the threshold used to obtain W i .The iterative restoration procedure can function well for rapidly degrading feature map M (k)  i by setting , where η is a hyperparameter.Next, we will further simplify. Let Similarly, we can obtain G l M i = g l * M i and update the gradient part at Eq. ( 14).14) could be revised as below.
arg min where ADMM solves Eq. ( 14) through the introduction of a dual variable B, the procedure is outlined below [32].
where ⊙ is the Hadamard product and then we can solve the Eq. ( 19) as follows.
where ρ is a hyperparameter.The overall CS-CSCGR-PKS algorithm for circular-array PACT reconstruction is shown below.

Results
To verify the benefits of the method proposed in this article, we first perform experimental verification using simulation data.The simulation environment is configured with 256 array elements uniformly distributed around the 10 centimeters diameter imaging area, a sampling frequency of 25 MHz, and a sound velocity of 1540 meters per second.We utilize an open-source photoacoustic image database of the human brain.The dataset consists of ten distinct digital photoacoustic images of the human brain, as depicted in Fig. 1.

Fig. 1. Digital human brain data
We utilize all 256 array components for full sampling, while the sparsity is 8.As a result, sparse reconstruction is dependent on the photoacoustic signals gathered by 32 array components that are evenly dispersed around the perimeter of the imaging region.In this paper, the peak signal-to-noise ratio (PSNR), root mean square deviation (RMSE), and structural similarity (SSIM) are employed to evaluate the algorithm's performance.Initially, this data set was used to train predetermined filters.After training the filter using a photoacoustic digital image of the human brain, 32 trained filters of size are obtained.As shown in Fig. 2. We evaluate our approach in comparison to a number of well used and successful methods, such as FBP, CS-DL, and CS-CSCGR.Figure 3 depicts both the original photoacoustic digital image of the human brain and the images that were restored through the aforementioned techniques.The images that are produced by the FBP reconstruction technique contain a large number of artefacts, which makes it challenging to employ them in a practical setting.Clear photoacoustic pictures that have a certain degree of detail restoration can be obtained through the use of sparse reconstruction techniques that are based on CS.When compared to CS-DL, both CS-CSCGR and CS-CSCGR-PKS improve the quality of the restored image, make the image sharper, and minimize the amount of speckle artefacts.It is difficult to tell the difference between CS-CSCGR and CS-CSCGR-PKS with the human eye; nevertheless, the absolute difference image (depicted in Fig. 4) reveals that CS-CSCGR-PKS demonstrates superior performance.achievable.The results of FBP reconstruction are essentially unavailable, and microscopic blood vessels in images reconstructed using CS-DL are nearly nonexistent.Figure 5(d) and (e) show that both CS-CSCGR and CS-CSCGR-PKS have outstanding sparse reconstruction performance, preserving the majority of image details, and CS-CSCGR-PKS achieves superior imaging.Simultaneously, we performed a quantitative analysis of the imaging outcomes.Table 1 displays the quantitative results of digital photoacoustic images of the human brain.We discovered that CS-CSCGR-PKS is the best of the three indicators.The iterative steps of the reconstruction algorithm based on compressive sensing are uniformly set to 50.The small animal imaging system used in this article (Fig. 6) is assembled based on a programmable ultrasound imaging system (Prodigy, S-sharp, Hong Kong, China), and the system schematic is shown in Fig. 6.The system mainly consists of a pulse laser module, a multi-channel data acquisition module, and a signal processing and timing control module.Convert the initial sound pressure signal into a photoacoustic image and test the algorithm proposed in this article.Concurrently, to demonstrate the efficacy of the algorithm, it was applied to photoacoustic imaging of the mouse torso.The imaging system is depicted in Fig. 6.Considering the comparatively deep penetration depth of 1064 nm laser, we employ 1064 nm laser to excite photoacoustic signals in subsequent studies.Using a laser with a wavelength of 1024 nm and a repetition rate of 20 Hz, we measured that the laser energy distribution density of the tissue specimen was 14.1 mJ/cm 2 , which is within the ANSI safety standard.An eight-way optical fiber equitably divides the laser beam into eight beams, and the annular irradiation of the imaging target area is achieved by optical path shaping.At the output end of the optical fiber, each laser beam should be designed as a narrow and elongated beam spot.Each fiber-emitted laser is collimated into a parallel beam through a convex lens (Thorlabs LA1131-YAG, f = 50.8)and then converged in a single direction through a flat cylindrical lens (Thorlabs LG1131-LU, f150) to form a high energy narrow beam spot.The circular device contains 256 array elements to collect photoacoustic signals.The imaging area is a 10-centimeter-diameter circular region.In the context of photoacoustic imaging, a multi-wavelength tunable optical parametric oscillation (OPO) pulse laser (SpitLight-600, Innolas, Munich, Germany) was employed (Fig. 7(c)).Additionally, a mouse clamping mechanism was devised with sole purpose of facilitating imaging of the mouse (Fig. 7(b)).In order to optimize the imaging quality, it is imperative to perform debugging procedures to verify the alignment between the cross-sectional profile of the ring array photoacoustic signal acquisition and the cross-sectional profile of the ring light path illumination.The usage of a specialized motion controller facilitates the axial displacement of the circular array probe and the imaging target, hence simplifying the processes of debugging and three-dimensional imaging.
We use the FBP, CS_DL, CS_CSCGR, and CS_CSC_PKS imaging algorithms for imaging based on 32 elements (1/8 of total 256 elements), and the label image is a fully sampled photoacoustic mouse torso image.The imaging results are depicted in Fig. 8. Figure 8(a) shows the position of the nude mouse to be imaged and the imaging cross-section.Figure 8(b)-8(c) show the images restored by the BP algorithm with 256 and 32 detecting elements, respectively.Under-sampling artifacts appear in the Fig. 8(c).Figure 8(d)-(f) show the images reconstructed by CS algorithm.They can greatly suppress the generation of under sampling artifacts and improve imaging quality.If a more accurate system matrix is used, the image reconstruction results of the CS method can be further improved.Based on the CS_CSCGR and CS_CSCGR_PKS reconstruction algorithms, we find that high-quality photoacoustic images can be obtained.Nonetheless, based on the absolute difference image (illustrated in Fig. 9), it is demonstrated that our algorithm produces the greatest quality reconstructed image.
To demonstrate the superiority of CS-CSCGR-PKS, the results of the mouse thoracic imaging are quantitatively analyzed in Table 2.We discovered that CS-CSCGR-PKS had the maximum score for each of the three metrics.It is consistent with photoacoustic digital imaging of the human brain.By comparing Fig. 9(a)-(d), it can be found that the light body image reconstructed based on the algorithm proposed in this paper is the closest in pixel intensity distribution to the light acoustic image reconstructed based on full ring array elements.Finally, the goal of achieving a sparse array imaging effect that is close to that of a full ring array imaging effect has been achieved.The second in vivo experiment was conducted on blood vessels in the brain of nude mouse.The imaging results are depicted in Fig. 10. Figure 10(a) shows the brain of the nude mouse.From the Fig. 11, it can be observed that the imaging results of the mouse brain are similar to mouse thoracic, based on CS_ CSCGR_ The photoacoustic image reconstructed by the PKS algorithm is closest to the label image.The absolute error image shows that the algorithm proposed in this paper has advantages in detail restoration, and compared to CS_ CSCGR requires fewer iteration steps to achieve the best imaging results.
We also conducted quantitative analysis on brain imaging of nude mouse.To demonstrate the superiority of CS-CSCGR-PKS, the results of mouse brain imaging are quantitatively analyzed in Table 3.We first calculate the histogram of the differential image as shown in Fig. 12, and label is the reconstructed image based on 256 elements.Figures 12(A)-(D) show the difference histograms between the reconstructed image and the label image using data from 32 transducer components by BP, CS-DL, CS-CSCGR, and CS-CSCGR-PKS, respectively.We find that the error pixels of   the reconstructed image based on CS-CSCGR-PKS are basically in the range of 0∼0.1, so it is more accurate.From Figs. 12(E)-(H), we can find that, for the images reconstructed with CS, the majority of the difference amplitudes is within the lower ranges.Figures 12 (I)-(L) show the similar results.
The imaging of human joints is a significant application of photoacoustic circular array imaging.Through the use of external contrast agents, it is able to perform morphological, microvascular, and functional imaging of joints, as well as molecular imaging.We also apply it to finger photoacoustic imaging.The laser beam is divided equitably into eight beams by eight optical fibers, and circular illumination of the imaging target area is achieved via optical path shaping.The collection of photoacoustic signals is carried out by an array of ultrasonic transducers arranged in a loop around them.Using a laser with a wavelength of 1024 nm and a repetition rate of 20 Hz, we measured that the laser energy distribution density of the tissue specimen was 15.2 mJ/cm 2 , which is within the ANSI safety standard.By installing a motor on the z-axis, a circular array can accomplish three-dimensional scanning.By scanning the left middle finger 100 times, we were able to obtain subcutaneous vascular imaging.
Similarly to photoacoustic digital human brain imaging, the BP and CS-CSCGR-PKS methods are used to reconstruct the finger image of each slice.After processing all cross-sectional images, a three-dimensional image of finger blood vessels is obtained by superimposing all scanned along the axial images of the finger.Figure 13 (B) depicts a 3D image of the finger's vasculature reconstructed using the BP method based on 256 elements.Figure 13 (C) and (D) depict 3D vascular images of digits reconstructed with BP and CS-CSCGR-PKS using 32 elements (1/8 of the total of 256 elements).To demonstrate the high-quality imaging capability of the CS-CSCGR-PKS method under scarce sampling conditions, we choose two out of one hundred representative images to illustrate the reconstruction effect.
The above results demonstrate that the CS-CSCGR-PKS algorithm proposed in this article is efficacious.Following is a quantitative examination of the outcome of the appeal.
Further, a local comparison is made by drawing a line along the cross section of the finger image.Figure 14 shows that the amplitude details of the photoacoustic image restored based on CS-CSCGR-PKS are closer to the target image.Figure 14

Discussion
In the field of rapid photoacoustic imaging, applying CS to sparse imaging is a highly efficient technique.In order to accomplish high-quality image restoration in the context of the sparse acquisition of photoacoustic signals in circular array photoacoustic tomography systems, we have designed a novel and effective sparse reconstruction architecture for photoacoustic images.CS-CSCGR-PKS applies convolutional sparse encoding in the gradient domain of images and introduces PKS.Compared to DL's reconstruction algorithm, it can avoid plaque aggregation artefacts, and gradient domain constraints can further suppress artefact generation.In order to accelerate the convergence of the reconstruction algorithm and reduce the number of iterations, this article implements PKS in the algorithm to accelerate convergence under certain known support conditions, thereby further reducing artefacts and obtaining higher quality photoacoustic images.This algorithm can restore high-quality photoacoustic images, as demonstrated by simulation, in vivo mouse and human finger experiments, and its efficacy has been confirmed.Research has demonstrated that CS-CSCGR-PKS has promising application prospects within the field of photoacoustic sparse imaging.On the basis of the algorithm proposed in this article, high-quality photoacoustic image can be accomplished using a small number of array components, which can be used for photoacoustic signal acquisition and effectively reduce the equipment's hardware cost.

Algorithm 1 :
Sparse reconstruction algorithm using the CS-CSCGR-PKS model Input: Photoacoustic signal y , measurement matrix K , preset filter i f Output: Initial sound pressure image u

Fig. 2 .
Fig. 2. Filters trained with photoacoustic digital human brain data set

Figure 5
Figure 5 depicts the blue dashed region of Fig. 3(a) to illustrate the reconstruction details.Multiple minuscule blood vessels are represented by the white arrows in the digital human brain.On the basis of the algorithm presented in this article, the most precise reconstruction results are

Figure 10 (
b)-8(c) show the images restored by the BP algorithm with 256 and 32 detecting elements, respectively.Figure10(d)-(f) show the images reconstructed by CS algorithm.At the same time, we also presented the difference images relative to the original image (Fig.11).

Fig. 12 .
Fig. 12.The amplitude histogram of the reconstructed image's difference image.(A)-(D) The distribution histogram of pixel error images in the digital human brain; (E)-(H) The distribution histogram of pixel error images in mouse torso; (I)-(L) The distribution histogram of pixel error images in the mouse brain.

Fig. 13 .
Fig. 13.Reconstructed images of the finger.(A) Photograph of the finger; (B)-(D) 3D images restored by FBP based on 256 elements, FBP based on 32 elements and CS-CSCGR-PKS based on 32 elements; (a)-(f) B-scan images along the white dashed lines in (B)-(D).

Fig. 14 .
Fig. 14.Local comparison of light absorption intensity along a defined straight line in the cross-section image of the finger.(A) The blue dashed line; (B) The green dashed line.
(A) shows the initial sound pressure amplitude distribution of the blue dashed line, whereas Fig. 14(B) shows the initial sound pressure amplitude distribution of the green dashed line.The results show that: (1) Compared to BP, reconstructed images based on CS-CSCGR-PKS are of higher quality; (2) The CS-CSCGR-PKS reconstruction image using only 32 transducers is basically equivalent to the BP image using all ultrasound unit data.