CN100361626C - 分布式心脏图像多维重建与交互可视化方法 - Google Patents

Abstract

一种采用分布式结构实现旋转扫描二维图像快速重建三维数据场心室功能评价以及实现重建结果的交互显示的方法，其步骤为：用户从客户端或浏览器载入欲处理的心脏二维图像序列；用户通过客户端或浏览器向服务器发送欲处理的心脏二维图像序列和处理请求；服务器接受来自客户端或浏览器的数据，并根据用户请求对心脏二维图像序列作二维或/和三维处理；服务器将处理完的数据发送给相应用户；本地用户在客户端接收和观察处理结果，通过键盘和鼠标控制观察方向绘制区域透明度；远程用户使用浏览器接收和观察处理结果，通过向服务器发送请求来改变观察方向、绘制区域、透明度。

Description

分布式心脏图像多维重建与交互可视化方法

技术领域

本发明属于医学图像成像与处理领域，特别涉及一种心脏图像多维重建与交互可视化方法。

背景技术

心脏超声成像是医学图像中一个发展迅速的领域，三维超声能够显示心脏的立体形态和解剖结构，为心脏疾病诊断以及辅助外科手术提供了一条新途径。因此，超声诊断仪以及相应的应用软件已应用于临床，医院对超声成像设备的需求不断增大。但目前医用超声成像设备的现状是：价格昂贵，单用户响应，功能不够完善。例如，荷兰Philip公司新近推出的Philip 7500实时三维超声成像***，价格高达两百多万人民币，仍然存在三维显示范围比较狭窄的缺点，在某些方面不能满足临床诊断的要求；德国TomTec公司的Compact 3D Cardiology Freehand***售价近20万人民币，但在三维重建的某些环节，还需要医生通过手工绘制出组织的轮廓，不仅耗时，而且重复性差。

对于心脏疾病的诊断，分布式多维重建与交互可视化***非常重要，但市场上尚没有此类***销售，也未见此类***或方法的文献报道。

发明内容

本发明的目的在于克服现有技术的不足，提供一种采用分布式结构实现旋转扫描二维图像快速重建三维数据场、心室功能评价以及实现重建结果的交互显示的方法，按照此种方法建立的***不仅能用于心脏疾病的辅助诊断，而且成本低，便于在各类医院推广。

本发明所述方法在***设计上采用分布式结构，***包括一个服务器、一个或多个客户端或/和浏览器，服务器、客户端和浏览器为软件，分别安装在不同的计算机中，但客户端和浏览器可以安装在同一台计算机中。客户端/服务器为CS模式，采用TCP作为传输协议；浏览器/服务器为BS模式，采用HTTP传输协议，图像序列采用XML描述，服务器收到XML文件后，需要将其转化为图像文件，浏览器也是如此。

本发明所述方法包括以下步骤：

用户从客户端或浏览器载入欲处理的心脏二维图像序列；

用户通过客户端或浏览器向服务器发送欲处理的心脏二维图像序列和处理请求；

服务器接受来自客户端或浏览器的数据，并根据用户请求对心脏二维图像序列作二维或/和三维处理；

服务器将处理完的数据发送给相应用户；

本地用户在客户端接收和观察处理结果，通过键盘和鼠标控制观察方向、绘制区域、透明度；

远程用户使用浏览器接收和观察处理结果，通过向服务器发送请求来改变观察方向、绘制区域、透明度。

本发明所述方法中，服务器对心脏二维图像序列作二维处理的内容为心室轮廓分割，服务器对心脏二维图像序列作三维处理的内容包括心室功能评价、旋转扫描插值和体绘制。

心室轮廓分割即确定心室区域，其步骤如下：

服务器记录用户选择的阈值t和用户选择的位于心室区域的初始像素点位置；

服务器将用户选择的初始像素点邻域的像素点放入一个顺序链表中：

①若链表非空，取出链表中的头元素y，计算y和初始像素点的灰度差的绝对值δ，若链表为空，则执行④；

②若δ≤t，则将该像素点y设为与初始像素点有相同的标识，找出y的邻域像素中既无标识、又不在链表中的元素，将这些元素加入链表，像素点在链表中的位置根据其与初始像素点的δ值按升序确定；

③若δ＞t，则将y标识为边界点，并返回步骤①；

④算法终止。

由上述算法所确定的边界点内的区域即为心室区域。

心室功能评价包括一个心动周期内的心室容积变化曲线、心室收缩末期容积、心室舒张末期容积和射血分数。获取心动周期内某一时相心室容积的步骤如下：

