CN100577103C - 使用两个图像的图像处理设备和方法 - Google Patents

Abstract

本发明旨在提供能够以高速产生差图像而不会出现失配的图像处理设备和方法。为了实现其，在图像处理设备和方法中，分别设置多个感兴趣区域以输入第一和第二图像，对于每个感兴趣区域计算表示第一和第二图像之间的失配量的偏移向量，对偏移向量执行过滤处理，内插过滤处理的偏移向量，根据内插的偏移向量配准第一和第二图像，以及在各配准的图像上的相应像素之间执行相减处理以获取差图像。

Description

使用两个图像的图像处理设备和方法

技术领域

本发明涉及图像处理设备和方法，其适用于在利用医学图像的CAD(计算机辅助诊断)技术等中作为处理如差处理的目标的两个图像的配准。

背景技术

近年来，在CAD领域，数字图像的使用得到发展。因此，通过使医学图像数字化，出现了在利用银盐摄影的传统诊断中难以进行的诊断形式的可能性。

更具体地说，在传统诊断中，在观察患者的状况期间在不同时间点拍摄的多个X射线照片被比较用于诊断的情况下，其上已分别显影X射线照片的胶片通常被悬挂在灯箱(schaukasten)上，并且悬挂的胶片被实际地比较及读取。

同时，在数字图像被用于诊断的情况下，关于一个患者的在不同时间点拍摄的两个数字图像要进行配准(registration)，以使一个数字图像的正常解剖结构与另一个数字图像的解剖结构相符，然后，对两个数字图像执行相减处理，由此产生并输出差图像。随后，将输出的差图像与该对两个原始数字图像相比较，由此可以更精确地掌握两个原始图像之间的变化。

在例如对应于USP 5359513并在下文中被称作文件1的日本专利申请公开No.H07-037074中公开了这样一种差图像产生方法。即，根据文件1中公开的产生方法，在不同时间点分别拍摄的两个胸部X射线图像要进行配准，并可产生差图像。这里，应当注意，这样的相减处理被称作时间相减处理。

随后，将参照图25描述如文件1中所公开的实现时间相减处理的设备的示意构成。

在图25中，首先，由图像输入单元1输入的一对医学数字图像通过预处理单元11进行密度校正处理，然后，输入到ROI(感兴趣区域)匹配单元12。在ROI匹配单元12中，通过计算互相关系数对于多个设置的ROI(感兴趣区域)执行匹配处理，并且对于每个ROI计算表示该对医学数字图像(两个图像)的位移量的偏移向量。

然后，在多项式内插单元13中，通过二元n级多项式对计算的偏移向量进行近似内插。随后，在配准单元5中，对两个图像的任意一个应用非线性失真。此外，在相减操作单元6中，在相应位置的像素之间执行相减操作，由此产生差信号。之后，在后处理单元7中，对差信号执行包括灰度处理等的后处理，并且处理的信号被输出到输出单元8。

此外，用于胸部X射线图像的时间相减技术是用于处理作为在不同时间点拍摄的人体的一部分的公用对象的第一和第二图像的技术。更具体地说，时间相减技术校正肺野的变形，该变形是由于各种因素比如对象的前后及左右移动、对象的呼吸、X射线管照射角度的改变等引起的，时间相减技术执行相减处理，然后提取包括变化的部分作为两个图像的差图像。通过应用上述时间相减技术，可以从第一和第二图像中提取仅对应于病变位置的变化的图像分量，而不会涉及一般的正常组织比如骨头、血管等的图像分量。因此，特别是在时间相减CAD技术中，可以临床期待早期病变检测、早期检测隐藏在正常组织比如肋骨及血管等之后的病变部位、防止病变的失察以及快速判读射线照片。

在任何情况下，时间相减技术的主要要素是用于校正在第一和第二图像之间出现的变形的配准技术。顺便地说，对应于美国专利公开No.2001048757并在下文中被称作文件2的日本专利申请公开No.2002-032735公开了常规的时间相减技术。更具体地说，在这样的常规时间相减技术中，执行如图26所示的处理。即，首先读取第一图像(原始图像或过去的图像)和第二图像(当前图像)，然后在第一图像的肺野区内均匀设置模板ROI。随后，在第二图像中，在对应于第一图像的模板ROI的位置设置搜索ROI。此时，在第二图像的搜索ROI中，搜索对应于第一图像的模板ROI的中心的位置，并且，从第一图像的模板ROI的中心到第二图像的搜索ROI中的相关位置的转变作为偏移向量被记录。

在实际记录偏移向量的情况下，模板ROI的中心坐标和从模板ROI的中心到搜索ROI中相关位置的转变被记录。典型地，在实现ROI的一致性(匹配)的情况下，匹配度被用作偏移向量的权重。然后，通常，在通过利用ROI的互相关实现匹配的情况下，相关系数本身如原样被用作偏移向量的权重。此外，在通过利用SSDA(序列类似性检测算法)实现匹配的情况下，通过计算并使残差的逆数归一化而获得的结果被用作偏移向量的权重。之后，通过利用获得的权重，对偏移向量执行利用多项式的内插，并且第二图像被卷曲(warp)到第一图像以获取差图像。

然而，在上述文件1中，当通过执行对于多个ROI的每个的匹配(ROI匹配)获得的偏移向量通过多项式近似被内插时，多项式的系数通过最小二乘法等被确定。为此，存在执行这样的处理需用较长时间的技术问题。

此外，在ROI匹配中，如果在对象中存在多个类似图案，则存在匹配精度的限制。因此，根据情况，在偏移向量中不可避免地会包括严重误差。在这样的情况下，如果偏移向量通过使用最小二乘法被内插，则包括的误差影响其它的偏移向量，由此整体上发生位移或失配。为此，存在在差图像中噪声分量增加的技术问题。

顺便提及，如下面所述，胸部简单X射线摄影包括信息量有别于其它的各种区域。

即，对于锁骨和身体边缘部分来说，尽管灰度信息是缺乏或简单的，但边缘信息是充分的。在肺野边缘部分，灰度信息和边缘信息均是充分的。此外，在心脏和横膈膜中，灰度信息和边缘信息均是简单或缺乏的。因此，例如在文件2中，当通过利用灰度信息或边缘信息缺乏的ROI执行匹配时，不一定能认为ROI的匹配被精确地执行，即使其偏移向量的权重较高，并且在这种情况下不可能判断偏移向量的信息是否被正确。为此，存在不利地影响随后的处理的技术问题。

此外，对于包括灰度信息、边缘信息等的纹理(texture)充分的ROI的偏移向量，希望执行精度高于纹理低的ROI的偏移向量的内插的内插，由此需要进行更精确的偏移向量内插。

发明内容

本发明被实现以解决上面的技术问题，其目的是提供能够以高速产生差图像同时相对减小失配的图像处理设备和方法，以及计算机使用的用来实现上述图像处理方法的程序。

用于实现上述目的的根据本发明的图像处理设备的特性在于：一种基于偏移向量将第二图像与第一图像相匹配的图像处理设备，包括：

匹配度获取单元，适于获取第一和第二图像之间相互对应的相应点的匹配度；

纹理估算单元，适于估算相应点外围的纹理；以及

偏移向量加权单元，适于基于匹配度和纹理的估算加权偏移向量。

用于实现上述目的的根据本发明的图像处理方法的特性在于：一种基于偏移向量将第二图像与第一图像相匹配的图像处理方法，包括：

匹配度获取步骤，用于获取第一和第二图像之间相互对应的相应点的匹配度；

纹理估算步骤，用于估算相应点外围的纹理；以及

偏移向量加权步骤，用于基于匹配度和纹理的估算加权偏移向量。

用于实现上述目的的根据本发明的图像处理设备的特性在于：一种输出第一图像和第二图像之间的差图像的图像处理设备，包括：

输入单元，适于输入第一图像和第二图像；

偏移向量计算单元，适于将多个感兴趣区域分别设置到第一图像和第二图像，并且对于每个感兴趣区域计算指示第一图像和第二图像之间的失配量的偏移向量；

过滤单元，适于对偏移向量执行过滤处理；

内插单元，适于对通过过滤单元进行过滤处理的偏移向量执行内插；

配准单元，适于根据通过内插单元内插的偏移向量配准第一图像和第二图像；

相减操作单元，适于在各配准的图像上的相应像素之间执行相减操作；以及

输出单元，适于输出通过相减操作单元获取的差图像。

用于实现上述目的的根据本发明的图像处理方法的特性在于：一种输出第一图像和第二图像之间的差图像的图像处理方法，包括：

输入步骤，输入第一图像和第二图像；

偏移向量计算步骤，将多个感兴趣区域分别设置到第一图像和第二图像，并且对于每个感兴趣区域计算指示第一图像和第二图像之间的失配量的偏移向量；

过滤步骤，对偏移向量执行过滤处理；

内插步骤，对在过滤步骤中进行过滤处理的偏移向量执行内插；

配准步骤，根据在内插步骤中内插的偏移向量配准第一图像和第二图像；

相减操作步骤，在各配准的图像上的相应像素之间执行相减操作；以及

输出步骤，输出在相减操作步骤中获取的差图像。

附图说明

图1是示出根据本发明的第一实施例的医学图像处理设备的功能构成的功能方框图；

图2A、2B、2C和2D是示出分别设置在肺野边缘、肺野中心、纵隔和横隔膜中的ROI的图像的示图；

