[0001] 本发明涉及一种图像处理方法，更具体地涉及一种用于足底压力(pedobarographical)统计分析的方法。

[0002] 在各个医疗领域广泛使用对图像进行处理以使得获得临床上有价值的数据的技术。不同医疗领域关心对呈现人体不同部分的图像进行处理。图像处理所使用的图像以及所期望获得的临床信息的特性，使得在不同医疗领域使用不同技术。

[0003] 足底压力术是二维医疗成像技术，在步态的站立期期间以及姿势任务期间，对作用于对象的足趾面的压力场进行处理和比较。足底压力术使用在各种临床应用中，包括：术后评估、关节炎监测、矫正术设计以及糖尿病神经病变检测和管理。

[0004] 典型地，通过适当地对从位于对象的足下方的压力传感器获取的数据进行转换，来获得足底压力图像数据。压力传感器被布置为提供各种输出，这些输出指示具体区域中足所施加的压力。压力传感器提供产生的图像数据内的像素值。压力传感器的数目、间隔以及大小确定了图像数据的分辨率。典型地，以这种方式获得的图像具有亚厘米分辨率。

[0005] 当前分析足底压力图像的方法采用定性视觉评估，或者各种子采样技术之一，子采样技术对多个分离部位中的图像进行处理。

[0006] 尽管图像的视觉评估是高效的，但典型地，这种评估是耗时的并具有高度主观性。此外，只有经过适当训练的个人才能进行高效的视觉评估。

[0007] 用于分析足底压力图像的现有子采样技术对原始足底压力数据进行分割，并且根据解剖区域对分割后的数据进行标记，从而实现共同标记的区域的自动或半自动隔离。然后在标记的区域内和/或在标记的区域之间进行统计测试。

[0008] 尽管比视觉评估更客观，但是，子采样技术具有许多缺点。假定足趾面是连续的，而非离散的，子采样破坏了数据并最终提供了连续足-地面相互作用的低分辨率视图。子采样技术典型地仅涉及足底压力图像的具体区域的比较而不是将足作为整体。

[0057] 本发明的实施例使用逐像素或逐体素的统计图像比较，来提供一种自动使分析足底压力数据的方法。

[0058] 逐像素比较在数学上与逐体素的比较相同，并因此这里使用术语“像素”来表示二维像素和三维体素。

[0059] 为了以下描述的目的，统计图像涉及其中每个像素体现单个统计值的图像。例如，统计图像内在具体位置处的具体像素值可以对多个图像上相同位置处的像素的统计分布进行参数化。

[0060] 为了获得足底压力数据，获得来自于位于对象的足下方的传感器的传感器值输出。例如，每个传感器可以与产生的足底压力数据的单个像素相关联，使得每个传感器的输出值确定相应像素的像素值。在优选实施例中，从每个传感器中获得多个传感器值，并且使用最大传感器值来确定相应像素的像素值，但是也可以使用其他技术。

[0061] 在本发明的一些实施例中，在压力传感器与像素之间不存在直接的一对一映射，执行一些预处理，以将传感器值与像素值相关。例如，如果使用大小为5.08mm×7.62mm的传感器，可以使用双线性插值或其他插值方法将压力传感器值重新采样到方卡笛尔(Cartesian)网格上，以使得在例如5.08mm×5.08mm的虚拟坐标处产生像素值。

[0062] 在本发明的一些实施例中，图像可以是三维压力图像，其中，第三维度是时间。

[0063] 期望比较多个图像。可以通过比较从单个对象获得的多个图像(对象内比较)，以及通过比较从不同对象获得的多个图像(对象间比较)，来获得临床有用数据。在逐像素的基础上比较来自不同对象的图像之前，需要将每个图像变换至公共模板图像，使得每个图像具有公共大小、形状以及定向。

[0064] 图1示出了对输入图像执行的处理。在步骤S1，接收到多个输入图像，每个图像包括多个像素，每个像素具有关联的像素值。在图2A中示出了左脚的示例输入图像，其中，较亮像素指示压力较高的区域。现在描述对输入图像的处理。将认识到，对每个输入图像进行类似处理。

[0065] 由于足趾面是连续和顺应(Compliant)的，相邻像素之间不应有压力的剧烈改变，使得对图像进行平滑不会导致有害的信息丢失，而是会移除输入图像内不期望的伪像，并提高信噪比。可以使用本领域技术人员公知的任何平滑方法来对图像进行平滑。

[0066] 在图1所示的本发明的实施例中，通过滤波(步骤S2)、阈值处理(步骤S3)以及形态学运算(步骤S4)的组合来对输入图像进行平滑。

[0067] 在步骤S2，例如利用三像素直径(～1.5cm)的对称卷积核对每个输入图像进行滤波，该卷积核比中心像素多2倍。图2B示出了向图2A的图像应用如上所述对称卷积核的结果。

