CN1637781A - 图像处理装置、图像处理方法以及图像处理*** - Google Patents

Abstract

本发明提供图像处理装置、图像处理方法以及图像处理***。输入第1、第2图像F1、F2(S110)。接着对第1、第2图像F1、F2进行对位(S120)。该对位检测在2张图像F1、F2之间存在的位置偏差，使用解剖学特征实施使一张的图像和另一张图像的对应像素位置一致的变形(Warping)处理。接着，使用对位完成图像X1、X2，进行独立成分分析处理(S130)。在主要成分判断中，在独立成分分析处理中的2个输出图像中，判断哪个是主要成分图像，哪个是经时变化成分图像(S140)。如果把一张判断为是主要成分图像，则另一张成为经时变化成分图像，然后输出主要成分图像和经时变化检测图像(S150)。

Description

图像处理装置、图像处理方法 以及图像处理***

技术领域

本发明涉及借助于CR和MRI等的医用诊断机器将不同时间拍摄的多张图像进行对位并进行差分处理的图像处理装置、图像处理方法以及图像处理***。

背景技术

近年，在医用图像诊断领域中的数字图像的利用得到发展。例如，使用半导体传感器拍摄X线图像的装置和以往的使用银盐照片的放射线照片***相比可以在极其宽的放射线曝露区域记录图像，进而具有在图像的保存和传送中容易构筑高效率***的实用的优点。

另外，通过数字化医用图像，有可能产生在以往的银盐照片中困难的诊断形态，即，以往在患者的经过观察等中在比较在不同的时间拍摄的X线图像时，一般是把胶片架在X线照片观察箱上进行比较读影。在该读影中，重要的是发现不同点，但为了查找2张胶片的不同点，需要具有能够非凡的洞察细部的观察力和丰富的经验。

在此如果使用数字图像，对不同的时间拍摄的2张数字图像进行对位使其与正常的解剖学构造一致，通过进行差分处理生成并输出差分图像，通过把该差分图像和成为基础的1对图像进行比较读影，可以更正确地把握图像间的变化。这是因为在生成的差分图像中只显示不同点，所以与通过胶片的读影相比可以容易地进行比较读影。这是支援被称为经时差分CAD(Computed Aided Diagnosis)的读影的技术。

有关生成这样的差分图像的处理方法被公开在美国专利5359513号中，通过将在不同的时间拍摄的2张胸部X线图像进行对位，可以生成差分图像。

图16是这样的差分图像的生成显示装置的构成图。第1以及第2图像是在不同的时间拍摄的特定部位的医用图像数据，但因为是在不同的时间拍摄的，所以两者的拍摄条件并不完全一致。例如在X线图像的情况下，X线的线量、管电压、管电流等的拍摄条件参数不一致，两者的图像浓度不同，另外摄影被拍摄体的位置和姿势也不同。

拍摄被拍摄体的2张图像数据在浓度补正单元1中被补正，使得两者的图像信号的浓度值分布大致相等，在对位单元2中使用解剖学构造的位置关系进行补正，在差分计算单元3中在对应的像素之间进行差分处理生成差分图像。而后，差分图像在显示单元4上和第1以及第2图像一同显示。

但是经时差分CAD在检测经时变化的医用图像处理方法中，并不是能完全满足的。差分处理的精度在浓度补正单元1中，被相等地补正两者的图像信号的浓度值分布的补正处理的精度左右。

作为例子，如果考虑从相同位置拍摄被拍摄体的2张图像，则假设第1图像在时刻t1摄影，第2图像在时刻t2摄影。如果为了简单把被拍摄体的姿势设置为一样，则第1图像和第2图像共用的部分和相似的部分非常多，而且细部和浓度必然不同。例如，如被拍摄体是患者的胸部，拍摄方法是X线摄影，则第1图像和第2图像横隔膜的位置等不同，另外，因为X线照射量不同所以拍摄图像的浓度不同。即使不是X线图像的情况下，当被拍摄体是人的情况下，可以认为不能完全同样地拍摄第1图像和第2图像。

如果把第1图像的浓度设置为x1(x，y)，把第2图像的浓度设置为x2(x，y)，把经时变化部分的浓度设置为Δ(x，y)，则它们有以下关系。

x1(x，y)＝s(x，y)

x2(x，y)＝a·s(x，y)+Δ(x，y)

s(x，y)是第1图像和第2图像共同的构造物的浓度，a是常数。因为第1、第2图像各自的摄影条件不同，所以常数并不限于1。在浓度补正单元1中推测常数a使第2图像的浓度与第1图像的浓度相等，如果把该推测结果设为a’，则差分处理结果h(x，y)为下式。

h(x，y)＝x2(x，y)-a’·x1(x，y)

       ＝a·s(x，y)+Δ(x，y)-a’·s(x，y)

       ＝(a-a’)·s(x，y)+Δ(x，y)

在正确地推测常数a的情况下，因为a’＝a，所以

h(x，y)＝Δ(x，y)

如果常数a的推测结果不正确，则

|(a-a’)·s(x，y)|＞|Δ(x，y)|

不能检测经时变化Δ(x，y)，经时变化Δ(x，y)还被埋入在左边的(a-a’)·s(x，y)中。另外，当在细部构造中没有不同的情况下，经时变化Δ(x，y)＝0，差分处理结果h(x，y)应该输出0。但是，常数a的推测结果如果不正确则为下式。

h(x，y)＝(a-a’)·s(x，y)

即，即使在被拍摄体中没有病变等的变化的情况下，医生也判断为有经时变化。该判断诱发医生的诊断错误。