图3A是示出图2A的图像的直方图的示图；

图3B是示出图2B的图像的直方图的示图；

图3C是示出图2C的图像的直方图的示图；

图3D是示出图2D的图像的直方图的示图；

图4是示出由加权处理单元90获取的偏移向量的权重的示例的示图；

图5是示出根据本发明的第一实施例的医学图像处理设备的操作的流程图；

图6是示出通过普通的胸部简单射线摄影检测的肺野边缘的典型示图；

图7是示出根据本发明的第二实施例的医学图像处理设备的操作的流程图；

图8是示出根据本发明的第三实施例的医学图像处理设备的操作的流程图；

图9A、9B、9C和9D是分别示出第三实施例的边缘检测结果的示图；

图10是示出根据本发明的第五实施例的图像处理设备的构成的方框图；

图11是示出能够实现根据本发明的第五实施例的图像处理设备的计算机系统的构成的示例的方框图；

图12是示出要由根据本发明的第五实施例的图像处理设备执行的整个处理的流程图；

图13是示出根据本发明的第五实施例的灰度转换特性曲线的示例的示图；

图14是示出根据本发明的第五实施例的偏移向量计算单元的详细构成的方框图；

图15是示出根据本发明的第五实施例的偏移向量计算单元的详细操作的流程图；

图16A和16B是用于描述根据本发明的第五实施例的ROI的设置的示图；

图17A、17B、17C、17D和17E是用于描述根据本发明的第五实施例的胸腔检测方法的示图；

图18是示出根据本发明的第五实施例的偏移向量的示例的示图；

图19是示出根据本发明的第五实施例的过滤单元的详细操作的流程图；

图20是示出根据本发明的第五实施例的偏移向量分量的示例的示图；

图21是示出根据本发明的第五实施例的内插的偏移向量的示例的示图；

图22是示出要由根据本发明的第六实施例的图像处理设备执行的整个处理的流程图；

图23是用于描述根据本发明的第七实施例的区域分割的示例的示图；

图24是用于描述根据本发明的第八实施例的样条内插的示图；

图25是示出常规图像处理设备的示图；以及

图26是示出常规图像处理方法的示图。

具体实施方式

下面将参照附图具体描述本发明的实施例。

<第一实施例>

首先，下文中将描述本发明的第一实施例。图1是示出根据本发明的第一实施例的医学图像处理设备的功能构成的功能方框图。顺便提及，应当注意，通过实现图1中所示功能的专用设备或通过使通用计算机执行后面描述的处理的控制程序，可实现根据本实施例的医学图像处理设备。然而，应当注意，通过硬件、软件或硬件和软件的组合均可实现图1中所示的每个功能块。

如图1中所示，根据本实施例的医学图像处理设备装备有图像输入单元10、模板ROI(感兴趣区域)设置单元20、搜索ROI匹配单元30、ROI纹理计算单元40、匹配度计算单元50、偏移向量计算单元60、偏移向量加权计算单元70和偏移向量内插单元80。

图像输入单元10从图像存储单元(未示出)读取两个数字图像(即第一和第二图像)。例如，图像输入单元10读取在不同时间点分别拍摄的同一对象的相同区域的两个图像。然后，模板ROI设置单元20在第一图像的肺野区域中均匀设置模板ROI，并且搜索ROI匹配单元30在对应于第一图像的模板ROI的第二图像的位置设置搜索ROI。这里，应当注意，搜索ROI被设置为大于模板ROI。

ROI纹理计算单元40计算由模板ROI设置单元20设置的模板ROI的纹理。这里，应当注意，包括在模板ROI中的纹理被不同地估算，并且模板ROI中像素值的变化被用作一个指标(index)。为此，在本实施例中，ROI纹理计算单元40形成模板ROI的直方图，并且将非零计数的像素值数C_i设置为纹理。这里，符号i表示模板ROI数。在任何情况下，通过执行这样的处理，设置在纵隔、心脏、横膈膜等区域中的模板ROI的纹理减少。相反，设置在肺野、肺野轮廓等中的模板ROI的纹理增加。

图2A、2B、2C和2D是示出分别设置在肺野边缘、肺野中心、纵隔和横隔膜中的模板ROI的图像(均具有12位的数据量)的示图，图3A是示出图2A的图像的直方图的示图，图3B是示出图2B的图像的直方图的示图，图3C是示出图2C的图像的直方图的示图，图3D是示出图2D的图像的直方图的示图。这里，在图3A-3D的每个图中，数据量被设置为8位。如图3A中所示，在设置在肺野边缘的模板ROI的直方图中，非零计数的像素值在大致从180延伸到230的约50宽度的范围内分布。如图3B中所示，在设置在肺野中心的模板ROI的直方图中，非零计数的像素值在大致从150延伸到180的约30宽度的范围内分布。一方面，如图3C中所示，在设置在纵隔的模板ROI的直方图中，非零计数的像素值在大致从240延伸到250的约10宽度的范围内分布。此外，如图3D中所示，在设置在横膈膜的模板ROI的直方图中，非零计数的像素值仅在大致从240延伸到256的约16宽度的范围内分布。因此，可以理解，设置在肺野边缘或肺野中心中的模板ROI的可靠性高于设置在纵隔或横膈膜中的模板ROI的可靠性。

当将模板ROI转换到第二图像中的搜索区时匹配度计算单元50计算互相关系数。这里，应当注意，互相关系数为最大的位置对应于模板ROI的中心。因此，在本实施例中，在模板ROI中心的互相关系数R_i被用作匹配度。

偏移向量计算单元60计算在第二图像中互相关系数为最大的位置和第一图像中模板ROI的中心之间的位移(或失配)，并且获取水平和垂直方向的偏移量作为偏移向量。

根据通过ROI纹理计算单元40获取的模板ROI的非零计数的像素值数C_i和通过匹配度计算单元50获取的最大互相关系数R_i，偏移向量加权计算单元70对于所有ROI的每个计算最大互相关系数R_i和模板ROI的纹理T_i的积，并且将通过方程式(1)获取的归一化结果设置为每个ROI的权重w_i。

$W_i=\frac{T_i\&CenterDot;R_i}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{T}_i\&CenterDot;R_i}...\left(1\right)$

这里，符号N为ROI的数目，并且在本实施例中满足T_i＝C_i。

在本实施例中，ROI纹理计算单元40、匹配度计算单元50和偏移向量加权计算单元70共同构成加权处理单元90(图1)。以这种连接，图4示出了通过加权处理单元90获取的偏移向量的权重的示例。这里，应当注意，图4中所示的示例从图26中所示的普通X射线摄影中被获得。如图4中所示，设置在胸腔域边缘附近的ROI的偏移向量具有较大的权重，设置在肺野中心、纵隔、心脏、腹部等附近的ROI的偏移向量具有较小的权重。

偏移向量内插单元80执行通过偏移向量加权计算单元70获取的偏移向量的权重和通过偏移向量计算单元60计算的偏移向量的多项式内插。因此，第一图像与第二图像的变形通过一个二元多项式被表示，第二图像通过随后的处理被卷绕(warp)到第一图像，然后，第一图像和第二图像之间的相减处理被执行以获取差图像。

随后，下文中将说明按上述构成的医学图像处理设备的操作。图5是示出根据本发明的第一实施例的医学图像处理设备的操作的流程图。

在本实施例中，首先，图像输入单元10将来自未示出的图像存储单元的第一图像读入存储器中(步骤S101)。然后，在第一图像中检测肺野的边缘，并且设置肺野区(步骤S102)。图6是示出通过普通的胸部简单射线摄影检测的肺野边缘的典型示图。在图6中，边缘的最上端、最下端、最左端和最右端的中心被看作肺野的中心，从最上端和最下端延伸距离Δh以及从最左端和最右端延伸距离Δw的矩形区域被看作肺野区10’。然而，如果按上述设置的矩形区超过图像的边界，则通过修整超出矩形区的部分而获取的结果被看作肺野区10’。

在肺野区被设置之后，模板ROI设置单元20通过在肺野区均匀地设置模板ROI中心设置模板ROI(步骤S103)。

随后，根据统计，ROI纹理计算单元40形成关于每一设置的ROI的图像直方图(步骤S104)。然后，ROI纹理计算单元40获得关于每一直方图的非零计数的像素值数C_i(步骤S105)。

之后，图像输入单元10将来自未示出的图像存储单元的第二图像读入存储器中(步骤S106)。随后，与步骤S102一样，在第二图像中检测肺野的边缘，并且获取肺野区的中心，由此设置肺野区(步骤S107)。

下一步，基于在步骤S102和S107的每个中检测的肺野边缘的最上端、最下端、最左端和最右端，在第一和第二图像之间的整个肺野的偏移量和垂直/水平方向的放大/缩小率被计算(步骤S108)。

随后，搜索ROI匹配单元30通过下面的方程式(2)将对应于步骤S103中设置的模板ROI的中心的位置设置在第二图像中。然后，搜索ROI匹配单元30根据设置的位置设置搜索ROI(步骤S109)。

x′＝r_h·x+ΔH

y′＝r_v·y+ΔV    ...(2)

这里，应当注意，(x，y)表示在第一图像中的肺野中心被设置为原点的情况下模板ROI的中心的坐标，(x’，y’)表示在第二图像中的肺野中心被设置为原点的情况下搜索ROI的中心的坐标，符号ΔH表示沿水平方向肺野中心的整个偏移量，符号ΔV表示沿垂直方向肺野中心的整个偏移量，符号r_h表示水平方向的放大/缩小率，符号r_v表示垂直方向的放大/缩小率。