[0068] 可以看出，滤波处理的边界效应是，围绕足的外围添加了像素，并因此人工地使足接触区膨胀。因此，在图1的步骤S2处输入图像进行滤波之后，在步骤S3处，向图2B所示获得的图像应用阈值，以使得产生图2C所示的图像。所应用的阈值移除了足外围的像素。阈值被设置为，移除与预定阈值压力之下的压力值相对应的任何像素，选择阈值压力以指示原始足接触区外围的像素的像素值。

[0069] 参照图2C，可以看出，在步骤S3处应用阈值导致，在横向指骨的低压区域中存在spur像素1。隔离的像素和spur像素易于影响统计分析，并且期望通过形态学开来从模板图像中移除这些像素。在步骤S4中执行这样的操作，并且获得图2D所示图像的产生。

[0070] 先前描述已经涉及对输入图像进行平滑，以移除输入图像内的任何伪像。对输入图像进行平滑提高了信噪比，并且有助于方法中以下参照图1更详细描述的后续步骤。然而，将认识到，平滑不是这里所述的方法的本质方面。

[0071] 为了能够适当比较图像，在步骤S5处，将输入图像与公共模板进行配准，以使输入图像在公共模板上对准。此外，如果在相对短时间段内获得对象内的足底压力数据(因此对象的足的尺寸没有改变)，与公共模板的对准允许对象内分析，而无需足底压力图像的其他变换。以下描述可以在本发明的实施例中使用的一种重新对准方法。

[0072] 重新对准过程包括以下操作的组合：将输入图像平移到给定位置；旋转图像；在旋转后图像的最后平移之前，进一步平移以使输入图像返回其初始位置。

[0073] 将三个一般参数作为向量q的元素，变换序列T(q)可以定义如下：

[0074]

[0075]

[0076]

[0077]

[0078]

[0079] 其中：

[0080] (xc，yc)是足质心的坐标；

[0081] q1和q2指示了分别沿着水平和垂直方向的图像平移；以及

[0082] q3指示图像旋转。

[0083] 变换序列(5)，如矩阵(2)所指示地将图像平移至原点，如矩阵(3)所示关于图像的质心旋转图像，以及如矩阵(2)的逆所示将旋转后的图像平移回至其原始位置。然后，如矩阵(4)所示，将旋转后的图像平移至任意位置。

[0084] 变换序列(5)涉及重新定位图像的像素。因此，如果图像由向量p表示，具有正确像素值p’的变换后的图像由方程(6)给出：

[0085] p′＝p(XT) (6)

[0086] 其中：

[0087] X是同质形式(即，添加全1的虚拟行，以允许矩阵乘)的J个像素坐标的(3×J)阵列。

[0088] 图3A示出了从以外展足姿态行走的对象获取的图像2。也示出了要与图像2对准的模板图像3。图3A示出了图像2的第一主轴4和次主轴5。角度θ表示图像2的次主轴5相对于参考轴6的定向。图3B示出了图像2与模板图像3的配准结果。

[0089] 图4A和4B与图3A和3B相对应，但是涉及从以并拢足姿态行走的对象获取的足图像。

[0090] 给出源图像p1以及要与源图像配准的模板图像p0，确定提供最优变换的向量q。(0)
利用足的几何形状，可以根据方程(7)进行的q、q 初始估计：

[0091]

[0092] 其中：

[0093] (xc0，yc0)是模板图p0中足的质心；

[0094] (xc1，yc1)是源图像p1中足的质心；

[0095] θ0是模板图像p0中足的次主轴相对于参考轴的定向；以及

[0096] θ1是源图像p1中足的次主轴相对于参考轴的定向。

[0097] 可以将主轴计算为压力图像(即，其中，像素值表示局部‘质量’)的‘惯性’矩阵的本征向量，或者根据二值图像(即，其中所有像素具有相等的‘质量’)的本征向量。次主轴是相对于其存在最小惯性运动的轴。

[0098] q(0)可以用于定义方程(5)所表示的变换序列，并且如方程(6)所示，该变换序列可以应用于源图像p1以提供变换后的源图像p1’。现在使用以下不相似度量ε来直接比较模板图像p0和变换后的源图像p1：

[0099] ε＝(p1′-p0)T(p1′-p0) (8)

[0100] 其中，ε是两个图像的元素之间的差值的平方和，其中，上标‘T’指示向量转置。

[0101] 现在配准目标能够正式规定为：

[0102]

[0103] 优化问题(9)在3D参数空间Ω中是无限制并无约束的。可以使用任何优化技术(0)来实现优化。方程(8)的目标函数是方程(7)所给出的初始估计q 附近的凸起(convex)。
因此，例如，可以使用拟牛顿最速下降梯度搜索来实现优化。

