(57)【要約】  圧力画像を処理する方法であり、該方法は、複数の第１の連続する圧力画像を受信するステップと、前記第１の連続する圧力画像の各々を処理して第２の連続する圧力画像をそれぞれ生成するステップと、前記第２の連続する圧力画像を処理することによって前記第１の圧力画像の特性を表す連続する統計画像を生成するステップと、を備える。【特許請求の範囲】
【請求項１】
圧力画像を処理する方法であって、該方法は、
複数の第１の連続する圧力画像を受信するステップと、
前記第１の連続する圧力画像の各々を処理して第２の連続する圧力画像をそれぞれ生成するステップと、
前記第２の連続する圧力画像を処理することによって前記第１の圧力画像の特性を表す連続する統計画像を生成するステップと、
を備える、圧力画像の処理方法。
【請求項２】
前記第１の連続する圧力画像を処理するステップは、前記第１の圧力画像に平滑化処理を加えるステップを備える、請求項１記載の方法。
【請求項３】
前記第１の圧力画像に平滑化処理を加えるステップは、受信した画像をフィルタ処理するステップを備える、請求項２記載の方法。
【請求項４】
前記第１の圧力画像に平滑化処理を加えるステップは、更に、前記フィルタ処理された画像に所定の閾値化処理を加えるステップを備える、請求項３記載の方法。
【請求項５】
前記第１の圧力画像に平滑化処理を加えるステップは、更に、前記所定の閾値化処理が加えられた画像にモルフォロジー処理を実行するステップを備える、請求項４記載の方法。
【請求項６】
前記第１の圧力画像を処理するステップは、前記第１の圧力画像をテンプレート画像にレジストレーション（位置合わせ）するステップを備える、請求項１−５のいずれかに記載の方法。
【請求項７】
前記第１の圧力画像をテンプレート画像にレジストレーションするステップは、前記第１の圧力画像の平行移動及び回転の少なくとも一つを実行するステップを備える、請求項６記載の方法。
【請求項８】
前記第１の圧力画像をテンプレート画像にレジストレーションするステップは、
初期画像変換パラメータを発生させるステップと、
前記初期画像変換パラメータを最適化するステップと、
を備える、請求項６又は７記載の方法。
【請求項９】
前記第１の圧力画像を処理するステップは、更に、前記第１の圧力画像をテンプレート画像に対して空間的にワーピングするステップを備える、請求項６,７又は８記載の方法。
【請求項１０】
前記空間的ワーピングステップは、前記第１の圧力画像に線形又は非線形の空間変換を加えるステップを備える、請求項９記載の方法。
【請求項１１】
前記空間ワーピングステップは、前記圧力画像をスケーリングするステップを備える、請求項９又は１０記載の方法。
【請求項１２】
前記スケーリングは２つの直交方向に個別に倍率を与える、請求項１１記載の方法。
【請求項１３】
前記第１の圧力画像の前記特性は一以上の所定の基準画像との比較である、請求項１−１２のいずれかに記載の方法。
【請求項１４】
前記第１の圧力画像の前記特性は前記第１の圧力画像の比較である、請求項１−１３のいずれかに記載の方法。
【請求項１５】
前記複数の第１の連続する圧力画像は一被検者から得られる、請求項１０−１４のいずれかに記載の方法。
【請求項１６】
前記複数の第１の連続する圧力画像は複数の被検者から得られる、請求項１０−１４のいずれかに記載の方法。
【請求項１７】
前記生成される連続する統計画像は前記又は各第１の圧力画像にほぼ等しい解像度を有する、請求項１−１６のいずれかに記載の方法。
【請求項１８】
前記生成される連続する統計画像はｔ−マップ、Ｆ−マップ又は他の統計マップである、請求項１−１７のいずれかに記載の方法。
【請求項１９】
前記統計画像を処理して更なる統計画像を生成するステップを更に備える、請求項１−１８のいずれかに記載の方法。
【請求項２０】
前記更なる統計画像を生成するステップは、前記統計画像に統計的推定を実行するステップを備える、請求項１９記載の方法。
【請求項２１】
前記更なる統計画像は、複数の第１の連続する圧力画像の特性の統計的有意性を表す段階的な連続統計推定画像である、請求項２０記載の方法。
【請求項２２】
前記統計的推定は、t検定により生成される前記統計画像のｔ値を処理して確率値を生成する、請求項２０又は２１記載の方法。
【請求項２３】
前記統計的推定は、前記統計画像の平滑度及び／又は前記統計画像の少なくとも一つの幾何学特性を考慮する、請求項２０,２１又は２２記載の方法。
【請求項２４】
前記統計的推定は、確率場理論、ボンフェローニ補正及び／又は非パラメトリック推定に基づく技術を含む、請求項２０−２３のいずれかに記載の方法。
【請求項２５】
前記第１の連続する圧力画像は足底面のペドバログラフィ画像である、請求項１−１４のいずれかに記載の方法。
【請求項２６】
前記第１の足底面ペドバログラフィ画像及び前記第２の足底面ペドバログラフィ画像の各々は足底面全体の画像である、請求項２５記載の方法。
【請求項２７】
各々が一以上の圧力画像の特性を表す複数の統計画像を生成するステップと、
入力圧力画像を受信するステップと、
前記入力圧力画像に基づいて前記複数の統計画像の一つを選択するステップと、
を更に備える、請求項２５又は２６記載の方法。
【請求項２８】
前記複数の統計画像の各々に対して識別データを格納するステップと、
前記複数の統計画像の前記選択された一つに対して識別データを出力するステップと、
を更に備える請求項２７に記載の方法。
【請求項２９】
前記第１の連続する圧力画像の各々は、第１の物体と第２の物体との間の圧力場の画像である、請求項１−２４のいずれかに記載の方法。
【請求項３０】
前記第１の連続する圧力画像は、被検者とマットレスとの間の圧力場の画像である、請求項２９記載の方法。
【請求項３１】
前記第１の連続する圧力画像は、被検者とシーティング装置との間の圧力場の画像である、請求項２９記載の方法。
【請求項３２】
請求項１−３１のいずれかに記載の方法を実行するように構成されたコンピュータプログラム。
【請求項３３】
請求項３２のコンピュータプログラムを担持するコンピュータ読み取り可能媒体。
【請求項３４】
プロセッサ読み取り可能な命令を格納するプログラムメモリと、前記プログラムメモリに格納された命令を読み出し、実行するプロセッサとを備え、
前記プログラム読み取り可能命令は、前記プロセッサを、請求項１−２５のいずれかに記載の方法を実行するように制御する命令を含む、コンピュータ装置。
【請求項３５】
複数の組の２次元圧力データを含む圧力データを表示する方法であって、該方法は、
前記圧力データの３次元体積表現を表示するステップを備え、前記３次元体積表現の第１及び第２の次元は空間データを表し、前記３次元体積表現の第３の次元は時間データを表す、圧力データの表示方法。
【請求項３６】
前記３次元体積表現をグラフィック描画するステップを更に備える、請求項３５記載の方法。
【請求項３７】
所定の圧力又は所定の範囲の圧力を表すデータから表面を計算するステップと、
計算された表面を表示するグラフィックデータを生成するステップと、
を備える、請求項３５又は３６記載の方法。
【請求項３８】
前記３次元体積表現は前記複数の圧力データから生成される３次元統計画像である、請求項３５,３６又は３７記載の方法。
【請求項３９】
請求項３５−３８のいずれかに記載の方法を実行するように構成されたコンピュータプログラム。
【請求項４０】
請求項３９のコンピュータプログラムを担持するコンピュータ読み取り可能媒体。
【請求項４１】
プロセッサ読み取り可能命令を格納するプログラムメモリと、前記プログラムメモリに格納された命令を読み出し、実行するプロセッサとを備え、
前記プログラム読み取り可能命令は、前記プロセッサを、請求項３５−３８のいずれかに記載の方法を実行するように制御する命令を含む、コンピュータ装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は画像処理に関し、特にはペドバログラフィカル統計分析方法に関する。
【背景技術】
【０００２】
臨床的に有用なデータが得られるように画像を処理する技術は種々の医学分野で広く使用されている。種々の医学分野は人体の種々の部分を表す画像の処理と関連する。使用される画像及び画像の処理により得られる臨床情報の種類は、様々な技術が様々な医学の領域で使用されるほどである。
【０００３】
ペドバログラフィは、歩行の立脚期中及びポーズタスク中の被検者の足の底面に作用する圧力場を処理し比較する２次元医用画像工学である。ペドバログラフィは、手術後評価、関節炎監視、矯正学デザイン及び糖尿病性神経障害の検出及び管理を含む種々の臨床応用に使用されている。
【０００４】
ペドバログラフィ画像データは、典型的には、被検者の足の真下に位置する圧力センサから得られるデータを適切に変換することにより得ている。圧力センサは特定の領域内において足で与えられる圧力を表す可変出力を提供するように配置される。これらの圧力センサは発生される画像データ内のピクセルの値を提供する。圧力センサの数、それらの間隔及びそれらの大きさが画像データの解像度(分解能)を決定する。このようにして得られる画像は典型的にはサブセンチメートルの解像度である。
【０００５】
ペドバログラフィ画像を解析する現在の方法は定性的視覚評価又は画像を複数の個別部分で処理する種々のサブサンプリング技術の一つのいずれかを用いている。
【０００６】
画像の視覚評価は有効であるが、典型的には時間がかかり、極めて主観的である。加えて、有効な視覚評価は適切に訓練された人が実行できるのみである。
【０００７】
ペドバログラフィ画像の解析に使用される既存のサブサンプリング技術は、生ペドバログラフィデータをセグメント化し、セグメント化されたデータに解剖学的領域に従ってラベル付けし、一般的にラベル付けされた領域の自動的又は半自動的な分離を可能にしている。次にラベル付けされた領域内及び／又は領域間で統計的検定を行うことができる。
【０００８】
サブサンプリング技術は、視覚評価より客観的であるが、多くの欠点を被る。足の底面が連続的で分離していないものとすると、サブサンプリングはデータを破壊し、最終的には連続的な足と地面の相互作用の低解像度ビューをもたらす。サブサンプリング技術は典型的には足全体ではなく特定の領域のペドバログラフィ画像の比較と関連するのみである。
【発明の概要】
【発明が解決しようとする課題】
【０００９】
本発明の目的は、上述した問題の一つ以上を除去もしくは軽減することにある。
【課題を解決するための手段】
【００１０】
本発明の第１の態様によれば、圧力画像を処理する方法が提供され、該方法は、複数の第１の連続する圧力画像を受信するステップと、前記第１の連続する圧力画像の各々を処理して第２の連続する圧力画像をそれぞれ生成するステップと、前記第２の連続する圧力画像を処理することによって前記第１の圧力画像の特性を表す連続する統計画像を生成するステップとを備える。
【００１１】
この態様の利点は、連続する統計画像が生成される結果、圧力画像のサブサンプルされた個別の部分にではなく、連続する圧力画像に統計的解析を実行することが可能になる点にある。
【００１２】
前記第１の圧力画像の処理は、前記第１の圧力画像に平滑化処理を加えるステップを備えることができる。
【００１３】
前記平滑化処理は、受信した画像をフィルタ処理するステップと、要すればフィルタ処理された画像に所定の閾値化（２値化）処理を加えるステップと、要すれば所定の閾値化処理が加えられた画像にモルフォロジー処理を実行するステップを備えることができる。モルフォロジー処理は、例えばモルフォロジーオープニング及びクロージングを含むことができる。
【００１４】
前記第１の圧力画像の処理は、前記第１の圧力画像をテンプレート画像にレジストレーション（位置合わせ）するステップを備えることができる。
【００１５】
前記第１の圧力画像のテンプレート画像へのレジストレーションは、前記第１の圧力画像の平行移動及び回転の少なくとも一つを実行するステップを備えることができる。前記第１の圧力画像のレジストレーションは、更に、初期レジストレーションパラメータを発生させるステップと、前記初期レジストレーションパラメータを最適化するステップとを備える。前記パラメータは線形又は非線形空間変換のいずれかを定義するものとし得る。線形変換は、例えばシャーリングアフィン変換、線形空間ワーピング及び一般射影変換を加えるステップを含むことができる。非線形変換は、例えばパラメータ化された非線形空間ワーピング変換を第１の圧力画像に加えるステップを含むことができる。
【００１６】
空間ワーピングは、圧力画像をスケーリングするステップを備えることができ、スケーリングはそれぞれの倍率を２つの直交方向に与えることができる。
【００１７】
第１の圧力画像の特性は一以上の所定の基準画像との比較とすることができる。
【００１８】
第１の圧力画像の特性は第１の圧力画像の比較とすることができる。
【００１９】
複数の第１の連続する圧力画像は一つの被検者又は複数の被検者から得ることができる。
【００２０】
データの生成は、前記第１の圧力画像の特性を表す連続画像を生成するステップを備えることができる。生成される画像は前記第１の圧力画像とほぼ同一の解像度を有するものとし得る。例えば、特性を比較する場合、生成された画像を比較画像とすることができる。比較画像は前記又は各第１の圧力画像にほぼ等しい解像度を有することができる。このような比較画像は前記又は各第１の圧力画像の便利な比較を可能にする。
【００２１】
前記生成される連続する統計画像は、複数の画像の一般線形統計モデルによって生成さ得るｔ−マップ、Ｆ−マップ又は他の統計マップとすることができる。
【００２２】
本方法は、更に、前記統計画像を処理して更なる統計画像を生成するステップを備えることができる。更なる統計画像の生成は、前記統計画像に統計的推定を実行するステップを備えることができる。更なる統計画像は、複数の第１の連続する圧力画像の特性の統計的有意性を表す段階的な連続する統計推定画像とすることができる。
【００２３】
前記統計的推定は、t検定により生成された統計画像のｔ値を処理し、確率値を生成することができる。
【００２４】
前記統計的推定は、統計画像の平滑度及び／又は統計画像の少なくとも一つの幾何学特性を考慮することができる。
【００２５】
前記統計的推定は、確率場理論、ボンフェローニ補正及び／又は非パラメトリック推定に基づく技術を含むことができる。
【００２６】
前記第１の連続する圧力画像は足底面ペドバログラフィ画像とすることができ、第１の足底面ペドバログラフィ画像及び第２の足底面ペドバログラフィ画像の各々は足底面全体の画像とすることができる。第１の圧力画像が足底面ペドバログラフィ画像である場合、本発明は、統計的解析を既存のサブサンプリング技術の場合のように足底面の個別の区分ではなく、連続する足底面に実行する利点を有する。
【００２７】
本方法は、更に、各々が一以上の圧力画像の特性を表す複数の統計画像を生成するステップと、入力圧力画像を受信するステップと、入力圧力画像に基づいて前記複数の統計画像の一つを選択するステップとを備えることができる。
【００２８】
前記第１の連続する圧力画像は、第１の物体と第２の物体との間の圧力場の画像とすることができる。
【００２９】
前記第１の連続する圧力画像は、被検者とマットレスとの間の圧力場の画像とすることができる。
【００３０】
前記第１の連続する圧力画像は、被検者とシーティング装置との間の圧力場の画像とすることができる。
【００３１】
本発明の第２の態様によれば、２次元圧力画像の連続的な時間サンプルを含む動的圧力データを表示する方法が提供され、該方法は、前記圧力データの３次元体積表現を表示するステップを備え、各体積要素は２つの空間的次元と１つの時間的次元からなる。
【００３２】
本発明の第３の態様によれば、複数の組の２次元圧力データを含む圧力データを表示する方法が提供され、該方法は、圧力データの３次元体積表現を表示するステップを備え、前記３次元体積表現の第１及び第２の次元は空間データを表し、前記３次元体積表現の第３の次元は時間データを表す。
【００３３】
本方法は、更に、前記３次元体積表現をグラフィック表示するステップを備えることができる。
【００３４】
本方法は。更に、所定の圧力又は所定の範囲の圧力を表すデータから表面を計算するステップと、決定された表面を表示するグラフィックデータを生成するステップとを備えることができる。
【００３５】
前記３次元体積表現は複数の圧力画像から生成される３次元統計画像とすることができる。
【００３６】
本発明は、任意の都合のよい手段、例えば適切な装置を用いて実施することができる。このような装置は、上述した方法を実行するようにプログラムされたコンピュータとすることができる。本発明は、上述した方法を実行するように構成されたコンピュータプログラムも提供する。このようなコンピュータプログラムは適切なキャリア媒体上に坦持して有形キャリア媒体（例えばハードディスクまたはＣＤ−ＲＯＭ）とすることができ、また通信信号のような無形のキャリア媒体とすることができる。
【００３７】
本発明の一つの態様に関連して記載される特徴は本発明の他の態様のいずれにも適用可能である点に留意されたい。
【００３８】
本発明の一実施例が、例証として、添付図面を参照して以下に記載される。
【図面の簡単な説明】
【００３９】
【図１】本発明の一実施例で実行される処理のフローチャートである。
【図２】図２Ａは図１の処理を用いて処理されるペドバログラフィ画像の一例である。  図２Ｂは図１の処理の一部分として実行されるフィルタ処理後の図２Ａのペドバログラフィ画像である。  図２Ｃは図１の処理の一部分として実行される閾値化処理後の図２Ｂの画像である。  図２Ｄは図１の処理の一部分として実行されるモルフォロジーオープニング処理後の図２Ｃの画像である。
【図３】図３Ａ及び図３Ｂは外転した足姿勢での歩行を表す足画像のレジストレーションを示す。
【図４】図４Ａ及び図４Ｂは内転した足姿勢での歩行を表す足画像のレジストレーションを示す。
【図５】図５Ａ及び図５Ｂは空間ワーピング処理によるレジストレーションを示す。
【図６】実験的に得られた一連の画像間のレジストレーション精度を示す画像である。
【図７】模範的な実験計画行列を示す。
【図８】図８Ａ〜図８Ｅは複数の統計的推定法を用いて生成される統計画像を示す。
【図９】図９Ａ〜図９Ｃは確率場理論を用いて生成される統計画像を示す。
【図１０】図１０Ａ〜図１０Ｃは、正常歩行、外転歩行及び内転歩行の種々の実験条件において得られた平均画像をそれぞれ示す。
【図１１】図１１Ａ〜図１１Ｂは、外転歩行と正常歩行を表す画像の統計的比較及び内転歩行と正常歩行を表す画像の統計的比較をそれぞれ示す。
【図１２】ここに記載する方法を個人の識別に利用する例を示すフローチャートである。
【図１３】等値面を用いるペドバログラフィの時空間体積表示である。
【図１４】ペドバログラフィデータの３次元表示を生成するために実行される処理のフローチャートである。
【図１５】３次元時空間統計画像の３次元表示である。
【発明を実施するための形態】
【００４０】
本発明の実施例は、ペドバログラフィデータを解析するプロセスを自動化するためにピクセル単位又はボクセル単位の統計画像比較を用いる。
【００４１】
ピクセル単位の比較はボクセル単位の比較と数学的に同一であるため、ここでは２次元ピクセル及び３次元ボクセルの両方を表すために「ピクセル」を用いる。
【００４２】
以下の説明のために、統計画像は、各ピクセルが一つの統計値を具現している画像に関するものである。例えば、統計画像内の特定の位置の特定のピクセル値は複数の画像の同じ位置のピクセルの統計的分布をパラメータ化したものとし得る。
【００４３】
ペドバログラフィデータを得るために、被検者の足の真下に位置するセンサから出力されるセンサ値を得る。例えば、各センサは、発生されるペドバログラフィデータの単一ピクセルと関連し、各センサの出力値がそれぞれのピクセルのピクセル値を決定するようにすることができる。好適実施例では、各センサから複数のセンサ値が得られ、そのうちの最大センサ値を用いてそれぞれのピクセルに対するピクセル値を決定するが、他の技術を用いることもできる。
【００４４】
本発明のいくつかの実施例では、圧力センサとピクセルとの間の１対１の直接マッピングを用いる代わりに、前処理を実行してセンサ値をピクセル値に関連づける。例えば、5.08mm×7.62mmの大きさのセンサが使用される場合、圧力センサ値を双線形補間又は他の補間法を用いて方形デカルト格子上に再サンプルして、例えば5.08mm×5.08mmの仮想座標にピクセル値を発生させることができる。
【００４５】
本発明のいくつかの実施例では、画像は第３次元が時間である３次元圧力画像にすることができる。
【００４６】
複数の画像を比較するのが望ましい。臨床的に有用なデータは、単一被検者から得られた複数の画像の比較（被検者内比較）によって及び異なる被検者から得られた複数の画像の比較（被検者間比較）によって得ることができる。異なる被検者からの画像をピクセル単位方式で比較する前に、各画像を、各画像が共通のサイズ、形状及び方向を有するように共通のテンプレート画像に変換する必要がある。
【００４７】
図１は入力画像に実行される処理を示す。ステップＳ１において、複数の入力画像が受信され、各画像は複数のピクセルを含み、各ピクセルは関連するピクセル値を有する。左足の入力画像の一例が図２Ａに示されており、より明るいピクセルがより高い圧力の領域を示している。入力画像の処理が以下に説明される。各入力画像は同様に処理される点に留意されたい。
【００４８】
足の底面は連続的で従順であるので、隣接ピクセル間の圧力にシャープな変化が存在しないため、画像の平滑化は情報の有害な損失を生じることなく入力画像内の望ましくないアーチファクトを除去し、信号対雑音比を増大する。画像は当業者に周知の任意の平滑化方法を用いて平滑化することができる。
【００４９】
図１に示す本発明の実施例では、入力画像はフィルタリング（ステップＳ２）、閾値化(ステップＳ３)及びモルフォロジー処理（ステップＳ４）の組み合わせによって平滑化される。
【００５０】
ステップＳ２において、各入力画像は、例えば３ピクセル径（〜1.5cm）の対称コンボリューションカーネルでフィルタ処理され、中心ピクセルを２倍にバイアスする。図２Ｂは図２Ａの画像に上記の対称コンボリューションカーネルを適用した結果を示す。
【００５１】
このフィルタ処理の副作用は、足の周辺の周りにピクセルを追加し、ひいては足の接触域が人工的に膨張することにあることが明らかである。それゆえ、図１のステップＳ２で入力画像がフィルタ処理されたら、ステップＳ３で図２Ｂに示される得られた画像に閾値化を提供して図２Ｃに示される画像を発生させる。適用する閾値は足の周囲のピクセルを除去する。閾値は、所定の閾圧力より低い圧力値に対応するピクセルを除去するように設定され、閾圧力はもとの足接触域の周囲のピクセルのピクセル値を示すように選択される。
【００５２】
図２Ｃを参照すると、ステップ３における閾値化の適用は外側指骨の低圧力領域に突出ピクセル１を生じる。孤立した突出ピクセルは統計的解析に影響を与えやすく、これらはテンプレート画像からモルフォロジーオープニングにより除去すべきである。このようなオープニングはステップＳ４で実行され、図２Ｄに示される画像を生成する。
【００５３】
以上の記載は入力画像内のアーチファクトを除去する入力画像の平滑化と関連する。入力画像の平滑化は信号対雑音比を高め、図１を参照して以下に更に詳細に記載される方法の順次のステップを容易にすることができる。しかし、平滑化はここに記載する方法の必須の要件ではない点に留意されたい。
【００５４】
画像を適切に比較することを可能にするために、入力画像はステップＳ５において、入力画像が共通のテンプレートに整列するように、共通テンプレートにレジストレーション（位置合わせ）される。更に、被検者内ペドバログラフィデータが比較的短い期間（被検者の足の寸法が変化しない期間）の間に得られる場合には、共通テンプレートへの整列はペドバログラフィ画像の更なる変換の必要なしに被検者内解析を可能にする効果がある。本発明の実施例で使用できる再整列の一つの方法が以下に記載される。
【００５５】
レジストレーションプロセスは、所定の位置への入力画像の平行移動、画像の回転、回転された画像の最終平行移動前に入力画像を初期位置に戻すための更なる平行移動の組み合わせを用いる。
【００５６】
ベクトルｑの要素として３つの一般化パラメータを取ると、変換シーケンスＴ(q)は以下のように定義することができる。
【数１】
2011506022.tif 000002
【００５７】
ここで、
(xc, yc)は足の重心の座標であり、
q1及びq2はそれぞれ画像の水平方向及び垂直方向の平行移動を示し、
q3は画像の回転を示す。
【００５８】
変換シーケンス(5)は、画像を行列(2)で示されるように原点に平行移動させ、画像を行列(3)で示されるように重心を中心に回転させ、回転された画像を行列(2)の逆行列で示される元の位置に平行移動させて戻す。次に、回転された画像は行列(4)で示されるように任意の位置に平行移動される。
【００５９】
変換シーケンス(5)は画像のピクセルの再配置と関連する。従って、画像がベクトルｐで表される場合、適切なピクセル値p’を有する変換された画像は方程式(6)で与えられる。
p’= p(XT)                                        (6)
ここで、
Xは同次型のJピクセル座標の(3xJ)行列である（即ち行列乗算を可能にするために１が付加されたダミー行を有する）
【００６０】
図３Ａは、外転した足姿勢で歩く被検者から得られる画像２を示す。画像２を整列させるべきテンプレート画像３も示されている。図３Ａは画像２の長主軸４及び短主軸５を示す。角度θは基準軸６に対する画像２の短主軸５の方向を示す。図３Ｂは画像２のテンプレート画像３へのレジストレーションの結果を示す。
【００６１】
図４Ａ及び４Ｂは図３Ａ及び３Ｂに対応するが、内転した足姿勢で歩く被検者から得られる足の画像に関連する。
【００６２】
ソース画像p1と、該ソース画像をレジストレーションすべきテンプレート画像p0が与えられると、最適な変換を提供するベクトルqが決定される。足の幾何形状を利用して、qの初期推定q(0)を方程式(7)に従って行うことができる。
【数２】
2011506022.tif 000003
【００６３】
ここで、
(xc0, yc0)はテンプレート画像p0内の足の重心であり、
(xc1, yc1)はソース画像p1内の足の重心であり、
θ0は基準軸に対するテンプレート画像p0内の足の短主軸の方向であり、
θ1は基準軸に対するソース画像p1内の足の短主軸の方向である。
【００６４】
主軸は、圧力画像の「慣性」行列の固有ベクトルとして計算することができ（即ち、ピクセル値が局所的「質量」を表す場合）、また２値画像から計算することができる（即ち、すべてのピクセルが等しい質量を有する場合）。短主軸はその軸を中心とする慣性モーメントが最小になる軸である。
【００６５】
q(0)は方程式(5)で表される変換シーケンスを定義するために使用でき、変換シーケンスを方程式(6)で示されるようにソース画像p1に適用して変換されたソース画像p1’を提供することができる。このとき、テンプレート画像p0と変換されたソース画像p1’を、相違関数εを用いて直接比較することができる。
【数３】
2011506022.tif 000004
ここで、εは２つの画像の要素間の二乗誤差の和であり、上付き文字「Ｔ」はベクトル転置を示す。
【００６６】
そして、レジストレーションの目標は形式的に、
【数４】
2011506022.tif 000005
として記載することができる。
【００６７】
最適化問題(9)は３Ｄパラメータ空間Ωにおいて無限で非拘束である。最適化は任意の最適化技術を用いて実行することができる。方程式(8)の目的関数は方程式(7)で与えられる初期推定q(0)の近くで凸関数である。従って、最適化は、例えば準ニュートン最急降下勾配探索法を用いて実行することができる。
【００６８】
画像レジストレーションはコンピュータ要求技術であるが、上記のように足幾何形状を利用すると、良好な初期近似を生成することができ、比較的簡単な最適化技術で満足なレジストレーションを達成することができる。満足な収束を小さなデータセットで与えるために簡単な勾配探索法が見出されている。とは言うものの、二乗誤差の和を用いる勾配探索法は広範囲の被検者及び／又は画一的な圧力プロファイルに対して有効でないかもしれない。なぜなら、方程式(8)の目的関数は初期近似の周囲の小領域においてのみ凸関数であることが確かめられているためである。スパース勾配最適化及び２値形状レジストレーションなどの他のアプローチの方が一般化されたレジストレーションアプローチに適切であり得る。
【００６９】
共通テンプレートに整列された一組の画像に対して、被検者内画像にピクセル単位の統計的解析を行うことができるが、被検者間の足の形状の差は、被検者間解析を行う前に、図１のステップＳ５で実行されるレジストレーションの一部として画像を空間的に正規化するのが望ましいことを意味している。一つの被検者からのソース画像を別の被検者から任意に選択されたテンプレート画像にレジストレーションすることにより空間ワーピングを達成することができる。主軸に沿って測定された足画像の幅及び長さを正規化するために５パラメータ非せん断(non-shearing)アフィン変換を用いることができる。このような変換を適用するために、方程式(5)で与えられる変換シーケンスは、
【数５】
2011506022.tif 000006
に変更される。
【００７０】
ここで、
【数６】
2011506022.tif 000007
及びU0，R及びU1は式(2),(3)及び(4)で与えられる。
【００７１】
ここで、ベクトルqは５つの要素を有し、最初の３つの要素は上述したとおりであり、２つの追加の要素q4及びq5は足画像の幅及び長さに対するそれぞれのスケーリングファクタ(倍率)を示す。方程式(10)の変換シーケンスは線形空間ワーピングの一形式を表し、足の幅と長さを等しくない量（q4≠q5の場合）でスケーリングすることができる。
【００７２】
行列(11)で与えられるRvは足を垂直方向まで回転させるため、行列(12) Sは幅をファクタq4倍に、長さをファクタq5倍にスケーリングすることができる。最適化を進めることができる良好な開始値は、
【数７】
2011506022.tif 000008
で与えられる。
ここで、
l0及びl1は長主軸５に沿って測定されたテンプレート画像p0及びソース画像p1の足の長さであり、
wo及びw1は短主軸４に沿って測定されたテンプレート画像p0及びソース画像p1の足の幅である。
【００７３】
q(0)を方程式(10)に使用してソース画像p1を方程式(6)に従って変換し、変換されたソース画像p1'を与えることができる。この場合、テンプレート画像p0と変換されたソース画像p1'を方程式(8)に従って直接比較することができるため、ベクトルqの値を方程式(9)に従って最適化することができる。
【００７４】
図５Ａは、方程式(10)を用いる変換前に９８kgの男性被検者から得られる足画像７を示す。図５Ｂは、方程式(10)及び(9)を用いた最適空間ワーピング後の同じ画像を示す。ここでは画像７は４７kgの女性被検者から得られたテンプレート画像８にレジストレーションされている。図５Ａ及び図５Ｂは、様々な足幾何形状を有する被検者からの画像をレジストレーションにより共通のテンプレート空間に変換できることを示している。
【００７５】
方程式(8)で示される二乗誤差の和を用いてベクトルqの値を最適化することができるが、もっと直感的な相違行列を用いて報告目的のレジストレーション性能指示を得ることができる。２つの画像がオーバラップしない程度は集合表記法を用いて次のように書き表せる。
【数８】
2011506022.tif 000009
ここで、
p0及びp1'は上述した画像であり、
|pi｜は集合pi内の非ゼロ要素の数であり、
2011506022.tif 000010
【００７６】
方程式(14)の分母は画像p0及び画像p1'の両方が占めるピクセルの総数を示す。方程式(14)の分子は一方又は他方の画像で覆われ、両方の画像で覆われないピクセルの数を示す。完全なレジストレーションはεの値が０％のときに与えられる。これは、方程式(14)の分子がゼロ値を有するとき生じる。ゼロオーバラップ（即ちεが１００％に等しい）は２つの集合間の共通集合がヌル集合であるときに生じる。方程式(14)の分子は２値画像であり、レジストレーション品質を確かめるために有用である。
【００７７】
【表１】
2011506022.tif 000011
【００７８】
表１は、ここに記載する方法の種々のステップの計算持続時間を示す。平滑化及び再サンプリングは無視できるほど急速に行われることがわかる。画像レジストレーションは１画像につき１秒程度を要した。総時間の大部分がレジストレーションに費やされる。被検者間レジストレーションの方が最適化ステップのサイズ制限のために被検者内レジストレーションより速かった。
【００７９】
レジストレーション精度は１０％程度であり、平均で、テンプレートピクセルとソースピクセルの約１０％がオーバラップしなかった。非オーバラップピクセルは、被検者内レジストレーションに対してはもっぱら足の周囲に見られた（以下に更に詳しく記載される図６参照）。被検者間レジストレーションに対しては、非オーバラップピクセルは中央アーチ部及び指骨部にも見られた。
【００８０】
図６の画像は、レジストレーションされた画像の一群に排他的論理和（ＸＯＲ）演算を加えることによって生成される。図６の画像は非オーバラップピクセルの発生頻度を示し、白いピクセルほどテンプレート画像とソース画像の共通部分に頻繁に含まれなかったことになる。
【００８１】
約９０％のレジストレーションのオーバラップ比はペドバログラフィ分解能に起因する非オーバラップピクセルの空間分布に起因し、解剖学的要因にも起因する。ペドバログラフィ分解能は５ｍｍ程度であるため、足の境界は５ｍｍの精度でトレースできるのみであり、厳重な制御対象物でも１ピクセルの厚さのＸＯＲ境界画像を生じる。これは、試行間のバリアビリティと相まって、不完全なレジストレーションオーバラップを説明する。被検者間の解剖学的差も同様に正規化後の中央アーチ部及びつま先部に観測されるバリアビリティを説明する（図６参照）。
【００８２】
上述したタイプのアフィン変換は全体的な足形状にのみ影響し、一般に完全な被検者間レジストレーションをもたらさない。性能を向上させるために非線形空間ワーピングを使用することができる。比較的高次で低周波数の１００パラメータディスクリートコサイン変換非線形ワーピングアルゴリズム及びロバストな進化的最適化アプローチを用いることによって、オーバラップを数パーセント高めることができる。尚、上述したレジストレーション方法の代わりに任意の適切なレジストレーション方法を使用してもよいことは当業者に明らかであると点に留意されたい。
【００８３】
上述のレジストレーションプロシージャの目的は、図１のステップＳ６〜Ｓ９で示されるように、統計的検定を直接画像にピクセルレベルで行うことを可能にすることにある。一連の統計的検定を使用できるが、選択される実際の検定は用途により決まる。
【００８４】
ステップＳ６において、一般線形モデルが定義される。定義された一般線形モデルは、各画像が関連する実験状態を記述する「計画」行列という数値行列により規定される。例えば、第１の状態（例えば手術前）からの５つの画像及び第２の状態（例えば手術後）からの５つの画像がある場合、計画行列は１０×２行列とし、２つの列がそれぞれ手術前及び手術後を表す。手術前列は最初の５行に１の値を有し、手術後列は最後の５行に１の値を有し、他の値はすべて０である。このようなモデルは２つの実験状態間を区別し、周知の２標本ｔ検定(two-sample t-test)を記述する。この簡単なモデルに加えて、一般線形モデルは、例えば１標本ｔ検定(one-sample t-test)、分散の解析(ANOVA(analysis of variance))(F tests)、分散の多変量解析(MANOVA(multivariate analysis of variance))、共分散の解析(ANCOVA(analysis of covariance))、共分散の多変量解析(MANCOVA(multivariate analysis of covariance))、固定効果解析及び混合効果解析の使用を含む任意の線形実験計画を可とする。
【００８５】
図７は計画行列のグラフ表示である。列は実験的要因を表し、行は実験的反復を表す。黒領域は０の値を有し、白領域は１の値を有する。図７の計画行列は、９つの被検者と、正常歩行、外転歩行及び内転歩行の３つの歩行状態を含む実験を示す。各状態の１０回の連続試行がすべての被検者により実行され、垂直軸に沿って示されるように、全部で２７０の画像が生じる。これは固定効果モデルの一例であり、この例では各被検者効果が行列の右側の被検者ブロックを通して固定効果として処理される。これは非無作為実験計画法の一例でもあり、この例では実験状態（最初の３列）は各被検者の３０試行の中で無作為に分散される試行ではなく単一ブロックの反復を含む。一般線形モデル化アプローチを使用するためにこのような計画行列をステップＳ６で指定しなければならない。
【００８６】
一般線形モデルは、例えば時間ドリフトなどの実験的関心因子でない因子をモデル化するために使用できる。これらの因子は迷惑因子と呼ばれる。図７の計画行列の右側に記された被検者因子はヌイサンス因子である。この例では、被検者間の差（例えば体重など）は関心因子でなく、正常歩行、外転歩行及び内転歩行と関連する結果のみが関心因子である。
【００８７】
一般線形モデルは、線形統計モデル化に対して最もフレキシブルなアプローチであるとここでは記載しているが、ステップＳ６において他の統計モデルを使用することもできる。
【００８８】
処理はステップＳ６からステップＳ７に移る。ステップＳ７において、実験状態を統計画像にマップするパラメータが推定される。例えば、２標本ｔ検定の場合には、これらのパラメータは２つの実験状態の平均圧力に対応する。この場合には、統計画像内の各ピクセルにつき２つのパラメータが存在するのに加えて各ピクセルにつき追加の分散パラメータが存在する。これらのパラメータは、以下に一例につき更に詳細に記載されるように、最小二乗法を用いて計算される。
【００８９】
処理はステップＳ７からＳ８へ移る。ステップＳ８において、統計画像がパラメータから生成される（一例につき以下に記載される）。統計画像において、各ピクセルは単一の統計値を表す。一般線形モデルは種々の統計量、例えばｔ統計量（いわゆるｔマップを生成する）、Ｆ統計量（いわゆるＦマップを生成する）及びχ2統計量（いわゆるχ2マップを生成する）の計算を実行する。
【００９０】
次に処理はステップＳ９に移り、統計的推定処理を用いて統計画像内の統計値及び／又は空間プロセスの統計的有意性を推定する。図８Ａ及び図９Ａの各々は各ピクセルが関連するｔ値を有する任意の統計画像を示す。多数のピクセルが存在し、多数のｔ値が存在するため、ｔ値自体が図８Ａ及び図９Ａに示す画像を形成する。
【００９１】
単一ｔ値は統計的有意性を搬送しない。同様に、各ピクセルがｔ値を表す統計画像は統計的有意性を搬送しない。統計的有意性を計算するために、統計的推定処理を用いて統計画像を構成するｔ値を関連する確率(‘p’)値に変換し、それによって図１のステップＳ１１に示すように、段階的な統計的推定画像を生成する。統計的推定画像及び対応するｐ値の解釈は使用する統計的推定処理のタイプに依存する。例えば、一つの解釈では、１％のｐ値（p=0.01）は、観測された実験的分散に基づいて、実験的治療が観測効果を発生した可能性は１％であることを意味する。１％のｐ値の別の解釈は、確率過程(ランダム過程)によりデータが発生された場合、この大きさのｔ値は多数の実験的反復の僅か１％に観測されることが予測されることを意味する。それゆえ、ｐ値が適切に小さいとき、この実験効果は確率過程に起因しなかったことが合理的に確かである。統計的有意性の典型的な臨界カットオフ(‘α’と呼ばれることもある)はｐ＝０．０５である。実験者設定α値は、帰無仮説が実際に真であるとき帰無仮説を棄却してしまう「偽陽性(false-positive)」としても知られている第１種過誤に対する保護レベルを指定する。
【００９２】
図８Ｂは、ｐ＝０．０５より小さいｐ値を有する図８Ａのピクセルを強調する２値画像を示し、ｐ閾値は補正されてないものとする。各ピクセルに対するｐ値はそれぞれのピクセルに対するｔ値に基づいて個別に計算されている。閾値処理は各ピクセルの個別処理のために多くのピクセルを持ち続ける。
【００９３】
単一のピクセルｔ値に対するこのようなｐ値の計算は単一の周知の統計方程式による簡単な１対１マッピングである。しかし、この処理は、多くの統計的検定が行われた事実も、画像データの空間特性も、隣接空間内のピクセルが相関動作を有することも無視するため、統計的に無効である。従って、統計画像に統計的推定を行う有効な方法は、上述した周知のｐ値計算より複雑である技術を必要とする。ｐ値を生成するために種々の方法を使用できるが、ここでは特に本発明の実施例に適用可能な３つのクラスの統計的推定：(1)ボンフェローニの補正、(2)確率場(ランダム場)理論及び(3)非パラメトリック推定：を記載する。 
【００９４】
これらの処理について記載する前に、これらの推定処理のすべての一般的な特性は多数の比較の問題の原因となることを理解することが重要である。(図８Ａ及び図９Ａに示すタイプの)統計画像の計算は、各ピクセルにつき１つの多数の統計的検定を含む。０．０５の無補正のα値の選択は１つのピクセルに対する偶発過程を保護するのみであるため、無補正のｐ閾値は一般に統計画像プロセスの有意性の過大評価を生じる。それゆえ有効なαを維持するために、多重比較を考慮しなければならない。
【００９５】
ボンフェローニ補正は多重比較を補正する最も簡単な方法である。モンフェローニ補正は、足表面全体に亘ってファミリー誤差率が所望の値のα（上述の例では５％）に維持されるようにｐ閾値を減少させる。ボンフェローニ補正は、
ｐcritical＝１−(１−α)1/k                                              (15)
として計算される。
【００９６】
ここで、Kは統計的検定の数又は、等価的に、統計画像内の非ゼロピクセルの数である。上述したように、図８Ｂは０．０５の無補正の閾値ｐ以下に下がるのに十分なほど大きなｔ値（この場合にはｔ＝１．７３）を有する図８Ａ内のピクセルを示す。図８Ｃはｔ＝４．４２のモンフェローニ補正閾値しきい値より大きいｔ値を有するピクセルを示す。ボンフェローニ補正後に有意になるピクセルは無補正の推定処理に対してよりかなり少なくなることがわかる。
【００９７】
無補正及びボンフェローニの両推定処理は、図８Ａに示す生統計画像のｔ値を超閾値クラスタ（αから計算される臨界ｔ閾値を越えるｔ値を有するピクセルのクラスタ）の群内に含む段階的な連続統計画像を生成することに注目すべきである。
【００９８】
図８Ｂ及び図８Ｃに示す統計的推定画像は統計的クラスタ間にギャップがあるという意味で段階的である点に留意されたい。ギャップは画像の不連続性（即ち、ｔ値の鋭い低下）を示すが、ｔ値は超閾値クラスタ内では連続的であり（即ちｔ値のピクセル単位の鋭い変化はなく）、従って「段階的−連続的」画像という。ボンフェローニ補正の場合には、実験者設定α＝０．０５であるものとすると、図８Ｂに示す統計的推定画像の解釈は、「各超閾値クラスタ内のピクセルが偶然に発生する確率は５％である」である。
【００９９】
上述のモンフェローニ補正の主な弱点は、データの空間特性を適切に考慮しない点にある。５００ピクセルの集団を処理する場合、ボンフェローニ補正はこれらのピクセルの空間的関係に関係なく同一であり、ピクセルは空間内にランダムに無限に散乱される。しかし、足圧力画像及び他の圧力画像は明らかにランダムな空間プロセスでなく、これらの画像は一貫した足形状の形態の空間秩序を有する。
【０１００】
確率場理論（ＲＦＴ）は、(i)多重比較に対してよりリアルで低保存性の補正を提供すること又は(ii)統計的有意性を個々のピクセルではなく超閾値クラスタに割り当てることのいずれかを活用するものである。これらのＲＦＴ処理はいずれも元の画像の２つの特性、画像平滑度及び空間幾何形状、に依存する。
【０１０１】
画像平滑度は一般線形モデルパラメータの推定中（図１のステップＳ７）に計算される残存画像から推定される。各残存画像Ｅi（その計算は方程式（35）に関して以下に記載される）を最初にピクセル標準誤差（ｓ）（その計算は方程式（34）に関して以下に記載される）により正規化して、正規化された残存画像Ｚi：
Ｚi＝Ｅi／ｓ                                              (16)
を生成する。
【０１０２】
次に、残存空間導関数の共分散行列∧を、
【数９】
2011506022.tif 000012
を用いて計算することができる。
【０１０３】
ここで、ｋは非ゼロピクセルのインデックスであり、ｘ及びｙは空間方向を示し、varは分散の測定値を与える関数であり、covは共分散の測定値を与える関数である。最後に、全体的な平滑度はＦＷＨＭ：
【数１０】
2011506022.tif 000013
により推定される。
【０１０４】
ここで、ＦＷＨＭは、確率場データとともに畳み込まれると、観測される残余と同一の平滑度を生成する、ガウスカーネルの「半値幅(full width at half maximum)」を表す。従って、ＦＷＨＭは画像平滑度の直接推定を提供する。
【０１０５】
ＲＦＴ推定処理は、ピクセルの接続性によりアルゴリズム的に決定し得る空間幾何形状にも依存する。２次元画像に対しては、ＲＦＴ関連空間幾何形状は「分解素子カウント(resolution element counts)」又は「リセルカウント(resel counts)」という３つのパラメータＲ0，Ｒ1及びＲ2により要約することができる。
Ｒ0＝Ｋ−（Ｅx＋Ｅy）+Ｆ                                  (19)
Ｒ1＝（Ｅx＋Ｅy−２Ｆ）／ＦＷＨＭ                         (20)
Ｒ2＝Ｆ／ＦＷＨＭ2                                        (21)
【０１０６】
ここで、Ｅx及びＥyはそれぞれｘ方向及びｙ方向に隣接するピクセルの対の数であり、Ｆは探索区域内の４ピクセル方形の数である。従って、Ｒパラメータは探索空間を次元数（本例では２）まで規定し、面積（Ｒ2）、周囲（Ｒ1）及びオイラー特性（Ｒ0）を測定する。これらの「リセルカウント」パラメータは平滑度を考慮した画像データセット内の独立した観測の有効数とみなすことができる。従って、場が平滑になるほど（ＦＷＨＭが大きくなるほど）、独立した観測の数が少なくなる。
【０１０７】
（ＦＷＨＭによる）画像平滑度及び（リセルカウントによる）画像幾何形状の計算後に、ＲＦＴベースの推定を進めることができる。上述したように、この推定は、２つの方法のひとつ、即ち(i)多重比較を考慮するためにＲＦＴ補正閾値を計算することによって、又は(ii)任意の実験者設定ｔ閾値に続いて超閾値クラスの統計的有意性を計算することによって、生じさせることができる。
【０１０８】
前者のＲＦＴ処理はボンフェローニ補正に類似し、臨界ｔ閾値ｈを計算し、閾値ｈを超えるｔ値を有するピクセルを統計的に有意とみなす。閾値ｈは次のように計算される。
【数１１】
2011506022.tif 000014
【０１０９】
ここで、ＲDはリセルカウントであり、ρDはオイラー特性密度であり、νは統計モデルの自由度の数であり（その計算は下記の方程式（36）で示される）、Γはガンマ関数である。ここで、P(tmax > h)は統計画像内の最大ｔ値が閾値ｈを越える理論的確率である。
【０１１０】
実験者設定αが与えられれば、パラメータｔを方程式（22）の右辺がαに等しくなるように繰り返し変化させることによって臨界閾値ｈが見つかる。この処理は、ボンフェローニ補正と比較して全体的に低い保存性の閾値を生じ、図８Ｄ及び図９Ｂに示されるような僅かに広い超閾値設定を生じる。α＝０．０５であるとすると、図８Ｄに示される統計的有意性は次の通りであり、観測画像の平滑度及び探索区域の空間幾何形状が与えられれば、これらの超閾値ピクセルが偶然生起しうる確率は５％である。
【０１１１】
上述したＲＦＴ推定処理の主な弱点は、統計的有意性がピクセル単位で割り当てられ、そのプロセスは超閾値ピクセルの空間クラスタリング（例えば図８Ｄに示されている）を活用していない。第２のＲＦＴベースの推定処理は任意の実験者設定閾値ｈを用い、次に所与の空間広がりで特定の超閾値クラスタが偶然発生する確率を計算する。最初に、エクスカーション集合ピクセル（Ｎ）及びクラスタ（ｍ）の予測数を、
E{N}＝Aρ0(h)                                      (26)
E{m} =Aρ2(h)                                      (27)
として計算する。
【０１１２】
ここで、AはＦＷＨＭ2により正規化された探索区域である。最大の超閾値クラスタがｎ以上のピクセルを含む確率［P(nmax>n)］は、
P(nmax>n)＝１−exp(−E{m}e-Bn)                       (28)
B＝Γ(2) E{m}/ E{N}                                  (29)
である。
【０１１３】
次に、これらの確率値を各超閾値クラスタについてそれらのサイズｋに基づいて計算し、図８Ｅ及び図９Ｃに示されるように、これらの確率値を統計的推定画像上にクラスタに隣接して表示することができる。α＝０．０５であるものとすると、統計的有意性の解釈は次の通りであり、「観測された画像平滑度、空間幾何形状、及び閾値ｈが与えられれば、同じサイズの超閾値統計クラスタが偶然発生する確率はｐ％である」と解釈される。ここで、ｐは特定のクラスタの計算されたｐ値である。
【０１１４】
第２のＲＦＴベースの推定処理に対する方程式はガウス場に対して導出されたものであるため、ｔ場に対しては高い自由度で近似値を与えるのみである。ｔ場及び他の統計分野に対するもっと正確な結果は、統計学の文献、例えば「Cao,J., 1999, “The size of the connected components of excursion sets ofχ2,t and F fields”,Advances in Applied Probability 31(3), 579-595」に与えられている。しかし、上記の方程式は高い自由度の場合に適切である。
【０１１５】
最後のクラスの推定処理は非パラメトリック推定である。上述したＲＦＴ処理の潜在的な弱点は、画像ランダム性の基礎的予想を記述するために統計パラメータに依存している点にある。具体的にいうと、ＲＦＴは統計分布を決定するためにピクセルレベルの「平均」及び「標準偏差」パラメータを用いている。このようなパラメトリックアプローチは強力であるが、種々の仮定、例えば残部のノーマル分布（方程式３５に含まれる）に最も臨界的に依存する。足圧力画像が幾何学的不一致の区域を含むとき（例えば、高いアーチ構造の被検者と低いアーチ構造の被検者を比較するとき）、高いアーチ構造の被検者に対して地面と接触しない区域において多くのゼロ観測が存在するためにこの仮定は破られる。
【０１１６】
この種の状態に対しては、非パラメトリック（ＮＰ）推定がより適切である。実験的統計画像（Ｅ画像）（それらの例は図８Ａ，９Ａに示されている）から出発して、ＮＰ処理はモンテカルロシミュレーションを用い、上述した一般線形モデルの生成中に供給される実験的パラメータの順序をランダムに入れ替えることによって一群の仮説統計画像（Ｈ画像）を計算する。Ｅ画像及び各Ｈ画像に対して、統計値、例えば最大の超閾値統計クラスタが計算される。これによりＨ画像及びＥ画像は計算された統計値の実験的確率分布を決定する。この場合、Ｅ画像の有意性は、その統計極値がランダムに生成されるＨ画像の値に匹敵するかどうかを決定することによって、この実験的確率分布から直接決定される。α＝０．０５であるものとすると、ＮＰベースの推定画像は次の通りに解釈され、即ち閾値ｈが与えられれば、同じサイズの超閾値統計クラスタが偶然発生する確率はｐ％である、と解釈される。
【０１１７】
ＮＰ処理に対しては、基礎的分散場を随意に平滑化してより滑らかな統計画像（図８Ａに示される画像とともに発生される）を生成することができ、過度に保存的なピクセル分散の推定の場合にはより大きな超閾値クラスタを生成しうるプロセスはこれらのクラスタからいくつかのピクセルを排除する。換言すると、ＮＰ推定処理はパラメータ化された分散に明確に依存しないため、推定処理の有効性に影響を与えることなく分散を平滑化することができる。
【０１１８】
次に上述した技術を用いて実行した実験について記載する。左足のペドバログラフィ記録を、０．５メートルのフットスキャン３Ｄシステム（ＲＳスキャン社（ベルギー）から市販されている）を用いて、９人の健康な被検者（５人の男性、４人の女性、年齢２９．８±７，７才）から収集した。３つの実験状態、即ち正常歩行足、外転歩行足及び内転歩行足、を調査した。他の歩行パラメータ（例えば歩行速度）は精密には制御しなかった。被検者は各状態の試行を連続的に１０回行い、全部で２７０（各被検者につき３０）のペドバログラフィ記録を生じさせた。各試行につき立脚期時系列から最大センサ値を抽出して一つの画像を形成した。任意のサイズ、形状及び向きの足画像を、上述した平滑化、再位置調整及び正規化技術を用いて変換した。統計的検定が実行された画像は５３×２２の格子内に含まれ、そのうちの７４２ピクセルが足に含まれる。上述したタイプの被検者間及び被検者内検定を実施した。
【０１１９】
【表２】
2011506022.tif 000015
【０１２０】
外転及び内転歩行における足の向き（主軸の向きで測定される）は、表２に示されているように、正常方向から約３０度相違している。小さい角度値のために非円形統計量であるものと仮定すると、ワンウェイＡＮＯＶＡは実験的処理の効果を確認した。
【０１２１】
図１のステップＳ６に関連して記載したように、一般線形統計モデル：
ｐik＝ｇijβjk＋ｈilγlk＋ｅik                    (30)
が生成され、ここでｐikは圧力記録ｉ内のｋ番のピクセルにおける圧力観測値である。ここで、ｉはＩ個の圧力記録のインデックスを示し、ｋは足を含むＫ個のピクセルのインデックスを示し、ｊはＪ個の実験因子のインデックスを示し、ｌはＬ（本例では９）人の被検者のインデックスを示す。誤差（ｅik）は、実験状態間及び被検者間で互いに独立で同等であり且つ正規分布するものと考えられるが、ピクセル間ではそうではない。説明変数ｇij及びｈilはそれぞれ処理効果及びヌイサンス因子（被検者効果）に対応する。βjk及びγlkは未知のパラメータである。
【０１２２】
図１のステップＳ７に関連して記載したように、画像を実験状態間で比較するためには、最初に方程式(30)の未知のパラメータを解く必要がある。方程式(30)のモデルは簡潔に、
Ｐ＝Ｇβ＋ｅ                                      (31)
と書き表せる。
【０１２３】
ここで、Ｐは平均補正された圧力データ行列であり、Ｇは図７に示される２値実験的計画行列であり、βはパラメータ行列（ヌイサンス因子を含む）、及びｅは誤差行列である。データＰはピクセル内で平均補正されるため、モデルは余分の定数項を必要としない。一般に、１の定数列をＧに加えることにより平均補正を不要にできるので、データＰを平均補正する必要はない。未知のパラメータの値を計算するために、βの最小二乗推定（ｂで示す）を、
【数１２】
2011506022.tif 000016
によって得る。ここで、Ｇ＋はＧの擬似反転である。
【０１２４】
図１のステップＳ８において、各ピクセルについてｔ統計値を次のように計算する。
【数１３】
2011506022.tif 000017
【０１２５】
ここで、ｃは（l*(J+L)）コントラストベクトルであり、ｂkはｂのｋ番の列であり、ｓkは各ピクセルの標準誤差推定であり、Ｅkは方形残余Ｅのｋ番の対角要素である。外転及び内転圧力が正常圧力を超える場合をそれぞれ示す２つのコントラストベクトル、[-1 1 0 0]及び[-1 01 0]、を試験した。ここで、０は８つの被検者ヌイサンス因子に対応する（１×８）ヌルベクトルである。上述したように、ｔ値自体が「ｔマップ」という画像を形成する。この統計画像は元のペドバログラフィ画像と同じ解像度を有し、元の足の２Ｄ形状を有する。しかし、ピクセル値が圧力値のみを表す元のペドバログラフィ画像と異なり、ｔマップのピクセルは実験的因子の一次結合の効果を表す。
【０１２６】
図１０Ａ，１０Ｂ及び１０Ｃはそれぞれ正常歩行、外転歩行及び内転歩行の実験状態において得られる平均画像を示す。図１０Ｂは、外転歩行においては、図１０Ａに示される正常歩行に比較して、第一中足骨頭及び第一趾の下部内側に増大した圧力を示す。図１０Ｃの画像データ表される内転歩行は外側中足骨頭の下部に増大した圧力を有するとともに外側指骨の下部及び縦アーチの下部に大きな接触面積を有する。これらの傾向の統計的有意性は図１１Ａ及び図１１Ｂに示されるｔマップにより確認される。図１１Ａから、外転歩行においては足内側部分に増大した圧力が存在することがわかるが、(図１１Ｂに示されるように)内転歩行においては足外側部分に増大した圧力が存在する。ノイズクラスタのサイズは小さいので、ＲＦＴベースの推定を行えば、周辺ノイズを除去できる。図１１Ａ及び１１Ｂに示される画像及び図８Ｅ及び９Ｃに示される画像はペドバログラフィ画像を比較する便利なメカニズムを提供する。
【０１２７】
ペドバログラフィデータの臨床／研究解析を実行するのに加えて、上述したペドバログラフィ画像解析方法はもっと一般的な応用を有する。例えば、本発明の一つの応用として、例えば空港で使用する個人識別方法がある。個人のデータベースを維持し、該データベースに、各個人について、各個人から取り出した複数のペドバログラフィ画像から上述の方法により生成された複数の処理画像及び統計画像を格納する。個人を識別するために、更にペドバログラフィ画像を取り出し、データベースに格納された複数の処理画像及び統計画像と比較して一致するかどうか決定する。
【０１２８】
図１２は個人を識別するように実行できる本発明によるプロセスを示すフローチャートであり、各個人の詳細を該個人の足のペドバログラフィ画像から取り出された統計画像と一緒に格納したデータベースが存在するものとする。
【０１２９】
図１２を参照すると、ステップＳ１２において個人の足の少なくとも一つのペドバログラフィ画像が得られる。圧力画像は、任意の適切な方法を用いて、例えば適切なセンサを含む圧力板の上を個人に歩いてもらうことにより得ることができる。少なくとも一つの適切な画像が得られたら、処理はステップＳ１３に移り、このステップにおいて、得られた画像を、データベース内のすべての被検者画像がレジストレーションされているテンプレートにレジストレーションする。画像のレジストレーションは上述したように行うことができる。図１に関連して記載したように、レジストレーション前に、得られた画像を平滑化して得られた画像からアーチファクトを除去することができる点に留意されたい。ステップＳ１３から、処理はステップＳ１４に移り、このステップにおいて、レジストレーションされた画像をデータベース内の各画像及び／又は統計画像及び必要に応じ各処理画像と比較して最も一致するものを決定する。得られた画像とデータベース内の画像との比較は、パターン認識技術、例えばｋ近接アルゴリズム（k-NN）を用いて実行することができる。あるいは、上述した統計手法を使用して一致の信頼性を決定することもできる。更に精密な反復一致を実行して探索結果を更に狭くすることができる点に留意されたい。
【０１３０】
図１２について記載したプロセスを、２４人の健康な若い被検者の標本サイズについて、各人につき１０回のデータベース繰り返しで実行し、１００％の識別精度を得た。
【０１３１】
別の応用としては、被検者とマットレス又は車椅子との間の圧力場を解析して被検者の動きを検出することができる。これは、本発明のリアルタイム実施で達成でき、この場合には圧力場の有意変化が被検者の動きを示し、患者が補助を必要としているかどうかを決定するのに役立つ。
【０１３２】
本発明は、更に、車両設計の分野で、例えば車のドアと車のフレームとの界面又はワイパー表面とフロントガラス面との界面の解析に使用することができる。一般に、本発明は、圧力場を解析する能力が必要とされる場合に適用可能である点に留意されたい。
【０１３３】
ここに記載のペドバログラフィデータの処理技術は３次元画像にも有効に適用できる。特に、３次元の時空間体積のデータの視覚化を可能にする技術が以下に記載され、その一例が図１３に示されている。図１４は図１３の視覚化を生じさせるために実行される処理を示す。
【０１３４】
図１４を参照すると、ステップＳ１５において、３次元の時空間圧力画像を受信する。このような３次元データは時間的に連続してサンプリングする任意のペドバログラフィハードウェアシステムにより生成され、典型的な市販のシステムは２Ｄ画像を５００Ｈｚで記録する。例えば、人間の歩行中の０．６秒の典型的な立脚期は約３００の２Ｄ画像を発生する。これらの２Ｄ画像を積み重ねて単一の３次元時空間体積を形成することができる。
【０１３５】
ステップＳ１６において、受信した３次元画像を３次元ガウスカーネルで畳み込むことによって平滑化する。ステップＳ１７において、平滑化された画像をアルゴリズム的に補間することによって各々が特定の圧力と関連する複数の等値面を計算する。３つの等値面が図１３に示されている。第１の等値面１０は１Ｎcm-2の圧力を有する区域を表し、第２の等値面１１は５Ｎcm-2の圧力を有する区域を表し、第３の等値面１２は１０Ｎcm-2の圧力を有する区域を表す。図１４のステップＳ１８において、各当表面に関連特性（例えば色、透明度及びテクスチャ）を付与する。ステップＳ１９において、明暗、カメラ及び他の環境パラメータを、所望の解釈可能特性を有する視覚化が提供されるように設定する。次に、ステップＳ２０において、この視覚化がレンダリングされ、ユーザに適切な形で表示される。図１３は単一の歩行試行からの圧力値の時空間体積の視覚化を示すが、同じ視覚化処理を上述した方法を用いて生成された３Ｄ統計画像に適用することができ、例えばその一例が図１５に示されている。
【０１３６】
特定の技術及び変換が例証として記載されているのみであり、これらの技術の順序は変更可能であり、他の適切な技術も同様に使用できることを理解されたい。例えば、統計解析技術は例証としてのみ示されている。実際上、本発明の範囲から外れることなく多くの変更を記載の実施例に対して加えることができる。
【図１】
2011506022.tif 000018
【図２Ａ】
2011506022.tif 000019
【図２Ｂ】
2011506022.tif 000020
【図２Ｃ】
2011506022.tif 000021
【図２Ｄ】
2011506022.tif 000022
【図３Ａ】
2011506022.tif 000023
【図３Ｂ】
2011506022.tif 000024
【図４Ａ】
2011506022.tif 000025
【図４Ｂ】
2011506022.tif 000026
【図５Ａ】
2011506022.tif 000027
【図５Ｂ】
2011506022.tif 000028
【図６】
2011506022.tif 000029
【図７】
2011506022.tif 000030
【図８Ａ−８Ｅ】
2011506022.tif 000031
【図９Ａ−９Ｃ】
2011506022.tif 000032
【図１０Ａ】
2011506022.tif 000033
【図１０Ｂ】
2011506022.tif 000034
【図１０Ｃ】
2011506022.tif 000035
【図１１Ａ】
2011506022.tif 000036
【図１１Ｂ】
2011506022.tif 000037
【図１２】
2011506022.tif 000038
【図１３】
2011506022.tif 000039
【図１４】
2011506022.tif 000040
【図１５】
2011506022.tif 000041
【国際調査報告】
2011506022.tif 000042
2011506022.tif 000043
2011506022.tif 000044
2011506022.tif 000045
2011506022.tif 000046
2011506022.tif 000047