图像处理的制作方法

专利名称：图像处理的制作方法
技术领域：
本发明涉及一种图像处理方法，更具体地涉及一种用于足底压力 (pedobarographical) ^if ^tlf StJTjfeo
背景技术：
在各个医疗领域广泛使用对图像进行处理以使得获得临床上有价值的数据的技 术。不同医疗领域关心对呈现人体不同部分的图像进行处理。图像处理所使用的图像以及 所期望获得的临床信息的特性，使得在不同医疗领域使用不同技术。足底压力术是二维医疗成像技术，在步态的站立期期间以及姿势任务期间，对作 用于对象的足趾面的压力场进行处理和比较。足底压力术使用在各种临床应用中，包括术 后评估、关节炎监测、矫正术设计以及糖尿病神经病变检测和管理。典型地，通过适当地对从位于对象的足下方的压力传感器获取的数据进行转换， 来获得足底压力图像数据。压力传感器被布置为提供各种输出，这些输出指示具体区域中 足所施加的压力。压力传感器提供产生的图像数据内的像素值。压力传感器的数目、间隔 以及大小确定了图像数据的分辨率。典型地，以这种方式获得的图像具有亚厘米分辨率。当前分析足底压力图像的方法采用定性视觉评估，或者各种子采样技术之一，子 采样技术对多个分离部位中的图像进行处理。尽管图像的视觉评估是高效的，但典型地，这种评估是耗时的并具有高度主观性。 此外，只有经过适当训练的个人才能进行高效的视觉评估。用于分析足底压力图像的现有子采样技术对原始足底压力数据进行分割，并且根 据解剖区域对分割后的数据进行标记，从而实现共同标记的区域的自动或半自动隔离。然 后在标记的区域内和/或在标记的区域之间进行统计测试。尽管比视觉评估更客观，但是，子采样技术具有许多缺点。假定足趾面是连续的， 而非离散的，子采样破坏了数据并最终提供了连续足-地面相互作用的低分辨率视图。子 采样技术典型地仅涉及足底压力图像的具体区域的比较而不是将足作为整体。

发明内容
本发明的目的是避免或消除一个或多个上述问题。根据本发明的第一方面，提供了一种处理压力图像的方法，该方法包括接收多个 第一连续压力图像；对每一个第一连续压力图像进行处理，以产生相应第二连续压力图像； 以及通过对所述第二连续压力图像进行处理，来产生表示第一压力图像的特性的连续统计 图像。第一方面的优点在于，产生连续统计图像，从而允许针对连续压力图像而不是压 力图像的离散、子采样部分，执行统计分析。所述处理可以包括向第一压力图像应用平滑操作。所述平滑可以包括对接收到的图像进行滤波；可选地，向滤波后的图像应用预定阈值；以及可选地，对已经应用了预定阈值的图像执行形态学运算。例如，形态学运算可 以包括形态学开和闭。所述处理可以包括将第一压力图像与模板图像配准。将第一压力图像与模板图像配准可以包括执行第一压力图像的平移和旋转中的 至少一个。配准第一压力图像还可以包括产生初始配准参数，以及对所述初始配准参数进 行优化。参数可以定义线性或非线性空间变换。线性变换可以包括，例如应用剪切仿射变 换、线性空间翘曲(wrap)以及一般投影变换等步骤。非线性变换可以包括，例如向第一压 力图像应用参数化非线性空间翘曲变换等步骤。空间翘曲可以包括对压力图像进行缩放，并且所述缩放可以沿两个彼此垂直的方 向应用相应缩放因子。第一压力图像的特性可以是与一个或多个预定参考图像的比较。第一压力图像的特性可以是第一压力图像的比较。可以从单个对象或多个对象获取多个第一连续压力图像。产生数据可以包括产生表示第一压力图像的特性的连续图像。产生的图像实质 上可以具有与第一压力图像相同的分辨率。例如，如果特性是比较，则产生的图像可以是比 较图像。比较图像可以具有实质上等于或第一压力图像或每一个第一压力图像的分辨率。 这样的比较图像实现了第一压力图像或每一个第一压力图像的适宜比较。产生的连续统计图像可以是t_图、F_图、或由多个图像的一般线性统计模型产生 的其他统计图。所述方法还可以包括对统计图像进行处理以产生另一统计图像。产生另一统计 图像可以包括对统计图像执行统计推断。另一统计图像可以是逐步连续统计推断图像，表 示多个第一连续压力图像的特性的统计显著性。统计推断可以对由t测试产生的统计图像的t值进行处理，并产生概率值。统计推断可以考虑统计图像的平滑度和/或统计图像的至少一个几何属性。统计推断可以包括基于随机场理论、Bonferroni校正、和/或非参数推断的技 术。第一连续压力图像可以是足趾面足底压力图像，并且第一足趾面足底压力图像和 第二足趾面足底压力图像中的每一个可以是整个足趾面的图像。在第一压力图像是足趾面 足底压力图像的情况下，本发明的优点在于，对连续足趾面而不是足趾面的离散部分(如 现有子采样技术的情况)进行统计分析。所述方法还可以包括产生多个统计图像，多个统计图像中的每一个表示一个或 多个压力图像的特性；接收输入压力图像；以及基于输入压力图像选择多个统计图像之
o第一连续压力图像可以是第一对象和第二对象之间的压力场的图像。第一连续压力图像可以是对象和垫之间的压力场的图像。第一连续压力图像可以是对象和座椅之间的压力场的图像。本发明的第二方面提供了一种用于显示动态压力数据的方法，这样的数据包括 二维压力图像的连续时间样本，所述方法包括显示所述压力数据的三维体(volumn)表 示，其中，每个体由两个空间维度和一个时间维度组成。
本发明的第三方面提供了一种用于显示压力数据的方法，所述压力数据包括多 个二维压力数据集合，所述方法包括显示压力数据的三维体表示，其中，三维体表示的第 一和第二维度表示空间数据，所述三维体表示的第三维度表示时间数据。所述方法还可以包括以图形方式呈现三维体表示。所述方法还可以包括根据所述数据对表示预定压力或压力范围的表面进行计 算，以及产生指示确定的表面的图形数据。三维体表示可以是根据多个压力数据产生的三维统计图像。本发明可以以任何适宜的方式(包括通过适合的装置)来实现。这样的装置可以 是被编程为执行上述方法的计算机。本发明还提供了一种被布置为执行上述方法的计算机 程序。可以在适当载体介质上执行这样的计算机程序，这样的适当载体介质可以是有形载 体介质(例如，硬盘或⑶-ROM)，或者备选地，诸如通信信号等无形载体介质。将认识到，本发明的一个方面的上下文中描述的特征可以在本发明任何其他方面 的上下文中应用。


