图像重建流程

Gate仿真

得到Root数据，用Analysis.C将数据转化成listmode数据

将listmode数据转化为3D的sinogram数据

将3Dsinogram用SSRB转化成2D的slice切片sinogram进行重建

	

本文作者：飞利浦资深CT专家，Colin Han，韩太林。感谢CT Clinical Evidence授权转载。
相对于MR，CT在参数及扫描方面并不是太难，但是CT的图像重建及各种算法则是非常难的，也是比较抽象的。这篇文章介绍CT图像重建算法等相关内容。
CT技术是CT诊断的基础，帮助医务工作者充分掌握CT技术是我们的责任和义务！
CT基础知识
Image Reconstruction Method
图像重建算法
在上一期内容CT原理1中,我们主要介绍了X线与物质的作用,也就是X线的衰减是如何发生的,正是由于这种衰减的存在,X线才可以被用于CT成像,那么探测器获得的衰减信号最终是如何被转换成CT图像的呢?这就涉及到图像的重建算法了,所以今天我们接着X线的衰减,继续介绍CT图像重建的原理。
1
重建算法的分类
CT重建算法共有3类,如下图。
(1)反投影法
(2)迭代重建算法
(3)解析法:包括滤波反投影法和傅里叶变换法
其中(3)解析法中的A滤波反投影法是在(1)反投影法的基础上发展起来的,通过加入滤波函数解决了图像锐利度的问题,如下图,没有加入滤波函数(Filter Function)时重建的图像是模糊的(红色),而加入滤波函数(Filter Function)后重建的图像是清晰的(黄色),由于这种方法(FBP)重建速度快,且图像质量比较高,成为了最常用的CT图像重建方法;而(2)迭代重建算法虽然很早就提出了(早于FBP法),但是由于其计算量非常大,重建速度慢,有赖于计算机性能的突破,所以这种方法的普及反而在FBP算法之后。
2
反投影法原理
反投影法是最简单的一种CT图像重建方法
CT图像重建的基础：相同的X射线强度，穿过不同物质，发生的衰减不同，利用这一规律可以区别人体内的不同物质。
CT图像重建的过程: 简单来说,CT扫描中X线穿透人体的每个层面内的结构可以被分成多个小立方体(被称为体素:Voxel),每个小立方体都对应一个单独的衰减信号,把这个信号输入到图像平面矩阵中相应的小格子(被称为像素:Pixel)中,把每一个体素的衰减信号都输入到相应的像素内,然后以不同的灰度反映出来,这就是CT图像重建的过程。
下面以4个体素为例对反投影法的原理进行简单说明:
假设四个体素的衰减系数分别为1、2、3、4,那么直接把相应体素的衰减系数直接输入到对应的像素内就可以生成图像了。是不是觉得太简单了呢?
但实际上,每个体素具体的衰减系数我们是不知道,通过扫描获得的是不同方向上,各体素衰减系数的总和,如下图,投影1-6(方向为从球管到探测器)获得的总衰减系数分别为:7、3、5、4、6、5,而每个体素具体的衰减系数是多少呢?接下来,我们利用已知的数据进行反投影(方向为从探测器到球管):
根据探测器获得的数据反过来对该方向上的每个体素进行相同的赋值,称之为反投影。
方向1反投影结果为:
方向2反投影结果为:
方向3反投影结果为:
方向4反投影结果为:
方向5反投影结果为:
方向6反投影结果为:
把这6次反投影的结果叠加到一起为:
叠加的过程导致每个像素的所得数据为3次自身衰减系数+另外三个体素的衰减系数(三个方向的三次反投影),所以需要减去一个相同的基数(四个体素的衰减系数之和),并除以2,才是确切的衰减系数.
由投影1+投影2或者投影4+投影5可得出4个体素总的衰减系数为10,那么减去基数后结果为:
然后再除以2结果为:
最终获得的这个像素矩阵就是对上述四个体素的准确还原.而这种方法就是CT图像重建中反投影法的基本原理,由于这个过程不需要对每一个体素的具体衰减系数进行求解,而是直接进行反投影获得叠加数据,然后再对获得叠加数据进行线性的调整就能实现对每个体素的还原,因此重建速度非常快,但是这种方法要求各个方向的投影数据必须是完备的,而且噪声与辐射剂量呈线性关系。
3
迭代重建算法原理
迭代重建算法是计算量非常大的一种CT图像重建方法
简单来说就是利用求解线性方程组的方式来重建图像，实际上是从—幅假设的初始图像出发，采用逐步逼近的方法，将理论投影值同实际测量投影值不断进行比较和迭代更新,直到最终获得最优解,最早被用在PET成像中,如下图[1]。
下面还是以4体素为例对迭代重建算法进行简单说明:
假设扫描获得的投影值结果为:
先对图像进行猜测,如每个像素都为10:
然后对这个猜测结果获得的投影值与实际扫描获得的投影值进行比较,获得校准系数(为实际测量的投影值与猜测所得的投影值之比),如下图:
然后根据本次校准系数进行一次迭代:
然后又获得新的校准系数,并进行二次迭代:
然后又获得了新的校准系数,并进行三次迭代:
......
以此类推,不断迭代,直到最终校准后的各体素衰减值对应的投影值与实际扫描获得的投影值差距为最小(Error=min),则满足了迭代终止条件,然后才根据最后校准获得的体素衰减值来重建图像.这种方法就是CT图像重建中迭代重建法的基本原理,由于这种方法需要对每一个体素所对应的具体衰减系数进行求解,而医学图像所用的矩阵至少为512x512,因此每层图像需要求解的未知数为262144个,所以这种方法的计算量非常之大,故重建速度慢,但是这种方法的噪声不再受到辐射剂量的限制,能够在辐射剂量较低的条件下，获得高质量图像。
重点回顾
CT图像重建的过程就是把每一个体素的衰减系数对应到每一个像素中去,反投影法不需要求解也能进行还原,重建速度快,而迭代重建需要先求得最优解,然后再进行图像重建,速度慢但在低剂量时也可获得高质量图像。
注意事项！
✔ 相同的X射线强度，穿过不同物质，发生的衰减不同，利用这一规律可以区别人体内的不同物质。
✔ 反投影法不需要对每一个体素的具体衰减系数进行求解,因此重建速度非常快,但是这种方法要求各个方向的投影数据必须是完备的,而且噪声与辐射剂量呈线性关系.
✔ 迭代重建需要对每一个体素所对应的具体衰减系数进行求解,因此重建速度慢,但能够在辐射剂量较低的条件下，获得高质量图像。
[1].Knesaurek K, Machac J, Vallabhajosula S, Buchsbaum MS. A new iterative reconstruction technique for attenuation correction in high-resolution positron emission tomography. Eur J Nucl Med. 1996;23(6):656‐661. doi:10.1007/BF00834527.

	