采用旋转扫描插值算法重建心室区域三维数据场；

采用逐层、逐行、逐体素扫描方法计算位于心室腔内的体素总数；

采用超声图像的标尺计算单位长度包含的像素个数，换算可得单位体积所含体素个数；

用心室腔内的体素总数除以单位体积体素个数得到心室容积。

采用上述步骤得到一个心动周期内不同时相所对应的心室容积，即可绘制一个心动周期内的心室容积变化曲线，获得心室收缩末期容积、舒张末期容积并计算出射血分数。

旋转扫描插值是根据旋转扫描得到的二维图像序列的空间分布特征，采用分区线性插值法重建三维数据场，得到心脏的三维结构图像。

体绘制是一种三维数据场可视化技术，将插值后得到的三维数据场作体绘制，用户可以根据需要从不同视角、不同位置观察重建后心脏的三维结构，此即交互可视化。

本发明具有以下有益效果：

1、***设计采用分布式结构，将重建三维数据场放在服务器完成，使得重建结果能同时为多个用户使用，并能响应不同用户的任务请求。

2、客户端/服务器为CS模式，浏览器/服务器为BS模式，前者效率高，后者支持远程操作，升级方便。

3、服务器的处理内容包括心室轮廓分割、心室功能评价、旋转扫描插值和体绘制，用户能够从不同视角、不同位置观察心脏的三维结构，进行心脏功能的定量评价，了解冠状动脉主干及大分支的特征(如管径、内膜、有无斑块等)，观察瓣膜功能和活动情况，以及确定起搏器电极的位置等，为临床诊断提供重要依据。

4、成本低，例如，由一台PC机安装服务器、三台PC机安装客户端组成的一个分布式***，成本只相当于一台德国Tomtec工作站的四分之一，一台Philip 7500的四十分之一。

5、有利于在各类医院、尤其是我国中小医院普及，造福于患者。

附图说明

图1是本发明所述分布式心脏图像多维重建与交互可视化方法的***结构图，图中，2D表示二维，3D表示三维；

图2是客户端的工作流程图；

图3是浏览器的工作流程图；

图4是服务器的工作流程图；

图5是心室轮廓分割流程图；

图6是心室功能评价流程图；

图7是原始图像同轴排列的第k层断面图；

图8是第k层像素插值示意图。

具体实施方式

本发明所述方法需要两台或两台以上具有图像处理能力的PC机，其中一台安装服务器软件，其它的安装客户端或/和浏览器软件。客户端/服务器为CS模式，采用TCP作为传输协议；浏览器/服务器为BS模式，采用HTTP传输协议，图像序列采用XML描述，服务器收到XML文件后，需要将其转化为图像文件，浏览器收到XML文件后，也需要将其转化为图像文件。

按照本发明所述方法，服务器、客户端、浏览器组成的***结构如图1所示。服务器接受客户端或浏览器的请求，负责二维和三维图像处理，并将处理结果返回客户端或浏览器。客户端负责载入二维图像，将用户需要处理的任务请求发送服务器，接受服务器处理结果，并显示、保存处理结果，而且客户端支持返回的三维数据场交互显示。浏览器载入二维图像，将用户请求发送服务器，接受服务器处理结果，显示并保存处理结果。

一、客户端的工作流程

客户端的工作流程如图2所示，具体步骤如下：

(1)根据用户的选择载入心脏二维图像序列。

(2)***初始化：

取得图像扫描类型

取得图像序列文件名

根据扫描角度排序

(3)初始化成功则显示图像，如果用户有图像序列的处理请求，则发送请求到服务器，并等待处理结果；初始化不成功则通知用户可能的错误原因。

(4)如果收到服务器的处理结果，则显示并保存处理结果。例如，如果收到重建后的三维数据场，则可以先保存该数据场，然后对该数据场作交互显示，用户通过键盘和鼠标控制观察方向、绘制区域、透明度等参数。

(5)如果上述过程出错，则通知用户可能的错误原因。

上述三维可视化采用shear-warp体绘制算法(该算法来自Lacroute P，Levoy M，Fastvolume rendering using a shear-warp factorization of the viewing transformation，Proc.SIGGRAPH’94，Orlando，Florida，July，1994：451-458)，基本流程如下：

(1)将重建后的三维数据场作游程长度编码；

(2)进行物质分类；

(3)赋予体素颜色和不透明度；

(4)对数据场作错切变换；

