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

Consider a SPECT imaging 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 $\boldsymbol{\lambda}$ and $\boldsymbol{\mu}$, respectively. Denote the $N$ elements of the activity and attenuation distribution vector, by $\{\lambda_{1},\ldots,\lambda_{N}\}$ and $\{\mu_{1},\ldots,\mu_{N}\}$, respectively, where $\lambda_{i}$ is the mean number of photons emitted from the $i^{\mathrm{th}}$ voxel per unit time, and $\mu_{i}$ is the attenuation coefficient of the $i^{\mathrm{th}}$ voxel. When a transmission scan is used to measure the attenuation distribution $\boldsymbol{\mu}$, then the inverse problem is to estimate $\boldsymbol{\lambda}$. However, in this paper, we are considering that a transmission scan is unavailable, so the inverse problem is jointly estimating $\boldsymbol{\lambda}$ and $\boldsymbol{\mu}$ 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 $j^{\mathrm{th}}$ detected event by the attribute vectors $\mathbf{A}_{j}$ and $\mathbf{\hat{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}}=\{\mathbf{\hat{A}}_{j},j=1,2,...,J\}$.

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

$\displaystyle g({\mathbf{A}})=\sum_{j=1}^{J}\delta({\mathbf{A}}-\mathbf{\hat{A}}_{j}).$

(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}(\mathbf{A})$ for $m=1,...M$. Here we assume that the binning functions represent non-overlapping hyperrectangle and encompass the full range of attribute space. In that case, the measurement in the $m^{\mathrm{th}}$ bin, denoted by $g_{m}$ is given by [36]

$\displaystyle g_{m}=\int g(\mathbf{A})b_{m}(\mathbf{A})d^{q}\mathbf{A}=\sum_{j=1}^{J}b_{m}(\mathbf{\hat{A}}_{j}),$

(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 previous studies [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 $\boldsymbol{\lambda}$ and $\boldsymbol{\mu}$. 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]

$$\mathrm{pr}(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})=\mathrm{pr}(T|\boldsymbol{\lambda},\boldsymbol{\mu})\prod_{j=1}^{J}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\boldsymbol{\lambda},\boldsymbol{\mu}),$$

(3)

where $T$ is the acquisition time required for detecting $J$ LM events. Taking the logarithm on both sides of the Eq. (3) yields the log-likelihood of the observed LM data, denoted by $\mathcal{L}(\boldsymbol{\lambda},\boldsymbol{\mu}|\hat{\mathcal{A}},T)$ and given by

$$\mathcal{L}(\boldsymbol{\lambda},\boldsymbol{\mu}|\hat{\mathcal{A}},T)=\sum_{j=1}^{J}\log\mathrm{pr}(\mathbf{\hat{A}}_{j}|\boldsymbol{\lambda},\boldsymbol{\mu})+\log\mathrm{pr}(T|\boldsymbol{\lambda},\boldsymbol{\mu}).$$

(4)

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

$$\mathrm{pr}(T|\boldsymbol{\lambda},\boldsymbol{\mu})=\frac{\beta^{J}T^{J-1}\exp(-\beta T)}{(J-1)!}.$$

(5)

Obtaining a similar direct analytic expression for $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\boldsymbol{\lambda},\boldsymbol{\mu})$ is complicated. To address this issue, note that any LM event is the result of a photon emitted from a voxel, travelling in a certain direction, and then, in some cases, scatterring 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{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})$, and the probability mass function of the path, $\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})$, can be derived, as described later. We thus decompose $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\boldsymbol{\lambda},\boldsymbol{\mu})$ as a mixture model over all possible paths. For this purpose, we use the following identity:

$$\mathrm{pr}(x|y=Y)=\sum_{z}\mathrm{pr}(x|y=Y,z=Z)\mathrm{Pr}(z=Z|y=Y),$$

(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{pr}(\mathbf{\hat{A}}_{j}|\boldsymbol{\lambda},\boldsymbol{\mu})$:

$$\mathrm{pr}(\mathbf{\hat{A}}_{j}|\boldsymbol{\lambda},\boldsymbol{\mu})=\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P},\boldsymbol{\lambda},\boldsymbol{\mu})\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu}).$$

(7)

