医学图像维显示

1、第二章冠状面、矢状面、横断面的提取及交互处理2.1引言通常人们都是通过CT或MRI(核磁共振成像)来获取医学图像的二维切片，但由于原始的CT切片有很多的局限性，例如，二维切片不能够清晰的反应出三维空间的内部组织结构；二维切片的提取所需时间较长；在存取时占用内存空间较大。因此，我们针对原始CT切片的众多局限性分析讨论，并对其改进。在获得医学图像横断面的序列切片后，我们引入了多平面重组法，进而获得了图像的冠状面及矢状面的序列二维切片。目前在医学图像显示领域，图像二维显示已发展的较为成熟，人们也已把目光放在更高维的显示方面。但是，二维显示却也是不可或缺的。对于医生来说，在研究手术方案时，提供一组二维

2、切片是必要的。针对这一方面我们提供了一个便于医生操作的二维显示交互平台。通过此平台，医生可以对二维切片进行变换，并对所显示的任意组织进行测量。2.2二维图像横断面获取在早期，计算机断层成像技术的方法主要是由投影来重建出该对象的二维图像。该投影重建方法在后来的医学图像成像研究中具有很大的影响。在国外，AllanM.Cormack等人早在50年代便开始探讨各种CT的原理。应用傅里叶变换，主要是利用对象截面投影和该截面图像之间的数学关系，重建出该截面图像，即傅里叶投影截面定理。60年代末期，英国EMI公司的实验中心建立了一套CT设备。于1972年春正式发表头部CT影像。随着CT设备的建立，计算机断层

3、成像技术便迅速发展起来，并在医疗诊断方面发挥着重要作用。现在断层图像重建技术发展更加成熟，人们可以重建出精度较高的二维图像。这些图像序列为本文进行二维切片三维显示所需体数据的获取打下很好的基础。2.2.1傅里叶变换在二维图像获取中的应用我们知道，二维图像的一维投影的傅里叶变换9与此二维图像的傅里叶变换的中心剖面相等，而傅里叶变换重建方法也是以此为基础的。通过将投影进行一定角度的旋转和傅里叶变换可以首先构造整个傅里叶变换的平面。若令f(x,y)代表一图像函数，则沿y轴在x轴上的投影示意图如图2.1所示。投影图像表达式为：1/17tg(x)f(x,y)dy（2.1）rr图2.1二维图像在x轴上的投

4、影示意图对投影图像函数作一维傅里叶变换：G(u)g(x)e2uxdxdy（2.2）取f(x,y)的傅里叶变换：F(u,v)f(x,y)e2(uxrv)dxdy（2.3）所以，F(u,0)f(x,y)e2uxdxdy（2.4）由（2.2）（2.4）式可知，F(u,0)G(u)（2.5）y其中F(u,0)为二维图像傅立叶变换中过y0的图像。2.2.2扩展傅里叶变换在二维图像的获取对上面所得到的结论我们进行扩展，讨论沿任意角度进行投影的一般情况。现在若将函数投影方向与y轴成一个角度，设该旋转角度为角。如图2.2所示。2/17图2.2二维图像在s轴上的投影示意图在这里，我们需要定义旋转坐标为：snsx

5、coysinstxsiyco（2.6）（2.7）即沿t轴在S轴上的投影表达式为：g(s,)f(x,y)dt（2.8）SOT坐标系与XOY坐标系的坐标变换关系如下：sscontsinsixscoy（2.9）现在对g(s,)作傅里叶变换为：G(R,)g(s,)e2Rds（2.10）由式（2.9）可知：snsxcoysi（2.11）所以展开后为：G(R,)osf(x,y)e2R(cysin)dxdy（2.12）vRsin若使展开式与投影的二维傅立叶变换相等，则先令指数项做某种代换得到：suRco（2.13）此时相当于点(u,v)位于与u轴成角的直线上，与原点的距离为R。则：3/17G(R,)F(u,

6、v)(2.14)即当变量u、v，和R，满足211式时，二维图像f(x,y)在与y轴成角的方向上的投影的傅里叶变换，等于该图像的二维傅里叶变换。采用不同的方向角插入到三维变换域中，建立一个傅立叶变换空间。在三维空间中的任意平面都可以被重建，则一个二维图像的重建也不失一般性。通过某一图像的不同方向角做投影，经过变换后可得图像中一条过原点的直线上的数据。最终可以获得其二维傅里叶变换，如图2.3所示。图2.3图像的二维傅里叶变换空间从上述论证可知，只要我们有足够多的投影傅立叶变换G(R,)，就能够获得旋转角度前的二维图像的傅立叶变换，通过傅立叶逆变换就可以获得旋转角度前的图像表达式f(x,y)。即：f

