Fisher information analysis of list-mode SPECT emission data for joint estimation of activity and attenuation distribution

Consider an SPECT system consisting of scintillation cameras that is detecting gamma-ray photons emitted by an object. This object consists of an activity distribution and an attenuation distribution, which parameterize the emission and attenuation of photons, respectively. Consider that the object being imaged is represented in a voxel basis consisting of N voxels, so that the activity and attenuation distributions are a set of N-dimensional vectors, denoted by λ and μ, respectively. Denote the N elements of the activity and attenuation distribution vector, by $\left\{{\lambda }_{1},\dots ,{\lambda }_{N}\right\}$ and $\left\{{\mu }_{1},\dots ,{\mu }_{N}\right\}$, respectively, where λ_i is themean number of photons emitted from the ith voxel per unit time and μ_i is the attenuation coefficient of the ith voxel. When a transmission scan is used to measure the attenuation distribution μ, then the inverse problem is to estimate λ. However, in this paper, we are considering that a transmission scan is unavailable, so the inverse problem is jointly estimating λ and μ from only the SPECT emission data.

For each gamma-ray photon emitted from the object that interacts with the scintillator, the position of interaction of the gamma-ray photon with the crystal and the energy deposited at the interaction site are estimated and recorded. Denote the true and estimated attributes of the jth detected event by the attribute vectors A_j and ${\hat{\mathbf{A}}}_{j}$, respectively. The attribute vector is a q-dimensional vector where q is the number of attributes in each LM event. We assume that the number of detected counts J is fixed, i.e. we have a preset–count system. Note that our analysis is general and can be easily extended to a preset–time system. Denote the full LM dataset of estimated attributes as $\hat{\mathcal{A}}=\left\{{\hat{\mathbf{A}}}_{j},j=1,2,\enspace \dots ,\enspace J\right\}$.

We can represent the LM data as an impulse-valued random process on attribute space [35], given by the generalized function:

$$\begin{equation}g\left(\mathbf{A}\right)=\sum _{j=1}^{J}\delta \left(\mathbf{A}-{\hat{\mathbf{A}}}_{j}\right).\end{equation} \tag{ 1 }$$

As mentioned above, this data can be binned and for completeness, we present the mathematical representation of the binned data. A bin can be defined as a hyperrectangle in the attribute space, denoted by the function b_m(A) for m = 1, ..., M. Here we assume that the binningfunctions represent non-overlapping hyperrectangle and encompass the full range of attribute space. In that case, the measurement in the mth bin, denoted by g_m is given by [36]

$$\begin{equation}{g}_{m}=\int g\left(\mathbf{A}\right){b}_{m}\left(\mathbf{A}\right)\enspace {\mathrm{d}}^{q}\mathbf{A}=\sum _{j=1}^{J}{b}_{m}\left({\hat{\mathbf{A}}}_{j}\right),\end{equation} \tag{ 2 }$$

where the integral is over the attribute space. We observe that the binned data is the result of an additional binning function applied to the LM data and thus intuitively is expected to lead to an information loss, as also mentioned above. This has also been observed in previousstudies [31, 32, 36]. Thus, we focus our analysis on jointly estimating the activity and attenuation distribution directly from the LM data.

More specifically, we intend to quantitatively determine if the LM data contains information to jointly estimate λ and μ. We use the widely used metric of Fisher information to quantify this information content [35, 37]. Deriving the Fisher information requires computing an expression for the likelihood. This is given by [30]

$$\begin{equation}\mathrm{p}\mathrm{r}\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)=\mathrm{p}\mathrm{r}\left(T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)\prod _{j=1}^{J}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right),\end{equation} \tag{ 3 }$$

where T is the acquisition time required for detecting J LM events. In the above equation, we have used the fact that the acquisition time and each of the list-mode events are independent of each other. Taking the logarithm on both sides of the equation (3) yields the log-likelihood of the observed LM data, denoted by $\mathcal{L}\left(\boldsymbol{\lambda },\boldsymbol{\mu }\vert \hat{\mathcal{A}},T\right)$ and given by

$$\begin{equation}\mathcal{L}\left(\boldsymbol{\lambda },\boldsymbol{\mu }\vert \hat{\mathcal{A}},T\right)=\sum _{j=1}^{J}\mathrm{log}\;\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)+\mathrm{log}\;\mathrm{p}\mathrm{r}\left(T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right).\end{equation} \tag{ 4 }$$

To quantify the Fisher information in the LM data for jointly estimating λ and μ, the log-likelihood must be differentiated with respect to the elements of λ and μ. This requiresobtaining an analytic expression for $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$ and pr(T|λ, μ). For a fixed number of detected counts J, the acquisition time T follows an Erlang distribution with shape J and rate β, where β is the mean rate of photons detected by the detector [38], so that

$$\begin{equation}\mathrm{p}\mathrm{r}\left(T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)=\frac{{\beta }^{J}{T}^{J-1}\mathrm{exp}\left(-\beta T\right)}{\left(J-1\right)!}.\end{equation} \tag{ 5 }$$

Obtaining a similar direct analytic expression for $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$ is complicated. To address this issue, note that any LM event is the result of a photon emitted from a voxel, traveling in a certain direction and then, in some cases, scattering in certain voxels, and finally being detected by the detector. In other words, any LM event is a result of a photon traveling within a discrete unit of space, which we refer to as a path. The concept of a path will be mathematically defined in the next section, but for now, suffice to say that a path is a discrete variable, denoted by $\mathbb{P}$. The expressions for the PDF of the LM event given the path, $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)$ and the probability mass function of the path, $\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$, can be derived, as described later. We thus decompose $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$ as a mixture model over all possible paths. For this purpose, we use the following identity:

$$\begin{equation}\mathrm{p}\mathrm{r}\left(x\vert y=Y\right)=\sum _{z}\mathrm{p}\mathrm{r}\left(x\vert y=Y,z=Z\right)\mathrm{Pr}\left(z=Z\vert y=Y\right),\end{equation} \tag{ 6 }$$

where x denotes a continuous random variable, and y and z denote discrete random variables. In the considered scenario, x, y and z correspond to the LM event attributes, emission and attenuation vectors and the path, respectively. Applying this identity yields the following mixture model for $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$:

$$\begin{equation}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)=\sum _{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P},\boldsymbol{\lambda },\boldsymbol{\mu }\right)\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right).\end{equation} \tag{ 7 }$$

The components of the mixture model are the probabilities that an LM event occurs given the photon traces a path and the weight of each component is the probability of the considered path. Because the LM event has already occurred, the probability of the event given thepath is independent of the activity and attenuation distribution, i.e. $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P},\boldsymbol{\lambda },\boldsymbol{\mu }\right)=\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)$, provided the probability of the path accounts for the emission and attenuation processes, as will be the case in our treatment. Using equation (4), we can rewrite the log-likelihood of the data given the activity and attenuation distribution in terms of this mixture-model decomposition as

$$\begin{equation}\mathcal{L}\left(\boldsymbol{\lambda },\boldsymbol{\mu }\vert \hat{\mathcal{A}},T\right)=\sum _{j=1}^{J}\mathrm{log}\sum _{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)+\mathrm{log}\enspace \mathrm{p}\mathrm{r}\left(T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right).\end{equation} \tag{ 8 }$$

