Tomographic reconstruction from Poisson distributed data: a fast and convergent EM-TV dual approach

Abstract

This paper focuses on tomographic reconstruction for nuclear medicine imaging, where a classical approach consists to maximize the likelihood of Poisson distributed data using the iterative Expectation Maximization algorithm. In this context and when the quantity of acquired data is low and produces low signal-to-noise ratio in the images, a step forward consists to incorporate a total variation prior on the solution into a MAP-EM formulation. This prior is not differentiable. The novelty of the paper is to propose a convergent and efficient numerical scheme to compute the MAP-EM optimizer, alternating regular maximum likelihood maximization steps and TV-denoising solved using the convex-duality principle of Fenchel-Rockafellar. The main theoretical result is the proof of stability and convergence of this scheme. We also present some numerical experiments where we compare the proposed algorithm with some other algorithms from the literature.

Appendices

Appendix A: Proof of Lemma 1

For a given $\varvec{\lambda }^{\prime }\in D$ and any $\varvec{\lambda }\in D$, we denote $K(\varvec{\lambda }\vert \varvec{\lambda }^{\prime })= F(\varvec{\lambda })-U(\varvec{\lambda }\vert \varvec{\lambda }^{\prime })$. Let us start with some results concerning $K(\varvec{\lambda }\vert \varvec{\lambda }^{\prime })$. We have:

$$\begin{aligned} K(\varvec{\lambda }\vert \varvec{\lambda }^{\prime })= \sum \limits _{i=1}^{I} \sum \limits _{j=1}^{J} t_{ij}\lambda _{j}- \sum \limits _{i=1}^{I} y_{i}\log \left( \sum \limits _{j=1}^{J} t_{ij}\lambda _{j}\right) - \sum \limits _{j=1}^{J} s_{j}\lambda _{j} +\sum \limits _{j=1}^{J} s_{j}\overline{\lambda ^{\prime }}_{j}\log (\lambda _{j}). \end{aligned}$$

Since $\sum \limits _{i=1}^{I} \sum \limits _{j=1}^{J} t_{ij}\lambda _{j} = \sum \limits _{j=1}^{J} s_{j}\lambda _{j}$, thus

$$\begin{aligned} K(\varvec{\lambda }\vert \varvec{\lambda }^{\prime })= \sum \limits _{j=1}^{J} s_{j}\overline{\lambda ^{\prime }}_{j}\log (\lambda _{j}) - \sum \limits _{i=1}^{I} y_{i}\log \left( \sum \limits _{j=1}^{J} t_{ij}\lambda _{j}\right) . \end{aligned}$$

(A1)

The partial derivatives of $K(\cdot \vert \varvec{\lambda }^{\prime })$ are:

$$\begin{aligned} \frac{\partial K(\cdot \vert \varvec{\lambda }^{\prime })}{\partial \lambda _{j}}(\varvec{\lambda }) = s_{j}\frac{\overline{\lambda ^{\prime }}_{j}}{\lambda _{j}}-\sum \limits _{i=1}^{I}\frac{t_{ij}y_{i}}{\sum \limits _{k=1}^{J} t_{ik}\lambda _{k}}. \end{aligned}$$

We recall that $\overline{\varvec{\lambda }^{\prime }}$ is the vector with entries $\overline{\lambda ^{\prime }}_{j} = \displaystyle \frac{\lambda ^{\prime }_{j}}{s_{j}}\sum _{i=1}^{I}\frac{t_{ij}y_{i}}{\sum _{k} t_{ik}\lambda ^{\prime }_{k}}$. This implies that for all $\varvec{\lambda }^{\prime }\in D$

$$\begin{aligned} \nabla K(\cdot \vert \varvec{\lambda }^{\prime })(\varvec{\lambda }^{\prime })=0. \end{aligned}$$

(A2)

The logarithm being a concave function, by the Jensen inequality we obtain for all $i\in \{1,\dots ,I\}$ that:

