Dose Images Reconstruction Based on X-ray-Induced Acoustic Computed Tomography

Abstract

The accurate reconstruction of the in vivo dose is a critical step in radiation therapy. X-ray-induced acoustic imaging is a promising technology for in vivo dose reconstruction, as it enables the nonradiative and noninvasive monitoring of radiation dose. However, current X-ray acoustic imaging methods suffer from several limitations, including high signal-to-noise ratio, poor imaging quality and massive loss of structural information. To address these limitations, we propose a dose image reconstruction method based on tensor sparse dictionary learning. Specifically, we combine tensor coding with compressed sensing data, extend two-dimensional dictionary learning to three-dimensional by using tensor product, and then utilize the spatial information of X-ray acoustic signal more efficiently. To reduce the artifacts of reconstruction images caused by spare sampling, we design the alternate iterative solution of the tensor sparse coefficient and tensor dictionary. In addition, we build the X-ray-induced acoustic dose images reconstruction system, simulate the X-ray acoustic signals based on patients’ information from Sichuan Cancer Hospital, and then create the simulated datasets. Compared to some typical state-of-the art imaging methods, the experimental results demonstrate that our method can significantly improve the quality of reconstructed images and the accuracy of dose distribution.

1. Introduction

X-ray acoustic imaging (XAI) is a novel biomedical imaging technique that builds up images based on acoustic signals generated by pulsed X-ray radiation. Figure 1 shows the basic schematic diagram of X-ray acoustic imaging. After X-ray excitation, irradiated tissues convert the deposited energy to heat and produce a temperature rise. The thermal expansion of tissues generates acoustic waves. The acoustic wave signals are recorded by the transducers arranged in a specified geometry, processed by the amplifier and data acquisition card, and finally form images. The basic mechanism of XAI is similar to that of photoacoustic imaging (PAI) [1,2]. PAI is a common biomedical imaging method, which is based on the transduction of light photons into acoustic waves, and has high contrast and is not easily affected by speckle artifacts [3,4,5,6,7,8]. However, due to the strong scattering of the laser in biological tissues, the depth of light-photon penetration is greatly limited. In contrast, X-ray photons travel in a straighter path and have shorter wavelengths [9]. There is less scattering and deeper penetrating in tissues. So, XAI has great potential in the field of biomedical imaging.

X-ray-induced acoustic computed tomography (XACT) has been proposed as an in vivo dosimetry method based on X-ray acoustic effects [10,11,12,13]. In radiation therapy, the X-ray acoustic signal is generated by the medical linear accelerator (LINAC), and its amplitude is proportional to the absorbed dose [1,14]. XACT uses X-ray acoustic signals detected by the sensors to reconstruct the dose distribution. Then, in vivo dosimetry can be performed. Previously, some researchers have conducted dose distribution reconstruction based on XACT in real soft tissue imitation experiments [15]. XACT has been demonstrated as an in vivo dose measurement tool during the simulation of prostate radiotherapy [16]. In the clinical application, there are only a few in vivo dosimetry tools that provide complete 3D dose distribution for patients. The most common tool is the electronic portal imaging device (EPID), which is commercially available [17,18,19,20]. However, the use of its MV photon beams leads to poor contrast and requires the delivery of additional ionizing radiation to the patient. By comparison, XACT has the advantage of measuring the dose with no radiation, noninvasively and in real time [21,22].

At present, reconstruction methods of photoacoustic images mainly include universal back projection [23], filtered back projection (FBP) [24], model-based iterative reconstruction [25] and learning-based reconstruction [26,27,28]. However, there are few methods to reconstruct X-ray-induced acoustic images: only FBP [29] and time reversal method [30]. By FBP, Kim et al. reconstructed X-ray excited acoustic images of a lead rod in a medical linear accelerator, and measured the X-ray dose distribution [1]. Wang et al. conducted real-time 3D dose images reconstruction during radiotherapy [30]. The iterative time reversal imaging technique is applied to the reconstruction of dose distribution in the simulated prostate radiotherapy plan [31]. However, these methods have some problems, such as high signal-to-noise ratio, low image quality and massive loss of the control information structure in the image.

In this paper, we present a new dose reconstruction method based on tensor sparse dictionary to improve the image quality. Our experimental datasets use 30 X-ray acoustic simulated datasets that depend on the nasopharyngeal tumor clinical cases provided by Sichuan Cancer Hospital. Compared to some state-of-the-art imaging methods, the experimental results show that our presented method can significantly enhance the quality of the reconstructed images and the accuracy of dose measurement. The contributions of this paper are as follows:

(1) Combining tensor sparse coding with compressed sensing data reconstruction theory, we propose a tensor joint sparse coding method to make more efficient use of the XACT signal spatial information.

(2) We present an alternate iterative solution of tensor sparse coefficient and tensor dictionary, which can effectively implement XACT signal data interpolation reconstruction and solve the problem of reconstruction image artifacts caused by undersampling.

(3) It is the first time that treatment planning system (TPS) information of nasopharyngeal carcinoma is used to construct the simulation workflow, which has guiding significance for clinical application.

The rest paper is organized as follows: Section 2 provides the related principles and our method in detail. In Section 3, we present the experimental datasets and results. Finally, Section 4 offers a discussion and conclusion.

2. Method

2.1. X-ray-Induced Acoustic Dose Images Reconstruction System with the Clinical Linear Accelerator

Figure 2a shows a three-dimensional model of XACT technology applied to the patient with nasopharyngeal carcinoma during radiotherapy, including the patient setup and ultrasound transducers array configuration. The equipment above the patient is the clinical linear accelerator system (LINAC) (Elekata Synergy Platform, Sweden). The red cylinder represents X-rays from the medical LINAC, radiating to the patient’s tumor area. A group of small yellow cylinders represent 2D ultrasound sensors array that are placed on the head and used to detect the generated X-ray acoustic signals. Considering the actual detection situation, the sensors are placed close to human tissue to reduce the attenuation of ultrasonic signals in the air. Figure 2b summarizes the steps taken from the X-ray excitation to the reconstruction of in vivo dose images during radiotherapy.

2.2. The Principle Background of XACT

In XACT, the duration of the X-ray pulse is shorter than the expected thermal constraint time. Therefore, in tissues and organs, the response to the spatial and temporal pulsed X-ray point sources can be used to approximately estimate the acoustic pressure $p(r,t)$ at location r and time t generated by the propagation of the acoustic waves in the homogeneous acoustic medium. The acoustic pressure $p(r,t)$ is expressed as [32]

$$\frac{{\partial }^{2}p(r,t)}{\partial t^2}-{v_s}^2{\nabla }^{2}p(r,t)=\frac{\beta {v_s}^2}{C_p}\frac{\partial S_s(r,t)}{\partial t}$$

(1)

where $\beta$ denotes the thermal coefficient of expansion, $v_s$ denotes the speed of sound, $v_s=(1/\sqrt{\rho k})$, $\rho$ denotes the density of tissues and organs, $C_p$ denotes the specific heat capacity at a constant pressure, r denotes the space coordinate position, t denotes the time and $S_s(r,t)$ is the heating source. Equation (1) can be expressed as the initial value problem of the homogeneous equation

$$\frac{{\partial }^{2}p(r,t)}{\partial t^2}-{v_s}^2{\nabla }^{2}p(r,t)=0$$

(2)

The initial conditions are as follows:

$${p(r,t)|}_{t=0}=\Gamma S_s(r,t),\frac{\partial p(r,t)}{\partial t}{|}_{t=0}=0$$

(3)

where $\Gamma$ is a Grueneisen parameter. We obtain Poisson-type integrals by using Green’s functions in free space, and then establish the analytical solution of the initial value problem in Equations (2) and (3), as follows:

$$p(r,t)=\frac{\Gamma }{4\pi v_s}\frac{\partial }{\partial t}{\int }_{S^{\prime }(r,t)}\frac{S_s(r^{\prime },t)}{|r-r^{\prime }|}d{S}^{\prime }(r,t)$$

(4)

where $S^{\prime }(r,t)$ denotes a spherical surface that is time dependent, $r^{\prime }$ denotes an absorbing point by X-ray excitation, $S_s(r^{\prime },t)$ denotes the heating source at location $r^{\prime }$ and time t, and $\left|r-r^{\prime }\right|$ denotes the distance between absorption point $r^{\prime }$ and propagation point r, and $\left|r-r^{\prime }\right|=v_{s}t$.

For the inversion method of Equation (4), the back projection method is more popular at present. If the surface S surrounding the optical absorber has a known pressure $p(r,t)$, the closed-form inversion Equation (4) is

$$S_s\left(r^{\prime }\right)=\frac{1}{\Gamma }{\int }_{\Omega }\frac{d\Omega }{\Omega }[2p(r,t)-2t\frac{\partial p(r,t)}{\partial t}]{|}_{t=|r-r^{\prime }|/v_s}$$

