单光子发射计算机断层成(SPECT)被广泛研究于医学和数学领域.在临床诊断中最重要的应用是SPECT投影数据的重建问题.在过去的几年中, 许多研究人员利用迭代和解析方法进行图像重建^[1-4].迭代算法, 能更好的利用成像系统模型包括各种物理因素和复杂的数据获取几何图形, 但与解析方法相比需要大量的计算.而解析方法, 则将投影数据视为完美的线积分, 并利用数学解析公式来有效的进行重建.我们知道, Radon变换及其反演公式在图像重建中起着重要的作用.一般来说, 光子衰减的投影数据由衰减Radon变换来表示, 定义为

$ \begin{eqnarray} [{\bf R}_{\mu}f](s, \theta)=\int^{\infty}_{-\infty}{\rm e}^{-[{\bf D}\mu]_{\theta}(s, t)}f(s\vec{\theta}+t\vec{\theta}^{\bot}){\rm d}t, %\eqno{(1.1)} \end{eqnarray} $

(1.1)

其中

$ \begin{eqnarray} [{\bf D}\mu]_{\theta}(s, t)=\int^{\infty}_{t}\mu_{\theta}(s, \tau){\rm d}\tau . %\eqno{(1.2)} \end{eqnarray} $

(1.2)

我们记$\vec{\theta}=(\cos\theta, \sin\theta)$, $\vec{\theta}^{\bot}=(-\sin\theta, \cos\theta)$.令$\mu_{\theta}$表示衰减系数.在SPECT中, $[R_{\mu}f](s, \theta)$和$f(x, y)$均为非负.反演公式在文献[5-7]中已得到充分的证实. Novikov和Natterer最近利用新复方法建立了反演公式.当衰减在体内是均匀的, 即$\mu_{\theta}$是常数, 衰减Radon变换就退化为指数型Radon变换, 各种反演公式可见文献[8-10].其中使用最广泛的重建方法是FBP算法, 即滤波背投影算法, 本文中也对其进行讨论.更进一步, 如果衰减可以忽略, 也就是说, $\mu_{\theta}=0$, 指数型Radon变换即为古典Radon变换, 反演公式可见文献[11-12].

在本文中, 我们用了两种不同的方法去得到反演公式.一种是利用经典的FBP算法^[12].我们得到了平行束几何扫描下的反演公式, 然后将结果拓展到扇形束投影.不同于其他方法, 我们直接在扇形束投影下计算公式中的偏导数, 所以在最终的结果中不包含偏导算子, 最终的反演公式也因此变得更加简洁, 更容易应用.另一种方法是利用对偶算子构建中间函数来得到反演公式.在这个方法中, 背投影和加权Hilbert变换对修复图像起了重要的作用.我们利用对偶算子构建中间函数, 通过这个中间函数和加权Hilbert变换我们可以更简单的得到反演公式.我们研究了在平行束几何的半扫描和对扇形光束几何的短扫描中获得的均匀衰减的数据取来重建图像.最后我们对有限加权的Hilbert变换进行了数值模拟^[9].

首先我们讨论研究指数型Radon变换以及它的反演公式.我们用下列形式定义指数型Radon变换

$ \begin{eqnarray} [{\bf T}_{\mu}f](\theta, s)=\int^{\infty}_{-\infty}{\rm e}^{\mu t}f(s\vec{\theta}+t\vec{\theta}^{\bot}){\rm d}t, %\eqno{(1.3)} \end{eqnarray} $

(1.3)

它的对偶算子, 通常叫做背投影运算, 定义为

$ \begin{eqnarray} [{\bf T}^{\sharp}_{\mu}g](\vec{r})=\int_{S^{1}}{\rm e}^{\mu\vec{r}\cdot\vec{\theta}^{\bot}}g(\theta, \vec{r}\cdot\vec{\theta}){\rm d}\theta, %\eqno{(1.4)} \end{eqnarray} $

(1.4)

其中$s=\vec{r}\cdot\vec{\theta}$.为简便起见, 在本文中, 我们令$p(\theta, s)=[{\bf T}_{\mu}f](\theta, s)$.

引理2.1  令$ f(\vec{r})\in C^{\infty}_{0}(\mathbb{R}^2)$, $ g\in C^{\infty}_{0}(S^1\times \mathbb{R}^1)$则

$ \begin{eqnarray} [({\bf T}^{\sharp}_{-\mu}g)\ast f](\vec{r})={\bf T}^{\sharp}_{-\mu}(g\ast {\bf T}_{\mu}f)(\vec{r}). %\eqno{(2.1)} \end{eqnarray} $

(2.1)

证  从已知条件出发, 我们有