To derive the expression for the Fisher information, analytic expressions for $\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$ and $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)$ must be derived. These are the topics of the next two sub-sections.

In this section, we derive the expression for $\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$. We first mathematically define a path. A path is a discrete unit of space that connects different voxels through which photon radiation propagates. A path can be described in terms of a set of sub-paths, where a sub-path describes the unit of space through which radiation propagates between two voxels. To describe the radiation transfer through a sub-path, we use an approach similar to the discrete-ordinates method for solving the equation of radiative transport [39]. Each sub-path is defined in terms of a start location and a finite angular range. First, assume that the directional coordinates are discretized by dividing the angular space of 4π radians into a finite number of solid-angle sub-domains, referred to as ordinates, each of size ΔΩ. Denote the 3D unitary direction vector by $\hat{\boldsymbol{s}}$ and discretized angular direction by ${\hat{\boldsymbol{s}}}_{k}$ associated with the kth instance of the ordinate. Then the indicator function of the kth angular ordinate ${\psi }_{k}\left(\hat{\boldsymbol{s}}\right)$ is defined as below:

$$\begin{equation}{\psi }_{k}\left(\hat{\boldsymbol{s}}\right)=\begin{cases}1,\quad \text{if}\;k\text{th}\;\text{ordinate}\;\text{contains}\;\text{the}\;\text{direction}\;\text{vector}\;\hat{\boldsymbol{s}}.\quad \hfill \\ 0,\quad \text{otherwise.}\quad \hfill \end{cases}\end{equation} \tag{ 9 }$$

A sub-path ${\mathbb{S}}_{ik}$ is now defined as a cone with a vertex at the center of the ith voxel and an angular range given by ${\psi }_{k}\left(\hat{\boldsymbol{s}}\right)$. In our analysis, we consider a path between voxels withindices i₀, i₁, ..., i_n. This path can be considered to consist of n sub-paths between i₀ to i₁, i₁ to i₂ and so on until i_n−1 to i_n. A schematic describing the above notations and other notations is presented in figure 2.

Zoom In Zoom Out Reset image size

The probability of a particular path $\mathbb{P}$, i.e. $\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$, is the ratio of the flux incident on the detector through the path $\mathbb{P}$ to the flux incident on the detector through all possible paths.Let ${\Phi}\left(\mathbb{P}\right)$ denote the flux of photons incident on the detector through a path $\mathbb{P}$. Then,

$$\begin{equation}\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)=\frac{{\Phi}\left(\mathbb{P}\right)}{{\sum }_{{\mathbb{P}}^{\prime }}{\Phi}\left({\mathbb{P}}^{\prime }\right)}.\end{equation} \tag{ 10 }$$

To derive the expression for ${\Phi}\left(\mathbb{P}\right)$, we first define a few other radiometric quantities. Let $\mathcal{Q}$ denote the radiant energy, the 3D vector r denote a location in space and E denote energy. The fundamental radiometric quantity we use to describe the photon transport is the photon distribution function $w\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)$, given by

$$\begin{equation}w\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{1}{E}\frac{{\partial }^{3}\mathcal{Q}}{\partial V\partial {\Omega}\partial E}.\end{equation} \tag{ 11 }$$

The quantity used to describe emission of photons is the source distribution function, denoted by ${\Xi}\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)$, and defined as

$$\begin{equation}{\Xi}\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{{\partial }^{3}{\Phi}}{\partial V\partial {\Omega}\partial E}.\end{equation} \tag{ 12 }$$

Finally, the radiant intensity, denoted by ${\Gamma}\left(\hat{\boldsymbol{s}}\right)$ is defined as

$$\begin{equation}{\Gamma}\left(\hat{\boldsymbol{s}}\right)={c}_{m}\iint w\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)\enspace {\mathrm{d}}^{3}\boldsymbol{r}\enspace \mathrm{d}E.\end{equation} \tag{ 13 }$$

Consider a point source at location r_s in the voxel indexed by i₀ emitting photons of energy E₀ at a constant rate ${\lambda }_{{i}_{0}}$. Assuming isotropic emission, the source distribution along the sub-path ${\mathbb{S}}_{{i}_{0},{k}_{0}}$ is given by

$$\begin{equation}{\Xi}\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{{\lambda }_{{i}_{0}}}{4\pi }\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}\right)\delta \left(E-{E}_{0}\right){\psi }_{{k}_{0}}\left(\hat{\boldsymbol{s}}\right),\end{equation} \tag{ 14 }$$

where δ(x) is the Dirac delta function. As photons travel from the voxel i₀ to i₁, a fraction of the photons scatter, leading to attenuation of photons along this path. The effect of the attenuation transform on the distribution function w is given by [35]

$$\begin{equation}\left[\mathcal{X}w\right]\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{1}{{c}_{m}}{\int }_{0}^{\infty }w\left(\boldsymbol{r}-\hat{\boldsymbol{s}}l\right)\mathrm{exp}\left[-{\int }_{0}^{l}\mu \left(\boldsymbol{r}-\hat{\boldsymbol{s}}{l}^{\prime },E\right)\mathrm{d}{l}^{\prime }\right]\mathrm{d}l,\end{equation} \tag{ 15 }$$

where μ(r, E) denotes the attenuation coefficient at location r and energy E, and c_m denotes the speed of light. Applying the attenuation transform to the source distribution function (equation (14)) yields

$$\begin{equation}\left[\mathcal{X}{\Xi}\right]\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{{\lambda }_{{i}_{0}}}{4\pi {c}_{m}}\delta \left(E-{E}_{0}\right){\psi }_{{k}_{0}}\left(\hat{\boldsymbol{s}}\right){\int }_{0}^{\infty }\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}-\hat{\boldsymbol{s}}l\right)\mathrm{exp}\left[-{\int }_{0}^{l}\mu \left(\boldsymbol{r}-\hat{\boldsymbol{s}}{l}^{\prime },E\right)\mathrm{d}{l}^{\prime }\right]\mathrm{d}l.\end{equation} \tag{ 16 }$$

Now, it can be shown that [35]

$$\begin{equation}{\int }_{l=0}^{\infty }\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}-\hat{\boldsymbol{s}}l\right)\mathrm{d}l=\frac{1}{\vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}{\vert }^{2}}\delta \left(\hat{\boldsymbol{s}}-{\hat{\boldsymbol{s}}}_{10}\right),\end{equation} \tag{ 17 }$$

where

$$\begin{equation}{\hat{\boldsymbol{s}}}_{10}=\frac{\boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}}{\vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}\vert }.\end{equation} \tag{ 18 }$$

Using the above relation and the sifting property of the delta function, the integral over l in equation (16) can be simplified to yield

$$\begin{equation}\left[\mathcal{X}{\Xi}\right]\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{{\lambda }_{{i}_{0}}}{4\pi {c}_{m}\vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}{\vert }^{2}}\mathrm{exp}\left[-{\int }_{0}^{\vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}\vert }\mu \left(\boldsymbol{r}-\hat{\boldsymbol{s}}{l}^{\prime },E\right)\enspace \mathrm{d}{l}^{\prime }\right]\delta \left(\hat{\boldsymbol{s}}-{\hat{\boldsymbol{s}}}_{10}\right)\delta \left(E-{E}_{0}\right){\psi }_{{k}_{0}}\left(\hat{\boldsymbol{s}}\right).\end{equation} \tag{ 19 }$$