$$\begin{aligned} \log \left( \sum \limits _{j=1}^{J} t_{ij}\lambda _{j}\right)= & {} \log \left( \sum \limits _{j=1}^{J} \frac{t_{ij}\lambda _{j}^{\prime }}{\sum _{k} t_{ik}\lambda _{k}^{\prime }}\frac{\lambda _{j}}{\lambda ^{\prime }_{j}}\right) + \log \left( \sum _{k} t_{ik}\lambda _{k}^{\prime }\right) \\\ge & {} \sum \limits _{j=1}^{J} \frac{t_{ij}\lambda _{j}^{\prime }}{\sum _{k} t_{ik}\lambda _{k}^{\prime }} \log \left( \frac{\lambda _{j}}{\lambda ^{\prime }_{j}}\right) + \log \left( \sum _{k} t_{ik}\lambda _{k}^{\prime }\right) \end{aligned}$$

thus which implies that for all $\varvec{\lambda },\varvec{\lambda }^{\prime }$,

$$\begin{aligned} K(\varvec{\lambda }\vert \varvec{\lambda }^{\prime })\le K(\varvec{\lambda }^{\prime }\vert \varvec{\lambda }^{\prime }). \end{aligned}$$

(A3)

Proof of Lemma 1

Let $\varvec{\lambda }^{*}$ be a minimum of F. From (A3), for all $\varvec{\lambda }\in D$, $K(\varvec{\lambda }\vert \varvec{\lambda }^{*})\le K(\varvec{\lambda }^{*}\vert \varvec{\lambda }^{*})$ thus $U(\varvec{\lambda }\vert \varvec{\lambda }^{*})\ge U(\varvec{\lambda }^{*}\vert \varvec{\lambda }^{*})$ and (16) is verified. Conversely, if $\varvec{\lambda }^{*}$ verifies (16) then $0\in \partial U(\cdot \,\vert \varvec{\lambda }^{*})(\varvec{\lambda }^{*}) = \partial F(\varvec{\lambda }^{*})-\nabla K(\cdot \vert \varvec{\lambda }^{*})(\varvec{\lambda }^{*})$. From (A2), $\nabla K(\cdot \vert \varvec{\lambda }^{*})(\varvec{\lambda }^{*}) =0$, thus $0\in \partial F(\varvec{\lambda }^{*})$ and $\varvec{\lambda }^{*}$ is a minimum of F. $\square$

Appendix B: Proof of Theorem 2

Proof

• (i) From the definition of K and from (A3) it follows that

  $$ \begin{array}{ll} F(M(\varvec{\lambda }^{\prime })) &{}= U(M(\varvec{\lambda }^{\prime })\vert \varvec{\lambda }^{\prime }) + K(M(\varvec{\lambda }^{\prime })\vert \varvec{\lambda }^{\prime }) \\ &{}\le U(M(\varvec{\lambda }^{\prime })\vert \varvec{\lambda }^{\prime }) + K(\varvec{\lambda }^{\prime }\vert \varvec{\lambda }^{\prime }). \end{array} $$

  Then from (18) and again the definition of K we obtain the result.

• (ii) If $\varvec{\lambda }^{*}$ is a fixed point of M then $U(M(\varvec{\lambda }^{*})\vert \varvec{\lambda }^{*}) = U(\varvec{\lambda }^{*}\vert \varvec{\lambda }^{*})$. Equation (19) follows from Definition 1. The reciprocal is obvious by the definition of M.

• (iii) This property immediately follows from (ii) and Lemma 1.

$\square$

Appendix C: Proof of Theorem 3

Proof