$ \begin{eqnarray*} g\ast {\bf T}_{\mu}f(\theta, s)&=&\int_{\mathbb{R}^1}g(\theta, s-\sigma){\bf T}_{\mu}f(\theta, \sigma){\rm d}\sigma\\ &=&\int_{\mathbb{R}^1}g(\theta, s-\sigma)\int_{\mathbb{R}^1}{\rm e}^{\mu t}f(\sigma\vec{\theta}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\sigma. \end{eqnarray*} $

因此

$ \begin{eqnarray*} {\bf T}_{-\mu} ^\sharp(g\ast {\bf T}_{\mu}f)(\vec{r})&=&\int_{S^{1}}{\rm e}^{-\mu \vec{r}\cdot\vec{\theta}^{\bot}}g\ast {\bf T}_{\mu}f(\theta, \vec{r}\cdot\vec{\theta}){\rm d}\theta\\ &=&\int_{S^{1}}{\rm e}^{-\mu \vec{r}\cdot\vec{\theta}^{\bot}}\int_{\mathbb{R}^1}\int_{\mathbb{R}^1}g(\theta, \vec{r}\cdot\vec{\theta}-\sigma){\rm e}^{\mu t}f(\sigma\vec{\theta}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\sigma {\rm d}\theta. \end{eqnarray*} $

记$\vec{j}=\sigma\vec{\theta}+t\vec{\theta}^{\bot}$, 则有$t=\vec{j}\cdot\vec{\theta}^{\bot}, \sigma=\vec{j}\cdot\vec{\theta}$, 更进一步有

$ \begin{eqnarray*} {\bf T}_{-\mu} ^\sharp(g\ast {\bf T}_{\mu}f)(\vec{r})&=&\int_{\mathbb{R}^{2}}f(\vec{j}){\rm d}\vec{j}\int_{S^{1}}g(\theta, (\vec{r}-\vec{j})\vec{\theta}){\rm e}^{-\mu(\vec{r}-\vec{j})\cdot\vec{\theta}^{\bot}}{\rm d}\theta\\ &=&\int_{\mathbb{R}^{2}}f(\vec{j})({\bf T}_{-\mu} ^\sharp g)(\vec{r}-\vec{j}){\rm d}\vec{j}\\ &=&({\bf T}_{-\mu} ^\sharp g)\ast f(\vec{r}). \end{eqnarray*} $

然后我们得到

$ ({\bf T}^{\sharp}_{-\mu}g)\ast f={\bf T}^{\sharp}_{-\mu}(g\ast {\bf T}_{\mu}f). $

证毕.

定理2.1  令$ f(\vec{r})\in C^{\infty}_{0}(\mathbb{R}^2)$, ${\bf T}_{\mu}f $是$ f(\vec{r})$的指数型Radon变换, 则

$ \begin{eqnarray} f(\vec{r})=\frac{1}{4{\pi}^2} \int_{S^{1}} {\rm e}^{-\mu \vec{r}\cdot\vec{\theta}^\perp} \int_{\mathbb{R}^1} \frac{\cos(\mu(\vec{r}\cdot\vec{\theta}-p))}{\vec{r}\cdot\vec{\theta}-p}\frac{\partial}{\partial p}({\bf T}_{\mu}f){\rm d}p{\rm d}\theta. %\eqno{(2.2)} \end{eqnarray} $

(2.2)

证  在引理2.1中, 利用(2.1)式得到一个反演公式, 如果我们找到一个$g$满足$T^{\sharp}_{-\mu}g=\delta$, 其中$\delta$是delta函数.在方程(1.4)中对$p$进行积分, 我们得到关于引理2.1的类比

$ g(\theta, p)=\frac{1}{4\pi^{2}}\frac{\partial}{\partial p}\left(\frac{\cos(\mu p)}{p}\right). $

因此, 方程(1.4)意味着我们必须证明等式

$ \begin{equation} \frac{1}{4\pi^{2}}\int_{S^{1}}{\rm e}^{-\mu \vec{r}\cdot\theta^\perp}\frac{\partial}{\partial p}\left(\frac{\cos(\mu p)}{p}\right)\bigg|_{p=\vec{r}\cdot\vec{\theta}}{\rm d}\theta=\delta(\vec{r}) . %\eqno{(2.3)} \end{equation} $

(2.3)

记

$ \begin{equation} \int_{\mathbb{R}^1}\frac{\cos(\mu p)}{p}{\rm e}^{{\rm i}ps}{\rm d}p=\varphi(s). %\eqno{(2.4)} \end{equation} $

(2.4)

显然, $\varphi^{\prime}(s)=\pi{\rm i}(\delta(s-\mu)+\delta(s+\mu))$.由于$\varphi(s)$是奇函数, 我们得到

$ \begin{equation} \varphi(s)=\left\{ \begin{array}{ll} -\pi{\rm i}, ~~ &s<-\mu, \\ 0, &|s|<\mu, \\ \pi{\rm i} , &s>\mu. \end{array} \right. \end{equation} $

(2.5)

利用上述方程, 我们得到等式

$ \frac{\partial}{\partial p}\left(\frac{\cos(\mu p)}{p}\right)=\frac{1}{2}\int_{|t|>\mu}|t|{\rm e}^{{\rm i}tp}{\rm d}t, $

代入到(2.3)式得到

$ \begin{eqnarray} &&\frac{1}{8\pi^{2}}\int_{|t|>\mu}|t|\int_{S^{1}}{\rm e}^{-\mu(\vec{r}\cdot\vec{\theta}^\perp)}{\rm e}^{{\rm i}t\vec{\theta}\cdot \vec{r}}{\rm d}\theta {\rm d}t \\ &=&\frac{1}{4\pi^{2}}\int^{\infty}_{\mu}t\int_{S^{1}}{\rm e}^{{\rm i}t\vec{\theta}\cdot \vec{x}}\cosh(\mu\vec{\theta}^\perp\cdot \vec{x}){\rm d}\theta {\rm d}t\\ &=&\frac{1}{4\pi^{2}}\int^{\infty}_{0}s\int_{S^{1}}{\rm e}^{{\rm i}\sqrt{s^{2}+\mu^{2}}\vec{\theta}\cdot \vec{x}}\cosh(\mu\vec{\theta}^\perp\cdot \vec{x}){\rm d}\theta{\rm d}s, \end{eqnarray} $

(2.6)

其中我们利用了变换$t^{2}=s^{2}+\mu^{2}$.

我们记(2.6)式的内部积分为$C(\mu)$.假设$C(\mu)$不依赖于$\mu$.那么, 在(2.6)式中令$\mu=0$, 我们注意到(2.6)式的右边是一个在柱坐标系中1的傅里叶逆变换的Delta函数.它还有待证明$C^{\prime}(\mu)=0$.令$|\vec{r}|=r$.则

$ \begin{eqnarray*} C(\mu)&=&\int^{2\pi}_{0}{\rm e}^{{\rm i}\sqrt{s^{2}+\mu^{2}}r\cos\theta}\cosh(\mu r\sin\theta){\rm d}\theta\\ &=&\int^{2\pi}_{0}\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &&+\int^{2\pi}_{0}\sin(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &=&\int^{2\pi}_{0}\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta) {\rm d}\theta\\ &&+\int^{\pi}_{0}\sin(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &&+\int^{2\pi}_{\pi}\sin(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &=&\int^{2\pi}_{0}\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta) {\rm d}\theta\\ &&+\int^{\pi}_{0}\sin(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &&-\int^{\pi}_{0}\sin(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &=&\int^{2\pi}_{0}\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &=&\int^{\pi}_{-\pi}\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta\\ &=&2\int^{\pi}_{0}\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\cosh(\mu r\sin\theta){\rm d}\theta. \end{eqnarray*} $

因此

$ \begin{eqnarray*} C^{\prime}(\mu)&=&4r\left[\int^{\frac{\pi}{2}}_{0}\sinh(\mu r\sin\theta)\cos(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\sin\theta {\rm d}\theta\right.\\ &&-\left.\int^{\frac{\pi}{2}}_{0}\cosh(\mu r\sin\theta)\sin(\sqrt{s^{2}+\mu^{2}}r\cos\theta)\frac{\mu}{\sqrt{s^{2}+\mu^{2}}}\cos\theta {\rm d}\theta\right]\\ &=&4r\left[\int^{1}_{0}\sinh(\mu r\sqrt{1-t^{2}})\cos(\sqrt{s^{2}+\mu^{2}}rt){\rm d}t\right.\\ &&-\left.\int^{1}_{0}\cosh(\mu rt)\sin(\sqrt{s^{2}+\mu^{2}}r\sqrt{1-t^{2}})\frac{\mu}{\sqrt{s^{2}+\mu^{2}}}{\rm d}t\right]. \end{eqnarray*} $

在最后一个方程右边的积分中积分, 我们得到$C^{\prime}(\mu)=0.$证毕.

下面的定理与定理2.1中的表达式稍有不同.

定理2.2  令$f(\vec{r})\in C^{\infty}_{0}(\mathbb{R}^{2})$, $[{\bf T}_{\mu}f](\theta, s)=p(\theta, s)$.则

$ \begin{eqnarray} f(\vec{r})=\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}} \left[\frac{\partial}{\partial s}\int_{\mathbb{R}^1}\frac{\cos(\mu(s-l))}{s-l}p(\theta, l){\rm d}l\right]|_{s=\vec{r}\cdot\vec{\theta}}{\rm d}\theta. %\eqno{(2.7)} \end{eqnarray} $

(2.7)

证  首先, 我们计算

$ \begin{eqnarray*} &&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}\vec{\theta}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}} \int_{\mathbb{R}^1}\frac{\cos(\mu(s-l))}{s-l}p(\theta, l){\rm d}l|_{s=\vec{r}\cdot\vec{\theta}}{\rm d}\theta\\ &=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}\vec{\theta}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}}\int_{\mathbb{R}^1} \frac{\cosh({\rm i}\mu(s-l))}{s-l}\int_{\mathbb{R}^1}f(l\vec{\theta}+t\vec{\theta}^{\perp}){\rm e}^{\mu t}{\rm d}t{\rm d}l|_{s=\vec{r}\cdot\vec{\theta}}{\rm d}\theta, \end{eqnarray*} $