In the considered path, Compton scattering occurs at some location r within voxel i₁. This operation can be described in terms of the scattering operator $\mathcal{K}$ as

$$\begin{equation}\left[\mathcal{K}w\right]\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)={\iint }_{4\pi }K\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}^{\prime },E,{E}^{\prime }\vert \boldsymbol{r}\right)w\left(\boldsymbol{r},{\hat{\boldsymbol{s}}}^{\prime },{E}^{\prime }\right)\enspace \mathrm{d}{{\Omega}}^{\prime }\enspace \mathrm{d}{E}^{\prime },\end{equation} \tag{ 20 }$$

where $K\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}^{\prime },E,{E}^{\prime }\vert \boldsymbol{r}\right)$ denotes the scattering kernel and is given by

$$\begin{equation}K\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}^{\prime },E,{E}^{\prime }\vert \boldsymbol{r}\right)={c}_{m}{n}_{\text{sc}}\left(\boldsymbol{r}\right)\frac{\partial {\sigma }_{\mathrm{s}\mathrm{c}}}{\partial {\Omega}}\delta \left\{E-{\left[\frac{1}{{E}^{\prime }}+\frac{1}{m{c}^{2}}\left(1-\mathrm{cos}\enspace \theta \right)\right]}^{-1}\right\},\end{equation} \tag{ 21 }$$

where n_sc(r) is the density of scatters and is related to the scattering coefficient by

$$\begin{equation}\mu \left(\boldsymbol{r}\right)={n}_{\text{sc}}\left(\boldsymbol{r}\right){\sigma }_{\text{sc}},\end{equation} \tag{ 22 }$$

where σ_sc is the scattering cross section. Also, θ is the angle at which the outgoing photon scatters relative to the incoming photon, so $\mathrm{cos}\enspace \theta =\hat{\boldsymbol{s}}\cdot {\hat{\boldsymbol{s}}}^{\prime }$. Finally, the differential scattering cross section $\frac{\partial {\sigma }_{\mathrm{s}\mathrm{c}}}{\partial {\Omega}}$ in equation (21) is given by the Klein–Nishina formula [17]. For notational simplicity, define ${K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}^{\prime },E\vert \boldsymbol{r}\right)$ by

$$\begin{equation}{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}^{\prime },E\vert \boldsymbol{r}\right)={n}_{\text{sc}}\left(\boldsymbol{r}\right){\left.\frac{\partial {\sigma }_{\mathrm{s}\mathrm{c}}}{\partial {\Omega}}\right\vert }_{\mathrm{cos}\theta =\hat{\boldsymbol{s}}\cdot {\hat{\boldsymbol{s}}}^{\prime }}.\end{equation} \tag{ 23 }$$

Also, denote the path integral between any two locations r_l and r_m by the function γ(r_l, r_m, E), i.e.

$$\begin{equation}\gamma \left({\boldsymbol{r}}_{l},{\boldsymbol{r}}_{m},E\right)={\int }_{0}^{\vert {\boldsymbol{r}}_{l}-{\boldsymbol{r}}_{m}\vert }\mu \left({\boldsymbol{r}}_{l}-t\frac{{\boldsymbol{r}}_{1}-{\boldsymbol{r}}_{m}}{\vert {\boldsymbol{r}}_{1}-{\boldsymbol{r}}_{m}\vert },E\right)\mathrm{d}t.\end{equation} \tag{ 24 }$$

Substituting the expression for the attenuation transformed source distribution function from equation (19) into equation (20) and using the sifting property of the delta function in angular coordinates yields

$$\begin{equation}\left[\mathcal{K}\mathcal{X}{\Xi}\right]\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{{\lambda }_{{i}_{0}}}{4\pi {c}_{m}\vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}{\vert }^{2}}K\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{10},E,{E}_{0}\vert \boldsymbol{r}\right)\mathrm{exp}\left[-\gamma \left(\boldsymbol{r},{\boldsymbol{r}}_{\mathrm{s}},{E}_{0}\right)\right]{\psi }_{{k}_{0}}\left({\hat{\boldsymbol{s}}}_{10}\right).\end{equation} \tag{ 25 }$$

We now integrate this distribution function over all possible locations within the ${i}_{1}^{\text{th}}$ voxel and over all possible energies. This yields the radiant intensity along direction $\hat{\boldsymbol{s}}$ due to photons traveling through the ${\mathbb{S}}_{{i}_{0},{k}_{0}}$ sub-path and scattering within the ${i}_{1}^{\text{th}}$ voxel. Denote this radiant intensity by ${{\Gamma}}_{{i}_{0}{i}_{1}{k}_{0}}\left(\hat{\boldsymbol{s}}\right)$. Substituting the distribution function from equation (25) into equation (13) yields

$$\begin{equation}{{\Gamma}}_{{i}_{0}{i}_{1}{k}_{0}}\left(\hat{\boldsymbol{s}}\right)=\int \frac{{\lambda }_{{i}_{0}}}{4\pi \vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}{\vert }^{2}}{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{10},{E}_{0}\vert \boldsymbol{r}\right)\mathrm{exp}\left[-\gamma \left(\boldsymbol{r},{\boldsymbol{r}}_{\mathrm{s}},{E}_{0}\right)\right]{\phi }_{{i}_{1}}\left(\boldsymbol{r}\right){\psi }_{{k}_{0}}\left({\hat{\boldsymbol{s}}}_{10}\right)\enspace {\mathrm{d}}^{3}\boldsymbol{r},\end{equation} \tag{ 26 }$$

where ϕ_i(r) is the voxel basis function defined as below:

$$\begin{equation}{\phi }_{i}\left(\boldsymbol{r}\right)=\begin{cases}1,\quad \boldsymbol{r}\enspace \text{lies}\;\text{within}\;\text{voxel}\;\text{i.}\quad \hfill \\ 0,\quad \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\quad \hfill \end{cases}\end{equation} \tag{ 27 }$$

Assuming that the functions $K\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{10},E,{E}_{0}\right)$ and $\mathrm{exp}\left[-\gamma \left(\boldsymbol{r},{\boldsymbol{r}}_{\mathrm{s}},{E}_{0}\right)\right]$ do not vary relatively within any location within voxel i₁, we can evaluate them when r is the center of the ${i}_{1}^{\text{th}}$ voxel and r_s is the center of the ${i}_{0}^{\text{th}}$ voxel. Denoting the center of the i₀ and i₁ voxels by r₀ and r₁, respectively, and denoting the direction vector joining r₀ and r₁ by ${\hat{\boldsymbol{s}}}_{c10}$, we obtain

