About one method of parallelization of calculations during the reconstruction of a tomographic image

TL;DR

This paper introduces a novel parallelization method for tomographic image reconstruction that leverages matrix symmetry and voxel numbering to split the problem into smaller, independent systems, enabling faster computations.

Contribution

The proposed method uniquely exploits matrix symmetry and voxel numbering to divide the system into smaller, independent systems for parallel processing in image reconstruction.

Findings

Matrices of the divided systems are four times smaller

The method reduces matrix dispersion

Enables parallel computation of tomographic reconstructions

Abstract

A new method for solving systems of linear algebraic equations of a special type arising in solving problems of image reconstruction has been proposed. This method, due to a certain symmetry of the matrix and the choice of the voxel numeration method for two-dimensional problems, allows us to divide the initial system of algebraic equations into two independent systems, which enables us to carry out the calculation in parallel. The dimension of the matrices of the resulting systems is 4 times less than the dimension of the original matrix, and these matrices are less dispersed.

Tables

Table 1. Table 1: CPU time and errors of the calculations

№	size of $A$	CPU time	error	CPU time for	error for
				our method	our method
1	$1024\times 2280$	1.324 sec	6.203e-8	0.227 sec	2.623e-9
2	$4096\times 9120$	73.434 sec	1.304e-6	3.233 sec	1.577e-7
3	$16384\times 36480$	4948.770 sec	4.498e-5	76.473 sec	2.361e-6

Equations

$$Af=p$$

$$w_{i,j}=w_{M-i+1,N-j+1}.$$

$$A^{(1)}{f}^{(1)}=p^{(1)},$$

$$A^{(2)}{f}^{(2)}=p^{(2)},$$

$$w_{i,j}^{(1)}=w_{i,j}-w_{M-i+1,j},w_{i,j}^{(2)}=w_{i,j}+w_{M-i+1,j},i=1,\dots ,\frac{M}{2},j=1,\dots ,\frac{N}{2},$$

$$p_i^{(1)}=p_i-p_{M-i+1},p_i^{(2)}=p_i+p_{M-i+1},i=1,\dots ,\frac{M}{2}.$$

$$f_j=\frac{f_j^{(2)}+f_j^{(1)}}{2},f_{N-j+1}=\frac{f_j^{(2)}-f_j^{(1)}}{2},j=1,\dots ,\frac{N}{2},$$

$$f_j^{(1)}=f_j-f_{N-j+1},f_j^{(2)}=f_j+f_{N-j+1},j=1,\dots ,N.$$

$$f_j^{(1)}=-f_{N-j+1}^{(1)},f_j^{(2)}=f_{N-j+1}^{(2)},j=1,\dots ,N,$$

$$A\left(\frac{f^{(1)}+f^{(2)}}{2}\right)=\frac{1}{2}A{f}^{(1)}+\frac{1}{2}A{f}^{(2)}=\frac{1}{2}{p}^{(1)}+\frac{1}{2}{p}^{(2)},$$

$$A{f}^{(1)}=p^{(1)},$$

$$A{f}^{(2)}=p^{(2)}.$$

