An automatic algorithm to exploit the symmetries of the system response matrix in positron emission tomography iterative reconstruction

Maximum likelihood estimation maximization (MLEM) (Vardi et al 1985) is one of the most popular algorithms used in iterative PET reconstruction. In its simplest form, the MLEM can be written for a system with N lines of response (LOR) and a field of view (FOV) of $V$ voxels as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle X^{q+1}=\frac{X^q}{M^\intercal\cdot 1} M^\intercal \cdot\left(\frac{P}{M \cdot X^q}\right) \label{eq:MLEM} \nonumber \end{align} \tag{ 1 }$$

where $X^q \in \!R^V$ is a vector containing the image at the qth iteration, $M \in \!R^{N\times V}$ is the SRM, $P \in \!R^N$ is a vector containing the coincidences detected in each LOR. The element $a, b$ of the SRM contains the probability that a photon pair emitted from the bth voxel of the FOV is detected in the ath LOR of the system (Qi and Leahy 2006). The set of values associated with the LOR a, $M_{a, 1}, M_{a, 2}, ..., M_{a, V}$, is often referred to as the tube of response (TOR) T_a. Since the number of nonnull element in each TOR is limited, the size of the SRM can be reduced exploiting its sparseness, i.e. storing only TOR's non-null elements. In order to exploit symmetries, it is more convenient to represent each voxel of the FOV using a coordinate triplet $(b_x, b_y, b_z)$ rather than the linear index b. A single entry Tⁱ of a TOR can be stored as a probability $T^i_0$ value followed by three voxel coordinates $(T^i_x, T^i_y, T^i_z)$ (see figure 1). Using this representation, the theoretical SRM memory footprint is  ≈$N_0\cdot(3\cdot S_{{index}}+S_{{value}})$ bytes where N₀ is the number of non-null elements in the whole matrix, S_index is the size of the variable representing each index of the triplet and $S_{{value}}$ is the size of the variable representing the probability. In this work, S_index and $S_{{value}}$ are set to 1 byte (unsigned char) and four bytes (float) respectively. Due to this choice, the size of the FOV is limited to 255 voxels per direction. The generalization of the SSA to other cases is trivial.

Zoom In Zoom Out Reset image size

Intuitively, two TORs can be considered symmetric if they are equal up to a geometric transformation. This means that their probability values are exactly equal and they can be mapped one into the other with a geometric transformation. However, in practice, due to small differences in their probabilities, the number of TORs that are exactly symmetric, according to this definition, is very limited. Therefore, a less strict condition has to be implemented in order to improve the compression efficacy of the SSA. In this work, two TORs l and m are said to be symmetric within a threshold t if:

• (i)  
  It is possible to transform the coordinates of all the elements of l into those of the m using a voxel space transformation.
• (ii)  
  All the probability values of the so-obtained TOR lie within a relative threshold t with respect to the original.

When a threshold t  =  0 is applied, the strict def for the IRIS SRMinition of symmetry is obtained. Only combinations of reflections, translations, and coordinates swap in the voxels space are exploited in this work. Using these hypotheses, the existence of a symmetry between two TORs can be reduced to three sufficient conditions:

• A  
  The number of non-null TOR elements n₀ of l and m is the same.
• B  
  Sorting the elements of the TORs l and m according to their probabilities, the following relation holds for each element $i=1 \ldots n_0$

  $$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle d(l^i_0,m^i_0)=\left|{l^i_0-m^i_0}\right| / \min{(l^i_0,m^i_0)} \leqslant t. \nonumber \end{align*}$$
• C  
  It is possible to sort l and m and to find $\vec {A}=(A_x, A_y, A_z)$ and $\vec {k}=(k_x, k_y, k_z)$ and z such that the following relation is true for each TOR entry i

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:symmetries} \vec {l^i}+\vec {A}\odot S_z(\vec {m^i})=\vec {k} \nonumber \end{align} \tag{ 2 }$$

with $A_x, A_y, A_z = \pm 1$, $\vec {k}$ a constant vector not depending on i, the $\odot$ is the element-wise product between vectors and S_z is the 'swap operator'. S_z transforms a vector into another containing its components permuted. Condition C is a convenient formula to characterize reflections, translations, and coordinate swap. For example, a pair of TORs symmetric up to a reflection along x axis can be described by a transformation with $A=(1, -1, -1)$, $\vec {k}=(0, 0, 0)$ and no coordinate swap (S_z equal to identity). A 90 degree rotation with respect to the z-axis can be implemented swapping $x, y$, i.e. using S_z such as $(x, y, z) \rightarrow (y, x, z)$, and setting $A=(-1, -1, -1)$ and $\vec {k}=(0, 0, 0)$. Since each component of $\vec {A}$ can assume only two values, i.e. $= \pm 1$ and the number of permutations of 3 indices is 6, 48 different symmetry classes are possible.