The components of the mixture model are the probabilities that a 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 the path is independent of the activity and attenuation distribution, i.e. $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P},\boldsymbol{\lambda},\boldsymbol{\mu})=\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})$, provided the probability of the path accounts for the emission and attenuation processes, as will be the case in our treatment. Using Eq. (4), we can rewrite the log-likelihood of the data given the activity and attenuation distribution in terms of this mixture-model decomposition as

$$\mathcal{L}(\boldsymbol{\lambda},\boldsymbol{\mu}|\hat{\mathcal{A}},T)=\sum_{j=1}^{J}\log\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})+\log\mathrm{pr}(T|\boldsymbol{\lambda},\boldsymbol{\mu}).$$

(8)

To derive the elements of the FIM, analytic expressions for $\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})$ and $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})$ must be derived. These are the topics of the next two sub-sections.

In this section, we derive the expression for $\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})$. 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\pi$ radians into a finite number of solid-angle sub-domains, referred to as ordinates, each of size $\Delta\Omega$ . Denote the 3-D unitary direction vector by $\hat{\boldsymbol{s}}$ and discretized angular direction by $\hat{\boldsymbol{s}}_{k}$ associated with the $k^{\mathrm{th}}$ instance of the ordinate. Then the indicator function of the $k^{\mathrm{th}}$ angular ordinate $\psi_{k}(\hat{\boldsymbol{s}})$ is defined as below:

$$\psi_{k}(\hat{\boldsymbol{s}})=\begin{cases}1,\quad\text{if $k^{\mathrm{th}}$ ordinate contains the direction vector $\hat{\boldsymbol{s}}$.}\\ 0,\quad\text{otherwise.}\end{cases}$$

(9)

A sub-path $\mathbb{S}_{ik}$ is now defined as a cone with a vertex at the center of the $i^{\mathrm{th}}$ voxel and an angular range given by $\psi_{k}(\hat{\boldsymbol{s}})$. In our analysis, we consider a path between voxels with indices $i_{0}$, $i_{1}$, …$i_{n}$. This path can be considered to consist of $n$ subpaths between $i_{0}$ to $i_{1}$, $i_{1}$ to $i_{2}$ and so on until $i_{n-1}$ to $i_{n}$. A schematic describing the above notations and other notations is presented in Fig. 2.

The probability of a particular path $\mathbb{P}$, i.e. $\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})$, 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(\mathbb{P})$ denote the flux of photons incident on the detector through a path $\mathbb{P}$. Then,

$$\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})=\frac{\Phi(\mathbb{P})}{\sum_{\mathbb{P}^{\prime}}\Phi(\mathbb{P}^{\prime})}.$$

(10)

To derive the expression for $\Phi(\mathbb{P})$, we first define a few other radiometric quantities. Let $\mathcal{Q}$ denote the radiant energy, the 3D vector $\boldsymbol{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(\boldsymbol{r},\hat{\boldsymbol{s}},E)$, given by

$$w(\boldsymbol{r},\hat{\boldsymbol{s}},E)=\frac{1}{E}\frac{\partial^{3}\mathcal{Q}}{\partial V\partial\Omega\partial E}.$$

(11)

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

$$\Xi(\boldsymbol{r},\hat{\boldsymbol{s}},E)=\frac{\partial^{3}\Phi}{\partial V\partial\Omega\partial E}.$$

(12)

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

$$\Gamma(\hat{\boldsymbol{s}})=c_{m}\int\int w(\boldsymbol{r},\hat{\boldsymbol{s}},E)d^{3}\boldsymbol{r}dE.$$

(13)

Consider a point source at location $\boldsymbol{r}_{\mathrm{s}}$ in the voxel indexed by $i_{0}$ emitting photons of energy $E_{0}$ 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

$$\Xi(\boldsymbol{r},\hat{\boldsymbol{s}},E)=\frac{\lambda_{i_{0}}}{4\pi}\delta(\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}})\delta(E-E_{0})\psi_{k_{0}}(\hat{\boldsymbol{s}}),$$

(14)

where $\delta(x)$ is the Dirac delta function. As photons travel from the voxel $i_{0}$ to $i_{1}$, 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]

$$[\mathcal{X}w](\boldsymbol{r},\hat{\boldsymbol{s}},E)=\frac{1}{c_{m}}\int_{0}^{\infty}w(\boldsymbol{r}-\hat{\boldsymbol{s}}l)\exp\left[-\int_{0}^{l}\mu(\boldsymbol{r}-\hat{\boldsymbol{s}}l^{\prime},E)dl^{\prime}\right]dl,$$