(5)

Here, $d\Omega$ represents the solid angle between the surface elements dS and $\Omega =4\pi$, corresponding to the whole surface S. Generally, the back projection method uses Equation (5) to reconstruct the light absorption distribution in a given region of interest (ROI). In fact, due to the limited number of sensor positions, Equation (5) must be discretized. If the distance between the ROI and the measuring point is significantly greater than the size of the ROI, $d\Omega /\Omega$ is approximately constant. Specifically, the discrete version of Equation (5) (all constant terms are deleted for simplicity) can be described as

$$S_s\left(r_j^{\prime }\right)=\sum _i[p(r_i,t_{ij})-t_{ij}\frac{\partial p(r_i,t_{ij})}{\partial t}]$$

(6)

where $r_i$ denotes the position of the i-th measurement point, $r_j^{\prime }$ denotes the position of the j-th ROI point, and $t_{ij}=\left|r_i-r_j^{\prime }\right|/v_s$. In both simulations and experiments, the discretization of the back projection algorithm described in Equation (6) is usually applied to two-dimensional reconstruction.

2.3. Three-Dimensional Discretization Solution of the Initial Acoustic Pressure Model

The reconstruction method derived in the previous section is based on the inversion of the discrete version of the forward model given in Equation (4). The pressure $p(r,t)$ generated by the X-ray-induced acoustic effect at a certain point r and at a certain moment t is described by Equation (4). If only the time distribution of pressure in an arbitrary unit is required, the constant term outside of the time derivative can be ignored. Then, the derivative can be used to approximately express Equation (4), as follows:

$$p(r,t)\approx \frac{I(r,t+\Delta t)-I(r,t-\Delta t)}{2\Delta t}$$

(7)

where $I(r,t)$ is the intermediate variable of the pressure of integration. It can be described as

$$I(r,t)={\int }_{S^{\prime }(r,t)}\frac{S_s\left(r^{\prime }\right)}{|r-r^{\prime }|}d{S}^{\prime }(r,t)$$

(8)

where $r^{\prime }$ denotes an absorbing point by X-ray excitation, $S_s\left(r^{\prime }\right)$ denotes the heating source at location $r^{\prime }$, $S^{\prime }(r,t)$ denotes a spherical surface that is time dependent, and $\left|r-r^{\prime }\right|$ denotes the distance between absorption point $r^{\prime }$ and propagation point r.

Considering the spherical coordinate system centered on the central position coordinate r, the surface element $d{S}^{\prime }(r,t)$ is given by

$$d{S}^{\prime }(r,t)={|r-r^{\prime }|}^{2}cos\varphi d\varphi d\theta$$

(9)

Thus, $I(r,t)$ is expressed as

$$I(r,t)={\int }_{\varphi }{\int }_{\theta }{S}_s\left(r^{\prime }\right)|r-r^{\prime }|cos\varphi d\varphi d\theta$$

(10)

Finally, we can discretize Equation (10) by calculating a set of points $S^{\prime }(r,t)$ on the position $r_l^{\prime }$ of the coordinate parameters $\varphi$ and $\theta$ of the spherical coordinate system. The equation is expressed as

$$I(r,t)\approx \Delta \varphi \Delta \theta \sum _l{S}_s\left(r^{\prime }\right)|r-r^{\prime }|cos{\varphi }_l$$

(11)

For simplicity, the constant term $\Delta \varphi \Delta \theta$ can be ignored when estimating pressure in any unit.

For each absorption point of the X-ray, we assume that the X-ray absorption is $S_s\left(r^{\prime }\right)\ne 0$ and limited to a finite three-dimensional region for each X-ray absorption point. Next, we create discrete regions of interest that are made up of regular grid points that cover the region (as shown in Figure 3). In our example, in order to facilitate the use of our experimental system for X-ray-induced acoustic reconstruction, the grid is composed of $\Delta xy$ that is a distance separated by $n_{xy}$ and y values, and $\Delta z$ that is a distance separated by $n_z$ and z. The point corresponding to the discrete integral in Equation (11) is depicted as a hollow circle in Figure 3a. By using trilinear interpolation, it is possible to represent the X-ray absorption of these grid points $S_s\left(r_l^{\prime }\right)$ as a function of the X-ray absorption of eight nearby grid points. As depicted in Figure 3b,c, eight adjacent points are represented as the red points a g in the 3D space and the four red dots in a box in 2D space, respectively. The absorption function is described as follows:

$$\begin{array}{c}\hfill S_s\left(r_l^{\prime }\right)\approx (1-\Delta x_a^{\prime })(1-\Delta y_a^{\prime })(1-\Delta z_a^{\prime })H_a+\Delta x_a^{\prime }(1-\Delta y_a^{\prime })(1-\Delta z_a^{\prime })H_b\\ \hfill +(1-\Delta x_a^{\prime })\Delta y_a^{\prime }(1-\Delta z_a^{\prime })H_c+(1-\Delta x_a^{\prime })(1-\Delta y_a^{\prime })\Delta z_a^{\prime }{H}_e\\ \hfill +\Delta x_a^{\prime }(1-\Delta y_a^{\prime })\Delta z_a^{\prime }{H}_f+(1-\Delta x_a^{\prime })\Delta y_a^{\prime }\Delta z_a^{\prime }{H}_g\\ \hfill +\Delta x_a^{\prime }\Delta y_a^{\prime }(1-\Delta z_a^{\prime })H_d+\Delta x_a^{\prime }\Delta y_a^{\prime }\Delta z_a^{\prime }{H}_h\end{array}$$

(12)

where $\Delta x_a^{\prime }=(x^{\prime }-x_a^{\prime })/\Delta xy$, $\Delta y_a^{\prime }=(y^{\prime }-y_a^{\prime })/\Delta xy$, $\Delta z_a^{\prime }=(z-z_a^{\prime })/\Delta z$ and $H_k=H(x_k^{\prime },y_k^{\prime },z_k^{\prime })$. Then, by the joint Formulas (7), (11) and (12), a linear combination of the X-ray absorption at the grid point can be used to express the sound pressure at various points and time periods. The equation is written as

$$P(r_i,t_j)=\sum _{k=1}^N{a}_k^{ij}H\left(r_k^{\prime }\right)$$

(13)

Equation (13) corresponds to the discrete forward model, which determines the pressure as the energy absorbed function in the discrete ROI, where H is expressed as $S_s$.

The position of a set of P sensors and the theoretical pressure of the instantaneous time I is a function of the optical adsorption at the grid point of Equation (12). As a result, we create a system of linear equations that may be written as a matrix as follows:

$$P=AH$$

(14)

where A is the model matrix, which is entirely dependent on the photoacoustic device measurement and the medium’s sound speed. Calculating the vector H corresponding to the X-ray absorption at the grid point is necessary for the photoacoustic reconstruction. We reduce the mean variance between the theoretical pressure in Equation (14) and the pressure that was actually measured at the same sensor position and instantaneous time, which is defined as

$$H_{sol}=\genfrac{}{}{0pt}{}{argmin}{H}{∥p_m-AH∥}^2$$

(15)

For our experimental setup in the simulated experiments, we must regularize the inversion process to obtain a representative and reasonable solution. In our example, we adopt the standard Tikhonov regularization, so Equation (15) is modified as follows:

$$H_{sol}=\genfrac{}{}{0pt}{}{argmin}{H}{∥p_m-AH∥}^2+{\lambda }^2{∥H∥}^2$$

(16)

Equation (16) can be solved by the sparse matrix algorithm that enables rapid reconstruction by take advantage of the sparsity of the matrix.

2.4. Dose Calculation

In an X-ray-induced acoustic system, when both the thermal and stress constraints are met, the acoustic emission amplitude produced by a single X-ray pulse in an X-ray-induced acoustic system is proportional to the specific amount of X-ray energy deposited in the tissues. According to Equations (3) and (4), the initial pressure $p_0\left(r\right)$ in this case is brought on by the X-ray excitation of tissues and X-ray energy deposition $A_p\left(r\right)$, as follows:

$$p_0\left(r\right)=\Gamma {\eta }_{th}{A}_p\left(r\right)$$

(17)

Here, we use the standard international (SI) unit Pa for the initial acoustic pressure, and the unit J/m³ for the X-ray energy deposition.

The generation of the initial acoustic signal pressure occurs naturally as a result of the dose deposition transformation during the process of X-ray energy absorption in tissues and organs. Grays and Pascals measure the dose and pressure using the standardized IU, respectively. The definitions of 1 Gy and 1 Pa, meanwhile, are 1 Gy = 1 J/kg and 1 Pa, respectively, and can be expressed as 1 N/m³ and 1 J/m³. The initial acoustic pressure in the target and the delivery dose have the following relationships if we assume that the X-ray pulse application time is very brief:

$$p_0\left(r\right)=\rho \Gamma {\eta }_{th}D\left(r\right)$$

(18)

where $D\left(r\right)$ kg/m³ denotes the dose deposition. SI unit kg/m³ denotes the mass density. Since the dose of instantaneous deposition and its corresponding thermal diffusion can be neglected, Equation (18) can be rewritten as

$$p_0\left(r\right)=\rho \Gamma D\left(r\right)$$

(19)

Hence, Equation (19) can be used to determine the initial acoustic pressure $p_0\left(r\right)$ produced in the tissues and organs, which is proportional to the dose delivered $D\left(r\right)$. The relative dose distribution during radiation therapy demonstrates the correlation between the initial pressure $p_0\left(r\right)$ and the local absorbed dose $D\left(r\right)$.

2.5. Reconstruction Method Based on the Tensor Dictionary Sparse Matrix

For Equation (16), we can use the least squares method for the inversion solution. However, in the case of few measured acoustic pressure data $p_m$, the reconstructed image will show many artifacts, so the sparse solution of Equation (16) is required. In this paper, we propose a dictionary learning method that uses the tensor dictionary to sparsely represent the initial acoustic pressure and then reconstruct the acoustic pressure images. We first perform the tensor dictionary sparse representation of H as

$$H^t={\mathbb{D}}^{\left(3\right)}\odot \beta$$

(20)

where ${\mathbb{D}}^{\left(3\right)}\in {\mathbb{R}}^{x\times y\times z}$ is the three-dimensional tensor dictionary, z is the depth of the third-order tensor dictionary and $\beta \in {\mathbb{R}}^{z\times n}$ is a two-dimensional sparse matrix (we assume that the size of the reconstructed image is $x\times y$, where $x=y$, and define $n=x=y$). The XACT reconstruction problem based on tensor dictionary learning can be expressed as follows:

$$\begin{array}{c}\hfill min\underset{z\in Z}{max}\left\{{∥{\beta }_z∥}_0\right\},\\ \hfill {s.t.H}^t={\mathbb{D}}^{\left(3\right)}\odot \beta ,\\ \hfill P=A{H}^t\end{array}$$

(21)

Equation (21) is the XACT reconstruction model based on the tensor format proposed in this paper. When there are very little measured acoustic pressure data P, we continue to consider the interference noise in the model as follows:

$$P=A{H}^t+\mathbb{F}$$

(22)

where $\mathbb{F}\in {\mathbb{R}}^{N_t\times N_m}$, $N_t$ denotes the number of sampling points for measuring the sound pressure signal, $N_m$ denotes the number of ultrasonic sensors, and $\mathbb{F}$ denotes the interference noise, which should obey the Poisson noise distribution. In XACT reconstruction, we should consider noise reduction and artifact suppression. According to the above sparse dictionary reconstruction model, the minimum and maximum problems of solving the 0 paradigm of the sparse matrix is a nondeterministic polynomial hard (NPH) problem, so we transform the problem of $min{max}_{z\in Z}\left\{\left|\right|{\beta }_z{\left|\right|}_0\right\}$ into the minimizing problem of the corresponding normal form $L_{p,q}$ for $\beta$, where the normal form $L_{p,q}$ of $\beta$ can be expressed as

$${∥\beta ∥}_{p,q}={\left({\sum _{j=1}^Z\left(\sum _{i=1}^n|{\beta }_{ij}{|}^p\right)}^{\frac{q}{p}}\right)}^{\frac{1}{q}}$$

(23)

Since the normal form $L_{2,1}$ is often used in sparse coding and data robustness analysis, this paper uses $L_{2,1}$ of the $\beta$ normal form as the target of the XACT reconstruction. The sparse reconstruction model can be rewritten as

$$\begin{array}{c}\hfill min{∥\beta ∥}_{2,1}\\ \hfill {s.t.H}^t={\mathbb{D}}^{\left(3\right)}\odot \beta ,\\ \hfill {∥p_m-A{H}^t∥}_{2,2}\le {\lambda }_s\end{array}$$

(24)

where ${\lambda }_s$ is the regularization parameter, which is used to realize the trade-off between the reconstruction image error and the noise of the measured data. The selection of ${\lambda }_s$ is related to the random distribution of $\mathbb{F}$. We combine the constraints and the equivalent function of sparse reconstruction to obtain the final model:

$$H^t⟵argmin{∥p_m-A{\mathbb{D}}^{\left(3\right)}\odot \beta ∥}_{2,2}+\frac{{\lambda }_s}{2}{\left|\right|\beta \left|\right|}_{2,1}$$