7、(x,y)F(u,v)e2(uxvy)dudv（2.15）f(x,y)G(R,)e2R(xcosysin)RdRd（2.16）由以上结论得出：我们可以采取新的方法来获得投影，即按给定的不同方向角经过原点抽取一个截面，再将这个截面作傅里叶逆变换。这一新方法的优点是，将原来在整个三维空域数据场中的重采样变为在对应的频域场中二维截面的抽取，即二维截面的采样，从而降低了采样计算的复杂度。2.2.3对扩展后的傅里叶变换在二维图像获取方式的分析与讨论随着医学图像在各个领域的深入研究，尤其是计算和测量方面所采用手段的日益完善，加大了系统所要处理的数据量。这就给很多设备带来了极大的压力。4/17虽然在研究中我

8、们已经降低了系统处理数据量，但是所花费的计算时间仍不能令人满意。其实上述算法在计算速度方面还具有一定的潜力，但并没有充分发挥出来。数据量的增加导致了存储量的增加这对绘制速度也产生一定的影响。通过分析讨论，我们还发现其存在另一问题。如图2.4所示，运用上述原理我们得到的图像类似于一张X光照片，尤其在医疗诊断方面，医生要求生成的图像具有更高的图像深度信息，且在图像前后遮挡关系方面也没有反映出来。这主要是由于线积分的结果与采样点到屏幕的距离无关；另一方面是和采样点的前后关系无关。目前这些二维图片已经很难满足医生在诊断中的需求，这些图像并不能清晰地反映三维空间的内部结构。图2.4由傅里叶截面定理所产生

9、的人脑横截面图像2.3改进后的二维图像横截面获取随着研究的深入，可视化系统所要处理的数据量日益扩大，针对以上问题，我们需要努力降低存储空间的需求。我们以医学图像为例，目前医学数据的体元规模都能达到512512100个，其为26MB字节，若赋以颜色值（R,G,B），则达到80MB字节左右。这对于很多图像处理工作站来说都是很重的负担。采用上述算法所得到的每一个点都必须由一个复数来表示，这就必须同时存储实部、虚部两部分数值。实际上我们最终得到的二维图像，都是直接定义在实数域中的，然而在中间过程中，却必须保留虚数，这就使得存储空间加大一倍。经过实验我们可以采用完全在实数域中定义的变换方法。5/172.

