Review 1. What is projection, Radon tranform? 2. What is sinogram, and its usage? 3. Direct back projection reconstruction method. 4. Center slice the

Size: px

Start display at page:

Download "Review 1. What is projection, Radon tranform? 2. What is sinogram, and its usage? 3. Direct back projection reconstruction method. 4. Center slice the"

辙翁
5 years ago
Views:

1 Lecture VIII CT Image Reconstruction(Ⅱ) 1

2 Review 1. What is projection, Radon tranform? 2. What is sinogram, and its usage? 3. Direct back projection reconstruction method. 4. Center slice theorem. 5. Fourier transform reconstruction. 2

3 Problems 1. 正弦图的两个坐标分别是什么? 2. 一点在正弦图上对应什么? 3. 正弦图有什么应用? 4. 什么是中心切片定理? 5. 什么是 Radon 空间, 物体空间向 Radon 空间的变化是什么过程? 反之是什么过程? 6. Radon 变换函数的两个变量是? 7. 什么是傅里叶变换重建, 这种重建方法存在什么问题? 3

4 Contents 1. Filtered back projection 滤波反投影 2. Convolution back projection 卷积反投影 3. R-L & S-L filter 4. CT reconstruction algorithm summary 5. Reconstruction process 6. An example 4

5 Problem of direct back projection b θ ( x, y) = g θ ( R)δ ( x cosθ + ysinθ R)dR b θ x 0, y 0 5 ( ) = g θ R ( ) = g θ ( R 0 ) π ( x, y) = b θ ( x, y)dθ = g θ x 0 cosθ + y 0 sinθ f b dθ = 0 π 0 g θ 如果系统输入是位于原点的一个 δ 函数 g θ ( R) = δ ( R) h b π ( r) = dθ = 0 π δ ( R)δ r cos( θ ϕ) R dr = δ r cos( θ ϕ) dθ δ θ π 2 +ϕ r cos ( θ ϕ ) θ π dθ = 1 0 rδ f x ( ) = 0 n ( )δ ( x 0 cosθ + y 0 sinθ R)dR ( R)δ ( x cosθ + ysinθ R)dR ( ) f ( ) δ x x n x n f b ( x, y) = f ( x, y) 1 r

6 Idea of removing artifacts crosssectional image (original) scan & projection back projection 2D filtering output image crosssectional image (original) scan & projection 1D filtering back projection Output image 6

7 Filtered backprojection reconstruction 1. 待重建图像可以从其二维傅里叶变换中重建出来 : (, ) ( ω, ω ) ( + ω ) j2πω1x 2y f x y F e d d = 2. 频域中极坐标形式, 即令 ω ω 1 2 2π 2π jω( xcosθ+ ysinθ) f ( x, y) F( cos, sin ) e d d 0 0 = 2π (, ) = (, ) ω = ωcos θ, ω = ωsinθ ω θ ω θ ω ω θ 3. 由中心切片定理, F( ωcos θ, ωsin θ) = P( ω, θ) (, ) (, + ) ( ωθ, + π) = P( ωθ, ) ( + ) 2π jω xcosθ ysinθ f x y d P e d 0 0 θ ω θ ω ω π 2π jω( xcosθ + ysinθ ) d P e d 0 0 = θ ω θ ω ω + π 2π jω( xcosθ + ysinθ ) d P e d 0 0 θ ω θ π ω ω 又因为, P, 故 π 2π jω( xcosθ+ ysinθ) (, ) (, ) f x y dθ P ω θ ω e dω = 0 (1) (2) (3) (4)

8 Filtered backprojection recon. deduction ( ) = dθ P( ω,θ ) ω e 2π jω R dω f x, y π π 0 ( ) = dθ [ P ω,θ 0 ω e 2π jω R dω]δ (xcosθ + ysinθ R)dR (5) Ø 对任意角度下的投影函数作一维傅立叶变换 ; ω Ø 乘上一维权重因子, 修正一维傅立叶变换的结果 ; Ø 对加权修正后的结果作一维傅立叶逆变换 ; Ø 用修正后的投影函数作直接反投影 ; Ø 改变投影角度, 重复上述步骤直至完成所有角度下的反投影滤波反投影重建的理论推证! 8

9 Filtered backprojection recon. 9

10 Convolution backprojection recon. g θ ( R) *C( R) g θ ( R) BP Interpolation ( x y) f, 空间域频域 1D IFT 1 1{ θ ( )} ω = 1{ θ ( )} ω C( R) = F { ω } ( ) 1 { ( )} ( ) 1 1 ω * 1 { ω} ( )* ( ) F g R F g R g θ R = F F gθ R = gθ R F = gθ R C R Key point: 1. 在频域中所做的滤波运算可以等效地在时域中用卷积完成 2. 从本质上讲, 卷积反投影与滤波反投影是一样的不同的是滤波反投影是将投影函数 g θ ( R) 变换到频域中, 然后用滤波函数 ω 对变换后的函数做滤波, 再反变换到空间域中作为修正过的投影函数而卷积反投影法是直接将 g θ ( R) 在空间域中进行卷积运算, 然后对卷积后的结果进行反投影 10