下一步，判断是否对所有ROI执行了稍后描述的步骤S111、S112、S113、S114和S115的处理(步骤S110)。如果判断对所有ROI执行了这些处理，则流程进到步骤S116。同时，如果判断存在未对其执行这些处理的ROI，则流程进到步骤S111。

在步骤S111中，判断是否对搜索范围内的所有位置执行了随后步骤S112中的处理。如果判断对所有位置执行了有关的处理，则流程进到步骤S113。同时，如果判断存在未对其执行有关处理的位置，则流程进到步骤S112。

在步骤S112中，匹配度计算单元50计算在模板ROI和搜索ROI中相应区域之间的互相关系数R_i。

在步骤S113中，匹配度计算单元50对于所有搜索位置找出最大互相关系数的位置。然后，偏移向量计算单元60根据在步骤S113中找出的位置计算偏移向量(步骤S114)。随后，对于所有搜索位置，通过利用在步骤S105中获取的模板ROI的纹理和在步骤S113中获取的最大互相关系数，根据方程式(1)分别计算偏移向量的权重(步骤S115)。

然后，偏移向量内插单元80通过利用所有偏移向量执行内插(步骤S116)。随后，第一图像和第二图像之间的相减处理被执行。

根据如上所述的第一实施例，要被给定到纹理充分的ROI的偏移向量的权重高于要被给定到纹理贫乏的ROI的偏移向量的权重，因此，与相关的背景技术相比，可以更精确地执行内插。为此，例如，可以增加在不同时间点分别拍摄的两个图像之间的配准的精度。

这里，应当注意，第一实施例中的处理顺序不限于图5的流程图中所示的处理顺序。即，本实施例的功能也可通过另一个过程来实现。此外，应当注意，方程式(1)可被应用在模板ROI的面积均一致的情况下。即，在ROI的面积各不相同的情况下，可通过使用方程式(3)来执行归一化。

$W_i=\frac{T_i\&CenterDot;R_i/S_i}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(T_i\&CenterDot;R_i/S_i)}...\left(3\right)$

这里，符号S_i为第i个ROI的面积，并且与方程式(1)一样，在本实施例中满足T_i＝C_i。

此外，除ROI的直方图的非零计数的像素值数C_i之外，像素值的离差、直方图的峰度等也可被加权并被添加到纹理中。

<第二实施例>

接下来，下文中将描述本发明的第二实施例。在第二实施例中，应当注意，功能块基本上与第一实施例中的功能块相同，仅ROI纹理计算单元40的功能与第一实施例中的不同。图7是示出根据本发明的第二实施例的医学图像处理设备的操作的流程图。

在本实施例中，在与第一实施例一样的对于第一图像设置ROI(步骤S103)之后，通过方程式(4)获取FFT(快速傅立叶变换)系数(步骤S201)。

$F(p,q)=\underset{m=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}I(m,n)e^{-j(2\mathrm{\&pi;}/M)\mathrm{pm}}{e}^{-j(2\mathrm{\&pi;}/N)\mathrm{qn}}...\left(4\right)$

其中p＝0，1，...，M-1，并且q＝0，1，...，N-1

这里，符号M表示模板ROI的水平像素数，符号N表示模板ROI的垂直像素数。

下一步，对于FFT系数，获取除p＝0，1，...，M/2和q＝0，1，...，N/2的频率分量之外的高频分量的绝对值的和SH_i(步骤S202)。这里，符号i表示ROI数。