$$\begin{equation}{{\Gamma}}_{{i}_{0}{i}_{1}{k}_{0}}\left(\hat{\boldsymbol{s}}\right)\approx \frac{{\lambda }_{{i}_{0}}}{4\pi }\mathrm{exp}\left[-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)\right]{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{c10},{E}_{0}\vert {\boldsymbol{r}}_{1}\right)\int \frac{1}{\vert \boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}{\vert }^{2}}{\psi }_{{k}_{0}}\left({\hat{\boldsymbol{s}}}_{10}\right){\phi }_{{i}_{1}}\left(\boldsymbol{r}\right)\enspace {\mathrm{d}}^{3}\boldsymbol{r}.\end{equation} \tag{ 28 }$$

Using the definition of ${\hat{\boldsymbol{s}}}_{10}$ from equation (18) and denoting |r − r_s| by R, we perform a change of variables $\boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{s}}=R{\hat{\boldsymbol{s}}}_{10}$ and replace r by ${R}^{2}\mathrm{d}R\mathrm{d}{\hat{\boldsymbol{s}}}_{10}$. Simplifying further yields

$$\begin{equation}{{\Gamma}}_{{i}_{0}{i}_{1}{k}_{0}}\left(\hat{\boldsymbol{s}}\right)\approx \frac{{\lambda }_{{i}_{0}}}{4\pi }\mathrm{exp}\left[-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)\right]{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{c10},{E}_{0}\vert {\boldsymbol{r}}_{1}\right)\int {\psi }_{{k}_{0}}\left({\hat{\boldsymbol{s}}}_{10}\right)\int {\phi }_{{i}_{1}}\left({\boldsymbol{r}}_{\mathrm{s}}+R{\hat{\boldsymbol{s}}}_{10}\right)\enspace \mathrm{d}R\enspace \mathrm{d}{\hat{\boldsymbol{s}}}_{10}.\end{equation} \tag{ 29 }$$

The integral over R is equal to the distance traversed by the ${\boldsymbol{r}}_{\mathrm{s}}+R{\hat{\boldsymbol{s}}}_{10}$ vector within the ${k}_{1}^{\text{th}}$ voxel, so this distance should vary with ${\hat{\boldsymbol{s}}}_{10}$. However, assuming that this variation is not substantial, we approximate it by the distance that is covered by the vector ${\boldsymbol{r}}_{0}+R{\hat{\boldsymbol{s}}}_{{k}_{0}}$ in the ${i}_{1}^{\text{th}}$ voxel, which we denote by ${{\Delta}}_{{i}_{1}}\left({\mathbb{S}}_{{i}_{0}{k}_{0}}\right)$. Note that ${\hat{\boldsymbol{s}}}_{{k}_{0}}$ is the discretized direction vector associated with the ${k}_{0}^{\text{th}}$ ordinate. Further, performing the integral over ${\hat{\boldsymbol{s}}}_{10}$ yields the following expression for${{\Gamma}}_{{i}_{0}{i}_{1}{k}_{0}}\left(\hat{\boldsymbol{s}}\right)$:

$$\begin{equation}{{\Gamma}}_{{i}_{0}{i}_{1}{k}_{0}}\left(\hat{\boldsymbol{s}}\right)\approx \frac{{\lambda }_{{i}_{0}}}{4\pi }\mathrm{exp}\left[-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)\right]{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{c10},{E}_{0}\vert {\boldsymbol{r}}_{1}\right){{\Delta}}_{{i}_{1}}\left({\mathbb{S}}_{{i}_{0}{k}_{0}}\right){\Delta}{\Omega}.\end{equation} \tag{ 30 }$$

To proceed further, for mathematical tractability, assume that this entire radiant intensity is concentrated at the center of the ${i}_{1}^{\text{th}}$ voxel, i.e. at location r₁ and divided uniformly over the sub-path from this voxel that includes the direction $\hat{\boldsymbol{s}}$. Denote this sub-path by ${\mathbb{S}}_{{i}_{1},{k}_{1}}$. Thus, the distribution function along this sub-path, denoted by ${w}_{1}\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)$, is given by

$$\begin{equation}{w}_{1}\left(\boldsymbol{r},\hat{\boldsymbol{s}},E\right)=\frac{{\lambda }_{{i}_{0}}}{4\pi {c}_{m}}\mathrm{exp}\left[-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)\right]K\left(\hat{\boldsymbol{s}},{\hat{\boldsymbol{s}}}_{c10},E,{E}_{0}\vert {\boldsymbol{r}}_{1}\right){{\Delta}}_{{i}_{1}}\left({\mathbb{S}}_{{i}_{0}{k}_{0}}\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{1}\right){\psi }_{{k}_{1}}\left(\hat{\boldsymbol{s}}\right).\end{equation} \tag{ 31 }$$

After this scattering operation, the photons along this path suffer from attenuation as they travel toward voxel i₂.

This series of operations continues until the path intersects with the detector. In the considered path, the last voxel where scattering occurs is i_n and the sub-path that connects the last voxel to the detector is denoted by ${\mathbb{S}}_{{i}_{n},{k}_{n}}$. Denote the coordinates on the front face of the detector by r_d and the normal unit vector to the detector plane by $\hat{\boldsymbol{n}}$. Then the distribution function at the face of the detector, denoted by ${w}_{d}\left({\boldsymbol{r}}_{d},\hat{\boldsymbol{s}},E\right)$ is given by

