CT重建笔记（四）——三维重建
人如果不思考不学习，天天刷短视频，跟咸鱼有什么区别？
平行的线积分数据（X射线变换）
平行光束图像重建的理论基础是中心切片定理（二维情形见我的博客https://leslielee.blog.csdn.net/article/details/134722426）。
三维平行光束的中心切片定理：
中心切片定理表示为：位于
Π
(
θ
)
\Pi(\theta)
Π(θ)的
F
(
w
x
,
w
y
,
w
z
)
F(w_x,w_y,w_z)
F(wx,wy,wz) 等于
P
(
w
u
,
w
v
)
P(w_u,w_v)
P(wu,wv)。
其中，
三维重建体
f
(
x
,
y
,
z
)
f(x,y,z)
f(x,y,z)在投影角度为
θ
\theta
θ的投影
p
(
u
,
v
)
p(u,v)
p(u,v)的傅里叶变换为
P
(
w
u
,
w
v
)
P(w_u,w_v)
P(wu,wv)；
三维重建体
f
(
x
,
y
,
z
)
f(x,y,z)
f(x,y,z)的傅里叶变换为
F
(
w
x
,
w
y
,
w
z
)
F(w_x,w_y,w_z)
F(wx,wy,wz)；
过原点且与探测器平面平行的平面为
Π
(
θ
)
\Pi(\theta)
Π(θ)。
基于中心切片定理，可以确定出哪些
θ
\theta
θ轨迹可以满足
F
(
w
x
,
w
y
,
w
z
)
F(w_x,w_y,w_z)
F(wx,wy,wz)的每一点都被测到。
比如圆轨迹，扫描过程中平行束始终与探测器法向量平行，与x-y平面的法向量垂直。但由于对称性，扫180度圆轨迹便可得到完整数据。
Orlov’s condition：如果单位球上的每个大圆都与
Ω
\Omega
Ω有交点，则扫描数据完整。
其中，平行束的方向为单位向量
θ
\theta
θ，
θ
\theta
θ的轨迹为
Ω
\Omega
Ω，大圆（great circle）是单位球上最大的圆周。
紧集（compact set）、完备集（perfect set）和完全集（complete set）的对比分析h. ttps://www.cnblogs.com/readalps/p/16671249.html
投影/反投影PSF 是指 将Dirac Delta函数先投影再反投影得到的响应。该PSF是空间移不变的。
b
=
f
∗
∗
∗
h
b=f***h
b=f∗∗∗h
其中，
h
h
h为投影/反投影PSF，
b
b
b为反投影，
f
f
f为重建体，
∗
∗
∗
***
∗∗∗为三维卷积。
2D情形下，360度扇束与360度平行束的
h
h
h相同，
h
=
1
/
x
2
+
y
2
h=1/\sqrt{x^2 + y^2}
h=1/x2+y2。
3D情形下，若
Ω
\Omega
Ω是整个单位球面，则平行束
h
=
1
/
(
x
2
+
y
2
+
z
2
)
h=1 / (x^2+y^2+z^2)
h=1/(x2+y2+z2)
疑问：两个脉冲响应相比，其中一个更快趋于0表明什么？
F
{
1
/
x
2
+
y
2
}
=
1
/
w
x
2
+
w
y
2
F\{ 1/\sqrt{x^2 + y^2} \} = 1/\sqrt{w_x^2+w_y^2}
F{1/x2+y2}=1/wx2+wy2
F
{
1
/
(
x
2
+
y
2
+
z
2
)
}
=
π
/
w
x
2
+
w
y
2
+
w
z
2
F\{ 1/(x^2 + y^2 + z^2) \} = \pi/\sqrt{w_x^2+w_y^2 + w_z^2}
F{1/(x2+y2+z2)}=π/wx2+wy2+wz2
推导这个推荐用deepseek,很方便
若已知三维斜坡滤波器
G
G
G，则对投影进行反投影再用
G
G
G进行三维滤波。
若想先滤波再反投影，二维滤波器可通过用
Π
(
θ
)
\Pi(\theta)
Π(θ)对三维滤波器截取得到。
平行的面积分数据（Radon变换）
射线源发出等间隔平行的平面束，被一维探测器所接收，探测器单元采集到的物体的面积分。其实面积分数据在实际中不常见。
三维平行面束的中心切片定理：
位于
L
(
θ
)
L(\theta)
L(θ)的
F
(
w
x
,
w
y
,
w
z
)
F(w_x,w_y,w_z)
F(wx,wy,wz) 等于
P
(
w
,
θ
ˉ
)
P(w,\bar{\theta})
P(w,θˉ)。
其中，
一维探测器的方向向量为
θ
ˉ
\bar{\theta}
θˉ；
与探测器平行且过原点的直线为
L
(
θ
ˉ
)
L(\bar{\theta})
L(θˉ)；
三维重建体
f
(
x
,
y
,
z
)
f(x,y,z)
f(x,y,z)的投影为
p
(
s
,
θ
ˉ
)
p(s,\bar{\theta})
p(s,θˉ)，其傅里叶变换为
P
(
w
,
θ
ˉ
)
P(w,\bar{\theta})
P(w,θˉ)；
三维重建体
f
(
x
,
y
,
z
)
f(x,y,z)
f(x,y,z)的傅里叶变换为
F
(
w
x
,
w
y
,
w
z
)
F(w_x,w_y,w_z)
F(wx,wy,wz)。
三维Radon反演公式：先求解投影二阶导数
d
2
p
(
s
,
θ
ˉ
)
d
s
2
\frac{d^2 p(s,\bar{\theta})}{ds^2}
ds2d2p(s,θˉ)，再进行反投影。
这里的反投影是将一维探测器上探元的导数值反投影至该探元对应的平面。即一个点扩展成一个面，而X射线变换中的反投影是将一个点扩展为一条线。
Radon变换反投影也可以分两步，先扩展出一条直线，再由直线上的每一点再扩展出一条直线。
锥形束数据
1983年Tuy发表的论文中有一个引理，几乎等同于锥形束的中心切片定理。
Tuy发表的论文还给出了锥形束数据足量条件，Tuy条件表述为：每一个与物体相交的平面都必须包含至少一个锥形束的焦点位置。
圆轨迹CBCT不满足Tuy条件，圆圈加直线或螺旋CBCT满足Tuy条件。
Feldkamp算法
将锥形束重建转化为扇形束重建。
Feldkamp算法：由Feldkamp，Davis，和 Kress 三人共同发表，用于重建圆轨迹CBCT。
远离圆轨道平面的区域伪影更严重，伪影表现为：图像数值较低，层间图像互串，物体边缘出现负值。
疑问：物体沿轴向是一个常数，则Feldkamp算法可精确重建，但如果沿轴向物体内部有不同缺陷呢？
锥形张角：焦点与中心列两端点形成的角度
Feldkamp算法
1、将投影乘以
c
o
s
α
cos\alpha
cosα，
α
\alpha
α为射线与中心射线的夹角
2、将处理后的投影逐行进行斜坡滤波
3、对滤波数据进行加权反投影，权重为射线长度
Grangeat算法
用三维Radon反演公式重建锥形束。
对投影进行线积分来构造物体的面积分投影。但这样得到是加权的Radon投影，权重为
1
/
r
1/r
1/r。
∫
L
c
o
n
e
b
e
a
m
t
r
a
n
s
f
o
r
m
=
1
r
R
a
d
o
n
t
r
a
n
s
f
o
r
m
\int_L conebeam \space transform = \frac{1}{r} Radon \space transform
∫Lconebeam transform=r1Radon transform
r
d
α
=
d
t
r d\alpha = dt
rdα=dt
联立两式得到：
∫
L
c
o
n
e
b
e
a
m
t
r
a
n
s
f
o
r
m
=
d
α
d
t
R
a
d
o
n
t
r
a
n
s
f
o
r
m
\int_L conebeam \space transform = \frac{d\alpha}{dt} Radon \space transform
∫Lconebeam transform=dtdαRadon transform
∂
(
∫
L
c
o
n
e
b
e
a
m
t
r
a
n
s
f
o
r
m
)
∂
α
=
∂
(
R
a
d
o
n
t
r
a
n
s
f
o
r
m
)
∂
t
\frac{\partial (\int_L conebeam \space transform)}{\partial \alpha} = \frac{\partial (Radon \space transform)}{\partial t}
∂α∂(∫Lconebeam transform)=∂t∂(Radon transform) 式1
三维Radon反演公式是对Radon投影进行二阶导再进行Radon反投影。由式1得到的结果相当于对Radon投影求了一阶导数，只需再求一阶导数便可进行反投影。
不早了，晚安，有空继续更新本文，(●'◡'●)