本发明就是为了解决上述那样的问题而完成的，其目的在于提供一种可以精确地检测出在时间系列上不同的2张图像间的经时变化的图像处理装置、图像处理方法以及图像处理***。

为了实现上述目的的本发明一个方式的图像处理方法包含：输入2张图像的图像输入步骤；对上述2张图像进行对位的对位步骤；使用在该位置对位步骤中所对位的图像进行独立成分分析处理的独立成分分析处理步骤。

另外，为了实现上述目的的本发明的一个方式的图像处理装置具备以下构成：即，输入2张图像的图像输入单元；对上述2张图像进行位置对位的对位单元；使用在该对位装置中所对位的图像进行独立成分分析处理的独立成分分析处理单元。

本发明的其它特征以及优点，根据结合附图的以下的说明将会明了，其中，在附图中对同样或者相似的部分附加相同的参考标号。

附图说明

附图包含在说明书中，构成说明书的一部分，表示本发明的一个实施方式，并与说明书所记述的内容一起用于说明本发明的原理。

图1是图像处理***的构成图。

图2是独立成分分析处理单元的构成图。

图3是实施例中的经时变化检测单元的构成图。

图4是实施例中的经时变化显示单元的构成图

图5是全体处理的流程图。

图6是全体处理的流程图。

图7是涉及对位的子程序的流程图。

图8是涉及独立成分分析处理的子程序的流程图。

图9是与主要成分判断处理有关使用相关值的子程序的流程图。

图10是与主要成分判断处理有关使用直方图的子程序的流程图。

图11是1组ROI设定的说明图。

图12是与移位向量有关的说明图。

图13是独立成分分析处理的说明图。

图14A、图14B是独立成分分析处理的一部分的说明图。

图15是经时变化强调处理的说明图。

图16是以往例子的说明图。

具体实施方式

下面将根据附图详细说明本发明的优选实施方式。

[实施例1]

图1是本发明的图像处理装置的构成图，从图像生成单元10输出的2张图像F1、F2经由经时变化检测单元20发送到经时变化显示单元30。在图像生成单元10中，在具有传感器111的摄像单元11中生成的图像输出经由存储单元12、选择单元13连接到经时变化检测单元20。

在经时变化检测单元20中，设置接收图像生成单元102的2张图像F1、F2的对位单元21，同时2张图像F1、F2直接被发送到经时变化显示单元30。在对位单元21中生成的对位图像X1、X2在被送到独立成分分析处理单元22的同时，发送到经时变化显示单元30。在独立成分分析处理单元22中得到的独立成分图像Y1、Y2、独立成分系数矩阵A被发送到主要成分判断单元23，在主要成分判断单元23中生成的主要成分图像、经时变化成分图像以及独立成分系数矩阵A被发送到经时变化显示单元30。

在经时变化显示单元30中，经时变化检测单元20的输出经由经时变化强调图像生成单元31、显示图像选择单元32与显示单元33连接。

在图像生成单元10中，借助于传感器111所得到的由摄像单元11生成的被拍摄体的拍摄图像被存储在存储单元12中，由选择单元13选择2张图像并输出数字摄影图像F1、F2。

作为传感器111的摄影装置只要输出数字图像信号即可，例如使用FPD(Flat Panel Detector)装置和CR(Computed Radiography)装置和CT(Computed Tomography)装置等的X线摄像装置，MRI(Magnetic Resonance Imaging)装置和US(Ultra Sonic)回波装置等。

存储单元12和选择单元13例如可以使用控制上述的控制摄像装置的计算机附带的硬盘等的外部存储装置，鼠标、键盘等的输入指示装置来实现。即，用户使用这些输入指示装置可以通过对话指定被保存在硬盘上的图像，也可以根据规定的规则自动地选择2张图像。

经时变化检测单元20输入来自图像生成单元10所输出的摄影图像F1、F2，输出后述的主要成分图像以及经时变化成分图像。在经时变化检测单元20中，已输入的2张摄影图像在对位单元21中进行用于补正相互位置偏移的对位，在后续的独立成分分析处理单元22中进行分析，输出构成图像的各成分，进而在主要成分判定单元23中生成并输出经时变化检测图像以及经时变化强调图像。

图2是经时变化检测单元20的独立成分分析处理单元22的构成图。该独立成分分析处理单元22由一维化处理单元221、分析前处理单元222、分离矩阵计算单元223、分析后处理单元224、二维化处理单元225、前处理信息存储单元226、二维信息存储单元227构成。

对位单元21、独立成分分析处理单元22、主要成分判断单元23可以作为软件实现，各单元可以作为构成在计算机上动作的软件的模块来实现。或者，作为其它的实现方法，也可以借助于硬件来安装对位单元21、独立成分分析处理单元22、主要成分判断单元23的一部分或者全部。

在经时变化显示单元30中，经时变化显示单元30进行输入的各图像的处理并输出，以适宜医生读影。经时变化强调图像生成单元31使用经时变化检测单元20生成的图像以及图像生成单元10生成的图像来生成经时变化强调图像。显示图像选择单元32选择经时变化强调图像生成单元31生成的图像，输入到后续的显示单元33。显示单元33显示显示图像选择单元32输出的图像。

为了实施这些功能，例如可以使用图3所示那样形态的构成。在图3中，计算机41经由网络42与1个或者多个图像生成装置43连接。作为该图像生成装置43，相当于上述图像生成单元10的各种摄影装置。或者，借助于未图示的摄影装置拍摄的图像被预先保持在文件服务器44中，即使把它和图像生成装置10同等处理，也可以完全同样使用。