令$l=\vec{\theta}\cdot\vec{w}$, $t=\vec{\theta}^{\bot}\cdot\vec{w}$, 我们得到

$ \begin{eqnarray*} &&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}\vec{\theta}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}}\int_{\mathbb{R}^1} \frac{\cosh({\rm i}\mu(s-l))}{s-l}\int_{\mathbb{R}^1}f(l\vec{\theta}+t\vec{\theta}^{\perp}){\rm e}^{\mu t}{\rm d}t{\rm d}l|_{s=\vec{r}\cdot\vec{\theta}}{\rm d}\theta\\ &=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}\vec{\theta}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}}\int_{\mathbb{R}^1} \frac{\cosh({\rm i}\mu(\vec{r}\cdot\vec{\theta}-\vec{\theta}\cdot\vec{\omega}))}{\vec{r}\cdot\vec{\theta}-\vec{\theta}\cdot\vec{\omega}} \int_{\mathbb{R}^1}f(\vec{w}){\rm e}^{\mu\vec{\theta}^{\perp}\cdot\vec{\omega}}{\rm d}\theta {\rm d}\vec{\omega}\\ &=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}\vec{\theta}{\rm e}^{-\mu(\vec{r}-\vec{\omega})\cdot\vec{\theta}^{\perp}}\int_{\mathbb{R}^1} \frac{\cosh({\rm i}\mu(\vec{r}-\vec{\omega})\cdot\vec{\theta})}{(\vec{r}-\vec{\omega})\cdot\vec{\theta}}\int_{\mathbb{R}^1}f(\vec{w}){\rm d}\theta {\rm d}\vec{\omega}, \end{eqnarray*} $