(5)根据用户选择的观察方向，将错切变换后的数据场进行切片合成得到中间图像；

(6)对中间图像作二维变换，得到最终的三维投影图像。

三维可视化也可以采用其他算法，例如光线投影算法。

二、浏览器的工作流程

浏览器的工作流程如图3所示，具体步骤与客户端大致相似。由于浏览器和服务器的通信协议为BS模式，用户请求采用HTTP传输，图像序列采用XML描述，因此重建后的三维数据场不发送回浏览器，而只根据用户请求中包含的可视化参数(如观察方向和位置)返回三维数据场的投影图像。这样，浏览器终端的用户仍然能够根据需要观察心脏的三维结构。三维数据场不通过XML发送的原因在于：(1)三维数据场本身数据量较大；(2)到XML的文件格式转换开销较大；(3)由前两个因素造成的传输滞后，影响使用。由服务器完成可视化，将投影结果作为二维图像传输到浏览器相对快捷。

三、服务器的工作流程

服务器的工作流程如图4所示，具体步骤如下：

(1)用户身份验证

如有用户请求，首先验证用户身份的合法性，若验证通过，则接受用户请求，并记录用户请求的消息ID；若验证未通过，则通知用户原因。

(2)接受用户请求及源数据

接收来自用户(客户端或浏览器)的请求和待处理心脏二维图像数据。

(3)对心脏二维图像序列作二维处理或/和三维处理

根据请求的消息ID调用图像处理模块进行二维处理或/和三维处理。

(4)发送处理结果

处理完毕后将相应结果发送回用户(客户端或浏览器)。

(5)异常通知

如果处理中发生异常，则通知用户原因。

服务器开设多个输入缓冲区(Buffer)，用于缓存来自用户的图像数据。设立一个工作缓冲区，用于处理数据，例如三维重建插值，数据场重建完成后，将其置换出工作缓冲区，并从输入缓冲区调入一组等待处理的数据。服务器采用多线程将结果发送回客户端或浏览器。

服务器对心脏二维图像序列作二维处理的内容为心室轮廓分割，服务器对心脏二维图像序列作三维处理的内容包括心室功能评价、旋转扫描插值和体绘制，下面分别予以说明。

1、心室轮廓分割

心室轮廓分割算法的流程如图5所示，具体原理和步骤如下：

用区域生长实现心室轮廓分割，其基本思想是将具有相似性质的像素集合起来构成区域。设I表示某图像的像素集合，则单个轮廓提取就是从某初始像素点(种子点)出发，在其邻域内找出所有与种子点性质(如灰度、颜色等)相同或相似的像素，合并到区域中，直到没有满足条件的像素为止。

心室轮廓分割算法表示为：

(1)初始化：记录用户选择的阈值t和用户选择的位于心室区域的初始像素点位置，给初始像素点一个标识，并将其邻域内的像素点放入一个顺序链表(Sequentially sortedlist)中。

(2)区域生长：

①若链表非空，取出链表中的头元素y，计算y和种子点的灰度差的绝对值δ；若链表为空，则执行④。

②若δ≤t，则将该像素点y设为与种子点有相同的标识。找出y的邻域像素中既无标识、又不在链表SSL中的元素，将这些元素加入链表，像素点在链表中的位置根据其与种子点的δ值按升序确定。

③若δ＞t，则将y标识为边界点，并返回步骤①。

④生长终止。

δ的定义如下：

δ＝|g(y)-g(s)|   (1)

其中，g(s)为初始像素点的灰度。链表中的初始元素为初始像素点s的邻域像素，且这些像素按其δ值排序。链表中的每个元素包含两个成员：该元素的坐标和δ值。

由上述算法所确定的边界点内的区域即为心室区域。

2、心室功能评价

心室功能评价的流程如图6所示，具体步骤如下：

(1)载入时相i的心脏图像

(2)心室轮廓分割

心室轮廓分割的步骤如上所述。

(3)三维插值重建

三维插值重建算法的原理如下：

假设n幅原始图像的大小均为w×h.则同一时相的原始二维图像在空间同轴排列。相邻两幅图像的夹角为δ＝π/n，且n幅原始图像的像素个数为n×w×h。为计算方便，假设图像宽度w为偶数。

原始图像的分布如图7所示，其特点为：中轴附近，像素点高度密集，越远离中轴，像素点越稀疏。如果用一个方格表示空间直角坐标系下的一个像素，可以看出，原点(中轴)附近，相邻两幅图像的点距小于一个相素大小，远离原点，相邻两幅图像的点距大于一个像素，甚至是数十倍像素；原始图像在柱坐标系中排列。