在图3中，经时变化检测单元20以及经时变化显示单元30可以用计算机41以及与之连接的***机器实现，作为一个例子经时变化检测单元20以及经时变化显示单元30可以作为显示装置45以及存储在硬盘46中的程序实现。这种情况下，通过未图示的用户输入，被保存在硬盘46中的程序经由总线47读出到RAM48中，CPU49顺序执行RAM48内的程序实现上述功能，处理的结果，把得到的图像显示在显示装置45上。在此作为显示装置45例如可以使用CRT监视器，液晶显示器等。

另外，实施本发明的程序并不必须限定被保存在硬盘46中，也可以存储在作为外部存储装置与计算机41连接的光磁盘50，或者经由网络连接的文件服务器44中。或者，本发明的经时变化检测单元20以及经时变化显示单元30也可以把一部分或者全部作为安装在计算机41中的专用硬件51来实现。

进而，本发明并不限定于与经由如图3所示的网络42连接，也可以采用图4所示的构成。在图4的构成中，计算机41和图像生成装置43直接连接，直接输入图像。另外，如上所述当图像生成装置43具备计算机41的情况下，在其中也可以通过上述的形态在摄影装置中实现。

图5、图6是本实施例中的全体处理流程图，图7是与全体处理流程图的对位有关的子程序的流程图，图8是与全体处理流程图的独立成分分析处理有关的子程序的流程图。图9是与全体处理流程图的主要成分判断处理有关使用相关值的子程序的流程图。图10是与全体处理流程图的主要成分判断处理有关使用直方图的子程序的流程图。

以下，参照图5～图10的各流程图说明各部分功能的动作。首先，输入第1图像F1和第2图像F2(步骤S110)。该图像在本实施例的情况下在图像生成单元10中由摄像单元11生成。由摄像单元11生成的图像被储存在储存用的存储单元12中，在患者单位、部位单位中按照时间系列存储。从被储存在存储单元12中的特定患者的特定部位的图像中，选择单元13选择2张图像F1、F2。这样，生成在步骤S110中的第1图像F1和第2图像F2。

接着，对第1图像F1和第2图像F2进行对位(步骤S120)。该对位对于被输入到对位单元21中的2张图像F1、F2，检测在2张图像F1、F2之间存在的位置偏移，使用解剖学特征检测2张图像间F1、F2的对应像素位置，实施变形(warping)处理，使得一张图像与另一张图像对齐。该对位的动作首先大致对第1图像F1和第2图像F2进行匹配(步骤S121)。其后，对第1图像F1和第2图像F2进行局部匹配(步骤S122)。通过将该匹配分为2个阶段来提高匹配精度。

在匹配时设定成为对位时的记号的关注区域ROI(Region ofinterest)。图11是展示已被设定的1组ROI的图，对第1图像F1设定样板ROI-TR，对第2图像F2设定检索ROI-SR。也可以对第2图像F2设定样板ROI-TR，对第1图像F1设定检索ROI-SR。有关用于对位的ROI的设定方法以及对位，在特开平7-37074号公报中详细叙述。

匹配通过把样板ROI-TR以及检索ROI-SR作为1组，对输入图像设定多个ROI的组，进行样板匹配。该样板匹配如图11所示，一边取得已设定的样板ROI-TR内的图像数据和检索ROI-SR内的图像数据的匹配一边使位置变化，求解一致度最高的位置，把和样板ROI-TR的设对位置的差作为2张图像F1、F2的局部位置偏差量进行输出。

在图11的情况下，因为第2图像F2与第1图像F1比较，向右下平行移动，所以偏移向量如图12所示。偏移向量意味着使图像上的某一点的坐标(x，y)移位到新的坐标(x’，y’)，所以这样称谓。

在步骤S124中，在已输入的多个移位向量中适用规定的插补模型，输出插补模型的参数。所谓插补模型的参数是用于唯一确定插补模型的常数。例如，如果把关于y＝ax²+bx+c的x的2次式作为插补模型，则(a，b，c)这3个常数是插补模型的参数。插补模型可以使用这样的二次多项式，也可以使用样条函数(spline)插补等其它的插补模型。插补模型的次数可以预先基于经时变化检测图像的画质评价等确定适宜的值，也可以进行适宜地确定使插补误差比规定的值小。

使用被计算出的插补模型的参数，进行把图像F1、F2中的一张与另一张对齐的处理(步骤S125)。因为该处理进行图像变形，所以被称为变形(Warping)处理。变形处理可以采用线性变换的图像变形，也可以采用非线性变换的图像变形，当在本实施例中进行仿射(affine)变换。通过该变形(Warping)处理，输出对位完成图像X1、X2。这样，对位处理(步骤S120)结束，但采用ROI设定和移位向量插补处理的对位计算也可以进行多次以提高精度。

以下，使用对位完成图像X1、X2，进行独立成分分析处理(步骤S130)。借助于独立成分分析处理单元22，用独立成分分析(ICA，Independent Component Analysis)对从对位单元21输出的对位完成图像X1、X2进行处理，生成主要成分图像和经时变化成分图像的2张图像。在此，所谓主要成分是与2张图像X1、X2共同的构造物，所谓经时变化成分是2张图像X1、X2之间的细部的不同。独立成分分析在主要成分和经时变化成分在统计上是独立的信号时，从图像X1、X2中分别抽出主要成分和经时变化成分。