(i) The fact that the sequence $(F(\varvec{\lambda }^{(n)}))$ is non-increasing is a direct consequence of Theorem 2(i). It is clear that $(\varvec{\lambda }^{(n)})$ is bounded. Indeed, if this would not be the case, a sub-sequence $(\varvec{\lambda }^{(n_{k})})$ such that $\lim _{k\rightarrow +\infty }\Vert \varvec{\lambda }^{(n_{k})}\Vert =+\infty$ may be extracted. Since F is coercive, the sequence $(F(\varvec{\lambda }^{(n_{k})}))$ would not be bounded either, which comes in contradiction with the fact that $(F(\varvec{\lambda }^{(n)}))$ is non-increasing and bounded below by the minimum of F. Let $(\varvec{\lambda }^{(n_{k})})$ be a convergent sub-sequence of $(\varvec{\lambda }^{(n)})$ with limit $\varvec{\lambda }^{*}$. Then the sub-sequence $(\varvec{\lambda }^{(n_{k}+1)})$ is also convergent and tends to $M(\varvec{\lambda }^{*})$. The sequence $(F(\varvec{\lambda }^{(n)}))$ being non-increasing and bounded below, it converges and

$$ \lim _{n\rightarrow +\infty } F(\varvec{\lambda }^{(n)}) =\lim _{k\rightarrow +\infty } F(\varvec{\lambda }^{(n_{k})}) = \lim _{k\rightarrow +\infty } F(\varvec{\lambda }^{(n_{k}+1)}), $$

thus $F(M(\varvec{\lambda }^{*}))=F(\varvec{\lambda }^{*})$. From Theorem 2(i) we then deduce that

$$\begin{aligned} \varvec{\lambda }^{*}\in \mathop {\mathrm {arg\,min\,}}\lbrace U(\varvec{\lambda }\vert \varvec{\lambda }^{*})\,:\, \varvec{\lambda }\in X\rbrace , \end{aligned}$$

and from the same theorem it results that $\varvec{\lambda }^{*}$ is a fixed point of M and $\varvec{\lambda }^{*}\in \mathop {\mathrm {arg\,min\,}}\lbrace F(\varvec{\lambda })\,:\, \varvec{\lambda }\in X\rbrace$. Thus

$$ \lim _{n\rightarrow +\infty } F(\varvec{\lambda }^{(n)}) = F(\varvec{\lambda }^{*}) = \min \lbrace F(\varvec{\lambda })\,:\, \varvec{\lambda }\in X\rbrace . $$

(ii) As an immediate consequence of the proof of (i) we have:

$$ F(\varvec{\lambda }^{*})=\lim _{k\rightarrow +\infty } F(\varvec{\lambda }^{(n_{k})}) =\lim _{n\rightarrow +\infty } F(\varvec{\lambda }^{(n)}) = \min \lbrace F(\varvec{\lambda })\,:\, \varvec{\lambda }\in X\rbrace . $$

(iii) If F is strictly convex there is an unique

$$\begin{aligned} \varvec{\lambda }^{*}\in \mathop {\mathrm {arg\,min\,}}\lbrace F(\varvec{\lambda })\,:\, \varvec{\lambda }\in X\rbrace . \end{aligned}$$

From (ii), any convergent sub-sequence of $(\varvec{\lambda }^{(n)})$ has to converge to $\varvec{\lambda }^{*}$, thus the sequence converges to the same limit. $\square$

Appendix D: Proof of Lemma 4

Proof

From the definition of the Fenchel-Legendre transform we have:

$$\begin{aligned} H^{*}(\varvec{p}) = \sup _{\varvec{u}} \left\{ \langle \varvec{u}, \varvec{p}\rangle - H(\varvec{u}) \right\} \\ = \sup _{\varvec{u}} \left\{ \langle \varvec{u}, \varvec{p}- \varvec{s}\rangle + \langle \widetilde{\varvec{\lambda }}\log \varvec{u}, \varvec{s}\rangle \right\} . \end{aligned}$$

If for some $j\in \{1,\dots ,J\}$, $p_{j}\ge s_{j}$ then