(15)

where $\mu(\boldsymbol{r},E)$ denotes the attenuation coefficient at location $\boldsymbol{r}$ and energy $E$, and $c_{m}$ denotes the speed of light. Applying the attenuation transform to the source distribution function (Eq. (14)) yields

$\displaystyle[\mathcal{X}\Xi](\boldsymbol{r},\hat{\boldsymbol{s}},E)=\frac{\lambda_{i_{0}}}{4\pi c_{\mathrm{m}}}\delta(E-E_{0})\psi_{k_{0}}(\hat{\boldsymbol{s}})\int_{0}^{\infty}\delta(\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}-\hat{\boldsymbol{s}}l)\exp\left[-\int_{0}^{l}\mu(\boldsymbol{r}-\hat{\boldsymbol{s}}l^{\prime},E)dl^{\prime}\right]dl.$

(16)

Now, it can be shown that [35]

$$\int_{l=0}^{\infty}\delta(\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}-\hat{\boldsymbol{s}}l)dl=\frac{1}{|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|^{2}}\delta\left(\hat{\boldsymbol{s}}-\hat{\boldsymbol{s}}_{10}\right),$$

(17)

where

$\displaystyle\hat{\boldsymbol{s}}_{10}=\frac{\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}}{|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|}.$

(18)

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

	$\displaystyle[\mathcal{X}\Xi](\boldsymbol{r},\hat{\boldsymbol{s}},E)$	$\displaystyle=\frac{\lambda_{i_{0}}}{4\pi c_{\mathrm{m}}|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|^{2}}\exp\left[-\int_{0}^{|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|}\mu(\boldsymbol{r}-\hat{\boldsymbol{s}}l^{\prime},E)dl^{\prime}\right]\times$
		$\displaystyle\delta\left(\hat{\boldsymbol{s}}-\hat{\boldsymbol{s}}_{10}\right)\delta(E-E_{0})\psi_{k_{0}}(\hat{\boldsymbol{s}}).$		(19)

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

$$[\mathcal{K}w](\boldsymbol{r},\hat{\boldsymbol{s}},E)=\int\int_{4\pi}K(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}^{\prime},E,E^{\prime}|\boldsymbol{r})w(\boldsymbol{r},\hat{\boldsymbol{s}}^{\prime},E^{\prime})d\Omega^{\prime}dE^{\prime},$$

(20)

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

$\displaystyle K(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}^{\prime},E,E^{\prime}|\boldsymbol{r})=c_{\mathrm{m}}n_{\mathrm{sc}}(\boldsymbol{r})\frac{\partial\sigma_{sc}}{\partial\Omega}\delta\left\{E-\left[\frac{1}{E^{\prime}}+\frac{1}{mc^{2}}(1-\cos\theta)\right]^{-1}\right\},$

(21)

where $n_{\mathrm{sc}}(\boldsymbol{r})$ is the density of scatterers and is related to the scattering coefficient by

$$\mu(\boldsymbol{r})=n_{\mathrm{sc}}(\boldsymbol{r})\sigma_{\mathrm{sc}},$$

(22)

where $\sigma_{\mathrm{sc}}$ is the scattering cross section. Also, $\theta$ is the angle at which the outgoing photon scatters relative to the incoming photon, so $\cos\theta=\hat{\boldsymbol{s}}\cdot\hat{\boldsymbol{s}}^{\prime}$. Finally, the differential scattering cross section $\frac{\partial\sigma_{sc}}{\partial\Omega}$ in Eq. (21) is given by the Klein-Nishina formula [17]. For notational simplicity, define $K_{\mathrm{mag}}(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}^{\prime},E|\boldsymbol{r})$ by

$$K_{\mathrm{mag}}(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}^{\prime},E|\boldsymbol{r})=n_{\mathrm{sc}}(\boldsymbol{r})\frac{\partial\sigma_{sc}}{\partial\Omega}\bigg{|}_{\cos\theta=\hat{\boldsymbol{s}}\cdot\hat{\boldsymbol{s}}^{\prime}}.$$

(23)

Also, denote the path integral between any two locations $\boldsymbol{r}_{l}$ and $\boldsymbol{r}_{m}$ by the function $\gamma(\boldsymbol{r}_{l},\boldsymbol{r}_{m},E)$, i.e.