A naïve implementation of the SSA would require checking all the possible TOR pairs in the SRM against the conditions A–C. However, this would lead to a computational complexity of  ∼N², with N equal to the number of TORs contained in the SRM. To reduce the computational burden of the SSA, the TORs are first grouped according to the number of non-null entries, i.e. voxels that had a non-null probability during the model computation, and then each group is stored in a separate file on disk. This reduce the memory demand of the SSA since each TOR group can be analyzed independently and allows for a straightforward multi-core implementation, where each thread searches the symmetries within a group. Using this method, condition A is automatically satisfied within the same group (see algorithm 1) and the complexity of the SSA is reduced to  ∼$G \cdot J^2$, where G is the number of groups and J is the average number of TOR per group, with both G and J $\ll N$. To check a pair of TORs against condition B, both the TORs are first sorted according to their probability values. This way, condition B can be checked element-wise (see algorithm 2). Finally, any TOR pair that meets condition B is validated against condition C. An algorithm to reduce the complexity of condition C is implemented. This solution allows to validate the condition pixel-wise, therefore reducing its complexity to the number of non-null elements of the involved TORs. The problem of checking condition C is that the values of $\vec {k}$ are not known a priori. A necessary condition to find $\vec {k}$ can be extracted by summing over the non-null entries on both sides of equation (2). This yields to equation (3), which leads to equation (4).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:KappaSum} (\sum\limits_{i=1...n_0}{\vec {l^i}})+\vec {A}\odot S_z (\sum\limits_{i=j...n_0}{\vec {m^{\,j}}})=n_0\cdot \vec {k} \nonumber \end{align} \tag{ 3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:KappaValues} \vec {k} =\frac{(\sum\limits_{i=1...n_0}{\vec {l^i}})+\vec {A}\odot S_z(\sum\limits_{j=1...n_0}{\vec {m^{\,j}}})}{n_0}. \nonumber \end{align} \tag{ 4 }$$

Algorithm 1. Implementation of condition A.

Given the TORs l, m each with $n_l, n_m$ non-null entries

1: if $n_l \neq n_m$ then

2:  No symmetry between l and m

3: end if

4: Check condition B

Algorithm 2. Implementation of condition B.

Given the TORs l, m fulfilling condition A

1: Sort all entries of l, m according to their probability values

2: for i in 1,n0 do

3:   if $d(l^i_0, m^i_0)>t$ then

4:    No symmetry between l and m

5:   end if

6: end for

7: Check condition C

Equation (4) gives a necessary but not sufficient condition for $(z, \vec {A}, \vec {k})$ to be a symmetry between two TORs l and m. For a given pair of TORs, 48 values of $\vec {k}$ are found, one for each combination of $\vec {A}$ and z. Some of the values found can be excluded by discarding those transformations that do not correspond to a $\vec {k}$ vector made of integer components. For all the other values, a direct test is needed. The direct test is performed transforming m according to $(z, \vec {A}, \vec {k})$ and then sorting l and m according to their coordinates. Since no tolerance is allowed in the coordinates of two symmetrical TORs, it is possible to check condition C pixel-wise as described in algorithm 3). If $(z, \vec {A}, \vec {k})$ transforms TOR l into TOR m within the threshold, then a symmetry is found. The result of the compression is a new SRM containing only the fundamental TORs (FT) and the transformations $(z, \vec {A}, \vec {k})$ to obtain those that are not fundamental. A non-fundamental TOR m can be retrieved by knowing the $(z, \vec {A}, \vec {k})$ values and the fundamental TOR l, by using the formula

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {m^i}=S^{-1}_z (\vec {A} \odot(\vec {k}- \vec {l^i})). \label{eq:Transformation} \nonumber \end{align} \tag{ 5 }$$

Algorithm 3. Implementation of condition C.

Given the TORs l, m fulfilling condition A,B

1: for Z* in Z1...Z6 do

2:   for A* in A1...A8 do

3:     Compute $\vec {K_*}$ using (4)

4:     if All $\vec {K_*}$ components are integer then

5:      Transform all mi using (5)

6:      Sort l and m TORs according to the pixel index

7:      for i in 1...n0 do

8:        if $l_{1, 2, 3}^{i} \neq m_{1, 2, 3}^{i}$ or $d(l^i_0, m^i_0) >t$ then

9:         $(Z_*, A_*, K_*)$ is not a symmetry of l, m

10:     end if

11:     end for

12:    Symmetry $(Z_*, A_*, K_*)$ found between l, m

13:  end if

14:  end for

15: end for