我们又有

$ \begin{eqnarray*} &&{\rm div}\frac{1}{4\pi^{2}}\int^{2\pi}_{0}\vec{\theta}{\rm e}^{-\mu(\vec{r}-\vec{\omega})\cdot\vec{\theta}^{\perp}}\int_{\mathbb{R}^1} \frac{\cosh({\rm i}\mu(\vec{r}-\vec{\omega})\cdot\vec{\theta})}{(\vec{r}-\vec{\omega})\cdot\vec{\theta}}{\rm d}\theta\\ &=&\frac{1}{4\pi^{2}}2\pi \div\frac{(\vec{r}-\vec{w})}{|\vec{r}-\vec{w}|^{2}} =\frac{1}{2\pi}\div\frac{(\vec{r}-\vec{w})}{|\vec{r}-\vec{w}|^{2}}. \end{eqnarray*} $

则我们得到

$ \begin{eqnarray*} f(\vec{r})&=&\iint_{\mathbb{R}^{2}}f(\vec{w})\delta(\vec{r}-\vec{w}){\rm d}\vec{w} =\frac{1}{2\pi}\iint_{\mathbb{R}^{2}}\div \frac{(\vec{r}-\vec{w})}{|\vec{r}-\vec{w}|^{2}}{\rm d}\vec{w}\\ &=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}} \left[\frac{\partial}{\partial s}\int_{\mathbb{R}^1}\frac{\cos(\mu(s-l))}{s-l}p(\theta, l){\rm d}l\right] \bigg|_{s=\vec{r}\cdot\vec{\theta}}{\rm d}\theta. \end{eqnarray*} $

证毕.

定理2.2给出了平行束扫描投影下的反演公式. 图 1显示了一种以圆形聚焦点轨迹的扇形束数据获取几何图形, 其中每个投影射线可用$(\beta, \sigma)$来表示.

首先我们引入扇形投影坐标, 记为$(\beta, \sigma)$.一个特定的投影射线是从焦点$ S $由角$\beta$和角$\sigma$确定.在本文中, 函数$f(x, y)$的扇形束均匀衰减的投影$D_{\mu}f$定义为

$ \begin{eqnarray} [D_{\mu}f](\beta, \sigma)=\int^{\infty}_{0}f(S+\tau\vec{\alpha}(\beta, \sigma)){\rm e}^{\mu\tau}{\rm d}\tau, %\eqno{(2.8)} \end{eqnarray} $

(2.8)

这里$\sigma\in[-\sigma_{m}, \sigma_{m}]$，$\vec{\alpha}(\beta, \sigma)$是$\mathbb{R}^2$中的单位向量, 表示从焦点到准直孔的方向, 如图 1所示. $\sigma_{m}\in(0, \pi/2)$表示扇形束的最大角.所以定义(1.3)和(2.8)可以由下列公式重新改写

$ \begin{eqnarray} p(\sigma+\beta, R\sin\sigma)={\rm e}^{-\mu R\cos\sigma}[D_{\mu}f](\beta, \sigma)=g(\beta, \sigma). %\eqno{(2.9)} \end{eqnarray} $

(2.9)

我们用$g(\beta, \sigma)$表示(2.8)式的修正投影.令$R$表示原点到焦点的距离, 则在平行束和扇形束几何下的坐标转化关系为

$ \begin{eqnarray} s=R\sin\sigma, \, \theta=\sigma+\beta. %\eqno{(2.10)} \end{eqnarray} $

(2.10)

焦点和$s$轴的距离为$R\cos\sigma$.令$r=|\vec{r}|$, $\vec{r}=(x, y)=r(\cos\varphi, \sin\varphi)$, $K(r, \varphi, \beta)$表示$ SP $的长度, $\sigma'(r, \varphi, \beta)$表示$ SO$与$SP$的夹角.显然以下关系式成立^[9]

$ \begin{equation} K(r, \varphi, \beta)=\sqrt{r^{2}+R^{2}+2rR\sin(\beta-\varphi)}, %\eqno{(2.11)} \end{equation} $

(2.11)

$ \begin{equation} \sigma'(r, \varphi, \beta)=\arctan\frac{r\cos(\beta-\varphi)}{R+r\sin(\beta-\varphi)}, %\eqno{(2.12)} \end{equation} $

(2.12)

$ \begin{equation} \vec{\alpha}(\beta, \sigma)=(\sin(\beta+\sigma), -\cos(\beta+\sigma)). %\eqno{(2.13)} \end{equation} $