$$\gamma(\boldsymbol{r}_{l},\boldsymbol{r}_{m},E)=\int_{0}^{|\boldsymbol{r}_{l}-\boldsymbol{r}_{m}|}\mu\left(\boldsymbol{r}_{l}-t\frac{\boldsymbol{r}_{1}-\boldsymbol{r}_{m}}{|\boldsymbol{r}_{1}-\boldsymbol{r}_{m}|},E\right)\mathrm{d}t.$$

(24)

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

$\displaystyle[\mathcal{K}\mathcal{X}\Xi](\boldsymbol{r},\hat{\boldsymbol{s}},E)$

$\displaystyle=\frac{\lambda_{i_{0}}}{4\pi c_{\mathrm{m}}|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|^{2}}K\left(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}_{10},E,E_{0}|\boldsymbol{r}\right)\exp\left[-\gamma(\boldsymbol{r},\boldsymbol{r}_{\mathrm{s}},E_{0})\right]\psi_{k_{0}}(\hat{\boldsymbol{s}}_{10}).$

(25)

We now integrate this distribution function over all possible locations within the $i_{1}^{\mathrm{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}}$ subpath and scattering within the $i_{1}^{\mathrm{th}}$ voxel. Denote this radiant intensity by $\Gamma_{i_{0}i_{1}k_{0}}(\hat{\boldsymbol{s}})$. Substituting the distribution function from Eq. (25) into Eq. (13) yields

$\displaystyle\Gamma_{i_{0}i_{1}k_{0}}(\hat{\boldsymbol{s}})=\int\frac{\lambda_{i_{0}}}{4\pi|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|^{2}}K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}_{10},E_{0}|\boldsymbol{r}\right)\exp\left[-\gamma(\boldsymbol{r},\boldsymbol{r}_{\mathrm{s}},E_{0})\right]\phi_{i_{1}}(\boldsymbol{r})\psi_{k_{0}}(\hat{\boldsymbol{s}}_{10})d^{3}\boldsymbol{r},$

(26)

where $\phi_{i}(\boldsymbol{r})$ is the voxel basis function defined as below:

$$\phi_{i}(\boldsymbol{r})=\begin{cases}1,\quad\boldsymbol{r}\ \text{lies within voxel i.}\\ 0,\quad\mathrm{otherwise.}\end{cases}$$

(27)

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

$\displaystyle\Gamma_{i_{0}i_{1}k_{0}}(\hat{\boldsymbol{s}})\approx\frac{\lambda_{i_{0}}}{4\pi}\exp[-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})]K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}_{c10},E_{0}|\boldsymbol{r}_{1}\right)\int\frac{1}{|\boldsymbol{r}-\boldsymbol{r}_{\mathrm{s}}|^{2}}\psi_{k_{0}}(\hat{\boldsymbol{s}}_{10})\phi_{i_{1}}(\boldsymbol{r})d^{3}\boldsymbol{r}.$

(28)

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

	$\displaystyle\Gamma_{i_{0}i_{1}k_{0}}(\hat{\boldsymbol{s}})\approx$	$\displaystyle\frac{\lambda_{i_{0}}}{4\pi}\exp[-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})]\times$
		$\displaystyle K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}_{c10},E_{0}|\boldsymbol{r}_{1}\right)\int\psi_{k_{0}}(\hat{\boldsymbol{s}}_{10})\int\phi_{i_{1}}(\boldsymbol{r}_{s}+R\hat{\boldsymbol{s}}_{10})dR~d\hat{\boldsymbol{s}}_{10}.$		(29)

The integral over $R$ is equal to the distance traversed by the $\boldsymbol{r}_{s}+R\hat{\boldsymbol{s}}_{10}$ vector within the $k_{1}^{\mathrm{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}^{\mathrm{th}}$ voxel, which we denote by $\Delta_{i_{1}}(\mathbb{S}_{i_{0}k_{0}})$. Note that $\hat{\boldsymbol{s}}_{k_{0}}$ is the discretized direction vector associated with the $k_{0}^{\mathrm{th}}$ ordinate. Further, performing the integral over $\hat{\boldsymbol{s}}_{10}$ yields the following expression for $\Gamma_{i_{0}i_{1}k_{0}}(\hat{\boldsymbol{s}})$:

$\displaystyle\Gamma_{i_{0}i_{1}k_{0}}(\hat{\boldsymbol{s}})\approx\frac{\lambda_{i_{0}}}{4\pi}\exp[-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})]K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}_{c10},E_{0}|\boldsymbol{r}_{1}\right)\Delta_{i_{1}}(\mathbb{S}_{i_{0}k_{0}})\Delta\Omega.$

