本站原创摘 要 CT图像重建理论是医学影像技术的基础,为了帮助医学影像学专业的学生更好地理解CT图像重建的思想,本文研究了基于MATLAB和C#的CT图像重建计算机模拟软件的制作.利用MATLAB进行编程,实现动态显示扫描过程及断层图像的重建.利用C#进行软件界面设计,处理好MATLAB与C#的接口,最终形成的软件可以满足医学院校影像学专业进行CT图像重建计算机模拟实验的要求.
关 键 词CT图像重建；MATLAB；C#；混合编程
中图分类号：TP391.41 文献标识码：A 文章编号：1671-7597（2014）15-0045-04
CT图像重建的思想和理论是医学影像物理学的重要内容,涉及较多的数学知识,内容抽象难学,教学难度大,需要开发直观形象的教学软件辅助教学.Microsoft Visual C#是Microsoft专门为使用.NET平台开发的一种强大的、面向组件的语言,可用于方便快捷的创建运行在.NET公共语言运行库（mon language runtime,CLR）上的Windows应用程
序[1].C#是微软.NET战略中核心的开发工具,它综合了Visual Basic的高效率和C++功能的强大性,具有良好的界面设计功能,可以很方便的建立应用程序的可视化界面[2].但其计算能力与MATLAB相比编程显得较为复杂.MATLAB是一种面向科学计算、数据可视化、系统仿真及交互式程序设计的高级语言[3].它以矩阵运算为基础,极少的代码即可实现复杂的功能.将C#界面设计能力和MATLAB的强大科学计算能力相结合,可以解决图像处理,工程计算等多个领域,具有良好的应用前景.基于MATLAB和C#的CT图像重建计算机模拟软件可运用于医学院校医学影像学专业的教学,学生可以通过本软件掌握CT图像重建的物理原理和数学方法.
1.软件设计与实现
CT图像重建是通过C#与MATLAB混合编程实现的.软件的设计目的是用两种X射线束（平行束和扇形束）对不同头模型（椭圆头模型和S-L头模型）进行扫描,用户可以选择不同的重建算法（直接反投影法,R-L滤波反投影法,S-L滤波反投影法）对头模型进行图像重建.此外用户还可以直观看到不同方法的扫描过程.软件功能图如图1所示.
图1 软件功能图
2.MATLAB程序优化
虽然MATLAB软件提供了大量的专业化的工具箱,但是MATLAB语言是一种解释性语言,其执行速度较慢,当用户选择不同的重建算法对头模型进行图像重建时,程序运行速度较慢.主要是因为在编程时并没有考虑到程序优化.MATLAB中通过tic函数和toc函数来查看一段程序的CPU使用时间,比较程序优化前后的运行速度.
2.1 使用MATLAB自带函数
MATLAB内部函数是由更底层的编程语言C构造的,它的执行速度显然快于自己编写的程序[4].这些函数都是经过严格测试的,运行不会出现错误.这样既可以避免自己编写函数繁琐的过程也可以加快运行速度.如在取整数的时候就可以直接使用MATLAB自带的fix函数、round函数.
2.2 大型矩阵的预先确定维数
大型矩阵动态扩展是非常耗费时间的,因为预先确定数组维数后不允许数组动态扩展,虽然数组的适应性降低,但其访问速度加快,访问变得非常容易,这样可以大大减少运行时间,所以用MATLAB的内部函数确定数组维数,然后再进行运算.比如数值数组用zeros（）,结构数组用repmat（）等,用repmat（）函数能获得连续的内存块,加快的程序运行速度[5].
2.3 自定义函数
在编译代码时会遇到多次使用到相同语句的时候,为了减少因为重复运行相同代码所消耗的时间,可以把程序中需要用到的语句写成一个函数,在之后碰到相同的情况下就可以直接调用已经定义好的函数,这样既缩短了MATLAB代码的长度,也加快了程序运行速度.例如程序会多次用到将反投影运算,为了避免代码重复,就可以将反投影运算单独写成一个函数,之后再遇到反投影运算时就可以直接调用该函数,加快代码运行速度.
2.4 循环向量化法
为了避免因多次使用while和for而减慢程序运行速度,可以将带有while和for循环的语句转化为向量操作.在很多情况下,优先考虑利用循环过程向量化来提高MATLAB代码的运行速度,这是MATLAB程序优化非常有效的方法.循环向量化解决重复执行动作的思想是通过MATLAB的矩阵操作完成的.向量是矩阵的一种特殊形式,如果运算对象是向量,那么运算结果也是向量,这就可以省去这其中循环结构的使用,大大提高运行速度.对于不能向量化的程序,应该尽量使内循环的次数多于外循环的次数[6].
例1 以椭圆头模型在平行束扫描下的直接反投影法为例,程序优化前后对比,具体代码如下：
程序优化前：
%N2为投影次数,N1为X（Y）方向上的像素数.
function zj（N2,N1）
tic
deta等于40/N1； %每个方向取值间隔
N3等于round（40*sqrt（2）/deta）+2； %每一次投影所得到的取样数
lo等于1.0； %lo为椭圆内部的密度值
d等于pi/N2； %R方向的角度变化
A等于10； %椭圆长半轴
B等于6； %椭圆短半轴
for j等于1：N3 %取样点变化
for i等于1：N2 %R方向的角度变化
rr等于A^2*cos（i*d）*cos（i*d）+B^2*sin（i*d）*sin（i*d）； %rr代表r^2
R等于j*deta-N3/2*deta； %将R轴的坐标原点移到最左边 if abs（R）<=sqrt（rr）
z（i,j）等于lo*2*A*B*sqrt（rr-R^2）/rr/8； %将得到的投影数据存储为矩阵
else
z（i,j）等于0；
end
end
end
f等于zeros（N1,N1）； %定义一个n*n矩阵,初始化为零
矩阵
for a等于1：N1
for b等于1：N1
for c等于1：N2
R0等于（a-1）*cos（c*d）+（b-1）*sin（c*d）-（N1-1）/2*（cos（c*d）+sin（c*d））；%任意一个像素（a,b）的R坐标
R1等于N3/2+R0；
n0等于fix（R1）； %整数部分
g等于R1-n0； %小数部分
f（a,b）等于f（a,b）+（1-g）*z（c,n0）+g*z（c,
n0+1）；
end
end
end
f等于f/N2；
imshow（f）
程序优化后：
%平行束X射线扫描获得椭圆头模型的投影数据
function pv等于radonel（N2,N3,deta）
A等于10；%椭圆长半轴
B等于6；%椭圆短半轴
a等于zeros（N2,N3）；
b等于zeros（N2,N3）；
pv等于zeros（N2,N3）；
m等于1：N2；
a等于A^2*（cos（m*pi/N2））.^2+B^2*（sin（m*pi/N2））.^2；
a等于repmat（a',1,N3）；
j等于1：N3；
b等于（（j-N3/2）*deta）.^2；
b等于repmat（b,N2,1）；
pv等于sqrt（（a-b））；
pv等于real（pv）；
pv等于2*A*B*pv./a；
%平行束扫描椭圆头模型投影值叠加
function ir等于iradonhy（N1,N2,N3,g）
a等于zeros（N1,N1）；
b等于zeros（N1,N1）；
R等于zeros（N1,N1）；
g1等于zeros（N1,N1）；
g2等于zeros（N1,N1）；
ir等于zeros（N1,N1）；
i等于1：N1；
j等于1：N1；
for m等于1：N2
a等于（i-1）*cos（m*pi/N2）；
a等于repmat（a',1,N1）；
b等于（j-1）*sin（m*pi/N2）；
b等于repmat（b,N1,1）；
c等于N3/2-（N1-1）*（cos（m*pi/N2）+sin（m*pi/N2））/2；
R等于a+b+c；
inte等于fix（R）；
dec等于R-inte；
g1等于g（m,inte）；
g1等于reshape（g1,N1,N1）；
g2等于g（m,inte+1）；
g2等于reshape（g2,N1,N1）；
ir等于ir+（1-dec）.*g1+dec.*g2；
end
ir等于ir/N2；
% 平行束X射线扫描椭圆头模型的直接反投影法重建图像
function TYZJP（N2,N1）
tic
deta等于40/N1；%为每个方向取值间隔
N3等于round（N1*sqrt（2）+2）； %每一次投影所得到的取样数
subplot（1,2,1,'replace'）；
xlim（[-10,10]）；
ylim（[-10,10]）；
title（'原图'）；
rectangle（'Position',[-10,-10,20,20],'FaceColor','k'）
rectangle（'Position',[-3,-5,6,10],'Curvature',[1,
1],'FaceColor','w'）
axis off；
daspect（[1,1,1]）
g等于radonel（N2,N3,deta）；
f等于iradonhy（N1,N2,N3,g）；
subplot（1,2,2）,imshow（f/10）,title（'重建后的图'）；
toc
若N2等于100；N1等于100,则程序优化前的输出结果为：Elapsed time is 0.583443 seconds,则程序优化后的输出结果为：Elapsed time is 0.243608 seconds.程序优化前后对比,通过将程序循环向量化,预先定义数组维数,将相同程序函数化和使用MATLAB自带的函数等优化方法有助于提高程序执行
效率.
3.C#和MATLAB混合编程
3.1 MATLAB编译环境的设置
为了能将MATLAB中的.m函数便以为动态链接库dll,对MATLAB的编译环境设置如下[7]： 在MATLAB的命令窗口（Command Window）键入“mex -setup”命令,并根据MATLAB的提示安装MATLAB编译器,过程如下：
>> mex -setup
Please choose your piler for building external interface （MEX） files：
Would you like mex to locate installed pilers [y]/n? Y
Select a piler：
[1] Digital Visual Fortran version 6.0 in C：＼Program Files＼Microsoft Visual Studio
[2] Lcc C version 2.4 in D：＼INSTALLPACKAGE＼MATLAB7＼sys＼lcc
[3] Microsoft Visual C/C++ version 6.0 in C：＼Program Files＼Microsoft Visual Studio
[0] None
Compiler： 3
Please verify your choices：
Compiler： Microsoft Visual C/C++ 6.0
Location： C：＼Program Files＼Microsoft Visual Studio
Are these correct?（[y]/n）： y
在MATLAB的命令窗口 （Command Window） 键入“mbuild -setup”命令后,并根据MATLAB的提示安装MATLAB编译器,过程如下：
>> mbuild -setup
Please choose your piler for building standalone MATLAB applications：
Would you like mbuild to locate installed pilers [y]/n? y
Select a piler：
[1] Lcc C version 2.4 in D：＼INSTALLPACKAGE＼MATLAB7＼sys＼lcc
[2] Microsoft Visual C/C++ version 6.0 in C：＼Program Files＼Microsoft Visual Studio
[0] None
Compiler： 2
Please verify your choices：
Compiler： Microsoft Visual C/C++ 6.0
Location： C：＼Program Files＼Microsoft Visual Studio
Are these correct?（[y]/n）： y
3.2 MATLAB创建COM组件
1）创建组件.
运行MATLAB 7.0,单击左下角的“start”,选择“MATLAB”,选择“MATLAB Builder for COM”.
单击File->New Project,输入信息,“Component name”项中填写组件名称“cton”,这时候会自动生成类“ctonclass”,在“Class name”项中填写类名称“f”,为了便于区分,将“Classes”中的“ctonclass”移除,再点击“Add>>”添加新类f,选择工程保存目录“Project directory”,在“Compiler options”中勾选“Create a singleton MCR”以及“Build dedug version：”点击OK,如果选择的工程所在目录并不存在,tool会提示该目录不存在,并询问是否创建该目录,选择Yes,tool就会创建工程到该目录.C成功创建COM组件[8].
2）添加m文件到COM组件.
在上图中文件夹f下的M-files,点击Add File按钮,在弹出对话框中添加的m文件选中,点击打开,这样就把所需的m文件添加到工程了.m文件必须是函数文件,即必须以function开头,否则会报错,如果是m脚本文件的话,只需要改为无输入输出参数的m函数即可.
选择Build->Com object,点击“Com object”,编译并注册组件,COM组件图如图2所示.
图2 COM组件图
3）打包组件及MCR.
选择Component→Package Component,点击“Com object”,选中“include MCR”,点击create,打包完成后将会在distrib中产生ct.exe的可执行文件与cton_1_0.dll注册文件等.
3.3 C#中调用MATLAB生成的动态链接库
1）C#中添加MATLAB生成的动态链接库.
新建一个C#项目,为体现C#组件界面的强大功能,建立主界面,主界面如图3.
图3 主界面
右击“解决方案资源管理器”的“引用”,选择“添加引用”[7],在弹出的窗体中通过点击“浏览”找到cton_1_0.dll,选择确定之后就可以添加到工程里面.
2）C#中调用MATLAB生成的动态链接库. 以椭圆头模型为例,双击“显示原模型”添加以下程序：
cton.f ff 等于 new f（）；
ff.TYTY（）；
编译运行后,点击“显示原模型”就会出现调用MATLAB组件所得到的图形,图4为显示原模型.
图4 显示原模型
4.软件发布
打开C#,新建项目,在弹出的窗口中点击“其他项目类型”,选择“Visual Studio解决方案”下的“空白解决方案”,将此项目命名为“cton”,确定后右击“解决方案资源管理器”,选择“添加新建项”,在弹出窗体中选择“安装与部署”下的“安装项目”,在“应用程序文件夹”下添加MATLAB创建组件中的cton_mrr文件,cton.ctf,cton.exe,cton_1_0.dll,MCRInstaller.exe,mwutil.dll以及C#中的.dll文件.然后点击“Innovate.exe”,右击添加快捷方式,去掉后缀,然后把这个快捷方式同样剪切放到“用户桌面”文件夹.
右击“cton”工程,依次点击“属性”-“系统必备”,选择“与我应用程序相同的位置下载系统必备组件”单选框,然后点击“确定”.将cton工程编译一下,进入其 debug或者release目录下就会发现Setup.exe 文件,双击该文件即可启动安装
程序.
5.结束语
该软件是基于C#与MATLAB混合编程实现的.充分利用C#与MATLAB各自的优点,较好地实现了CT图像重建的模拟.该软件操作简单,界面友好,实用项强.此外,笔者通过对此项目的研究,熟悉并运用各种重建算法的研究和设计、C#与MATLAB的混合编程、图像的处理,提升了笔者的创新实验能力.
项目基金
海南医学院大学生创新性实验项目[HYCX201211810071]和国家大学生创新创业训练计划项目[201211810071].