11 R-L Filter R-L 滤波函数是 1971 年由印度数学家 G.N. Ramachandran 和 A.V. Lakeshmnarayanan. 提出的, 以后 H ω 表示 R-L 滤波函数 ; 用 R 表示对应的卷积函数 ( ) ( ) R L 1. 设计 R-L 滤波函数的基本出发点是认为实际的二维图像函数总有一个频率上限, 所以频域中的滤波函数 ω 可以乘以一个宽度为 ω0 ω0 的窗函数 Q( ω) h R L H R L ( ω) (, ) Q ( ω ) 1 ω < ω0 = 0 ω ω0 ω 0 +ω 0 频带有限的滤波函数 ω 式中, ω 0 为截止频率, 即对投影的傅里叶变换进行滤波的是 : H ( ω) ω Q( ω) R L = 11

12 R-L filter 2. 卷积函数可以写成 ω ( t) 0 = H R L (ω ) e j2πωt dω h RL sinc x = ω 0 = ω 0 2 sin ( π x) π x 2sinc( 2ω 0 t) sinc 2 ( ω 0 t) 12

13 R-L filter 3. 对于离散化的数字图像已经设定最高频率为, 根据 Nyquist 采样定理选择空间域中的采样率为 2ω 0, 或者采样间隔为 δ = 1/2ω0 令上式中 t = nδ = n/2ω0, 可得离散化的 R-L 卷积函数 h R L ( nδ ) ω 0 1 n = 0 2 4δ = 0 n = Even 1 n = Odd π n δ 如果探测器板 200mm 宽, 有 256 个探测器, 则 δ =200/256; 然后计算每个探测器对应的 h R-L, 再进行傅里叶变换计算 H R-L n=(i-256/2) i=

14 S-L filter function S-L 滤波函数是美国学者 L.A. Shepp 和 B.F. Logan 于 1974 年提出的 1. 在频域中不使用矩形函数去截断 ω, 而是采用平滑的窗函数去约束滤波例如, 可以令窗函数为 2. S-L 滤波函数表示为 H S L W ( ) ( ω) ω = sinc Q( ω) 2ω0 ω ω = ω sinc Q( ω) 2ω0 14

15 3. 对应的卷积函数为 S-L filter function h S L ( t) 0 = ( πω t) 2 ( ω t) 1 4ω 1 4ω tsin 2 2 π 同样假设采样间隔 δ = 1/2ω t = nδ = n/2ω 0 0带入上式, 可得离散化的 S-L 卷积函数 0 h S L ( nδ ) = πδ 2 ( 4n 1) Key point: 1. S-L 卷积函数的优点是图像分辨率高且平滑, 适于脑部扫描 2. S-L 卷积函数重建的图像在高频响应方面不如 R-L 卷积函数好, 因为 S-L 函数在高频段偏离了理想的滤波函数 ω 3. 滤波器为滤波核与窗函数的乘积, 窗函数可以为 Sinc 或 Hanning 等 15

16 Analysis reconstruction summary No Method Features Merit&Drawback Notes 1 傅里叶切片重建先傅里叶变换, 再在频缺点 : 域内插值, 然后 2D 逆傅 1. 频域内插值和坐标变换计里叶变换算量大 2. 高频信号缺失导致失真 2 直接反投影先反投影, 在空间域插值和坐标变换, 然后在频域修正, 最后 2D 傅里叶逆变换 3 滤波反投影先对一维投影函数做傅里叶变换, 进行滤波处理修正, 再反投影缺点 : 1. 两次二维傅里叶变换花费时间相当可观 2. 图像模糊, 需费时修正优点 : 1. 两次一维傅里叶变换, 速度快, 准确度高图像模糊, 商用 CT 不使用商用 CT 几乎不使用多数商用 CT 都采用该方法, 边扫描边运算 4 卷积反投影直接对投影函数进行修优点 : 避免了傅里叶变换正后反投影缺点 : 在空间域需要大量的积分 ( 加法和乘法运算 ) 5... 与滤波反投影类似 16

17 Interpolation Key point: 1. 对于反投影, 插值是必需的 2. 插值可以在被滤波后的投影采样上或者重建图像上 3. 像素驱动反投影节省计算量 4. 在反投影之前, 滤波投影常此采用高阶插值方案预插值到一个更精细的采样栅格中, 反投影时就可采用简单的最近邻或线性插值 17

18 Reconstruction process 在角度 Theta 下, 获得投影数据 p(t,theta) 对 p(t,theta) 补零, 得到 p (t,theta) 对 p (t,theta) 进行傅里叶变换, 得到 P (Omega,Theta) P (Omega,Theta) 乘以斜变滤波器 H R-L, 得到 G(Omega, Theta) 对 G(Omega, Theta) 进行反傅里叶变换获得滤波投影 g(t, Theta) 将 g(t, Theta) 反投影累加到图像函数 f(x,y) 上 Last time Success 18

19 Conclusion 1. Direct backprojection mathematic analysis? 2. Direct backprojection process 3. Filtered backprojection 4. R-L filter & S-L filter 5. Reconstruction process 19

20 Questions??? 20

21 Example 已知一个物体的投影数据 ( 即正弦图, 三个投影角度 0 o,45 o,90 o ), 卷积函数为 (k1, k2, k3), 卷积修正后 45 o 时的反投影数据为 A45-I45, 利用卷积反投影重建方法, 计算重建矩阵 RA-RI 提示 : 计算卷积时, 扩展区填零 21

ii

i ii iii iv Abstract This senior project is to use compute simulation to accomplish analysis and synthesis of Cam. The object of these focuses on three major partsthe first one is to establish the mathematical