[0104] 尽管图像配准是一种计算量很大的技术，但是利用如上所述的足几何形状允许产生良好的初始近似，以实现了具有相对简单优化技术的成功配准。已经发现，简单的梯度搜索技术在小数据集合中提供令人满意的收敛。即，使用误差平方和的梯度搜索对于较宽范围的对象和/或对于传统压力分布图不是有效的，特别地，这是由于已经发现方程(8)的目标函数仅在围绕初始近似的小区域中是凸起的。诸如稀疏梯度优化和二值形状批准等其他方法可以更适合于一般的配准方法。

[0105] 根据与公共模板对准的图像集合，能够对对象内图像进行逐像素统计分析，但是对象的足之间的形状差异意味着，期望在进行任何对象间分析之前，对图像进行空间归一化作为在图1的步骤S5处执行的配准的一部分。通过将来自于一个对象的源图像与从另一对象任意选择的模板图像进行配准，来实现空间翘曲。五参数非剪切仿射变换可以用于对沿着主轴测量的足图像的宽度和长度进行归一化。为了应用这样的变换，将方程(5)给出的变换序列修改为：

[0106]

[0107] 其中，

[0108]

[0109]

[0110] U0、R和U1分别由(2)、(3)和(4)给出

[0111] 这里，向量q有5个元素，前三个元素如上所述，两个附加元素q4和q5指示针对足图像的宽度和长度的相应缩放因子。方程(10)的变换序列表示线性空间翘曲的形式，能够将足的宽度和长度按缩放不相等的量(如果q4≠q5)。

[0112] 矩阵(11)给出的Rv将足旋转到垂直方向，使得矩阵(12)S可以将宽度按缩放因子q4，并将长度按缩放因子q5。开始进行优化的良好初始值q(0)由以下方程给出：

[0113]

[0114] 其中：

[0115] l0和l1是分别沿着次主轴5测量的、模板图像p0和源图像p1的足长度；以及[0116] w0和w1是分别沿着第一主轴4测量的、模板图像p0和源图像p1的足宽度。

[0117] 可以在方程(10)中使用q(0)来根据方程(6)对源图像p1进行变换，以使给出变换后的源图像p1’。现在根据方程(8)直接比较模板图像p0和变换后的源图像p1，使得根据方程(9)优化向量q的值。

[0118] 图5A示出了使用方程(10)变换之前从98kg男性对象获取的足图像7。图5B示出了使用方程(10)和(9)的优化空间翘曲之后的同一图像。这里，将图像7与从47kg女性对象获取的模板图像8进行配准。图5A和5B示出了如何经由配准将来自于具有不同足几何形状的对象的图像变换到公共模板空间。

[0119] 尽管使用由方程(8)表示的误差平方和来优化向量q的值，更直观的不相似矩阵可以用于获得配准性能的指示，以便进行报告。可以使用集合表示法来表达两个图像未交叠的程度：

[0120]

[0121] 其中，

[0122] p0和p1’是如上所述的图像；

[0123] |pi|是集合pi中非零元素的数目；

[0124] ⊕是逻辑XOR运算符；以及

[0125] ∨是逻辑OR运算符。

[0126] 方程(14)的分母指示由图像p0和图像p1’共同拥有的像素总数。方程(14)的分子指示由一个图像或另一个图像(而不是二者)覆盖的像素数目。当ε的值为0％时，给出理想配准。这会在方程(14)的分子具有零值时出现。当两个集合之间的交集是空(NULL)集合时，发生零交叠(即，ε等于100％)。方程(14)的分子是二值图像，并本身用于确认配准质量。

[0127]数据处理阶段 持续时间(s) 精确度(％)
平滑 0.006 N.A.
对象内配准 1.009 8.22
对象间配准 0.721 14.71
统计测试 2.855 N.A.
重采样 0.018 N.A.

[0128] 表1

[0129] 表1示出了这里所述方法的各个步骤的计算持续时间。可以看出，可忽略地、快速地进行平滑和重新采样。图像配准需要每图像处理1秒的量级。在配准上花费大部分总时间。由于优化步长大小的约束，对象间配准比对象内配准更快。

[0130] 配准精度是10％的量级，这指示平均起来近似10％的模板和源像素不交叠。对于对象内配准，仅在足外围上发现非交叠像素(如可以从以下进一步描述的图6中看出)。对于对象间配准，还在中间足弓和指骨的区域中找到非交叠像素。

[0131] 通过在配准图像集合上应用异或(XOR)运算来产生图6的图像。图6的图像示出了非交叠像素出现的频率，像素越白，在源图像的模板的交集中越不会频繁出现非交叠像素。

