CT重建笔记(五)—2D平行束投影公式
写的又回去了,因为我发现我理解不够透彻,反正想到啥写啥,尽量保证内容质量好简洁易懂
2D平行束投影公式
p ( s , θ ) = ∫ ∫ f ( x , y ) δ ( x c o s θ + y s i n θ − s ) d x d y p(s,\theta)=\int \int f(x,y)\delta(x cos\theta + ysin\theta - s) dxdy p(s,θ)=∫∫f(x,y)δ(xcosθ+ysinθ−s)dxdy式1
p ( s , θ ) = ∫ f ( s c o s θ − t s i n θ , s s i n θ + t c o s θ ) d t p(s,\theta)=\int f(scos\theta - tsin\theta,s sin\theta + t cos\theta)dt p(s,θ)=∫f(scosθ−tsinθ,ssinθ+tcosθ)dt式2
p ( s , θ ) = ∫ f ( s θ ⃗ + t θ ⃗ ) d t p(s,\theta) = \int f(s\vec \theta + t \vec \theta)dt p(s,θ)=∫f(sθ+tθ)dt式3
p ( s , θ ) = ∫ f θ ( s , t ) d t p(s,\theta)=\int f_\theta(s,t)dt p(s,θ)=∫fθ(s,t)dt式4
这些公式来自医学图像重建1.5节
,这些都是2D平行束投影公式(我偷懒没写积分的上下限)。
s s s表示探测器像素单元的坐标, θ \theta θ表示投影角度。
f ( x , y ) f(x,y) f(x,y)为物体截面的线衰减系数。
x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s是平面直线方程的一种写法(Hesse Normal Form,见https://leslielee.blog.csdn.net/article/details/145670396
),表示投影角度为 θ \theta θ投影至探测器坐标 s s s处的射线。
式1
δ ( z ) \delta(z) δ(z)是一个广义函数,具有筛选的性质(见https://leslielee.blog.csdn.net/article/details/144859730
,但由于我的数学水平不够,这篇写的比较差劲)
给定 s , θ s,\theta s,θ后,若固定 y y y,则式1
中位于里面的积分可写作:
∫ f ( x , y ˉ ) δ ( x c o s θ ˉ + y ˉ s i n θ ˉ − s ˉ ) d x = f ( x ~ , y ˉ ) \int f(x,\bar y)\delta(x cos \bar \theta + \bar ysin \bar \theta - \bar s) dx = f(\tilde x,\bar y) ∫f(x,yˉ)δ(xcosθˉ+yˉsinθˉ−sˉ)dx=f(x~,yˉ)
其中, x ~ c o s θ ˉ + y ˉ s i n θ ˉ = s ˉ \tilde x cos \bar \theta + \bar ysin \bar \theta = \bar s x~cosθˉ+yˉsinθˉ=sˉ,即 ( x ~ , y ˉ ) (\tilde x, \bar y) (x~,yˉ)在直线 x c o s θ ˉ + y s i n θ ˉ = s ˉ x cos \bar \theta + ysin \bar \theta = \bar s xcosθˉ+ysinθˉ=sˉ上。
给定 s , θ s,\theta s,θ后,遍历所有的 y y y(式1
中位于外面积分的作用),便可实现将位于直线 x c o s θ ˉ + y s i n θ ˉ = s ˉ x cos \bar \theta + ysin \bar \theta = \bar s xcosθˉ+ysinθˉ=sˉ上的 f ( x , y ) f(x,y) f(x,y)进行求和。
因此,式1
得到的是位于直线 x c o s θ + y s i n θ = s x cos \theta + ysin \theta = s xcosθ+ysinθ=s上的 f ( x , y ) f(x,y) f(x,y)的和。
若令 θ ⃗ = ( c o s θ , s i n θ ) \vec \theta = (cos \theta, sin \theta) θ=(cosθ,sinθ), x ⃗ = ( x , y ) \vec x = (x,y) x=(x,y),则式1
可得到向量写法:
p ( s , θ ) = ∫ ∫ f ( x , y ) δ ( x ⃗ ⋅ θ ⃗ − s ) d x d y p(s,\theta)=\int \int f(x,y)\delta(\vec x \cdot \vec \theta - s) dxdy p(s,θ)=∫∫f(x,y)δ(x⋅θ−s)dxdy
式2
( x , y ) (x,y) (x,y)要位于直线 x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s上才可计算线积分。因此, x = s c o s θ − t s i n θ , y = s s i n θ + t c o s θ x = s cos\theta - t sin\theta, y=s sin\theta + t cos \theta x=scosθ−tsinθ,y=ssinθ+tcosθ必然已经将 ( x , y ) (x,y) (x,y)约束至直线 x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s上。
定义 t t t,令 x = s c o s θ − t s i n θ , y = s s i n θ + t c o s θ x = s cos\theta - t sin\theta, y=s sin\theta + t cos \theta x=scosθ−tsinθ,y=ssinθ+tcosθ,将 x , y x,y x,y带入直线方程中得到:
( s c o s θ − t s i n θ ) c o s θ + ( s s i n θ + t c o s θ ) s i n θ = s (s cos\theta - t sin\theta) cos\theta + (s sin\theta + t cos \theta) sin\theta = s (scosθ−tsinθ)cosθ+(ssinθ+tcosθ)sinθ=s
进一步化简可得:
s = s s=s s=s
因此, x = s c o s θ − t s i n θ , y = s s i n θ + t c o s θ x = s cos\theta - t sin\theta, y=s sin\theta + t cos \theta x=scosθ−tsinθ,y=ssinθ+tcosθ 等效于 x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s,即将 ( x , y ) (x,y) (x,y)约束至直线 x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s上。
(那是如何才能想到这么表示 x , y x,y x,y呢,对此我咨询了元宝
,元宝
告我跟旋转矩阵相关。)
x = s c o s θ − t s i n θ , y = s s i n θ + t c o s θ x = s cos\theta - t sin\theta, y=s sin\theta + t cos \theta x=scosθ−tsinθ,y=ssinθ+tcosθ 可写成向量矩阵形式:
(见https://leslielee.blog.csdn.net/article/details/135566902
)
[ c o s θ − s i n θ s i n θ c o s θ ] [ s t ] = [ x y ] \begin{bmatrix} cos\theta & -sin\theta \\ sin\theta & cos\theta \end{bmatrix} \begin{bmatrix} s \\ t \end{bmatrix}=\begin{bmatrix} x \\ y \end{bmatrix} [cosθsinθ−sinθcosθ][st]=[xy]
x , y x,y x,y是 s , t s,t s,t逆时针旋转得到的。
那么, s , t s,t s,t是 x , y x,y x,y顺时针旋转得到的,可得到表达式:
[ c o s θ s i n θ − s i n θ c o s θ ] [ x y ] = [ s t ] \begin{bmatrix} cos\theta & sin\theta \\ -sin\theta & cos\theta \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}=\begin{bmatrix} s \\ t \end{bmatrix} [cosθ−sinθsinθcosθ][xy]=[st]
这样我们便可明白, t = − x s i n θ + y c o s θ t = -xsin\theta + ycos\theta t=−xsinθ+ycosθ,式2
给的 t t t并不是没有限制的。
直线 − x s i n θ + y c o s θ = t -xsin\theta + ycos\theta = t −xsinθ+ycosθ=t与 x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s垂直,因为 ( − s i n θ , c o s θ ) ⋅ ( c o s θ , s i n θ ) = 0 (-sin\theta,cos\theta)\cdot (cos\theta,sin\theta) = 0 (−sinθ,cosθ)⋅(cosθ,sinθ)=0。
给定 s , θ s,\theta s,θ后,遍历所有 t t t便可取到直线 x c o s θ + y s i n θ = s x cos\theta + ysin\theta = s xcosθ+ysinθ=s上的所有 ( x , y ) (x,y) (x,y)点。
式3
令 θ ⃗ ⊥ = ( − s i n θ , c o s θ ) \vec \theta^{\perp} = (-sin\theta,cos\theta) θ⊥=(−sinθ,cosθ)
将 s θ ⃗ + t θ ⃗ ⊥ s\vec \theta + t \vec \theta^{\perp} sθ+tθ⊥展开便得到: ( s c o s θ − t s i n θ , s s i n θ + t c o s θ ) (s cos\theta - t sin\theta, s sin\theta + t cos\theta) (scosθ−tsinθ,ssinθ+tcosθ)
而函数 f f f是一个二元函数,自变量是 x , y x,y x,y,因此有:
s c o s θ − t s i n θ = x s cos\theta - t sin\theta = x scosθ−tsinθ=x
s s i n θ + t c o s θ = y s sin\theta + t cos\theta = y ssinθ+tcosθ=y
这不就是表示 ( x , y ) (x,y) (x,y)是由 ( s , t ) (s,t) (s,t)逆时针旋转 θ \theta θ得到的。
因此,式3
是式2
的向量写法。
式4
令 f θ ( s , t ) = f ( s c o s θ − t s i n θ , s s i n θ + t c o s θ ) f_\theta (s,t) = f(scos\theta - tsin\theta,s sin\theta + t cos\theta) fθ(s,t)=f(scosθ−tsinθ,ssinθ+tcosθ),便得到式4
。
式1
表示射线源与探测器同时逆时针旋转(射线旋转),物体不动。式2
表示物体顺时针旋转,射线不动。