参照附图，通过示例描述本发明的实施例，在附图中图1是在本发明实施例中执行的处理的流程图；图2A是使用图1的处理对足底压力图像进行处理的示例；图2B是在执行作为图1的处理的一部分的滤波操作之后的图2A的足底压力图 像；图2C是在执行作为图1的处理的一部分的阈值处理操作之后的图2B的图像；图2D是在执行作为图1的处理的一部分的形态学开运算之后的图2C的图像；图3A和3B示出了以外展的足姿态行走的足图像的配准；图4A和4B示出了以并拢的足姿态行走的足图像的配准；图5A和5B示出了经由空间翘曲处理的配准；图6是示出了通过实验获得的图像序列的配准精确度的图；图7是示例实验设计矩阵；图8A至8E示出了使用多个统计推断方法产生的统计图像；图9A至9C示出了使用随机场理论产生的统计图像；图10A至10C分别示出了在正常行走、外展行走以及并拢行走的不同实验条件下 获得均值图像；图11A和11B分别示出了表示外展行走与正常行走的获得图像的统计比较以及表 示并拢行走与正常行走的获得图像的统计比较；图12是示出了这里在个人识别中描述的方法的使用的流程图；图13是使用等值面的足底压力数据的空时体可视化；图14是执行以产生足底压力数据的体可视化的处理的流程图；图15是三维空时统计图像的体可视化。
具体实施例方式本发明的实施例使用逐像素或逐体素的统计图像比较，来提供一种自动使分析足 底压力数据的方法。逐像素比较在数学上与逐体素的比较相同，并因此这里使用术语“像素”来表示二 维像素和三维体素。为了以下描述的目的，统计图像涉及其中每个像素体现单个统计值的图像。例如， 统计图像内在具体位置处的具体像素值可以对多个图像上相同位置处的像素的统计分布 进行参数化。为了获得足底压力数据，获得来自于位于对象的足下方的传感器的传感器值输 出。例如，每个传感器可以与产生的足底压力数据的单个像素相关联，使得每个传感器的输 出值确定相应像素的像素值。在优选实施例中，从每个传感器中获得多个传感器值，并且使 用最大传感器值来确定相应像素的像素值，但是也可以使用其他技术。在本发明的一些实施例中，在压力传感器与像素之间不存在直接的一对一映射， 执行一些预处理，以将传感器值与像素值相关。例如，如果使用大小为5.08mmX7.62mm 的传感器，可以使用双线性插值或其他插值方法将压力传感器值重新采样到方卡笛尔 (Cartesian)网格上，以使得在例如5. 08mmX5. 08mm的虚拟坐标处产生像素值。在本发明的一些实施例中，图像可以是三维压力图像，其中，第三维度是时间。期望比较多个图像。可以通过比较从单个对象获得的多个图像(对象内比较)，以 及通过比较从不同对象获得的多个图像(对象间比较)，来获得临床有用数据。在逐像素的 基础上比较来自不同对象的图像之前，需要将每个图像变换至公共模板图像，使得每个图 像具有公共大小、形状以及定向。图1示出了对输入图像执行的处理。在步骤S1，接收到多个输入图像，每个图像包 括多个像素，每个像素具有关联的像素值。在图2A中示出了左脚的示例输入图像，其中，较 亮像素指示压力较高的区域。现在描述对输入图像的处理。将认识到，对每个输入图像进 行类似处理。由于足趾面是连续和顺应(Compliant)的，相邻像素之间不应有压力的剧烈改 变，使得对图像进行平滑不会导致有害的信息丢失，而是会移除输入图像内不期望的伪像， 并提高信噪比。可以使用本领域技术人员公知的任何平滑方法来对图像进行平滑。在图1所示的本发明的实施例中，通过滤波(步骤S2)、阈值处理(步骤S3)以及 形态学运算(步骤S4)的组合来对输入图像进行平滑。在步骤S2，例如利用三像素直径( 1. 5cm)的对称卷积核对每个输入图像进行滤波， 该卷积核比中心像素多2倍。图2B示出了向图2A的图像应用如上所述对称卷积核的结果。可以看出，滤波处理的边界效应是，围绕足的外围添加了像素，并因此人工地使足 接触区膨胀。因此，在图1的步骤S2处输入图像进行滤波之后，在步骤S3处，向图2B所示 获得的图像应用阈值，以使得产生图2C所示的图像。所应用的阈值移除了足外围的像素。 阈值被设置为，移除与预定阈值压力之下的压力值相对应的任何像素，选择阈值压力以指 示原始足接触区外围的像素的像素值。参照图2C，可以看出，在步骤S3处应用阈值导致，在横向指骨的低压区域中存在 spur像素1。隔离的像素和spur像素易于影响统计分析，并且期望通过形态学开来从模板图像中移除这些像素。在步骤S4中执行这样的操作，并且获得图2D所示图像的产生。先前描述已经涉及对输入图像进行平滑，以移除输入图像内的任何伪像。对输入 图像进行平滑提高了信噪比，并且有助于方法中以下参照图1更详细描述的后续步骤。然 而，将认识到，平滑不是这里所述的方法的本质方面。为了能够适当比较图像，在步骤S5处，将输入图像与公共模板进行配准，以使输 入图像在公共模板上对准。此外，如果在相对短时间段内获得对象内的足底压力数据(因 此对象的足的尺寸没有改变)，与公共模板的对准允许对象内分析，而无需足底压力图像的 其他变换。以下描述可以在本发明的实施例中使用的一种重新对准方法。重新对准过程包括以下操作的组合将输入图像平移到给定位置；旋转图像；在 旋转后图像的最后平移之前，进一步平移以使输入图像返回其初始位置。将三个一 般参数作为向量q的元素，变换序列T(q)可以定义如下
其中(xc, y。)是足质心的坐标；
q1和q2指示了分别沿着水平和垂直方向的图像平移；以及
q3指示图像旋转。 变换序列(5)，如矩阵(2)所指示地将图像平移至原点，如矩阵(3)所示关于图像 的质心旋转图像，以及如矩阵(2)的逆所示将旋转后的图像平移回至其原始位置。然后，如 矩阵(4)所示，将旋转后的图像平移至任意位置。变换序列(5)涉及重新定位图像的像素。因此，如果图像由向量p表示，具有正确 像素值P’的变换后的图像由方程(6)给出p ‘ = p (XT)(6)
其中X是同质形式(即，添加全1的虚拟行，以允许矩阵乘)的J个像素坐标的(3XJ) 阵列。图3A示出了从以外展足姿态行走的对象获取的图像2。也示出了要与图像2对准 的模板图像3。图3A示出了图像2的第一主轴4和次主轴5。角度0表示图像2的次主 轴5相对于参考轴6的定向。图3B示出了图像2与模板图像3的配准结果。图4A和4B与图3A和3B相对应，但是涉及从以并拢足姿态行走的对象获取的足 图像。给出源图像Pl以及要与源图像配准的模板图像P(l，确定提供最优变换的向量q。 利用足的几何形状，可以根据方程(7)进行的q、q(°)初始估计 其中(xc0, yc0)是模板图pQ中足的质心；(xcl, ycl)是源图像Pl中足的质心；e^是模板图像P(1中足的次主轴相对于参考轴的定向；以及6工是源图像Pl中足的次主轴相对于参考轴的定向。可以将主轴计算为压力图像(即，其中，像素值表示局部‘质量’)的‘惯性’矩阵 的本征向量，或者根据二值图像(即，其中所有像素具有相等的‘质量’)的本征向量。次主 轴是相对于其存在最小惯性运动的轴。q(°)可以用于定义方程(5)所表示的变换序列，并且如方程(6)所示，该变换序列 可以应用于源图像Pl以提供变换后的源图像Pl’。现在使用以下不相似度量￡来直接比 较模板图像PQ和变换后的源图像Pi e = (Pl' -P0)T(P1' -p0)(8)其中，￡是两个图像的元素之间的差值的平方和，其中，上标‘T’指示向量转置。现在配准目标能够正式规定为
(9)优化问题(9)在3D参数空间Q中是无限制并无约束的。可以使用任何优化技术 来实现优化。方程⑶的目标函数是方程(7)所给出的初始估计q(°)附近的凸起(convex)。 因此，例如，可以使用拟牛顿最速下降梯度搜索来实现优化。尽管图像配准是一种计算量很大的技术，但是利用如上所述的足几何形状允许产 生良好的初始近似，以实现了具有相对简单优化技术的成功配准。已经发现，简单的梯度搜 索技术在小数据集合中提供令人满意的收敛。即，使用误差平方和的梯度搜索对于较宽范 围的对象和/或对于传统压力分布图不是有效的，特别地，这是由于已经发现方程(8)的目 标函数仅在围绕初始近似的小区域中是凸起的。诸如稀疏梯度优化和二值形状批准等其他 方法可以更适合于一般的配准方法。
根据与公共模板对准的图像集合，能够对对象内图像进行逐像素统计分析，但是 对象的足之间的形状差异意味着，期望在进行任何对象间分析之前，对图像进行空间归一 化作为在图1的步骤S5处执行的配准的一部分。通过将来自于一个对象的源图像与从另 一对象任意选择的模板图像进行配准，来实现空间翘曲。五参数非剪切仿射变换可以用于 对沿着主轴测量的足图像的宽度和长度进行归一化。为了应用这样的变换，将方程(5)给 出的变换序列修改为U^RySR'y'RU^其中，
U。、R和分别由(2)、(3)和(4)给出
这里，向量q有5个元素，前三个元素如上所述，两个附加元素q4和q5指示针对足
图像的宽度和长度的相应缩放因子。方程(10)的变换序列表示线性空间翘曲的形式，能够 将足的宽度和长度按缩放不相等的量(如果q4卓q5)。 矩阵(11)给出的Rv将足旋转到垂直方向，使得矩阵(12)S可以将宽度按缩放因 子q4，并将长度按缩放因子q5。开始进行优化的良好初始值q 由以下方程给出
一 Xcl
■VC0 一 yci 线一鸽
w0/w,
ik
“(13)其中10和“是分别沿着次主轴5测量的、模板图像P(l和源图像Pl的足长度；以及w0和Wl是分别沿着第一主轴4测量的、模板图像P(l和源图像Pl的足宽度。可以在方程(10)中使用q(°)来根据方程(6)对源图像？1进行变换，以使给出变换 后的源图像P/。现在根据方程(8)直接比较模板图像P(l和变换后的源图像Pl，使得根据 方程(9)优化向量q的值。图5A示出了使用方程(10)变换之前从98kg男性对象获取的足图像7。图5B示 出了使用方程(10)和(9)的优化空间翘曲之后的同一图像。这里，将图像7与从47kg女 性对象获取的模板图像8进行配准。图5A和5B示出了如何经由配准将来自于具有不同足
11几何形状的对象的图像变换到公共模板空间。尽管使用由方程(8)表示的误差平方和来优化向量q的值，更直观的不相似矩阵 可以用于获得配准性能的指示，以便进行报告。可以使用集合表示法来表达两个图像未交 叠的程度 其中，pQ和Pl，是如上所述的图像；|Pi是集合Pi中非零元素的数目； 是逻辑X0R运算符；以及V是逻辑OR运算符。方程(14)的分母指示由图像pQ和图像Pl’共同拥有的像素总数。方程(14)的分 子指示由一个图像或另一个图像(而不是二者)覆盖的像素数目。当￡的值为0%时，给出 理想配准。这会在方程(14)的分子具有零值时出现。当两个集合之间的交集是空(NULL) 集合时，发生零交叠(即，￡等于100% )。方程(14)的分子是二值图像，并本身用于确认 配准质量。 表1表1示出了这里所述方法的各个步骤的计算持续时间。可以看出，可忽略地、快速 地进行平滑和重新采样。图像配准需要每图像处理1秒的量级。在配准上花费大部分总时 间。由于优化步长大小的约束，对象间配准比对象内配准更快。配准精度是10%的量级，这指示平均起来近似10%的模板和源像素不交叠。对于 对象内配准，仅在足外围上发现非交叠像素(如可以从以下进一步描述的图6中看出)。对 于对象间配准，还在中间足弓和指骨的区域中找到非交叠像素。通过在配准图像集合上应用异或(X0R)运算来产生图6的图像。图6的图像示出 了非交叠像素出现的频率，像素越白，在源图像的模板的交集中越不会频繁出现非交叠像
o
12
近似90%的配准交叠比归因于由于足底压力分辨率而引起的非交叠像素的空间 分布，并且还归因于解剖学因素。由于足底压力分辨率是5mm量级，预期仅以5mm的精度来 绘制足的周长，甚至刚性控制对象应获得具有1像素厚度的X0R周长图像。这与实验间的 变化性相联系说明了不理想的配准交叠。类似地，对象间解剖学差异能够说明，归一化之后 在中间足弓和脚趾区域中观测到的变化性(如可以从图6中看出)。上述类型的仿射变换仅影响大体足形状，并且一般不能获得理想对象间配准。非 线性空间翘曲可以用于改善性能。通过使用相对高次低频100参数离散余弦变换，非线性 翘曲算法和鲁棒演进优化方法例如能够将交叠增加几个百分点。实际上，将认识到除了上 述那些配准方法以外，可以使用对本领域技术人员而言显而易见的任何适合的配准方法。上述配准过程的目的在于，如图1的步骤S6至S9所示，允许在像素级别对图像直 接进行统计测试。能够使用大量统计测试，所选的实际测试取决于应用。在步骤S6处，定义一般线性模型。定义的一般线性模型由被称作‘设计矩阵’的 数字矩阵指定，该矩阵描述了与每个图像相关联的实验条件。例如，如果有来自于第一条件 (例如，手术前)的5个图像，以及来自于第二条件(例如，手术后)的5个图像，则设计矩 阵是10 X 2矩阵，其中，两个列分别表示手术前和手术后。手术前列在其前5行中具有1值， 手术后列在其后5行中具有1值，并且所有其他值为零。这样的模型在两个实验条件之间 进行正确区分，并且描述公知的2-样本t测试。除了这种简单模型以外，一般线性模型允 许任意线性实验设计，包括使用例如，1样本t测试、方差分析(ANOVA) (F测试)、多元方差 分析(MAN0VA)、协方差分析(ANC0VA)、多元协方差分析(MANC0VA)、固定效果分析以及混合 效果分析。图7是示例设计矩阵的图形表示。列表示实验因素，行表示实验重复。黑色区域 具有零值，白色区域具有1值。图7的设计矩阵示出了设计9个对象和3个行走条件(正 常行走、外翻(或外展)行走、以及反向(或并拢)行走)的实验。所有对象在每个条件下 执行10次连续试验，如沿着垂直轴所示，总共获得270个图像。这是固定效果模型的示例， 通过矩阵右侧上的对象块，每个对象效果被视为固定效果。这也可以是非随机实验设计的 示例，其中，实验条件(前3列)包含单个重复块，而不具有在每对象的30次试验之间随机 分散的试验。必须在步骤S6处指定诸如上述的设计，以使用一般线性建模方法。一般线性模型可以用于模型因素，模型因素不具有例如时间漂移的实验影响。这 些被称作多余(nuisance)因素。图7中矩阵右侧上所列的对象因素是多余因素。在该示 例中，对象之间的差异(例如，体重)不对仅与所关心的正常、外翻以及反向行走相关联效 果有影响。一般线性模型是一种最灵活的线性统计建模方法，因此在这里对其进行描述，但 是在步骤S6也可以使用其他统计模型。处理从步骤S6前进至步骤S7。在步骤S7处，对将实验条件映射至统计图像的参 数进行估计。在2-样本t测试的情况下，例如，这些参数与两个实验条件的平均压力相对 应。在这种情况下，对于统计图像中的每一像素应当有两个参数，加上针对每个像素的附加 方差参数。使用如以下参照示例进一步描述的最小二乘方法来计算参数。处理从步骤S7前进至步骤S8。在步骤S8处，根据参数产生统计图像(以下参照 示例进行描述)。在统计图像中，每个像素表示单个统计值。一般线性模型允许计算各种统计量，包括t统计(从而产生所谓的t图)、F统计(从产生所谓的F图)，以及x 2统计 (产生所谓的x2图)。然后，处理前进至步骤S9，并且使用统计推断过程来评估统计值的统计显著性和 /或统计图像中的空间处理。图8A和9A分别示出了任意统计图像，其中，每个像素具有关 联的t值。由于存在多个像素，存在多个t值，因此t值本身形成图8A和9A所示的图像。单个t值不传达统计显著性。类似地，其中每个像素表示t值的统计图像不传达 统计显著性。为了计算统计显著性，使用统计推断过程来将构成统计图像的t值变换成关 联的概率(‘P’)值，从而产生如图1的步骤S11处所示的逐步统计推断图像。统计推断 图像和相应P值的解释取决于所采用的统计推断过程的类型。例如，对于一种解释，的 p值(P = 0.01)意味着，基于所观测的实验方差，有的机会，使得实验处理产生观测的 效果。的P值的另一种解释是，如果随机过程产生数据，则应当期望，仅在的大量实 验重复中观测到这种量级的t值。因此，当P值适当小时，可以合理地确信该实验效果不是 由随机过程产生的。对于统计显著性的典型的临界截止(有时被称作‘ CI ’)为p = 0. 05。 实验员-设置a值指定对于I类型统计误差的保护等级，也被称作‘错误肯定’，在空假设 实际为真时拒绝空假设。图8B示出了二值图像，假设未校正的p阈值，该二值图像突出了图8A中具有小于 P = 0. 05的p值的像素。基于针对相应像素的t值来单独计算针对每个像素的p值。可以 看出，由于每个像素的单独处理，阈值处理过程保留了许多像素。对于单个像素t值，以这种方式计算p值是经由单个公知的统计方程的直接一对 一映射。但是该过程在统计上是无效的，这是由于其既忽略了已经进行许多统计测试的事 实又忽略了图像数据的空间属性；相邻空间中的像素具有相关行为。因此，对于统计图像进 行统计推断的有效方法涉及比上述公知的P值计算更复杂的技术。尽管可以使用各种方法 来产生P值，但是这里描述具体适用于本发明的实施例的三类统计推断(l)B0nferr0ni校 正，(2)随机场理论，以及(3)非参数推断。在描述这些过程之前，重要的是认识到，所有这些推断过程的一般特性在于，它们 应当解决多次比较的问题。计算统计图像(图8A和9A所示类型的统计图像)涉及，大量 统计测试，对于每个图像像素进行一次测试。选择未校正的0.05的a值仅保护免于针对 一个像素的机会过程，因此，未校正的P阈值通常导致统计图像处理的显著性的过高估计。 因此应当考虑多次比较，以便保留有效a。Bonferroni校正是用于对多次比较进行校正的最简单方法。Bonferroni校正减 小了 P阈值，使得类似误差率判断族在整个足面上保持在期望的a值(在所述示例中为 5% )。Bonferroni校正计算为Pctitical = l"(l-a)1/K(15)其中，K是统计测试的数目，或者相当于统计图像中非零像素的数目。如上所述， 图8B示出了图8A中具有足够大t值以致于低于0. 05的未校正p阈值的那些像素(在这 种情况下t = 1. 73)。图8C示出了具有比t = 4. 24的Bonferroni校正阈值大的t值的那 些像素。可以看出，在与未校正推断过程有关的Bonferroni校正之后，更少的像素达到显著性。值得注意的是，未校正推断过程和Bonferroni推断过程均产生逐步连续统计推
接着，可以使用以下方程来计算残差空间导数的协方差矩阵A
(17) 其中，k是非零像素的索引，其中，x和y指示空间方向，var是提供方差测量的函 数，COV是提供协方差测量的函数。最后，全局平滑度由FWHM来估计 这里，FWHM表示高斯核的“半峰全宽”，当与随机场数据卷积时，高斯核将产生与观 测到的残差相同的平滑度。因此，FWHM提供图像平滑度的直接估计。RFT推断过程也取决于可以由像素连通性在算法上定义的空间几何形状。对于二 维图像而言，RFT相关空间几何形状由被称作‘分辨率元素计数’或‘resel计数’的三个参 数R0、礼和R2来概括R0 = K-(EX+Ey)+F(19)= (Ex+Ey-2F) /FWHM (20)R2 = F/FWHM2(21)
断图像，这些图像包含在超阈值聚类(具有大于根据a计算的临界t阈值的t值的像素聚 类)的集合中的图8A所示的原始统计图像的t值。将认识到，图8B和8C中所示的统计推断图像在统计聚类之间存在间隙的意义上 是逐步的。间隙表示图像不连续(即，t值的剧烈下降)，但t值在超阈值聚类内是连续的 (即，在像素到像素的t值中没有剧烈改变)，因此称作“逐步连续”图像。在Bonferroni 校正的情况下，并且假定实验员-设置a =0.05，图8B所示的统计推断图像的一种解释在 于，偶然出现每个超阈值聚类中的像素的概率为5%。上述Bonferroni校正的主要缺点在于，该Bonferroni校正没有适当考虑数据的 空间属性。如果对500个像素的集合进行处理，Bonferroni校正应当是相同的，而不管那 些像素的空间关系(像素在空间上可能随机分散至无限远)如何。然而，很明显，足压力图 像和其他压力图像不是随机空间过程；它们具有连续足形状形式的空间顺序。随机场理论(RFT)利用下述事实，以⑴为多次比较提供更实际且不保守校正，或 者(ii)向超阈值聚类而不是各个单独像素分配统计显著性。这两种RFT过程均取决于原 始图像的两个特性图像平滑度和空间几何形状。根据在一般信息模型参数的估计(图1的步骤S7)期间计算的残差图像，来估计 图像平滑度。首先通过像素标准误差(s)(以下参照方程(34)像素标准误差(s)进行描 述)对每个残差图像&(以下参照方程(35)对￡1的计算进行描述)进行归一化，以产生归 一化的残差图像Zi:
6 11
/k
15 其中，RD是resel计数，P D是欧拉特性密度，v是统计模型中自由度的数目(在以 下方程(36)处示出了 V的计算)，以及r是伽马函数。这里P(t_>h)是理论概率，在统 计图像中最大t值大于阈值h。给定实验员-设置a，通过迭代改变参数t找到临界阈值h，使得方程(22)的右 侧等于a。该过程获得与Bonferroni校正相比一般较小的保守阈值，产生如图8D和98所 示略微较宽的超阈值集合。假定a = 0. 05，图8D所示的统计推断图像的统计显著性如下 给定观测图像平滑度和搜索区域的空间几何形状，偶尔出现这些超阈值像素有5%的概率。上述RFT推断过程的主要缺点在于，统计显著性是逐像素分配的，其过程没有利 用超阈值像素的空间聚类(例如，图8D所示)。第二基于RFT的推断过程采用任意实验 员-设置阈值h，然后计算偶尔出现特定超阈值聚类(在给定其空间范围的情况下)的概 率。首先，预期数目的偏移集合像素(N)和聚类(m)计算如下 其中，艮和民是分别沿着x和y方向的相邻像素对的数目，F是搜索区域中4像 素方块的数目。R参数因此定义了高达维数(在这种情况下是2维)的搜索空间，并且测量 区域(R2)、周长(礼)以及欧拉特性(R。)。这些‘resel计数’参数可以视为图像数据集合中 有效数目的独立观测(用于平滑度)。因此，场越平滑(FWHM越大)，独立观测的数目就越 少。在计算(经由FWHM)图像平滑度以及计算(经由resel计数)图像几何形状之后， 可以进行基于RFT的推断。如上所述，这可以采用两种方式之一来进行，(i)通过计算RFT 校正的阈值以多次比较，或者(ii)通过遵照任意实验员-设置t阈值来计算超阈值聚类的 统计显著性。先前RFT过程与Bonferroni校正的类似之处在于，计算临界t阈值h，具有大于阈 值h的t值的像素被视作在统计上是显著的。阈值h计算如下
E{m}=Ap2(h)(27)其中，A是由FWHM2归一化的搜索区域。最大超阈值聚类包含n或更多像素的概率 [P(nmax > n)]是P (nmax > n) = 1-exp (_E {m} e_Bn)(28)B = T (2) E {m} /E {N}(29)然后基于超阈值聚类的大小k来计算每个超阈值聚类的这些概率值，这些概率值 被指示为仅次于如图8E和9C所示统计推断图像上的聚类。假定a =0.05，统计显著性的 解释如下“给定观测图像平滑度、空间几何形状以及阈值h，偶尔出现相同大小的超阈值 统计聚类的概率为”。这里P是特定聚类的计算的P值。注意，为高斯场推导出第二基于RFT的推断过程的方程，从而仅为具有高自 由度的t场给出近似解。例如在统计文献[Cao，J.，1999. Thesize of the connected components of excursion set of x2, t and F fields. Advances in Applied Probability 31 (3) 579-595]中给出了用于t场和其他统计场的更精确结果。然而，上述方程对于高自由 度的情况是适合的。最后一类推断过程是非参数推断。上述RFT过程的缺点在于，它们依赖于统计参 数来描述图像随机性的基本期望。具体地，RFT使用像素级别‘均值’和‘标准偏差’参数来 定义统计分布。这样的参数方法功能强，但是依赖于各种假定，包括残差的最严格正态分布 (如方程35中所包括的)。当足压力图像包含有几何差异的区域(例如将足弓高的对象与 足弓低的对象进行比较)时，会违背这种假定，这是由于对于足弓高的对象而言，在足弓没 有接触地面的区域中存在许多空观测。对于这些情况，非参数(NP)推断更适合。以实验统计图像(‘E-图像’)开始，图 8A、9A中示出了实验统计图像的示例，NP过程采用蒙特卡罗模拟，通过随机改变上述一般 线性模型的创建期间应用的实验参数，来计算假定统计图像(‘H图像’)集合。对于E-图 像和每个H-图像，计算统计值，例如，最大超阈值统计阈值的大小。H-图像和E-图像从而 映射出所计算的统计量的实验概率分布。然后通过确定将E-图像的统计值与随机产生的 H-图像值进行比较的极端程度，来根据该实验概率分布直接确定E-图像统计量的显著性。 假定a =0.05，基于NP的推断图像可以解释如下给定阈值h，偶尔出现相同大小的超阈 值统计聚类的概率为P%。对于NP过程，可选地，可以对基本方差场进行平滑，以产生更平滑的统计图像(如 在图8A所示的图像上所发生的)，能够在过度保守像素方差估计的情况下产生更大超阈值 聚类的过程应当排除来自于这些聚类的特定像素。换言之，由于NP推断过程没有显式地取 决于参数化方差，因此能够在不影响推断过程有效性的前提下对方差进行平滑。现在描述使用上述技术执行的实验。使用0. 5米足扫描3D系统(可从RS Scan, Belgium获得)，从9个健康对象(5个男性，4个女性，年龄29. 8士7. 7岁)收集左足足底 压力记录。研究三个实验条件正常足、外展足以及并拢足。不对其他步态参数(例如，行 走速度)进行精确控制。对象在每个条件下执行10次连续试验，获得总共270个足底压力 记录(每个对象30个记录)。从站立期时间序列中提取最大传感器值，可以形成针对次试 验的单个图像。使用上述平滑、重新对准以及归一化技术来对具有任意大小、形状和定向的 足图像进行变换。对其执行统计测试的图像包含在包含53X22网格中，该网格的742个像素包含了整个足。执行上述类型的对象间和对象内测试。
表2如主轴定向所测量的外展和并拢行走中的足定向与如表2所示近似30度的正常 定向不同。假定非循环统计，由于较小角度值，单向AN0VA证实了实验处理的效果。如上所述，参照图1的步骤S6，产生一般线性统计模型Pik = giJ ^ jk+h^ y lk+eik(30)其中，pik是压力记录i中第k个像素处的压力观测。这里，i是I个压力记录的 索引，k是包含足的K个像素的索引，j是J个实验因素的索引，以及1是L个(在该示例 中为9个)对象的索引。假定误差(eik)在条件和对象上是独立、相等以及正态分布的，而 在像素上并非如此。解释变量和分别与处理效果和多余因素(对象效果)相对应。
禾口 Ylk未知参数。如上所述，参照图1的步骤S7，为了比较实验条件之间的图像，必须首先求解出方 程(30)中的未知参数。这可以使用线性代数容易进行。方程(30)的模型可以精确表达 为P = G3+e(31)其中，P是均值校正的压力数据矩阵，G是图7所示的二值实验设计矩阵，0是参 数矩阵(包括多余因素)，以及e是误差矩阵。在像素内对数据P进行均值校正，使得该模 型不需要多余的常数项。通常，由于可以将常数列添加至G以使得平均校正变得不必要，因 此不需要对数据P进行平均校正。为了计算未知参数的值，0 (表示为‘b’)的最小二乘估 计通过以下方程来获得b = G+P = (GtG) _1GtP(32)其中，G+是G的伪倒置。在图1的步骤S8处，针对每个像素的t统计计算如下
cbttk =—
sk.(33)S2k ^^-Ag^YC7
V x ;(34)E = (P-Gb)T (P-Gb)(35)v = I-rank(G)(36)其中，c是(1X(J+L))对比度向量，bk是b的第k个列，sk是针对每个像素的标准 偏差估计，以及Ek是平方残差E的第k个对角元素。测试两个对比度向量[_1 1 0 0]和 [-1 0 1 0]，分别表示外展和并拢压力超过正常压力的情况；这里0是与8个对象多余因素 相对应的(1X8)零向量。如上所述，t值本身形成被称为‘t图’的图像。该统计图像具有
18与原始足底压力图像相同的分辨率，并实际上具有原始足的2D形状。然而，与其像素值仅 表示压力值的原始足底压力图像不同，t图的像素代表实验因素的线性组合的效果。图10A、10B和10C分别示出了在正常行走、外展行走以及并拢行走实验条件中获 得的均值图像。可以看出，与图10A所示的正常行走相比，图10B示出了外翻行走中第一跖 骨头和拇指下平均增加的压力。由图10C的图像表示的并拢行走与侧部跖骨头下的增加的 压力相关联，并且还与侧部指骨下的横向足弓下的更多接触相关联。图11A和11B所示的t 图证实了这些趋势的统计显著性。在图11A中可以看出，外展行走中足的中部上有增加的 压力，而并拢行走中足的侧部下有增加的压力(如图11B所示)。如果进行基于RFT的推 断，由于噪声聚类的大小较小，可以移除外围‘噪声’。可以看出，图11A和11B所示的图像 以及图8E和9C所示的图像提供了一种用于比较足底压力图像的适宜机制。除了执行足底压力数据的临床/实验室分析以外，将认识到，上述用于足底压力 图像分析的方法具有更一般的应用。例如，本发明的一个应用是，作为个人识别的方法，以 例如在机场使用。可以保留个人数据库，该数据库为每个个人存储了所各处理后的图像，以 及使用如上所述方法根据从该个人获取的多个足底压力图像而产生的统计图像。为了识别 个人，获得其他足底压力图像，并与数据库中多个处理后的图像和统计图像进行比较，以确 定是否存在匹配。图12是示出了可以使用本发明来执行以识别个人的处理的流程图，假定存在存 储有个人详情以及从那些个人的足的足底压力图像获取的相应统计图像的数据库。参照图12，在步骤S12，获得个人的一只足或两只足之一的至少一个足底压力图 像。可以使用任何适合的方法来获得压力图像，例如，通过让个人在包含适当传感器的压力 板上行走。一旦获得了至少一个适合的图像，处理前进至步骤S13，在步骤S13，将获得的图 像与模板图像进行配准，该模板图像已与数据库中所有对象图像配准。可以如上所述进行 图像的配准。将认识到，参照附图1所述，在配准之前，可以对获得的图像进行平滑，以从获 得的图像中移除伪像。处理从步骤S13前进至步骤S14，在步骤S14，将配准的图像与数据 库中的每个图像和/或统计图像进行比较，如果需要，与处理后的图像进行比较，以确定最 接近的匹配。可以使用模式识别技术(例如，k-最邻近算法(k-NN))来执行获得的图像与 数据库中图像的比较。备选地，可以采用上述统计过程来确定匹配的置信度。将认识到，可 以对返回的匹配执行经一步细化，以进一步缩小搜索结果。已经对样本大小为24个健康年轻对象执行参照图12描述的处理，对10个数据库 分别重复进行，并且获得100%的识别精确度。其他应用能够分析对象与垫或轮椅之间的压力场，以检测对象的运动。这可以在 本发明的实时实现方式中进行，其中，压力场的显著变化指示对象的运动，并从而有助于确 定病人是否需要帮助。本发明还可以使用在车辆设计领域，例如，分析车门与车框之间的接触面，或者雨 刷表面与挡风玻璃之间的接触面。实际上，将认识到一般而言，本发明能够被期望用于分析 压力场的应用。这里描述的用于处理足底压力数据的技术还可以有用地应用于三维图像。具体 地，现在描述实现了数据的三维空时体的可视化的技术，在图13中示出了这种可视化的示 例。图14示出了被执行以产生图13的可视化的处理。
参照图14，在步骤S15，接收到三维空时压力图像。这样的三维数据由任何足底压 力硬件系统产生，足底压力硬件系统在时间上进行连续采样，并且典型的商业系统以500Hz 记录2D图像。例如，人类行走期间0.6秒的典型站立期间获得大约300个2D图像。对这 些2D图像进行堆叠以形成单个三维空时体。在步骤S16，通过将三维图像与三维高斯核进行卷积，来平滑接收到的图像。在步 骤S17，通过对平滑后的图像进行算法插值来计算大量等值面，每个等值面与特定压力相关 联。在图13中示出了三个等值面。第一等值面10表示表示lNcm—2压力的区域，第二等值面 11表示5NcnT2压力的区域，第三等值面12表示lONcm—2压力的区域。在图14的步骤S18处， 给予每个等值面关联的属性(例如，颜色、透明度以及纹理)。在步骤S19，可以设置照明、 摄像机以及其他环境参数，以使得提供具有期望解释属性的可视化。然后在步骤S20处，以 适当的形式向用户呈现和显示该可视化。尽管图13示出了来自于单次行走试验的压力值 的空时体的可视化，但是相同的可视化过程可以应用于使用上述方法产生的3D统计图像， 在图15中示出了这种可视化的示例。应理解，仅通过示例提供特定技术和变型，可以改变技术的顺序，以及也可以使用 其他适合的技术。例如，统计分析技术仅用于示例性的。实际上，在不背离本发明的范围的 前提下可以对描述实施例进行各种修改。
权利要求
一种处理压力图像的方法，该方法包括接收多个第一连续压力图像；对所述第一连续压力图像中的每一个进行处理，以产生相应第二连续压力图像；以及通过对所述第二连续压力图像进行处理，来产生表示第一压力图像的特性的连续统计图像。
2.根据权利要求1所述的方法，其中，处理所述第一连续压力图像包括向所述第一压 力图像应用平滑操作。
3.根据权利要求2所述的方法，其中，向所述第一压力图像应用所述平滑操作包括对 接收到的图像进行滤波。
4.根据权利要求3所述的方法，其中，向所述第一压力图像应用所述平滑操作包括向 滤波后的图像应用预定阈值。
5.根据权利要求5所述的方法，其中，向所述第一压力图像应用所述平滑操作包括对 已经应用了预定阈值的图像执行形态学运算。
6.根据前述权利要求中任一项所述的方法，其中，处理所述第一压力图像包括将所 述第一压力图像与模板图像配准。
7.根据权利要求6所述的方法，其中，将所述第一压力图像与模板图像配准包括执行 第一压力图像的平移和旋转中的至少一个。
8.根据权利要求6或7中任一项所述的方法，其中，将所述第一压力图像与模板图像配 准包括产生初始图像变换参数，以及对所述初始图像变换参数进行优化。
9.根据权利要求6、7或8所述的方法，其中，处理所述第一压力图像还包括参照模板 图像对第一压力图像进行空间翘曲。
10.根据权利要求9所述的方法，其中，所述空间翘曲包括向所述第一压力图像应用 线性或非线性空间变换。
11.根据权利要求9或10所述的方法，其中，所述空间翘曲包括缩放所述压力图像。
12.根据权利要求11所述的方法，其中，所述缩放沿着两个彼此垂直的方向应用相应 缩放因子。
13.根据前述权利要求中任一项所述的方法，其中，所述第一压力图像的所述特性是与 一个或多个预定参考图像的比较。
14.根据前述权利要求中任一项所述的方法，其中，所述第一压力图像的所述特性是所 述第一压力图像的比较。
15.根据前述权利要求中任一项所述的方法，其中，所述多个第一压力图像是从单个对 象获取的。
16.根据权利要求1至14中任一项所述的方法，其中，所述多个第一压力图像是从多个 对象获取的。
17.根据前述权利要求中任一项所述的方法，其中，所产生的连续统计图像具有实质上 与所述第一压力图像或每个所述第一压力图像的分辨率相同的分辨率。
18.根据前述权利要求中任一项所述的方法，其中，所产生的连续统计图像是t-图、F_图、或其他统计图。
19.根据前述权利要求中任一项所述的方法，还包括 处理所述统计图像，以产生另一统计图像。
20.根据权利要求19所述的方法，其中，所述处理以产生另一统计图像包括对所述统 计图像执行统计推断。
21.根据权利要求20所述的方法，其中，所述另一统计图像是逐步连续统计推断图像， 表示多个第一连续压力图像的特性的统计显著性。
22.根据权利要求20或21所述的方法，其中，所述统计推断对由t测试产生的所述统 计图像的t值进行处理，并产生概率值。
23.根据权利要求20、21或22所述的方法，其中，所述统计推断考虑所述统计图像的平 滑度和/或所述统计图像的至少一个几何属性。
24.根据权利要求20至23中任一项所述的方法，其中，所述统计推断包括基于随机 场理论、Bonferroni校正、和/或非参数推断的技术。
25.根据前述权利要求中任一项所述的方法，其中，所述第一连续压力图像是足趾面足 底压力图像。
26.根据权利要求25所述的方法，其中，第一足趾面足底压力图像和第二足趾面足底 压力图像中的每一个是整个足趾面的图像。
27.根据权利要求25或26所述的方法，还包括产生多个统计图像，所述多个统计图像中的每一个表示一个或多个压力图像的特性； 接收输入压力图像；以及基于所述输入压力图像，选择所述多个统计图像之一。
28.根据权利要求27所述的方法，还包括存储所述多个所述统计图像中的每一个的标识数据；以及 输出所述多个统计图像中所选的一个的标识数据。
29.根据权利要求1至24中任一项所述的方法，其中，所述第一连续压力图像中的每一 个是第一对象和第二对象之间的压力场的图像。
30.根据权利要求29所述的方法，其中，所述第一连续压力图像中的每一个是对象和 垫之间的压力场的图像。
31.根据权利要求29所述的方法，其中，所述第一连续压力图像中的每一个是对象和 座椅装置之间的压力场的图像。
32.—种被配置为执行根据前述权利要求中任一项所述的方法的计算机程序。
33.一种承载根据权利要求25所述的计算机程序的计算机可读介质。
34.一种计算机装置，包括存储器，存储处理器可读指令；以及处理器，被配置为读取和执行存储在所述程序存储器中的指令； 其中，所述处理器可读指令包括控制处理器以执行根据权利要求1至25中任一项所 述的方法的指令。
35.一种用于显示包括多个二维压力数据集合在内的压力数据的方法，所述方法包括显示所述压力数据的三维体表示，其中，所述三维体表示的第一和第二维度表示空间 数据，以及所述三维体表示的第三维度表示时间数据。
36.根据权利要求35所述的方法，还包括以图形方式呈现所述三维体表示。
37.根据权利要求35或36所述的方法，还包括 根据表示预定压力或压力范围的所述数据，来计算表面； 产生对所确定的表面加以指示的图形数据。
38.根据权利要求35、36或37所述的方法，其中，三维体表示是根据多个所述压力数据 而产生的三维统计图像。
39.一种被配置为执行根据权利要求35至38中任一项所述的方法的计算机程序。
40.一种承载根据权利要求39所述的计算机程序的计算机可读介质。
41.一种计算机装置，包括存储器，存储处理器可读指令；以及处理器，被配置为读取和执行存储在所述程序存储器中的指令； 其中，所述处理器可读指令包括控制处理器以执行根据权利要求35至38中任一项所 述的方法的指令。
全文摘要
一种处理压力图像的方法，该方法包括接收多个第一连续压力图像；对所述第一连续压力图像中的每一个进行处理，以产生相应第二连续压力图像；以及通过对所述第二连续压力图像进行处理，来产生表示第一压力图像的特性的连续统计图像。
文档编号G06T7/00GK101933044SQ200880121950
公开日2010年12月29日 申请日期2008年12月19日 优先权日2007年12月21日
发明者托德·科林·保陶基, 约翰·扬尼斯·古勒马斯, 罗宾·休·克朗普顿 申请人:优莱企业有限公司