要显示心脏的三维解剖结构，需要：A、通过插值得到空间各向同性的体素(voxel)；B、同时将柱坐标系变换为空间直角坐标系。

插值后的体数据为w×w×h。如图8，对于第k(0≤k＜h)层断面某一点P(i，j，k)，根据方程：

$\left\{\begin{array}{c}i=r\cos \mathrm{\θ}\\ j=r\sin \mathrm{\θ}\\ r=\sqrt{i^2+j^2}\\ z=k\end{array}\right.---\left(2\right)$

计算该点的柱坐标(r，θ，z)后，若r＞w/2，则该点灰度为0，若r≤w/2，则分两种情形确定该点的灰度。

①中轴附近的像素灰度计算

中轴线附近的像素分布太密集，对此范围采用最临近插值法。为此需要确定该范围(邻域)的大小。若用δw/2表示δ＝π/n对应的弧长，w/2表示半径，1表示单位像素大小，r_c表示相邻两幅图像间点距小于1个像素的邻域半径，则r_c由下式确定：

(δw/2)/(w/2)＝1/r_c    (3)

由上式得出r_c＝1/δ，对r_c取整得[r_c]，当r＜[r_c]时，采用最临近插值法确定像素灰度。

②远离中轴的像素灰度计算

当r≥[r_c]时，空间直角坐标系下的体素点灰度应该采用以下方法计算：

A、根据(r，θ，k)确定与P点相邻的两幅原始图像f_i，f_i+1以及该点在两幅图像中的对应位置(w/2+r，k)，由于w/2+r不见得是像素大小的整数倍，因此可以找到与该位置相邻的四个像素点f_i([w/2+r]，k)，f_i([w/2+r]+1，k)，f_i+1([w/2+r]，k)，f_i+1([w/2+r]+1，k)，如图8。

B、根据下式计算P点的灰度：

f_p＝δ₁t_i+1+(1-δ₁)t_i    (4)

其中，[]表示取整运算，

t_i＝(1-α₁)f_i(w/2+[r]，k)+α₁f_i(w/2+[r]+1，k)，

t_i+1＝(1-α₁)f_i+1(w/2+[r]，k)+α₁f_i+1(w/2+[r]+1，k).δ₁＝θ-iδ，α₁＝r-[r]。

C、确定P点灰度的另一种方法是：除了用P点所在的第k层像素外，还考虑与P点相邻的两层像素，即第k-1和第k+1层像素，每层如步骤B中考虑四个点，这样就可以计算出三个中间灰度f_k-1，p、f_k，p和f_k+1，p，P点的最终灰度由这三个中间灰度的加权来计算，即

f_p＝λf_k-1，p+λf_k+1，p+(1-2λ)f_k，p    (5)

其中，0≤λ＜1/2为常数，在计算中根据需要选定，λ的取值越大，则第k-1和第k+1层像素对P点的灰度贡献越大，反之越小。如果k＝0或k＝h-1，则取λ＝0。

将上述算法命名为旋转扫描分区线性插值法，简称分区线性插值法，其步骤如下：

①初始化：读取原始图像的幅数n，图像宽度w和高度h，根据公式(3)计算出r_c；

②计算w×w×h三维数据场中所有点(体素)的灰度，用P(i，j，k)表示三维数据场中的点：

for(k＝0；k＜h；k++)

   for(j＝-w/2；j＜w/2；j++)

      for(i＝-w/2；i＜w/2；i++){

$\mathrm{if}(r=\sqrt{i^2+j^2}<[r_c\left]\right),$ 采用最邻近插值法计算P点的灰度；

         else，采用公式(4)或(5)计算P点的灰度；

      }

③算法终止。

由上述分区线性插值法即可重建由原始图像序列确定的三维数据场。

(4)心室体素个数计算

将分割后的心室二维图像按时相分组，同一时相下不同角度的图像同轴排列，采用上述分区线性插值法得到心室的三维数据场。

下面以右心室为例说明如何计算心室体素。根据上述插值算法重建的心室数据场有如下特征：①右心室区域位于正方形的一个内切圆之内；②由于右心室是一个腔体结构，因此灰度非常低；③圆内右心室区域外的区域灰度接近或等于255。

对于三维数据场的任一体素(voxel)，上述三个特征足可以区分该体素是否属于右心室腔。设定一个灰度阈值，对三维数据场作逐层、逐行、逐体素扫描，灰度小于此阈值的体素视为在心室区域内，这样可以得出位于右心室腔内的体素个数。

(5)时相i容积计算

计算容积的依据是：图像中某区域的面积可以通过计算该区域的像素个数来完成，三维数据场中某区域的容积则通过计算该区域内包含的体素个数实现。例如：在某次诊断中，通过图像获取时超声仪上的标尺，计算得到10cm包含196.31个像素宽度，可以算得10×10×10cm³的数据场包含约7565426个体素。

一旦得到右心室腔包含的体素个数num，则右心室的容积V由下式计算：

V＝1000num/E_num    (6)

其中，num为右心室包含的体素个数，E_num为1000毫升(ml)的立方体包含的体素个数，V的单位为毫升(ml)。

(6)如果i大于一个心动周期内的时相个数m，则绘制容积变化曲线，否则继续计算下一时相的容积；

(7)完成一个心动周期内不同时相所对应的心室容积计算后，即可绘制该周期内的心室容积变化曲线，得出心室收缩末期容积ESV、舒张末期容积EDV并计算出射血分数RVEF：

RVEF＝(EDV-ESV)/EDV    (7)

至此，心室功能评价完成。

3.旋转扫描插值和体绘制

旋转扫描插值采用上述心室功能评价中的分区线性插值法，体绘制采用上述客户端工作流程中的shear-warp体绘制算法。

Claims (3)

1.一种分布式心脏图像多维重建与交互可视化方法，其特征在于该方法采用分布式结构***，所述***包括一个服务器、一个或多个客户端或/和浏览器，服务器、客户端和浏览器分别安装在不同的计算机中，或客户端和浏览器安装在同一台计算机中，服务器安装在与客户端和浏览器不同的计算机中，服务器开设多个输入缓冲区，用于缓存来自用户的图像数据，设立一个工作缓冲区，用于处理数据，采用多线程将结果发送回客户端或浏览器，心脏图像多维重建与交互可视化包括以下步骤：

用户从客户端或浏览器载入欲处理的心脏二维图像序列；

用户通过客户端或浏览器向服务器发送欲处理的心脏二维图像序列和处理请求；

服务器接受来自客户端或浏览器的数据，并根据用户请求对心脏二维图像序列作二维或/和三维处理，服务器对心脏二维图像序列作二维处理的内容为心室轮廓分割，服务器对心脏二维图像序列作三维处理的内容包括心室功能评价、旋转扫描插值和体绘制；

心室轮廓分割的步骤如下：

服务器记录用户选择的阈值t和用户选择的位于心室区域的初始像素点位置；

服务器将用户选择的初始像素点邻域的像素点放入一个顺序链表中：

①若链表非空，取出链表中的头元素y，计算y和初始像素点的灰度差的绝对值δ，若链表为空，则执行④，

②若δ≤t，则将该像素点y设为与初始像素点有相同的标识，找出y的邻域像素中既无标识、又不在链表中的元素，将这些元素加入链表，像素点在链表中的位置根据其与初始像素点的δ值按升序确定，

③若δ＞t，则将y标识为边界点，并返回步骤①，

④算法终止；

服务器将处理完的数据发送给相应用户；

本地用户在客户端接收和观察处理结果，通过键盘和鼠标控制观察方向、绘制区域、透明度；

远程用户使用浏览器接收和观察处理结果，通过向服务器发送请求来改变观察方向、绘制区域、透明度。

2.根据根据权利要求1所述的分布式心脏图像多维重建与交互可视化方法，其特征在于客户端/服务器为CS模式，采用TCP作为传输协议，浏览器/服务器为BS模式，采用HTTP传输协议。

3.根据权利要求1或2所述的分布式心脏图像多维重建与交互可视化方法，其特征在于心室功能评价包括一个心动周期内的心室容积变化曲线、心室收缩末期容积、心室舒张末期容积和射血分数，

获取心动周期内某一时相心室容积的步骤如下：

采用旋转扫描插值算法重建心室区域三维数据场，

采用逐层、逐行、逐体素扫描方法计算位于心室腔内的体素总数，

采用超声图像的标尺计算单位长度包含的像素个数，进而得到单位体积所含体素个数，

用心室腔内的体素总数除以单位体积体素个数得到心室容积；

采用上述步骤得到一个心动周期内不同时相所对应的心室容积，即可绘制一个心动周期内的心室容积变化曲线，获得心室收缩末期容积、舒张末期容积并计算出射血分数。

Images