${\mathrm{SH}}_i=\underset{p=0}{\overset{M/2}{\mathrm{\&Sigma;}}}\underset{q=N/2+1}{\overset{N}{\mathrm{\&Sigma;}}}\left|F(p,q)\right|+\underset{P=M/2+1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{q=0}{\overset{N}{\mathrm{\&Sigma;}}}\left|F(p,q)\right|...\left(5\right)$

随后，如方程式(6)所示，通过将左上区的频率分量的绝对值加到和SH_i中，获取整个频率分量的绝对值的和SA_i(步骤S203)。

${\mathrm{SA}}_i=\underset{p=0}{\overset{M}{\mathrm{\&Sigma;}}}\underset{q=0}{\overset{N}{\mathrm{\&Sigma;}}}\left|F(p,q)\right|...\left(6\right)$

之后，通过方程式(7)获取比率值Ra_i，并且该获取的值被设置为ROI的纹理(步骤S204)。

${\mathrm{Ra}}_i=\frac{{\mathrm{SH}}_i}{{\mathrm{SA}}_i}\&times;100\%...\left(7\right)$

在任何情况下，在常规的胸部简单射线摄影中，肺野边缘的ROI、肺野中心的ROI、纵隔的ROI和横膈膜的ROI的各高频分量比率Ra如下所示：

肺野边缘的ROI：20.73％

肺野中心的ROI：23.81％

纵隔的ROI：    6.22％

横膈膜的ROI：  3.33％

如上所述，设置在肺野边缘或肺野中心的ROI包括较高的纹理。

之后，与第一实施例一样，步骤S106和下面步骤中的处理被执行。然而，当通过利用方程式(1)执行归一化时，偏移向量的权重在T_i＝Ra_i的条件下被计算。

如上所述，在第二实施例中，ROI的FFT(快速傅立叶变换)系数被获取，并且高频分量与整个分量的比率被设置为ROI的纹理。因此，与第一实施例一样，与相关的背景技术相比，可以更精确地执行内插，由此可以增加在不同时间点分别拍摄的两个图像之间的配准的精度。

顺便提及，在本实施例中，ROI的图像可通过使用DCT(离散余弦变换)、小波变换等而非FFT来变换，以将高频分量与整个频率分量的比率设置为纹理。此外，在上述说明中，小于预定值的频率p和q的每一个被设置为低频。然而，如果(p+q)或(p²+q²)^1/2小于预定值，则频率p和q的每一个可被设置为低频，以便计算高频分量比率Ra。此外，当和SH_i或SA_i被计算时，每个频率分量的平方和可代替绝对值的和被使用。而且，可以获取ROI的低频分量与整个频率分量的比率Ra_L，并且随后将“1-Ra_L”设置为纹理。

<第三实施例>

随后，下文中将说明本发明的第三实施例。在第三实施例中，应当注意，功能块基本上与第一实施例中的功能块相同，仅ROI纹理计算单元40的功能与第一和第二实施例中的不同。图8是示出根据本发明的第三实施例的医学图像处理设备的操作的流程图。

在本实施例中，在与第一实施例一样的对于第一图像设置ROI(步骤S103)之后，如方程式(8)所示的水平苏贝尔算子(Sobel operator)被乘到ROI，从而计算位置(i，j)的图像的水平边缘强度(intensity)bx(i，j)(步骤S301)。然后，如方程式(8)所示的垂直苏贝尔算子被乘到ROI，从而计算位置(i，j)的图像的垂直边缘强度by(i，j)(步骤S302)。

                 -1 0  1

水平苏贝尔算子   -2 0  2

                 -1 0  1

                 -1 -2 -1

垂直苏贝尔算子   0  0  0

                 1  2  1    ...(8)

随后，根据位置(i，j)的图像的水平边缘强度bx(i，j)和垂直边缘强度by(i，j)计算位置(i，j)的图像的梯度的强度g(i，j)(步骤S303)。

g²(i，j)＝bx²(i，j)+by²(i，j)  ...(9)

之后，判断是否对所有像素位置执行了稍后描述的步骤S305、S306、S307和S308的处理(步骤S304)。如果判断对所有像素位置执行了这些处理，则流程进到步骤S309。同时，如果判断存在未对其执行这些处理的像素位置，则流程进到步骤S305。

在步骤S305中，判断梯度的强度g(i，j)是否大于预定阈值。如果判断强度g(i，j)大于预定阈值，则流程进到步骤S306。同时，如果判断强度g(i，j)等于或小于预定阈值，则流程返回到步骤S304。

在步骤S306中，如果bx(i，j)＞by(i，j)，则判断该边缘为水平边缘。此外，还判断水平边缘强度bx(i，j)是否为最大。如果满足bx(i，j)＞by(i，j)并且水平边缘强度bx(i，j)为最大，则流程进到步骤S308。同时，如果不满足bx(i，j)＞by(i，j)或者水平边缘强度bx(i，j)不为最大，则流程进到步骤S307。

在步骤S307中，如果by(i，j)＞bx(i，j)，则判断该边缘为垂直边缘。此外，还判断垂直边缘强度by(i，j)是否为最大。如果满足by(i，j)＞bx(i，j)并且垂直边缘强度by(i，j)为最大，则流程进到步骤S308。同时，如果不满足by(i，j)＞bx(i，j)或者垂直边缘强度by(i，j)不为最大，则流程返回到步骤S304。

然后，在步骤S308中，位置(i，j)被检测为边缘(步骤S304)，并且流程返回到步骤S304。

此外，在步骤S309中，所检测边缘的数量与ROI的面积的比率(边缘比)P_i被计算，并且比率P_i被设置为ROI的纹理。这里，符号i表示ROI数。

例如，在一般的胸部简单射线摄影中，肺野边缘的ROI、肺野中心的ROI、纵隔的ROI和横膈膜的ROI的检测结果分别在图9A、9B、9C和9D中被示出。顺便提及，这些ROI的各边缘比P_i如下所示：

肺野边缘的ROI：6.66％

肺野中心的ROI：6.30％

纵隔的ROI：    0.00％

横膈膜的ROI：  0.00％

如上所述，对于设置在纵隔或横膈膜中的ROI，不能执行满意的匹配。

之后，与第一实施例一样，步骤S106和下面步骤中的处理被执行。然而，当利用方程式(1)执行归一化时，偏移向量的权重在T_i＝P_i的条件下被计算。在这种情况下，由于设置在纵隔或横膈膜中的模板ROI的偏移向量的权重较低，因此计算这样的权重的处理被省略，从而在通过内插方程式计算偏移向量中不使用相关的权重。

如上所述，在第三实施例中，在每个ROI中边缘检测被执行，并且边缘数量与ROI面积的比率被设置为ROI的纹理。因此，与第一实施例一样，与相关的背景技术相比，可以更精确地执行内插，从而可以增加在不同时间点分别拍摄的两个图像之间的配准的精度。

顺便提及，在执行边缘检测以计算本实施例的ROI的纹理时，Prewitt法、Roberts法、Canny法等均可代替苏贝尔(Sobel)算子被使用。

此外，在第一至第三实施例中，在通过偏移向量加权计算单元70计算偏移向量的权重的情况下，ROI的纹理可被归一化，而非通过使用方程式(3)进行归一化，以使归一化纹理和匹配度的加权和被设置为偏移向量的权重。

<第四实施例>

接下来，在下文中将说明本发明的第四实施例。在第四实施例中，应当注意，功能块基本上与第一实施例中的功能块相同，但仅偏移向量加权计算单元70与第一实施例中的不同。

在本实施例中，偏移向量加权计算单元70检测胸腔(或胸廓)，对设置在胸腔附近的ROI给定较大的权重，对其它ROI给定较小的权重。然后，偏移向量加权计算单元70将给定权重和匹配度的归一化结果设置给偏移向量。

在这样的处理中，例如，其中存在胸腔边缘的模板ROI可被判断为设置在胸腔附近的ROI。此外，通过计算模板ROI的中心和胸腔边缘之间的水平及垂直距离，可以基于计算距离的较小的一个判断在胸腔附近是否存在相关的边缘。另外，通过计算模板ROI的中心和胸腔的最近边缘之间的距离，可以基于计算的距离判断在胸腔附近是否存在相关的边缘。

根据上述第四实施例，可以具有与第一至第三实施例中的效果一样的效果。

<第五实施例>

图10是示出根据本发明的第五实施例的图像处理设备的构成的方框图。

顺便提及，应当注意，图像处理设备的各构成部件可通过专用硬件或通过在通用计算机上操作程序来实现。在后面的情况下，当CPU执行相关程序的模块时，图10中所示的各构成部件可被实现。

在下文中，将说明能够实现根据第五实施例的图像处理设备的计算机系统的构成的一个示例。

图11是示出能够实现根据本发明的第五实施例的图像处理设备的计算机系统的构成的一个示例的方框图。

用作根据第五实施例的图像处理设备的计算机2000可通过网络1001与用于产生医学X射线图像的图像产生设备1000以及文件服务器1002相连。然而，当然也可以单独构成计算机2000。

计算机2000包括加速器2001、硬盘2002、CPU 2003、RAM 2004和ROM 2005，并且与诸如磁光盘2007、鼠标2008、键盘2009、打印机2010和显示设备2011之类的各种外围设备相连。这里应当注意，这样的构成部件通过总线2006彼此相互连接。

CPU 2003控制通过总线2006相连的各构成部件，以实现第五实施例中的图像处理设备。加速器2001实现各种图像处理功能，并且还通过与CPU 2003合作实现各种处理。硬盘2002存储涉及实现第五实施例的各种处理的控制程序和诸如要被处理的图像数据等之类的数据。

RAM 2004用作各种数据的工作区和出栈(pull-off)区，ROM2005存储各种数据，比如控制程序、各种参数等。通常被称作MO的磁光盘2007存储控制程序、比如要被处理的图像数据的各种数据等。鼠标2008或键盘2009用作输入处理执行指令和各种数据的输入设备。这里，除鼠标2008和键盘2009之外，也可使用另外的点击设备比如笔等。

打印机2010打印各种数据，比如要被处理的图像数据等。这里，应当注意，可以使用各种方法作为打印机2010的打印方法，比如喷墨打印方法、激光束打印方法、热转印打印方法等。显示设备2011显示用来执行各种处理的操作屏幕，还显示各种处理结果。这里，应当注意，作为显示设备2011，可以使用CRT、LCD(液晶显示器)等。

此外，计算机2000可通过由接口(未示出)连接的网络1001将图像数据传送到外部设置的图像产生设备1000和文件服务器1002，或从外部设置的图像产生设备1000和文件服务器1002接收图像数据。

在这样的构成中，用于实现根据本发明的图像处理设备的功能的程序被存储在例如硬盘2002或通过网络1001连接的文件服务器1002中。然后，通过使用输入设备比如鼠标2008、键盘2009等响应用户的指示，程序被读取并被写入计算机2000的RAM 2004中。因此，CPU 2003顺序读取并执行程序，从而能够实现根据本发明的图像处理设备的功能。

下文中，将参照图12中所示的流程图详细说明图10中所示的各构成部件的操作。

图12是示出要由根据本发明的第五实施例的图像处理设备执行的整个处理的流程图。

(步骤S100)

首先，响应预定的输入指示，作为相减处理的目标的多个图像(时间序列图像)被输入到图像输入单元1。然后，图像输入单元1对输入的图像执行缩减(reduction)处理，并且将处理的图像输出到偏移向量计算单元2。顺便提及，从其提供要被输入到图像输入单元1的图像的图像输入设备对应于存储介质(比如为直接或间接地与计算机2000相连的硬盘2002、磁光盘2007等)或图像产生设备1000。

此外，用于将时间序列图像输入到图像输入单元1的预定的输入指示通过操作相关图像处理设备的用户的操作，或通过控制相关图像处理设备的控制器(未示出)等被给出。

顺便提及，例如，如果水平方向的像素数和垂直方向的像素数按照要被应用到输入图像的缩减处理中的缩减率分别被设置以具有1/4×1/4尺寸，则从增加处理效率以保持差图像的必要分辩率的方面来说，该尺寸是所希望的。然而，缩减率不一定是上述值，即，当然也可使用另外的缩减率。

此外，可以将图像输出到后面描述的图像处理单元9，而不用执行缩减处理(即，具有不变的尺寸)。通过这样的处理，可以在保持高分辩率的同时，将作为要被解释或读取的目标的图像与其中变化被强调的差图像相比较。

这里，尽管时间序列图像通常表示针对同一患者并且在不同时间点拍摄的一组图像，但在本实施例中，为了简化说明，仅一对第一和第二图像(IM1和IM2)被看作时间序列图像。然而，即使图像数为三个或更多，从三个或更多个图像中选择一对图像作为一对，并且仅对每对图像应用后面所述的处理。

在图10中，第一图像IM1和第二图像IM2共同构成要从图像输入单元1被输入的时间序列图像。在本实施例中，假定第一图像IM1为晚些拍摄的图像，第二图像IM2为第一图像IM1被拍摄之前拍摄的图像。顺便提及，在下面的说明中，假定图像IM1和IM2为胸部X射线正面图像。然而，本发明不限于此。即，在不脱离其目的的情况下，本发明当然也可适用于其它种类的图像。

此外，例如，第一图像IM1和第二图像IM2为通过数字射线摄影系统产生的图像，所述数字射线摄影系统采用使用平板检测器或光激励荧光体的CR(计算机射线摄影)等，并且由这些图像所指示的图像数据的特性与拍摄对象时所获取的相关X射线量的对数成比例。

换句话说，当解释由数字射线摄影系统拍摄的医学图像时，通常执行非线性灰度转换处理，以便符合已被传统使用的银盐胶片的特性。然而，在本发明中，使用在上述灰度转换处理被执行之前所获取的图像数据。

一方面，第一图像IM1和第二图像IM2被输出到图像处理单元9。因此，图像处理单元9对这些图像执行上述的灰度转换处理，产生最适于解释的灰度转换图像，然后，将产生的图像输出到输出单元8。这里，要由图像处理单元9执行的灰度转换处理可以是具有如图13中所示的转换特性曲线的处理，其类似于传统的银盐胶片。此外，也可附加地执行诸如钝化掩膜(unsharp mask)处理等之类的频率加重处理。

在任何情况下，本发明也适用于通过图像处理单元9灰度转换处理被事先执行的情况，灰度转换图像分别与处理前的图像相关或相关联，相关的图像被存储在预定的存储设备(例如硬盘2002)中，此后，当第一图像IM1和第二图像IM2被指定时，存储的图像与指定的图像一起被读取。

(步骤S200)

偏移向量计算单元2计算并输出表示输入的第一图像IM1和第二图像IM2之间的对应像素的物理关系的偏移向量。这里，将参照图14说明偏移向量计算单元2的详细构成，并且将参照图15说明偏移向量计算单元2的详细操作。

即，图14是示出根据本发明的第五实施例的偏移向量计算单元的详细构成的方框图，图15是示出根据本发明的第五实施例的偏移向量计算单元的详细操作的流程图。

(步骤S21)

从图像输入单元1输入的第一图像IM1和第二图像IM2的每个的尺寸通过图像缩减单元21被缩减到预定尺寸。在本实施例中，每个图像的水平和垂直方向的尺寸被缩减到1/4×1/4，以产生第一缩减图像IM1和第二缩减图像IM2。然后，这些缩减图像分别被输入到胸腔(胸廓)检测单元22。同时，尺寸未被缩减的第一图像IM1和第二图像IM2分别被输入到第二ROI设置单元25。

因此，应当注意，尺寸为原始图像尺寸1/16×1/16的图像被输入到胸腔检测单元22，尺寸为原始图像尺寸1/4×1/4的图像被输入到第二ROI设置单元25。

(步骤S22)

胸腔检测单元22分析第一缩减图像IM1和第二缩减图像IM2，检测在各图像中对象的交点，将检测的交点设置为参考点，然后，基于设置的参考点输出在两个图像之间的对象的大致的失配(或位移)量。

图16A示出根据后面描述的方法在第一缩减图像IM1中被检测的参考点Lm11-Lm18，图16B示出根据后面描述的方法在第二缩减图像IM2中被检测的参考点Lm21-Lm28。

这些参考点基于包括在对象图像中的固有特性被确定。例如，在所示的胸部X射线正面图像中，肺的顶点(Lm11、Lm15、Lm21、Lm25)、CP(肋膈)角(Lm14、Lm18、Lm24、Lm28)、胸腔(或胸廓)的外边缘(Lm12、Lm13、Lm16、Lm17、Lm22、Lm23、Lm26、Lm27)被检测，并且包含这些参考点的矩形区R1和R2然后被检测。

更具体地说，CP(肋膈)角对应于胸腔的边缘外侧与隔膜的阴影相交的部分(即对应于图16A的第一图像中的Lm14和Lm18的部分)。

在本实施例中，胸腔检测单元22设置多个矩形区(分析区)，用于从要被处理的图像中检测参考点，然后，根据包括在多个矩形区中的图像数据实际检测多个参考点。图17A示出对于胸部正面图像用于检测参考点的矩形区被设置的状态。下面，将说明如何检测参考点。

然后，胸腔检测单元22通过在垂直方向累积要被处理的图像数据产生轮廓(profile)数据PM，然后将产生的轮廓数据PM存储在内部存储器比如RAM 2004等中。图17B示出轮廓数据PM的示例。在图17B中，像素值变为最大的位置的水平坐标Mx被设置为图17A中的中心线M的水平坐标。

随后，胸腔检测单元22将图17A中所示的多个分析区Rt1-Rt8、Rr1-Rr4和Rc1-Rc8分别设置在基于中心线M近似对称的位置。这里，假定各区的尺寸基于要被拍摄的对象的平均尺寸被预定。

之后，对于分析区Rt1-Rt8的每个，胸腔检测单元22通过沿水平方向累积要被处理的图像数据产生轮廓数据Pt，并且将产生的轮廓数据Pt存储在内部存储器比如RAM 2004等中。图17C示出轮廓数据Pt的示例。在图17C中，线yt对应于胸部正面图像中肺野的上边缘，并且胸腔检测单元22检测其中轮廓数据Pt第一个负变化的位置作为垂直方向的上边缘的候选位置yt。

随后，对每个区执行与上面一样的处理，并且对于位于中心线M两侧的四个区Rt1-Rt4和四个区Rt5-Rt8的每个执行二维内插，由此图16A中所示的点Lm11和Lm15被确定。此外，通过对第二缩减图像IM2执行类似的处理，图16A中所示的点Lm21和Lm25被确定。

下面，胸腔检测单元22解析分析区Rc1-Rc8，从而确定图16A中所示的点Lm14和Lm18。各区的图像数据沿水平方向被累积，按如上所述的同样方式，轮廓数据Pc被产生，并且产生的轮廓数据Pc被存储在内部存储器比如RAM 2004等中。图17D示出轮廓数据Pc的示例。在图17D中，线yc对应于肺野和横膈膜之间的边界。然后，胸腔检测单元22分析轮廓数据Pc，并且确定平均亮度级变化最大的边界yc。在这样的确定中，轮廓数据Pc的主差值首先被分析，并且仅需检测分析值变化最大的位置。

对每个区执行同样的处理，并且其中边界yc具有最大值的区被检测。然后，检测的区进一步被分成精细区(点)，并且边界yc具有最大值的点被设置成肺野的最低边缘(Lm14和Lm18)。

顺便提及，如果区域从肺野脱离，则图17D的轮廓数据Pc不会很大地变化。由于该原因，胸腔检测单元22将预定阈值与上述的主差值进行比较。然后，如果不存在超过阈值的差值，则认为在相关区域的肺野和横膈膜之间不存在边界，从而该区从要被处理的目标中被排除。

通过对第一缩减图像IM1和第二缩减图像IM2执行类似的处理，图16A中所示的点Lm14、Lm18、Lm24和Lm28被检测。

下面，胸腔检测单元22通过沿垂直方向累积区域Rr1中的图像数据产生轮廓数据Pr，并且将产生的轮廓数据Pr存储在内部存储器比如RAM 2004等中。图17E示出轮廓数据Pr的示例。如图17E中所示，在区域Rr1和Rr2中，胸腔检测单元22将在最左侧轮廓数据Pr具有最大值处的坐标xc设置为点Lm11和Lm13的临时水平坐标值。同时，在区域Rr3和Rr4中，胸腔检测单元22将在最右侧轮廓数据Pr具有最大值处的坐标xc设置为点Lm26和Lm27的临时水平坐标值。

在检测在该处轮廓数据具有最大值的位置的情况中，如果封闭相关位置的轮廓数据的值不在预定范围内，则相关的最大值被排除作为对象和背景之间的边界。在避免错误地检测对象的边界方面，希望这样做。

此外，胸腔检测单元22将各区域的中心点的垂直坐标设置为点Lm14、Lm18、Lm24和Lm28的临时垂直坐标。然后，根据由此获取的临时参考点Lm14’、Lm18’、Lm24’和Lm28’设置更小的区域，并且对这些区域执行相同的处理，从而确定最终参考点Lm14、Lm18、Lm24和Lm28。

顺便提及，本发明不限于上述方法。即，例如对应于USP 5790690的日本专利申请公开No.H08-335271所公开的其它方法也可适用。

这里，在日本专利申请公开No.H08-335271公开的方法中，多个一维轮廓数据从图像数据中被捕获并被分析以检测胸腔(或胸廓)。此外，上述方法的特性在于，基于从部分轮廓数据获取的特征点，考虑肺野的整个形状检测整个肺野。

(步骤S23)

胸腔检测单元22确定区域R1和R2，其针对由上述方法检测的参考点被限定并且作为要被进行相减处理的目标。这里，应当注意，在下文中，区域R1和R2也分别被称作相减处理目标区R1和R2。图16A示出分别关于第一缩减图像IM1和第二缩减图像IM2被设置的区域R1和R2。如图16A中所示，每个区域被确定为最大矩形区，其包括各自图像中检测的参考点。然后，胸腔检测单元22将区域R1的左上坐标和右下坐标分别输出给第一ROI设置单元23和第二ROI设置单元25。

随后，胸腔检测单元22计算各区域R1和R2的中心坐标(xc1，yc1)和(xc2，yc2)，通过方程式(10)计算第一缩减图像IM1和第二缩减图像IM2之间的整体偏移(或位移)量G，并且将计算的失配量G输出到第一ROI设置单元23。

$G=\left[\begin{array}{c}{G}_x\\ G_y\end{array}\right.==\left[\begin{array}{c}\mathrm{xc}2-\mathrm{xc}1\\ \mathrm{yc}2-\mathrm{yc}1\end{array}\right]...\left(10\right)$

而且，胸腔检测单元22通过方程式(11)，根据图16A中所示的区域R1和R2的左上和右下坐标计算包括在第一缩减图像IM1和第二缩减图像IM2中的对象的放大率变化M，并且将计算的放大率变化M输出到第一ROI设置单元23。

$M=\left[\begin{array}{c}{M}_x\\ M_y\end{array}\right]=\left[\begin{array}{c}\frac{x_{L2}-x_{U2}}{x_{L1}-x_{U1}}\\ \\ \frac{y_{L2}-y_{U2}}{y_{L1}-y_{U1}}\end{array}\right]...\left(11\right)$

(步骤S24)

根据通过胸腔检测单元22输入的相减处理目标区R1的坐标、整体偏移量G和放大率变化M，第一ROI设置单元23关于第一缩减图像IM1和第二缩减图像IM2设置多个ROI。

为了简化起见，图16B仅从关于第一缩减图像IM1和第二缩减图像IM2设置的ROI中示出左上的三个ROI。在图16B中，尺寸为Wt×Ht的矩形ROI关于第一缩减图像IM1被设置，尺寸为Ws×Hs的矩形ROI相对于第二缩减图像IM2被设置。这里，关于第一缩减图像IM1被设置的ROI被称作模板ROI，关于第二缩减图像IM2被设置的ROI被称作搜索ROI。

第一ROI设置单元23排列第一缩减图像IM1中的第一模板ROI的中心，以便与对应于相减处理目标区R1的左上方的位置重叠。随后，第一ROI设置单元23以水平间隔Δx和垂直间隔Δy将模板ROI设置在覆盖整个区域R1的范围内。

下面，根据设置在第一缩减图像IM1中的模板ROI的位置、从胸腔检测单元22输入的整体偏移量G和放大率变化M，第一ROI设置单元23设置搜索ROI。即，如果假定对应的第n个模板ROI的中心坐标为(x_n ^t1，y_n ^t1)，搜索ROI的中心坐标为(x_n ^s1，y_n ^s1)，则第n个搜索ROI的中心坐标通过方程式(12)和方程式(13)被计算。

x_n ^s1＝x_n ^t1+mod(n，C)M_xΔx+G_x    (n＝0，...，N-1)

                                              ...(12)

y_n ^s1＝y_n ^t1+floor(n/C)M_yΔy+G_y   (n＝0，...，N-1)

                                              ...(13)

这里，应当注意，符号N表示模板ROI和对应的搜索ROI的数目，并且数目N根据可被设置在第一缩减图像IM1的差处理目标区R1中的模板ROI的数目被确定。符号C表示水平方向的ROI的数目，符号mod()表示模运算，符号floor表示“地板”函数。

模板ROI的尺寸、搜索ROI的尺寸、水平间隔Δx和垂直间隔Δy被事先确定并被存储在例如内部存储器，比如偏移向量计算单元2的RAM 2004等中。顺便提及，也可根据要被进行相减处理的一种对象事先选择适当的值作为这些值。

例如，在目标为胸部正面图像的情况下，最好大约将模板ROI的尺寸设置为25×20(mm)，将搜索ROI的尺寸设置为25×25(mm)，将间隔设置为3(mm)。然而，本发明不限于这些值，即，其它值当然也可适用于本发明。而且，ROI的形状不限于矩形，即，其它形状当然也可适用于ROI。

第一ROI设置单元23将按上面设置的模板ROI和搜索ROI的位置和尺寸输出到第一匹配单元24。

(步骤S25)

第一匹配单元24对包括在输入的模板ROI和对应的搜索ROI中的第一缩减图像IM1和第二缩减图像IM2的数据执行匹配处理，计算匹配度为最高处的位置作为关于每一组ROI的偏移向量，然后输出计算的结果。

顺便提及，应当注意，在匹配处理中可使用各种公知的方法。例如，可使用序列类似性检测算法和互相关法等，其在由东京大学出版社出版的“图像分析手册”一书中被描述。在第五实施例中，第一匹配单元24和第二匹配单元26通过互相关法计算偏移向量。然而，由于在上述文件中已描述了互相关法的详细内容，因此这里对其的说明将被省略。

图18示出了通过对于每一组模板ROI和搜索ROI执行匹配处理以获取偏移向量以及将获取的偏移向量叠加在第一缩减图像IM1上所获取的图像。而且，尽管图18通过示例方式示出25个偏移向量，但偏移向量的实际数目取决于模板ROI的尺寸、设置的间隔和相减处理目标区的尺寸。每一偏移向量指示包括在第一缩减图像IM1中的对象的结构对应于相应的第二缩减图像IM2的部分，即，每一偏移向量指示在这两个图像之间的对象的失配量。

第一匹配单元24将作为第一偏移向量V_n ¹(n＝0，...，24)的偏移向量输出到第二ROI设置单元25。

(步骤S26)

根据从第一匹配单元24输入的第一偏移向量V_n ¹和从胸腔检测单元22预先输入的相减处理目标区R1，第二ROI设置单元25将模板ROI和搜索ROI分别设置到第一图像IM1和第二图像IM2。

在第一ROI设置单元23的处理中，模板ROI和搜索ROI根据相减处理目标区的位置和尺寸被设置。然而，在第二ROI设置单元25中，ROI根据第一偏移向量V_n ¹被设置。

首先，根据事先通过图像缩减单元21所使用的缩减倍率来缩减第一图像IM1和第二图像IM2，第二ROI设置单元25将相减处理目标区R1的坐标和从第一匹配单元24输入的第一偏移向量V_n ¹转换为第一图像IM1或第二图像IM2中的值。

在第五实施例中，第一缩减图像IM1和第二缩减图像IM2分别等于垂直方向和水平方向的第一图像IM1和第二图像IM2的1/4。因此，第一偏移向量V_n ¹和相减处理目标区R1的坐标被分别放大四倍，由此获取转换的第一偏移向量V_n ¹’和相减处理目标区R1的转换坐标R1’。这里，应当注意，相减处理目标区R1的转换坐标R1’也被称作相减处理目标区R1’。

接下来，以与第一ROI设置单元23所使用的方式相同的方式，第二ROI设置单元25关于第一图像IM1设置模板ROI。此时，相减处理目标区为相减处理目标区R1’，并且ROI的尺寸和设置间隔仅相对于由第一ROI设置单元23设置的各值被放大四倍。然而，从增加处理效率和抑制由于不必要的匹配出现的误差的观点来看，最好将ROI的尺寸设置为小于第一模板ROI中设置的值。

随后，第二ROI设置单元25根据方程式(14)和(15)关于第二图像IM2排列搜索ROI。即，如果假定对应的第n个模板ROI的中心坐标为(x_n ^t2，y_n ^t2)，搜索ROI的中心坐标为(x_n ^s2，y_n ^s2)，并且转换的第n个第一偏移向量为v_n ¹’＝(x_n ¹’，y_n ¹’)，则第n个搜索ROI的中心坐标位置通过方程式(14)和(15)被计算。

x_n ^s2＝x_n ^t2+x_n ¹′(n＝0，...，N-1)     ...(14)

y_n ^s2＝y_n ^t2+y_n ¹′(n＝0，...，N-1)     ...(15)

第二ROI设置单元25将按上面设置的模板ROI和搜索ROI的位置和尺寸输出到第二匹配单元26。

(步骤S27)

第二匹配单元26对包括在输入的模板ROI和对应的搜索ROI中的第一图像IM1和第二图像IM2的数据执行匹配处理，计算匹配度为最高处的位置作为对于每组ROI的偏移向量V_n ²，然后将计算的结果输出到过滤(filter)单元3(图10)。这里，第二匹配单元26的匹配处理与第一匹配单元24的匹配处理相同，因此对其的说明将被省略。

顺便提及，在第五实施例中，第一匹配单元24和第二匹配单元26的每个均通过互相关法计算偏移向量。然而，第一匹配单元24和第二匹配单元26的任何一个也可使用另外的方法。例如，如果第一匹配单元24使用序列类似性检测算法而第二匹配单元26使用互相关法，则有可能在将保持一定程度的匹配的精度同时缩短整个处理时间。

随后，将参照图12和19中所示的流程图说明过滤单元3的详细操作。

(步骤S300)

过滤单元3对输入的第二偏移向量V_n ²执行过滤处理，以消除在匹配中混合的误差分量，并且输出误差分量已从中被消除的第三偏移向量V_n ³。在下文中，将参照图19中所示的流程图说明过滤单元3的详细操作。

图19是示出根据本发明的第五实施例的过滤单元的详细操作的流程图。

(步骤S31)

过滤单元3对于输入的偏移向量确定过滤处理的过滤处理间隔。在本实施例中，处理间隔等于整个图像，即，处理间隔等于包括如图18中所示的所有25个偏移向量的间隔。

(步骤S32)

过滤单元3初始化要被用于后面所述处理的每一偏移向量的权重系数。在本实施例中，通过第二匹配单元26在每一ROI中计算的互相关系数CC_n被用作初始值。

(步骤S33)

过滤单元3对输入的第二偏移向量V_n ²＝(x_n ²，y_n ²)的水平分量x_n ²执行过滤处理。顺便提及，在下列步骤S331-S335中将说明过滤处理的详细内容。

(步骤S331)

过滤单元3仅将水平分量x_n ²暂时存储在内部存储器比如RAM2004中。图20示出了此时偏移向量的水平分量x_n ²。然后，过滤单元3通过使用先前输入的互相关系数CC_n作为权重对偏移向量分量执行加权平均处理(或加权平均值处理)，并且将处理结果x_n ²’暂时存储在内部存储器比如RAM 2004中。

(步骤S332)

过滤单元3通过方程式(16)计算在加权平均处理之前的偏移向量分量和在加权平均处理之后的偏移向量分量之间的残差r_n。

r_n＝x_n ²-x_n ²’      ...(16)

而且，残差的中值绝对偏差MAD通过方程式(17)被计算。这里，应当注意在方程式(17)中符号“median(x)”为“x”的中值。

MAD＝median(|r_n|)           ...(17)

(步骤S333)

过滤单元3根据方程式(18)计算加权平均处理的权重系数w_n。

$w_n=\left\{\begin{array}{cc}{\mathrm{CC}}_n{(1-{\left(\frac{r_n}{\mathrm{\&kappa;MAX}}\right)}^2)}^2& \left|r_n\right|<\mathrm{\&kappa;MAD}\\ 0& \left|r_n\right|\&GreaterEqual;\mathrm{\&kappa;MAD}\end{array}\right....\left(18\right)$

这里，符号k表示根据残差的分布先前确定的值。在本实施例中，例如，k＝6。然而，本发明不限于此，也可使用其它的值。

(步骤S334)

过滤单元3通过使用在先前步骤中更新的权重系数w_n对偏移向量的水平分量x_n ²再次执行加权平均处理。

(步骤S335)

过滤单元3判断步骤S332和其后步骤中的处理是否被执行了等于预定的次数T的次数。如果判断相关的处理被执行了预定的次数T(即步骤S335中的“是”)，则流程进到下一步骤。同时，如果判断相关的处理未被执行预定的次数(即步骤S335中的“否”)，则流程返回到步骤S332。这里，预定的次数T可根据过滤处理的处理效率和效果的平衡被确定。然而，最好将T设置为2上下。

(步骤S34)

过滤单元3对偏移向量的垂直分量y_n ²执行如上所述的相同处理。这里，由于步骤S34中处理的内容与涉及水平分量x_n ²的步骤S33中的处理内容基本上相同，因此对其的说明将被省略。

(步骤S35)

过滤单元3通过上述方法将对于每一水平和垂直方向已被进行过滤处理的第三偏移向量V_n ³＝(x_n ³，y_n ³)输出到样条内插单元4。

然后，将再次描述图12。

(步骤S400)

样条内插单元4通过使用公知的三次样条内插对第三偏移向量V_n ³＝(x_n ³，y_n ³)执行内插处理，并且将内插处理中获取的第四偏移向量V_n ⁴＝(x_n ⁴，y_n ⁴)输出到配准单元5。在该处理中，样条内插单元4对第三偏移向量V_n ³的条的每一水平和垂直分量应用三次样条内插。

顺便提及，可通过例如1990年IEEE计算机学会出版的G.Wolberg所著“数字图像卷绕(Digital Image Warping)”等中所述的方法执行三次样条内插，并且该方法是公知的，因此对其的详细描述将被省略。

图21示出了第三偏移向量V_n ³的某一条的水平分量的状态。在图21中，符号K表示在执行内插前的水平方向的偏移向量数目，符号L表示在执行内插后的偏移向量数目。顺便提及，可事先确定一值作为数L，通过该值可在后面描述的差图像产生中实现充分的图像质量。例如，如果数L为数K的百倍，就足够了。

(步骤S500)

配准单元5通过利用从样条内插单元4输入的第四偏移向量V_n ⁴使第二图像IM2变形，产生已与第一图像IM1配准的卷绕图像IM3，并且随后将产生的卷绕图像IM3输出到相减操作单元6。

即，第四偏移向量V_n ⁴指示分别示出在第一图像IM1和第二图像IM2上的公共构成之间的精确的失配(位移)量。由于该原因，通过反向地将第四偏移向量V_n ⁴应用到卷绕的第二图像IM3，第一图像上对应的位置被计算，并且通过根据内插处理执行抽样可确定像素值。

在任何情况下，由于在上述文件中已描述了这样一种方法的详细内容，因此对其的详细说明将被省略。顺便提及，作为一种在执行重抽样情况下的内插方法，考虑到差图像的图像质量，最好使用例如双线性内插法等。

(步骤S600)

相减操作单元6在输入的第一图像IM1的像素和变形的第二图像IM3的像素之间执行相减以产生差图像IMS，并且将产生的差图像IMS输出到后处理单元7。

(步骤S700)

后处理单元7对输入的差图像IMS执行灰度转换处理以产生其像素值已被转换到适于显示的像素范围内的差图像IMS’，并且将产生的差图像IMS’输出到输出单元8。这里，灰度转换可根据相减处理之前图像和相减处理之后图像的精度来确定。

例如，如果在相减处理被执行之前的图像的像素值的精度为无符号的12位，则通过相减处理产生具有带符号的13位大小的差图像。这里，如果输出单元8可显示的像素值的大小为无符号的8位，则其必须将13位大小转换为8位大小。因此，对对于对应于X射线图像的对象X射线照射量具有线性灰度特性的图像从输出单元8被输出。

或者，差图像的直方图被计算，并且，基于对应于直方图中的模值的像素值，预定的大小可被线性地转换为8位大小。顺便提及，应当注意，在输入和输出范围之间的这种转换函数不一定总是线性的。例如，可使用特性曲线为如图13中所示的函数。

(步骤S800)

输出单元8以适于解释内插的形式显示输入的灰度转换的差图像IMS’连同对原始的第一图像IM1和第二图像IM2进行适合的灰度转换的图像。例如，第一图像IM1和第二图像IM2以及差图像IMS’被对照地显示。

输出单元8例如为显示设备，比如CRT监视器、LCD等。然而，显示设备不限于这样的电子显示设备，即，输出单元8也可以是将图像输出为硬拷贝的图像输出设备，比如激光成像器等。

而且，在本发明中，输出单元8可以不一定为显示设备，即，与计算机连接的硬盘、网络输入/输出设备等也可被应用于输出单元8。

如上所述，根据第五实施例，在过滤单元3中，指示在输入的第一图像IM1和第二图像IM2之间的对应像素的位置关系的偏移向量通过使用权重系数被进行过滤处理，所述权重系数根据不易受离群值(outlier)影响的中值(即，根据对偏移向量的加权平均处理被执行之前和之后的状态之间的差的中值绝对偏差)来确定，并且获取的偏移向量通过样条内插被内插，由此执行配准。

因此，在后阶段的样条内插中，可以抑制由匹配所引起的误差的影响，由此可以实现不易受部分离群值影响的高精度的配准。

根据本实施例，可以保证配准的精度以及实现处理的效率，由此可以产生高质量的差图像。

<第六实施例>

在第五实施例中，第二图像IM2仅被变形一次。然而，本发明不限于此。即，如第六实施例中将说明的，可以使第二图像IM2变形多次并根据这样的次数改变偏移向量的过滤处理的处理内容。

图22是示出要由根据本发明的第六实施例的图像处理设备执行的整个处理的流程图。

顺便提及，在第六实施例中，整个图像处理设备的构成与图10中所示的图像处理设备的构成相同，并且控制单元(未示出)控制整个处理。

与第五实施例中的图12的流程图相比，在第六实施例的图22的流程图中增加了步骤S900。因此，控制单元(未示出)判断在步骤S500中的卷绕处理是否被执行了预定次数。如果判断卷绕处理未被执行预定次数(即步骤S900中的“否”)，则流程返回到步骤S200。与此同时，如果判断卷绕处理被执行了预定次数(即步骤S900中的“是”)，则流程进到步骤S600。

此外，在第六实施例中，在第二和接下来执行的卷绕处理中，涉及步骤S200的详细操作的图15中的步骤S21-S25的处理被省略。然后，在步骤S26的第二ROI设置中，通过使用第三偏移向量V_n ³＝(x_n ³，y_n ³)，该ROI通过方程式(19)和(20)被设置。

x_n ^s2＝x_n ^t2+x_n ³  (n＝0，...，N-1)      ...(19)

y_n ^s2＝y_n ^t2+y_n ³  (n＝0，...，N-1)      ...(20)

在这种情况下，如先前所述，不用说偏移向量的尺寸不根据图像的缩减尺寸而改变。

而且，在第六实施例中，在步骤S300的过滤处理中，过滤单元3根据方程式(21)计算过滤处理的权重系数w_n。

$w_n=\left\{\begin{array}{cc}{\mathrm{CC}}_n{(1-{\left(\frac{r_n}{{\mathrm{\&kappa;}}_i\mathrm{MAX}}\right)}^2)}^2& \left|r_n\right|<{\mathrm{\&kappa;}}_i\mathrm{MAD}\\ 0& \left|r_n\right|\&GreaterEqual;{\mathrm{\&kappa;}}_i\mathrm{MAD}\end{array}\right....\left(21\right)$

这里，符号i表示在相应处理中卷绕的次数。即，相对于中值绝对偏差MAD的系数k根据重复处理的次数被改变。此时，最好根据的重复处理的次数的增加使系数k变小。

例如，在第一重复处理中，与第五实施例一样，系数k被设置为6，并且在第二重复处理和后面的重复处理中，该系数被设置为更小，比如4或2。因此，在计算权重系数的情况下，离群值的容差变得更小。

通常，随着卷绕被重复，要被进行配准的两个图像之间的失配(位移)量变得较小。因此，根据这一点，如果在过滤处理中离群值的容差被设置得较小，则能够实现更有效的过滤处理。

如上所述，根据第六实施例，除了上述第五实施例的效果以外，还可以实现比第五实施例中的过滤处理更有效的过滤处理。

<第七实施例>

在上述第五和第六实施例中，过滤单元3通过方程式(18)计算对于在整个图像上延伸的偏移向量的权重系数。然而，通过计算对于图像的每部分的权重系数，也可以执行过滤处理。

在第七实施例中，图像输入单元1通过分析第一图像IM1和第二图像IM2确定多个局部区，然后将确定的局部区的位置信息输出到过滤单元3。图23是用于描述根据本发明的第七实施例的区域分割的示例的示图。在图23中，符号A1-A6分别表示基本上对应于根据对象的结构的下列区域的矩形区域。即，区域A1和A2对应于非对象部分，区域A3和A4对应于肺野部分，区域A5对应于纵隔部分，区域A6对应于腹腔部分。

通过根据除区域A1和A2以外的每一区域改变方程式(18)的系数k的值，过滤单元3计算权重系数。例如，系数k＝6被设置到区域A3和A4的肺野部分，系数k＝4被分别设置到区域A5的纵隔部分和区域A6的腹腔部分。

在区域A5的纵隔部分和区域A6的腹腔部分，原始图像的对比度通常较低。为此，即使相关系数本身较高，偏移向量也可能变化。因此，在过滤单元3中，通过使离群值的容差对于这些区域变小，可以有力地抑制偏移向量中的噪声分量。

顺便提及，过滤处理可以不根据上述的对象的固有区域被分割，而可以对于偏移向量的垂直和水平线的一个被独立地执行。

如上所述，根据第七实施例，除了上述第五实施例的效果之外，通过使过滤单元3对于要被处理的图像中的特性部分区的每个执行过滤处理，以实现更适合和更有效的过滤处理。

<第八实施例>

在上述第五和第六实施例中，对于在整个图像上延伸的偏移向量执行一个三次样条内插。然而，也可以对于多个分割的区域的每个执行样条内插。

这里，区域可根据如图23中所示的构造被分割。此外，如图24中所示，可以通过使用关于每四点的偏移向量分量执行三次样条内插，并且输出一部分处理区域。

在图24中，使用偏移向量的水平分量x₁ ²至x₄ ²执行三次样条内插，并且内插结果so₂的分量x₂ ²和x₃ ²之间的结果被给定作为该间隔的最终结果。然后，随着按一个抽样移动窗口，相同的处理被执行，由此所有间隔的内插结果被产生并被输出。顺便提及，对于数据两端仅必须使用常规结果。

如上所述，根据第八实施例，除了上述第五实施例的效果以外，还可以抑制内插结果的振动，即使在图像等中的偏移向量的变化较大(因为其包括许多噪声分量)等情况下，由此可以保证配准的精度。

<第九实施例>

在上述第五-第八实施例中，互相关系数CC_n被用作方程式(18)所示的加权平均处理中的权重系数。然而，也可采用其它方法。

即，在第九实施例中，偏移向量计算单元2计算包括在涉及偏移向量计算的ROI中的图像数据的纹理，并且将计算的纹理输出到过滤单元3。这里，尽管可使用各种纹理，但在本实施例中ROI中像素值的变化被用作指标(即纹理的复杂度)。即，在本实施例中，偏移向量计算单元2形成模板ROI的直方图，并且将频率为非零的图像像素值数

设置为纹理。

典型地，在胸部正面图像中，纵隔部分、心脏部分、横膈膜部分等的对比度较低，由此这些部分的纹理较低。相反，设置在肺野部分和其周围部分的模板ROI的纹理较高。因此，很有可能在低纹理ROI中计算的偏移向量包括较多误差，由此最好降低偏移向量的过滤处理的权重。

由于该原因，在第九实施例中，过滤单元3根据方程式(22)和(23)计算加权平均处理的权重系数w_n。

$w_n=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_n{(1-{\left(\frac{r_n}{\mathrm{kMAX}}\right)}^2)}^2& \left|r_n\right|<\mathrm{kMAD}\\ 0& \left|r_n\right|\&GreaterEqual;\mathrm{kMAD}\end{array}\right....\left(22\right)$

${\mathrm{\&Phi;}}_n=\frac{{\mathrm{\&phi;}}_n\&CenterDot;{\mathrm{CC}}_n}{\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_n\&CenterDot;{\mathrm{CC}}_n}...\left(23\right)$

这里，应当注意，方程式(23)计算的值φ_n表示归一化的纹理。

如由方程式(18)所示的，用于偏移向量的计算的互相关系数被用作第五实施例中权重的指标(或量度)。然而，如上所述，互相关系数指示较高值，即使对于纹理贫乏的部分，由此有可能较大的权重被给定到包括较大误差的偏移向量。

与此同时，根据第九实施例，ROI中的纹理(即纹理的复杂度)被用作权重的指标(或量度)，由此较小的权重被给定到位于比如相关信息量较小的纵隔部分、腹腔部分等部分中的偏移向量。因此，可以增加过滤处理的精度。

顺便提及，本实施例中的纹理不限于在上述直方图中频率为非零的像素值数。即，也可使用其它的指标(或量度)。例如，在要被处理的ROI中的图像数据被进行频率变换的情况下，可使用频率变换的图像数据的高频分量的比率。此外，也可使用像素值的离差、直方图的峰度等。

如上所述，本发明的实施例被详细说明。顺便提及，本发明可被应用于例如系统、设备、方法、程序、存储介质等。更具体地说，本发明可被应用于由多个仪器组成的系统或应用于由单一仪器组成的设备。

在任何情况下，本发明也适用于这样一种情况，即用于实现上述实施例的功能的软件的程序(即在上述实施例中，对应于后附流程图的程序代码)被直接或间接地提供给系统或设备，并且系统或设备中的计算机根据提供的实现上述实施例的程序操作各种设备。

就此而言，由于上述实施例的功能通过计算机来实现，因此安装在相应计算机中的程序代码本身实现本发明。即，用于实现本发明的功能处理的计算机程序本身也包括在本发明的概念中。

在这种情况下，可使用目标代码、要由解释器执行的程序、要被提供给OS的脚本等，只要其具有程序功能。

这里，作为用于提供程序的存储介质，可使用例如软盘、硬盘、光盘、磁光盘、MO、CD-ROM、CD-R、CD-RW、磁带、非易失性存储卡、ROM、DVD(DVD-ROM、DVD-R)等。

此外，作为提供程序的方法，存在一种方法，其通过使用客户计算机的浏览器与因特网上的主页连接，并将本发明的计算机程序本身或者包括自动安装功能连同计算机程序的压缩文件下载到比如硬盘等的记录介质中。此外，还存在一种方法，其将构成本发明的程序的程序代码分成多个文件，并且从不同的主页下载各文件。即，用于利用计算机将实现本发明的处理的功能的程序文件下载给多个用户的WWW服务器包括在本发明的范围内。

而且，可以加密本发明的程序，将加密的程序存储在比如CD-ROM等的存储介质中，将获得的存储介质分发给用户，使已满足预定条件的用户能够通过因特网从主页下载用于解密加密程序的密钥信息，使该用户能够将解密的程序安装到适当的计算机中，并因此实现本发明的功能。

而且，本发明不仅包括通过执行由计算机读取的程序代码实现上述实施例的功能的情况，也包括运行在计算机上的OS等根据程序代码的指令执行部分或所有实际处理的情况，由此上述实施例的功能通过该处理被实现。

此外，在从存储介质读取的程序被一次写入设置在功能扩展板(其被插在计算机或与计算机连接的功能扩展单元中)中的存储器中，并且随后设置在功能扩展板或功能扩展单元中的CPU等根据程序指令执行部分或所有实际处理的情况下，也可实现上述实施例的功能。

本申请对2004年6月8日申请的日本专利申请No.2004-170231和2005年1月28日申请的日本专利申请No.2005-021827要求优先权，其整个内容以参照的方式被包含在这里。

Claims (12)

1、一种图像处理设备，其输出第一图像和第二图像之间的差图像，包括：

输入单元，适于输入第一图像和第二图像；

偏移向量计算单元，适于将多个感兴趣区域分别设置到第一图像和第二图像，并且对于每个感兴趣区域计算指示第一图像和第二图像之间的失配量的偏移向量；

过滤单元，适于对偏移向量执行过滤处理；

内插单元，适于对通过所述过滤单元进行过滤处理的偏移向量执行内插；

配准单元，适于根据通过所述内插单元内插的偏移向量配准第一图像和第二图像；

相减操作单元，适于在各配准的图像上的相应像素之间执行相减操作；以及

输出单元，适于输出通过所述相减操作单元获取的差图像。

2、根据权利要求1的图像处理设备，其中，所述输入单元输入对于对象均具有X射线照射量的线性灰度特性的图像，作为第一图像和第二图像。

3、根据权利要求1的图像处理设备，其中，所述偏移向量计算单元通过互相关法计算感兴趣区域之间的偏移向量。

4、根据权利要求3的图像处理设备，其中，所述过滤单元根据执行加权平均处理之前和之后的状态之间的差的中值绝对偏差，通过执行预定次数的加权平均处理对偏移向量执行过滤处理。

5、根据权利要求4的图像处理设备，其中，所述过滤单元根据由所述偏移向量计算单元计算的互相关系数对偏移向量执行加权平均处理。

6、根据权利要求4的图像处理设备，其中，所述过滤单元根据感兴趣区域中纹理的复杂度对偏移向量执行加权平均处理。

7、根据权利要求1的图像处理设备，其中，所述过滤单元以预定间隔单位对偏移向量执行过滤处理。

8、根据权利要求1的图像处理设备，其中，所述内插单元对由所述过滤单元进行过滤处理的偏移向量执行样条内插。

9、根据权利要求1的图像处理设备，其中，所述内插单元以预定间隔单位对由所述过滤单元进行过滤处理的偏移向量执行样条内插。

10、根据权利要求1的图像处理设备，还包括控制单元，适于控制要被重复预定次数的所述偏移向量计算单元、所述过滤单元、所述内插单元和所述配准单元的处理。

11、根据权利要求10的图像处理设备，其中，在计算关于加权平均处理的权重系数的情况下，所述过滤单元基于由所述控制单元重复处理的次数改变过滤处理的处理内容。

12、一种图像处理方法，其输出第一图像和第二图像之间的差图像，包括：

输入步骤，输入第一图像和第二图像；

偏移向量计算步骤，将多个感兴趣区域分别设置到第一图像和第二图像，并且对于每个感兴趣区域计算指示第一图像和第二图像之间的失配量的偏移向量；

过滤步骤，对偏移向量执行过滤处理；

内插步骤，对在所述过滤步骤中进行过滤处理的偏移向量执行内插；

配准步骤，根据在所述内插步骤中内插的偏移向量配准第一图像和第二图像；

相减操作步骤，在各配准的图像上的相应像素之间执行相减操作；以及

输出步骤，输出在所述相减操作步骤中获取的差图像。

Images