(2.13)

定理2.3  令$f(\vec{r})\in C^{\infty}_{0}(\mathbb{R}^{2})$, $\hat{s}=\vec{r}\cdot\vec{\theta}$ $\hat{t}=\vec{r}\cdot\vec{\theta}^{\perp}$, 在扇形束投影下, 则

$ \begin{eqnarray} f(\vec{r})&=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\hat{t}}\int^{\sigma_{m}}_{-\sigma_{m}}(\frac{\mu\cos\mu\hat{s}}{K\sin(\sigma'-\sigma)}p^{fan}_{s}(\beta, \sigma) -\frac{\mu\sin\mu\hat{s}}{K\sin(\sigma'-\sigma)}p^{fan}_{c}(\beta, \sigma)\\ &&-\frac{\cos\mu\hat{s}}{(K\sin(\sigma'-\sigma))^{2}}p^{fan}_{c}(\beta, \sigma) -\frac{\sin\mu\hat{s}}{(K\sin(\sigma'-\sigma))^{2}}p^{fan}_{s}(\beta, \sigma)){\rm d}\sigma {\rm d}\beta, \end{eqnarray} $

(2.14)

其中

$ p^{fan}_{c}(\beta, \sigma)=R\cos\sigma\cos(\mu R\sin\sigma)p^{fan}(\beta, \sigma), \\ p^{fan}_{s}(\beta, \sigma)=R\cos\sigma\sin(\mu R\sin\sigma)p^{fan}(\beta, \sigma). $

证  我们重新改写(2.7)式如下

$ \begin{eqnarray} f(\vec{r})&=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}} \left[\frac{\partial}{\partial s}\int_{\mathbb{R}^1}\frac{\cos(\mu(s-l))}{s-l}p(\theta, l){\rm d}l\right]{\rm d}\theta\\ &=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}}\int_{\mathbb{R}^1} \left(-\frac{\mu\sin(\mu(s-l))}{(s-l)}-\frac{\cos(\mu(s-l))}{(s-l)^{2}}\right)p(\theta, l){\rm d}l{\rm d}\theta\\ &=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\hat{t}}\int_{\mathbb{R}^1}\left(\frac{\mu(\cos\mu s\sin\mu l-\sin\mu s\cos\mu l)}{(s-l)}\right.\\ &&\left.-\frac{\cos\mu s\cos\mu l+\sin\mu s\sin\mu l}{(s-l)^{2}}\right)p(\theta, l){\rm d}l{\rm d}\theta. \end{eqnarray} $

(2.15)

为了简便起见, 我们令$p_{c}(\theta, l)=\cos(\mu l)p(\theta, l)$, $p_{s}(\theta, l)=\sin(\mu l)p(\theta, l)$, 则(2.15)式可表达为如下

$ \begin{eqnarray} f(\vec{r})&=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\hat{t}}\int_{\mathbb{R}^1}\left(\frac{\mu\cos\mu s}{s-l}p_{s}(\theta, l)\right.\\ &&\left.-\frac{\mu\sin\mu s}{s-l}p_{c}(\theta, l)-\frac{\cos\mu s}{(s-l)^{2}}p_{c}(\theta, l)-\frac{\sin\mu s}{(s-l)^{2}}p_{s}(\theta, l)\right){\rm d}l{\rm d}\theta. \end{eqnarray} $

(2.16)

利用坐标转化, 得到

$ |J|=R\cos\sigma, $

所以在扇形束下重建可表达为

$ \begin{eqnarray} f(\vec{r})&=&\frac{1}{4\pi^{2}}\int^{2\pi}_{0}{\rm e}^{-\mu\hat{t}}\int^{\sigma_{m}}_{-\sigma_{m}} \bigg(\frac{\mu\cos\mu\hat{s}}{K\sin(\sigma'-\sigma)}p^{fan}_{s}(\beta, \sigma) -\frac{\mu\sin\mu\hat{s}}{K\sin(\sigma'-\sigma)}p^{fan}_{c}(\beta, \sigma)\\ &&-\frac{\cos\mu\hat{s}}{K(\sin(\sigma'-\sigma))^{2}}p^{fan}_{c}(\beta, \sigma) -\frac{\sin\mu\hat{s}}{K(\sin(\sigma'-\sigma))^{2}}p^{fan}_{s}(\beta, \sigma)\bigg){\rm d}\sigma {\rm d}\beta. \end{eqnarray} $

(2.17)

证毕.

定理2.3为在扇形束投影下的的反演公式, 在最终的结果中并不含偏导算子, 所以这个反演公式更容易去实现.

定理2.4  令$f(\vec{r})\in C^{\infty}_{0}(\mathbb{R}^2)$, 则