[0132] 近似90％的配准交叠比归因于由于足底压力分辨率而引起的非交叠像素的空间分布，并且还归因于解剖学因素。由于足底压力分辨率是5mm量级，预期仅以5mm的精度来绘制足的周长，甚至刚性控制对象应获得具有1像素厚度的XOR周长图像。这与实验间的变化性相联系说明了不理想的配准交叠。类似地，对象间解剖学差异能够说明，归一化之后在中间足弓和脚趾区域中观测到的变化性(如可以从图6中看出)。

[0133] 上述类型的仿射变换仅影响大体足形状，并且一般不能获得理想对象间配准。非线性空间翘曲可以用于改善性能。通过使用相对高次低频100参数离散余弦变换，非线性翘曲算法和鲁棒演进优化方法例如能够将交叠增加几个百分点。实际上，将认识到除了上述那些配准方法以外，可以使用对本领域技术人员而言显而易见的任何适合的配准方法。

[0134] 上述配准过程的目的在于，如图1的步骤S6至S9所示，允许在像素级别对图像直接进行统计测试。能够使用大量统计测试，所选的实际测试取决于应用。

[0135] 在步骤S6处，定义一般线性模型。定义的一般线性模型由被称作‘设计矩阵’的数字矩阵指定，该矩阵描述了与每个图像相关联的实验条件。例如，如果有来自于第一条件(例如，手术前)的5个图像，以及来自于第二条件(例如，手术后)的5个图像，则设计矩阵是10×2矩阵，其中，两个列分别表示手术前和手术后。手术前列在其前5行中具有1值，手术后列在其后5行中具有1值，并且所有其他值为零。这样的模型在两个实验条件之间进行正确区分，并且描述公知的2-样本t测试。除了这种简单模型以外，一般线性模型允许任意线性实验设计，包括使用例如，1样本t测试、方差分析(ANOVA)(F测试)、多元方差分析(MANOVA)、协方差分析(ANCOVA)、多元协方差分析(MANCOVA)、固定效果分析以及混合效果分析。

[0136] 图7是示例设计矩阵的图形表示。列表示实验因素，行表示实验重复。黑色区域具有零值，白色区域具有1值。图7的设计矩阵示出了设计9个对象和3个行走条件(正常行走、外翻(或外展)行走、以及反向(或并拢)行走)的实验。所有对象在每个条件下执行10次连续试验，如沿着垂直轴所示，总共获得270个图像。这是固定效果模型的示例，通过矩阵右侧上的对象块，每个对象效果被视为固定效果。这也可以是非随机实验设计的示例，其中，实验条件(前3列)包含单个重复块，而不具有在每对象的30次试验之间随机分散的试验。必须在步骤S6处指定诸如上述的设计，以使用一般线性建模方法。

[0137] 一般线性模型可以用于模型因素，模型因素不具有例如时间漂移的实验影响。这些被称作多余(nuisance)因素。图7中矩阵右侧上所列的对象因素是多余因素。在该示例中，对象之间的差异(例如，体重)不对仅与所关心的正常、外翻以及反向行走相关联效果有影响。

[0138] 一般线性模型是一种最灵活的线性统计建模方法，因此在这里对其进行描述，但是在步骤S6也可以使用其他统计模型。

[0139] 处理从步骤S6前进至步骤S7。在步骤S7处，对将实验条件映射至统计图像的参数进行估计。在2-样本t测试的情况下，例如，这些参数与两个实验条件的平均压力相对应。在这种情况下，对于统计图像中的每一像素应当有两个参数，加上针对每个像素的附加方差参数。使用如以下参照示例进一步描述的最小二乘方法来计算参数。

[0140] 处理从步骤S7前进至步骤S8。在步骤S8处，根据参数产生统计图像(以下参照示例进行描述)。在统计图像中，每个像素表示单个统计值。一般线性模型允许计算各种2
统计量，包括t统计(从而产生所谓的t图)、F统计(从产生所谓的F图)，以及χ 统计
2
(产生所谓的χ 图)。

[0141] 然后，处理前进至步骤S9，并且使用统计推断过程来评估统计值的统计显著性和/或统计图像中的空间处理。图8A和9A分别示出了任意统计图像，其中，每个像素具有关联的t值。由于存在多个像素，存在多个t值，因此t值本身形成图8A和9A所示的图像。