当2张图像F1、F2是医用图像的情况下，例如主要成分图像相当于正常解剖学构造，经时变化成分图像相当于非正常的解剖学构造。在不同时间拍摄同一患者的同一部位，当在该拍摄间隔期间内产生非正常的解剖学变化的情况下，在经时变化成分图像中出现该变化。

例如，图像F1、F2是胸部图像，假设第1图像F是发生癌症前的拍摄图像，第2图像F2是发生癌症后的拍摄图像。如果假设是在对位单元21中经对位的完成图像X1、X2，则主要成分图像是图像X1和图像X2共用的胸部构造的图像，经时变化成分图像是癌症的图像。因为胸部构造和癌是独立的事项，所以通过独立成分分析可以抽出2个事项。

主要成分图像s(x，y)和经时变化成分图像Δ(x，y)如果使用经对位的第1图像X1＝x1(x，y)和经对位的第2图像X2＝x2(x，y)，则存在以下关系。

x1(x，y)＝a₁₁·s(x，y)+a₁₂·Δ(x，y)

x2(x，y)＝a₂₁·s(x，y)+a₂₂·Δ(x，y)  ......(1)

通过对独立成分分析，输出s(x，y)，Δ(x，y)，系数a₁₁，a₁₂，a₂₁、a₂₂，该输出结果在后述的主要成分判断单元23中使用。在以往的经时差分CAD中的差分处理中，在浓度补正单元1中，需要正确推测系数a₁₁和a₂₁，但在独立成分分析中因为自动地计算这些值，所以浓度补正处理是不需要的。

在此，为了以后的说明，把独立成分系数矩阵(行列)A定义为下式(2)([数1])。

[数1]

$A=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right)$

对象图像是胸部图像，第1对位完成图像X1是癌症发生前的图像，如果第2完成图像X2是癌症发生后的图像，则在式(1)中，可以认为a₁₂＝0，成为以下的式(3)。

x1(x，y)＝a₁₁·s(x，y)

x2(x，y)＝a₂₁·s(x，y)+a₂₂·Δ(x，y)  ......(3)

通过独立成分分析处理，输出作为胸部构造的主要成分图像s(x，y)，和相当于癌的经时变化成分图像Δ(x，y)。不需要对象图像是胸部图像，可以是另一部位的图像。另外，经时变化成分图像是表示经时变化的图像，并不限于癌，把被检测体的异常区域作为经时变化成分图像输出。图13是展示表示在独立成分分析处理单元22中的输入图像和输出图像关系的说明图。

独立成分分析处理单元22输出主要成分图像s(x，y)和经时变化成分图像Δ(x，y)的2张图像，而不知道该图像中的哪一张是主要成分图像，哪一张是经时变化成分图像。因此，在后续的主要成分判断部分23中，判断哪一张是主要成分图像s(x，y)。

所谓独立成分分析处理是从相互独立的信号混合存在的观测信号中，在没有信号独立性以外的预见性知识的情况下抽出原信号的数据分离方法，是在声音和图像处理的领域中受到关注的新的技术。不管原信号和结合系数是否都未知，根据重复学习推测结合系数，可以从观测信号中分离并抽出原信号。

即，在把N个独立的原信号设置为S(t)＝[s₁(t)，s₂(t)，...，s_N(t)]^T(t＝0，1，2，......)，把N个观测信号记述为X(t)＝[x₁(t)，x₂(t)，...，x_N(t)]^T(t＝0，1，2，......)时，在观测信号可以记述为X(t)＝A·S(t)的情况下，不使用原信号的独立性和观测信号X(t)以外的任何信息，推测原信号S(t)。在此，A是N·N的实数矩阵，[  ]^T是矩阵转置。

如果找到矩阵A的逆矩阵W，则原信号通过S(t)＝WX(t)这一计算可以复原。在独立成分分析中，假设各原信号s₁(t)是相互在统计上独立的信号，来计算矩阵W。即，进行Y(t)＝WX(t)这一计算，计算W使得得到的各信号Y(t)＝[y₁(t)，y₂(t)，...，y_N(t)]^T(t＝1，2，......)成为相互在统计上独立的信号。矩阵W因为从作为混合原信号S(t)的信号的X(t)中分离原信号，所以被称为分离矩阵。

求分离矩阵W的方法有一边逐次学习一边求的方法和一并求解的方法。在此，说明逐次学习的方法。另外，测定是否在统计上独立的尺度有相互信息量和非高斯性等，而在此采用相互信息量。相互信息量具有非负的值，只在相互统计上独立的情况下，是具有0这一性质的尺度。因而，查找信号Y(t)的各成分间的相互信息量尽可能接近0那样的矩阵W。

学习方法用以下步骤进行。

步骤(1)：进行分离矩阵W的初始化。

步骤(2)：计算信号Y(t)＝WX(t)。

步骤(3)：根据更新式W←W+ηΔW更新分离矩阵W。

η是学习系数，ΔW是更新的变化部分，如果把相互信息量设置为I(Y)，则有下式。

(ΔW)_ij∝-{(I(Y))/W_ij}W^TW

这是与相互信息量I(Y)有关的最急下降方向。

步骤(4)：设ε是收敛条件，在满足(ΔW)_ij＜ε为止，重复步骤(2)～(3)。

在此，补足相互信息量。如果预先把信号Y(t)的同时分布概率密度函数设置为p(Y)＝p(y₁，y₂，...，y_N)，则相互信息量I(Y)用下式定义。

[数2]

$I\left(Y\right)=\∫p\left(Y\right)\log \frac{p\left(Y\right)}{{\mathrm{\Π}}_{i=1}^{N}p\left(y_i\right)}\mathrm{dY}$

