CT图像重建计算机模拟实验研究

２０１４年第１５期总第１５９期　Ｓ＿ＬＩＣｏＮ　ＶＡＬＬＥＹ　ＣＴ图像重建计算机模拟实验研究　秦闻聆，何治艳，戴榕　（海南医学院，海南海口５７１　１９９）　摘　要　ｃＴ图像重建理论是医学影像技术的基础，为了帮助医学影像学专业的学生更好地理解ｃＴ图像重建的思想，　本文研究了基于ＭＡＴＬＡＢ和ｃ＃的ｃＴ图像重建计算机模拟软件的制作。利用ＭＡＴＬＡＢ进行编程，实现动态显示扫描过程　及断层图像的重建。利用ｃ＃进行软件界面设计，处理好ＭＡＴＬＡＢ与ｃ＃的接口，最终形成的软件可以满足医学院校影　像学专业进行ｃＴ图像重建计算机模拟实验的要求。　关键词ｃＴ图像重建；ＭＡＴＬＡＢ；ｃ＃；混合编程　中图分类号：ＴＰ３９１．４１　文献标识码：Ａ　文章编号：１６７１—７５９７（２０１４）１　５－００４５—０４　ｃＴ图像重建的思想和理论是医学影像物理学的重要内　容，涉及较多的数学知识，内容抽象难学，教学难度大，需　２　ＭＡＴＬＡＢ程序优化　虽然ＭＡＴＬＡＢ软件提供了大量的专业化的工具箱，但是　ＭＡＴＬＡＢ语言是一种解释性语言，其执行速度较慢，当用户选　择不同的重建算法对头模型进行图像重建时，程序运行速度较　慢。主要是因为在编程时并没有考虑到程序优化。ＭＡＴＬＡＢ中通　过ｔｉｃ函数和ｔｏｃ函数来查看一段程序的ＣＰＵ使用时间，比较　程序优化前后的运行速度。　要开发直观形象的教学软件辅助教学。Ｍｉｃｒｏｓｏｆｔ　Ｖｉｓｕａｌ　Ｃ＃　是Ｍｉｃｒｏｓｏｆｔ专门为使用．ＮＥＴ平台开发的一种强大的、面向　组件的语言，可用于方便快捷的创建运行在．ＮＥＴ公共语言运　行库（ｃｏｍｍｏｎ　ｌａｎｇｕａｇｅ　ｒｕｎｔｉｍｅ，ＣＬＲ）上的Ｗｉｎｄｏｗｓ应用程　序Ⅲ。ｃ＃是微软．ＮＥＴ战略中核心的开发工具，它综合了　Ｖｉｓｕａｌ　Ｂａｓｉｃ的高效率和ｃ＋＋功能的强大性，具有良好的界面　设计功能，可以很方便的建立应用程序的可视化界面　。但其　计算能力与ＭＡＴＬＡＢ相比编程显得较为复杂。ＭＡＴＬＡＢ是一种面　向科学计算、数据可视化、系统仿真及交互式程序设计的高级　语言　Ｊ。它以矩阵运算为基础，极少的代码即可实现复杂的功　能。将ｃ＃界面设计能力和ＭＡＴＬＡＢ的强大科学计算能力相结合，　可以解决图像处理，工程计算等多个领域，具有良好的应用前景。　基于ＭＡＴＬＡＢ和ｃ＃的ＣＴ图像重建计算机模拟软件可运用于医学　院校医学影像学专业的教学，学生可以通过本软件掌握ｃＴ图像　重建的物理原理和数学方法。　２．１使用ＭＡＴＬＡＢ自带函数　ＭＡＴＬＡＢ内部函数是由更底层的编程语言ｃ构造的，它的执　行速度显然快于自己编写的程序　。这些函数都是经过严格测　试的，运行不会出现错误。这样既可以避免自己编写函数繁琐　的过程也可以加快运行速度。如在取整数的时候就可以直接使　用ＭＡＴＬＡＢ自带的ｆｉｘ函数、ｒｏｕｎｄ函数。　２．２大型矩阵的预先确定维数　大型矩阵动态扩展是非常耗费时间的，因为预先确定数组　维数后不允许数组动态扩展，虽然数组的适应性降低，但其访　１软件设计与实现　ｃＴ图像重建是通过ｃ＃与ＭＡＴＬＡＢ混合编程实现的。软件的　设计目的是用两种ｘ射线束（平行柬和扇形束）对不同头模型　（椭圆头模型和ｓ—Ｌ头模型）进行扫描，用户可以选择不同的重　建算法（直接反投影法，Ｒ—Ｌ滤波反投影法，ｓ—Ｌ滤波反投影法）　对头模型进行图像重建。此外用户还可以直观看到不同方法的　扫描过程。软件功能图如图ｌ所示。　问速度加快，访问变得非常容易，这样可以大大减少运行时间，　所以用ＭＡＴＬＡＢ的内部函数确定数组维数，然后再进行运算。比　如数值数组用ｚｅｒｏｓ０，结构数组用ｒｅｐｍａｔＯ等，用ｒｅｐｍａｔ０　函数能获得连续的内存块，加快的程序运行速度　。　２．３自定义函数　在编译代码时会遇到多次使用到相同语句的时候，为了减　少因为重复运行相同代码所消耗的时间，可以把程序中需要用　到的语句写成一个函数，在之后碰到相同的情况下就可以直接　调用已经定义好的函数，这样既缩短了ＭＡＴＬＡＢ代码的长度，也　加快了程序运行速度。例如程序会多次用到将反投影运算，为　了避免代码重复，就可以将反投影运算单独写成一个函数，之　后再遇到反投影运算时就可以直接调用该函数，加快代码运行　速度。　２．４循环向量化法　为了避免因多次使用ｗｈｉｌｅ和ｆｏｒ而减慢程序运行速度，　可以将带有ｗｈｉｌｅ和ｆｏｒ循环的语句转化为向量操作。在很多　情况下，优先考虑利用循环过程向量化来提高ＭＡＴＬＡＢ代码的运　行速度，这是ＭＡＴＬＡＢ程序优化非常有效的方法。循环向量化解　决重复执行动作的思想是通过ｇ￣ＴＬＡＢ的矩阵操作完成的。向量　是矩阵的一种特殊形式，如果运算对象是向量，那么运算结果　也是向量，这就可以省去这其中循环结构的使用，大大提高运　行速度。对于不能向量化的程序，应该尽量使内循环的次数多　Ｙ　Ｒ＆Ｄ　于外循环的次数　。　例１以椭圆头模型在平行束扫描下的直接反投影法为例，　程序优化前后对比，具体代码如下：　程序优化前：　％Ｎ２为投影次数，Ｎｌ为ｘ（Ｙ）方向上的像素数。　ｆｕｎｃｔｉｏｎ　ｚｊ（Ｎ２，Ｎ１）　ｔｉＣ　ｄｅｔａ＝４０／Ｎ１：　％每个方向取值间隔　Ｎ３＝ｒｏｕｎｄ（４０＊ｓｑｒｔ（２）／ｄｅｔａ）＋２；　％每一次投影所得到的　取样数　Ｉｏ＝ｉ．０：　％１ｏ为椭圆内部的密度值　ｄ＝ｐｉ／Ｎ２：　％Ｒ方向的角度变化　Ａ＝１０：　％椭圆长半轴　Ｂ＝６：　％椭圆短半轴　ｆｏｒ　ｊ＝ｌ：Ｎ３　％取样点变化　ｆｏｒ　ｉ＝ｌ：Ｎ２　％Ｒ方向的角度变化　ｒｒ＝Ａ　２＊ｃｏｓ（ｉ＊ｄ）＊ＣＯＳ（ｉ＊ｄ）＋Ｂ　２＊ｓｉｎ（ｉ＊ｄ）＊ｓｉｎ（ｉ　ｄ）；％ｒｒ代表ｒ　２　Ｒ＝ｊ＊ｄｅｔａ—Ｎ３／２＊ｄｅｔａ：　％将Ｒ轴的坐标原点移　到最左边　ｉｆ　ａｂｓ（Ｒ）＜＝ｓｑｒｔ（ｒｒ）　ｚ（ｉ，Ｊ）＝１ｏ　２　Ａ　Ｂ　ｓｑｒｔ（ｒｒ—Ｒ　２）／ｒｒ／８：％将　得到的投影数据存储为矩阵　ｅ１ｓｅ　ｚ（ｉ，ｊ）＝Ｏ：　ｅｎｄ　ｅｎｄ　ｅｎｄ　．ｆ＝ｚｅｒｏｓ（，　Ｎ１，Ｎ１）：　％定义一个ｎ＊ｎ矩阵，初始化为零　矩阵　扎一　一毗一．，　～．，　一一～　…　．，＇蔷．　．，　一一一　一一　一ｆｏｒ　ａ＝ｌ：Ｎｌ　ｆｏｒ　ｂ＝ｌ：Ｎ１　ｆｏｒ　ｃ＝１：Ｎ２　Ｒ０＝（ａ一１咖～　）＊ＣＯＳ（ｃ＊ｄ）＋（ｂ一１）＊ｓｉｎ（ｃ＊ｄ）一（Ｎｌ．　．　１）／２＊（ＣＯＳ（ｃ＊ｄ）＋ｓｉｎ（ｃ＊ｄ））：％任意一个像素（ａ，ｂ）的Ｒ坐标　ＲＩ＝Ｎ３／２＋ＲＯ：　ｎＯ＝ｆｉｘ（Ｒ１）：％整数部分　ｇ＝Ｒ１一ｎＯ：　％小数部分　ｆ（ａ，ｂ）＝ｆ（ａ，ｂ）＋（１一ｇ）　ｚ（Ｃ，ｎＯ）　ｇ　ｚ（Ｃ　ｎ０＋１）：　ｅｎｄ　ｅｎｄ　ｅｎｄ　ｆ＝ｆ／Ｎ２：　ｉｍｓｈｏｗ（ｆ）　程序优化后：　％平行束Ｘ射线扫描获得椭圆头模型的投影数据　ｆｕｎｃｔｉｏｎ　ｐｖ＝ｒａｄｏｎｅｌ（Ｎ２，Ｎ３，ｄｅｔａ）　Ａ＝１０：％椭圆长半轴　Ｂ＝６：％椭圆短半轴　ａ＝ｚｅｒｏｓ（Ｎ２，Ｎ３）：　ｂ＝ｚｅｒｏｓ（Ｎ２，Ｎ３）：　ｐｖ＝ｚｅｒｏｓ（Ｎ２，Ｎ３）：　ｍ＝ｌ：Ｎ２：　ａ＝Ａ　２：Ｉｃ（ＣＯＳ（ｍ＊ｐｉ／Ｎ２））．　２＋Ｂ　２　（ｓｉｎ（ｍ＊ｐｉ／Ｎ２））．　２：　～枷一　　一～咖一一．，　心吼　一　吼叭．莩　狲　一州一～　唧　一　一　，一　　一一～一　一　　～圳Ｉ至．　圳．，．，　．，　．，　值　福一缉．　，　ｌＯ，２０，２０］，’ＦａｃｅＣｏｌｏｒ’，’ｋ’）　ｒｅｃｔａｎｇｌｅ（’Ｐｏｓｉｔｉｏｎ’，［一３，一５，６，ｉ０］，’Ｃｕｒｖａｔｕｒｅ’，［Ｉ，　１］，’ＦａｃｅＣｏｌｏｒ’，’ｗ’）　ａｘｉＳ　ｏｆｆ：　ｄａｓｐｅｃｔ（［１，１，１］）　ｇ＝ｒａｄｏｎｅｌ（Ｎ２，Ｎ３，ｄｅｔａ）：　ｆ＝ｉｒａｄｏｎｈｙ（Ｎ１，Ｎ２，Ｎ３，ｇ）：　ｓｕｂｐｌｏｔ（１，２，２），ｉｍｓｈｏｗ（ｆ／ｌＯ），ｔｉｔｌｅ（’重建后的图’）：　ｔｏｃ　若Ｎ２＝１００：ＮＩ＝１００，则程序优化前的输出结果为：Ｅｌａｐｓｅｄ　ｔｉｍｅ　ｉｓ　０．５８３４４３　ｓｅｃｏｎｄｓ，则程序优化后的输出结果为：　Ｅｌａｐｓｅｄ　ｔｉｍｅ　ｉｓ　０．２４３６０８　ｓｅｃｏｎｄｓ。程序优化前后对比，通　　一一　　一．，ｈ川～　．．　一．二