CT重建笔记(二)
大部分内容来自wiki与《医学图像重建入门》
X-ray变换
X-ray变换也称ray变换、John变换,由Fritz John在1938,与Radon变换在2D情形下是一致的。
John变换是线积分,Radon变换是超平面积分,2D情形时超平面为直线。
X-ray变换是John方程(一种超双曲型偏微分方程)的解。
John方程:
∂
η
i
,
ϵ
j
2
u
(
ϵ
,
η
)
−
∂
η
j
,
ϵ
i
2
u
(
ϵ
,
η
)
=
0
;
i
,
j
=
1
,
2
,
3
\partial_{\eta_i, \epsilon_j}^2 u(\epsilon, \eta) - \partial_{\eta_j, \epsilon_i}^2 u(\epsilon, \eta) = 0; i,j=1,2,3
∂ηi,ϵj2u(ϵ,η)−∂ηj,ϵi2u(ϵ,η)=0;i,j=1,2,3
X-ray变换:
u
(
ϵ
,
η
)
=
∫
R
f
(
ϵ
+
t
(
η
−
ϵ
)
)
d
t
u(\epsilon, \eta) = \int_R f(\epsilon + t(\eta-\epsilon))dt
u(ϵ,η)=∫Rf(ϵ+t(η−ϵ))dt
其中,
ϵ
,
η
\epsilon, \eta
ϵ,η分别是射线源坐标与探测器单元坐标。
[1] John's equation. https://en.wikipedia.org/wiki/John%27s_equation
[2] X-ray transform. https://en.wikipedia.org/wiki/X-ray_transform
[3] John’s Equation-based Consistency Condition and Corrupted Projection Restoration in Circular Trajectory Cone Beam CT
超平面
超平面是n维欧氏空间中,余维数(codimension)为1的子空间。即超平面是n维空间中的n-1维的子空间。
超平面是平面中的直线、空间中的平面的推广。
若
W
W
W是向量空间
V
V
V的一个子空间,
W
W
W在
V
V
V的余维数是商空间
V
/
W
V/W
V/W的维数:
c
o
d
i
m
(
W
)
=
d
i
m
(
V
/
W
)
codim(W) = dim(V/W)
codim(W)=dim(V/W)
若
V
V
V是有限维的,有:
c
o
d
i
m
(
W
)
=
d
i
m
(
V
/
W
)
=
d
i
m
(
V
)
−
d
i
m
(
W
)
codim(W) = dim(V/W) = dim(V) - dim(W)
codim(W)=dim(V/W)=dim(V)−dim(W)
令
F
F
F为域(如:
R
R
R),则n维空间
F
n
F^n
Fn中的超平面为:
Σ
i
=
1
n
a
i
x
i
=
b
\Sigma_{i=1}^n a_i x_i = b
Σi=1naixi=b
其中,系数
a
i
a_i
ai不全为0
[1] 余维数. https://zh.wikipedia.org/wiki/餘維數
[2] 超平面. https://zh.wikipedia.org/wiki/超平面
[3] 域. https://www.bananaspace.org/wiki/域
[4] John's equation. https://en.wikipedia.org/wiki/John%27s_equation
[5] X-ray transform. https://en.wikipedia.org/wiki/X-ray_transform
Radon变换
Radon变换的复数版为Penrose变换;
Radon变换与Fourier变换的关系:中心切片定理;
Radon变换
R
\mathcal{R}
R的对偶变换为在超平面上反投影
R
∗
\mathcal{R}^*
R∗,
Radon变换的逆变换为DHB算法(求导、希尔伯特变换、反投影)。
Laplace算子
Δ
\Delta
Δ与
∂
s
2
\partial^2_s
∂s2都具有旋转不变性,可以得到:
R
(
Δ
f
)
=
∂
s
2
(
R
f
)
\mathcal{R} (\Delta f) = \partial^2_s (Rf)
R(Δf)=∂s2(Rf)
R
∗
(
∂
s
2
g
)
=
Δ
(
R
∗
g
)
\mathcal{R}^* (\partial^2_s g) = \Delta (\mathcal{R}^* g)
R∗(∂s2g)=Δ(R∗g)
这个性质叫Intertwining property
[1] What does "regularity condition" mean. https://stats.stackexchange.com/questions/654975/what-does-regularity-condition-mean
[2] Radon transform. https://en.wikipedia.org/wiki/Radon_transform
[3] Invariant differential operator. https://en.wikipedia.org/wiki/Invariant_differential_operator
[4] 部分经典反演的推导(附矩阵形式). https://www.cnblogs.com/HaHeHyt/p/18236125
斜坡滤波器的拆分
q
(
s
)
=
p
(
s
)
∗
h
(
s
)
q(s)=p(s) * h(s)
q(s)=p(s)∗h(s)
q
(
s
)
=
d
p
(
s
)
d
s
∗
−
1
s
q(s) = \frac{dp(s)}{ds} * \frac{-1}{s}
q(s)=dsdp(s)∗s−1
Q
(
w
)
=
P
(
w
)
×
H
(
w
)
Q(w) = P(w) \times H(w)
Q(w)=P(w)×H(w)
H
(
w
)
=
∣
w
∣
=
i
w
×
s
g
n
(
w
)
−
i
H(w) = |w| = iw \times \frac{sgn(w)}{-i}
H(w)=∣w∣=iw×−isgn(w)
其中,
s
s
s为探测器坐标,
w
w
w为角频率
由傅里叶微分性质可得:在频域与
i
w
iw
iw相乘,相当于在空间域微分;
由符号函数的傅里叶变换以及对偶性可得:在频域与
s
g
n
(
w
)
−
i
\frac{sgn(w)}{-i}
−isgn(w)相乘,相当于在时域与
1
/
s
1/s
1/s卷积(希尔伯特变换)。
傅里叶变换
傅里叶变换:
X
(
j
w
)
=
∫
−
∞
∞
x
(
t
)
e
−
j
w
t
d
t
X(jw) = \int_{-\infty}^{\infty}x(t) e^{-jwt}dt
X(jw)=∫−∞∞x(t)e−jwtdt
逆变换:
x
(
t
)
=
1
2
π
∫
−
∞
∞
X
(
j
w
)
e
j
w
t
d
w
x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(jw) e^{jwt} dw
x(t)=2π1∫−∞∞X(jw)ejwtdw
微分性质:
若
x
(
t
)
↔
X
(
w
)
x(t) \leftrightarrow X(w)
x(t)↔X(w),则
d
x
(
t
)
d
t
↔
i
w
X
(
w
)
\frac{dx(t)}{dt} \leftrightarrow iwX(w)
dtdx(t)↔iwX(w)
积分性质:
若
x
(
t
)
↔
X
(
w
)
x(t) \leftrightarrow X(w)
x(t)↔X(w),则
∫
−
∞
t
x
(
τ
)
d
τ
↔
1
i
w
X
(
w
)
+
π
X
(
0
)
δ
(
w
)
\int_{-\infty}^t x(\tau) d\tau \leftrightarrow \frac{1}{iw}X(w)+\pi X(0) \delta(w)
∫−∞tx(τ)dτ↔iw1X(w)+πX(0)δ(w)
对偶性:
若
x
(
t
)
↔
X
(
w
)
x(t) \leftrightarrow X(w)
x(t)↔X(w),则
X
(
t
)
↔
x
(
−
w
)
X(t) \leftrightarrow x(-w)
X(t)↔x(−w)
[1] The Step Function and the Signum Function. https://www.thefouriertransform.com/pairs/step.php
[2] Cauchy principal value. https://en.wikipedia.org/wiki/Cauchy_principal_value
[3] Sign function. https://en.wikipedia.org/wiki/Sign_function
[4] 符号函数的傅里叶变换. https://www.tutorialspoint.com/fourier-transform-of-signum-function
不同的书里,傅里叶变换写法不同,有的会写
2
π
2\pi
2π
符号函数的定义也不同,有的定义
s
g
n
(
0
)
=
0
sgn(0)=0
sgn(0)=0,有的定义
s
g
n
(
0
)
=
1
sgn(0)=1
sgn(0)=1
不过这个不影响对理论的理解
ROI重建
内部重建问题(Interior Problem):探测器不够大,无法把整个物体摄入其内,但可把 ROI 完全看到。
问题1:
R
O
I
<
F
O
V
ROI<FOV
ROI<FOV
可以通过DBH算法精确重建,求导与反投影是局部运算,尽管希尔伯特变换无法应用于截断数据,但可以用有限的希尔伯特反变换。
问题2:
R
O
I
=
F
O
V
ROI=FOV
ROI=FOV,探测器一端在物体外面
利用复变函数中的解析开拓理论求解,但需要事先已知
F
O
V
FOV
FOV中的一个区域(可以很小)。
如果复变函数 h ( z ) h(z) h(z)在区域内解析(即在区域内每一点都无穷阶可导),只要知道很小的子区域中 h ( z ) h(z) h(z)的值,便可唯一地解出该区域内其它点的值。
解析函数的定义:在定义域内每一点的泰勒级数皆逐点收敛的光滑函数。
解析函数具有唯一性。
解析延拓是将解析函数的定义域扩大的方法。
解析函数
f
(
z
)
,
z
∈
D
f(z), z\in D
f(z),z∈D,如果存在解析函数
F
(
z
)
,
z
∈
G
,
G
⊂
D
F(z), z\in G, G\subset D
F(z),z∈G,G⊂D,使得
F
(
z
)
=
f
(
z
)
,
∀
z
∈
D
F(z)=f(z), \forall z\in D
F(z)=f(z),∀z∈D,则称
F
(
z
)
F(z)
F(z)是
f
(
z
)
f(z)
f(z)从
D
D
D到
G
G
G的解析延拓。
举例
物体横截面
f
(
x
,
y
)
f(x,y)
f(x,y),半径为
r
r
r的圆区域为已知,半径为1的圆区域为FOV。
以求解
f
(
x
,
0
)
f(x,0)
f(x,0)为例,即过横截面的直线
y
=
0
y=0
y=0为例,其余直线同理。
将
y
=
0
y=0
y=0代入Radon逆变换得到:
f
(
x
,
0
)
=
1
4
π
∫
0
2
π
d
θ
∫
−
∞
∞
1
s
−
x
c
o
s
θ
∂
p
(
s
,
θ
)
∂
s
d
s
f(x,0)=\frac{1}{4\pi} \int_0^{2\pi} d\theta \int_{-\infty}^\infty \frac{1}{s-xcos\theta} \frac{\partial p(s,\theta)}{\partial s} ds
f(x,0)=4π1∫02πdθ∫−∞∞s−xcosθ1∂s∂p(s,θ)ds
将该式拆分为两项,一项为
s
s
s在
[
−
1
,
1
]
[-1,1]
[−1,1]积分,一项为
s
s
s在
[
−
1
,
1
]
[-1,1]
[−1,1]之外的区域积分
第一项可用测得到数据计算出,第二项并未测得,无法计算,将其记为
h
(
x
)
h(x)
h(x)
h
(
x
)
=
1
4
π
∫
0
2
π
d
θ
∫
∣
s
∣
>
1
1
s
−
x
c
o
s
θ
∂
p
(
s
,
θ
)
∂
s
d
s
,
∣
x
∣
<
1
h(x)=\frac{1}{4\pi} \int_0^{2\pi} d\theta \int_{|s|>1} \frac{1}{s-xcos\theta} \frac{\partial p(s,\theta)}{\partial s} ds, |x|<1
h(x)=4π1∫02πdθ∫∣s∣>1s−xcosθ1∂s∂p(s,θ)ds,∣x∣<1
由于
∣
s
∣
>
1
,
∣
x
∣
<
1
|s|>1, |x|<1
∣s∣>1,∣x∣<1,使得
1
/
(
s
−
x
c
o
s
θ
)
1/(s-xcos\theta)
1/(s−xcosθ)无穷阶可导
假设
∂
p
(
s
,
θ
)
∂
s
\frac{\partial p(s,\theta)}{\partial s}
∂s∂p(s,θ)连续
此时,用复变量
z
=
x
+
i
y
z=x+iy
z=x+iy替换实变量
x
x
x,
h
(
z
)
,
∣
z
∣
<
1
h(z), |z|<1
h(z),∣z∣<1就是一个复解析函数
h
(
x
)
h(x)
h(x)又可写作:
h
(
x
)
=
f
(
x
,
0
)
−
1
4
π
∫
0
2
π
d
θ
∫
∣
s
∣
≤
1
1
s
−
x
c
o
s
θ
∂
p
(
s
,
θ
)
∂
s
d
s
h(x)=f(x,0) - \frac{1}{4\pi} \int_0^{2\pi} d\theta \int_{|s|\leq 1} \frac{1}{s-xcos\theta} \frac{\partial p(s,\theta)}{\partial s} ds
h(x)=f(x,0)−4π1∫02πdθ∫∣s∣≤1s−xcosθ1∂s∂p(s,θ)ds
第一项中,
f
(
x
,
0
)
f(x,0)
f(x,0)在
∣
x
∣
≤
r
|x|\leq r
∣x∣≤r是已知的
第二项是可计算出的
利用解析开拓,
h
(
z
)
h(z)
h(z)在区域
∣
z
∣
<
1
|z|<1
∣z∣<1的值就确定了,从而
h
(
x
)
h(x)
h(x)在区域
∣
x
∣
<
1
|x|<1
∣x∣<1也确定了。从而可以计算出
f
(
x
,
0
)
,
∣
x
∣
<
1
f(x,0), |x|<1
f(x,0),∣x∣<1
[1] 闭合形式的解. https://mathworld.wolfram.com/Closed-FormSolution.html
[2] J. Zhao, Z. Chen, L. Zhang and D. Wu, "CT reconstruction method from truncated projection based on FSM analytic continuation," 2015 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), San Diego, CA, USA, 2015, pp. 1-4, doi: 10.1109/NSSMIC.2015.7582222.
Revision of Complex Analysis; Analytic Continuation; Residues
复变函数 f ( z ) = u ( z ) + i v ( z ) , z = x + i y f(z) = u(z)+iv(z), z=x+iy f(z)=u(z)+iv(z),z=x+iy
f
(
z
)
f(z)
f(z)可导,则可推出Cauchy-Riemann条件:
u
x
=
v
y
,
u
x
=
−
u
y
u_x = v_y, u_x=-u_y
ux=vy,ux=−uy
f
(
z
)
,
z
∈
D
f(z), z\in D
f(z),z∈D解析是指:
f
′
(
z
)
f'(z)
f′(z)在整个定义域
D
D
D上都存在。
解析函数也称全纯函数、正则函数(regular)。
Cauchy-Goursat定理:
如果
f
(
z
)
f(z)
f(z)解析,则其在闭曲线上的积分为0,即
∮
f
(
z
)
d
z
=
0
\oint f(z) dz=0
∮f(z)dz=0
可得出:
∫
a
b
f
(
z
)
d
z
\int_a^b f(z)dz
∫abf(z)dz与路径无关
Cauchy积分公式:
如果
f
(
z
)
f(z)
f(z)解析,则
f
(
z
0
)
=
1
2
π
i
∮
f
(
z
)
/
(
z
−
z
0
)
d
z
=
0
f(z_0) = \frac{1}{2\pi i} \oint f(z)/(z-z_0) dz=0
f(z0)=2πi1∮f(z)/(z−z0)dz=0
f
(
n
)
(
z
0
)
=
n
!
2
π
i
∮
f
(
z
)
/
(
z
−
z
0
)
n
+
1
d
z
=
0
f^{(n)}(z_0) = \frac{n!}{2\pi i} \oint f(z)/(z-z_0)^{n+1} dz=0
f(n)(z0)=2πin!∮f(z)/(z−z0)n+1dz=0
f
(
z
)
f(z)
f(z)在曲线
C
C
C包围的区域内解析,则其所有阶导数也在该区域解析。
整函数(entire):复平面上除
∞
\infty
∞外,处处解析的函数
Liouville定理:有界整函数是常函数
Taylor定理:解析函数有一个关于 z z z收敛的展开, z z z在其定义域内。
孤立奇点(isolated singularity):若
f
(
z
)
f(z)
f(z)在
z
0
z_0
z0不解析,但在
z
0
z_0
z0的某一去心邻域内解析,则称
z
0
z_0
z0是
f
(
z
)
f(z)
f(z)的孤立奇点。
(邻域的缩写是nbhd?)
Laurent级数:如果
f
(
z
)
f(z)
f(z)的孤立奇点为
z
0
z_0
z0,则有
f
(
z
)
=
σ
n
=
0
∞
a
n
(
z
−
z
0
)
n
+
b
1
/
(
z
−
z
0
)
+
b
2
/
(
z
−
z
0
)
2
+
.
.
.
f(z) = \sigma_{n=0}^\infty a_n (z-z_0)^n + b_1/(z-z_0) + b_2/(z-z_0)^2+...
f(z)=σn=0∞an(z−z0)n+b1/(z−z0)+b2/(z−z0)2+...
根据系数
b
n
b_n
bn可将孤立奇点分类:
simple pole
b
n
=
0
,
f
o
r
n
>
1
b_n = 0, for\space n > 1
bn=0,for n>1
pole of order N
b
n
=
0
,
f
o
r
n
>
N
b_n = 0, for\space n > N
bn=0,for n>N
essential singularity 无穷多幂级数
1
/
(
z
−
z
0
)
1/(z-z_0)
1/(z−z0)
将复变量
z
z
z写作极坐标形式
z
=
R
e
i
θ
z=Re^{i\theta}
z=Reiθ,得到:
l
n
z
=
l
n
R
+
i
θ
lnz = lnR + i\theta
lnz=lnR+iθ
由
z
=
R
e
i
θ
=
R
e
i
θ
+
i
2
π
n
z=Re^{i\theta} = Re^{i\theta+i2\pi n}
z=Reiθ=Reiθ+i2πn,得到:
l
n
z
=
l
n
R
+
i
θ
+
i
2
π
n
lnz = lnR + i\theta + i2\pi n
lnz=lnR+iθ+i2πn
所以一个
z
z
z对应多个值。
恒等定理(Identity theorem):如果两个函数在同一区域内解析,对同一区域内的曲线积分值相同,则两个函数在该区域内处处相等。
Abel定理:幂级数在
x
0
x_0
x0收敛,则满足
∣
x
∣
<
∣
x
0
∣
|x|<|x_0|
∣x∣<∣x0∣的任意
x
x
x处,幂级数收敛;幂级数在
x
0
x_0
x0发散,则满足
∣
x
∣
>
∣
x
0
∣
|x| >|x_0|
∣x∣>∣x0∣的任意
x
x
x处,幂级数发散。
幂级数在
x
0
x_0
x0收敛,则
l
i
m
n
→
∞
a
n
x
0
n
=
0
lim_{n \rightarrow \infty} a_n x_0^n = 0
limn→∞anx0n=0
Abel定理表明一定存在一个
R
R
R,使幂级数的收敛区间为
(
−
R
,
R
)
(-R,R)
(−R,R),收敛半径为
R
R
R
解析延拓:
给定一个以
z
1
z_1
z1为中心点且具有有限收敛半径的幂级数,则幂级数至少表示了一个在其收敛域内解析的函数
f
1
(
z
)
f_1(z)
f1(z)。
将
f
1
(
z
)
f_1(z)
f1(z)在点
z
2
z_2
z2展开为新的级数,则其收敛域可能超过原收敛域。
f
2
(
z
)
f_2(z)
f2(z)是
f
1
(
z
)
f_1(z)
f1(z)在新区域的解析延拓。
举例:
幂级数
Σ
n
=
0
∞
z
n
\Sigma_{n=0}^\infty z^n
Σn=0∞zn在
∣
z
∣
<
1
|z|<1
∣z∣<1内收敛,代表一个解析函数
f
1
(
z
)
f_1(z)
f1(z);在
∣
z
∣
>
1
|z|>1
∣z∣>1内发散。
f
1
(
z
)
f_1(z)
f1(z)在对原收敛域内任一点(如
i
/
2
i/2
i/2)的函数值以及各阶导数,从而构造出泰勒级数:
Σ
n
=
0
∞
f
1
(
n
)
(
i
/
2
)
(
z
−
i
/
2
)
n
\Sigma_{n=0}^\infty f_1^{(n)}(i/2)(z-i/2)^n
Σn=0∞f1(n)(i/2)(z−i/2)n,该级数在
∣
z
−
i
/
2
∣
<
5
/
2
|z-i/2|<\sqrt{5}/2
∣z−i/2∣<5/2收敛,代表一个解析函数
f
2
(
z
)
f_2(z)
f2(z)
f
1
(
z
)
=
1
/
(
1
−
z
)
,
g
1
:
∣
z
∣
<
1
f_1(z) = 1/(1-z), g_1: |z|<1
f1(z)=1/(1−z),g1:∣z∣<1
f
2
(
z
)
=
1
/
(
1
−
z
)
,
g
2
:
∣
z
−
i
/
2
∣
<
5
/
2
f_2(z) = 1/(1-z), g_2: |z-i/2|<\sqrt{5}/2
f2(z)=1/(1−z),g2:∣z−i/2∣<5/2
f 1 ( z ) f_1(z) f1(z)是 f 2 ( z ) f_2(z) f2(z)在 g 1 g_1 g1内的解析延拓, f 2 ( z ) f_2(z) f2(z)是 f 1 ( z ) f_1(z) f1(z)在 g 2 g_2 g2内的解析延拓。
[1] Revision of Complex Analysis; Analytic Continuation; Residues. https://www2.ph.ed.ac.uk/~mevans/amm/lecture02.pdf
[2] 复变函数的柯西-黎曼条件 - 光能丰的文章 - 知乎. https://zhuanlan.zhihu.com/p/21827856
[3] 复变多值函数的黎曼面 (Riemann surface)、分支点 (branch point) 与割线 (branch cut) - 東雲正樹的文章 - 知乎. https://zhuanlan.zhihu.com/p/422338793
[4] 复旦大学-幂级数. https://math.fudan.edu.cn/_upload/article/files/8a/c7/b6bc60b64933acee041a351adda9/66ceeb73-64e3-4a1a-b59c-a341e9914deb.pdf
解析延拓还得学一下,推荐课程:北京大学 数学物理方法 吴崇试
希尔伯特变换
两次希尔伯特变换:
H
(
H
f
)
=
−
f
H(Hf) = -f
H(Hf)=−f
实值函数与其希尔伯特变换是正交的
F
(
G
f
)
=
−
i
⋅
s
g
n
(
w
)
⋅
(
F
f
)
F(Gf) = -i \cdot sgn(w) \cdot (Ff)
F(Gf)=−i⋅sgn(w)⋅(Ff)
希尔伯特变换是一个全通滤波器,相位移动
±
9
0
o
\pm 90^o
±90o
DBH算法中,求导、希尔伯特变换、反投影可以互相交换位置
有限希尔伯特变换
假定 f ( t ) = 0 , ∣ t ∣ ≥ 1 ; f ( t ) ≠ 0 , ∣ t ∣ < 1 f(t)=0, |t|\geq1; f(t)\neq 0, |t|<1 f(t)=0,∣t∣≥1;f(t)=0,∣t∣<1
g
=
H
f
=
∫
−
1
1
f
(
τ
)
1
π
(
t
−
τ
)
d
τ
g = Hf = \int_{-1}^1 f(\tau) \frac{1}{\pi (t-\tau)} d\tau
g=Hf=∫−11f(τ)π(t−τ)1dτ
g
g
g未必只在有限区间内非零
而下面两种有限希尔伯特反变换仅用到
g
g
g在
[
−
1
,
1
]
[-1,1]
[−1,1]区间的值:
f
(
t
)
=
1
π
t
−
1
t
+
1
∫
−
1
1
τ
+
1
τ
−
1
g
(
τ
)
t
−
τ
d
τ
f(t) = \frac{1}{\pi} \sqrt \frac{t-1}{t+1} \int_{-1}^1 \sqrt \frac{\tau+1}{\tau-1} \frac{g(\tau)}{t-\tau} d\tau
f(t)=π1t+1t−1∫−11τ−1τ+1t−τg(τ)dτ
f ( t ) = 1 π 1 − t 2 ∫ − 1 1 f ( τ ) d τ + 1 π ∫ − 1 1 1 − τ 2 τ − t g ( τ ) d τ f(t) = \frac{1}{\pi \sqrt{1-t^2}} \int_{-1}^1 f(\tau)d\tau + \frac{1}{\pi} \int_{-1}^1 \frac{\sqrt{1-\tau^2}}{\tau - t} g(\tau) d\tau f(t)=π1−t21∫−11f(τ)dτ+π1∫−11τ−t1−τ2g(τ)dτ