$$\begin{align}\hfill {w}_{d}\left({\boldsymbol{r}}_{d},\hat{\boldsymbol{s}},E\right)& =\frac{{\prod }_{m=0}^{n-1}{{\Delta}}_{{i}_{m+1}}\left({\mathbb{S}}_{{i}_{m}{k}_{m}}\right)}{\vert {\boldsymbol{r}}_{d}-{\boldsymbol{r}}_{n}{\vert }^{2}}\mathrm{exp}\left\{-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)-\dots \gamma \left({\boldsymbol{r}}_{d},{\boldsymbol{r}}_{n},{E}_{n}\right)\right\}\hfill \\ \hfill & \quad {\times}\prod _{m=0}^{n-1}\left\{{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left({\hat{\boldsymbol{s}}}_{c,m+2,m+1},{\hat{\boldsymbol{s}}}_{c,m+1,m},{E}_{m}\vert {\boldsymbol{r}}_{m+1}\right)\right\}\delta \left(E-{E}_{n}\right){\psi }_{{k}_{n}}\left(\hat{\boldsymbol{s}}\right).\hfill \end{align} \tag{ 32 }$$

Now the plane of the front face of the detector is given by $\delta \left(p-{\boldsymbol{r}}_{d}\cdot \hat{\boldsymbol{n}}\right)$, where p is the perpendicular distance from the origin to the detector plane. Let the sensitivity of the collimator to a photon emitted from location r in direction $\hat{\boldsymbol{s}}$ be denoted by $t\left(\boldsymbol{r},\hat{\boldsymbol{s}}\right)$. Then, the flux detected through the considered path is given by

$$\begin{align}\hfill {\Phi}\left(\mathbb{P}\right)& ={c}_{m}\iiint t\left({\boldsymbol{r}}_{n},\hat{\boldsymbol{s}}\right)\left(\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{s}}\right)\delta \left(p-{\boldsymbol{r}}_{d}\cdot \hat{\boldsymbol{n}}\right){w}_{d}\left({\boldsymbol{r}}_{d},\hat{\boldsymbol{s}},E\right)\enspace {\mathrm{d}}^{3}{\boldsymbol{r}}_{d}\enspace \mathrm{d}{\Omega}\enspace \mathrm{d}E\hfill \\ \hfill & =\frac{{\lambda }_{{i}_{0}}}{4\pi }{\iint }_{2\pi }\int t\left({\boldsymbol{r}}_{n},\hat{\boldsymbol{s}}\right)\left(\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{s}}\right)\delta \left(p-{\boldsymbol{r}}_{d}\cdot \hat{\boldsymbol{n}}\right)\mathrm{exp}\left\{-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)-\dots \gamma \left({\boldsymbol{r}}_{d},{\boldsymbol{r}}_{n},{E}_{n}\right)\right\}\hfill \\ \hfill & \quad {\times}\enspace \frac{{\prod }_{m=0}^{n-1}{{\Delta}}_{{i}_{m+1}}\left({\mathbb{S}}_{{i}_{m}{k}_{m}}\right)}{\vert {\boldsymbol{r}}_{d}-{\boldsymbol{r}}_{n}{\vert }^{2}}\prod _{m=0}^{n-1}\left\{{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left({\hat{\boldsymbol{s}}}_{c,m+2,m+1},{\hat{\boldsymbol{s}}}_{c,m+1,m},{E}_{m}\vert {\boldsymbol{r}}_{m+1}\right)\right\}\hfill \\ \hfill & \quad {\times}\enspace \delta \left(E-{E}_{n}\right){\psi }_{{k}_{n}}\left(\hat{\boldsymbol{s}}\right)\enspace {\mathrm{d}}^{3}{\boldsymbol{r}}_{d}\enspace \mathrm{d}{\Omega}\enspace \mathrm{d}E.\hfill \end{align} \tag{ 33 }$$

Perform a change of variables by replacing r_d − r_n by $R\hat{\boldsymbol{s}}$, so that ${\mathrm{d}}^{3}{r}_{d}={R}^{2}\mathrm{d}R\mathrm{d}\hat{\boldsymbol{s}}$. Next, integrating over E and Ω using the sifting property of the delta function yields

$$\begin{align}\hfill {\Phi}\left(\mathbb{P}\right)& =\frac{{\lambda }_{{i}_{0}}}{4\pi }\iint \delta \left(p-{\boldsymbol{r}}_{n}\cdot \hat{\boldsymbol{n}}-R\hat{\boldsymbol{s}}\cdot \hat{\boldsymbol{n}}\right)t\left({\boldsymbol{r}}_{n},\hat{\boldsymbol{s}}\right)\left(\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{s}}\right)\mathrm{exp}\left\{-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)-\dots \gamma \left({\boldsymbol{r}}_{n}+R\hat{\boldsymbol{s}},{\boldsymbol{r}}_{n},{E}_{n}\right)\right\}\hfill \\ \hfill & \quad {\times}\enspace \prod _{m=0}^{n-1}\left\{{{\Delta}}_{{i}_{m+1}}\left({\mathbb{S}}_{{i}_{m}{k}_{m}}\right)\right\}\prod _{m=0}^{n-1}\left\{{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left({\hat{\boldsymbol{s}}}_{c,m+2,m+1},{\hat{\boldsymbol{s}}}_{c,m+1,m},{E}_{m}\vert {\boldsymbol{r}}_{m+1}\right)\right\}\hfill \\ \hfill & \quad {\times}\enspace {\psi }_{{k}_{n}}\left(\hat{\boldsymbol{s}}\right)\enspace \mathrm{d}\hat{\boldsymbol{s}}\enspace \mathrm{d}R.\hfill \end{align} \tag{ 34 }$$

The above expression formalizes the radiation through a path for a general SPECT system. We now derive the specific form of this expression for an SPECT system with a parallel-hole collimator. This collimator allows only photons that are incident in a small range of angles around $\hat{\boldsymbol{n}}$ to pass through the collimator. Thus, assuming $\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{s}}\approx 1$ when $t\left(\boldsymbol{r},\hat{\boldsymbol{s}}\right){ >}0$, the integral over R yields