It is worth noticing that the storage saved depends on the TOR size but not on the symmetry detected.

The IRIS PET/CT is a novel pre-clinical system for mice and rats featuring a full ring PET and a high resolution CT system (Belcari et al 2017). The PET component of the scanner consists of 8 heads arranged in an octagonal ring. Each head comprises two modules, each based on a LYSO matrix of 702 crystals of 1.6 mm  ×  1.6 mm  ×  12 mm directly coupled to a 64 anodes PMT (Hamamatsu H8500). Each head can acquire coincidences with the three heads facing it, therefore the scanner can acquire coincidences from 12 head pairs and features $(702\cdot 2) ^2\cdot 12=23\, 654\, 592$ LORs. The detector implements a component based normalization procedure (Belcari et al 2017).

To show that the SSA can be operated both in an automatic and in a semi-automatic way, two types of experiments were performed in this study. In the first experiment, an SRM, named hereafter SRM_O (Original), including only the TORs corresponding to the head pairs shown in figure 2 was used. In this configuration, the other TORs are obtained by rotating the image, e.g. the TORs belonging to the pair 3–7 can be obtained by rotating the image by $90^{\circ}$ and using the corresponding TORs of the pair 1–5. In this work, a 3rd order b-spline interpolation is used to implement the image rotation. The manual exploitation of this set of rotational symmetries, allows to reduce the SRM computation time and its storage by a factor 6. SRM_O is used in this work as a 'gold standard', i.e. all the reconstructed images are compared to those obtained with SRM_O. In the second experiment, an SRM, named hereafter SRM_C (Complete), that contains all the TORs of the scanner was used. In both cases, the two models were the result of a multi-ray Siddon with $4\times 4$ integration points on the crystal surface and 8 in the crystal depth (Moehrs et al 2008), with the FOV sampled into pixels of 0.855 mm  ×  0.855 mm  ×  1.71 mm. This setup differs considerably from the size of the standard reconstructed image used in pre-clinical routine (0.42 mm  ×  0.42 mm  ×  0.855 mm). The large pixel size was necessary in order to limit the size of SRM_O and being able to run reconstruction without using the symmetries on the 64 GB RAM workstation available for this work. To show that the SSA does not affect the quality of the reconstructed images, a ¹⁸FDG filled NEMA image quality (IQ) phantom (National Electrical Manufacturers Association 2008), with an initial activity of 37 MBq (100 μCi), was acquired for 20 min with the IRIS PET scanner. The IQ was reconstructed performing 100 MLEM iterations, using the models mentioned above and correcting for decay time, dead-time and using detector normalization (Belcari et al 2017). Two analysis were performed on the reconstructed images of the IQ: first, the absolute values of the images reconstructed with the compressed models and those obtained with SRM_O were compared. The differences between the images reconstructed with different models were evaluated using a relative difference metric: for each voxel i and iteration q the following value was computed

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta I^{q}_i= \frac{C^q_i-O^q_i}{O^q_i} \nonumber \end{align} \tag{ 6 }$$

where $O^q_i$ and $C^q_i$ are the values of the ith voxel at the qth iteration obtained with SRM_O and with one of the compressed models (C), respectively. All the relative values reported in this paper are expressed in decimal representation rather than in percentages, e.g. 0.01 instead of 1%. The results of this analysis were resumed in three figures of merit: the maximum relative difference between two images $M_q=\max\nolimits_i(|\Delta I^{q}|)$, the average relative difference $A_q={avg}_i(|\Delta I^{q}|)$ and the standard deviation of the relative difference $S_q={std}_i(|\Delta I^{q}|)$. As a second analysis, an evaluation of the image quality according to the NEMA standard NU-2008 for small animal was performed (National Electrical Manufacturers Association 2008). In order to assess the image quality, the NEMA prescribes to investigate two figures of merit: the recovery coefficients and the image uniformity. For studying the uniformity, a cylindrical region of interest of diameter 22.5 mm and height 10 mm (VOI) was drawn in the uniform region of the phantom. The image uniformity was then evaluated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathrm{IU}=\frac{\mathrm{std(VOI)}}{\mathrm{mean(VOI)}} \nonumber \end{align} \tag{ 7 }$$

where $\mathrm{std(VOI)}$ and $\mathrm{mean(VOI)}$ are the standard deviation and mean of the VOI region, respectively. The recovery coefficient (RC) of the nth rod was evaluated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathrm{RC}_n = \frac{\max (\mathrm{ROD}_n)}{\mathrm{mean(VOI)}} \nonumber \end{align} \tag{ 8 }$$

where $\mathrm{ROD}_n$ is a cylindrical region centered on the nth rod, with two times the diameter of the rod and with a height of 10 mm. The idea behind this analysis is that RCs are correlated with the spatial resolution of the system, whereas the image uniformity, is intended to quantify the spatial noise perceived in an individual image.