英文名称为Algebraic reconstruction technique：即代数重建算法
代数重建技术(ART)是一种用于计算机层析成像的迭代重建技术。它从一系列的角度投影(正弦图)中重建一幅图像。Gordon、Bender和Herman首次证明了它在图像重建中的应用。而这种方法被称为数值线性代数中的Kaczmarz方法。相对于其他重建方法(如滤波反投影)，ART的一个优点是，将先验知识(已知的约束条件)纳入重建过程是相对容易的。该算法的实质是用迭代法求解线性方程组的解。
clc;clear all;close all;N = 180; %图像大小N2 = N^2;I = phantom(N);% 产生头模型图像theta = linspace(0,180,181);theta = theta(1:180);%% = = = = = = 产生投影数据 = = = = = = %%P_num = 260; %探测器通道个数P = ParallelBeam(theta ,N ,P_num); %产生投影数据%P = radon(I,theta);%% = = = = = = 获取投影矩阵 = = = = = =%%delta = 1;% 网格大小，角度增量[W_ind,W_dat] = SystemMatrix(theta,N,P_num,delta);%% = = = = = = 进行ART迭代 = = = = = = %%F = zeros(N2,1);  %初始图像向量lambda = 0.25;  %松弛因子c = 0;  %迭代计数器irt_num = 5;  while(c    for j = 1:length(theta)        for i = 1:1:P_num            % 取得一条射线所穿过的网格编号和长度            u = W_ind((j-1)*P_num + i,:);  % 编号            v = W_dat((j-1)*P_num + i,:);  % 长度                        if any(u) == 0                continue;            end            %恢复投影矩阵中与这一条射线对应的行向量 w            w = zeros(1,N2);            ind = u > 0;            w(u(ind))=v(ind);            % 图像进行一次迭代            PP = w * F;   %前向投影            C = (P(i,j)-PP)/sum(w.^2) * w';  % 修正项                       F = F + lambda * C;        end    end    F(F<0) = 0;   % 小于0的像素值置为0    c = c+1;endF = reshape(F,N,N)'; % 转换成N x N的图像矩阵%% = = = = = = 仿真结果显示 = = = = = = %%figure(1);imshow(I);xlabel('(a)180x180头模型图像');figure(2);%A = imadjust(F);imshow(F,[]);xlabel('(b)ART算法重建的图像');     
往期回顾>>>>>>
【模式识别】Matlab指纹识别【优化求解】A*算法解决三维路径规划问题  matlab自动识别银行卡号【优化问题】基于栅格地图——遗传算法的机器人最优路径规划【优化求解】模拟退火遗传实现带时间窗的车辆路径规划问题【数学建模】Matlab实现SEIR模型分享到朋友圈获取6个赞即可获取完整代码