$$\{\begin{array}{l}{A}^{(1)}{f}^{(1)}=p^{(1)},\\ -A^{(1)}{f}^{(1)}=-p^{(1)},\end{array}\iff A^{(1)}{f}^{(1)}=p^{(1)}.$$

$$\{\begin{array}{l}{A}^{(2)}{f}^{(2)}=p^{(2)},\\ A^{(2)}{f}^{(2)}=p^{(2)},\end{array}\iff A^{(2)}{f}^{(2)}=p^{(2)}.$$

$${\bar{f}}_j=\frac{{\bar{f}}_j^{(2)}+{\bar{f}}_j^{(1)}}{2},{\bar{f}}_{N-j+1}=\frac{{\bar{f}}_j^{(2)}-{\bar{f}}_j^{(1)}}{2},j=1,\dots ,\frac{N}{2},$$

$$\Vert f{\Vert }_2^2=\sum _{j=1}^N{f}_j^2=\sum _{j=1}^{N/2}\left({\left(\frac{f_j^{(2)}+f_j^{(1)}}{2}\right)}^2+{\left(\frac{f_j^{(2)}-f_j^{(1)}}{2}\right)}^2\right)=$$

$$\frac{1}{2}\sum _{j=1}^{N/2}(f_j^{(1)}{)}^2+\frac{1}{2}\sum _{j=1}^{N/2}(f_j^{(2)}{)}^2=\frac{1}{2}\Vert f^{(1)}{\Vert }_2^2+\frac{1}{2}\Vert f^{(2)}{\Vert }_2^2.$$

$$\left(\begin{array}[]{cccccc}1&3&5&7&9&1\\ 2&4&6&8&3&7\\ 7&3&8&6&4&2\\ 1&9&7&5&3&1\end{array}\right)\left(\begin{array}[]{c}f_{1}\\ f_{2}\\ f_{3}\\ f_{4}\\ f_{5}\\ f_{6}\end{array}\right)=\left(\begin{array}[]{c}5\\ 6\\ 8\\ 7\end{array}\right).$$

$$\left(\begin{array}[]{ccc}0&-6&-2\\ -5&1&-2\end{array}\right)\left(\begin{array}[]{c}f_{1}^{(1)}\\ f_{2}^{(1)}\\ f_{3}^{(1)}\end{array}\right)=\left(\begin{array}[]{c}-2\\ -2\end{array}\right),$$

$$\left(\begin{array}[]{ccc}2&12&12\\ 9&7&14\end{array}\right)\left(\begin{array}[]{c}f_{1}^{(2)}\\ f_{2}^{(2)}\\ f_{3}^{(2)}\end{array}\right)=\left(\begin{array}[]{c}12\\ 14\end{array}\right).$$

$$f_1^{(1)}=14x,f_2^{(1)}=10x,f_3^{(1)}=-30x+1,\text{where}x\in R,$$

$$f_1^{(2)}=84y,f_2^{(2)}=80y,f_3^{(2)}=-94y+1,\text{where}y\in R.$$

$$f=\left(\begin{array}[]{c}7x+42y\\ 5x+40y\\ -15x-47y+1\\ 15x-47y\\ -5x+40y\\ -7x+42y\end{array}\right).$$

$$f^{(1)}={\left(0.3512,0.2508,0.2475\right)}^T,\Vert f^{(1)}{\Vert }_2=0.4975,$$

$$f^{(2)}={\left(0.3542,0.3373,0.6036\right)}^T,\Vert f^{(2)}{\Vert }_2=0.7769,$$

$$f={\left(0.3527,0.2941,0.4256,0.1781,0.0433,0.0015\right)}^T,\Vert f{\Vert }_2=\mathrm{0.6523.}$$

$$\det (A)=\det (A^{(1)})\cdot \det (A^{(2)}).$$

$$B^{(1)}=(A^{(1)}{)}^T\cdot A^{(1)},B^{(2)}=(A^{(2)}{)}^T\cdot A^{(2)}.$$

Full text

Keywords. system matrix, algorithmic computational speedup, cluster computing, tomosynthesis, image reconstruction.

1 Introduction

An important method of non-invasive imaging widely used in medical diagnostics is x-ray computed tomography. Image reconstruction is performed both by direct and iterative methods. Iterative methods are preferred when using a limited set of projections, for example, in digital tomosynthesis. Digital tomosynthesis is a new method of three-dimensional reconstruction of an object using projections at limited viewing angles. Currently, the method of simultaneous algebraic reconstruction (SART) is a generally accepted iterative reconstruction method [1]. The disadvantage of SART is its high computational complexity and cost. In the works [2, 3, 4], the possibility of improving the SART image reconstruction speed was investigated. The use of parallel methods to speed up the SART reconstruction was proposed in [5, 6, 7, 8, 9, 10, 11].

In [12] the method of partitioning the reconstruction region into an uneven voxel grid was considered. This allows one to split the original three-dimensional problem into a chain of two-dimensional independent tasks.

In this paper, we propose a method that, by virtue of a certain symmetry of the matrix and the choice of the numbering method for voxels for two-dimensional problems, allows us to split the initial system of algebraic equations into two independent systems, which makes it possible to perform a parallel calculation. The size of the matrices of the resulting systems is 4 times less than the size of the original matrix, and these matrices are less dispersed.

2 Methods

Let $f=(f_{1},\ldots,f_{N})\in R^{N}$ be a discrete representation of an unknown image ($N$ is the total number of image voxels) to be reconstructed, $p=(p_{1},\ldots,p_{M})\in R^{M}$ be the measured projection data ($M$ is the total number of rays), and $A=(w_{i,j})$ be a known system matrix, whose element $w_{i,j}$ is the weighting factor that represents the contribution of the $j$th voxel of the $i$th ray integral. Hereinafter we shall assume that the numbers $N$ and $M$ are even (as a rule, it takes place in practice). The image reconstruction problem can be formulated as a system of linear equations:

[TABLE]

Suppose that during the motion of the radiator, the irradiation is performed in positions that are pairwise symmetrical about the axis OZ. let the number of radiation is $K$. The reconstruction region is divided into layer-by-layer voxels, and each of the $j$ voxels is characterized by the parameter $f_{j}$. However, voxels are numbered in a special order (according to the vertical "snake" rule) as shown in Fig. 1 for the case $N=32$.

In this case, the matrix $А$ is symmetric in the sense

[TABLE]

Consider two systems of linear algebraic equations that are independent of each other:

[TABLE]

[TABLE]

where the elements of matrices $A^{(1)}$ and $A^{(2)}$ depend on the elements of the matrix A, and are found by the formulas:

[TABLE]

and the elements of vectors $p^{(1)}$ and $p^{(2)}$ depend on the elements of vector $p$:

[TABLE]

Theorem 1. Let the matrix of system (1) satisfy the symmetry conditions (2), whereas $f^{(1)}$ and $f^{(2)}$ are the solutions of systems (3) and (4) respectively, then $f=(f_{1},\ldots,f_{N})$, where

[TABLE]

is a solution to equation (1).

Proof. We introduce the notation

[TABLE]

It’s obvious that

[TABLE]

From equation (1), by virtue of the equalities $f=\frac{f^{(1)}+f^{(2)}}{2}$ and $p=\frac{p^{(1)}+p^{(2)}}{2}$, we get

[TABLE]

where $p^{(1)}$ and $p^{(2)}$ are extended for all $i=1,\ldots,{N}$ so that the symmetry conditions are satisfied (6).

By virtue of the symmetry condition (2) for matrix $A$, vector $Af^{(1)}$ satisfies the first symmetry condition from (6), whereas vector $Af^{(2)}$ satisfies the second symmetry condition from (6). Since the vector representation as the sum of two vectors satisfying the symmetry conditions (6) is unique, then (7) splits into two systems:

[TABLE]

[TABLE]

System (8) can be represented as

[TABLE]

And system (9) can be represented as

[TABLE]

Theorem 1 is proved.

The size of the matrices $A^{(1)}$ and $A^{(2)}$ is 4 times less than the size of the original matrix $A$. In addition, the matrices $A^{(1)}$ and $A^{(2)}$ are less dispersed. The independence of these matrices makes it possible to find solutions of systems (3) and (4) by the method of parallel calculations. Obviously, if the number of arithmetic operations required for solving system (1) is proportional to $N^{\alpha}\cdot M^{\beta}$ ($\alpha>0,\,\beta>0$), then the time spent on solving system (1) using the same method with the help of (5) will be $2^{\alpha+\beta}$ times less.

Theorem 2. Let $\bar{f}^{(1)}=(\bar{f}^{(1)}_{1},\ldots,\bar{f}^{(1)}_{\frac{N}{2}})$ - be the normal pseudo-solution of system (3), and $\bar{f}^{(2)}=(\bar{f}^{(2)}_{1},\ldots,\bar{f}^{(2)}_{\frac{N}{2}})$ be the normal pseudo-solution of system (4), then

[TABLE]

is a normal pseudo-solution of equation (1).

The proof of Theorem 2 follows from the following obvious equalities:

[TABLE]

[TABLE]

Example 1. Consider the case $N=6$, $M=4$

[TABLE]

In order to solve system (13) consider the following two systems:

[TABLE]

[TABLE]

The general solution of systems (14) and (15) respectively has the form:

[TABLE]

[TABLE]

From (16) and (17) based on Theorem 1, we obtain the general solution of system (13)

[TABLE]

Normal pseudo-solutions, up to four decimal places, of systems (14) and (15) respectively, have the form:

[TABLE]

[TABLE]

Normal pseudo solutions, up to four decimal places, for the initial system (13) have the form:

[TABLE]

For this example, the validity of Theorem 2 is confirmed.

Numerous calculations for square matrices show the validity of the following equality

[TABLE]

If the matrix $A$ is not square and satisfies symmetry conditions in the sense of (2), then the matrix $B=A^{T}\cdot A$ will also be symmetric in the sense of (2). In addition, the following equalities hold true:

[TABLE]

Example 2. The emitter rotates about a fixed axis parallel to the detector matrix. An object is fixed between the emitter and detector matrix. The lower boundary of the object is at the distance of $h_{e}=1$ m from the emitter and at a distance of $h_{m}=0.25$ m from the detector matrix. The size of the matrix is $L_{m}\times L_{m}$ ($L_{m}=0.43$ m) and the resolution is $N_{p}\times N_{p}$ ($N_{p}=1024$ pixels). The spot size of an x-ray tube is $a_{e}=7\cdot 10^{-4}$ m. When rotating the emitter the angle of inclination changes from $-\gamma$ to $+\gamma$, where $\gamma=30^{\circ}$ (see Fig. 2).

The passage of the beam through voxels of the object is computationally modeled, thus building the system of equations using system of matrix $A$ and the right-hand side $p$ of the system (1). In this example, the reconstruction of the Shepp-Logan phantom is modeled once the object is non-uniformly divided into a voxel grid, as shown in [12]. This allows to obtain a series of two-dimensional independent problems. In example 2 reconstruction of one of the $2D$ - problem using the method proposed is considered. In the case when $N=32\times 32$ for the $2D$ - object restoration by the direct method 12 exposures from points located evenly between the extreme left position and center are sufficient as well as at the points symmetric to them relative to the central axis (24 positions in total). Thus we get matrix A with the size $1024\times 2280$. With each doubling of the number of partitions of the object in each direction, the number of positions doubles, and the size of the matrix $A$ increases 16 times. For example, if $N=64\times 64$ then the exposure is produced from 48 positions and the size of matrix $A$ is $4096\times 9120$. While in the case when $N=128\times 128$, the exposure is produced from 96 positions with matrix $A$ size $16384\times 36480$. The calculations obtained for these systems of equations are given in Table 1. There are three dimensions for matrix $A$ system, the calculation time of the direct method, the errors of the direct method, the calculation time of the method proposed in this paper, and the errors of the proposed method. Table 1 shows that the advantages of the proposed method are in the reduced calculation time by an order of magnitude and also in the calculation accuracy. Fig. 3 shows the reconstruction results of the Schepp-Logan phantom for dividing an object into a) $32\times 32$, b) $64\times 64$ and c) $128\times 128$ voxels.