Zoom In Zoom Out Reset image size

The sizes of SRM_O and SRM_C are respectively 52 GB and 346 GB. SRM_O and SRM_C were compressed using different thresholds ranging from 0.01 to 2 and discarding condition B, which is equivalent to set $t=\infty$. Each model produced by the SSA was named according to the model compressed and the threshold used, e.g. the name of the model obtained compressing SRM_O with t  =  0.01 is SRM$_{\rm O_{0.01}}$. The results achieved in terms of compression factors and number of fundamental TORs at different thresholds are shown in tables 1 and 2.

Table 1. Size and compression factors (CF) of SRM_O as obtained at different thresholds.

Threshold	SRMO size	Compression factor	# of fundamental TORs
—	52 GB	—	3 942 432
0.01	2.4 GB	22	141 468
0.05	1.3 GB	40	72 786
0.5	1.1 GB	47	57 945
1	1.0 GB	52	56 975
2	1.0 GB	52	56 344
$\infty$	984 MB	54	55 449

Table 2. Size and compression factors (CF) of SRM_C as obtained at different thresholds.

Threshold	SRMC size (GB)	Compression factor	# of Fundamental TORs
—	346	—	23 654 592
0.01	11	31	605 957
0.05	3.6	96	194 819
0.5	1.9	182	106 745
1	1.8	192	101 090
2	1.8	192	98 145
$\infty$	1.7	204	94 015

The SSA was implemented in C++ and was run on an Intel(R) Xeon(R) @ 3.50 GHz. The time needed to compress a SRM depends on the threshold used and reaches a maximum when $t=\infty$ is used. The maximum compression time was found to be 3000 s for SRM_O and 15 000 s for SRM_C. These numbers have to be compared with the time needed to compute SRM_O and SRM_C, which are 42 000 s and 250 000 s, respectively. All these results refer to single core implementations of the SSA. To quantify the overhead due to the implementation of the symmetries detected by the SSA in the reconstruction software, two multi-core MLEM algorithms were implemented: the first does not make use of the SSA detected symmetries, the second implements equation (5) before projection and retro-projection operations. The reconstructions were run on the same workstation used for the compression tests. The time needed to perform a single MLEM iteration using any of the models compressed with the SSA was 67 s, regardless of the threshold used to produce the model. The time needed to perform the a single MLEM iteration using SRM_O, i.e. without the use of SSA detected symmetries, was 56 s.

The plots of M_q, A_q, S_q versus the number of iterations (q) for some of the models produced by compressing SRM$_{\mathrm{O}}$ and SRM$_{\mathrm{C}}$ are reported in figure 3. As an example, the distribution of values of $\Delta I^{100}$ are shown in figure 4 for four models obtained compressing SRM$_{\mathrm{O}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The maximum recovery coefficients MRC obtained within the first 100 iterations and the iteration where the maximum was obtained are reported in table 3. The results of all the models compressed using a finite threshold are not reported as they were exactly the same as those obtained with SRM$_{\mathrm{O}}$. An IU of 6.75% at the 100th MLEM iteration was measured for all the models.

Table 3. Maximum recovery coefficients of the NEMA IQ phantom obtained using SRM$_{\mathrm{O}}$, SRM$_{\mathrm{O}_{\infty}}$ and SRM$_{\mathrm{C}_{\infty}}$.

Rod / MRC (mm)	SRM$_{\mathrm{O}}$ (iteration)	SRM$_{\mathrm{O}_{\infty}}$	SRM$_{\mathrm{C}_{\infty}}$ (iteration)
1	0.21  ±  0.03 (100)	0.21  ±  0.03 (100)	0.21  ±  0.03 (100)
2	0.50  ±  0.07 (100)	0.50  ±  0.07 (100)	0.51  ±  0.07 (100)
3	0.81  ±  0.06 (47)	0.80  ±  0.06 (47)	0.81  ±  0.06 (46)
4	0.93  ±  0.05 (25)	0.92  ±  0.05 (25)	0.93  ±  0.05 (25)
5	0.91  ±  0.05 (19)	0.91  ±  0.05 (19)	0.91  ±  0.05 (18)

An experiment using the model of the YAP-PET small animal scanner was performed. The model of the YAP-PET was computed as described in Motta et al (2002) and its size was 5.9 GB with the FOV segmented into $128\times 128 \times 80$ voxels of $0.375 \times 0.375 \times 0.5$ mm³. The maximum compression obtained with the SSA using $t=\infty$ is 92. This compression factor is comparable with that provided by the authors using manual implementation of reflections and translations (∼80).