如上所述，I(Y)＝0，当y₁，y₂，...，y_N相互独立的情况下，即当

[数3]

$p\left(Y\right)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}p\left(y_i\right)$

时，为I(Y)＝0。因而，通过最小化相互信息量I(Y)，可以使Y的各成分独立。在步骤(3)中，用最急下降法求相互信息量I(Y)的最小值。

以下，如果考虑用分离矩阵W表现相互信息量I(Y)，则如步骤(3)所示，在相互信息量I(Y)减小的方向上更新分离矩阵W时需要。Y的熵(entropy)H(Y)如下：

H(Y)＝-∫p(Y)log p(Y)dY

因为根据Y＝WX，有下式

p(Y)＝p(X)/det(W)

dY＝det(W)dX

所以Y的熵H(Y)为下式。

H(Y)＝H(X)+log|det(W)|

相互信息量I(Y)利用熵H(Y)可得：

[数4]

$I\left(Y\right)=-H\left(Y\right)+\underset{j=1}{\overset{N}{\mathrm{\Σ}}}H\left(y_j\right)$

所以结果如下：

[数5]

$I\left(W\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}H\left(y_j\right)-H\left(X\right)-\log \left|\det \left(W\right)\right|$

可以用分离矩阵W的函数表现相互信息量I(Y)。

以下，具体地计算在步骤(3)中的更新变化部分ΔW。

[数6]

${\left(\mathrm{\ΔW}\right)}_{\mathrm{ik}}\≅-\frac{\∂I\left(W\right)}{{\∂W}_{\mathrm{ik}}}{W}^{T}W$

$=-\frac{\∂({\mathrm{\Σ}}_{j=1}^{N}H\left(y_j\right)-H\left(X\right)-\log |\det \left(W\right)\left|\right)}{{\∂W}_{\mathrm{ik}}}{W}^{T}W$

$=-\frac{\∂\left({\mathrm{\Σ}}_{j=1}^{N}h\left(y_j\right)\right)-\∂\left(\log \right|\det \left(W\right)\left|\right)}{{\∂W}_{\mathrm{ik}}}{W}^{T}W$

在此，最后的式子右边的第1项如果参照“S.Haykin，NeuralNetworks，Prentice Hall，1999.”，则可以用下式来近似。

[数7]

在此(y_i)被称为活性化函数，是非线形函数，例如是如下的函数。

(y_i)＝y_i ²

另外，也可以是以下函数。

(y_i)＝(1/2)y_i ⁵+(2/7)y_i ⁷+(15/2)y_i ⁹+(2/15)y_i ¹¹

进而，det(W)在把W的j行k列的要素设置为W_ik时，使用余因子W_jk，可以如下进行余因子分解。

[数8]

$\det \left(W\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{jk}}{\hat{W}}_{\mathrm{jk}}$

$j=\mathrm{1,2},...,N$

因而，基于右边第2项的对数的W_jk的微分如果归纳为：

[数9]

$\frac{\∂}{{\∂W}_{\mathrm{ik}}}\log \left|\det \left(W\right)\right|$

$=\frac{1}{\det \left(W\right)}\frac{\∂}{{\∂W}_{\mathrm{ik}}}\det \left(W\right)$

$=\frac{{\hat{W}}_{\mathrm{ik}}}{\det \left(W\right)}={\left(W_{\mathrm{ik}}^T\right)}^{-1}=W_{\mathrm{ik}}^{-T}$

则ΔW如下式：

ΔW_jk＝{(W^-T)_jk-(y_i)x_k}W^TW

如果把它用矩阵表现，则如下式：

ΔW＝{(W^-T)-(Y)X^T}W^TW

在此，如果使用Y^T＝X^TW^T这一关系，则为：

ΔW＝{(W^-T)-(Y)X^T}W^TW

   ＝{I-(Y)X^TW^T}W

   ＝{(I)-(Y)Y^T}W

得到在步骤(3)中的具体的更新变化部分ΔW。另外，上式的I为单位矩阵。

这样，对位完成图像X1、X2在一维化处理单元221中，从二维数据形式向一维数据形式变换。变换后的完成图像X1、X2如下式表示。另外，n是图像的像素数。

X(n)＝[x₁(n)，x₂(n)]^T    (n＝1，2，......，N)

x₁(n)和x₂(n)的构成方法如图14A所示可以是从图像的左角顺序取出排列的方法，如图14B所示也可以是Z形排列的方法，但x₁(n)和x₂(n)的构成方法必须相同。在变换中使用的排列方法被保持在二维信息存储单元227中。该信息在分析处理后从一维数据形式向二维数据形式变换时被参照。

在变换为一维数据形式后，X(n)为了在分析前处理单元222中高精度地进行独立成分分析处理而被整形。在整形中所使用的信息被保持在前处理信息存储单元226中，在进行分析处理后的整形逆变换时使用。在本实施例中，分别减去直流成分值使x₁(n)的平均值和x₂(n)的平均值等于0。被减去的直流成分值信息保持在前处理信息存储单元226中。

经整形的X(n)被输入到成为独立成分分析处理中心的分离矩阵计算单元223。分离矩阵计算单元223计算分离矩阵W，输出独立成分系数矩阵A和独立成分Y(n)。

该分离矩阵W用前面叙述过的学习步骤(1)～(4)计算。

步骤(1)：进行分离矩阵W的初始化(步骤S131)。在本实施方式中，如以下那样进行初始化。