[0142] 单个t值不传达统计显著性。类似地，其中每个像素表示t值的统计图像不传达统计显著性。为了计算统计显著性，使用统计推断过程来将构成统计图像的t值变换成关联的概率(‘p’)值，从而产生如图1的步骤S11处所示的逐步统计推断图像。统计推断图像和相应p值的解释取决于所采用的统计推断过程的类型。例如，对于一种解释，1％的p值(p＝0.01)意味着，基于所观测的实验方差，有1％的机会，使得实验处理产生观测的效果。1％的p值的另一种解释是，如果随机过程产生数据，则应当期望，仅在1％的大量实验重复中观测到这种量级的t值。因此，当p值适当小时，可以合理地确信该实验效果不是由随机过程产生的。对于统计显著性的典型的临界截止(有时被称作‘α’)为p＝0.05。实验员-设置α值指定对于I类型统计误差的保护等级，也被称作‘错误肯定’，在空假设实际为真时拒绝空假设。

[0143] 图8B示出了二值图像，假设未校正的p阈值，该二值图像突出了图8A中具有小于p＝0.05的p值的像素。基于针对相应像素的t值来单独计算针对每个像素的p值。可以看出，由于每个像素的单独处理，阈值处理过程保留了许多像素。

[0144] 对于单个像素t值，以这种方式计算p值是经由单个公知的统计方程的直接一对一映射。但是该过程在统计上是无效的，这是由于其既忽略了已经进行许多统计测试的事实又忽略了图像数据的空间属性；相邻空间中的像素具有相关行为。因此，对于统计图像进行统计推断的有效方法涉及比上述公知的p值计算更复杂的技术。尽管可以使用各种方法来产生p值，但是这里描述具体适用于本发明的实施例的三类统计推断：(1)Bonferroni校正，(2)随机场理论，以及(3)非参数推断。

[0145] 在描述这些过程之前，重要的是认识到，所有这些推断过程的一般特性在于，它们应当解决多次比较的问题。计算统计图像(图8A和9A所示类型的统计图像)涉及，大量统计测试，对于每个图像像素进行一次测试。选择未校正的0.05的α值仅保护免于针对一个像素的机会过程，因此，未校正的p阈值通常导致统计图像处理的显著性的过高估计。因此应当考虑多次比较，以便保留有效α。

[0146] Bonferroni校正是用于对多次比较进行校正的最简单方法。Bonferroni校正减小了p阈值，使得类似误差率判断族在整个足面上保持在期望的α值(在所述示例中为5％)。Bonferroni校正计算为：

[0147] pctitical＝1-(1-α)1/K (15)

[0148] 其中，K是统计测试的数目，或者相当于统计图像中非零像素的数目。如上所述，图8B示出了图8A中具有足够大t值以致于低于0.05的未校正p阈值的那些像素(在这种情况下t＝1.73)。图8C示出了具有比t＝4.24的Bonferroni校正阈值大的t值的那些像素。可以看出，在与未校正推断过程有关的Bonferroni校正之后，更少的像素达到显著性。

[0149] 值得注意的是，未校正推断过程和Bonferroni推断过程均产生逐步连续统计推断图像，这些图像包含在超阈值聚类(具有大于根据α计算的临界t阈值的t值的像素聚类)的集合中的图8A所示的原始统计图像的t值。

[0150] 将认识到，图8B和8C中所示的统计推断图像在统计聚类之间存在间隙的意义上是逐步的。间隙表示图像不连续(即，t值的剧烈下降)，但t值在超阈值聚类内是连续的(即，在像素到像素的t值中没有剧烈改变)，因此称作“逐步连续”图像。在Bonferroni校正的情况下，并且假定实验员-设置α＝0.05，图8B所示的统计推断图像的一种解释在于，偶然出现每个超阈值聚类中的像素的概率为5％。

[0151] 上述Bonferroni校正的主要缺点在于，该Bonferroni校正没有适当考虑数据的空间属性。如果对500个像素的集合进行处理，Bonferroni校正应当是相同的，而不管那些像素的空间关系(像素在空间上可能随机分散至无限远)如何。然而，很明显，足压力图像和其他压力图像不是随机空间过程；它们具有连续足形状形式的空间顺序。

[0152] 随机场理论(RFT)利用下述事实，以(i)为多次比较提供更实际且不保守校正，或者(ii)向超阈值聚类而不是各个单独像素分配统计显著性。这两种RFT过程均取决于原始图像的两个特性：图像平滑度和空间几何形状。

[0153] 根据在一般信息模型参数的估计(图1的步骤S7)期间计算的残差图像，来估计图像平滑度。首先通过像素标准误差(s)(以下参照方程(34)像素标准误差(s)进行描述)对每个残差图像Ei(以下参照方程(35)对Ei的计算进行描述)进行归一化，以产生归一化的残差图像Zi：

[0154]

[0155] 接着，可以使用以下方程来计算残差空间导数的协方差矩阵Λ：

[0156]

[0157] 其中，k是非零像素的索引，其中，x和y指示空间方向，var是提供方差测量的函数，cov是提供协方差测量的函数。最后，全局平滑度由FWHM来估计：

[0158]

