JP2011506022A - 画像処理 - Google Patents

画像処理 Download PDF

Info

Publication number
JP2011506022A
JP2011506022A JP2010538893A JP2010538893A JP2011506022A JP 2011506022 A JP2011506022 A JP 2011506022A JP 2010538893 A JP2010538893 A JP 2010538893A JP 2010538893 A JP2010538893 A JP 2010538893A JP 2011506022 A JP2011506022 A JP 2011506022A
Authority
JP
Japan
Prior art keywords
image
pressure
statistical
images
data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
JP2010538893A
Other languages
English (en)
Inventor
コリン パタキー トッド
ヒュー クロムプトン ロビン
ヤニス ゴウラーマス ジョン
Original Assignee
ユライブ エンタープライゼズ リミテッド
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by ユライブ エンタープライゼズ リミテッド filed Critical ユライブ エンタープライゼズ リミテッド
Publication of JP2011506022A publication Critical patent/JP2011506022A/ja
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/103Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/1036Measuring load distribution, e.g. podologic studies
    • A61B5/1038Measuring plantar pressure during gait
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/50Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/32Determination of transform parameters for the alignment of images, i.e. image registration using correlation-based methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/40Analysis of texture
    • G06T7/41Analysis of texture based on statistical description of texture
    • G06T7/44Analysis of texture based on statistical description of texture using image operators, e.g. filters, edge density metrics or local histograms
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20036Morphological image processing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Molecular Biology (AREA)
  • Dentistry (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pathology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Biomedical Technology (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Probability & Statistics with Applications (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Image Processing (AREA)
  • Measuring And Recording Apparatus For Diagnosis (AREA)
  • Image Analysis (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Measuring Fluid Pressure (AREA)
  • Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)

Abstract

圧力画像を処理する方法であり、該方法は、複数の第1の連続する圧力画像を受信するステップと、前記第1の連続する圧力画像の各々を処理して第2の連続する圧力画像をそれぞれ生成するステップと、前記第2の連続する圧力画像を処理することによって前記第1の圧力画像の特性を表す連続する統計画像を生成するステップと、を備える。

Description

本発明は画像処理に関し、特にはペドバログラフィカル統計分析方法に関する。
臨床的に有用なデータが得られるように画像を処理する技術は種々の医学分野で広く使用されている。種々の医学分野は人体の種々の部分を表す画像の処理と関連する。使用される画像及び画像の処理により得られる臨床情報の種類は、様々な技術が様々な医学の領域で使用されるほどである。
ペドバログラフィは、歩行の立脚期中及びポーズタスク中の被検者の足の底面に作用する圧力場を処理し比較する2次元医用画像工学である。ペドバログラフィは、手術後評価、関節炎監視、矯正学デザイン及び糖尿病性神経障害の検出及び管理を含む種々の臨床応用に使用されている。
ペドバログラフィ画像データは、典型的には、被検者の足の真下に位置する圧力センサから得られるデータを適切に変換することにより得ている。圧力センサは特定の領域内において足で与えられる圧力を表す可変出力を提供するように配置される。これらの圧力センサは発生される画像データ内のピクセルの値を提供する。圧力センサの数、それらの間隔及びそれらの大きさが画像データの解像度(分解能)を決定する。このようにして得られる画像は典型的にはサブセンチメートルの解像度である。
ペドバログラフィ画像を解析する現在の方法は定性的視覚評価又は画像を複数の個別部分で処理する種々のサブサンプリング技術の一つのいずれかを用いている。
画像の視覚評価は有効であるが、典型的には時間がかかり、極めて主観的である。加えて、有効な視覚評価は適切に訓練された人が実行できるのみである。
ペドバログラフィ画像の解析に使用される既存のサブサンプリング技術は、生ペドバログラフィデータをセグメント化し、セグメント化されたデータに解剖学的領域に従ってラベル付けし、一般的にラベル付けされた領域の自動的又は半自動的な分離を可能にしている。次にラベル付けされた領域内及び/又は領域間で統計的検定を行うことができる。
サブサンプリング技術は、視覚評価より客観的であるが、多くの欠点を被る。足の底面が連続的で分離していないものとすると、サブサンプリングはデータを破壊し、最終的には連続的な足と地面の相互作用の低解像度ビューをもたらす。サブサンプリング技術は典型的には足全体ではなく特定の領域のペドバログラフィ画像の比較と関連するのみである。
本発明の目的は、上述した問題の一つ以上を除去もしくは軽減することにある。
本発明の第1の態様によれば、圧力画像を処理する方法が提供され、該方法は、複数の第1の連続する圧力画像を受信するステップと、前記第1の連続する圧力画像の各々を処理して第2の連続する圧力画像をそれぞれ生成するステップと、前記第2の連続する圧力画像を処理することによって前記第1の圧力画像の特性を表す連続する統計画像を生成するステップとを備える。
この態様の利点は、連続する統計画像が生成される結果、圧力画像のサブサンプルされた個別の部分にではなく、連続する圧力画像に統計的解析を実行することが可能になる点にある。
前記第1の圧力画像の処理は、前記第1の圧力画像に平滑化処理を加えるステップを備えることができる。
前記平滑化処理は、受信した画像をフィルタ処理するステップと、要すればフィルタ処理された画像に所定の閾値化(2値化)処理を加えるステップと、要すれば所定の閾値化処理が加えられた画像にモルフォロジー処理を実行するステップを備えることができる。モルフォロジー処理は、例えばモルフォロジーオープニング及びクロージングを含むことができる。
前記第1の圧力画像の処理は、前記第1の圧力画像をテンプレート画像にレジストレーション(位置合わせ)するステップを備えることができる。
前記第1の圧力画像のテンプレート画像へのレジストレーションは、前記第1の圧力画像の平行移動及び回転の少なくとも一つを実行するステップを備えることができる。前記第1の圧力画像のレジストレーションは、更に、初期レジストレーションパラメータを発生させるステップと、前記初期レジストレーションパラメータを最適化するステップとを備える。前記パラメータは線形又は非線形空間変換のいずれかを定義するものとし得る。線形変換は、例えばシャーリングアフィン変換、線形空間ワーピング及び一般射影変換を加えるステップを含むことができる。非線形変換は、例えばパラメータ化された非線形空間ワーピング変換を第1の圧力画像に加えるステップを含むことができる。
空間ワーピングは、圧力画像をスケーリングするステップを備えることができ、スケーリングはそれぞれの倍率を2つの直交方向に与えることができる。
第1の圧力画像の特性は一以上の所定の基準画像との比較とすることができる。
第1の圧力画像の特性は第1の圧力画像の比較とすることができる。
複数の第1の連続する圧力画像は一つの被検者又は複数の被検者から得ることができる。
データの生成は、前記第1の圧力画像の特性を表す連続画像を生成するステップを備えることができる。生成される画像は前記第1の圧力画像とほぼ同一の解像度を有するものとし得る。例えば、特性を比較する場合、生成された画像を比較画像とすることができる。比較画像は前記又は各第1の圧力画像にほぼ等しい解像度を有することができる。このような比較画像は前記又は各第1の圧力画像の便利な比較を可能にする。
前記生成される連続する統計画像は、複数の画像の一般線形統計モデルによって生成さ得るt−マップ、F−マップ又は他の統計マップとすることができる。
本方法は、更に、前記統計画像を処理して更なる統計画像を生成するステップを備えることができる。更なる統計画像の生成は、前記統計画像に統計的推定を実行するステップを備えることができる。更なる統計画像は、複数の第1の連続する圧力画像の特性の統計的有意性を表す段階的な連続する統計推定画像とすることができる。
前記統計的推定は、t検定により生成された統計画像のt値を処理し、確率値を生成することができる。
前記統計的推定は、統計画像の平滑度及び/又は統計画像の少なくとも一つの幾何学特性を考慮することができる。
前記統計的推定は、確率場理論、ボンフェローニ補正及び/又は非パラメトリック推定に基づく技術を含むことができる。
前記第1の連続する圧力画像は足底面ペドバログラフィ画像とすることができ、第1の足底面ペドバログラフィ画像及び第2の足底面ペドバログラフィ画像の各々は足底面全体の画像とすることができる。第1の圧力画像が足底面ペドバログラフィ画像である場合、本発明は、統計的解析を既存のサブサンプリング技術の場合のように足底面の個別の区分ではなく、連続する足底面に実行する利点を有する。
本方法は、更に、各々が一以上の圧力画像の特性を表す複数の統計画像を生成するステップと、入力圧力画像を受信するステップと、入力圧力画像に基づいて前記複数の統計画像の一つを選択するステップとを備えることができる。
前記第1の連続する圧力画像は、第1の物体と第2の物体との間の圧力場の画像とすることができる。
前記第1の連続する圧力画像は、被検者とマットレスとの間の圧力場の画像とすることができる。
前記第1の連続する圧力画像は、被検者とシーティング装置との間の圧力場の画像とすることができる。
本発明の第2の態様によれば、2次元圧力画像の連続的な時間サンプルを含む動的圧力データを表示する方法が提供され、該方法は、前記圧力データの3次元体積表現を表示するステップを備え、各体積要素は2つの空間的次元と1つの時間的次元からなる。
本発明の第3の態様によれば、複数の組の2次元圧力データを含む圧力データを表示する方法が提供され、該方法は、圧力データの3次元体積表現を表示するステップを備え、前記3次元体積表現の第1及び第2の次元は空間データを表し、前記3次元体積表現の第3の次元は時間データを表す。
本方法は、更に、前記3次元体積表現をグラフィック表示するステップを備えることができる。
本方法は。更に、所定の圧力又は所定の範囲の圧力を表すデータから表面を計算するステップと、決定された表面を表示するグラフィックデータを生成するステップとを備えることができる。
前記3次元体積表現は複数の圧力画像から生成される3次元統計画像とすることができる。
本発明は、任意の都合のよい手段、例えば適切な装置を用いて実施することができる。このような装置は、上述した方法を実行するようにプログラムされたコンピュータとすることができる。本発明は、上述した方法を実行するように構成されたコンピュータプログラムも提供する。このようなコンピュータプログラムは適切なキャリア媒体上に坦持して有形キャリア媒体(例えばハードディスクまたはCD−ROM)とすることができ、また通信信号のような無形のキャリア媒体とすることができる。
本発明の一つの態様に関連して記載される特徴は本発明の他の態様のいずれにも適用可能である点に留意されたい。
本発明の一実施例が、例証として、添付図面を参照して以下に記載される。
本発明の一実施例で実行される処理のフローチャートである。 図2Aは図1の処理を用いて処理されるペドバログラフィ画像の一例である。 図2Bは図1の処理の一部分として実行されるフィルタ処理後の図2Aのペドバログラフィ画像である。 図2Cは図1の処理の一部分として実行される閾値化処理後の図2Bの画像である。 図2Dは図1の処理の一部分として実行されるモルフォロジーオープニング処理後の図2Cの画像である。 図3A及び図3Bは外転した足姿勢での歩行を表す足画像のレジストレーションを示す。 図4A及び図4Bは内転した足姿勢での歩行を表す足画像のレジストレーションを示す。 図5A及び図5Bは空間ワーピング処理によるレジストレーションを示す。 実験的に得られた一連の画像間のレジストレーション精度を示す画像である。 模範的な実験計画行列を示す。 図8A〜図8Eは複数の統計的推定法を用いて生成される統計画像を示す。 図9A〜図9Cは確率場理論を用いて生成される統計画像を示す。 図10A〜図10Cは、正常歩行、外転歩行及び内転歩行の種々の実験条件において得られた平均画像をそれぞれ示す。 図11A〜図11Bは、外転歩行と正常歩行を表す画像の統計的比較及び内転歩行と正常歩行を表す画像の統計的比較をそれぞれ示す。 ここに記載する方法を個人の識別に利用する例を示すフローチャートである。 等値面を用いるペドバログラフィの時空間体積表示である。 ペドバログラフィデータの3次元表示を生成するために実行される処理のフローチャートである。 3次元時空間統計画像の3次元表示である。
本発明の実施例は、ペドバログラフィデータを解析するプロセスを自動化するためにピクセル単位又はボクセル単位の統計画像比較を用いる。
ピクセル単位の比較はボクセル単位の比較と数学的に同一であるため、ここでは2次元ピクセル及び3次元ボクセルの両方を表すために「ピクセル」を用いる。
以下の説明のために、統計画像は、各ピクセルが一つの統計値を具現している画像に関するものである。例えば、統計画像内の特定の位置の特定のピクセル値は複数の画像の同じ位置のピクセルの統計的分布をパラメータ化したものとし得る。
ペドバログラフィデータを得るために、被検者の足の真下に位置するセンサから出力されるセンサ値を得る。例えば、各センサは、発生されるペドバログラフィデータの単一ピクセルと関連し、各センサの出力値がそれぞれのピクセルのピクセル値を決定するようにすることができる。好適実施例では、各センサから複数のセンサ値が得られ、そのうちの最大センサ値を用いてそれぞれのピクセルに対するピクセル値を決定するが、他の技術を用いることもできる。
本発明のいくつかの実施例では、圧力センサとピクセルとの間の1対1の直接マッピングを用いる代わりに、前処理を実行してセンサ値をピクセル値に関連づける。例えば、5.08mm×7.62mmの大きさのセンサが使用される場合、圧力センサ値を双線形補間又は他の補間法を用いて方形デカルト格子上に再サンプルして、例えば5.08mm×5.08mmの仮想座標にピクセル値を発生させることができる。
本発明のいくつかの実施例では、画像は第3次元が時間である3次元圧力画像にすることができる。
複数の画像を比較するのが望ましい。臨床的に有用なデータは、単一被検者から得られた複数の画像の比較(被検者内比較)によって及び異なる被検者から得られた複数の画像の比較(被検者間比較)によって得ることができる。異なる被検者からの画像をピクセル単位方式で比較する前に、各画像を、各画像が共通のサイズ、形状及び方向を有するように共通のテンプレート画像に変換する必要がある。
図1は入力画像に実行される処理を示す。ステップS1において、複数の入力画像が受信され、各画像は複数のピクセルを含み、各ピクセルは関連するピクセル値を有する。左足の入力画像の一例が図2Aに示されており、より明るいピクセルがより高い圧力の領域を示している。入力画像の処理が以下に説明される。各入力画像は同様に処理される点に留意されたい。
足の底面は連続的で従順であるので、隣接ピクセル間の圧力にシャープな変化が存在しないため、画像の平滑化は情報の有害な損失を生じることなく入力画像内の望ましくないアーチファクトを除去し、信号対雑音比を増大する。画像は当業者に周知の任意の平滑化方法を用いて平滑化することができる。
図1に示す本発明の実施例では、入力画像はフィルタリング(ステップS2)、閾値化(ステップS3)及びモルフォロジー処理(ステップS4)の組み合わせによって平滑化される。
ステップS2において、各入力画像は、例えば3ピクセル径(〜1.5cm)の対称コンボリューションカーネルでフィルタ処理され、中心ピクセルを2倍にバイアスする。図2Bは図2Aの画像に上記の対称コンボリューションカーネルを適用した結果を示す。
このフィルタ処理の副作用は、足の周辺の周りにピクセルを追加し、ひいては足の接触域が人工的に膨張することにあることが明らかである。それゆえ、図1のステップS2で入力画像がフィルタ処理されたら、ステップS3で図2Bに示される得られた画像に閾値化を提供して図2Cに示される画像を発生させる。適用する閾値は足の周囲のピクセルを除去する。閾値は、所定の閾圧力より低い圧力値に対応するピクセルを除去するように設定され、閾圧力はもとの足接触域の周囲のピクセルのピクセル値を示すように選択される。
図2Cを参照すると、ステップ3における閾値化の適用は外側指骨の低圧力領域に突出ピクセル1を生じる。孤立した突出ピクセルは統計的解析に影響を与えやすく、これらはテンプレート画像からモルフォロジーオープニングにより除去すべきである。このようなオープニングはステップS4で実行され、図2Dに示される画像を生成する。
以上の記載は入力画像内のアーチファクトを除去する入力画像の平滑化と関連する。入力画像の平滑化は信号対雑音比を高め、図1を参照して以下に更に詳細に記載される方法の順次のステップを容易にすることができる。しかし、平滑化はここに記載する方法の必須の要件ではない点に留意されたい。
画像を適切に比較することを可能にするために、入力画像はステップS5において、入力画像が共通のテンプレートに整列するように、共通テンプレートにレジストレーション(位置合わせ)される。更に、被検者内ペドバログラフィデータが比較的短い期間(被検者の足の寸法が変化しない期間)の間に得られる場合には、共通テンプレートへの整列はペドバログラフィ画像の更なる変換の必要なしに被検者内解析を可能にする効果がある。本発明の実施例で使用できる再整列の一つの方法が以下に記載される。
レジストレーションプロセスは、所定の位置への入力画像の平行移動、画像の回転、回転された画像の最終平行移動前に入力画像を初期位置に戻すための更なる平行移動の組み合わせを用いる。
ベクトルqの要素として3つの一般化パラメータを取ると、変換シーケンスT(q)は以下のように定義することができる。
Figure 2011506022
ここで、
(xc, yc)は足の重心の座標であり、
q1及びq2はそれぞれ画像の水平方向及び垂直方向の平行移動を示し、
q3は画像の回転を示す。
変換シーケンス(5)は、画像を行列(2)で示されるように原点に平行移動させ、画像を行列(3)で示されるように重心を中心に回転させ、回転された画像を行列(2)の逆行列で示される元の位置に平行移動させて戻す。次に、回転された画像は行列(4)で示されるように任意の位置に平行移動される。
変換シーケンス(5)は画像のピクセルの再配置と関連する。従って、画像がベクトルpで表される場合、適切なピクセル値p’を有する変換された画像は方程式(6)で与えられる。
p’= p(XT) (6)
ここで、
Xは同次型のJピクセル座標の(3xJ)行列である(即ち行列乗算を可能にするために1が付加されたダミー行を有する)
図3Aは、外転した足姿勢で歩く被検者から得られる画像2を示す。画像2を整列させるべきテンプレート画像3も示されている。図3Aは画像2の長主軸4及び短主軸5を示す。角度θは基準軸6に対する画像2の短主軸5の方向を示す。図3Bは画像2のテンプレート画像3へのレジストレーションの結果を示す。
図4A及び4Bは図3A及び3Bに対応するが、内転した足姿勢で歩く被検者から得られる足の画像に関連する。
ソース画像p1と、該ソース画像をレジストレーションすべきテンプレート画像p0が与えられると、最適な変換を提供するベクトルqが決定される。足の幾何形状を利用して、qの初期推定q(0)を方程式(7)に従って行うことができる。
Figure 2011506022
ここで、
(xc0, yc0)はテンプレート画像p0内の足の重心であり、
(xc1, yc1)はソース画像p1内の足の重心であり、
θ0は基準軸に対するテンプレート画像p0内の足の短主軸の方向であり、
θ1は基準軸に対するソース画像p1内の足の短主軸の方向である。
主軸は、圧力画像の「慣性」行列の固有ベクトルとして計算することができ(即ち、ピクセル値が局所的「質量」を表す場合)、また2値画像から計算することができる(即ち、すべてのピクセルが等しい質量を有する場合)。短主軸はその軸を中心とする慣性モーメントが最小になる軸である。
q(0)は方程式(5)で表される変換シーケンスを定義するために使用でき、変換シーケンスを方程式(6)で示されるようにソース画像p1に適用して変換されたソース画像p1’を提供することができる。このとき、テンプレート画像p0と変換されたソース画像p1’を、相違関数εを用いて直接比較することができる。
Figure 2011506022
ここで、εは2つの画像の要素間の二乗誤差の和であり、上付き文字「T」はベクトル転置を示す。
そして、レジストレーションの目標は形式的に、
Figure 2011506022
として記載することができる。
最適化問題(9)は3Dパラメータ空間Ωにおいて無限で非拘束である。最適化は任意の最適化技術を用いて実行することができる。方程式(8)の目的関数は方程式(7)で与えられる初期推定q(0)の近くで凸関数である。従って、最適化は、例えば準ニュートン最急降下勾配探索法を用いて実行することができる。
画像レジストレーションはコンピュータ要求技術であるが、上記のように足幾何形状を利用すると、良好な初期近似を生成することができ、比較的簡単な最適化技術で満足なレジストレーションを達成することができる。満足な収束を小さなデータセットで与えるために簡単な勾配探索法が見出されている。とは言うものの、二乗誤差の和を用いる勾配探索法は広範囲の被検者及び/又は画一的な圧力プロファイルに対して有効でないかもしれない。なぜなら、方程式(8)の目的関数は初期近似の周囲の小領域においてのみ凸関数であることが確かめられているためである。スパース勾配最適化及び2値形状レジストレーションなどの他のアプローチの方が一般化されたレジストレーションアプローチに適切であり得る。
共通テンプレートに整列された一組の画像に対して、被検者内画像にピクセル単位の統計的解析を行うことができるが、被検者間の足の形状の差は、被検者間解析を行う前に、図1のステップS5で実行されるレジストレーションの一部として画像を空間的に正規化するのが望ましいことを意味している。一つの被検者からのソース画像を別の被検者から任意に選択されたテンプレート画像にレジストレーションすることにより空間ワーピングを達成することができる。主軸に沿って測定された足画像の幅及び長さを正規化するために5パラメータ非せん断(non-shearing)アフィン変換を用いることができる。このような変換を適用するために、方程式(5)で与えられる変換シーケンスは、
Figure 2011506022
に変更される。
ここで、
Figure 2011506022
及びU0,R及びU1は式(2),(3)及び(4)で与えられる。
ここで、ベクトルqは5つの要素を有し、最初の3つの要素は上述したとおりであり、2つの追加の要素q4及びq5は足画像の幅及び長さに対するそれぞれのスケーリングファクタ(倍率)を示す。方程式(10)の変換シーケンスは線形空間ワーピングの一形式を表し、足の幅と長さを等しくない量(q4≠q5の場合)でスケーリングすることができる。
行列(11)で与えられるRvは足を垂直方向まで回転させるため、行列(12) Sは幅をファクタq4倍に、長さをファクタq5倍にスケーリングすることができる。最適化を進めることができる良好な開始値は、
Figure 2011506022
で与えられる。
ここで、
l0及びl1は長主軸5に沿って測定されたテンプレート画像p0及びソース画像p1の足の長さであり、
wo及びw1は短主軸4に沿って測定されたテンプレート画像p0及びソース画像p1の足の幅である。
q(0)を方程式(10)に使用してソース画像p1を方程式(6)に従って変換し、変換されたソース画像p1'を与えることができる。この場合、テンプレート画像p0と変換されたソース画像p1'を方程式(8)に従って直接比較することができるため、ベクトルqの値を方程式(9)に従って最適化することができる。
図5Aは、方程式(10)を用いる変換前に98kgの男性被検者から得られる足画像7を示す。図5Bは、方程式(10)及び(9)を用いた最適空間ワーピング後の同じ画像を示す。ここでは画像7は47kgの女性被検者から得られたテンプレート画像8にレジストレーションされている。図5A及び図5Bは、様々な足幾何形状を有する被検者からの画像をレジストレーションにより共通のテンプレート空間に変換できることを示している。
方程式(8)で示される二乗誤差の和を用いてベクトルqの値を最適化することができるが、もっと直感的な相違行列を用いて報告目的のレジストレーション性能指示を得ることができる。2つの画像がオーバラップしない程度は集合表記法を用いて次のように書き表せる。
Figure 2011506022
ここで、
p0及びp1'は上述した画像であり、
|pi|は集合pi内の非ゼロ要素の数であり、
Figure 2011506022
方程式(14)の分母は画像p0及び画像p1'の両方が占めるピクセルの総数を示す。方程式(14)の分子は一方又は他方の画像で覆われ、両方の画像で覆われないピクセルの数を示す。完全なレジストレーションはεの値が0%のときに与えられる。これは、方程式(14)の分子がゼロ値を有するとき生じる。ゼロオーバラップ(即ちεが100%に等しい)は2つの集合間の共通集合がヌル集合であるときに生じる。方程式(14)の分子は2値画像であり、レジストレーション品質を確かめるために有用である。
Figure 2011506022
表1は、ここに記載する方法の種々のステップの計算持続時間を示す。平滑化及び再サンプリングは無視できるほど急速に行われることがわかる。画像レジストレーションは1画像につき1秒程度を要した。総時間の大部分がレジストレーションに費やされる。被検者間レジストレーションの方が最適化ステップのサイズ制限のために被検者内レジストレーションより速かった。
レジストレーション精度は10%程度であり、平均で、テンプレートピクセルとソースピクセルの約10%がオーバラップしなかった。非オーバラップピクセルは、被検者内レジストレーションに対してはもっぱら足の周囲に見られた(以下に更に詳しく記載される図6参照)。被検者間レジストレーションに対しては、非オーバラップピクセルは中央アーチ部及び指骨部にも見られた。
図6の画像は、レジストレーションされた画像の一群に排他的論理和(XOR)演算を加えることによって生成される。図6の画像は非オーバラップピクセルの発生頻度を示し、白いピクセルほどテンプレート画像とソース画像の共通部分に頻繁に含まれなかったことになる。
約90%のレジストレーションのオーバラップ比はペドバログラフィ分解能に起因する非オーバラップピクセルの空間分布に起因し、解剖学的要因にも起因する。ペドバログラフィ分解能は5mm程度であるため、足の境界は5mmの精度でトレースできるのみであり、厳重な制御対象物でも1ピクセルの厚さのXOR境界画像を生じる。これは、試行間のバリアビリティと相まって、不完全なレジストレーションオーバラップを説明する。被検者間の解剖学的差も同様に正規化後の中央アーチ部及びつま先部に観測されるバリアビリティを説明する(図6参照)。
上述したタイプのアフィン変換は全体的な足形状にのみ影響し、一般に完全な被検者間レジストレーションをもたらさない。性能を向上させるために非線形空間ワーピングを使用することができる。比較的高次で低周波数の100パラメータディスクリートコサイン変換非線形ワーピングアルゴリズム及びロバストな進化的最適化アプローチを用いることによって、オーバラップを数パーセント高めることができる。尚、上述したレジストレーション方法の代わりに任意の適切なレジストレーション方法を使用してもよいことは当業者に明らかであると点に留意されたい。
上述のレジストレーションプロシージャの目的は、図1のステップS6〜S9で示されるように、統計的検定を直接画像にピクセルレベルで行うことを可能にすることにある。一連の統計的検定を使用できるが、選択される実際の検定は用途により決まる。
ステップS6において、一般線形モデルが定義される。定義された一般線形モデルは、各画像が関連する実験状態を記述する「計画」行列という数値行列により規定される。例えば、第1の状態(例えば手術前)からの5つの画像及び第2の状態(例えば手術後)からの5つの画像がある場合、計画行列は10×2行列とし、2つの列がそれぞれ手術前及び手術後を表す。手術前列は最初の5行に1の値を有し、手術後列は最後の5行に1の値を有し、他の値はすべて0である。このようなモデルは2つの実験状態間を区別し、周知の2標本t検定(two-sample t-test)を記述する。この簡単なモデルに加えて、一般線形モデルは、例えば1標本t検定(one-sample t-test)、分散の解析(ANOVA(analysis of variance))(F tests)、分散の多変量解析(MANOVA(multivariate analysis of variance))、共分散の解析(ANCOVA(analysis of covariance))、共分散の多変量解析(MANCOVA(multivariate analysis of covariance))、固定効果解析及び混合効果解析の使用を含む任意の線形実験計画を可とする。
図7は計画行列のグラフ表示である。列は実験的要因を表し、行は実験的反復を表す。黒領域は0の値を有し、白領域は1の値を有する。図7の計画行列は、9つの被検者と、正常歩行、外転歩行及び内転歩行の3つの歩行状態を含む実験を示す。各状態の10回の連続試行がすべての被検者により実行され、垂直軸に沿って示されるように、全部で270の画像が生じる。これは固定効果モデルの一例であり、この例では各被検者効果が行列の右側の被検者ブロックを通して固定効果として処理される。これは非無作為実験計画法の一例でもあり、この例では実験状態(最初の3列)は各被検者の30試行の中で無作為に分散される試行ではなく単一ブロックの反復を含む。一般線形モデル化アプローチを使用するためにこのような計画行列をステップS6で指定しなければならない。
一般線形モデルは、例えば時間ドリフトなどの実験的関心因子でない因子をモデル化するために使用できる。これらの因子は迷惑因子と呼ばれる。図7の計画行列の右側に記された被検者因子はヌイサンス因子である。この例では、被検者間の差(例えば体重など)は関心因子でなく、正常歩行、外転歩行及び内転歩行と関連する結果のみが関心因子である。
一般線形モデルは、線形統計モデル化に対して最もフレキシブルなアプローチであるとここでは記載しているが、ステップS6において他の統計モデルを使用することもできる。
処理はステップS6からステップS7に移る。ステップS7において、実験状態を統計画像にマップするパラメータが推定される。例えば、2標本t検定の場合には、これらのパラメータは2つの実験状態の平均圧力に対応する。この場合には、統計画像内の各ピクセルにつき2つのパラメータが存在するのに加えて各ピクセルにつき追加の分散パラメータが存在する。これらのパラメータは、以下に一例につき更に詳細に記載されるように、最小二乗法を用いて計算される。
処理はステップS7からS8へ移る。ステップS8において、統計画像がパラメータから生成される(一例につき以下に記載される)。統計画像において、各ピクセルは単一の統計値を表す。一般線形モデルは種々の統計量、例えばt統計量(いわゆるtマップを生成する)、F統計量(いわゆるFマップを生成する)及びχ2統計量(いわゆるχ2マップを生成する)の計算を実行する。
次に処理はステップS9に移り、統計的推定処理を用いて統計画像内の統計値及び/又は空間プロセスの統計的有意性を推定する。図8A及び図9Aの各々は各ピクセルが関連するt値を有する任意の統計画像を示す。多数のピクセルが存在し、多数のt値が存在するため、t値自体が図8A及び図9Aに示す画像を形成する。
単一t値は統計的有意性を搬送しない。同様に、各ピクセルがt値を表す統計画像は統計的有意性を搬送しない。統計的有意性を計算するために、統計的推定処理を用いて統計画像を構成するt値を関連する確率(‘p’)値に変換し、それによって図1のステップS11に示すように、段階的な統計的推定画像を生成する。統計的推定画像及び対応するp値の解釈は使用する統計的推定処理のタイプに依存する。例えば、一つの解釈では、1%のp値(p=0.01)は、観測された実験的分散に基づいて、実験的治療が観測効果を発生した可能性は1%であることを意味する。1%のp値の別の解釈は、確率過程(ランダム過程)によりデータが発生された場合、この大きさのt値は多数の実験的反復の僅か1%に観測されることが予測されることを意味する。それゆえ、p値が適切に小さいとき、この実験効果は確率過程に起因しなかったことが合理的に確かである。統計的有意性の典型的な臨界カットオフ(‘α’と呼ばれることもある)はp=0.05である。実験者設定α値は、帰無仮説が実際に真であるとき帰無仮説を棄却してしまう「偽陽性(false-positive)」としても知られている第1種過誤に対する保護レベルを指定する。
図8Bは、p=0.05より小さいp値を有する図8Aのピクセルを強調する2値画像を示し、p閾値は補正されてないものとする。各ピクセルに対するp値はそれぞれのピクセルに対するt値に基づいて個別に計算されている。閾値処理は各ピクセルの個別処理のために多くのピクセルを持ち続ける。
単一のピクセルt値に対するこのようなp値の計算は単一の周知の統計方程式による簡単な1対1マッピングである。しかし、この処理は、多くの統計的検定が行われた事実も、画像データの空間特性も、隣接空間内のピクセルが相関動作を有することも無視するため、統計的に無効である。従って、統計画像に統計的推定を行う有効な方法は、上述した周知のp値計算より複雑である技術を必要とする。p値を生成するために種々の方法を使用できるが、ここでは特に本発明の実施例に適用可能な3つのクラスの統計的推定:(1)ボンフェローニの補正、(2)確率場(ランダム場)理論及び(3)非パラメトリック推定:を記載する。
これらの処理について記載する前に、これらの推定処理のすべての一般的な特性は多数の比較の問題の原因となることを理解することが重要である。(図8A及び図9Aに示すタイプの)統計画像の計算は、各ピクセルにつき1つの多数の統計的検定を含む。0.05の無補正のα値の選択は1つのピクセルに対する偶発過程を保護するのみであるため、無補正のp閾値は一般に統計画像プロセスの有意性の過大評価を生じる。それゆえ有効なαを維持するために、多重比較を考慮しなければならない。
ボンフェローニ補正は多重比較を補正する最も簡単な方法である。モンフェローニ補正は、足表面全体に亘ってファミリー誤差率が所望の値のα(上述の例では5%)に維持されるようにp閾値を減少させる。ボンフェローニ補正は、
critical=1−(1−α)1/k (15)
として計算される。
ここで、Kは統計的検定の数又は、等価的に、統計画像内の非ゼロピクセルの数である。上述したように、図8Bは0.05の無補正の閾値p以下に下がるのに十分なほど大きなt値(この場合にはt=1.73)を有する図8A内のピクセルを示す。図8Cはt=4.42のモンフェローニ補正閾値しきい値より大きいt値を有するピクセルを示す。ボンフェローニ補正後に有意になるピクセルは無補正の推定処理に対してよりかなり少なくなることがわかる。
無補正及びボンフェローニの両推定処理は、図8Aに示す生統計画像のt値を超閾値クラスタ(αから計算される臨界t閾値を越えるt値を有するピクセルのクラスタ)の群内に含む段階的な連続統計画像を生成することに注目すべきである。
図8B及び図8Cに示す統計的推定画像は統計的クラスタ間にギャップがあるという意味で段階的である点に留意されたい。ギャップは画像の不連続性(即ち、t値の鋭い低下)を示すが、t値は超閾値クラスタ内では連続的であり(即ちt値のピクセル単位の鋭い変化はなく)、従って「段階的−連続的」画像という。ボンフェローニ補正の場合には、実験者設定α=0.05であるものとすると、図8Bに示す統計的推定画像の解釈は、「各超閾値クラスタ内のピクセルが偶然に発生する確率は5%である」である。
上述のモンフェローニ補正の主な弱点は、データの空間特性を適切に考慮しない点にある。500ピクセルの集団を処理する場合、ボンフェローニ補正はこれらのピクセルの空間的関係に関係なく同一であり、ピクセルは空間内にランダムに無限に散乱される。しかし、足圧力画像及び他の圧力画像は明らかにランダムな空間プロセスでなく、これらの画像は一貫した足形状の形態の空間秩序を有する。
確率場理論(RFT)は、(i)多重比較に対してよりリアルで低保存性の補正を提供すること又は(ii)統計的有意性を個々のピクセルではなく超閾値クラスタに割り当てることのいずれかを活用するものである。これらのRFT処理はいずれも元の画像の2つの特性、画像平滑度及び空間幾何形状、に依存する。
画像平滑度は一般線形モデルパラメータの推定中(図1のステップS7)に計算される残存画像から推定される。各残存画像Ei(その計算は方程式(35)に関して以下に記載される)を最初にピクセル標準誤差(s)(その計算は方程式(34)に関して以下に記載される)により正規化して、正規化された残存画像Zi
i=Ei/s (16)
を生成する。
次に、残存空間導関数の共分散行列∧を、
Figure 2011506022
を用いて計算することができる。
ここで、kは非ゼロピクセルのインデックスであり、x及びyは空間方向を示し、varは分散の測定値を与える関数であり、covは共分散の測定値を与える関数である。最後に、全体的な平滑度はFWHM:
Figure 2011506022
により推定される。
ここで、FWHMは、確率場データとともに畳み込まれると、観測される残余と同一の平滑度を生成する、ガウスカーネルの「半値幅(full width at half maximum)」を表す。従って、FWHMは画像平滑度の直接推定を提供する。
RFT推定処理は、ピクセルの接続性によりアルゴリズム的に決定し得る空間幾何形状にも依存する。2次元画像に対しては、RFT関連空間幾何形状は「分解素子カウント(resolution element counts)」又は「リセルカウント(resel counts)」という3つのパラメータR0,R1及びR2により要約することができる。
0=K−(Ex+Ey)+F (19)
1=(Ex+Ey−2F)/FWHM (20)
2=F/FWHM2 (21)
ここで、Ex及びEyはそれぞれx方向及びy方向に隣接するピクセルの対の数であり、Fは探索区域内の4ピクセル方形の数である。従って、Rパラメータは探索空間を次元数(本例では2)まで規定し、面積(R2)、周囲(R1)及びオイラー特性(R0)を測定する。これらの「リセルカウント」パラメータは平滑度を考慮した画像データセット内の独立した観測の有効数とみなすことができる。従って、場が平滑になるほど(FWHMが大きくなるほど)、独立した観測の数が少なくなる。
(FWHMによる)画像平滑度及び(リセルカウントによる)画像幾何形状の計算後に、RFTベースの推定を進めることができる。上述したように、この推定は、2つの方法のひとつ、即ち(i)多重比較を考慮するためにRFT補正閾値を計算することによって、又は(ii)任意の実験者設定t閾値に続いて超閾値クラスの統計的有意性を計算することによって、生じさせることができる。
前者のRFT処理はボンフェローニ補正に類似し、臨界t閾値hを計算し、閾値hを超えるt値を有するピクセルを統計的に有意とみなす。閾値hは次のように計算される。
Figure 2011506022
ここで、RDはリセルカウントであり、ρDはオイラー特性密度であり、νは統計モデルの自由度の数であり(その計算は下記の方程式(36)で示される)、Γはガンマ関数である。ここで、P(tmax > h)は統計画像内の最大t値が閾値hを越える理論的確率である。
実験者設定αが与えられれば、パラメータtを方程式(22)の右辺がαに等しくなるように繰り返し変化させることによって臨界閾値hが見つかる。この処理は、ボンフェローニ補正と比較して全体的に低い保存性の閾値を生じ、図8D及び図9Bに示されるような僅かに広い超閾値設定を生じる。α=0.05であるとすると、図8Dに示される統計的有意性は次の通りであり、観測画像の平滑度及び探索区域の空間幾何形状が与えられれば、これらの超閾値ピクセルが偶然生起しうる確率は5%である。
上述したRFT推定処理の主な弱点は、統計的有意性がピクセル単位で割り当てられ、そのプロセスは超閾値ピクセルの空間クラスタリング(例えば図8Dに示されている)を活用していない。第2のRFTベースの推定処理は任意の実験者設定閾値hを用い、次に所与の空間広がりで特定の超閾値クラスタが偶然発生する確率を計算する。最初に、エクスカーション集合ピクセル(N)及びクラスタ(m)の予測数を、
E{N}=Aρ0(h) (26)
E{m} =Aρ2(h) (27)
として計算する。
ここで、AはFWHM2により正規化された探索区域である。最大の超閾値クラスタがn以上のピクセルを含む確率[P(nmax>n)]は、
P(nmax>n)=1−exp(−E{m}e-Bn) (28)
B=Γ(2) E{m}/ E{N} (29)
である。
次に、これらの確率値を各超閾値クラスタについてそれらのサイズkに基づいて計算し、図8E及び図9Cに示されるように、これらの確率値を統計的推定画像上にクラスタに隣接して表示することができる。α=0.05であるものとすると、統計的有意性の解釈は次の通りであり、「観測された画像平滑度、空間幾何形状、及び閾値hが与えられれば、同じサイズの超閾値統計クラスタが偶然発生する確率はp%である」と解釈される。ここで、pは特定のクラスタの計算されたp値である。
第2のRFTベースの推定処理に対する方程式はガウス場に対して導出されたものであるため、t場に対しては高い自由度で近似値を与えるのみである。t場及び他の統計分野に対するもっと正確な結果は、統計学の文献、例えば「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」に与えられている。しかし、上記の方程式は高い自由度の場合に適切である。
最後のクラスの推定処理は非パラメトリック推定である。上述したRFT処理の潜在的な弱点は、画像ランダム性の基礎的予想を記述するために統計パラメータに依存している点にある。具体的にいうと、RFTは統計分布を決定するためにピクセルレベルの「平均」及び「標準偏差」パラメータを用いている。このようなパラメトリックアプローチは強力であるが、種々の仮定、例えば残部のノーマル分布(方程式35に含まれる)に最も臨界的に依存する。足圧力画像が幾何学的不一致の区域を含むとき(例えば、高いアーチ構造の被検者と低いアーチ構造の被検者を比較するとき)、高いアーチ構造の被検者に対して地面と接触しない区域において多くのゼロ観測が存在するためにこの仮定は破られる。
この種の状態に対しては、非パラメトリック(NP)推定がより適切である。実験的統計画像(E画像)(それらの例は図8A,9Aに示されている)から出発して、NP処理はモンテカルロシミュレーションを用い、上述した一般線形モデルの生成中に供給される実験的パラメータの順序をランダムに入れ替えることによって一群の仮説統計画像(H画像)を計算する。E画像及び各H画像に対して、統計値、例えば最大の超閾値統計クラスタが計算される。これによりH画像及びE画像は計算された統計値の実験的確率分布を決定する。この場合、E画像の有意性は、その統計極値がランダムに生成されるH画像の値に匹敵するかどうかを決定することによって、この実験的確率分布から直接決定される。α=0.05であるものとすると、NPベースの推定画像は次の通りに解釈され、即ち閾値hが与えられれば、同じサイズの超閾値統計クラスタが偶然発生する確率はp%である、と解釈される。
NP処理に対しては、基礎的分散場を随意に平滑化してより滑らかな統計画像(図8Aに示される画像とともに発生される)を生成することができ、過度に保存的なピクセル分散の推定の場合にはより大きな超閾値クラスタを生成しうるプロセスはこれらのクラスタからいくつかのピクセルを排除する。換言すると、NP推定処理はパラメータ化された分散に明確に依存しないため、推定処理の有効性に影響を与えることなく分散を平滑化することができる。
次に上述した技術を用いて実行した実験について記載する。左足のペドバログラフィ記録を、0.5メートルのフットスキャン3Dシステム(RSスキャン社(ベルギー)から市販されている)を用いて、9人の健康な被検者(5人の男性、4人の女性、年齢29.8±7,7才)から収集した。3つの実験状態、即ち正常歩行足、外転歩行足及び内転歩行足、を調査した。他の歩行パラメータ(例えば歩行速度)は精密には制御しなかった。被検者は各状態の試行を連続的に10回行い、全部で270(各被検者につき30)のペドバログラフィ記録を生じさせた。各試行につき立脚期時系列から最大センサ値を抽出して一つの画像を形成した。任意のサイズ、形状及び向きの足画像を、上述した平滑化、再位置調整及び正規化技術を用いて変換した。統計的検定が実行された画像は53×22の格子内に含まれ、そのうちの742ピクセルが足に含まれる。上述したタイプの被検者間及び被検者内検定を実施した。
Figure 2011506022
外転及び内転歩行における足の向き(主軸の向きで測定される)は、表2に示されているように、正常方向から約30度相違している。小さい角度値のために非円形統計量であるものと仮定すると、ワンウェイANOVAは実験的処理の効果を確認した。
図1のステップS6に関連して記載したように、一般線形統計モデル:
ik=gijβjk+hilγlk+eik (30)
が生成され、ここでpikは圧力記録i内のk番のピクセルにおける圧力観測値である。ここで、iはI個の圧力記録のインデックスを示し、kは足を含むK個のピクセルのインデックスを示し、jはJ個の実験因子のインデックスを示し、lはL(本例では9)人の被検者のインデックスを示す。誤差(eik)は、実験状態間及び被検者間で互いに独立で同等であり且つ正規分布するものと考えられるが、ピクセル間ではそうではない。説明変数gij及びhilはそれぞれ処理効果及びヌイサンス因子(被検者効果)に対応する。βjk及びγlkは未知のパラメータである。
図1のステップS7に関連して記載したように、画像を実験状態間で比較するためには、最初に方程式(30)の未知のパラメータを解く必要がある。方程式(30)のモデルは簡潔に、
P=Gβ+e (31)
と書き表せる。
ここで、Pは平均補正された圧力データ行列であり、Gは図7に示される2値実験的計画行列であり、βはパラメータ行列(ヌイサンス因子を含む)、及びeは誤差行列である。データPはピクセル内で平均補正されるため、モデルは余分の定数項を必要としない。一般に、1の定数列をGに加えることにより平均補正を不要にできるので、データPを平均補正する必要はない。未知のパラメータの値を計算するために、βの最小二乗推定(bで示す)を、
Figure 2011506022
によって得る。ここで、GはGの擬似反転である。
図1のステップS8において、各ピクセルについてt統計値を次のように計算する。
Figure 2011506022
ここで、cは(l*(J+L))コントラストベクトルであり、bkはbのk番の列であり、skは各ピクセルの標準誤差推定であり、Ekは方形残余Eのk番の対角要素である。外転及び内転圧力が正常圧力を超える場合をそれぞれ示す2つのコントラストベクトル、[-1 1 0 0]及び[-1 01 0]、を試験した。ここで、は8つの被検者ヌイサンス因子に対応する(1×8)ヌルベクトルである。上述したように、t値自体が「tマップ」という画像を形成する。この統計画像は元のペドバログラフィ画像と同じ解像度を有し、元の足の2D形状を有する。しかし、ピクセル値が圧力値のみを表す元のペドバログラフィ画像と異なり、tマップのピクセルは実験的因子の一次結合の効果を表す。
図10A,10B及び10Cはそれぞれ正常歩行、外転歩行及び内転歩行の実験状態において得られる平均画像を示す。図10Bは、外転歩行においては、図10Aに示される正常歩行に比較して、第一中足骨頭及び第一趾の下部内側に増大した圧力を示す。図10Cの画像データ表される内転歩行は外側中足骨頭の下部に増大した圧力を有するとともに外側指骨の下部及び縦アーチの下部に大きな接触面積を有する。これらの傾向の統計的有意性は図11A及び図11Bに示されるtマップにより確認される。図11Aから、外転歩行においては足内側部分に増大した圧力が存在することがわかるが、(図11Bに示されるように)内転歩行においては足外側部分に増大した圧力が存在する。ノイズクラスタのサイズは小さいので、RFTベースの推定を行えば、周辺ノイズを除去できる。図11A及び11Bに示される画像及び図8E及び9Cに示される画像はペドバログラフィ画像を比較する便利なメカニズムを提供する。
ペドバログラフィデータの臨床/研究解析を実行するのに加えて、上述したペドバログラフィ画像解析方法はもっと一般的な応用を有する。例えば、本発明の一つの応用として、例えば空港で使用する個人識別方法がある。個人のデータベースを維持し、該データベースに、各個人について、各個人から取り出した複数のペドバログラフィ画像から上述の方法により生成された複数の処理画像及び統計画像を格納する。個人を識別するために、更にペドバログラフィ画像を取り出し、データベースに格納された複数の処理画像及び統計画像と比較して一致するかどうか決定する。
図12は個人を識別するように実行できる本発明によるプロセスを示すフローチャートであり、各個人の詳細を該個人の足のペドバログラフィ画像から取り出された統計画像と一緒に格納したデータベースが存在するものとする。
図12を参照すると、ステップS12において個人の足の少なくとも一つのペドバログラフィ画像が得られる。圧力画像は、任意の適切な方法を用いて、例えば適切なセンサを含む圧力板の上を個人に歩いてもらうことにより得ることができる。少なくとも一つの適切な画像が得られたら、処理はステップS13に移り、このステップにおいて、得られた画像を、データベース内のすべての被検者画像がレジストレーションされているテンプレートにレジストレーションする。画像のレジストレーションは上述したように行うことができる。図1に関連して記載したように、レジストレーション前に、得られた画像を平滑化して得られた画像からアーチファクトを除去することができる点に留意されたい。ステップS13から、処理はステップS14に移り、このステップにおいて、レジストレーションされた画像をデータベース内の各画像及び/又は統計画像及び必要に応じ各処理画像と比較して最も一致するものを決定する。得られた画像とデータベース内の画像との比較は、パターン認識技術、例えばk近接アルゴリズム(k-NN)を用いて実行することができる。あるいは、上述した統計手法を使用して一致の信頼性を決定することもできる。更に精密な反復一致を実行して探索結果を更に狭くすることができる点に留意されたい。
図12について記載したプロセスを、24人の健康な若い被検者の標本サイズについて、各人につき10回のデータベース繰り返しで実行し、100%の識別精度を得た。
別の応用としては、被検者とマットレス又は車椅子との間の圧力場を解析して被検者の動きを検出することができる。これは、本発明のリアルタイム実施で達成でき、この場合には圧力場の有意変化が被検者の動きを示し、患者が補助を必要としているかどうかを決定するのに役立つ。
本発明は、更に、車両設計の分野で、例えば車のドアと車のフレームとの界面又はワイパー表面とフロントガラス面との界面の解析に使用することができる。一般に、本発明は、圧力場を解析する能力が必要とされる場合に適用可能である点に留意されたい。
ここに記載のペドバログラフィデータの処理技術は3次元画像にも有効に適用できる。特に、3次元の時空間体積のデータの視覚化を可能にする技術が以下に記載され、その一例が図13に示されている。図14は図13の視覚化を生じさせるために実行される処理を示す。
図14を参照すると、ステップS15において、3次元の時空間圧力画像を受信する。このような3次元データは時間的に連続してサンプリングする任意のペドバログラフィハードウェアシステムにより生成され、典型的な市販のシステムは2D画像を500Hzで記録する。例えば、人間の歩行中の0.6秒の典型的な立脚期は約300の2D画像を発生する。これらの2D画像を積み重ねて単一の3次元時空間体積を形成することができる。
ステップS16において、受信した3次元画像を3次元ガウスカーネルで畳み込むことによって平滑化する。ステップS17において、平滑化された画像をアルゴリズム的に補間することによって各々が特定の圧力と関連する複数の等値面を計算する。3つの等値面が図13に示されている。第1の等値面10は1Ncm-2の圧力を有する区域を表し、第2の等値面11は5Ncm-2の圧力を有する区域を表し、第3の等値面12は10Ncm-2の圧力を有する区域を表す。図14のステップS18において、各当表面に関連特性(例えば色、透明度及びテクスチャ)を付与する。ステップS19において、明暗、カメラ及び他の環境パラメータを、所望の解釈可能特性を有する視覚化が提供されるように設定する。次に、ステップS20において、この視覚化がレンダリングされ、ユーザに適切な形で表示される。図13は単一の歩行試行からの圧力値の時空間体積の視覚化を示すが、同じ視覚化処理を上述した方法を用いて生成された3D統計画像に適用することができ、例えばその一例が図15に示されている。
特定の技術及び変換が例証として記載されているのみであり、これらの技術の順序は変更可能であり、他の適切な技術も同様に使用できることを理解されたい。例えば、統計解析技術は例証としてのみ示されている。実際上、本発明の範囲から外れることなく多くの変更を記載の実施例に対して加えることができる。

Claims (41)

  1. 圧力画像を処理する方法であって、該方法は、
    複数の第1の連続する圧力画像を受信するステップと、
    前記第1の連続する圧力画像の各々を処理して第2の連続する圧力画像をそれぞれ生成するステップと、
    前記第2の連続する圧力画像を処理することによって前記第1の圧力画像の特性を表す連続する統計画像を生成するステップと、
    を備える、圧力画像の処理方法。
  2. 前記第1の連続する圧力画像を処理するステップは、前記第1の圧力画像に平滑化処理を加えるステップを備える、請求項1記載の方法。
  3. 前記第1の圧力画像に平滑化処理を加えるステップは、受信した画像をフィルタ処理するステップを備える、請求項2記載の方法。
  4. 前記第1の圧力画像に平滑化処理を加えるステップは、更に、前記フィルタ処理された画像に所定の閾値化処理を加えるステップを備える、請求項3記載の方法。
  5. 前記第1の圧力画像に平滑化処理を加えるステップは、更に、前記所定の閾値化処理が加えられた画像にモルフォロジー処理を実行するステップを備える、請求項4記載の方法。
  6. 前記第1の圧力画像を処理するステップは、前記第1の圧力画像をテンプレート画像にレジストレーション(位置合わせ)するステップを備える、請求項1−5のいずれかに記載の方法。
  7. 前記第1の圧力画像をテンプレート画像にレジストレーションするステップは、前記第1の圧力画像の平行移動及び回転の少なくとも一つを実行するステップを備える、請求項6記載の方法。
  8. 前記第1の圧力画像をテンプレート画像にレジストレーションするステップは、
    初期画像変換パラメータを発生させるステップと、
    前記初期画像変換パラメータを最適化するステップと、
    を備える、請求項6又は7記載の方法。
  9. 前記第1の圧力画像を処理するステップは、更に、前記第1の圧力画像をテンプレート画像に対して空間的にワーピングするステップを備える、請求項6,7又は8記載の方法。
  10. 前記空間的ワーピングステップは、前記第1の圧力画像に線形又は非線形の空間変換を加えるステップを備える、請求項9記載の方法。
  11. 前記空間ワーピングステップは、前記圧力画像をスケーリングするステップを備える、請求項9又は10記載の方法。
  12. 前記スケーリングは2つの直交方向に個別に倍率を与える、請求項11記載の方法。
  13. 前記第1の圧力画像の前記特性は一以上の所定の基準画像との比較である、請求項1−12のいずれかに記載の方法。
  14. 前記第1の圧力画像の前記特性は前記第1の圧力画像の比較である、請求項1−13のいずれかに記載の方法。
  15. 前記複数の第1の連続する圧力画像は一被検者から得られる、請求項10−14のいずれかに記載の方法。
  16. 前記複数の第1の連続する圧力画像は複数の被検者から得られる、請求項10−14のいずれかに記載の方法。
  17. 前記生成される連続する統計画像は前記又は各第1の圧力画像にほぼ等しい解像度を有する、請求項1−16のいずれかに記載の方法。
  18. 前記生成される連続する統計画像はt−マップ、F−マップ又は他の統計マップである、請求項1−17のいずれかに記載の方法。
  19. 前記統計画像を処理して更なる統計画像を生成するステップを更に備える、請求項1−18のいずれかに記載の方法。
  20. 前記更なる統計画像を生成するステップは、前記統計画像に統計的推定を実行するステップを備える、請求項19記載の方法。
  21. 前記更なる統計画像は、複数の第1の連続する圧力画像の特性の統計的有意性を表す段階的な連続統計推定画像である、請求項20記載の方法。
  22. 前記統計的推定は、t検定により生成される前記統計画像のt値を処理して確率値を生成する、請求項20又は21記載の方法。
  23. 前記統計的推定は、前記統計画像の平滑度及び/又は前記統計画像の少なくとも一つの幾何学特性を考慮する、請求項20,21又は22記載の方法。
  24. 前記統計的推定は、確率場理論、ボンフェローニ補正及び/又は非パラメトリック推定に基づく技術を含む、請求項20−23のいずれかに記載の方法。
  25. 前記第1の連続する圧力画像は足底面のペドバログラフィ画像である、請求項1−14のいずれかに記載の方法。
  26. 前記第1の足底面ペドバログラフィ画像及び前記第2の足底面ペドバログラフィ画像の各々は足底面全体の画像である、請求項25記載の方法。
  27. 各々が一以上の圧力画像の特性を表す複数の統計画像を生成するステップと、
    入力圧力画像を受信するステップと、
    前記入力圧力画像に基づいて前記複数の統計画像の一つを選択するステップと、
    を更に備える、請求項25又は26記載の方法。
  28. 前記複数の統計画像の各々に対して識別データを格納するステップと、
    前記複数の統計画像の前記選択された一つに対して識別データを出力するステップと、
    を更に備える請求項27に記載の方法。
  29. 前記第1の連続する圧力画像の各々は、第1の物体と第2の物体との間の圧力場の画像である、請求項1−24のいずれかに記載の方法。
  30. 前記第1の連続する圧力画像は、被検者とマットレスとの間の圧力場の画像である、請求項29記載の方法。
  31. 前記第1の連続する圧力画像は、被検者とシーティング装置との間の圧力場の画像である、請求項29記載の方法。
  32. 請求項1−31のいずれかに記載の方法を実行するように構成されたコンピュータプログラム。
  33. 請求項32のコンピュータプログラムを担持するコンピュータ読み取り可能媒体。
  34. プロセッサ読み取り可能な命令を格納するプログラムメモリと、前記プログラムメモリに格納された命令を読み出し、実行するプロセッサとを備え、
    前記プログラム読み取り可能命令は、前記プロセッサを、請求項1−25のいずれかに記載の方法を実行するように制御する命令を含む、コンピュータ装置。
  35. 複数の組の2次元圧力データを含む圧力データを表示する方法であって、該方法は、
    前記圧力データの3次元体積表現を表示するステップを備え、前記3次元体積表現の第1及び第2の次元は空間データを表し、前記3次元体積表現の第3の次元は時間データを表す、圧力データの表示方法。
  36. 前記3次元体積表現をグラフィック描画するステップを更に備える、請求項35記載の方法。
  37. 所定の圧力又は所定の範囲の圧力を表すデータから表面を計算するステップと、
    計算された表面を表示するグラフィックデータを生成するステップと、
    を備える、請求項35又は36記載の方法。
  38. 前記3次元体積表現は前記複数の圧力データから生成される3次元統計画像である、請求項35,36又は37記載の方法。
  39. 請求項35−38のいずれかに記載の方法を実行するように構成されたコンピュータプログラム。
  40. 請求項39のコンピュータプログラムを担持するコンピュータ読み取り可能媒体。
  41. プロセッサ読み取り可能命令を格納するプログラムメモリと、前記プログラムメモリに格納された命令を読み出し、実行するプロセッサとを備え、
    前記プログラム読み取り可能命令は、前記プロセッサを、請求項35−38のいずれかに記載の方法を実行するように制御する命令を含む、コンピュータ装置。
JP2010538893A 2007-12-21 2008-12-19 画像処理 Pending JP2011506022A (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
GBGB0725094.7A GB0725094D0 (en) 2007-12-21 2007-12-21 Image processing
PCT/GB2008/004185 WO2009081113A2 (en) 2007-12-21 2008-12-19 Image processing

Publications (1)

Publication Number Publication Date
JP2011506022A true JP2011506022A (ja) 2011-03-03

Family

ID=39048646

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2010538893A Pending JP2011506022A (ja) 2007-12-21 2008-12-19 画像処理

Country Status (8)

Country Link
US (1) US20100296752A1 (ja)
EP (1) EP2232439A2 (ja)
JP (1) JP2011506022A (ja)
CN (1) CN101933044A (ja)
AU (1) AU2008339679A1 (ja)
CA (1) CA2710143A1 (ja)
GB (1) GB0725094D0 (ja)
WO (1) WO2009081113A2 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013183809A (ja) * 2012-03-06 2013-09-19 Terumo Corp 足底圧計測装置及びその処理方法
JP2016502061A (ja) * 2012-12-04 2016-01-21 ゲナント ヴェルスボールグ インゴ シトーク 熱処理監視システム
JP2020018366A (ja) * 2018-07-30 2020-02-06 花王株式会社 歩行分析方法

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012016370A1 (en) * 2010-08-02 2012-02-09 Peking University Representative motion flow extraction for effective video classification and retrieval
US8754929B1 (en) * 2011-05-23 2014-06-17 John Prince Real time vergence control for 3D video capture and display
CN102846320A (zh) * 2011-06-30 2013-01-02 鸿富锦精密工业(深圳)有限公司 步姿校正系统及方法
EP2803057A4 (en) * 2012-01-13 2015-07-08 Enhanced Surface Dynamics Inc SYSTEM AND METHODS FOR RISK MANAGEMENT ANALYSIS OF A PRESSURE DETECTION SYSTEM
US20140324400A1 (en) * 2013-04-30 2014-10-30 Marquette University Gesture-Based Visualization System for Biomedical Imaging and Scientific Datasets
CN104082905B (zh) * 2014-06-18 2016-02-17 杭州华亭科技有限公司 多功能智能鞋垫及步态相似性检测方法
CN106154874A (zh) * 2015-04-16 2016-11-23 中兴通讯股份有限公司 一种智能鞋垫及其工作方法
CN105310654A (zh) * 2015-08-06 2016-02-10 跑动(厦门)信息科技有限公司 一种检测足内旋的方法和用于检测足内旋的智能鞋垫
US10661738B2 (en) * 2018-01-02 2020-05-26 Wipro Limited Method, system, and device for controlling internal systems within a vehicle based on user preferences
CN110840459B (zh) * 2019-11-19 2022-07-22 京东方科技集团股份有限公司 人体平衡能力获取方法及系统、计算机设备及介质
CN111407232A (zh) * 2020-03-31 2020-07-14 湖北民族大学 基于足底压力分布的足部运动特征提取方法及系统
CN112529782B (zh) * 2021-02-18 2021-06-01 棉捷(北京)网络科技有限公司 压力数据3d成像方法、装置、终端设备和可读存储介质
CN116576994B (zh) * 2023-05-12 2024-05-10 爱梦睡眠(珠海)智能科技有限公司 一种基于压电传感器的在离床辅助判断装置和方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5550928A (en) * 1992-12-15 1996-08-27 A.C. Nielsen Company Audience measurement system and method
DE19648272A1 (de) * 1996-11-21 1998-05-28 Emitec Emissionstechnologie Verfahren und Vorrichtung zur Bestimmung einer Zelldichte eines Wabenkörpers, insbesondere für einen Abgaskatalysator
JP3114668B2 (ja) * 1997-10-03 2000-12-04 日本電気株式会社 物体検出・背景除去方法、装置およびプログラムを記録した記録媒体
US7085401B2 (en) * 2001-10-31 2006-08-01 Infowrap Systems Ltd. Automatic object extraction
US7280710B1 (en) * 2002-05-24 2007-10-09 Cleveland Clinic Foundation Architecture for real-time 3D image registration
KR20050065543A (ko) * 2002-09-12 2005-06-29 엔라인 코포레이션 복소 영상의 획득 및 처리 시스템 및 방법
US7536043B2 (en) * 2003-08-18 2009-05-19 Siemens Medical Solutions Usa, Inc. Flow representation method and system for medical imaging
WO2006025963A2 (en) * 2004-07-16 2006-03-09 New York University Method, system and storage medium which includes instruction for analyzing anatomical structures
US8971984B2 (en) * 2005-04-04 2015-03-03 Hypermed Imaging, Inc. Hyperspectral technology for assessing and treating diabetic foot and tissue disease
EP1897495A4 (en) * 2005-06-28 2011-01-05 Konica Minolta Med & Graphic METHOD FOR THE DETECTION OF ANOMALIC SHADE CANDIDATES, DEVICE FOR DETECTING ANOMALIC SHADE CANDIDATES
US8040519B2 (en) * 2006-05-23 2011-10-18 Hitachi Medical Corporation Biological optical measurement apparatus
KR101150987B1 (ko) * 2007-07-16 2012-06-08 삼성전자주식회사 화상처리장치 및 그 제어방법

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013183809A (ja) * 2012-03-06 2013-09-19 Terumo Corp 足底圧計測装置及びその処理方法
JP2016502061A (ja) * 2012-12-04 2016-01-21 ゲナント ヴェルスボールグ インゴ シトーク 熱処理監視システム
US11013237B2 (en) 2012-12-04 2021-05-25 Ingo Stork Genannt Wersborg Heat treatment monitoring system
JP2020018366A (ja) * 2018-07-30 2020-02-06 花王株式会社 歩行分析方法
JP7096096B2 (ja) 2018-07-30 2022-07-05 花王株式会社 歩行分析方法及び歩行分析装置

Also Published As

Publication number Publication date
EP2232439A2 (en) 2010-09-29
WO2009081113A3 (en) 2009-10-29
AU2008339679A1 (en) 2009-07-02
CA2710143A1 (en) 2009-07-02
WO2009081113A2 (en) 2009-07-02
GB0725094D0 (en) 2008-01-30
US20100296752A1 (en) 2010-11-25
CN101933044A (zh) 2010-12-29

Similar Documents

Publication Publication Date Title
JP2011506022A (ja) 画像処理
CA2602218C (en) Method and system for characterization of knee joint morphology
US9561004B2 (en) Automated 3-D orthopedic assessments
US8126242B2 (en) Computer program products and methods for detection and tracking of rheumatoid arthritis
JP6741305B2 (ja) 脊椎骨姿勢推定装置
EP3114645B1 (en) Apparatus and method for generating and using a subject-specific statistical motion model
Ballester et al. Segmentation and measurement of brain structures in MRI including confidence bounds
Zhang et al. MPF-net: An effective framework for automated cobb angle estimation
US20120249551A1 (en) alignment of shapes of body parts from images
Reyes et al. Automatic digital biometry analysis based on depth maps
US20240277283A1 (en) Detecting spinal shape from optical scan
Davatzikos Measuring biological shape using geometry-based shape transformations
Strickland et al. Development of subject-specific geometric spine model through use of automated active contour segmentation and kinematic constraint-limited registration
Hanik et al. Bi-invariant dissimilarity measures for sample distributions in Lie groups
Loeckx et al. Non-rigid image registration using a statistical spline deformation model
Langs et al. Model-based erosion spotting and visualization in rheumatoid arthritis
Camara et al. Accuracy assessment of global and local atrophy measurement techniques with realistic simulated longitudinal data
Bui et al. Volterra series modelling and compensation of non-linear distortions caused by susceptibility difference artefacts related to the presence of ferromagnetic implants in magnetic resonance imaging
Hanik et al. Check for updates Bi-invariant Two-Sample Tests in Lie Groups for Shape Analysis Data from the Alzheimer's Disease Neuroimaging Initiative
Zeng et al. Colon surface registration using Ricci flow
Boisvert Articulated statistical shape models of the Spine
Mohr et al. A method for the automatic segmentation of autoradiographic image stacks and spatial normalization of functional cortical activity patterns

Legal Events

Date Code Title Description
A711 Notification of change in applicant

Free format text: JAPANESE INTERMEDIATE CODE: A711

Effective date: 20111101

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20111102

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20111102