[数10]

$W=\left(\begin{array}{cc}1& 0\\ 1& 1\end{array}\right)$

步骤(2)：计算独立成分Y(t)＝WX(t)(步骤S132)。

步骤(3)：按照更新式更新分离矩阵(步骤S133)。更新式如上所述为：

W←W+ηΔW

ΔW＝{I+(Y)Y^T}W

(y_i)＝(1/2)y_i ⁵+(2/7)y_i ⁷+(15/2)y_i ⁹+(2/15)y_i ¹¹

η＝1

步骤(4)：直到满足(ΔW)_ij＜ε为止，重复步骤(2)～(3)(步骤S134)。关于收敛条件ε在本实施例中假设为0.0001。另外，在本实施方式中把学习系数η设定为常数，但也可以是使η和更新次数一同减小那样的更新方法。

步骤(5)：得到独立成分Y(n)＝[y₁(n)，y₂(n)]^T(n＝1，2，...，N)(步骤S135)。另外，通过计算A＝W^-1，得到以下的独立成分系数矩阵A。

[数11]

$A=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right)$

由分离矩阵计算单元223输出的独立成分Y(n)使用被保持在前处理信息存储单元226中的信息，在分析后处理单元224中进行整形逆变换。在本实施方式的情况下，通过把被保持在前处理信息存储单元226中的直流成分值相加到独立成分Y(n)来进行整形逆变换。

在分析后处理单元224中进行了整形逆变换的独立成分Y(n)使用被保持在二维信息存储单元227中的信息，在二维化处理单元225中将一维数据形式变换为二维数据形式，作为独立成分图像Y1、Y2进行输出。

在此，如果说道和以往的经时差分CAD的构成的不同时，则在以往的经时差分CAD中，把第1图像和第2图像作为输入，作为差分处理的结果输出1张图像。在本实施例中，把第1图像和第2图像作为输入，作为独立成分分析处理的结果输出2张图像。另外，作为被检测体的异常区域的经时变化成分图像Δ(x，y)虽然是和以往例的差分处理的结果相似的图像，但并不一样。

在主要成分判定单元23中，对独立成分分析处理单元22中的2个输出图像，判断哪个是主要成分图像，哪个是经时变化成分图像(步骤S140)。

在本实施方式中，使用相互相关值判断主要成分图像和经时变化图像的判断。通过在独立成分分析处理单元4中的处理结果，得到以下的关系式。

x1(x，y)＝a₁₁·y1(x，y)+a₁₂·y2(x，y)

x2(x，y)＝a₂₁·y1(x，y)+a₂₂·y2(x，y)  ......(4)

判断y1(x，y)和y2(x，y)的哪个是主要成分图像。经时变化成分图像与主要成分图像相比因为信号值非常小，所以在式(4)的第1式中，比较左边x1(x，y)和第1项a₁₁·y1(x，y)的相关值R₁₁的绝对值，和左边x1(x，y)和第2项a₁₂·y2(x，y)的相关值R₁₂的绝对值，把大的一方判断为主要成分图像。即，如果

|R₁₁|＞|R₁₂|

则判断为y1(x，y)是主要成分图像，如果

|R₁₁|＜|R₁₂|

则判断为y2(x，y)是主要成分图像(步骤S141a，S142a)。当相关的绝对值|R₁₁|和|R₁₂|大致相等的情况下，判断为独立成分分析处理单元22不能正常地分离主要成分图像和经时变化成分图像。

作为相关值的计算方法，虽然在本实施例中使用了相互相关值，但并不限定于此，也可以使用其它的相关值计算方法。在y1(x，y)和y2(x，y)中，如果一个被判断为主要成分图像，则另一个就为经时变化成分图像。这样，主要成分判断单元23输出主要成分图像和经时变化检测图像(步骤S150)。

经时变化强调图像生成单元31使用经时变化检测单元20生成的图像，生成经时变化检测图像和经时变化强调图像(步骤S160)。经时变化检测图像h1(x，y)、h2(x，y)把经时变化成分图像和其对应的系数相乘而生成。如果假设y2(x，y)是经时变化成分图像，则生成以下式子。

h1(x，y)＝a₁₂·y2(x，y)

h2(x，y)＝a₂₂·y2(x，y)

进而，如以下那样生成经时变化强调图像w1(x，y)，w2(x，y)，k1(x，y)，k2(x，y)，v1(x，y)，v2(x，y)。其中，c是充分大的常数。

w1(x，y)＝c·h1(x，y)

w2(x，y)＝c·h2(x，y)

k1(x，y)＝y1(x，y)+c·h1(x，y)

k2(x，y)＝y2(x，y)+c·h2(x，y)

v1(x，y)＝x1(x，y)+c·h1(x，y)

v2(x，y)＝x2(x，y)+c·h2(x，y)

经时变化强调图像w1(x，y)，w2(x，y)是强调经时变化检测图像的图像。因为经时变化部分是微小的变化，所以即使不知道有没有经时变化的情况下，如果看经时变化强调图像，则可以一目了然地判断是否有经时变化。

经时变化强调图像k1(x，y)、k2(x，y)是同样地强调经时变化，是重叠写在主要成分图像上的图像。因为在主要成分图像上重叠显示经时变化部分，所以可以判断经时变化部分位于主要成分部分的哪个位置。

经时变化强调图像v1(x，y)、v2(x，y)因为同样在被拍摄体的拍摄图像上重叠显示经时变化部分，所以可以判断经时变化部分处于被拍摄体的哪个位置上。