$$\begin{align}\hfill {\Phi}\left(\mathbb{P}\right)& =\frac{{\lambda }_{{i}_{0}}}{4\pi }\int t\left({\boldsymbol{r}}_{n},\hat{\boldsymbol{s}}\right)\mathrm{exp}\left\{-\gamma \left({\boldsymbol{r}}_{1},{\boldsymbol{r}}_{0},{E}_{0}\right)-\dots \gamma \left({\boldsymbol{r}}_{{d}_{0}},{\boldsymbol{r}}_{n},{E}_{n}\right)\right\}\prod _{m=0}^{n-1}\left\{{{\Delta}}_{{i}_{m+1}}\left({\mathbb{S}}_{{i}_{m}{k}_{m}}\right)\right\}\hfill \\ \hfill & \quad {\times}\enspace \prod _{m=0}^{n-1}\left\{{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left({\hat{\boldsymbol{s}}}_{c,m+2,m+1},{\hat{\boldsymbol{s}}}_{c,m+1,m},{E}_{m}\vert {\boldsymbol{r}}_{m+1}\right)\right\}{\psi }_{{k}_{n}}\left(\hat{\boldsymbol{s}}\right)\enspace \mathrm{d}\hat{\boldsymbol{s}},\hfill \end{align} \tag{ 35 }$$

where ${\boldsymbol{r}}_{d0}={\boldsymbol{r}}_{n}+\left(p-{\boldsymbol{r}}_{n}\cdot \hat{\boldsymbol{n}}\right)\hat{\boldsymbol{s}}$ is the coordinate on the detector plane where the gamma-ray photon is incident. Further, assuming the angular ordinates in the object space to be small compared to the angular range allowed by the collimator, we replace $\hat{\boldsymbol{s}}$ with ${\hat{\boldsymbol{s}}}_{{k}_{n}}$, where ${\hat{\boldsymbol{s}}}_{{k}_{n}}$ denotes the discretized direction vector corresponding to the ${k}_{n}^{\text{th}}$ ordinate. Then the integration over $\hat{\boldsymbol{s}}$ can be performed to yield

$$\begin{align}\hfill {\Phi}\left(\mathbb{P}\right)& =\frac{{\lambda }_{{i}_{0}}{\Delta}{\Omega}}{4\pi }\mathrm{exp}\left\{-\gamma \left({\boldsymbol{r}}_{0},{\boldsymbol{r}}_{1},{E}_{0}\right)-\dots \gamma \left({\boldsymbol{r}}_{dc0},{\boldsymbol{r}}_{n},{E}_{n}\right)\right\}t\left({\boldsymbol{r}}_{n},{\hat{\boldsymbol{s}}}_{{k}_{n}}\right)\prod _{m=0}^{n-1}{{\Delta}}_{{i}_{m+1}}\left({\mathbb{S}}_{{i}_{m}{k}_{m}}\right)\hfill \\ \hfill & \quad {\times}\enspace \prod _{m=0}^{n-1}{K}_{\mathrm{m}\mathrm{a}\mathrm{g}}\left({\hat{\boldsymbol{s}}}_{c,m+2,m+1},{\hat{\boldsymbol{s}}}_{c,m+1,m},{E}_{m}\vert {\boldsymbol{r}}_{m+1}\right),\hfill \end{align} \tag{ 36 }$$

where ${\boldsymbol{r}}_{dc0}={\boldsymbol{r}}_{n}+\left(p-{\boldsymbol{r}}_{n}\cdot \hat{\boldsymbol{n}}\right){\hat{\boldsymbol{s}}}_{{k}_{n}}$. To simplify the expression for ${\Phi}\left(\mathbb{P}\right)$, note that in our analysis, we need to only consider the dependence of this expression on the emission and attenuation coefficients. Quantities that are not dependent on these parameters in the above expression are denoted by ${\Lambda}\left(\mathbb{P}\right)$. Further, to simplify notation, we define ${s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)$ as

$$\begin{equation}{s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)={\Lambda}\left(\mathbb{P}\right)\mathrm{exp}\left[-\sum _{m=0}^{n}\gamma \left({\boldsymbol{r}}_{m},{\boldsymbol{r}}_{m+1},\boldsymbol{\mu }\left({E}_{m}\right)\right)\right]\prod _{m=1}^{n}{\mu }_{{k}_{m}}\left({E}_{m-1}\right).\end{equation} \tag{ 37 }$$

This term actually denotes the effective sensitivity for path $\mathbb{P}$ to the detector. Further, denote the activity in the starting voxel of the path, i.e. ${\lambda }_{{i}_{0}}$, by $\lambda \left(\mathbb{P}\right)$. Equation (36) can then be rewritten as

$$\begin{equation}{\Phi}\left(\mathbb{P}\right)={s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\lambda \left(\mathbb{P}\right).\end{equation} \tag{ 38 }$$

The attenuation coefficient in equation (37) is a function of energy, which does not lend itself to Fisher information matrix analysis with respect to the attenuation coefficient. To address this issue, based on the energy spectrum for common radionuclides such as Technetium-99m (Tc99m), we model the energy dependence of the attenuation coefficient [40] as a linear function of energy, i.e.

$$\begin{equation}\boldsymbol{\mu }\left({E}_{m}\right)=a\boldsymbol{\mu }+b.\end{equation} \tag{ 39 }$$

Here μ denotes the attenuation coefficient at the photo-peak energy i.e. μ = μ(E₀). Also, a and b denote constants that can be computed from the spectrum. Finally, substituting the expression of ${\Phi}\left(\mathbb{P}\right)$ in equation (10) yields:

$$\begin{equation}\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)=\frac{\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}{{\sum }_{{\mathbb{P}}^{\prime }}\lambda \left({\mathbb{P}}^{\prime }\right){s}_{\mathrm{e}\mathrm{ff}}\left({\mathbb{P}}^{\prime }\right)}.\end{equation} \tag{ 40 }$$

This explicit modeling of probability of path in terms of λ and μ is then used to express the log-likelihood and perform the Fisher information analysis, as we describe later.

The term $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)$ denotes the PDF of the measured attribute vector ${\hat{\mathbf{A}}}_{j}$ given the photon followed a particular path $\mathbb{P}$. This attribute vector, as mentioned above, consists of the position of interaction of the gamma-ray photon with the crystal, denoted by ${\hat{\boldsymbol{r}}}_{j}$ and the energydeposited by the photon in the crystal, denoted by ${\hat{E}}_{j}$. To model the randomness in estimating ${\hat{\mathbf{A}}}_{j}$ due to the uncertainty in the estimation process and the finite energy and spatial resolution of the detectors, we first write $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)$ as

$$\begin{equation}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)=\int \mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert {\mathbf{A}}_{j}\right)\mathrm{p}\mathrm{r}\left({\mathbf{A}}_{j}\vert \mathbb{P}\right)\mathrm{d}{\mathbf{A}}_{j}.\end{equation} \tag{ 41 }$$

To obtain the expression for $\mathrm{p}\mathrm{r}\left({\mathbf{A}}_{j}\vert \mathbb{P}\right)$, consider a path that connects several voxels, as in the section above. Consider a photon that propagates exactly between the center of these voxels before reaching the detector. Denote the energy of the photon at the end of the path by${E}_{\mathbb{P}}$ and the location where the photon interacts with the detector by ${\boldsymbol{r}}_{\mathbb{P}}$. Assuming that the paths are discretized finely enough, we can assume that all the photons within this path will have approximately the same energy and interact with the detector at the same location. Thus, we can write

$$\begin{equation}\mathrm{p}\mathrm{r}\left({\mathbf{A}}_{j}\vert \mathbb{P}\right)=\mathrm{p}\mathrm{r}\left(\boldsymbol{r},E\vert \mathbb{P}\right)\approx \delta \left({\mathbf{A}}_{j}-{\mathbf{A}}_{j}\left(\mathbb{P}\right)\right),\end{equation} \tag{ 42 }$$

where ${\mathbf{A}}_{j}\left(\mathbb{P}\right)=\left\{{\boldsymbol{r}}_{\mathbb{P}},{E}_{\mathbb{P}}\right\}$.

Next, to derive the expression for the term $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert {\mathbf{A}}_{j}\right)$, assume that the finite spatial and energy resolution of the detector and the uncertainty in the estimation of the LM attributes can be modeled by a Gaussian distribution. This is a very reasonable assumption for most SPECT systems [31, 41]. Then, we can write

$$\begin{equation}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert {\mathbf{A}}_{j}\right)=\frac{1}{{\left(2\pi \right)}^{2}\vert {\boldsymbol{K}}_{\mathrm{d}\mathrm{e}\mathrm{t}}\vert }\mathrm{exp}\left[-\frac{{\left({\hat{\mathbf{A}}}_{j}-{\mathbf{A}}_{j}\right)}^{{\dagger}}{\boldsymbol{K}}_{\mathrm{d}\mathrm{e}\mathrm{t}}^{-1}\left({\hat{\mathbf{A}}}_{j}-{\mathbf{A}}_{j}\right)}{2}\right],\end{equation} \tag{ 43 }$$

where K_det denotes the covariance matrix quantifying the variances and co-variances in the energy and position estimates and where |K_det| denotes the determinant of the matrix K_det. Substituting equations (43) and (42) into equation (41) and using the sifting property of the delta function yields

$$\begin{equation}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)=\frac{1}{{\left(2\pi \right)}^{2}\vert {\boldsymbol{K}}_{\mathrm{d}\mathrm{e}\mathrm{t}}\vert }\mathrm{exp}\left[-\frac{{\left({\hat{\mathbf{A}}}_{j}-{\mathbf{A}}_{j}\left(\mathbb{P}\right)\right)}^{{\dagger}}{\boldsymbol{K}}_{\mathrm{d}\mathrm{e}\mathrm{t}}^{-1}\left({\hat{\mathbf{A}}}_{j}-{\mathbf{A}}_{j}\left(\mathbb{P}\right)\right)}{2}\right].\end{equation} \tag{ 44 }$$

We now have the expressions for $\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$ (equation (40)) and $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)$ (equation (44)) as required to formalize the likelihood of the LM data. We proceed to deriving the expression for the elements of the Fisher information matrix.

The general expression for the elements of an FIM is given by