[0159] 这里，FWHM表示高斯核的“半峰全宽”，当与随机场数据卷积时，高斯核将产生与观测到的残差相同的平滑度。因此，FWHM提供图像平滑度的直接估计。

[0160] RFT推断过程也取决于可以由像素连通性在算法上定义的空间几何形状。对于二维图像而言，RFT相关空间几何形状由被称作‘分辨率元素计数’或‘resel计数’的三个参数R0、R1和R2来概括：

[0161] R0＝K-(Ex+Ey)+F (19)

[0162] R1＝(Ex+Ey-2F)/FWHM (20)

[0163] R2＝F/FWHM2 (21)

[0164] 其中，Ex和Ey是分别沿着x和y方向的相邻像素对的数目，F是搜索区域中4像素方块的数目。R参数因此定义了高达维数(在这种情况下是2维)的搜索空间，并且测量区域(R2)、周长(R1)以及欧拉特性(R0)。这些‘resel计数’参数可以视为图像数据集合中有效数目的独立观测(用于平滑度)。因此，场越平滑(FWHM越大)，独立观测的数目就越少。

[0165] 在计算(经由FWHM)图像平滑度以及计算(经由resel计数)图像几何形状之后，可以进行基于RFT的推断。如上所述，这可以采用两种方式之一来进行，(i)通过计算RFT校正的阈值以多次比较，或者(ii)通过遵照任意实验员-设置t阈值来计算超阈值聚类的统计显著性。

[0166] 先前RFT过程与Bonferroni校正的类似之处在于，计算临界t阈值h，具有大于阈值h的t值的像素被视作在统计上是显著的。阈值h计算如下：

[0167]

[0168]

[0169]

[0170]

[0171] 其中，RD是resel计数，ρD是欧拉特性密度，v是统计模型中自由度的数目(在以下方程(36)处示出了v的计算)，以及Γ是伽马函数。这里P(tmax＞h)是理论概率，在统计图像中最大t值大于阈值h。

[0172] 给定实验员-设置α，通过迭代改变参数t找到临界阈值h，使得方程(22)的右侧等于α。该过程获得与Bonferroni校正相比一般较小的保守阈值，产生如图8D和9B所示略微较宽的超阈值集合。假定α＝0.05，图8D所示的统计推断图像的统计显著性如下：给定观测图像平滑度和搜索区域的空间几何形状，偶尔出现这些超阈值像素有5％的概率。

[0173] 上述RFT推断过程的主要缺点在于，统计显著性是逐像素分配的，其过程没有利用超阈值像素的空间聚类(例如，图8D所示)。第二基于RFT的推断过程采用任意实验员-设置阈值h，然后计算偶尔出现特定超阈值聚类(在给定其空间范围的情况下)的概率。首先，预期数目的偏移集合像素(N)和聚类(m)计算如下：

[0174] E{N}＝Aρ0(h) (26)

[0175] E{m}＝Aρ2(h) (27)

[0176] 其中，A是由FWHM2归一化的搜索区域。最大超阈值聚类包含n或更多像素的概率[P(nmax＞n)]是：

[0177] P(nmax＞n)＝1-exp(-E{m}e-Bn) (28)

[0178] B＝Γ(2)E{m}/E{N} (29)

[0179] 然后基于超阈值聚类的大小k来计算每个超阈值聚类的这些概率值，这些概率值被指示为仅次于如图8E和9C所示统计推断图像上的聚类。假定α＝0.05，统计显著性的解释如下：“给定观测图像平滑度、空间几何形状以及阈值h，偶尔出现相同大小的超阈值统计聚类的概率为p％”。这里p是特定聚类的计算的p值。

[0180] 注意，为高斯场推导出第二基于RFT的推断过程的方程，从而仅为具有高自由度的t场给出近似解。例如在统计文献[Cao，J.，1999.Thesize of the connected2
components of excursion set of x，t and F fields.Advances in Applied Probability
31(3)579-595]中给出了用于t场和其他统计场的更精确结果。然而，上述方程对于高自由度的情况是适合的。

[0181] 最后一类推断过程是非参数推断。上述RFT过程的缺点在于，它们依赖于统计参数来描述图像随机性的基本期望。具体地，RFT使用像素级别‘均值’和‘标准偏差’参数来定义统计分布。这样的参数方法功能强，但是依赖于各种假定，包括残差的最严格正态分布(如方程35中所包括的)。当足压力图像包含有几何差异的区域(例如将足弓高的对象与足弓低的对象进行比较)时，会违背这种假定，这是由于对于足弓高的对象而言，在足弓没有接触地面的区域中存在许多空观测。