The proposed method is applicable for an arbitrary symmetric motion geometry of the source and detector, which ensures the symmetry of the reconstruction matrix. For example, in paper [13] it was considered the problem of recovering images of nuclear waste in a cylindrical container using tomographic Gamma Scanning. Each reconstructed layer was divided into $72$ voxels using the polar coordinate system, and the voxel numbering was chosen in the order shown in Fig. 4a. If we redesign the voxel numbering by our proposed method, as shown in Fig. 4b, then the system of linear algebraic equations from [13] can be divided into two independent systems with $36$ unknowns and the number of equations two times smaller than the original one.

3 Conclusion

A new method for solving systems of linear algebraic equations of a special type arising in solving problems of image reconstruction has been proposed. This method, due to a certain symmetry of the matrix and the choice of voxel numbering method for two-dimensional problems, allows us to divide the initial system of algebraic equations into two independent systems, which makes it possible to carry out the calculation in parallel. The size of the matrices of the resulting systems is 4 times less than the size of the original matrix, and these matrices are less dispersed.

The proposed method is applicable for an arbitrary symmetric motion geometry of the source and detector, which ensures the symmetry of the reconstruction matrix.

The combined use of this method and efficient algorithms, such as those proposed in [6, 11, 12, 13], for reconstructing a tomographic image will increase the speed of these algorithms at least 4 times, while maintaining the accuracy of the reconstruction.

This parallelization method can also be applied when solving multidimensional problems of mathematical physics considered in symmetric domains, the approximation of which leads to matrices symmetric in the sense of (2).

4 Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

5 Acknowledgements

The research results outlined in this paper were obtained with financial support from the Ministry of Education and Science of the Russian Federation, as part of the execution of the project entitled "Establishment and creation of hightech production for manufacturing of digital X-ray complex with the function of tomographic image synthesis", in accordance with the Government resolution of the Russian Federation №218 from 09.04.2010.