(25)

For the solution calculation of Equation (25), we define the sparse operator ⊙ as follows:

$$H_{xy}^t=\sum _{z=1}^Z{\mathbb{D}}_{xyz}^{\left(3\right)}{\beta }_{zy}$$

(26)

The detailed training process of the tensor sparse dictionary model for Equation (25) is shown in Algorithm 1. The training data are composed of different densities of XACT sensor signals collected, described as:

$$X=P_m\left\{\begin{array}{c}{P_m}^{high}\hfill \\ {P_m}^{low}\hfill \end{array}\right.$$

where ${P_m}^{high}$ denotes the high density signal, and ${P_m}^{low}$ denotes the low density signal.

Algorithm 1: The training process of the tensor sparse dictionary model for Equation (25).

Input: training data X, parameters A, ${\lambda }_s$, the maximum number of iterations $Miter$, $num$.  Output: the dictionary $\mathbb{D}$, the sparse coefficient $\beta$.  1 Initialize the dictionary $\mathbb{D}$, sparse coefficient ${\beta }_0=0\in R^{r\times n\times k}$, parameters A, ${\lambda }_s$, Lagrange multiplier $\xi$, and when $\eta >1$, make $B_1={\beta }_0,t_1=1$.  2 for$iter=1$ to $Miter$do$\phantom{(}$  3  for$p=1$ to $num$do$\phantom{(}$  4     $L_p={\eta }^p({\sum }_l^k{∥{\mathbb{D}}^{\Gamma }\mathbb{D}∥}_{2,2})$  5     $\nabla f\left(B_p\right)=2({\mathbb{D}}^T\ast \mathbb{D}\ast B_p-{B_p}^T\ast X)$  6     ${\beta }_p=Pro{x}_{\xi /L_P}(B_P-\nabla f\left(B_P\right)/L_p)$  7     $t_{p+1}=\frac{(1+\sqrt{1+4t_p^2})}{2}$  8     $B_{p+1}={\beta }_p+\left(\frac{t_{p-1}}{t_{p+1}}\right)({\beta }_p-{\beta }_{p-1})$  9  end  10   $\phantom{(}$$\widehat{\mathbb{D}}=fft(\mathbb{D},\left[\right],3),\widehat{\beta }=fft(\beta ,\left[\right],3)$  11   for$l=1$ to K do  12     $\phantom{(}$ Solve $\lambda$ by Conjugate Gradient Method (where $\lambda$ is the parameter to be solved in Lagrange dual algorithm)  13     ${\widehat{\mathbb{D}}}^{\left(l\right)}=({\widehat{X}}^{\left(l\right)}{\widehat{\beta }}^{{\left(t\right)}^H}){({\widehat{X}}^{\left(l\right)}{\widehat{X}}^{{\left(t\right)}^H}+diag\left(\lambda \right))}^{-1},l\in \left\{1,2,\cdots ,k\right\}$  14   end  15   $\phantom{(}$$\mathbb{D}=ifft(\widehat{\mathbb{D}},\left[\right],3)$  16 end$\phantom{(}$

In Algorithm 1, the $3\sim 9$ steps are the training process of solving the sparse coefficient $\beta$ with the specified tensor dictionary $\mathbb{D}$. In step 4, the superscript symbol $\Gamma$ is the conjugate transpose of the tensor, k is the number of training datasets, and $\eta$ is the variation coefficient of the iterative solution model. In step 5, $\nabla f$ is the gradient defined in the tensor space. In step 6, $Prox$ is soft threshold operator. We use the soft threshold algorithm to solve ${\beta }_p$. The $10\sim 15$ steps are the training process of solving tensor dictionary $\mathbb{D}$ with the specified sparse coefficient $\beta$, where $diag$ denotes returning a diagonal matrix.

Finally, based on the previous derivation of physical mechanism, we can combine Equation (25) with dose Equation (19) to establish the dose equation of the entire X-ray-induced acoustic system. It is expressed as

$$p_0\left(r\right)=\rho \Gamma D\left(r\right)=A{H}^{t=0}$$

(27)

By Equation (27), we can obtain the dose image $D\left(r\right)$.

3. Experiment

3.1. X-ray Acoustic Signal Datasets

Our experimental datasets are based on 30 simulated X-ray acoustic signal datasets of gross tumor volume (GTV) target areas of nasopharyngeal carcinoma. The step of dataset generation includes the collection of clinic nasopharyngeal carcinoma datasets and XA signal simulation.

3.1.1. Clinic Nasopharyngeal Carcinoma Datasets

After ethical institutional review board approval (SCCHE-02-2022-003), 30 nasopharyngeal carcinoma datasets from Sichuan Cancer Hospital were used to XA signal simulation. The patients were treated with volume modulated arc therapy. The radiotherapy planning system is Monaco (Elekta medical system). The patients CT scans are shown in Figure 4, and the size of every slice of data is $512\times 512$.

According to CT scans of the tumor target area contoured and the prescription dose formulated by the doctor, the radiation therapy physicists design treatment plans. Figure 5 shows contoured GTV target areas and the dose distribution from the TPS information of a nasopharyngeal carcinoma case.

3.1.2. Acoustic Signal Simulation

Based on the dose image reconstruction system (Figure 2a), we simulate X-ray acoustic signals of nasopharyngeal carcinoma GTV. Figure 6 shows the simulation workflow. At first, we use a kind of Monte Carlo software called Primo to simulate the medical LINAC. According to the patients’ CT scans and treatment planning, the initial dose distribution is simulated. Then, based on the X-ray acoustic effect theory, the initial pressure signal can be calculated, and the sound pressure time function can be obtained by calculating the LINAC pulse time function. Finally, we use the K-wave toolbox in MATLAB (MathWorks Inc., Natick, MA, USA), and obtain the simulated acoustic signals. During radiotherapy, the tumor target area is the most important therapeutic part, so the GTV target areas are set as the acoustic source and only simulated. In the paper, based on the nasopharyngeal carcinoma datasets provided by the Hospital, we build 30 simulated acoustic signals datasets of GTV.

3.2. Evaluation Metrics

To assess the reconstruction performance of the proposed method, we use mean-square error ($MSE$), peak signal-to-noise ratio ($PNSR$) and structural similarity ($SSIM$) as evaluation criteria.

$MSE$ is the average squared difference between the real image and the noisy image. It can be written as

$$MSE=\frac{1}{mn}\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}{[I(i,j)-K(i,j)]}^2$$

(28)

where I is the real image of size $m\times n$, and K is a noisy image.

$PSNR$ is the ratio of the peak signal to average noise. In other words, it is the ratio of the peak signal to MSE and is defined as

$$PSNR=10\times {log}_{10}\left(\frac{MA{X}_I^2}{MSE}\right)$$

(29)

where $MA{X}_I^2$ is the maximum possible pixel value in the picture.

$SSIM$ is used to evaluate the similarity of image X and image Y, expressed as

$$SSIM=\frac{(2{\mu }_x{\mu }_y+C_1)(2{\sigma }_{xy}+C_2)}{({\mu }_x^2+{\mu }_y^2+C_1)({\sigma }_x^2+{\sigma }_y^2+C_2)}$$

(30)

where ${\mu }_x$ and ${\mu }_y$ denote the mean values of X and Y, respectively. ${\sigma }_x^2$ and ${\sigma }_y^2$ denote the variance of X and Y, respectively. ${\sigma }_{xy}$ denotes the covariance of X and Y.

3.3. Result

3.3.1. The Dose Image Reconstruction

In the experiment, all methods use the simulated X-ray acoustic signals datasets of GTV to reconstruct the dose image. Figure 7 shows the qualitative results of the 3D dose distribution of GTV reconstructed by our method. We can see that the maximum values of X direction and Y direction are different in every case because the pixel spacing of every patient case is different. However, the pixel interval of the z direction of every case is the same, 3 mm. Due to the different number of slices in different patient, the maximum value of the z direction is also different.

Table 1. Quantitative results of 3D dose images reconstruction in Figure 7.

Case	MSE	PSNR	SSIM
Case 1	0.0005	35.7850	0.8041
Case 2	0.0009	34.4623	0.7602
Case 3	0.0011	30.9631	0.7035
Case 4	0.0005	36.5276	0.8168

shows the quantitative results of our method, which correspond to the four results shown in Figure 7. As can be seen in Table 1, the higher PSNR and SSIM values and lower MSE value indicate a better image quality.

3.3.2. Comparison with the Exiting Methods

Our experiment introduce two model architectures that are the most representative methods, namely FBP and iteration time inversion (ITR), as the baseline. FBP is often applied to XACT imaging modality and 2D /3D dose image reconstruction [30]. ITR is used to reconstruct wave-based X-ray acoustic 3D imaging, which can effectively reduce the artifacts caused by incomplete data, noise and limited sampling [31]. We randomly select four cases for visual comparisons of the reconstruction performance in the two models and our model, as shown Figure 8. For better visual effect, we crop the images and preserve the feature areas. The first column is the initial dose distribution examples, and these are original images. The second, third, and fourth columns are the dose images reconstructed by FBP, ITR and our method, respectively. As can be seen in Figure 8, artifacts obviously appear around GTV when using FBP and ITR, but there is nothing in our method. In addition, the morphological and textural details of the reconstructed images using FBP and ITR show significant differences from the initial dose images, especially at the edges of the contour. However, the results of our method are similar to the initial dose images in terms of visual effect.

Figure 9, Figure 10 and Figure 11 represent the statistical results of the measurements of patients datasets. In Figure 9, we can see that the $MSE$ values of our method are lowest, so our method is better than the other two methods in terms of the performance of the $MSE$ indicator. As can be seen in Figure 10 and Figure 11, the $PSNR$ and $SSIM$ values of our method are higher than other methods, and it demonstrates that our method displays the best performance for both evaluation indicators.

To characterize the spatial accuracy of the imaging technique, we compared the initial doses with the reconstructed dose by FBP, ITR and our method on the vertical and horizontal profiles. As shown in Figure 12a–d, along the vertical and horizontal white dashed lines, we can obtain the X profiles (shown in Figure 12e–g) and Y profiles (shown in Figure 12h–j), respectively. Due to orders of magnitude differences between the reconstructed dose of the three methods, the values are normalized. From Figure 12e–j, we can see that the reconstructed dose distribution from our method is closer to the initial dose distribution, with less error and more accuracy than FBP and ITR.

4. Discussion and Conclusions

Recent research studies have shown the feasibility of using XACT imaging as a dosimetry technique during radiotherapy [30]. Dose reconstruction methods primarily include FBP and ITR. FBP has gained popularity in XACT due to its simplicity and computational efficiency. However, the X-ray-induced acoustic signal is a time-series signal with both positive and negative intensities, and the radial projection sum is unable to quantify the dose. Furthermore, the X-ray excited acoustic signal is mixed with a lot of redundant noise, and FBP is prone to artifacts in the reconstructed images. ITR is a wave equation-based method, which can perform the spherical characteristic of the sound signal propagation and reduce the artifacts. However, when there are few detectors measuring acoustic signals or the geometric structure of the detectors are incomplete, the obtained sound pressure data are decreased, and the images still lie in some artifacts. Given the above problems, we comprehensively make use of the theory that the amplitude of acoustic signal excited by X-ray is proportional to the deposition of X-ray energy and the physical characteristics of acoustic wave propagation. Meanwhile, we combine tensor sparse coding with compressed sensing data reconstruction. By means of the tensor product, the two-dimensional dictionary learning is extended to three-dimensional to utilize the spatial information of XACT signals more effectively. Then, to solve the noise problem of the reconstructed image caused by spare sampling, we design an alternate iteration solution between the tensor sparse coefficient and the tensor dictionary. The iterative shrinkage threshold algorithm is used to solve the tensor sparse coefficient, and the Lagrangian dual method is used to solve the tensor dictionary. This can effectively perform XACT signal data interpolation reconstruction and then suppress the artifacts generated.

In the experiment, XACT simulation was carried out based on the clinic datasets of nasopharyngeal carcinoma. The biological anatomy of nasopharyngeal carcinoma has the characteristics of complex structure, a great variability in the shape and size, and a high similar intensity value to adjacent tissues, which is a great challenge to the dose reconstruction. However, this is closer to the actual situation, which has guiding significance for further clinical trials of XACT technology.

Compared to some state-of-the-art methods, the experimental results show that the proposed method can significantly improve the quality of reconstructed images and the accuracy of dose distribution. It is helpful for X-ray excited acoustic imaging to become an in vivo dose measurement technology. In this paper, we build an XACT simulation system, but there is still a error between the experimental results and the treatment planning dose because the model parameter setting is ideal and simple. If the physical density of the human body is accurately measured and the loss of acquired X-ray acoustic signals is reduced, the accuracy of dose distribution will be greatly improved. In addition, the experiments only use simulated data and do not use experimental data to verify. The results show some limitation. At present, experimental XACT imaging research works are under study [23]. Therefore, in the future, we intend to use X-ray acoustic experimental data to improve the model and perform the in vivo dose monitoring in clinical radiotherapy with X-ray excitation acoustic imaging technology.