图15所示的说明图展示经时变化强调图像与对位完成图像和经时变化检测图像和主要成分图像有什么样关系。

在图像显示选择单元32中，诊断者选择经时变化强调图像生成单元31输出的图像。即，选择多个或者只选择1个主要成分图像、经时变化成分图像、经时变化检测图像、经时变化强调图像、对位完成图像。在显示单元33中，同时显示图像显示选择单元32选择的图像。

[实施例2]

作为实施例2说明主要成分判断单元采用直方图判别的方式。对y1(x，y)和y2(x，y)的两者计算直方图，把直方图的分布幅度小的一方判断为经时变化成分图像(步骤S141b，S142b)。这是因为虽然主要成分图像的像素值从最小值到最大值分布全部的值，但经时变化成分图像只分布在某一窄范围内的缘故。在y1(x，y)和y2(x，y)中，如果判断为一方是经时变化成分图像，则另一方为主要成分图像。这样，主要成分判断单元23输出主要成分图像和经时变化检测图像(步骤S150)。

[实施例3]

作为实施例3说明采用高通滤波器等进行强调图像生成的方式。

在本实施例3中，如以下那样生成经时变化强调图像w1(x，y)，w2(x，y)，k1(x，y)，k2(x，y)，v1(x，y)，v2(x，y)。

w1(x，y)＝u(x，y)*h1(x，y)

w2(x，y)＝u(x，y)*h2(x，y)

k1(x，y)＝y1(x，y)+u(x，y)*h1(x，y)

k2(x，y)＝y2(x，y)+u(x，y)*h2(x，y)

v1(x，y)＝x1(x，y)+u(x，y)*h1(x，y)

v2(x，y)＝x2(x，y)+u(x，y)*h2(x，y)

在此，u(x，y)是强调高频成分的高通滤波器函数，运算符*表示卷积计算。

经时变化强调图像w1(x，y)，w2(x，y)，k1(x，y)，k2(x，y)，v1(x，y)，v2(x，y)强调经时变化，可以判断经时变化部分在哪个位置上。

在实施方式1、2中，通过把经时变化检测图像放大常数倍生成经时变化强调图像，但在本实施例3中使用强调高频的高通滤波器来生成。生成方法并不限于这些，还可以使用非线性滤波器等各种滤波器。

在图像显示选择单元32中，选择经时变化强调图像生成单元31输出的图像。即，选择多个或者只选择1个主要成分图像、经时变化成分图像、经时变化检测图像、经时变化强调图像、对位完成图像。虽然在实施例1、2中诊断者选择了所希望的图像，但在本实施例中，图像显示选择单元32用规定的规则选择图像。所谓该规定的规则例如是“初诊的患者因为需要详细信息所以选择上述全部图像”，“经过观察中的患者因为掌握关注区域所以只选择经时变化检测图像”等。

[实施例4]

作为实施例4说明显示形态的图像自动选择、顺序显示等的其它方式。在本实施例4中，作为测量是否是统计上独立的尺度使用相互信息量，但也可以使用4次积累量和更高次的统计量。另外，作为独立成分分析计算方法，还可以使用Fast ICA和JADE。这些算法在以前的甘利俊一、村田升合著，独立成分分析、临时另册·数理科学“SGC图书18”，科学研究社(2002)中有解说。

在本实施例中，虽然求分离矩阵W的方法使用一边逐次学习一边求解的方法，但也可以是一并求的方法。一边逐次学习一边求的方法在解到收敛为止需要时间，但有用少量的运算用计算区域就可完成的优点。另一方面，一并求的方法虽然求解的时间不长，但存在需要大量的运算用计算区域的问题。采用哪种方法的判断依赖于图像尺寸的大小和计算机的能力。

分离矩阵W由记载于实施例1中的学习步骤(1)～(4)来计算。在步骤(5)中，得到[数11]的独立成分系数矩阵A。

在分析前处理单元222中减去直流成分值后，也可以为了使X(n)无相关化而实施白色化处理。如果实施白色化处理，则分离矩阵W具有以少的步骤数收敛的趋势。但是，这仅限于在图像F1、F2的S/N比大的情况。

用分离矩阵计算单元223输出的独立成分Y(n)使用被保持在前处理信息存储单元226中的信息，用分析后处理单元224进行整形逆变换。在本实施方式4的情况下，通过把被保持在前处理信息存储单元226中的直流成分值相加到独立成分Y(n)，进行整形逆变换。

在分析后处理单元224中进行了整形逆变换的独立成分Y(n)使用被保持在二维信息存储单元227中的信息，在二维化处理单元225中从一维数据形式变换为二维数据形式，作为独立成分图像Y1、Y2输出。

和实施例1一样经由步骤S140、步骤S141a、S142a，在步骤S150中主要成分判断部分23输出主要成分图像和经时变化检测图像。

经时变化强调图像生成单元31在步骤S160中，和实施例1一样进行卷积运算计算生成下式。

h1(x，y)＝a₁₂·y2(x，y)

h2(x，y)＝a₂₂·y2(x，y)

进而，和实施例3一样如以下那样生成经时变化强调图像w1(x，y)，w2(x，y)，k1(x，y)，k2(x，y)，v1(x，y)，v2(x，y)。

w1(x，y)＝u(x，y)*h1(x，y)

w2(x，y)＝u(x，y)*h2(x，y)

k1(x，y)＝y1(x，y)+u(x，y)*h1(x，y)

k2(x，y)＝y2(x，y)+u(x，y)*h2(x，y)

v1(x，y)＝x1(x，y)+u(x，y)*h1(x，y)