(30)

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

$\displaystyle w_{1}(\boldsymbol{r},\hat{\boldsymbol{s}},E)=\frac{\lambda_{i_{0}}}{4\pi c_{\mathrm{m}}}\exp\left[-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})\right]K\left(\hat{\boldsymbol{s}},\hat{\boldsymbol{s}}_{c10},E,E_{0}|\boldsymbol{r}_{1}\right)\Delta_{i_{1}}(\mathbb{S}_{i_{0}k_{0}})\delta(\boldsymbol{r}-\boldsymbol{r}_{1})\psi_{k_{1}}(\hat{\boldsymbol{s}}).$

(31)

After this scattering operation, the photons along this path suffer from attenuation as they travel towards voxel $i_{2}$.

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 subpath 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 $\boldsymbol{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}(\boldsymbol{r}_{d},\hat{\boldsymbol{s}},E)$ is given by

	$\displaystyle w_{d}(\boldsymbol{r}_{d},\hat{\boldsymbol{s}},E)$	$\displaystyle=\frac{\prod_{m=0}^{n-1}\Delta_{i_{m+1}}(\mathbb{S}_{i_{m}k_{m}})}{|\boldsymbol{r}_{d}-\boldsymbol{r}_{n}|^{2}}\exp\{-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})-\ldots\gamma(\boldsymbol{r}_{d},\boldsymbol{r}_{n},E_{n})\}\times$
		$\displaystyle\prod_{m=0}^{n-1}\{K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}}_{c,m+2,m+1},\hat{\boldsymbol{s}}_{c,m+1,m},E_{m}|\boldsymbol{r}_{m+1}\right)\}\delta(E-E_{n})\psi_{k_{n}}(\hat{\boldsymbol{s}}).$		(32)

Now the plane of the front face of the detector is given by $\delta(p-\boldsymbol{r}_{d}\cdot\hat{\boldsymbol{n}})$, 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 $\boldsymbol{r}$ in direction $\hat{\boldsymbol{s}}$ be denoted by $t(\boldsymbol{r},\hat{\boldsymbol{s}})$. Then, the flux detected through the considered path is given by

	$\displaystyle\Phi(\mathbb{P})=c_{\mathrm{m}}\int\int\int t(\boldsymbol{r}_{n},\hat{\boldsymbol{s}})(\hat{\boldsymbol{n}}\cdot\hat{\boldsymbol{s}})\delta(p-\boldsymbol{r}_{d}\cdot\hat{\boldsymbol{n}})w_{d}(\boldsymbol{r}_{d},\hat{\boldsymbol{s}},E)~d^{3}\boldsymbol{r}_{d}~d\Omega~dE$
	$\displaystyle=\frac{\lambda_{i_{0}}}{4\pi}\int\int_{2\pi}\int t(\boldsymbol{r}_{n},\hat{\boldsymbol{s}})(\hat{\boldsymbol{n}}\cdot\hat{\boldsymbol{s}})\delta(p-\boldsymbol{r}_{d}\cdot\hat{\boldsymbol{n}})\times$
	$\displaystyle\exp\{-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})-\ldots\gamma(\boldsymbol{r}_{d},\boldsymbol{r}_{n},E_{n})\}\frac{\prod_{m=0}^{n-1}\Delta_{i_{m+1}}(\mathbb{S}_{i_{m}k_{m}})}{|\boldsymbol{r}_{d}-\boldsymbol{r}_{n}|^{2}}\times$
	$\displaystyle\prod_{m=0}^{n-1}\{K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}}_{c,m+2,m+1},\hat{\boldsymbol{s}}_{c,m+1,m},E_{m}|\boldsymbol{r}_{m+1}\right)\}\delta(E-E_{n})\psi_{k_{n}}(\hat{\boldsymbol{s}})~d^{3}\boldsymbol{r}_{d}~d\Omega~dE.$		(33)

