正投影与反投影
1. 正投影
正投影（即Radon变换）
\[p(s,\theta)=\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)\delta(x\cos\theta+y\sin\theta - s)dxdy\]
2. 反投影
\[b(x,y)=\int\limits_{0}^{\pi}p(s,\theta)|_{s=x\cos\theta+y\sin\theta}d\theta=\frac{1}{2}\int\limits_{0}^{2\pi}p(x\cos\theta+y\sin\theta,\theta)d\theta\]
\(\colorbox{Orange}{注意上述反投影过程并没有得到f(x,y)，而是b(x,y)，反投影运算不是正投影的逆运算。}\)
3. 矩阵简化示意
将被投影区域离散化后，平行束在各角度探测器上的正投影的过程可以表示为：
\[\vec{P}=A\vec{F}\]
其中，投影值的列向量\(\vec{P}\)元素个数为\((L_{\theta}\times K_{s})\)，被投影域的值向量\(\vec{F}\)元素个数为\((M_{column}\times N_{row})\)，为简单起见，认为投影矩阵\(A\)的元素为非0即1。
直接反投影的矩阵表示：
\[\vec{B}=A^T\vec{P}\]
反投影过程\(A^T\)只是正投影过程\(A\)的转置，并非\(A^{-1}\)（实际上\(A\)也不一定有\(A^{-1}\)），考虑到简化处理\(A\)中元素非1即0，反投影得到的\(\vec{B}\)只是\(\vec{P}\)中元素的简单相加（直接反投影操作）。
\(\colorbox{Orange}{实际上考虑到投影线与体素相交的线段长度不一，上述矩阵A的元素并非0即1}\)
\(\colorbox{Orange}{反投影矩阵在离散计算时有时却通常非0即1（即只要体素在反投影线上即可，不考虑相交线段长度）}\)
\(\colorbox{Orange}{因此反投影矩阵严格上也不是正投影矩阵的转置，此处仅简化考虑，具体参见重建的迭代算法}\)
4. 正投影的逆变换
那理论上如何基于\(p(s,\theta)\)得到\(f(x,y)\)呢，利用傅立叶变换。
Step1:
\[p(s,\theta)=\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)\delta(x\cos\theta+y\sin\theta - s)dxdy\]
两侧同时关于\(s\)进行傅立叶变换
\[S(\omega,\theta)=\int\limits_{-\infty}^{+\infty}p(s,\theta)e^{-j\omega s}ds=\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)\delta(x\cos\theta+y\sin\theta - s)dxdye^{-j\omega s}ds\]
交换右侧积分次序：
\[S(\omega,\theta)=\int\limits_{-\infty}^{+\infty}p(s,\theta)e^{-j\omega s}ds=\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}[\int\limits_{-\infty}^{+\infty}f(x,y)\delta(x\cos\theta+y\sin\theta - s)e^{-j\omega s}ds]dxdy\]
利用\(\delta(x)\)函数性质，等式右侧等于：
\[S(\omega,\theta)=\int\limits_{-\infty}^{+\infty}p(s,\theta)e^{-j\omega s}ds=\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)e^{-j\omega(x\cos\theta + y\sin\theta)}dxdy= \int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)e^{-j[(\omega \cos\theta)x +(\omega\sin\theta)y]}dxdy\]
由于的\(f(x,y)\)二维傅立叶变换\(F(u,v)\)为:
\[F(u,v) = \int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)e^{-j2\pi(ux+vy)}dxdy\]
对比可知：
\[S(\omega,\theta) = F(\omega \cos\theta,\omega \sin\theta)=F(u,v)|_{u=w\cos\theta,v=w\sin\theta}=F_{polar}(\omega,\theta)\]
(忽略对\(2\pi\)的处理，实际操作中保持正、逆变换的自洽一致即可)
即\(F(u,v)\)平面的一条线(极径\(r =w\)，极角为\(\theta\)的一条线)，也即中心切片定理
Step2
再对上述进行二维傅立叶逆变换，即可得到\(f(x,y)\)
\[F(u,v) = \int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}f(x,y)e^{-j2\pi(ux+vy)}dxdy\]
\[f(x,y) = \int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}F(u,v)e^{j2\pi(ux+vy)}dudv\]
（上述公式形式上未考虑\(2\pi\)系数）
具体来说：
\[f(x,y) = \int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}F(u,v)|_{u=w\cos\theta,u=w\sin\theta}e^{j2\pi(ux+vy)}dudv\]
进行坐标变换：
\[u=\omega\cos\theta,v=\omega\sin\theta\]
\[dudv = \omega d\omega d\theta\]
可得
\[f(x,y) = \int\limits_{0}^{2\pi}\int\limits_{0}^{+\infty}F_{polar}(\omega,\theta)e^{j2\pi\omega(x\cos\theta +y\sin\theta)}\omega d\omega d\theta\]
考虑到
\[F_{polar}(\omega,\theta + \pi) = F_{polar}(-\omega,\theta)\]
则：
\[f(x,y) = \int\limits_{0}^{\pi}\int\limits_{-\infty}^{+\infty}F_{polar}(\omega,\theta)e^{j2\pi\omega(x\cos\theta +y\sin\theta)}|\omega| d\omega d\theta=\int\limits_{0}^{\pi}\int\limits_{-\infty}^{+\infty}S(\omega,\theta)e^{j2\pi\omega(x\cos\theta +y\sin\theta)}|\omega| d\omega d\theta\]
令：
\[s= x\cos\theta +y\sin\theta\]
\[q(s,\theta)=\int\limits_{-\infty}^{+\infty}S(\omega,\theta)e^{j2\pi\omega s}|\omega| d\omega\]
\(\colorbox{Orange}{最终得到原投影域的值为：}\)
\[f(x,y)=\int\limits_{0}^{\pi}\int\limits_{-\infty}^{+\infty}S(\omega,\theta)e^{j2\pi\omega(x\cos\theta +y\sin\theta)}|\omega| d\omega d\theta=\int\limits_{0}^{\pi}q(s,\theta)|_{s=x\cos\theta+y\sin\theta}d\theta\]
5. 总结
对比上述过程：
\[b(x,y)=\int\limits_{0}^{\pi}p(s,\theta)|_{s=x\cos\theta+y\sin\theta}d\theta\]
\[f(x,y)=\int\limits_{0}^{\pi}q(s,\theta)|_{s=x\cos\theta+y\sin\theta}d\theta\]
其中：
\[p(s,\theta)=\int\limits_{-\infty}^{+\infty}S(\omega,\theta)e^{j2\pi\omega s}d\omega\]
\[q(s,\theta)=\int\limits_{-\infty}^{+\infty}S(\omega,\theta)e^{j2\pi\omega s}|\omega| d\omega=\int\limits_{-\infty}^{+\infty}|\omega|S(\omega,\theta)e^{j2\pi\omega s} d\omega\]
\[S(\omega,\theta)=\int\limits_{-\infty}^{+\infty}p(s,\theta)e^{-j\omega s}ds\]
\(S(\omega,\theta)\)为原始投影数据\(p(s,\theta)\)的直接傅立叶变换（再逆变换回来得到的仍为原始数据\(p(s,\theta)\)），若在逆变换前进行\(|\omega|\)滤波后再进行傅立叶逆变换，则得到\(q(s,\theta)\)。如果直接使用投影数据\(p(s,\theta)\)进行反投影操作，则得到是\(b(x,y)\)，若用\(q(s,\theta)\)进行直接反投影操作，则得到的为最终所求的\(f(x,y)\)
获取\(b(x,y)\)的直接反投影过程等价于：
- 傅立叶变换
- 对1中结果傅立叶逆变换
- 反投影
上述前两步实际上是抵消的，等效于只有第3步
获取\(f(x,y)\)的滤波反投影过程等价于
- 傅立叶变换
- 对1中结果进行滤波\(|\omega|\)
- 对2中结果傅立叶逆变换
- 反投影
\(\colorbox{green}{总之，理论上反投影操作本身并不妨碍还原得到原始f(x,y)的值}\)
\(\colorbox{green}{关键在于不能直接用投影数据直接反投影，对投影数据进行加工后再进行反投影,理论上可以还原得到f(x,y)}\)
理论上通过\(|\omega|\)滤波反投影过程可以还原得到\(f(x,y)\)，但实际过程中由于离散化等因素，不可能完全还原得到\(f(x,y)\)，因此也发展出各种类型的滤波核。
Ref