$$\begin{equation}{F}_{q{q}^{\prime }}=-{\left\langle \frac{{\partial }^{2}\mathcal{L}\left(\boldsymbol{\lambda },\boldsymbol{\mu }\vert \hat{\mathcal{A}},T\right)}{\partial {\theta }_{q}\partial {\theta }_{{q}^{\prime }}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)},\end{equation} \tag{ 45 }$$

where θ_q and θ_q' denote the parameters we intend to estimate and thus in our case are the activity-attenuation coefficients in the qth and q'th voxels of the object, and where $\mathcal{L}\left(\boldsymbol{\lambda },\boldsymbol{\mu }\vert \hat{\mathcal{A}},T\right)$ is the log-likelihood of the observed LM data (equation (8)). Substituting the expression for $\mathrm{Pr}\left(\mathbb{P}\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)$ from equation (40) into equation (8), and further using equation (5) yields

$$\begin{equation}\mathcal{L}\left(\boldsymbol{\lambda },\boldsymbol{\mu }\vert \hat{\mathcal{A}},T\right)=\sum _{j=1}^{J}\mathrm{log}\left[\sum \mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\right]+\left(J-1\right)\mathrm{log}\enspace T-\beta T-\mathrm{log}\left(J-1\right)!.\end{equation} \tag{ 46 }$$

Here β, which is the mean rate of photons detected, is equivalently the total radiant flux over all possible paths. Thus, β is given by

$$\begin{equation}\beta =\sum _{{\mathbb{P}}^{\prime }}\lambda \left({\mathbb{P}}^{\prime }\right){s}_{\mathrm{e}\mathrm{ff}}\left({\mathbb{P}}^{\prime }\right).\end{equation} \tag{ 47 }$$

Using equation (47) to substitute for β and differentiating the log-likelihood with respect to the activity in the qth voxel λ_q yields