Perform a change of variables by replacing $\boldsymbol{r}_{d}-\boldsymbol{r}_{n}$ by $R\hat{\boldsymbol{s}}$, so that $d^{3}r_{d}=R^{2}dRd\hat{\boldsymbol{s}}$. Next, integrating over $E$ and $\Omega$ using the sifting property of the delta function yields

	$\displaystyle\Phi(\mathbb{P})=\frac{\lambda_{i_{0}}}{4\pi}\int\int\delta(p-\boldsymbol{r}_{n}\cdot\hat{\boldsymbol{n}}-R\hat{\boldsymbol{s}}\cdot\hat{\boldsymbol{n}})t\left(\boldsymbol{r}_{n},\hat{\boldsymbol{s}}\right)(\hat{\boldsymbol{n}}\cdot\hat{\boldsymbol{s}})\times$
	$\displaystyle\exp\{-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})-\ldots\gamma(\boldsymbol{r}_{n}+R\hat{\boldsymbol{s}},\boldsymbol{r}_{n},E_{n})\}\prod_{m=0}^{n-1}\{\Delta_{i_{m+1}}(\mathbb{S}_{i_{m}k_{m}})\}\times$
	$\displaystyle\prod_{m=0}^{n-1}\{K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}}_{c,m+2,m+1},\hat{\boldsymbol{s}}_{c,m+1,m},E_{m}|\boldsymbol{r}_{m+1}\right)\}\psi_{k_{n}}(\hat{\boldsymbol{s}})~d\hat{\boldsymbol{s}}~dR.$		(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 a 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(\boldsymbol{r},\hat{\boldsymbol{s}})>0$, the integral over $R$ can be performed to yield

	$\displaystyle\Phi(\mathbb{P})=\frac{\lambda_{i_{0}}}{4\pi}\int t\left(\boldsymbol{r}_{n},\hat{\boldsymbol{s}}\right)\exp\{-\gamma(\boldsymbol{r}_{1},\boldsymbol{r}_{0},E_{0})-\ldots\gamma(\boldsymbol{r}_{d_{0}},\boldsymbol{r}_{n},E_{n})\}\times$
	$\displaystyle\prod_{m=0}^{n-1}\{\Delta_{i_{m+1}}(\mathbb{S}_{i_{m}k_{m}})\}\prod_{m=0}^{n-1}\{K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}}_{c,m+2,m+1},\hat{\boldsymbol{s}}_{c,m+1,m},E_{m}|\boldsymbol{r}_{m+1}\right)\}\psi_{k_{n}}(\hat{\boldsymbol{s}})~d\hat{\boldsymbol{s}},$		(35)

where $\boldsymbol{r}_{d0}=\boldsymbol{r}_{n}+(p-\boldsymbol{r}_{n}\cdot\hat{\boldsymbol{n}})\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}^{\mathrm{th}}$ ordinate.. Then the integration over $\hat{\boldsymbol{s}}$ can be performed to yield

	$\displaystyle\Phi(\mathbb{P})=\frac{\lambda_{i_{0}}\Delta\Omega}{4\pi}\exp\{-\gamma(\boldsymbol{r}_{0},\boldsymbol{r}_{1},E_{0})-\ldots\gamma(\boldsymbol{r}_{dc0},\boldsymbol{r}_{n},E_{n})\}t\left(\boldsymbol{r}_{n},\hat{\boldsymbol{s}}_{k_{n}}\right)$
	$\displaystyle\prod_{m=0}^{n-1}\Delta_{i_{m+1}}(\mathbb{S}_{i_{m}k_{m}})\prod_{m=0}^{n-1}K_{\mathrm{mag}}\left(\hat{\boldsymbol{s}}_{c,m+2,m+1},\hat{\boldsymbol{s}}_{c,m+1,m},E_{m}|\boldsymbol{r}_{m+1}\right),$		(36)

where $\boldsymbol{r}_{dc0}=\boldsymbol{r}_{n}+(p-\boldsymbol{r}_{n}\cdot\hat{\boldsymbol{n}})\hat{\boldsymbol{s}}_{k_{n}}$. To simplify the expression for $\Phi(\mathbb{P})$, 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(\mathbb{P})$. Further, to simplify notation, we define $s_{\mathrm{eff}}(\mathbb{P})$ as

$$s_{\mathrm{eff}}(\mathbb{P})=\Lambda(\mathbb{P})\exp\left[-\sum_{m=0}^{n}\gamma(\boldsymbol{r}_{m},\boldsymbol{r}_{m+1},\boldsymbol{\mu}(E_{m}))\right]\prod_{m=1}^{n}\mu_{k_{m}}(E_{m-1}).$$