10、3.1引入Hartley变换对其改进设一维函数f(x)，如果存在与其对应的连续Hartley变换，那么它也是一个一维函数，记为H(u)，其定义为：H(u)f(x)cas(ux)dx（2.17）上式中，cas()函数的定义为：cas()cos()sin()（2.18）由上所诉，我们可以将其推广到多维情况，即：H(u)f(x)cas(ux)dx（2.19）在实际进行数值计算时，使用更多的是Hartley变换的离散形式。三维空域离散数据场可以通过三维离散Hartley变换而得到三维频域离散数据场。其定义为：H(u,u,u)2N31N3d123Nx11x21x31f(x,x,x)cas(2(3d123

11、ux11N1ux2N22ux33)（2.20）N32.3.2改进后的傅里叶变换众所周知，任何一实函数f都可以分解为一个偶函数feven和一个奇函数fodd之和，即：ffevenfodd（2.21）我们可以根据傅里叶变换的线性性质得到：F(f)F(feven)F(fodd)（2.22）（2.24）又因为偶像函数、奇函数的傅里叶变换分别是纯实函数和纯虚函数，所以，假设：F(f)Re(F(f)jIm(F(f)（2.23）（这里，Re(F(f)，Im(F(f)分别是F(f)的实部和虚部），有：Re(F(f)F(f)evenjIm(F(f)F(fodd)那么：6/17H(f)Re(F(f)Im(F(f)

12、F(feven)jF(fodd)F(fevenjfodd)（2.25）由此可见，对实函数而言，改进后的傅里叶变换可以看成是傅里叶变换的线性组合，其本身就是傅里叶变换，故此傅里叶投影截面定理对改进后的傅里叶变换仍是有效的。2.4冠状面、矢状面的断层图像获取通过改进后的傅里叶变换，我们成功的提取了一序列断层图象。但对于本文所要求的，需要分别提取三个不同面的序列二维断层图象，因此，在提取完横截面的基础上，我们还需要对物体的冠状面、矢状面的序列断层图象加以提取。在这里我们提出一种新的方法，该方法不仅能够成功提取出所需的序列断层图象，同时有利于在交互式浏览方面的操作。2.4.1基于序列断层医学图像形成三

13、维数据场在引入多平面重组法前，我们需要通过一组序列断层图象来建立三维数据场。当我们成功提取到一组序列断层图像(CT、MRI、US、PET和SPECT等)后,需要对图像进行一些预处理,例如,格式上的转换（从医学图像格式到便于计算机处理的图像格式）、滤波、位置校正和二维图像剪辑等。然后沿着坐标系中横坐标的方向将这一组序列断层图像一张切片挨着一张切片的重组在一起，形成一长方体。每一张切片中的每一象素是作为体积中的一个体素,若体积数据三维方向采样率不一致,那么我们就需要进行图像插值处理,使得每一体素成为一立方体。经过以上的步骤，我们就对序列断层医学图像形成了三维数据场。在生成三维数据场后，还需要对三维

14、数据进行滤波处理，例如：三维高斯滤波、非线性滤波2等。最后我们要对各项同性的数据进行体素重组，进而实现多平面重组。其步骤如图2.5所示：7/17图2.5医学图像三维数据场形成流程2.4.2多平面重组法当由一组序列断层医学图像形成三维数据场后,体数据是各向同性的,我们在这里可以用一平面方程:axbyczd0（2.26）来切割体数据,从而获得该截面内的体素分布图像。如指定截面方程为:axd0（2.27）byd0（2.28）czd0（2.29）通过方程（2.27）、（2.28）、（2.29）即可分别获得符合医生所要求的冠状面、矢状面、横断面的断层图像。如图2.6所示。图2.6冠状面、矢状面、横断面的

15、二维断层图象2.5冠状面、矢状面、横断面的交互平台冠状面、矢状面、横断面的交互操作平台，主要是方便医生对于一些二维图像数据进行处理及分析，医生可通过一些功能操作来获取所要了解的脑内部组织信息。操作平台主要功能：图像数据的读取操作，可支持加载多种数据格式，如DICOM、BMP、JPEG、IMO、IMP等，图像参数设置；图像数据的平移；图像数据8/17的旋转；图像数据的测量，用户可测得脑内部组织结构的长度、面积等，所测得的结果可通过测量结果区显示。2.5.1冠状面、矢状面、横断面的图像数据读取冠状面、矢状面、横断面的图像数据读取的主要思路是，设置并初始化读取，显示其环境；分析图像数据的具体格式；根

16、据图像数据的格式定义接收图像文件数据的数组并向操作系统申请内存；将图像文件数据读入数组；显示图像数据；图像数据的处理；将内存的数据，按一定的格式和顺序存取。（1）首先我们需要分析一下图像的数据格式。图像数据格式分析是十分重要的，我们需要了解图像文件的种类。由于目前国际都是以DICOM格式为标准，那么我们就针对DICOM格式来进行分析。简单介绍一下DICOM标准。DICOM(医学数字图像和通信)标准，是由美国放射学院(ACR)和美国国家电器制造商协会(NEMA)组成的联合委员会共同制定的，目的是为多种图像设备提供连通性和互操作性，统一各种数字化影像设备的数据格式和数据传输标准，以便于在这些个别系

17、统之间的图像通信和交换，它有效地利用了各项设备的有效资源5。DICOM标准中涵盖了医学数字图像的采集、归档、通信、显示及查询等几乎所有信息交换的协议；以开放互联的架构和面向对象的方法定义了一套包含各种类型的医学诊断图像及其相关的分析、报告等信息的对象集。(2)根据图像数据的格式定义接收图像文件数据的C+数组，并向操作系统申请图像数据的内存空间。对于DICOM、BMP、JPEG、TIFF等标准格式的影像文件，其本身一般不能存储医学影像处理所需的一些信息，需要通过人为设置；对于RAW、IMG等格式文件更是如此。我们还需设置影像的长、宽、数据类型等；而对于DICOM，IMO等格式的图像文件，因其文件

18、内部已完整地保存了所有需要的信息，因此不需要人为设置其它信息。通过以上对医学图像数据存储的分析，我们设置了一个图像参数设置的界面，其中包括图像数据的大小，如图像的长度、宽度等；图像的数据类型，如char型或short型等；图像的像素间距；图像层与层之间的距离；图像数据的排列顺序等。具体的图像参数设置如图2.7所示。9/17图2.7图像参数设置(3)用C+语言的fread()函数将图像的头文件及图像的文件内容读入到C+数组并进行显示，如有需要，可进行一系列的处理。(4)将暂时存放在内存的图像数据，按所需的格式并按一定的顺序使用C+语言的fwrite()函数存回磁盘。该操作平台支持常见的图像文件格

19、式，如：DICOM、BMP、JPEG、TIFF、RAW、IMO、IMG。我们可以通过该操作平台的相应菜单项来读取二维图像数据。可参见图2.8所示。图2.8二维医学图像数据读取界面2.5.2冠状面、矢状面、横断面的图像数据显示对DCIOM文件的计算机显示，需要先提取图像数据的各项参数和图像像素数据值。我们可以先设两个短整型的变量a、b，用以存放读入文件中数据标签的两个16位二进制。然后需要判断DCIOM的编码方式，如果是BGIENDIAN则需要高低位字节交换。然后查得DCIOM的第六部分的数据字典数据标签得出文件的压缩方式。本文中文件采用的是未压缩的方式。以后顺次读出图像数据的各项参数。最后根据

20、(7FE，00010)读出图像的像素数据。计算机通常显示每像素用8位表示的灰度图像。而DICOM医学文件一般每个像素由16位或12位来表示的，因此我们需要进行转换。由于CT图像的灰度级比较大，一般都在40%级以上，而一般的显示器很难一次动态的显示图像的全部10/17信息。为了逼真地显示出医学图像的全部信息。我们通过得到窗位与窗宽，采用图像数据窗口变换的方法。所谓图像数据的窗口变换是指通过限定一个数据观察窗口，将窗口区域的图像数据线性的转换到显示器的最大显示范围。高于或低于窗口上下限的图像数据则分别设置为最亮或最暗的显示值29。这样通过动态地调节窗口的窗宽(图像数据的显示范围)和窗位(显示的图像

21、数据的中心值)，可以依次看到医学图像的全部细节。然后通过下面的公式计算缩放比例，对所有的图像数据进行变换。灰度级小于O的，默认为0；大于255的默认为255。本实验数据的窗宽设为255，窗位设为127.5。如果可以动态调节动态窗口的窗宽与窗位，能够观察到医学图像的全部细节。最后一部分就是图像数据。我们只要完成有关的数据结构和图像数据的转换，就可以显示位图了。但是要注意每一行的字节数必须是4的倍数，不是要补齐。它们存储图像数据的方式也不一样，BMP文件是从下到上，从左到右的。而DICOM图像数据是从上至下、从左到右的。如图2.9所示，为冠状面、矢状面、横断面的图像数据显示。1x（2.30）图2.

22、9冠状面、矢状面、横断面的图像数据显示2.5.3冠状面、矢状面、横断面的图像数据平移对于二维医学图像数据的平移主要原理就是把图像数据中的每一点都按规定好的数据量来进行水平或垂直的移动。先令(x,y)为平移前图像上的一点，设图像数据水平方向移动的距离为t，00 x垂直方向移动的距离为t，则图像数据平移后点(x,y)坐标将变成(x,y)，图y0011像数据平移前后两者的关系表达式为:xxt0y1x0ty采用矩阵形式表示为：11/17y01txy0y01txy10 x10tx11y010011对（2.31）式的矩阵求其逆得：x10tx00y110011即：xxt1xy0 x1ty（2.31）（2.3

23、2）（2.33）通过运用以上公式，可得出图像数据平移后所显示的每一点都能够在原图像数据中找到其对应点。2.5.4冠状面、矢状面、横断面的图像数据旋转图像数据的旋转主要是对图像数据做绕坐标原点为中心的，以角为旋转角度的图像数据旋转。如图2.10所示。先令坐标原点为(x,y)，则图像数据旋转前可用公式（2.34）表示为：00y0rsina()sxrcoa()0（2.34）经旋转角度为旋转后，其坐标变成(x,y)，则图像数据旋转后可用公式11（2.35）表示为：11ysin()cos()0y0y1ysin()xrcos(a)rcos(a)sin()rsin(a)sin()xcos()ysin()00

24、yrsin(a)rsin(a)cos()rcos(a)sin()xsin()ycos()00采用矩阵形式表式为：xcos()sin()0 x11010011对（2.36）式做逆运算得：sxco()sin()0 x0s0co()011001112/17（2.35）（2.36）（2.37）图2.10二维图像数据的旋转2.5.5冠状面、矢状面、横断面的图像数据测量本节主要实现的是面积（病灶面积）、长度（病灶长度）等的测量。二维图像数据的测量主要是关于点的选取，该点在这里我们称为关键点，之后获取关键点在坐标系中的对映位置。找到对映位置的坐标后，我们可以推算出点与点之间的距离。关键点的标记，如图2.11

25、所示。图2.11关键点标记因此，在这里我们需要引入MITK，通过MITK平台来实现这一操作步骤。目前在医学图像研究中VTK和ITK被研究人员广泛应用。其中VTK(visualizationToolkit)是一套进行数据可视化的开发工具包,其完全采用面向对象的设计思想来设计开发。ITK(insightsegmentationandregistrationToolkit)是一个提供医学影像分割与配准的算法平台。但由于一些原因，使用VTK、ITK或使用VTK+ITK,都存在着一定的问题，使它们在应用上具有局限性。因此，由中国科学院自动化研究所开发了一套集成化的医学影像处理与分析的算法平台MITK（M

26、edicalImagingToolKit的缩写）。其主要目的是为医学影像领域提供一套整合了医学图像分割、配准和可视化等功能的,具有一致接口的、可复用的、灵活高效的算法开发工具。对于两点距离测量问题，我们借助于MITK平台所提供的mitkwidgetModel13/17类、mitkManipulator类、mitkobserver类来实现测量功能。mitkwidgetModel类和mitkxmageModel(即二维影像模型)一样，其共同的父类都是mitkModel，是一种模型，也可以绘制并显示出来。该类借鉴了BrookshireD.Conner等在1992年首先提出的3Dwidgets56，即

27、“一种封装了三维几何形状及行为的实体”。其中的“三维几何形状”即指3DWidgets三身的外观，“行为”则指3DWidgets对三维场景中其他物体的控制和对其他物体的信息的显示。mitkwidgetModel类有两个子类:分别是mitkwidgetMode12D和mitkwidgetMode13D。mitkwidgetMode12D类有7个子类，可以绘制二维的矩形、多边形、直线、角度、椭圆等；mitkwidgetMode13D有3个子类，可以在三维场景中绘制直线、平面、角度等，两者经过绘制显示后均可返回指定的信息。因为mitkwidgetMode12D本身也是一种模型，它不能控制mitk工ma

28、geModel模型达到测量的目的，为了实现对mitkwidgetModel的控制和测量，必须添加mitkManipulator类，由它来截获mitklmageView的命令，然后驱动mitkwidgetModel对mitklmageModel进行控制或测量。mitkManiPulator类为交互式地添加所需的测量模型、挑选绘制好的模型进行缩放、移动、删除等操作提供了一个平台，但用户必须自行实现其三个函数(_onMouseDown()、_onMouseMove、_onMouseup()才能实现交互添加测量模型的功能。为了得到mitkwidgetModel或mitklmageModel的一些几何或

29、物理参数，我们还须用到mitkobserver类，该类是一个抽象的类，必须实现其Update()函数，在该函数内把mitkwidgetModel或mitklmageModel的更新数据保存或展示给用户。要注意的是，创建mitkobserver类的对象后必须把其加入mitkwidgetModel对象中进行关联。1、两点距离的测量通过以上对测量的分析，在这里我们借助于mitkLinewidgetMode12D类，医生可以通过交互平台，点击鼠标左键来选取测量的点。我们根据从mitklmageView返回的这两个点的坐标，结合影像像素的间距，可算出两点之间的距离。本实验取两点一是颅脑的上表皮点，令一点为喉部，测得这两点的长度为，如图2.12所示。14/17图2.12两点间的测量2、规则图形面积的测量借助于mitkRectwidgetMode12D类，用户交互式在二维影像中用鼠标左键点击两个点(这两个点为矩形的两个对角点)画出一个矩形，我们根据从mitkImageView返回的这两点的坐标，结合影像像素的间距，可算出矩形四条边的长度，再计算出矩形所含的区域的面积，如图2.13所示。图2.13规则图形面积测量3、不规则图形面积