表格 4.2 不同精度的重构结果
original Raum Sektor
2倍精度 Sektor
4倍精度
SART
原始图 Raum Sektor
2倍精度 Sektor
4倍精度
SMART
MLEM
FDK
5 重构方法的评估
这一章将实现章节2.2中介绍的变化重构法。目前有两种标准：“最大准确性（maximale Dichtegenauigkeit）”和“梯度（Gradient）”，来对重构后的体素进行评估。好的体素用λj = 0标记，坏的体素用λj = 1标记。
5.1 梯度（Gradient）
根据经验，被检测物体边界上的体素，重构质量总是很差，尤其是对于小角度范围的圆形轨迹。所以要首先把这些边界提取出来，再次重构，来提高整体的重构质量。在本范文中对边界的提取是利用“Sobel-Operator”[11]。“Sobel-Operator”是在图像处理中经常用到的边界识别过滤器。对于二维的图像有如下结构，分别对x-和y-方向进行处理。
(5.1)
x-1, y-1 x, y-1 x+1, y-1
x-1, y x,y x+1, y
x-1, y+1 x, y+1 x+1, y+1
图 5.1 3×3的二维图像
图5.1描述了二维图像中任意一个像素（x,y）和它的8个相邻像素。像素（x,y）的明暗值用B（x,y） 表示。利用公式5.1能计算出像素（x,y）的水平和竖直的梯度。
于是像素（x,y）的梯度模为：
(5.2)
一个三维重构体积中的体素拥有在x-,y-和z-方向上三个离散的梯度，为此需要参照“Sobel-Operator”开发一个三维的边界识别过滤器，如图5.2。
图5.2 三维边界提取结构
图5.3 3×3×3 的三维体积
以图5.3的中心（x,y,z）为参照点，对应的x-方向上的梯度为：
(5.3)
对应的y-，z-方向上的梯度为：
(5.4)
(5.5)
同样，体素（x，y，z） 的梯度模为：
(5.6)
表格5.1是“空心球”Raum和Sektor的重构图，对其运用公式5.3到5.6，得到的沿着箭头所示的梯度值曲线描绘在图5.4中。
表格 5.1 “空心球”重构图
原始图 Raum Sektor
2倍精度
图 5.4 “空心球”的梯度值曲线
从图5.4可以看出，边界上的梯度值远大于其它的区域，也就是说被检测物体的边界被正确的识别出来。具体在编程中还需要定义一个门槛值0.3，所有体素的梯度值小于0.3，用λj = 0 标记，其余的用λj = 1 标记。被λj = 1标记的体素将再次做重构。表格5.2描绘了相应的四种算法的梯度图。
表格 5.2 “空心球”的梯度图
重构
(Raum) 重构
(Sektor 2倍精度) 梯度
(Raum) 二进制化
(Raum) 梯度 (Sektor) 二进制化
(Sektor)
SART
(门槛值 = 0,3)
SMART
(门槛值 = 0,4)
MLEM
(门槛值 = 0,4)
FDK
(门槛值 = 13)
下面将对另一个检测物体“Small Piece”进行模拟，表格5.3是对应的梯度图。
表格 5.3 “Small Piece”的梯度图
重构
(Raum) 重构
(Sektor 2倍精度) 梯度
(Raum) 二进制化
(Raum) 梯度 (Sektor)
二进制化
(Sektor)
SART
(门槛值 = 0,06)
SMART
(门槛值 = 0,05)
MLEM
(门槛值 = 0,05)
FDK
(门槛值 = 2,6)
5.2 最大准确性（maximale Dichtegenauigkeit）
评估标准“梯度”不足之处在于，所求的某个体素的梯度值只与周围相邻的26个体素有关。因此有可能把一个体素评估为“好”，当这个体素周围相邻的体素都为“坏”。这个问题可以通过第二种评估标准“最大准确性（maximale Dichtegenauigkeit）”得以解决。
“最大准确性（maximale Dichtegenauigkeit）”主要是基于体素值偏差的计算，即公式5.7。如果重构后的体素值越接近理想值，这个偏差(fid)就会越小[9]。再通过公式5.8的计算，就能得出最大准确性(gjd)，实现对重构质量的评价。当重构的体素质量越差，其体素值偏差就越大，对应的最大准确性的值也就越小。
理想值vdef_d下第j个体素值的偏差： (5.7)
理想值vdef_d下第j个体素值的最大准确性: (5.8)
N = 被重构物体的体素个数
P = 投影像素的个数
D = 重构体积中不同材料的个数
vj = 第j 个体素值
vdef_d = 第d种材料的理想体素值
pi = 在投影中第i个像素的明暗值
wijF = 第j个体素在第i条正向投影光线下的加权系数
wijB = 第j个体素在第i条逆向投影光线下的加权系数
接下来就是编程的问题，如何将“最大正确性”这一评估标准用“OpenGL”在显卡中实现。首先必须定义每种材料的理想体素值，通常人们是取这种材料重构后的平均体素值作为它的理想体素值。对于图5.5中的例子，它有4种不同材料组成（包含空气）。而对于每种材料的最大正确性，是呈高斯分布的。
图 5.5 体素j的偏差和最大正确性
如果要实现不同精度的最大正确性评估，需要将公式5.7做如下变化，原理详见第四章。
(5.9)
但是上述公式还有所缺陷。当实际测得的投影值Pi与重构后Raum和Sektor的投影值Σwinvn之间的偏差过大时，会导致重构体素值vj与之前定义的材料体素值vdef_d之间的偏差对结果起不了作用。为此在实践中要增加一个附加的限制条件：
(5.10)
其中 是每两种定义的材料体素值的平均差。 例如，图5.5中定义的四种材料体素值分别为： vdef_1 = 0； vdef_2 = 0.009； vdef_3 = 0.018 和 vdef_4 = 0.027， 于是 = 0.009.
公式5.9 的编程与重构很相似。更确切地说只用到投影和逆投影，这些步骤在 [1] 中都有详细的介绍。公式5.10左边部分其实就是第四章中不同精度下的校正这一步骤。至于算数运算加，减，乘，除以及平方可以用“OpenGL Shading Language”很容易实现。加权系数wij的确定是通过一个常量“1”代替“vn”向接收器做投影。图5.6评估标准“最大正确性”的程序设计结构图。
图 5.6 评估标准“最大正确性”的程序设计结构图
表格 5.4 “空心球”最大准确性评估结果
重构
(Raum) 重构
(Sektor 2倍精度) 最大准确性 (Raum) 最大准确性
(Sektor)
SART
表格5.4是对“空心球”做最大准确性评估后的结果。从中可以看出，被检测物体的边界能很好的被提取出来。最大准确性的值在边界上是最小的。图5.7描绘了沿着表格5.4中箭头方向的最大准确性曲线。
图 5.7 “空心球”最大准确性曲线-SART
评估标准是用来鉴别“坏”的重构体素。因此还需定义一个门槛值，例如图5.7中值为5000的这条临界线。所有的体素，当它们的最大准确性的数值小于5000，用λj = 1 标记，其余的用λj = 0 标记。被λj = 1标记的体素将再次做重构。
表格 5.5 “空心球”的最大准确性评估结果
重构
(Raum) 重构
(Sektor
2倍精度) 最大准确性
(Raum) 二进制化
(Raum) 最大准确性
(Sektor) 二进制化
(Sektor)
SART
(门槛值 = 5000)
SMART
(门槛值 = 2000)
MLEM
(门槛值 = 2000)
FDK
(门槛值 = 1,5)
从表格5.5看出，对FDK-算法的评估效果不是很好。沿着箭头方向将FDK的最大准确性曲线描绘在图5.8中，很难找到一条分界线，来区分重构体素的“好”与“坏”。
图 5.8 “空心球”最大准确性曲线-FDK
原因是FDK采用的是分析的方法，被重构的体素不是像重复叠代法那样逐次接近理想值，而是所有的X光照片一次性向重构体积做逆投影，所以相对来说不是很稳定，容易产生大的偏差。FDK-算法根本上也不适用于变化重构法。
下面将对另一个检测物体“Small Piece”进行模拟，结果显示在下表中。
表格 5.6 “Small Piece”的最大准确性评估结果
重构
(Raum) 重构
(Sektor
2倍精度) 最大准确性
(Raum) 二进制化
(Raum) 最大准确性
(Sektor) 二进制化
(Sektor)
SART
(门槛值 = 45000)
SMART
(门槛值 = 35000)
MLEM
(门槛值 = 35000)
FDK
(门槛值 = 6)
5.3 两种评估标准的结合
两种评估标准“梯度”和“最大准确性”应该一起使用，以更好地对重构体素进行鉴定。值得一提的是，在评估过程之前应该先用传统重构法做几次重构（一般为3次），待重构体素值相对稳定后再做评估，效果比较好。下图是完整的程序流程图，包含两种评估标准和重构体素值的标准化。
图 5.9 完整程序流程图
6 比较和评价
6.1 变化重构法的计算量和内存消耗
为了比较变化重构法的计算量和内存消耗，将对7个不同精度和不同角度圆形轨迹的物体做三维重构。首先在表格6.1中对每个被检测的物体进行字母排列，在图6.1中对程序的各部分进行编号，以便能将比较结果清楚地展现出来。
表格 6.1 被检测物体字母排列
体积精度 接收器精度 投影张数 材料种数 字母排列
64×64×64 64×64 64 2 A
128×128×128 128×128 54 2 B
125 2 C
168×128 164 2 D
507×512 235 2 E
256×256×256 168×128 164 2 F
507×512 235 2 G
表格 6.2 变化重构法计算量和内存消耗的比较
Prüfkörper aus Tabelle 6.1
A B C D E F G
图6.1中每个步骤的计算量 (秒) 1 1,4 6,4 22,1 15 34,5 232,5 243,9
2 0,5 1 2,3 29 63,1 131 213,7
3 0,03 0,06 0,125 0,204 3,07 0,187 4,07
4 1,2 3 6,6 60,2 163,3 271,2 304,6
5 0,03 0,03 0,032 0,032 0,032 0,032 0,032
6 0,085 0,78 0,781 0,64 0,672 6,6 5,5
7 3,4 11,8 31,2 105,5 265,3 647,1 779,2
内存消耗 (MB) 73 180 193 203 889 967 1664
从上表中可以发现,与传统重构法相比,变化重构法的计算量非常大(消耗的时间多)。其中评估过程（“2”到“6”）占了很大的比重。在评估过程中，大部分运行时间在“2”和“4”上，用来计算体素值偏差（公式5.7）。随着材料种数的增加，步骤6的运行时间也会相应增加。
6.2 变化重构法的评价
本节中将对两个检测物体分别做小角度圆形轨迹的传统重构和变化重构。Sektor的精度是Raum的两倍。每种方法的测试结果都列在下面的表格中，从中可以很容易的发现变化重构法的优势或者改进。
表格 6.3 “Platine”146°圆形轨迹重构
重构次数 3 重构次数 6 重构次数 9 重构次数 20
传统 SART Raum-
重构
Sektor-
重构
变化 SART Raum-
重构
Raum-
“坏”的体素
Sektor-
重构
Sektor-
“坏”的体素
表格 6.4 “Small Piece”135°圆形轨迹重构-SART
重构次数 3 重构次数 6 重构次数 9 重构次数 20
传统 SART 重构
变化 SART 重构
“坏”的体素
重构
“坏”的体素
表格6.3和6.4中的传统重构图比较模糊，主要原因是圆形轨迹的角度范围小，所有X光照片上的信息量不足以把被检测的物体三维的重构出来。利用第五章中介绍的评估标准可以把“坏”的体素识别出来，并通过再次重构，减少“坏”的体素的个数，提高整体重构质量。比如，不管是Raum还是Sektor，经过20次变化重构后，表面显得很均匀，边界也很清晰。这也是为什么要引进变化重构法的原因。
下面是用变化重构法对另两种算法（SMART和MLEM）的测试，通过对比也能得到与SART相同的结论。
表格 6.5 “Small Piece”135°圆形轨迹重构-SMART
重构次数 3 重构次数 6 重构次数 9 重构次数 20
传统 SMART 重构
变化 SMART 重构
“坏”的体素
重构
“坏”的体素
表格 6.6 “Small Piece”135°圆形轨迹重构-MLEM
重构次数 3 重构次数 6 重构次数 9 重构次数 20
传统 MLEM 重构
变化 MLEM 重构
“坏”的体素
重构
“坏”的体素
7 总结
本范文主要是对三维重构检测做了两方面的研发。
重构体素值标准化的引入，能使其与实际的物理体素值相匹配。
在实际应用中需要被检测的物体用变化重构法实现不同精度的重构。为此要在Raum中另外定义一个精度高的Sektor, 然后两者同步执行重构的四个步骤。根据变化重构法的原理，只有被鉴别为“坏”的体素才能再次重构。用来鉴别的标准有两种：“梯度”和“最大准确性”。变化重构法的优势在于，当没有足够的X光照片来清晰的重构检测的物体时，采用变化重构法连续做多次（一般20次）的重构，依然能使重构图接近原始的物体。但是这种方法需要大量的运算，比较费时。
参考文献
[1] Wei Chen, GPU-basierte Beschleunigung der dreidimensionalen Rekonstruktionstechnik SART, Masterarbeit am Institut für Innovations-Transfer an der Fachhochschule Hannover, 2007
[2] Wei Chen, GPU-basierte Beschleunigung der dreidimensionalen Rekonstruktionstechnik SART, Ergänzung der Masterarbeit am Institut für Innovationstransfer an der Fachhochschule Hannover, 2007
[3] Aninash C. Kak, Malcolm Slaney, Principles of Computerized Tomographic Imaging, pp. 285-292, ISBN 089871494X, 2001
[4] Cristian Badea, Richard Gordon, Experiments with the nonlinear and chaotic behavior of the multiplicative algebraic reconstruction technique (MART) algorithm for computed tomography IOP Publishing Ltd, 2004
[5] Georg Steidl, Bildrekonstruktion für PET/SPECT-System, Universität Erlangen, 2006
[6] L. Feldkamp, L. Davis, and J. Kress, Practical cone beam algorithm, J. Opt. Soc. Am., pp. 612-619, 1984
[7] Stefan Speer, Simulation von PET-Daten und Rekonstruktion unter Berücksichtigung der Dämpfungskarte, 2004
[8] Fastest Fourier Transformation in the West, FFTW, hhtp://www.fftw.org, 23.8.2009
[9] A. Frost, Iterative Rekonstruktionsverfahren, inklusive Algebraische Rekon- struktionstechnik ART unter Berücksichtigung von a-priori Information, Institut für Innovations-Transfer an der Fachhochschule Hannover, 2007
[10] http://de.wikipedia.org/wiki/Sobel-Operator, 30.8.2009
[11] Wolfgang Abmayr, Einführung in die digitale Bildverarbeitung, Fachhochschule München, 1994
[12] Ning Ye, Anwendung der linearen Regression zur Schätzung von Materialdichten in Röntgenrekonstruktion, Praktikumsbericht am Institut für Innovations-Transfer an der Fachhochschule Hannover, 2009