(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(\mathbb{P})$. Eq. (36) can then be rewritten as

$$\Phi(\mathbb{P})=s_{\mathrm{eff}}(\mathbb{P})\lambda(\mathbb{P}).$$

(38)

The attenuation coefficient in Eq. (37) is a function of energy, which does not lend itself to Fisher information matrix analysis. 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 coefficent [40] as a linear function of energy, i.e.

$$\boldsymbol{\mu}(E_{m})=a\boldsymbol{\mu}+b.$$

(39)

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

$$\mathrm{Pr}(\mathbb{P}|\boldsymbol{\lambda},\boldsymbol{\mu})=\frac{\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}{\sum_{\mathbb{P^{\prime}}}\lambda(\mathbb{P^{\prime}})s_{\mathrm{eff}}(\mathbb{P^{\prime}})}.$$

(40)

This explicit modeling of probability of path in terms of $\boldsymbol{\lambda}$ and $\boldsymbol{\mu}$ is then used to express the log-likelihood and Fisher information analysis.

The term $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})$ denotes the PDF of the measured attribute vector $\mathbf{\hat{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 energy deposited by the photon in the crystal, denoted by $\hat{E}_{j}$. To model the randomness in estimating $\mathbf{\hat{A}}_{j}$ due to the uncertainty in the estimation process and the finite energy and spatial resolution of the detectors, we first write $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})$ as

$$\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})=\int\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbf{A}_{j})\mathrm{pr}(\mathbf{A}_{j}|\mathbb{P})d\mathbf{A}_{j}.$$

(41)

To obtain the expression for $\mathrm{pr}(\mathbf{A}_{j}|\mathbb{P})$, 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

$$\mathrm{pr}(\mathbf{A}_{j}|\mathbb{P})=\mathrm{pr}(\boldsymbol{r},E|\mathbb{P})\approx\delta(\mathbf{A}_{j}-\mathbf{A}_{j}({\mathbb{P}})),$$

(42)

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

Next, to derive the expression for the term $\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbf{A}_{j})$, 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 [41, 31]. Then, we can write

$$\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbf{A}_{j})=\frac{1}{(2\pi)^{2}|\boldsymbol{K}_{\mathrm{det}}|}\exp\left[-\frac{(\mathbf{\hat{A}}_{j}-\mathbf{A}_{j})^{\dagger}\boldsymbol{K}_{\mathrm{det}}^{-1}(\mathbf{\hat{A}}_{j}-\mathbf{A}_{j})}{2}\right],$$

(43)

where $\boldsymbol{K}_{\mathrm{det}}$ denotes the covariance matrix quantifying the variances and co-variances in the energy and position estimates, and where $|\boldsymbol{K}_{\mathrm{det}}|$ denotes the determinant of the matrix $\boldsymbol{K}_{\mathrm{det}}$. Substituting Eq. (43) and (42) into Eq. (41), and using the sifting property of the delta function yields

$$\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})=\frac{1}{(2\pi)^{2}|\boldsymbol{K}_{\mathrm{det}}|}\exp\left[-\frac{(\mathbf{\hat{A}}_{j}-\mathbf{A}_{j}(\mathbb{P}))^{\dagger}\boldsymbol{K}_{\mathrm{det}}^{-1}(\mathbf{\hat{A}}_{j}-\mathbf{A}_{j}(\mathbb{P}))}{2}\right].$$

(44)

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

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

$$F_{qq^{\prime}}=-\left\langle\frac{\partial^{2}\mathcal{L}(\boldsymbol{\lambda},\boldsymbol{\mu}|\hat{\mathcal{A}},T)}{\partial\theta_{q}\partial\theta_{q^{\prime}}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})},$$

(45)

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

$\displaystyle\mathcal{L}(\boldsymbol{\lambda},\boldsymbol{\mu}|\hat{\mathcal{A}},T)=\sum_{j=1}^{J}\log\left[\sum\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\right]+(J-1)\log T-\beta T-\log(J-1)!.$

(46)

Now, note that $\beta$, which is the mean rate of photons detected, is equivalently the total radiant flux over all possible paths. Thus, $\beta$ is given by

$$\beta=\sum_{\mathbb{P^{\prime}}}\lambda(\mathbb{P^{\prime}})s_{\mathrm{eff}}(\mathbb{P^{\prime}}).$$