[0182] 对于这些情况，非参数(NP)推断更适合。以实验统计图像(‘E-图像’)开始，图8A、9A中示出了实验统计图像的示例，NP过程采用蒙特卡罗模拟，通过随机改变上述一般线性模型的创建期间应用的实验参数，来计算假定统计图像(‘H图像’)集合。对于E-图像和每个H-图像，计算统计值，例如，最大超阈值统计阈值的大小。H-图像和E-图像从而映射出所计算的统计量的实验概率分布。然后通过确定将E-图像的统计值与随机产生的H-图像值进行比较的极端程度，来根据该实验概率分布直接确定E-图像统计量的显著性。
假定α＝0.05，基于NP的推断图像可以解释如下：给定阈值h，偶尔出现相同大小的超阈值统计聚类的概率为p％。

[0183] 对于NP过程，可选地，可以对基本方差场进行平滑，以产生更平滑的统计图像(如在图8A所示的图像上所发生的)，能够在过度保守像素方差估计的情况下产生更大超阈值聚类的过程应当排除来自于这些聚类的特定像素。换言之，由于NP推断过程没有显式地取决于参数化方差，因此能够在不影响推断过程有效性的前提下对方差进行平滑。

[0184] 现在描述使用上述技术执行的实验。使用0.5米足扫描3D系统(可从RS Scan，Belgium获得)，从9个健康对象(5个男性，4个女性，年龄：29.8±7.7岁)收集左足足底压力记录。研究三个实验条件：正常足、外展足以及并拢足。不对其他步态参数(例如，行走速度)进行精确控制。对象在每个条件下执行10次连续试验，获得总共270个足底压力记录(每个对象30个记录)。从站立期时间序列中提取最大传感器值，可以形成针对次试验的单个图像。使用上述平滑、重新对准以及归一化技术来对具有任意大小、形状和定向的足图像进行变换。对其执行统计测试的图像包含在包含53×22网格中，该网格的742个像素包含了整个足。执行上述类型的对象间和对象内测试。

[0185]正常 外展 并拢
均值 0 -32.0 +26.3
标准偏差 (1.8) (7.9) (10.3)

[0186] 表2

[0187] 如主轴定向所测量的外展和并拢行走中的足定向与如表2所示近似30度的正常定向不同。假定非循环统计，由于较小角度值，单向ANOVA证实了实验处理的效果。

[0188] 如上所述，参照图1的步骤S6，产生一般线性统计模型：

[0189] pik＝gijβjk+hilγlk+eik (30)

[0190] 其中，pik是压力记录i中第k个像素处的压力观测。这里，i是I个压力记录的索引，k是包含足的K个像素的索引，j是J个实验因素的索引，以及l是L个(在该示例中为9个)对象的索引。假定误差(eik)在条件和对象上是独立、相等以及正态分布的，而在像素上并非如此。解释变量gij和hil分别与处理效果和多余因素(对象效果)相对应。βjk和γlk未知参数。

[0191] 如上所述，参照图1的步骤S7，为了比较实验条件之间的图像，必须首先求解出方程(30)中的未知参数。这可以使用线性代数容易进行。方程(30)的模型可以精确表达为：

[0192] P＝Gβ+e (31)

[0193] 其中，P是均值校正的压力数据矩阵，G是图7所示的二值实验设计矩阵，β是参数矩阵(包括多余因素)，以及e是误差矩阵。在像素内对数据P进行均值校正，使得该模型不需要多余的常数项。通常，由于可以将常数列添加至G以使得平均校正变得不必要，因此不需要对数据P进行平均校正。为了计算未知参数的值，β(表示为‘b’)的最小二乘估计通过以下方程来获得：

[0194] b＝G+P＝(GTG)-1GTP (32)

[0195] 其中，G+是G的伪倒置。在图1的步骤S8处，针对每个像素的t统计计算如下：

[0196]

[0197]

[0198] E＝(P-Gb)T(P-Gb) (35)

[0199] v＝I-rank(G) (36)

[0200] 其中，c是(1×(J+L))对比度向量，bk是b的第k个列，sk是针对每个像素的标准偏差估计，以及Ek是平方残差E的第k个对角元素。测试两个对比度向量：[-1 1 0 0]和[-1 0 1 0]，分别表示外展和并拢压力超过正常压力的情况；这里0是与8个对象多余因素相对应的(1×8)零向量。如上所述，t值本身形成被称为‘t图’的图像。该统计图像具有与原始足底压力图像相同的分辨率，并实际上具有原始足的2D形状。然而，与其像素值仅表示压力值的原始足底压力图像不同，t图的像素代表实验因素的线性组合的效果。