$$\begin{equation}\frac{\partial \mathcal{L}}{\partial {\lambda }_{q}}=\sum _{j=1}^{J}\frac{{\sum }_{{\mathbb{P}}_{q}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}-T\sum _{{\mathbb{P}}_{q}}{s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right),\end{equation} \tag{ 48 }$$

where the summation in the numerator is only over the paths that start from voxel q, denoted by ${\mathbb{P}}_{q}$. Similarly, differentiating the log-likelihood (equation (46) with respect to the attenuation coefficient in the qth voxel, i.e. μ_q, yields

$$\begin{equation}\frac{\partial \mathcal{L}}{\partial {\mu }_{q}}=\sum _{j=1}^{J}\frac{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\left[-{{\Delta}}_{q}\left(\mathbb{P}\right)+\frac{{\zeta }_{q}\left(\mathbb{P}\right)}{{\mu }_{q}}\right]}{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}-T\sum _{\mathbb{P}}\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\left[-{{\Delta}}_{q}\left(\mathbb{P}\right)+\frac{{\zeta }_{q}\left(\mathbb{P}\right)}{{\mu }_{q}}\right],\end{equation} \tag{ 49 }$$

where ${\zeta }_{q}\left(\mathbb{P}\right)$ and ${{\Delta}}_{q}\left(\mathbb{P}\right)$ are the number of scatter events occurring in the qth voxel in the considered path and the distance that the considered path covers in the qth voxel, respectively. To derive the FIM elements, we must differentiate equations (48) and (49) further with respect to the activity and attenuation coefficients in q'th voxel and then average over the observed LM data.The derivations are detailed in appendix A and the final expressions are as below:

$$\begin{align}\hfill {\left\langle \frac{{\partial }^{2}\mathcal{L}}{\partial {\mu }_{{q}^{\prime }}\partial {\mu }_{q}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)}& =-\left\langle \sum _{j=1}^{J}\frac{\left\{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\left[-{{\Delta}}_{q}\left(\mathbb{P}\right)+\frac{{\zeta }_{q}\left(\mathbb{P}\right)}{{\mu }_{q}}\right]\right\}}{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}\right.\hfill \\ \hfill & \quad {\times}{\left.\frac{\left\{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\left[-{{\Delta}}_{{q}^{\prime }}\left(\mathbb{P}\right)+\frac{{\zeta }_{{q}^{\prime }}\left(\mathbb{P}\right)}{{\mu }_{{q}^{\prime }}}\right]\right\}}{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)},\hfill \end{align} \tag{ 50 }$$

$$\begin{equation}{\left\langle \frac{{\partial }^{2}\mathcal{L}}{\partial {\lambda }_{{q}^{\prime }}\partial {\lambda }_{q}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)}=-{\left\langle \sum _{j=1}^{J}\frac{\left\{{\sum }_{{\mathbb{P}}_{q}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\right\}\left\{{\sum }_{{\mathbb{P}}_{{q}^{\prime }}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\right\}}{{\left\{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\right\}}^{2}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)},\end{equation} \tag{ 51 }$$

$$\begin{align}\hfill {\left\langle \frac{{\partial }^{2}\mathcal{L}}{\partial {\mu }_{{q}^{\prime }}\partial {\lambda }_{q}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)}& =-\left\langle \sum _{j=1}^{J}\frac{\left\{{\sum }_{{\mathbb{P}}_{q}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\right\}}{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}\right.\hfill \\ \hfill & \quad {\times}{\left.\frac{\left\{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\left[-{{\Delta}}_{{q}^{\prime }}\left(\mathbb{P}\right)+\frac{{\zeta }_{{q}^{\prime }}\left(\mathbb{P}\right)}{{\mu }_{{q}^{\prime }}}\right]\right\}}{{\sum }_{\mathbb{P}}\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)},\hfill \end{align} \tag{ 52 }$$

$$\begin{equation}{\left\langle \frac{{\partial }^{2}\mathcal{L}}{\partial {\lambda }_{{q}^{\prime }}\partial {\mu }_{q}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)}={\left\langle \frac{{\partial }^{2}\mathcal{L}}{\partial {\mu }_{q}\partial {\lambda }_{{q}^{\prime }}}\right\rangle }_{\left(\hat{\mathcal{A}},T\vert \boldsymbol{\lambda },\boldsymbol{\mu }\right)}.\end{equation} \tag{ 53 }$$

Since we cannot simplify these expressions further, we use Monte Carlo integration to evaluate these expressions from simulated LM data, thus yielding the elements of the FIM. The inverse of the FIM yields the Cramér–Rao bound (CRB), which is the lower bound on thevariance of any unbiased estimate of the activity and attenuation coefficients from the SPECT emission data. Thus, using the CRB, we can quantify the information content of the SPECT emission data on the task of jointly estimating activity and attenuation.

To evaluate the information content in the LM emission data for joint reconstruction, we simulated a clinical 2D parallel-hole SPECT system configuration similar to GE's Optima NM/CT 640 SPECT/CT system. Details of this SPECT system geometry are given in section 3.3. The geometric sensitivity of the detector to different paths, the finite extent and bore diameters of the collimator and the finite energy and spatial resolution of the detectors were all modeled. For generating the LM data, we simulated the photon transport via a Monte Carlo-basedapproach in MATLAB. The scattering was modeled using the Klein–Nishina formula, which was normalized for a 2D system so that all the scattering was in plane. While simulating the photon transport, we considered photons that scattered at most once or twice, depending on the experimental setup. This was done to reduce the computational requirements in the FIM computation code, as we discuss in the next section. For each detected photon, the 1D estimated position of interaction of the incident gamma-ray photon with the scintillator, the estimated energy of the gamma-ray photon and the angular orientation of the detector were recorded in LM format.

For computing the FIM using the proposed path-based approach, software was developed in C programming language that used GPU acceleration using CUDA parallel computing platform. The first step was to implement the path-based formalism to describe the radiation transfer through different paths, as described in section 2.2. This formalism was validated by comparing it with the results obtained using the Monte Carlo approach described above (section 3.1) through several studies. The results obtained with the two approaches, for example, the photon flux, the energy spectrum and also the profiles of the projection data, were found to match. We do not show these results since our focus in this paper is on studying the information content of LM data.

The validated path-based formalism was then used to develop software to compute the FIM terms, as described by equations (50)–(53). For each detected LM event, we considered all the possible unscattered and scattered paths that the photon could have taken, regardless of the energy of the photon. In other words, photons were not labeled a priori as scattered or unscattered, since even a photon that falls in the photopeak (PP) window may have scattered. As is convention, PP window was characterized as having a width equal to twice the full-width half maximum (FWHM) centered at the photopeak energy. Further, in all our experiments, the energy window between 90 keV and the lower threshold of the PP window was considered as the scatter window. The FIM elements were computed at the true value of the parameters. These FIM terms were used to compute the CRB, which was then used to compute the lower bound on the standard deviations for the activity and attenuation coefficients of the different voxels.

A major challenge in the FIM computation was the large computational and memory requirements. These computational requirements were addressed using various algorithmic strategies. For example, to reduce the memory requirements, certain quantities that did not require large computation times were pre-computed and stored. This included quantitiessuch as the radiological path between the different voxels, distance covered by a sub-path inside a voxel and sensitivity of the collimator as a function of the angular and spatial voxel index. Further the code was parallelized and executed on graphics processing units (GPUs) for faster execution. More specifically, quantities that require sum over paths $\mathbb{P}$, such as $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)$ and $\mathrm{p}\mathrm{r}\left({\hat{\mathbf{A}}}_{j}\vert \mathbb{P}\right)\lambda \left(\mathbb{P}\right){s}_{\mathrm{e}\mathrm{ff}}\left(\mathbb{P}\right)\left[-{{\Delta}}_{q}\left(\mathbb{P}\right)+\frac{{\zeta }_{q}\left(\mathbb{P}\right)}{{\mu }_{q}}\right]$, were computed in parallel and summed using reduction algorithm. A more detailed procedure to compute the FIM terms, including a pseudocode of the GPU-based implementation, are presented in appendix B.

However, even after these computational optimizations, we observed that the code took substantial time to execute. A major reason for this was that the code scales as the number of paths, P, which in turn scales as the square of the number of voxels and as the number of angular samples used to define a path. Further, if we consider photons that have scattered M times, the number of path scales as P^M. To reduce this computational expenses and conduct experiments within a practical time limit, we considered phantoms with 32 × 32 pixels. Further,for most of the experiments, we modeled photons that scattered at most once; although in a subset of experiments, photons that scattered twice were modeled. In particular, in phantoms that had activity only within a single pixel, we could computationally model the dual-scattered photons in the FIM computation code within a practical time limit. Thus, in summary, the code was implemented to model single and dual-scattered events.

In the experiments, our objectives were to use the FIM computation framework to study the CRB of SPECT LM data for joint estimation of activity and attenuation distribution in simulation studies. The emission source was assumed to be Tc-99m, one of the most commonly used SPECT tracers, emitting photons at 140 keV. Photons were acquired at 64 angular positions spaced uniformly over 360°. A low-energy high-resolution parallel-hole collimator, with specifications similar to GE Optima NM/CT 640 SPECT/CT scanner was simulated. The system yielded a resolution of 7.8 mm at 10 cm depth. The scintillation detector had an intrinsic resolution of 4 mm and a length of 35 mm. Further, the energy resolution of the scintillation detector was set to 10% at 140 keV.

The system was used to first image a set of synthetic phantoms. For the synthetic phantoms, circular-shaped field of view with a diameter of 35 cm and 32 × 32 spatial pixels were considered. Two attenuation map configurations were considered: a uniform attenuation map with a constant attenuation value of 0.15 cm⁻¹ inside the field of view (figure 3(a)) and a non-uniform attenuation map (figure 3(b)) that simulated the cardiac region with an attenuation coefficient of 0.151 cm⁻¹ and 0.056 cm⁻¹ in background and lung region, respectively.Three activity distributions were considered, namely phantom with radiotracer uptake confined to a single pixel (figure 4(a)), uptake in multiple isolated pixels (figure 4(d)) and uptake over a donut-shaped region (figure 4(g)), referred to as single-pixel, multi-pixel and donut phantom, respectively. Experiments with these synthetic phantoms were conducted to gain an understanding of how activity distribution in uniform and non-uniform attenuation maps affect the information content of scattered photons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two anthropomorphic phantoms, namely the cardiac and the brain phantom, were also imaged to provide more clinical realism to our studies. The cardiac activity (figure 5(a)) and attenuation (figure 5(b)) phantoms were generated using a 2D slice of extended cardiac and torso (XCAT) phantom [42]. The brain activity (figure 5(d)) and attenuation (figure 5(e)) phantoms were generated using a 2D slice of the Zubal phantom [43]. For the study with thecardiac phantom, the SPECT system parameters were similar to as described above. For the study with the brain phantom, we simulated a DaTScan SPECT study [44] by modeling a system that was imaging ioflupane (I-123) tracer emitting photons with energy of 159 keV and had a circular field-of-view with a diameter of 30 cm.

Zoom In Zoom Out Reset image size

LM data for these phantoms were generated using the simulated SPECT system. In the first set of experiments, photons that scattered more than once were discarded.The FIM for the LM data was computed and used to determine the CRB for the activity and attenuation estimates in the different pixels. The first experiment was conducted with 4 × 10⁵ detected LM events, which is a clinically realistic count level in cardiac SPECT studies for a 2D slice [45]. Next, the detected photon counts were reduced from 4 × 10⁵ to 40 000 to study the CRB at different count levels. We also computed the CRB of the activity and attenuation coefficients when events only within the PP window were considered. This study was conducted to assess the impact of including the scattered photons on the CRB of these coefficients. Next, we studied the effect of varying the energy resolution of the system.

In the second set of experiments, we considered photons that scattered at most twice. These experiments were conducted for the single-pixel activity phantom with non-uniformattenuation (figure 3(b)). This dual scattering was modeled in the FIM code. Our objective was to assess the impact of including dual-scattered photons on estimating the attenuation coefficient. In the first experiment, we varied the number of detected LM events and computed the CRB for attenuation coefficients. Our second experiment investigated the impact of change in attenuation coefficient on the CRB. The rationale behind this experiment was that an increase in the true value of the attenuation coefficient causes an increase in the number ofscattered photons, but also increases the proportion of dual-scattered photons. An important question is whether even in these scenarios, the scattered photons provide information to estimate the attenuation coefficients. In this experiment, we varied the attenuation coefficient in the lung region of the non-uniform phantom map (figure 3(b)).