$$\begin{aligned} \lim \limits _{u_{j}\rightarrow +\infty } \left( u_{j}(p_{j}-s_{j})+s_{j}\widetilde{\lambda }_{j}\log (u_{j})\right) =+\infty . \end{aligned}$$

By taking all other entries of $\varvec{u}$ as ones, it can be seen that $H^{*}(\varvec{p})$ is infinite. In the case ${\varvec{p}}<\varvec{s}$ component-wise, the resolution of the Euler equation associated to the optimization problem shows that

$$\begin{aligned} H^{*}({\varvec{p}}) = \left\langle \frac{\varvec{s}\widetilde{\varvec{\lambda }}}{\varvec{s}- {\varvec{p}} }, {\varvec{p}} \right\rangle - H\left( \frac{\varvec{s}\widetilde{\varvec{\lambda }}}{\varvec{s}- {\varvec{p}} }\right) . \end{aligned}$$

(D4)

Finally, we have

$$\begin{aligned} H^{*}({\varvec{p}})= & {} -\left\langle \frac{\varvec{s}\widetilde{\varvec{\lambda }}}{\varvec{s}- {\varvec{p}} },\varvec{s}- {\varvec{p}} \right\rangle +\left\langle \log \frac{\varvec{s}\widetilde{\varvec{\lambda }}}{\varvec{s}- {\varvec{p}}},\varvec{s}\widetilde{\varvec{\lambda }}\right\rangle \\\nonumber= & {} \left\langle \log (\varvec{s}\widetilde{\varvec{\lambda }})-\varvec{1},\varvec{s}\widetilde{\varvec{\lambda }}\right\rangle - \left\langle \log (\varvec{s}- {\varvec{p}}),\varvec{s}\widetilde{\varvec{\lambda }}\right\rangle \\= & {} C- \left\langle \log (\varvec{s}- {\varvec{p}}),\varvec{s}\widetilde{\varvec{\lambda }}\right\rangle , \end{aligned}$$

with $C=\left\langle \log (\varvec{s}\widetilde{\varvec{\lambda }})-\varvec{1},\varvec{s}\widetilde{\varvec{\lambda }}\right\rangle$ independent from $\varvec{p}$. $\square$

Appendix E: Proof of Lemma 5

Proof The Fenchel-Legendre transform of the total variation G is (see [32]):

$$\begin{aligned} G^{*}(\varvec{p}) = \delta _{\alpha K}(\varvec{p}) , \end{aligned}$$

(E5)

where $K = \lbrace \textrm{div}_{d} \varvec{\varphi }{:} \varvec{\varphi }\in Z \text { such that } \vert \varphi _{j}\vert \le 1, \ j=1,\dots ,J \rbrace ,$ and $\delta$ is the indicator function,

$$ \delta _{\alpha K}(\varvec{p}) = \left\{ \begin{array}{cl} 0 &{} \text { if } \varvec{p}\in \alpha K \\ +\infty &{} \text {otherwise} \end{array} \right. . $$

Thus the solution $\varvec{p}^{*}$ of the dual problem can be expressed as $\varvec{p}^{*} = -\alpha \textrm{div}_{d} \varvec{\varphi }^{*}$, and the conclusion follows now from Lemma 4.

Appendix F: Proof of Lemma 6

Proof

If $\varvec{\varphi }\in S$, we have $\Vert \textrm{div}_{d} \varvec{\varphi }\Vert _{\infty } \le 4$ and then for $\alpha < s_{\textrm{min}}/4$ we obtain

$$\begin{aligned} \varvec{s}+\alpha \textrm{div}_{d} \varvec{\varphi }\ge s_{\textrm{min}} - 4 \alpha > 0. \end{aligned}$$

The function h is thus finite and continuous on the compact set S. h is thus bounded and $\varvec{\varphi }^{*}$ exists. The positivity of $\varvec{u}^{*}$ is obvious from the previous inequality. $\square$