$ \begin{eqnarray} [{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]=2\iint_{\mathbb{R}^{2}}\frac{\cosh(\mu|\vec{r}-\vec{w}|)}{|\vec{r}-\vec{w}|}f(\vec{w}){\rm d}\vec{w}. %\eqno{(2.18)} \end{eqnarray} $

(2.18)

证  我们从${\bf T}^{\sharp}_{\mu}$的定义出发

$ \begin{eqnarray*} [{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]&=&\int_{S^{1}}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\bot}}[{\bf T}_{\mu}f](\theta, s){\rm d}\theta\\ &=&\int_{S^{1}}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\bot}} \int_{\mathbb{R}^{1}}{\rm e}^{\mu t}f(s\vec{\theta}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta, \end{eqnarray*} $

利用$s=\vec{r}\cdot\vec{\theta}$和$\vec{r}=(\vec{r}\cdot\vec{\theta})\vec{\theta}+(\vec{r}\cdot\vec{\theta}^{\bot})\vec{\theta}^{\bot}$, 我们得到

$ \begin{eqnarray*} &&\int_{S^{1}}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\bot}} \int_{\mathbb{R}^{1}}{\rm e}^{\mu t}f(s\vec{\theta}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta\\ &=&\int_{S^{1}}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\bot}} \int_{\mathbb{R}^{1}}{\rm e}^{\mu t}f(\vec{r}+(t-\vec{r}\cdot\vec{\theta}^{\bot})\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta\\ &=&\int_{S^{1}}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\bot}} \int_{\mathbb{R}^{1}}{\rm e}^{\mu(t+\vec{r}\cdot\vec{\theta}^{\bot})}f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta\\ &=&\int_{S^{1}}\int_{\mathbb{R}^{1}}{\rm e}^{\mu t}f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta\\ &=&\int_{S^{1}}\int^{\infty}_{0}{\rm e}^{\mu t}f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta+ \int_{S^{1}}\int^{0}_{-\infty}{\rm e}^{\mu t}f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta, \end{eqnarray*} $

其中在内部积分中, 我们做了替换$t=t-\vec{r}\cdot\vec{\theta}^{\bot}$.在最后一个等式的第二个积分中, 令$t\rightarrow-t, \theta\rightarrow\theta+\pi$, 我们得到