v2(x，y)＝x2(x，y)+u(x，y)*h2(x，y)

在实施例1、2中，是通过把经时变化检测图像放大常数倍生成经时变化强调图像，但在本实施例4中和实施例3一样，使用强调高频的高通滤波器。

在显示单元33中，显示图像显示选择单元32选择的图像。在实施例1、2、3中，虽然同时显示图像显示选择单元32选择的图像，但在本实施方式4中，每次1张顺序地显示。显示方法并不限定于此，也可以使用其它的方法，例如眼镜方式的显示装置，显示右眼和左眼不同的图像。

如果采用本发明的图像处理方法以及处理装置，则通过独立成分分析处理，可以高精度地进行在检测经时变化时重要的图像浓度补正，可以精确地生成检测经时变化的图像。

检测出经时变化的图像因为表示肿瘤、炎症等病状的变化，所以医生和技师通过观看精确地生成的检测图像，可以容易确认经时变化。如果诊断行为长时间地在诊断中消耗使集中力下降，则漏看病状变化的可能性高，但也有降低该可能性的效果，抑制诊断结果的差异程度，通过这些相乘效果提高诊断精度。

进而，本发明的目的当然也可以向装置或者***提供存储有实现实施方式的装置或者***功能的软件程序代码的存储介质，通过该装置或者***的计算机(CPU或者MPU等)读出被存储在存储介质上的程序代码执行来实现。

这种情况下，从存储介质读出的程序代码自身实现实施方式的功能，存储有该程序代码的存储介质以及该程序代码构成本发明。

作为用于提供程序代码的存储介质，可以使用ROM、软盘(注册商标)、硬盘、光盘、光磁盘、CD-ROM、CD-R、磁带、非易失性存储卡等。

另外，计算机通过执行读出的程序代码，不仅可以实现实施方式的功能，而且根据该程序代码的指示，在计算机上运行的OS等进行实际处理的一部分或者全部，通过该处理实现本实施方式1～4的功能，这当然也包含在本发明的实施方式中。

进而，从存储介质读出的程序代码在写入在被***在计算机上的功能扩展板和在与计算机连接的功能扩展单元中具备的存储器中后，根据该程序代码的指示，在该功能扩展板和在功能扩展单元中具备的CPU等进行实际处理的一部分或者全部，通过该处理实现实施方式的功能，这当然也包含在本发明的实施方式中。

当在这样的程序或者存储有该程序的存储介质中适用本发明的情况下，该程序例如由与上述的图5～图10所示的流程图对应的程序代码构成。

如上所述，如果采用本发明，则可以提供精确地检测在时间系列上不同的2张图像间的经时变化的图像处理装置、图像处理方法以及图像处理***。

Claims (10)

1、一种图像处理方法，其特征在于包括：

输入2张图像的图像输入步骤；

对上述2张图像进行对位的对位步骤；

使用在该位置对位步骤中所对位的图像进行独立成分分析处理的独立成分分析处理步骤。

2、根据权利要求1所述的方法，其特征在于还包括：

进行用于把在上述独立成分分析步骤中输出的处理结果区分为主要成分图像和经时变化成分图像的判断的主要成分判断步骤。

3、根据权利要求1所述的方法，其特征在于：

上述独立成分分析处理步骤具备：

把在上述对位步骤中进行了对位的图像从二维数据形式变换为一维数据形式的一维处理步骤；

从在该一维处理步骤中得到的一维数据中减去直流成分的分析前处理步骤；

从在该分析前处理步骤中得到的一维数据中计算分离矩阵、独立成分系数矩阵以及独立成分的分离矩阵计算步骤；

在该分离矩阵计算步骤中计算出的独立成分上加上上述直流成分的分析后处理步骤；

把在分析后处理步骤中加上了直流成分的独立成分变换为二维数据形式的二维处理步骤。

4、根据权利要求2所述的方法，其特征在于：

上述主要成分判断步骤调查从上述独立成分分析处理中所输出的处理结果与在上述独立成分分析处理步骤中使用的经上述对位的图像和另一张图像的2张图像中的一张的相似性，判断上述主要成分图像。

5、根据权利要求2所述的方法，其特征在于：

上述主要成分判断步骤将从上述独立成分分析处理输出的处理结果进行直方图处理，把该直方图的分布宽度小的一方判断为上述经时变化成分图像。

6、根据权利要求2所述的方法，其特征在于还包括：

通过使借助于滤波处理强调了在上述主要成分判断步骤中判断的上述经时变化成分图像的图像，与在上述图像输入步骤中所输入的图像、或者在上述独立成分分析处理步骤中使用的经上述对位的图像、或者在上述主要成分判断步骤中判断的上述主要成分图像重合，生成强调了经时变化部分的图像的经时变化强调图像生成步骤。

7、根据权利要求1所述的方法，其特征在于：

上述2张图像是在不同的时刻拍摄的同一部位的图像。

8、根据权利要求1所述的方法，其特征在于：

上述对位步骤使用解剖学特征检测上述2张图像的对应像素位置，实施把一方图像和另一方图像对齐的变形(Warping)处理。

9、一种图像处理装置，其特征在于包括：

输入2张图像的图像输入单元；

对上述2张图像进行对位的对位单元；

使用在该对位装置中所对位的图像进行独立成分分析处理的独立成分分析处理单元。

10、一种图像处理***，其特征在于：

由权利要求9所述的图像处理装置经由网络与图像生成单元、显示单元、文件服务器中的至少1个连接而组成。

Images