(47)

Using Eq. (47) to substitute for $\beta$, and differentiating the log-likelihood with respect to the activity in the $q^{\mathrm{th}}$ voxel $\lambda_{q}$ yields

$$\frac{\partial\mathcal{L}}{\partial\lambda_{q}}=\sum_{j=1}^{J}\frac{\sum_{\mathbb{P}_{q}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}-T{\sum_{\mathbb{P}_{q}}s_{\mathrm{eff}}(\mathbb{P})},$$

(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 (Eq. (46) with respect to the attenuation coefficient in the $q^{\mathrm{th}}$ voxel, i.e. $\mu_{q}$, yields

	$\displaystyle\frac{\partial\mathcal{L}}{\partial\mu_{q}}=$	$\displaystyle\sum_{j=1}^{J}\frac{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\left[-\Delta_{q}(\mathbb{P})+\frac{\zeta_{q}(\mathbb{P})}{\mu_{q}}\right]}{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}\times$
		$\displaystyle-T{\sum_{\mathbb{P}}\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}\left[-\Delta_{q}(\mathbb{P})+\frac{\zeta_{q}(\mathbb{P})}{\mu_{q}}\right],$		(49)

where $\zeta_{q}(\mathbb{P})$ and $\Delta_{q}(\mathbb{P})$ are the number of scatter events occurring in the $q^{\mathrm{th}}$ voxel in the considered path and the distance that the considered path covers in the $q^{\mathrm{th}}$ voxel, respectively. To derive the FIM elements, we must differentiate Eqs. (48) and (49) further with respect to the activity and attenuation coefficients in some other $q^{\prime}{\mathrm{th}}$ voxel, and then average over the observed LM data . The derivations are detailed in Appendix A, and the final expressions are as below:

	$\displaystyle\left\langle\frac{\partial^{2}\mathcal{L}}{\partial\mu_{q^{\prime}}\partial\mu_{q}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})}$	$\displaystyle=-\left\langle\sum_{j=1}^{J}\frac{\left\{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\left[-\Delta_{q}(\mathbb{P})+\frac{\zeta_{q}(\mathbb{P})}{\mu_{q}}\right]\right\}}{{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}}\times\right.$
		$\displaystyle\quad\left.\frac{\left\{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\left[-\Delta_{q^{\prime}}(\mathbb{P})+\frac{\zeta_{q^{\prime}}(\mathbb{P})}{\mu_{q^{\prime}}}\right]\right\}}{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})},$		(50)
	$\displaystyle\left\langle\frac{\partial^{2}\mathcal{L}}{\partial\lambda_{q^{\prime}}\partial\lambda_{q}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})}$	$\displaystyle=-\left\langle\sum_{j=1}^{J}\frac{\{\sum_{\mathbb{P}_{q}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\}\{\sum_{\mathbb{P}_{q^{\prime}}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\}}{\{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\}^{2}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})},$		(51)
	$\displaystyle\left\langle\frac{\partial^{2}\mathcal{L}}{\partial\mu_{q^{\prime}}\partial\lambda_{q}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})}$	$\displaystyle=-\left\langle\sum_{j=1}^{J}\frac{\left\{\sum_{\mathbb{P}_{q}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\right\}}{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}\times\right.$
		$\displaystyle\quad\left.\frac{\left\{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})\left[-\Delta_{q^{\prime}}(\mathbb{P})+\frac{\zeta_{q^{\prime}}(\mathbb{P})}{\mu_{q^{\prime}}}\right]\right\}}{\sum_{\mathbb{P}}\mathrm{pr}(\mathbf{\hat{A}}_{j}|\mathbb{P})\lambda(\mathbb{P})s_{\mathrm{eff}}(\mathbb{P})}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})},$		(52)
	$\displaystyle\left\langle\frac{\partial^{2}\mathcal{L}}{\partial\lambda_{q^{\prime}}\partial\mu_{q}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})}$	$\displaystyle=\left\langle\frac{\partial^{2}\mathcal{L}}{\partial\mu_{q}\partial\lambda_{q^{\prime}}}\right\rangle_{(\hat{\mathcal{A}},T|\boldsymbol{\lambda},\boldsymbol{\mu})}.$		(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 Cramer-Rao bound (CRB), which is the lower bound on the variance 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.