$ \begin{eqnarray*} [{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]&=&\int_{S^{1}}\int^{\infty}_{0}{\rm e}^{\mu t}f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta +\int_{S^{1}}\int^{\infty}_{0}{\rm e}^{-\mu t}f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta\\ &=&\int_{S^{1}}\int^{\infty}_{0}({\rm e}^{\mu t}+{\rm e}^{-\mu t})f(\vec{r}+t\vec{\theta}^{\bot}){\rm d}t{\rm d}\theta\\ &=&2\iint_{\mathbb{R}^{2}}\cosh(\mu|\vec{z}|)f(\vec{r}+\vec{z})\frac{1}{|\vec{z}|}{\rm d}\vec{z}, \end{eqnarray*} $

其中我们令$\vec{z}=t\vec{\theta}^{\bot}, t=|\vec{z}|$, 做变量替换$\vec{w}=\vec{r}+\vec{z}$, 最终我们得到

$ [{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]=2\iint_{\mathbb{R}^{2}}\frac{\cosh(\mu|\vec{r}-\vec{w}|)}{|\vec{r}-\vec{w}|}f(\vec{w}){\rm d}\vec{w}. $

证毕.

为了简便起见, 在本文中我们只讨论扇形束短扫描的重建问题.如图 2所示, $y\in(-1, 1)$, $\beta(y)=\arcsin(y/R)$我们假设$R\sin\sigma_{m}=1$, 则图 2中内部的圆盘是单位圆盘$D$.

定理2.5  令$f(\vec{r})\in C^{\infty}_{0}(\mathbb{R}^{2})$, 如果$f(x, y)$在这条直线上的支集是$[-L(y), L(y)]$, 则在扇形束投影下我们得到如下反演公式

$ \begin{eqnarray} f(x, y)&=&\frac{1}{2\pi^{2}}\int^{L(y)}_{-L(y)}\sqrt{\frac{L^{2}(y)-x^{2}}{L^{2} (y)-\tau^{2}}}\bigg[\frac{1}{\tau-x}+k_{\mu, L(y)}(\tau, x)\bigg] \\ &&\times \int^{\pi/2-\beta(y)}_{\beta(y)-\pi/2}{\rm e}^{-\mu\vec{r}\cdot\vec{\alpha}^{\perp}} g(\beta, \sigma')\frac{R\cos\sigma'(r, \varphi, \beta)}{K(r, \varphi, \beta)}{\rm d}\beta {\rm d}\tau, \end{eqnarray} $

(2.19)

其中$k_{\mu, L(y)}(\tau, x)$是一个解析的函数.

证  函数$h(t)$在$I_{q}=(-q, q)$上的加权Hilbert变换$H_{\mu}(s)$定义为

$ \begin{eqnarray} H_{\mu}(s)=\frac{1}{\pi}\int^{q}_{-q}\frac{\cosh(\mu(s-t))}{s-t}h(t){\rm d}t.%\eqno{(2.20)} \end{eqnarray} $

(2.20)

在文献[13]中已经证明

$ \begin{eqnarray} h(t)=\frac{1}{\pi}\int^{q}_{-q}\bigg[\frac{1}{s-t}+k_{\mu, q}(s, t)\bigg]H_{\mu}(s)\sqrt{\frac{q^{2}-t^{2}}{q^{2}-s^{2}}}{\rm d}s, %\eqno{(2.21)} \end{eqnarray} $

(2.21)

其中$k_{\mu, q}(s, t)$是一个解析的函数且满足$k_{0, q}(s, t)\equiv0$.我们知道, 偏导和背投影是局部算子.如果$f(x, y)$在有限区域$[-L(y), L(y)]$中有支集, 如图 3所示.

利用定理2.4, 可以从$[{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]$中得到$f(x, y)$的反演公式, 我们可以得到

$ \frac{1}{2\pi}[{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]= \frac{1}{\pi}\iint_{\mathbb{R}^{2}}\frac{\cosh(\mu|\vec{r}-\vec{w}|)}{|\vec{r}-\vec{w}|}f(\vec{w}){\rm d}\vec{w}, $

所以

$ \begin{eqnarray} f(x, y)=\frac{1}{2\pi^{2}}\int^{L(y)}_{-L(y)}\sqrt{\frac{L^{2}(y)-x^{2}}{L^{2} (y)-\tau^{2}}}\bigg[\frac{1}{\tau-x}+k_{\mu, L(y)}(\tau, x)\bigg] [{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]{\rm d}\tau. %\eqno{(2.22)} \end{eqnarray} $

(2.22)

现在我们可以计算$[{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]$.如图 3所示, 我们将平行光束几何下对$\theta$的积分转化为扇形束下对$\beta$的积分.我们得到

$ \begin{eqnarray*} [{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]&=&\int^{\frac{\pi}{2}}_{-\frac{\pi}{2}}{\rm e}^{-\mu\vec{r}\cdot\vec{\theta}^{\perp}}p(\theta, s){\rm d}\theta =\int^{\pi/2-\beta(y)}_{\beta(y)-\pi/2}{\rm e}^{-\mu\vec{r}\cdot\vec{\alpha}^{\perp}} g(\beta, \sigma')\frac{R\cos\sigma'(r, \varphi, \beta)}{K(r, \varphi, \beta)}{\rm d}\beta. \end{eqnarray*} $

对于给定的$\vec{r}$, $\vec{\theta}$, $\beta$, 满足下列微分关系

$ \begin{equation} {\rm d}\theta={\rm d}(\beta+\sigma'(r, \varphi, \beta))= \frac{R\cos\sigma'(r, \varphi, \beta)}{K(r, \varphi, \beta)}{\rm d}\beta. %\eqno{(2.23)} \end{equation} $

(2.23)

最终我们得到

$ \begin{eqnarray*} f(x, y)&=&\frac{1}{2\pi^{2}}\int^{L(y)}_{-L(y)}\sqrt{\frac{L^{2}(y)-x^{2}} {L^{2}(y)-\tau^{2}}}\bigg[\frac{1}{\tau-x}+k_{\mu, L(y)}(\tau, x)\bigg]\\ &&\times\int^{\pi/2-\beta(y)}_{\beta(y)-\pi/2}{\rm e}^{ -\mu\vec{r}\cdot\vec{\alpha}^{\perp}} g(\beta, \sigma')\frac{R\cos\sigma'(r, \varphi, \beta)}{K(r, \varphi, \beta)} {\rm d}\beta {\rm d}\tau, \end{eqnarray*} $

其中$k_{\mu, L(y)}(\tau, x)$是一个解析函数.证毕.

定理2.4和定理2.5利用对偶算子和加权Hilbert变换来得到反演公式.

定理2.5的关键是有限加权Hilbert变换的反演.我们用$q=1$来进行模拟. (2.21)式中$k_{\mu, 1}(s, t)$内核可以被预计算^[14].这里我们用文献[9]中的数值模拟, 比较容易实现.利用(2.21)式, 当$\mu=0$, 我们可以用下列方程将$h(t)$和$H_{\mu}(s)$联系起来.

$ \begin{equation} h(t)+\int^{1}_{-1}\Psi_{\mu}(t, p)h(p){\rm d}p=\int^{1}_{-1}\frac{H_{\mu}(s)}{\pi(s-t)}\sqrt{\frac{1-t^{2}}{1-s^{2}}}{\rm d}s, %\eqno{(3.1)} \end{equation} $

(3.1)

其中

$ \begin{equation} \Psi_{\mu}(t, p)=\int^{1}_{-1}\frac{A_{\mu}(s-p)}{\pi^{2}(s-t)}\sqrt{\frac{1-t^{2}}{1-s^{2}}}{\rm d}s %\\ =\int^{1}_{-1}\frac{A_{\mu}(s-p)-A_{\mu}(t-p)}{\pi^{2}(s-t)}\sqrt{\frac{1-t^{2}}{1-s^{2}}}{\rm d}s, \end{equation} $

(3.2)

$ \begin{equation} A_{\mu}(t)=(\cosh(\mu t)-1)/t. %\eqno{(3.3)} \end{equation} $

(3.3)

注意到$[A_{\mu}(s-p)-A_{\mu}(t-p)]/(s-t)$是一个光滑的函数, 因此$\Psi_{\mu}(t, p)$核定义了$L^{2}_{w}(I_{1})$上的一个紧积分算子, 记为${\bf \Psi}$.我们把(3.1)式写为如下表达式

$ \begin{equation} [({\bf I}+{\bf \Psi})h](t)=\int^{1}_{-1}\frac{H_{\mu} (s)}{\pi(s-t)}\sqrt{\frac{1-t^{2}}{1-s^{2}}}{\rm d}s. %\eqno{(3.4)} \end{equation} $

(3.4)

假设$({\bf I}+{\bf \Psi})^{-1}$存在, 在SPECT中$\mu$取实值, 用如下反演公式可完全重建原始函数$h(t)$

$ \begin{equation} h(t)=[({\bf I}+{\bf \Psi})^{-1}\int^{1}_{-1}\frac{H_{\mu}(s)} {\pi(s-t)}\sqrt{\frac{1-t^{2}}{1-s^{2}}}{\rm d}s](t). %\eqno{(3.5)} \end{equation} $

(3.5)

方程(3.4)实际上构造了第二种类型的Fredholm算子.同时注意到当$\mu\rightarrow 0$时$\Psi_{\mu}(t, p)$是正的且一致趋于0.因此${\bf \Psi}$可以是一个小$\mu$的收缩映射.在SPECT图像中, 事实上$f(x, y)\geq0$且在$D$中有支集, 因此对于所有的$\mu$来说(3.4)式的反演公式很有可能存在.假设在$(-1, 1)\times(-1, 1)$上$[{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]$存在, 结合(2.22)和(3.4)式, 我们推导出积分方程

$ \begin{equation} [({\bf I}+{\bf \Psi})f](x, y)=\int^{1}_{-1} \frac{[{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f]}{2\pi^{2}(\tau-x)} \sqrt{\frac{1-x^{2}}{1-\tau^{2}}}{\rm d}\tau. \end{equation} $

(3.6)

方程(3.4)是第二种类型的Fredholm积分方程.文献[15]和[16]中讨论的基于矩阵求逆算法是稳定的方式去数值解决(3.4)式.利用对平均采样数据的梯形规则, 现在我们用以下四个步骤来实现(3.5)式的数值.

第一步  初始化

(1) 紧支集区间是[-1, 1], 也就是说，$q=1$.

(2) 采样间隔为$\Delta t=1/M$, 本文中当采取两个一维信号时取M=512, SPECT数据时取M=128.

(3) 支集中的采样点是$t_{m}=(m+0.5)\Delta t$, $-M\leq m\leq M-1$.

第二步  计算$\Psi_{\mu}(t, p)$内核

$ \begin{equation} A_{\mu}(t)=\left\{ \begin{array}{ll} (\cosh(\mu t)-1)/t, ~~&\mid t\mid\geq0.5\Delta t, \\ 0, &\mid t\mid<0.5\Delta t. \end{array} \right. \end{equation} $

(3.7)

$ \begin{equation} M_{\mu}(s, t, p)=\left\{ \begin{array}{ll} \frac{A_{\mu}(s-p)-A_{\mu}(t-p)}{\pi^{2}(s-t)}, ~~&\mid s-t\mid\geq0.5\Delta t, \\[2mm] 0.5(\mu/\pi), &\mid s-t\mid<0.5\Delta t. \end{array} \right. \end{equation} $

(3.8)

$ \begin{equation} \Psi_{\mu}(t_{m}, t_{n})=\Delta t\sqrt{1-t^{2}_{m}}\times\sum\limits^{M-1}\limits_{j=-M}[M_{\mu}(t_{j}, t_{m}, t_{n})/\sqrt{1-t^{2}_{j}}]. %\eqno{(3.9)} \end{equation} $

(3.9)

第三步  执行积分变换

$ \begin{equation} \hat{h}(t_{m})=\int^{1}_{-1}\frac{H_{\mu}(s)}{\pi(s-t_{m})}\sqrt{\frac{1-t^{2}_{m}}{1-s^{2}}}{\rm d}s =\sum\limits^{M-1}\limits_{n=-M}\frac{H_{\mu}(t_{n})}{\pi(m-n)}\sqrt{\frac{1-t^{2}_{m}}{1-t^{2}_{n}}}. \end{equation} $

(3.10)

在实现中, 被积函数在奇异点的数据未被采样.

第四步  求逆矩阵求$h$

$ ({\bf I}+{\bf \Psi})^{-1}=\{\Psi_{\mu}(t_{m}+t_{n})+\delta_{m, n}\}^{-1}, $

其中$\delta_{m, n}$是Kroneckerm delta函数.

$ \begin{equation} h(t_{m})=\sum\limits^{M-1}\limits_{n=-M}\{\Psi(t_{m}, t_{n})+\delta_{m, n}\}^{-1}\hat{h}(t_{n}). %\eqno{(3.11)} \end{equation} $

(3.11)

在步骤2中计算$\Psi(t_{m}, t_{n})$并不包含任何奇异点; 因此, 所有的数值计算都是稳定的.当将这个数值过程应用到图像重建问题(2.22)时, 我们令$H_{\mu}(s)=\frac{1}{2\pi}[{\bf T}^{\sharp}_{-\mu}{\bf T}_{\mu}f](x=s, y)$并且对于一个固定的$y$找到$f(x, y)=h(t=x)$.这个逐行的有限逆过程是对所有$y$坐标的重复.逆矩阵$({\bf I}+{\bf \Psi})^{-1}$仅需被计算一次.

在SPECT模拟中, 目标函数被选为修正后的Shepp-Logan幻影, 如图 4所示. 图 5是通过无噪音的数据来重建图像.

我们可以看到定理2.5更容易实现数值模拟, 而定理2.3更容易进行计算.