[0201] 图10A、10B和10C分别示出了在正常行走、外展行走以及并拢行走实验条件中获得的均值图像。可以看出，与图10A所示的正常行走相比，图10B示出了外翻行走中第一跖骨头和拇指下平均增加的压力。由图10C的图像表示的并拢行走与侧部跖骨头下的增加的压力相关联，并且还与侧部指骨下的横向足弓下的更多接触相关联。图11A和11B所示的t图证实了这些趋势的统计显著性。在图11A中可以看出，外展行走中足的中部上有增加的压力，而并拢行走中足的侧部下有增加的压力(如图11B所示)。如果进行基于RFT的推断，由于噪声聚类的大小较小，可以移除外围‘噪声’。可以看出，图11A和11B所示的图像以及图8E和9C所示的图像提供了一种用于比较足底压力图像的适宜机制。

[0202] 除了执行足底压力数据的临床/实验室分析以外，将认识到，上述用于足底压力图像分析的方法具有更一般的应用。例如，本发明的一个应用是，作为个人识别的方法，以例如在机场使用。可以保留个人数据库，该数据库为每个个人存储了所各处理后的图像，以及使用如上所述方法根据从该个人获取的多个足底压力图像而产生的统计图像。为了识别个人，获得其他足底压力图像，并与数据库中多个处理后的图像和统计图像进行比较，以确定是否存在匹配。

[0203] 图12是示出了可以使用本发明来执行以识别个人的处理的流程图，假定存在存储有个人详情以及从那些个人的足的足底压力图像获取的相应统计图像的数据库。

[0204] 参照图12，在步骤S12，获得个人的一只足或两只足之一的至少一个足底压力图像。可以使用任何适合的方法来获得压力图像，例如，通过让个人在包含适当传感器的压力板上行走。一旦获得了至少一个适合的图像，处理前进至步骤S13，在步骤S13，将获得的图像与模板图像进行配准，该模板图像已与数据库中所有对象图像配准。可以如上所述进行图像的配准。将认识到，参照附图1所述，在配准之前，可以对获得的图像进行平滑，以从获得的图像中移除伪像。处理从步骤S13前进至步骤S14，在步骤S14，将配准的图像与数据库中的每个图像和/或统计图像进行比较，如果需要，与处理后的图像进行比较，以确定最接近的匹配。可以使用模式识别技术(例如，k-最邻近算法(k-NN))来执行获得的图像与数据库中图像的比较。备选地，可以采用上述统计过程来确定匹配的置信度。将认识到，可以对返回的匹配执行经一步细化，以进一步缩小搜索结果。

[0205] 已经对样本大小为24个健康年轻对象执行参照图12描述的处理，对10个数据库分别重复进行，并且获得100％的识别精确度。

[0206] 其他应用能够分析对象与垫或轮椅之间的压力场，以检测对象的运动。这可以在本发明的实时实现方式中进行，其中，压力场的显著变化指示对象的运动，并从而有助于确定病人是否需要帮助。

[0207] 本发明还可以使用在车辆设计领域，例如，分析车门与车框之间的接触面，或者雨刷表面与挡风玻璃之间的接触面。实际上，将认识到一般而言，本发明能够被期望用于分析压力场的应用。

[0208] 这里描述的用于处理足底压力数据的技术还可以有用地应用于三维图像。具体地，现在描述实现了数据的三维空时体的可视化的技术，在图13中示出了这种可视化的示例。图14示出了被执行以产生图13的可视化的处理。

[0209] 参照图14，在步骤S15，接收到三维空时压力图像。这样的三维数据由任何足底压力硬件系统产生，足底压力硬件系统在时间上进行连续采样，并且典型的商业系统以500Hz记录2D图像。例如，人类行走期间0.6秒的典型站立期间获得大约300个2D图像。对这些2D图像进行堆叠以形成单个三维空时体。

[0210] 在步骤S16，通过将三维图像与三维高斯核进行卷积，来平滑接收到的图像。在步骤S17，通过对平滑后的图像进行算法插值来计算大量等值面，每个等值面与特定压力相关-2联。在图13中示出了三个等值面。第一等值面10表示表示1Ncm 压力的区域，第二等值面-2 -2
11表示5Ncm 压力的区域，第三等值面12表示10Ncm 压力的区域。在图14的步骤S18处，给予每个等值面关联的属性(例如，颜色、透明度以及纹理)。在步骤S19，可以设置照明、摄像机以及其他环境参数，以使得提供具有期望解释属性的可视化。然后在步骤S20处，以适当的形式向用户呈现和显示该可视化。尽管图13示出了来自于单次行走试验的压力值的空时体的可视化，但是相同的可视化过程可以应用于使用上述方法产生的3D统计图像，在图15中示出了这种可视化的示例。

[0211] 应理解，仅通过示例提供特定技术和变型，可以改变技术的顺序，以及也可以使用其他适合的技术。例如，统计分析技术仅用于示例性的。实际上，在不背离本发明的范围的前提下可以对描述实施例进行各种修改。