JP4385139B2 - 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラム - Google Patents
変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラム Download PDFInfo
- Publication number
- JP4385139B2 JP4385139B2 JP2007556824A JP2007556824A JP4385139B2 JP 4385139 B2 JP4385139 B2 JP 4385139B2 JP 2007556824 A JP2007556824 A JP 2007556824A JP 2007556824 A JP2007556824 A JP 2007556824A JP 4385139 B2 JP4385139 B2 JP 4385139B2
- Authority
- JP
- Japan
- Prior art keywords
- displacement
- phase
- phase singularity
- procedure
- singularity
- 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.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/002—Measuring arrangements characterised by the use of optical techniques for measuring two or more coordinates
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Image Analysis (AREA)
Description
【技術分野】
【０００１】
本発明は変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチングプログラムに係り、特に、変位前後の位相特異点から変位を検出する変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラムに関する。
【背景技術】
【０００２】
近年、レーザースペックルパターンやランダムドットパターンなどのように空間的にランダムな構造をもつパターンやテクスチャを指標として物体の微小変位を非接触で計測する技術は非破壊検査や材料強度試験など産業応用上の重要な位置を占めている。パターンやテクスチャを用いて微小変位を検出する際の信号処理技術として相関法が知られている。
【０００３】
図３８、図３９は相関法の動作を説明するため図を示す。図３８（Ａ）は場の変位、図３８（Ｂ）相関法の概念図、図３８（Ｃ）は本測定法の概念図を示している。また、図３９（Ａ）は物体Ａの一方向への変位、図３９（Ｂ）は物体Ａの回転変位を説明するための図を示している。
【０００４】
相関法を用いて変位を検出する場合、図３８（Ａ）に矢印で示すように場の変位が位置により異なるときには図３８（Ｂ）に示すように各領域に分けて計算する必要があり、各点での相関値を決定するごとに２次元相関積分計算または２次元フーリエ変換・逆変換を行う必要があったため、１ピクセル以下での変位検出の分解能を上げようとしてピクセル内部での相関値の計算点数を増やすと計算時間が膨大になっていた。
【０００５】
また、従来の相関法では、図３９（Ａ）に示すように一方向の変位に対しては相関をとることにより変位量を求めることができるが、図３９（Ｂ）に示すように回転に対しては適用できなかった。
【０００６】
なお、位相情報を利用する相関関数の計算方法としては、位相限定相関法が知られている（特許文献１参照）。この相関法は、複素関数である合成スペクトルの振幅を一定化又は対数関数等により抑制し、位相情報のみからなる振幅限定複素合成スペクトルを作成してこれをフーリエ逆変換して相互相関関数を求めるものである。
【先行技術文献】
【特許文献１】
特許第３０３５６５４号
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかるに、従来の変位検出方法では、計算時間が膨大になるとともに、回転変位の検出が容易に行えないなどの問題点があった。
【０００８】
本発明は上記の点に鑑みてなされたもので、位相特異点の構造を特徴付ける特徴パラメータから所望の位相特異点を特定することにより変位前後の位相特異点を確実に特定し、その変位を容易、かつ、確実に、さらには、回転変位を含む種々の変位を検出できる変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
本発明は、取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、位相特異点の位相構造から所定の特徴パラメータを取得し、特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを有することを特徴とする。
【００１０】
このとき、取得される特徴パラメータは、位相特異点近傍の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度、位相特異点近傍の強度分布の形状の離心率、位相特異点近傍の強度の勾配であることを特徴とする。
【００１１】
また、特定された位相特異点の変位前後の位置を特徴パラメータに基づいて特定し、特定された変位前後の位相特異点の位置に基づいて位相特異点の変位を測定することを特徴とする。このとき、変位前後の位相特異点の位置の座標の差分の二乗平方根を位相特異点の変位量として算出することを特徴とする。
【００１２】
複数組の変位前後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる点を回転中心として求め、変位前後の位相特異点の位置と回転中心とからなる三角形の回転中心を挟む辺の角度を回転角として求めることを特徴とする。
【００１３】
さらに、本発明は、取得した画像の強度信号をフーリエ変換し、フーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行い、擬似位相情報を取得し、取得した擬似位相情報に基づいて位相特異点を取得することを特徴とする。
【００１４】
このとき、フーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振幅特性を持つラゲールガウスビームフィルタ特性によりフィルタリングすることを特徴とする。また、画像の強度信号の平均値を画像の各強度信号から減算した画像を元画像とすることを特徴とする。
【００１５】
さらに、フーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒルベルト変換によりフィルタリングを行い、擬似位相情報を取得することを特徴とする。
【００１６】
また、本発明は、変位前に特定された位相特異点と変位後に特定された位相特異点とをマッチングを行うポイントマッチング処理方法であって、変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とする。
【００１７】
また、このとき、本発明は、変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる位相特異点の組を変位前と変位後とで互いにマッチングする特徴点として特定することを特徴とする。
【００１８】
さらに、本発明は、変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
【００１９】
【数１】
が最小となる位相特異点の組を変位前後でマッチングする位相特異点の組として抽出することを特徴とする。
【００２０】
なお、位相特異点は、画像の強度情報に基づいて取得される位相特異点であることを特徴とし、位相特異点は特徴パラメータとして位置情報、離心率、及び、渦度、並びに、実軸と虚軸とが交わる角度を有することを特徴とする。
【発明の効果】
【００２１】
本発明によれば、位相特異点の位相構造から所定の特徴パラメータを取得し、特徴パラメータに基づいて位相特異点を特定することにより、変位前後の各位相特異点を一意に特定することができ、これによって、変位後の位相特異点を追跡できるなどの効果を得られる。
【図面の簡単な説明】
【００２２】
【図１】本発明の一実施例のシステム構成図である。
【図２】変位検出プログラムの処理フローチャートである。
【図３】位相特異点取得処理の処理フローチャートである。
【図４】位相特異点取得処理の動作説明図である。
【図５】位相特異点取得処理の動作説明図である。
【図６】位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
【図７】位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
【図８】位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
【図９】位相特異点特定処理の処理フローチャートである。
【図１０】離心率ｅ、渦度ωを取得する動作を説明するための図である。
【図１１】離心率ｅ、渦度ωを取得する動作を説明するための図である。
【図１２】離心率ｅ、渦度ωを取得する動作を説明するための図である。
【図１３】全位相特異点の位置の差を求めｘ方向，ｙ方向の変位をヒストグラムに表した図である。
【図１４】回転変位の検出動作を説明するための図である。
【図１５】シミュレーション実験を説明するための図である。
【図１６】回転前後の位相特異点の数をヒストグラムで表した図である。
【図１７】回転前後の位相特異点の回転角θiをヒストグラムで表した図である。
【図１８】回転前後の角度Δθの差をヒストグラムで表した図である。
【図１９】回転前後の対応する位相特異点の離心率ｅをヒストグラムで表した図である。
【図２０】回転前後の対応する位相特異点の渦度ωの差をヒストグラムで表した図である。
【図２１】ポイントマッチング処理の動作説明図である。
【図２２】ポイントマッチング処理のフローチャートである。
【図２３】マッチングアルゴリズムの処理動作を説明するための図である。
【図２４】マッチングアルゴリズムの処理動作を説明するための図である。
【図２５】マッチングアルゴリズムの処理動作を説明するための図である。
【図２６】実験・検証に用いたＬＧフィルタの振幅特性図である。
【図２７】マッチング処理の動作説明図である。
【図２８】マッチング処理の動作説明図である。
【図２９】実験・検証結果を説明するための図である。
【図３０】実験・検証結果を説明するための図である。
【図３１】実験・検証結果を説明するための図である。
【図３２】実験・検証結果を説明するための図である。
【図３３】実験・検証結果を説明するための図である。
【図３４】実験・検証結果を説明するための図である。
【図３５】実験・検証結果を説明するための図である。
【図３６】点間総和距離の差によりマッチングを行う動作を説明するための図である。
【図３７】点間総和距離の差によりマッチングを行う動作を説明するための図である。
【図３８】相関法の動作を説明するため図である。
【図３９】相関法の動作を説明するため図である。
【発明を実施するための形態】
【００２３】
〔システム構成〕
図１は本発明の一実施例のシステム構成図を示す。
【００２４】
本実施例の変位検出装置１０１は、画像検出部１１１及び信号処理部１１２から構成されている。
【００２５】
画像検出部１１１は、光源１２１及び撮像装置１２２から構成されている。画像検出部１１１は、光源１２１から測定対象１０２に光を照射し、撮像装置１２２により測定対象１０２の表面画像を撮像する。
【００２６】
画像検出部１１１で撮像された画像は、信号処理部１１２に供給される。信号処理部１１２は、例えば、コンピュータシステムであり、フレームメモリ１３１、処理部１３２、フーリエ変換部１３３、ＲＯＭ１３４、ハードディスクドライブ１３５、ディスプレイ１３６，入力装置１３７から構成されている。
【００２７】
フレームメモリ１３１は、画像検出部１１１からの画像をフレーム毎に記憶する。処理部１３２は、ハードディスクドライブ１３５に予めインストールされた変位検出プログラムに基づいて測定対象１０２の微小変位を検出するための処理を実行する。
【００２８】
フーリエ変換部１３３は、画像検出部１１１で撮像された画像に対してフーリエ変換を行う。メモリ１３４は、処理部１３２の作業用記憶領域として用いられる。ハードディスクドライブ１３５は、変位検出プログラムや画像検出部１１１で検出された画像データ、あるいは、変位検出プログラムでの途中結果、変位検出結果などのデータが記憶される。
【００２９】
ディスプレイ１３６は、ＬＣＤ、ＣＲＴなどから構成されており、画像、変換画像、途中結果、変位検出結果などの処理部１３２での処理結果を表示する。入力装置１３７は、キーボード、マウスなどから構成されており、データ入力や各種コマンドの入力するために用いられる。
【００３０】
本実施例の変位検出装置１０１は、インストールされた変位検出プログラムに基づいて測定対象１０２から取得した画像を解析して、その変位を検出する。
【００３１】
〔変位検出プログラム〕
図２は変位検出プログラムの処理フローチャートを示す。
【００３２】
処理部１３２は、まず、変位検出プログラムに従って、ステップＳ１−１で変位前後の対象画像をハードディスクドライブ１３５からメモリ１３４などに読み出す。
【００３３】
次に処理部１３２は、ステップＳ１−２で、読み出された変位前後の対象画像に対して各々位相特異点を取得するための位相特異点取得処理を実行する。
【００３４】
次に処理部１３２は、ステップＳ１−３で変位前後の対象画像から取得した位相特異点が持つ、実部と虚部とのゼロクロス角度、離心率、渦度などの特徴パラメータにより位相特異点を特定する位相特異点特定処理を実行する。
【００３５】
次に処理部１３２は、ステップＳ１−４で、実部と虚部とのゼロクロス角度、離心率、渦度などの位相特異点に特有の特徴パラメータにより変位前後で対応する位相特異点を特定し、変位前後の位相特異点の位置から対象画像の変位を検出する変位検出処理を実行する。
【００３６】
位相特異点によって、対象画像の変位を検出することで対象画像の変位を正確に特定することが可能となる。
【００３７】
〔複素解析信号取得処理〕
次にステップＳ１−２の位相特異点取得処理について説明する。
【００３８】
図３は位相特異点取得処理の処理フローチャート、図４、図５は位相特異点取得処理の動作説明図を示す。
【００３９】
ステップＳ１−２の位相特異点取得処理は、取得した画像の強度信号をフーリエ変換し、フーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行い、擬似位相情報を取得し、取得した擬似位相情報に基づいて位相特異点を取得する処理である。
【００４０】
処理部１３２は、まず、ステップＳ２−１で前処理を行う。ステップＳ２−１の前処理は、対象画像の全体の強度情報の平均値を求め、求めた全体画像の強度情報の平均値を対象画像の強度情報から減算する処理を行う処理であり、この処理によって、画像データ強度情報中の直流成分が除去される。直流成分が除去されることにより、位相特異点が強調される。
【００４１】
次に、処理部１３２は、ステップＳ２−２でフーリエ変換部１３３を制御して、強度情報に対して２次元高速フーリエ変換を行い、空間周波数スペクトルを得る。
【００４２】
例えば、図４において、画像検出部１１１で得られた物体のテクスチャの強度情報２１１の関数２０１を
【００４３】
【数２】
とする。この強度情報を２次元フーリエ変換することにより、空間周波数スペクトル
【００４４】
【数３】
を求めることができる。なお、図４において２１２は空間周波数スペクトルの振幅分布、２１３は位相分布を示している。空間周波数スペクトルは、実部のみをもつ強度情報をフーリエ変換したものなので、周波数原点に対して複素共役の対称性、すなわちエルミート対称性をもつ。ＦＦＴを関数Ｆとすると、例えば、式２０２のように表せる。
【００４５】
そこで、処理部１３２は、ステップＳ２−３で実部を基に虚部を作り出して複素解析信号を得るために、図４に２２１〜２２５で示すような分布を持つ、周波数領域上の複素共役の対称性を崩すようなフィルタリング処理を行う。フィルタリング処理では、例えば、ヒルベルト変換、スパイラル位相フィルタ、ＬＧ（ラゲールガウス）フィルタなどが用いられる。なお、図４において２２１はヒルベルト変換の振幅分布、２２２はスパイラル位相フィルタの位相分布、２２４はＬＧフィルタの振幅分布、２２５はＬＧフィルタの位相分布を示している。
【００４６】
ヒルベルト変換の演算式は、例えば、図４、２０３で表すことができる。また、スパイラル位相フィルタの演算式は、例えば、図４、２０４で表すことができる。ＬＧフィルタの振幅透過分布２２４の演算式は、図４、２０５（Ａ）で表すことができる。ＬＧフィルタは、中心を対称に複素共役を位相上でπだけずらすことにより，ドーナツ状の振幅透過率分布が得られる。なお、ＬＧフィルタの位相分布２２５の演算式は、例えば、図４、２０５（Ａ）に角度φを導入することにより図４、２０５（Ｂ）に示すように表される。
【００４７】
以上のように、演算式２０２に演算式２０３、２０４、２０５などを乗算することにより、本実施例で必要となる解析信号を得ることができる。また、エラーの原因となる高周波数成分を除去でき、直流成分を自動的、かつ、連続的にゼロにできる。さらに、ドーナツ状の強度分布により低周波数成分を取り除くことにより、位相特異点の数を増やすことができる。
【００４８】
さらに、本実施例ではランダムなパターンの強度情報を用いて解析を行っており、測定対象１０２の構造などによってスペックルという丸い斑点のようなパターンが多数存在することになる。このスペックルは、本実施例で指標としている位相特異点に対応関係がある。
【００４９】
このスペックルのサイズによりフーリエ変換したときの周波数領域に位相特異点の情報を多く含む周波数が変わる。例えば、スペックルのサイズが小さいと位相特異点の情報を含む周波数成分は高くなり、逆にサイズが大きいと位相特異点の情報を含む周波数成分は低くなる。つまり、位相特異点の情報を高精度に求めるにはドーナツ状の強度分布をもつＬＧフィルタが最適であり、ドーナツ状の幅のサイズをスペックルのサイズに合わせて調節することにより、より位相特異点の情報を高精度に求めることができる。一方、ローパスフィルタでは、低周波数成分を強調し、位相特異点の情報をうまく取り出すことができない。
【００５０】
以上のように、本実施例における解析信号の生成過程ではＬＧフィルタを用いることが最適となる。
【００５１】
なお、図４において２１４はフィルタリング処理後の振幅分布、２１５はフィルタリング処理後の位相分布を示しており、その演算式は、例えば、図４、２０６で表すことができる。
【００５２】
なお、２２１、２２３に示すようなヒルベルト変換またはスパイラル位相フィルタをかけても同種の解析信号を得ることはできるが、ヒルベルト変換またはスパイラル位相フィルタを使うためには予め直流成分を除去しておかなければならない。
【００５３】
処理部１３２は、ステップＳ２−４でフィルタをかけた空間周波数スペクトルをフーリエ逆変換して複素解析信号
【００５４】
【数４】
を得る。なお、図４において２１６はフーリエ逆変換後の振幅分布、２１７はフーリエ逆変換後の位相特性を示しており、その演算式は、例えば、図４、２０７で表すことができる。
【００５５】
図５（Ａ）はフーリエ逆変換後の振幅分布を拡大したものであり、図５（Ｂ）は図５（Ａ）の所定の部分を更に拡大したものである。
【００５６】
この複素解析信号の擬似位相情報から図５（Ｂ）に示すような位相特異点Ｐを取得することができる。
【００５７】
次に処理部１３２は、ステップＳ２−５で位相特異点の位置を１ピクセルより小さい分解能であるサブピクセルの分解能で決定する処理を行う。
【００５８】
図６、図７、図８は位相特異点の位置をサブピクセルで決定する動作を説明するための図を示す。
【００５９】
図６（Ａ）は擬似位相情報、図６（Ｂ）は位相特異点周辺画素の拡大図、図６（Ｃ）は実部、図６（Ｄ）は虚部を示している。また、図７（Ａ）は位相特異点近傍の擬似位相情報、図７（Ｂ）は位相特異点近傍の擬似位相情報の虚部、図７（Ｃ）は位相特異点近傍の擬似位相情報の実部、図７（Ｄ）は位相特異点近傍の擬似位相情報を直線補間したもの、図７（Ｅ）は位相特異点近傍の擬似位相情報の虚部を直線補間したもの、図７（Ｆ）は位相特異点近傍の擬似位相情報の実部を直線補間したものを示している。
【００６０】
擬似位相は(-π，π)までの値をとり、図６、図７において黒から白になるにしたがって値が大きくなる。
【００６１】
位相特異点Ｐは、その点を囲む閉経路上の位相勾配の積分が±２πの整数倍になり、また、複素解析信号の実部と虚部が共にゼロとなるという性質を持っている。
【００６２】
図６（Ｂ）、図７（Ｃ）に示すような互いに隣接しあう４画素を結ぶ正方形閉経路に対して位相勾配の積分が±２πの整数倍になるという性質を用いて位相特異点の存在している位置を決定していた。このため、１ピクセル以上の分解能をもって位相特異点の位置を決定することができない。
【００６３】
本実施例では、複素解析信号の実部と虚部が共にゼロとなるという性質を用いて１ピクセルの分解能より小さいサブピクセルの分解能で位相特異点の位置を決定するようにしている。本実施例では、図７（Ｂ）、図７（Ｃ）に示す特異点近傍の離散化された実部と虚部を平面で線形補間することにより図７（Ｅ）、図７（Ｆ）に示すように実部及び虚部を共に滑らかな構造として、実部がゼロとなる線分Ｌre0、虚部がゼロとなる線分Ｌim0を求める。図７（Ｄ）に示すように位相特異点の位置は実部がゼロとなる線分Ｌre0と虚部がゼロとなる線分Ｌim0の交わる点を位相特異点Ｐ0として設定することによりサブピクセルの分解能をもって位相特異点を決定することができる．
図８（Ａ）に示すように位相勾配の積分が±２πの整数倍になる性質を用いる方法では、１ピクセルの面積を持つ正方形の閉経路内のどこかに位相特異点が存在することはわかってもその位置を特定することができなかった。これに対して図８（Ｂ）に示すように複素解析信号の実部と虚部が共にゼロとなるという性質を用いて実部と虚部の内挿と実部と虚部のゼロの交差点から位相特異点を求めることにより１ピクセルより小さいサブピクセルの分解能をもって位相特異点の位置を決定できている。
【００６４】
〔位相特異点特定処理〕
図９は位相特異点特定処理の処理フローチャートを示す。
【００６５】
ステップＳ１−３の位相特異点特定処理は、位相特異点の位相構造から所定の特徴パラメータを取得し、特徴パラメータに基づいて位相特異点を特定する処理である。
【００６６】
処理部１３２は、ステップＳ３−１で位相特異点を特定する特徴パラメータとしてゼロクロス角を取得する。ゼロクロス角は、図７（Ｄ）、図８（Ｂ）に示す実部がゼロとなる線分Ｌre0と虚部がゼロとなる線分Ｌim0との交差する角度Δθである。
【００６７】
また、処理部１３２は、ステップＳ３−２で位相特異点を特定する特徴パラメータとして離心率ｅを取得し、ステップＳ３−３で渦度ωを取得する。
【００６８】
図１０、図１１、図１２は離心率ｅ、渦度ωを取得する動作を説明するための図を示す。図１０（Ａ）は擬似位相情報、図１０（Ｂ）は擬似位相情報の位相特異点近傍の振幅分布、図１０（Ｃ）は位相分布を示している。
【００６９】
図１０（Ｂ）に示すように振幅分布は位相特異点近傍に近づくにつれて暗くなっていて、これは値がゼロに近づいていることを示している。
【００７０】
また、図１１は特異点近傍の振幅分布の構造図を示している。図１１（Ａ）は振幅を等高線で表示したものを示しており、図１１（Ｂ）は三次元表示したものを示す。
【００７１】
図１１（Ａ）に示すように位相特異点近傍の等振幅の分布は略楕円形状をしている。一般に位相特異点近傍の振幅分布は略楕円形をしており、その形状は位相特異点に応じて異なっており、位相特異点を特定する特徴パラメータとして用いることができる。この楕円形状の特徴を二次曲線の特徴を表す値である離心率によって示すことによって、位相特異点を特定することが可能となる。
【００７２】
例えば、図１２に示すような楕円の場合、一般に数式で、
【００７３】
【数５】
で表せる。図１２に示す楕円の離心率ｅはａ：長軸半径，ｂ：短軸半径，ｆ：焦点距離とすると、離心率ｅは
【００７４】
【数６】
で表させる。
【００７５】
図１１（Ｂ）に示される振幅を三次元で表示した図を見ると位相特異点近傍の形状は楕円形の円錐のような構造をしていることが分かる。したがって、本実施例では、各位相特異点を特定する特徴パラメータの一つとして位相特異点近傍の渦度を求める。渦度とは微小部分の局所的な回転の強さと回転の方向を表す量である。図１１（Ａ）では振幅の等高線の間隔が狭いほど渦度は大きくなり、図１１（Ｂ）では三次元表示した位相特異点近傍の勾配（変化の度合い）が急になるほど渦度は大きくなる。すなわち、渦度は位相特異点近傍の振幅の大きさに依存する。
【００７６】
渦度ωは次式で表される。
【００７７】
【数７】
以上のように位相特異点は位置（ｘ、ｙ），ゼロクロス角Δθ，離心率e，渦度ωの特徴パラメータを持っていて、これらの特徴パラメータを用いて各位相特異点を一意に特定することができる。このとき、複数の特徴パラメータを用いることにより確実に位相特異点を特定することができる。
【００７８】
〔直線変位検出〕
次にステップＳ１−４の位相特異点変位検出処理について説明する。
【００７９】
位相変位検出処理では、変位前後の位相特異点を上記の特定方法により特定することにより変位量を検出する。
【００８０】
このとき、まず、変位前後で位相特異点をそれぞれ一意に特定するためのパラメータとして前述のように求めたゼロクロス角Δθ，離心率ｅ，渦度ωを用いて位相特異点を特定する。その方法は変位前の各パラメータをΔθ1、ｅ1、ω1、変位後の各パラメータをΔθ2、ｅ2、ω2とする。
【００８１】
上記パラメータ（Δθ1、ｅ1、ω1）、（Δθ2、ｅ2、ω2）のそれぞれの差（Δθ2−Δθ1、ｅ2−ｅ1、ω2−ω1）の絶対値がある値より小さい時に同一の位相特異点であるとみなす。各位相特異点をそれぞれ一意に特定できるため、変位前後で各位相特異点がどこに移動したのか分かる。このとき、変位前の位相特異点の位置（ｘ1、ｙ1）と変位後の位相特異点の位置（ｘ2、ｙ2）の差から変位量を求めることができる。
【００８２】
実際の変位量ΔＬは、Δｘ＝（ｘ2−ｘ1）、Δｙ＝（ｙ2−ｙ1）とすると、
【００８３】
【数８】
で表すことができる。
【００８４】
なお、図１３は、全位相特異点の位置の差を求めｘ方向，ｙ方向の変位をヒストグラムに表した図を示す。図１３（Ａ）はｘ方向の変位、図１３（Ｂ）はｙ方向の変位を示している。
【００８５】
図１３に示すように全位相特異点の位置がｘ方向に略平行に変位したことがわかる。
【００８６】
〔回転変位検出〕
また、本実施例では、各位相特異点を回転変位前後で特定できるため、回転の検出も可能である。
【００８７】
図１４は回転変位の検出動作を説明するための図を示す。
【００８８】
図１４（Ａ）に示すように回転前の２組の位相特異点Ｐ11、Ｐ12のｘｙ座標をそれぞれ
（ｘ1、ｙ1）、（ｘ2、ｙ2）
とし，回転後の対応する位相特異点Ｐ21、Ｐ22のｘｙ座標をそれぞれ
（ｘ'1、ｙ'1）、（ｘ'2、ｙ'2）
とすると、回転の中心
（ｘc、ｙc）
は図１４（Ａ）に示すように回転前後の位相特異点Ｐ11、Ｐ12と位相特異点Ｐ21、Ｐ22の垂直二等分線
【００８９】
【数９】
の交わる点となる。
【００９０】
一般的に回転の中心を求めるには２組の回転前後の位相特異点が分かれば十分であるが、実際は多数の位相特異点が存在するので回転前後の各位相特異点から得られる垂直二等分線が図１４（Ｂ）に示すに必ずしも１点では交わらない。このような場合には、回転の中心（ｘc、ｙc）と各垂直二等分線との距離の和が最小となるような回転の中心（ｘc、ｙc）を最小二乗法により求めるようにする。
【００９１】
以上により回転中心を決定できる。
【００９２】
次に図１４（Ｃ）に示すように回転前後の各位相特異点Ｐ11、Ｐ12と回転中心（ｘc、ｙc）からなる三角形に対して余弦定理を用いることによって、回転角θを求めることができる。
【００９３】
以上により、本実施例によれば、回転変位の検出を行うことが可能となる。
【００９４】
〔シミュレーション実験結果〕
実験で得た強度情報をコンピュータ内部で９０度だけ回転させて，回転させる前後の強度情報を用いてシミュレーション実験を行った。
【００９５】
図１５はシミュレーション実験を説明するための図を示す。なお、図１５は、画素数が１０２４×１０２４で、回転させる前後の強度情報からＬＧフィルタを用いて２つの解析信号を得たときの振幅分布と位相分布を表示したものである。図１５（Ａ）は回転変位前の擬似位相情報、図１５（Ｂ）は回転変位前の擬似位相情報の振幅分布、図１５（Ｃ）は回転変位前の擬似位相情報の位相分布、図１５（Ｄ）は回転変位後の擬似位相情報、図１５（Ｅ）は回転変位後の擬似位相情報の振幅分布、図１５（Ｆ）は回転変位後の擬似位相情報の位相分布を示している。
【００９６】
この回転前後の位相特異点の数は両方とも正の位相特異点が322個，負の位相区異点321個の合計643個で同数の位相特異点が存在していた。図１５（Ｂ）及び図１５（Ｅ）に示される振幅分布からも明らかなように回転前の三角で囲った２つのリング状の振幅が回転後には回転の中心（全画素の中心）に対して９０°回転している。
【００９７】
次に位相特異点を特定するための特徴パラメータであるRe=0とIm=0の直線の交わる角度：Δθ，離心率（eccentricity）：ｅ，渦度（vorticity）：ωを回転前後の全位相特異点に対して求める。回転前後で対応する位相特異点を探す条件は以下の条件を用いた。
【００９８】
条件１：Re=0とIm=0の直線の交わる角度Δθの差が±0.5度以下
条件２：回転前後の離心率ｅの差が±0.005以下
条件３：回転前後の渦度ωの差が±0.05以下
回転前の位相特異点に対応する回転後の位相特異点が１対１の対応になっているのかを調べた。
【００９９】
図１６は回転前の位相特異点に対応する回転後の位相特異点の数をヒストグラムで表した図を示す。
【０１００】
回転前後でほぼ全ての位相特異点が一意に特定できた。
シミュレーションで得られた測定結果の回転の中心は
（ｘc、ｙc）＝（511.5、511.5）
となった。これは設定した回転の中心
（ｘc、ｙc）＝（511.5、511.5）
と一致する。
【０１０１】
図１７は回転前の位相特異点に対応する回転後の位相特異点と回転の中心から求めたそれぞれの回転角θiをヒストグラムで表した図を示す。
【０１０２】
シミュレーション結果である回転角のヒストグラムは±0.001度の精度で全位相特異点が９０度回転していることを示している。シミュレーション実験では、測定対象１０２を９０度回転させており、図１７に示す測定結果も９０度回転していることを示しているので、この測定結果から本発明が回転に対しても有効に計測できることを確認できる。
【０１０３】
図１８は回転前後の対応する位相特異点の角度Δθの差をヒストグラムで表した図、図１９は回転前後の対応する位相特異点の離心率ｅの差をヒストグラムで表した図、図２０は回転前後の対応する位相特異点の渦度ωの差をヒストグラムで表した図を示す。
【０１０４】
図１８、図１９、図２０に示すようにゼロ付近に固まっており、誤差が少なく、位相特異点を確実に特定できていることがわかる。
【０１０５】
なお、上記実施例では、位相特異点を用いて回転変位を検出する処理について説明したが縮尺、移動、回転を簡単な処理により検出することも可能である。
【０１０６】
〔ポイントマッチング処理〕
ポイントマッチング処理は、任意の複数点の位置関係と、その点の空間構造の類似性を用いて変位前のパターンから変位後のパターンを検出する処理である。
【０１０７】
すなわち、位相構造が類似する点を候補点として、その候補点を含んだパターン内の点間の距離が最小となる点を変位後のパターンとして決定するというものである。
【０１０８】
以下にポイントマッチング処理の手順を説明する。
【０１０９】
まず、変位前後の２つの強度分布画像に対して各々前述の位相特異点を抽出するためのフィルタ演算を行うことにより位相特異点を抽出する。抽出した位相特異点近傍にある実部及び虚部を線形補間することによって、サブピクセルの精度で全位相特異点の位置を特定する。更に、位置が特定された各位相特異点に対して離心率ｅ、実軸と虚軸との角Δθ、渦度ωを求め、ラベリングを行う。
【０１１０】
図２１は、ポイントマッチング処理の動作説明図を示す。図２１（Ａ）は変位前の位相特異点の状態、図２１（Ｂ）は変位後の位相特異点の状態を示している。
【０１１１】
まず、図２１（Ａ）に示す位相特異点から黒丸で示す位相特異点を選択する。選択する位相特異点はその特徴パラメータである離心率ｅ、実軸と虚軸との角Δθ、渦度ωが他の多くの位相特異点とは明らかにことなるものを抽出することが望ましい。これによって、変位後の位相特異点を特定しやすくなる。
【０１１２】
図２１（Ａ）に示す変位前の画像から選択した位相特異点の特徴パラメータと位置と図２１（Ｂ）に示す変位後の画像から選択した位相特異点の特徴パラメータと位置とから図２１（Ｂ）に黒丸で示すように図２１（Ａ）に黒丸で示す位相特異点に対応する位相特異点を検出する。
【０１１３】
位相特異点の特徴パラメータ、離心率e、渦度ω、実軸と虚軸との角Δθは、変位前後で同じであるか、または、ほとんど変化しない。このため、変位前後の位相特異点を正確に特定することによって移動、回転、縮尺などを正確に測定することができる。
【０１１４】
図２２はポイントマッチング処理のフローチャートを示す。
【０１１５】
まず、処理部１３２は、ステップＳ４−１で、図３とともに説明した位相特異点取得処理によって変位前後の位相特異点を取得し、ラベリングする。
【０１１６】
次に、処理部１３２は、ステップＳ４−２で変位前の各位相特異点にほぼ同じ、あるいは、所定の範囲内に近似した特徴パラメータを持つ位相特異点を変位後の位相特異点から探索する。
【０１１７】
次に処理部１３２は、ステップＳ４−３で変位前の位相特異点から所定の２点を選択し、それを組み合わせて組み合わせＭとする。次に、処理部１３２は、ステップＳ４−４で変数ｋに「３」が代入される。
【０１１８】
次に、処理部１３２は、ステップＳ４−５で変数ｋが定数ｎより小さいか否かを判定する。なお、定数ｎは、探索する位相特異点の数である。
【０１１９】
次に処理部１３２は、ステップＳ４−６で選択された複数の位相特異点から次の探索対象となる位相特異点を抽出する。
【０１２０】
次に処理部１３２は、ステップＳ４−７で新たに抽出された位相特異点に対して特徴パラメータが近似する変位後の位相特異点を探索し、組み合わせＭとして追加する。
【０１２１】
次に処理部１３２は、ステップＳ４−８で組み合わせＭの位相特異点の特徴パラメータに対して下記の式（１）によって計算を行い、閾値以下のものを変位前の位相特異点と変位後の位相特異点との有力な組み合わせの候補として残す。
【０１２２】
【数１０】
式（１）においてＴｉとＰw(i)は変位前と変位後の点を表している。なお、ＴｉはＰw(i)と対応している。
Ｔ＝｛Ｔ１，Ｔ２，Ｔ３，・・・, Ｔｎ｝
は変位前の位相特異点の組み合わせを表している。そして、
Ｐ＝｛Ｐ１，Ｐ２，Ｐ３，・・・, Ｐｎ｝
は変位後の位相特異点の組み合わせを表している。
【０１２３】
また、｜Ｔi1−Ｔi2｜は、点Ｔi1と点Ｔi2との間の距離，αは縮尺のスケール(度合い)を表している。さらに、変数ｗはＴ→Ｐの対応関係を表している。
【０１２４】
なお、上記式（１）の関数をコスト関数（cost function）と呼ぶ。
【０１２５】
処理部１３２は、ステップＳ４−９で変数ｋに（ｋ＋１）を代入し、ステップＳ４−５に戻る。処理部１３２は、上記ステップＳ４−６〜Ｓ４−９をステップＳ４−５で変数ｋが定数ｎより大きくなるまで繰り返す。
【０１２６】
次に処理部１３２は、ステップＳ４−１０で残った組み合わせＭの位相特異点の特徴パラメータ間の距離が最小となるものをマッチングパターンとして抽出する。
【０１２７】
次に、処理部１３２は、ステップＳ４−１１で抽出した特徴パラメータ間の距離が最小となる位相特異点の変換パラメータ、例えば、平行移動(shift),回転(rotation)、縮尺(scale)を計算する。
【０１２８】
上記ポイントマッチング処理により、目標はパターン内の複数の点の位置関係において変位前後で最も良く一致するものを見つけ出すことができる。
【０１２９】
次に、上記マッチングアルゴリズムの処理動作を、図面を用いてさらに詳細に説明する。
【０１３０】
図２３、図２４、図２５はマッチングアルゴリズムの処理動作を説明するための図を示す。
【０１３１】
図２３（Ａ）は変位前の位相特異点の中から複数の位相特異点の組を選んだものであり、これをＴ空間とし、ｔi（ｉ＝１、．．．,ｎ）、Ｔi（Ｘi,Ｙi,ｅi,ωi,Δθi）で表す。図２３（Ｂ）は変位後の位相特異点でＰ空間とし、ｐj（ｊ＝１、．．．,ｋ）、Ｐj（Ｘj,Ｙj,ｅj,ωj,Δθj）で表される。
【０１３２】
まず、図２３（Ａ）に黒丸で示すように、変位前の位相特異点から任意の個数の位相特異点を選ぶ。
【０１３３】
変位前のパターン内のそれぞれの位相特異点ｔiの特徴パラメータである離心率ｅti、ゼロクロス角Δθti、渦度ωtiと変位後の位相特異点ｐjの特徴パラメータである離心率ｅpj、ゼロクロス角Δθpj、渦度ωpjとの差を閾値と比較する。それは、例えば、次式で表せる。
【０１３４】
【数１１】
なお、式（２）で２番目の分母の｜ωti＋ωpj｜は正規化するために付け加えたものである。実際のプログラム内ではこの形で使われている。また、ｅthは離心率の閾値、Δθthはゼロクロス角の閾値、ωthは渦度の閾値を示している。
【０１３５】
次にＴ空間の位相特異点の組とＰ空間の位相特異点の組とのマッチングを行う。
【０１３６】
図２４において、実線で示す丸印は正の位相特異点を示しており、破線で示した丸印は負の位相特異点を点線で示している。
【０１３７】
まず、図２４（Ａ）に示すＴ空間の２個の位相特異点の極性、正か負か符号と、その位置関係とに基づいて図２４（Ｂ）に破線の直線で示すようにＰ空間に存在しうる可能性のある２個の位相特異点の組を探索する。
【０１３８】
このマッチングでは移動(shift)，回転(rotation)，だけでなく縮尺(scale)も考慮しているので、Ｔ空間での初めの２つの組をＰ空間の中から２つを選ぶ組み合わせは、この例では正の位相特異点と負の位相特異点の組み合わせ全てであり、正が７個、負が６個なので４２通り存在する。
【０１３９】
なお、ここで、正の位相特異点と負の位相特異点というのは特徴パラメータの渦度(vorticity)の値が正のものを正の位相特異点、また、負のものを”負の位相特異点としている。
【０１４０】
この得られた位相特異点の組み合わせに対して，位相特異点の特徴パラメータの条件に合わない組み合わせを切り捨てる。このとき、各特徴パラメータは実験環境により変化しやすいことを考え、条件を緩くとることが好ましい。ここで言う条件とは、式（２）の閾値ｅth、Δθth、ωthのことである。
【０１４１】
次に、ステップＳ４−３で空間Ｔ内の２つの位相特異点ｔi1、ｔi2を選ぶ。この位相特異点ｔi1、ｔi2は少なくとも１つのマッチング候補点を持っている。位相特異点ｔi1はｎ1個のマッチング候補の点を、位相特異点ｔi2はｎ2個のマッチング候補点を持つ。したがって、変位後の空間Ｐ内の候補点の組み合わせはｎ1ｎ２通りあることがわかる。それをテンポラルマッチング（temporal matching）と呼ぶことにする。
【０１４２】
次に、図２５のようにＴ空間の位相特異点を１つ増やし、Ｐ空間の中から３点目として可能性のある点を探す。これにより３点目を含む組み合わせが得られる。なお、特徴パラメータを用いて明らかに間違った点は除外される。このような操作をＴ空間で選択される位相特異点がｎ個になるまで繰り返す。
【０１４３】
つまり、新しいマッチングポイントを加えつつ、上記テンポラルマッチングを最終的にすべての位相特異点のマッチング点が決定されるまで繰り返す。
【０１４４】
この操作を行うためには、新しく追加されていく点の地理的な位置関係が変位前と変位後で類似性があるかどうか調べる必要がある。
【０１４５】
式（１）では、テンポラルマッチングの類似性を調べると同時に距離の閾値は間違って検出されたマッチングポイントを除外する操作も行っていることになる。
【０１４６】
最終的にで得られた位相特異点の組み合わせＭの中から式（１）の評価式によりＴ空間の全位相特異点間の距離と対応するＰ空間の位相特異点間の空間的な距離及び特徴点の距離の差の総和が最小になるものを選び出してマッチングを行っている。
【０１４７】
さらに、図面を用いて説明を続ける。
【０１４８】
図２６、図２７はマッチング処理の動作説明図を示す。図２６（Ａ）〜（Ｃ）、図２７（Ａ）〜（Ｃ）は変位前の位相特異点の分布、すなわち、Ｔ空間における位相特異点の分布、図２６（Ｄ）〜（Ｆ）、図２７（Ｄ）〜（Ｆ）は変位後の位相特異点の分布、すなわち、Ｐ空間における位相特異点の分布を示す。
【０１４９】
まず、図２６（Ａ）、（Ｄ）に示すようにＴ空間及びＰ空間のそれぞれの位相特異点の特徴パラメータを計算し、注目ポイントを選択する。図２６（Ａ）では、ｔ1、ｔ2、ｔ4の組み合わせを選択した。
【０１５０】
次に、図２６（Ｂ）、（Ｅ）に示すようにＴ空間内のそれぞれの点の特徴パラメータと似たＰ空間の候補点Ｐ1〜Ｐ4を探索する。
【０１５１】
次に、図２６（Ｃ）、（Ｅ）に示すようにＴ空間内の２つの位相特異点ｔ1、ｔ2を選択して、Ｐ空間内の対応する２つの位相特異点ｐ1、ｐ2を選択し、組み合わせＭとする。
【０１５２】
次に、図２７（Ａ）、（Ｄ）に示すように既に選択した位相特異点以外で候補となる位相特異点ｔ3、ｔ4及びそれに対応するＰ空間内の位相特異点ｐ3、ｐ4を探索する。
【０１５３】
次に、図２７（Ｂ）、（Ｆ）に示すようにＴ空間の位相特異点ｔ1、ｔ2、ｔ4とＰ空間の位相特異点ｐ1、ｐ2、ｐ4の特徴パラメータを式（１）に代入して、式（１）の値を算出する。同様に、Ｔ空間の注目位相特異点ｔ1、ｔ2、ｔ4とＰ空間の位相特異点ｐ1、ｐ2、ｐ3の特徴パラメータを式（１）に代入して、式（１）の値を算出する。
【０１５４】
算出結果のうち閾値より小さい位相特異点の組を抽出する。例えば、ここでは、組み合わせＭ１（ｐ1、ｐ2、ｐ4）及びＭ2（ｐ1、ｐ2、ｐ3）が抽出されたとする。抽出された位相特異点の組Ｍ1、Ｍ2をマッチングした組み合わせとして抽出する。
【０１５５】
次に図２７（Ｃ）、（Ｇ）に示すようにマッチングした組Ｍ1、Ｍ2のうち最もコスト、すなわち、式（１）の値が低い組Ｍ1をマッチングした組として抽出する。
【０１５６】
以上は、説明を簡単にするために注目する位相特異点の数を３点とした場合について説明したが、３点以上の位相特異点に注目して同様な処理を行うことも可能である。
【０１５７】
次に、ステップＳ４−１１で算出される変換パラメータについて説明する。すなわち、位相特異点の移動が平行移動、回転、縮尺の３つのパラメータに分解できることについて説明する。
【０１５８】
一般に、流体力学の観点からどんな複雑な動きもそれを平行移動(transition)，回転(rotation)，縮尺(scale)で表すことができることが知られており、変位前後座標の回転関係式は次式で表すことができる。
【０１５９】
【数１２】
ここで、（ｘi、ｙi）は変位前の位相特異点の座標、（ｘi'、ｙi'）は変位後の位相特異点の座標、αは縮尺、θは回転角である。
【０１６０】
【数１３】
は変位前の位相特異点の重心、
【０１６１】
【数１４】
は変位後の位相特異点の重心を表している。
まず、はじめに、マッチングを表すのに適した関数について説明する。
【０１６２】
この関数はそれぞれ変位前後で対応している位相特異点の差を足し合わせた形で書ける。
【０１６３】
【数１５】
ここで、Ｃ＝αcosθ、Ｓ＝αsinθを表している。Ｅの値の最小値を見つけるためにはそれぞれのパラメータをまず計算しなければならない。例えば、
【０１６４】
【数１６】
などである。
【０１６５】
ここで、
【０１６６】
【数１７】
とおくと、
【０１６７】
【数１８】
から次式が得られる。
【０１６８】
【数１９】
これを解くと次のようになる。
【０１６９】
【数２０】
同様にして
【０１７０】
【数２１】
から次式が得られる。
【０１７１】
【数２２】
これを解くと次式が得られる。
【０１７２】
【数２３】
ＳとＣの定義より以下のパラメータが得られる。
【０１７３】
【数２４】
このように、位相特異点の移動は、平行移動Δｘ、Δｙ、回転θ、縮尺αの３つのパラメータで表せることがわかる。この３つのパラメータを変換パラメータと呼び、ステップＳ４−１１で算出される。
【０１７４】
〔実験・検証〕
次に、上記マッチング処理方法の実験・検証結果について説明する。
【０１７５】
図２８は実験・検証に用いられたＬＧフィルタの振幅特性図、図２９〜図３５は実験・検証結果を説明するための図を示す。図２９は実験画像を示す図であり、図２９（Ａ）は変位前の画像、図２９（Ｂ）は変位後の画像を示している。図３０は図２９に示す画像から抽出された位相特異点を示す画像であり、図２９（Ａ）は変位前の位相特異点の分布、図２９（Ｂ）は変位後の位相特異点の分布を示している。図３１は選択した位相特異点を拡大した図であり、図３１（Ａ）は変位前、図３１（Ｂ）は変位後の拡大図を示している。
【０１７６】
まず、今回の実験では、対象としてみどりふぐを使用し、高速ＣＣＤカメラを用いて秒間３０枚、連続撮影した画像を用いてマッチング処理を検証している。
【０１７７】
撮影した画像に対してＬＧフィルタ演算を行い、位相特異点を得る。このとき、ＬＧフィルタのフィルタ幅は、ω2＝３０としている。なお、ＬＧフィルタは、spiral phase filterにバンドパスフィルタの特性を加えたものであり、図２８に示すような振幅特性とした。なお、ＬＧフィルタの振幅幅ωを広くとれば、出てくる位相特異点の数は多くなり、逆に狭くすれば少なくなる。
【０１７８】
図２９に示す変位前後の画像に対して今回提案したマッチングアルゴリズムを適用した例について説明する。図２９に示す画像の変位前後の位相特異点の分布は図３０に示すようになる。図３０（Ａ）に示す変位前の位相特異点から黒丸で示される６つの位相特異点を選択してマッチングを行った。
【０１７９】
図３０（Ｂ）において、黒丸で示す位相特異点は変位前の黒丸で示す位相特異点に対応しており、実際にマッチングアルゴリズムにより検出されたパターンである。
【０１８０】
図３１は変位前に選んだ位相特異点の組み合わせと、ポイントマッチングアルゴリズムによって検出された変位後の位相特異点の組み合わせを拡大図示したものである。このように変位前，変位後のパターンを見比べると移動，回転，縮尺を組み合わせた変形が起こっているがマッチングアルゴリズムによって検証された。
【０１８１】
異なる２枚の画像に対して異なったパターンに着目して、本実施例のポイントマッチング処理を採用した場合の結果についても検証した。
【０１８２】
図３２は異なる２枚の画像に対して異なったパターンに着目して、本実施例のポイントマッチング処理を採用した場合の変位前後の位相特異点の分布を示す図であり、図３２（Ａ）は変位前の分布、図３２（Ｂ）は変位後の分布を示している。図３３は図３２に示した位相特異点の変位を示す図であり、図３３（Ｂ）は選択された位相特異点を拡大した図である。
【０１８３】
この検証例においても図３２（Ａ）に示すように変位前の４つの位相特異点が変位後、図３２（Ｂ）に示す位置に移動し、図３３に示すように注目した４つの位相特異点が変位していることが検証できた。
【０１８４】
さらに同じ作業をふぐの連続撮影画像１００枚に対して行った。
【０１８５】
図３４は連続した画像について実験・検証例を行った場合の画像を示す図である。図３４（Ａ）は最初の画像、図３４（Ｂ）は最後の画像を示している。図３５は図３４に示す画像から６つの位相特異点に注目し、その軌跡を追跡したものである。
【０１８６】
このように、本実施例のポイントマッチング処理を用いることにより、位相特異点を連続して追跡できる。
【０１８７】
さらに、任意の位相特異点の組み合わせに対して位相特異点が持つ位相構造と他の位相特異点との距離とによって対応する位相特異点の組を検出し、それらの位相特異点の組を追跡することによってランダムな動きを追跡できることがわかる。特徴パラメータが変化してしまう理由で計測不可能であった大きな変位を伴った物体の運動変位計測が可能となった。
【０１８８】
また、本実施例では、位相特異点の組のうち式（１）が最小となる組み合わせを検出することにより、マッチングを行ったが、位相特異点の点間総和距離の差を算出し、変位前のパターンと変位後のパターンとで点間総和距離の差が最小となる位相特異点を検出することによりマッチングを行うようにしてもよい。
【０１８９】
図３６、図３７は点間総和距離の差によりマッチングを行う動作を説明するための図を示す。
【０１９０】
図３６（Ａ）、図３７（Ａ）は変位前の位相特異点の分布、図３６（Ｂ）、図３７（Ｂ）は変位後の位相特異点の分布を示す。
【０１９１】
まず、変位前の画像から位相特異点ｐ1〜ｐ4を抽出するとともに、変位後の画像から位相特異点ｐ1'、ｐ2'、ｐ3'、ｐ4'を抽出し、変位前の位相特異点ｐ1〜ｐ4と変位後の位相特異点ｐ1'、ｐ2'、ｐ3'、ｐ4'との点間総和距離の差Ｓ1を求める。点間総和距離の差Ｓ1は、
Ｓ1＝（｜ｐ1ｐ3｜−｜ｐ1'ｐ3'｜）2＋（｜ｐ1ｐ2｜−｜ｐ1'ｐ2'｜）2
＋（｜ｐ2ｐ3｜−｜ｐ2'ｐ3'｜）2＋（｜ｐ3ｐ4｜−｜ｐ3'ｐ4'｜）2
＋（｜ｐ1ｐ4｜−｜ｐ1'ｐ4'｜）2＋（｜ｐ2ｐ4｜−｜ｐ2'ｐ4'｜）2
で求められる。
ここで、｜ｐiｐj｜は、点ｐiと点ｐjとの間の距離を表している。
【０１９２】
また、変位後の画像から位相特異点ｐ1'、ｐ2'、ｐ3'、ｐ4'とは異なる位相特異点ｐ1"、ｐ2"、ｐ3"、ｐ4"を抽出し、変位前の位相特異点ｐ1〜ｐ4と変位後の位相特異点ｐ1"、ｐ2"、ｐ3"、ｐ4"との点間総和距離の差Ｓ2を求める。点間総和距離の差Ｓ2は、
Ｓ2＝（｜ｐ1ｐ3｜−｜ｐ1"ｐ3"｜）2＋（｜ｐ1ｐ2｜−｜ｐ1"ｐ2"｜）2
＋（｜ｐ2ｐ3｜−｜ｐ2"ｐ3"｜）2＋（｜ｐ3ｐ4｜−｜ｐ3"ｐ4"｜）2
＋（｜ｐ1ｐ4｜−｜ｐ1"ｐ4"｜）2＋（｜ｐ2ｐ4｜−｜ｐ2"ｐ4"｜）2
で求められる。
【０１９３】
このとき、例えば、Ｓ2＜Ｓ1であれば、変位前の画像から位相特異点ｐ1〜ｐ4は、位相特異点ｐ1"、ｐ2"、ｐ3"、ｐ4"に対応すると判定する。
【０１９４】
なお、ここでは、説明を簡単にするために４点について説明したが、２点以上であれば、点間総和距離の差によりマッチングを行うことが可能であることは言うまでもない。
【０１９５】
本実施例のポイントマッチング処理は、複数の位相特異点の位置関係を用いるので縮尺、移動そして、回転により変位前後で完全にその位置関係が一致しなくてもポイントマッチングを行うことで一番近い組を選び出すことができるので、本発明はあらゆる変位を計測することができる。
【０１９６】
なお、本実施例では、特徴点として位置情報、離心率、渦度、実部がゼロとなる線分Ｌre0と虚部がゼロとなる線分Ｌim0との交叉角度であるゼロクロス角を特徴パラメータとして持つ位相特異点を例に説明したが、特徴パラメータはこれに限定されるものではなく、例えば、位置情報だけであってもよい。
【０１９７】
以上、本実施例によれば、干渉計なしで空間的なランダムパターンの強度情報を位相情報に対応づけ、位相情報の中に存在する全位相特異点の特徴を基に変位前後で各位相特異点がどこに移動したか特定することにより、計測物体の全体の変位計測に対しても、物体の局所領域における変位計測のいずれに対しても適用することができ、さらに場の回転の計測をも実現した変位検出装置を提供することできる。
【０１９８】
本国際出願は、２００６年２月１日に出願した日本国特許出願２００６−０２４８４６号、及び、２００６年１１月８日に出願した日本国特許出願２００６−３０３２３８号に基づく優先権を主張するものであり、日本国特許出願２００６−０２４８４６号、及び、日本国特許出願２００６−３０３２３８号の全内容を本国際出願に援用する。
【符号の説明】
【０１９９】
１０１ 変位検出装置、１０２ 測定対象
１１１ 画像検出部、１１２ 信号処理部
１２１ 光源、１２２ 撮像装置
１３１ フレームメモリ、１３２ 処理部、１３３ フーリエ変換部
１３４ ＲＯＭ
１３５ ハードディスクドライブ、１３６ ディスプレイ、１３７ 入力装置
【０００１】
本発明は変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチングプログラムに係り、特に、変位前後の位相特異点から変位を検出する変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラムに関する。
【背景技術】
【０００２】
近年、レーザースペックルパターンやランダムドットパターンなどのように空間的にランダムな構造をもつパターンやテクスチャを指標として物体の微小変位を非接触で計測する技術は非破壊検査や材料強度試験など産業応用上の重要な位置を占めている。パターンやテクスチャを用いて微小変位を検出する際の信号処理技術として相関法が知られている。
【０００３】
図３８、図３９は相関法の動作を説明するため図を示す。図３８（Ａ）は場の変位、図３８（Ｂ）相関法の概念図、図３８（Ｃ）は本測定法の概念図を示している。また、図３９（Ａ）は物体Ａの一方向への変位、図３９（Ｂ）は物体Ａの回転変位を説明するための図を示している。
【０００４】
相関法を用いて変位を検出する場合、図３８（Ａ）に矢印で示すように場の変位が位置により異なるときには図３８（Ｂ）に示すように各領域に分けて計算する必要があり、各点での相関値を決定するごとに２次元相関積分計算または２次元フーリエ変換・逆変換を行う必要があったため、１ピクセル以下での変位検出の分解能を上げようとしてピクセル内部での相関値の計算点数を増やすと計算時間が膨大になっていた。
【０００５】
また、従来の相関法では、図３９（Ａ）に示すように一方向の変位に対しては相関をとることにより変位量を求めることができるが、図３９（Ｂ）に示すように回転に対しては適用できなかった。
【０００６】
なお、位相情報を利用する相関関数の計算方法としては、位相限定相関法が知られている（特許文献１参照）。この相関法は、複素関数である合成スペクトルの振幅を一定化又は対数関数等により抑制し、位相情報のみからなる振幅限定複素合成スペクトルを作成してこれをフーリエ逆変換して相互相関関数を求めるものである。
【先行技術文献】
【特許文献１】
特許第３０３５６５４号
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかるに、従来の変位検出方法では、計算時間が膨大になるとともに、回転変位の検出が容易に行えないなどの問題点があった。
【０００８】
本発明は上記の点に鑑みてなされたもので、位相特異点の構造を特徴付ける特徴パラメータから所望の位相特異点を特定することにより変位前後の位相特異点を確実に特定し、その変位を容易、かつ、確実に、さらには、回転変位を含む種々の変位を検出できる変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
本発明は、取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、位相特異点の位相構造から所定の特徴パラメータを取得し、特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを有することを特徴とする。
【００１０】
このとき、取得される特徴パラメータは、位相特異点近傍の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度、位相特異点近傍の強度分布の形状の離心率、位相特異点近傍の強度の勾配であることを特徴とする。
【００１１】
また、特定された位相特異点の変位前後の位置を特徴パラメータに基づいて特定し、特定された変位前後の位相特異点の位置に基づいて位相特異点の変位を測定することを特徴とする。このとき、変位前後の位相特異点の位置の座標の差分の二乗平方根を位相特異点の変位量として算出することを特徴とする。
【００１２】
複数組の変位前後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる点を回転中心として求め、変位前後の位相特異点の位置と回転中心とからなる三角形の回転中心を挟む辺の角度を回転角として求めることを特徴とする。
【００１３】
さらに、本発明は、取得した画像の強度信号をフーリエ変換し、フーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行い、擬似位相情報を取得し、取得した擬似位相情報に基づいて位相特異点を取得することを特徴とする。
【００１４】
このとき、フーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振幅特性を持つラゲールガウスビームフィルタ特性によりフィルタリングすることを特徴とする。また、画像の強度信号の平均値を画像の各強度信号から減算した画像を元画像とすることを特徴とする。
【００１５】
さらに、フーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒルベルト変換によりフィルタリングを行い、擬似位相情報を取得することを特徴とする。
【００１６】
また、本発明は、変位前に特定された位相特異点と変位後に特定された位相特異点とをマッチングを行うポイントマッチング処理方法であって、変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とする。
【００１７】
また、このとき、本発明は、変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる位相特異点の組を変位前と変位後とで互いにマッチングする特徴点として特定することを特徴とする。
【００１８】
さらに、本発明は、変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
【００１９】
【数１】
が最小となる位相特異点の組を変位前後でマッチングする位相特異点の組として抽出することを特徴とする。
【００２０】
なお、位相特異点は、画像の強度情報に基づいて取得される位相特異点であることを特徴とし、位相特異点は特徴パラメータとして位置情報、離心率、及び、渦度、並びに、実軸と虚軸とが交わる角度を有することを特徴とする。
【発明の効果】
【００２１】
本発明によれば、位相特異点の位相構造から所定の特徴パラメータを取得し、特徴パラメータに基づいて位相特異点を特定することにより、変位前後の各位相特異点を一意に特定することができ、これによって、変位後の位相特異点を追跡できるなどの効果を得られる。
【図面の簡単な説明】
【００２２】
【図１】本発明の一実施例のシステム構成図である。
【図２】変位検出プログラムの処理フローチャートである。
【図３】位相特異点取得処理の処理フローチャートである。
【図４】位相特異点取得処理の動作説明図である。
【図５】位相特異点取得処理の動作説明図である。
【図６】位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
【図７】位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
【図８】位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
【図９】位相特異点特定処理の処理フローチャートである。
【図１０】離心率ｅ、渦度ωを取得する動作を説明するための図である。
【図１１】離心率ｅ、渦度ωを取得する動作を説明するための図である。
【図１２】離心率ｅ、渦度ωを取得する動作を説明するための図である。
【図１３】全位相特異点の位置の差を求めｘ方向，ｙ方向の変位をヒストグラムに表した図である。
【図１４】回転変位の検出動作を説明するための図である。
【図１５】シミュレーション実験を説明するための図である。
【図１６】回転前後の位相特異点の数をヒストグラムで表した図である。
【図１７】回転前後の位相特異点の回転角θiをヒストグラムで表した図である。
【図１８】回転前後の角度Δθの差をヒストグラムで表した図である。
【図１９】回転前後の対応する位相特異点の離心率ｅをヒストグラムで表した図である。
【図２０】回転前後の対応する位相特異点の渦度ωの差をヒストグラムで表した図である。
【図２１】ポイントマッチング処理の動作説明図である。
【図２２】ポイントマッチング処理のフローチャートである。
【図２３】マッチングアルゴリズムの処理動作を説明するための図である。
【図２４】マッチングアルゴリズムの処理動作を説明するための図である。
【図２５】マッチングアルゴリズムの処理動作を説明するための図である。
【図２６】実験・検証に用いたＬＧフィルタの振幅特性図である。
【図２７】マッチング処理の動作説明図である。
【図２８】マッチング処理の動作説明図である。
【図２９】実験・検証結果を説明するための図である。
【図３０】実験・検証結果を説明するための図である。
【図３１】実験・検証結果を説明するための図である。
【図３２】実験・検証結果を説明するための図である。
【図３３】実験・検証結果を説明するための図である。
【図３４】実験・検証結果を説明するための図である。
【図３５】実験・検証結果を説明するための図である。
【図３６】点間総和距離の差によりマッチングを行う動作を説明するための図である。
【図３７】点間総和距離の差によりマッチングを行う動作を説明するための図である。
【図３８】相関法の動作を説明するため図である。
【図３９】相関法の動作を説明するため図である。
【発明を実施するための形態】
【００２３】
〔システム構成〕
図１は本発明の一実施例のシステム構成図を示す。
【００２４】
本実施例の変位検出装置１０１は、画像検出部１１１及び信号処理部１１２から構成されている。
【００２５】
画像検出部１１１は、光源１２１及び撮像装置１２２から構成されている。画像検出部１１１は、光源１２１から測定対象１０２に光を照射し、撮像装置１２２により測定対象１０２の表面画像を撮像する。
【００２６】
画像検出部１１１で撮像された画像は、信号処理部１１２に供給される。信号処理部１１２は、例えば、コンピュータシステムであり、フレームメモリ１３１、処理部１３２、フーリエ変換部１３３、ＲＯＭ１３４、ハードディスクドライブ１３５、ディスプレイ１３６，入力装置１３７から構成されている。
【００２７】
フレームメモリ１３１は、画像検出部１１１からの画像をフレーム毎に記憶する。処理部１３２は、ハードディスクドライブ１３５に予めインストールされた変位検出プログラムに基づいて測定対象１０２の微小変位を検出するための処理を実行する。
【００２８】
フーリエ変換部１３３は、画像検出部１１１で撮像された画像に対してフーリエ変換を行う。メモリ１３４は、処理部１３２の作業用記憶領域として用いられる。ハードディスクドライブ１３５は、変位検出プログラムや画像検出部１１１で検出された画像データ、あるいは、変位検出プログラムでの途中結果、変位検出結果などのデータが記憶される。
【００２９】
ディスプレイ１３６は、ＬＣＤ、ＣＲＴなどから構成されており、画像、変換画像、途中結果、変位検出結果などの処理部１３２での処理結果を表示する。入力装置１３７は、キーボード、マウスなどから構成されており、データ入力や各種コマンドの入力するために用いられる。
【００３０】
本実施例の変位検出装置１０１は、インストールされた変位検出プログラムに基づいて測定対象１０２から取得した画像を解析して、その変位を検出する。
【００３１】
〔変位検出プログラム〕
図２は変位検出プログラムの処理フローチャートを示す。
【００３２】
処理部１３２は、まず、変位検出プログラムに従って、ステップＳ１−１で変位前後の対象画像をハードディスクドライブ１３５からメモリ１３４などに読み出す。
【００３３】
次に処理部１３２は、ステップＳ１−２で、読み出された変位前後の対象画像に対して各々位相特異点を取得するための位相特異点取得処理を実行する。
【００３４】
次に処理部１３２は、ステップＳ１−３で変位前後の対象画像から取得した位相特異点が持つ、実部と虚部とのゼロクロス角度、離心率、渦度などの特徴パラメータにより位相特異点を特定する位相特異点特定処理を実行する。
【００３５】
次に処理部１３２は、ステップＳ１−４で、実部と虚部とのゼロクロス角度、離心率、渦度などの位相特異点に特有の特徴パラメータにより変位前後で対応する位相特異点を特定し、変位前後の位相特異点の位置から対象画像の変位を検出する変位検出処理を実行する。
【００３６】
位相特異点によって、対象画像の変位を検出することで対象画像の変位を正確に特定することが可能となる。
【００３７】
〔複素解析信号取得処理〕
次にステップＳ１−２の位相特異点取得処理について説明する。
【００３８】
図３は位相特異点取得処理の処理フローチャート、図４、図５は位相特異点取得処理の動作説明図を示す。
【００３９】
ステップＳ１−２の位相特異点取得処理は、取得した画像の強度信号をフーリエ変換し、フーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行い、擬似位相情報を取得し、取得した擬似位相情報に基づいて位相特異点を取得する処理である。
【００４０】
処理部１３２は、まず、ステップＳ２−１で前処理を行う。ステップＳ２−１の前処理は、対象画像の全体の強度情報の平均値を求め、求めた全体画像の強度情報の平均値を対象画像の強度情報から減算する処理を行う処理であり、この処理によって、画像データ強度情報中の直流成分が除去される。直流成分が除去されることにより、位相特異点が強調される。
【００４１】
次に、処理部１３２は、ステップＳ２−２でフーリエ変換部１３３を制御して、強度情報に対して２次元高速フーリエ変換を行い、空間周波数スペクトルを得る。
【００４２】
例えば、図４において、画像検出部１１１で得られた物体のテクスチャの強度情報２１１の関数２０１を
【００４３】
【数２】
とする。この強度情報を２次元フーリエ変換することにより、空間周波数スペクトル
【００４４】
【数３】
を求めることができる。なお、図４において２１２は空間周波数スペクトルの振幅分布、２１３は位相分布を示している。空間周波数スペクトルは、実部のみをもつ強度情報をフーリエ変換したものなので、周波数原点に対して複素共役の対称性、すなわちエルミート対称性をもつ。ＦＦＴを関数Ｆとすると、例えば、式２０２のように表せる。
【００４５】
そこで、処理部１３２は、ステップＳ２−３で実部を基に虚部を作り出して複素解析信号を得るために、図４に２２１〜２２５で示すような分布を持つ、周波数領域上の複素共役の対称性を崩すようなフィルタリング処理を行う。フィルタリング処理では、例えば、ヒルベルト変換、スパイラル位相フィルタ、ＬＧ（ラゲールガウス）フィルタなどが用いられる。なお、図４において２２１はヒルベルト変換の振幅分布、２２２はスパイラル位相フィルタの位相分布、２２４はＬＧフィルタの振幅分布、２２５はＬＧフィルタの位相分布を示している。
【００４６】
ヒルベルト変換の演算式は、例えば、図４、２０３で表すことができる。また、スパイラル位相フィルタの演算式は、例えば、図４、２０４で表すことができる。ＬＧフィルタの振幅透過分布２２４の演算式は、図４、２０５（Ａ）で表すことができる。ＬＧフィルタは、中心を対称に複素共役を位相上でπだけずらすことにより，ドーナツ状の振幅透過率分布が得られる。なお、ＬＧフィルタの位相分布２２５の演算式は、例えば、図４、２０５（Ａ）に角度φを導入することにより図４、２０５（Ｂ）に示すように表される。
【００４７】
以上のように、演算式２０２に演算式２０３、２０４、２０５などを乗算することにより、本実施例で必要となる解析信号を得ることができる。また、エラーの原因となる高周波数成分を除去でき、直流成分を自動的、かつ、連続的にゼロにできる。さらに、ドーナツ状の強度分布により低周波数成分を取り除くことにより、位相特異点の数を増やすことができる。
【００４８】
さらに、本実施例ではランダムなパターンの強度情報を用いて解析を行っており、測定対象１０２の構造などによってスペックルという丸い斑点のようなパターンが多数存在することになる。このスペックルは、本実施例で指標としている位相特異点に対応関係がある。
【００４９】
このスペックルのサイズによりフーリエ変換したときの周波数領域に位相特異点の情報を多く含む周波数が変わる。例えば、スペックルのサイズが小さいと位相特異点の情報を含む周波数成分は高くなり、逆にサイズが大きいと位相特異点の情報を含む周波数成分は低くなる。つまり、位相特異点の情報を高精度に求めるにはドーナツ状の強度分布をもつＬＧフィルタが最適であり、ドーナツ状の幅のサイズをスペックルのサイズに合わせて調節することにより、より位相特異点の情報を高精度に求めることができる。一方、ローパスフィルタでは、低周波数成分を強調し、位相特異点の情報をうまく取り出すことができない。
【００５０】
以上のように、本実施例における解析信号の生成過程ではＬＧフィルタを用いることが最適となる。
【００５１】
なお、図４において２１４はフィルタリング処理後の振幅分布、２１５はフィルタリング処理後の位相分布を示しており、その演算式は、例えば、図４、２０６で表すことができる。
【００５２】
なお、２２１、２２３に示すようなヒルベルト変換またはスパイラル位相フィルタをかけても同種の解析信号を得ることはできるが、ヒルベルト変換またはスパイラル位相フィルタを使うためには予め直流成分を除去しておかなければならない。
【００５３】
処理部１３２は、ステップＳ２−４でフィルタをかけた空間周波数スペクトルをフーリエ逆変換して複素解析信号
【００５４】
【数４】
を得る。なお、図４において２１６はフーリエ逆変換後の振幅分布、２１７はフーリエ逆変換後の位相特性を示しており、その演算式は、例えば、図４、２０７で表すことができる。
【００５５】
図５（Ａ）はフーリエ逆変換後の振幅分布を拡大したものであり、図５（Ｂ）は図５（Ａ）の所定の部分を更に拡大したものである。
【００５６】
この複素解析信号の擬似位相情報から図５（Ｂ）に示すような位相特異点Ｐを取得することができる。
【００５７】
次に処理部１３２は、ステップＳ２−５で位相特異点の位置を１ピクセルより小さい分解能であるサブピクセルの分解能で決定する処理を行う。
【００５８】
図６、図７、図８は位相特異点の位置をサブピクセルで決定する動作を説明するための図を示す。
【００５９】
図６（Ａ）は擬似位相情報、図６（Ｂ）は位相特異点周辺画素の拡大図、図６（Ｃ）は実部、図６（Ｄ）は虚部を示している。また、図７（Ａ）は位相特異点近傍の擬似位相情報、図７（Ｂ）は位相特異点近傍の擬似位相情報の虚部、図７（Ｃ）は位相特異点近傍の擬似位相情報の実部、図７（Ｄ）は位相特異点近傍の擬似位相情報を直線補間したもの、図７（Ｅ）は位相特異点近傍の擬似位相情報の虚部を直線補間したもの、図７（Ｆ）は位相特異点近傍の擬似位相情報の実部を直線補間したものを示している。
【００６０】
擬似位相は(-π，π)までの値をとり、図６、図７において黒から白になるにしたがって値が大きくなる。
【００６１】
位相特異点Ｐは、その点を囲む閉経路上の位相勾配の積分が±２πの整数倍になり、また、複素解析信号の実部と虚部が共にゼロとなるという性質を持っている。
【００６２】
図６（Ｂ）、図７（Ｃ）に示すような互いに隣接しあう４画素を結ぶ正方形閉経路に対して位相勾配の積分が±２πの整数倍になるという性質を用いて位相特異点の存在している位置を決定していた。このため、１ピクセル以上の分解能をもって位相特異点の位置を決定することができない。
【００６３】
本実施例では、複素解析信号の実部と虚部が共にゼロとなるという性質を用いて１ピクセルの分解能より小さいサブピクセルの分解能で位相特異点の位置を決定するようにしている。本実施例では、図７（Ｂ）、図７（Ｃ）に示す特異点近傍の離散化された実部と虚部を平面で線形補間することにより図７（Ｅ）、図７（Ｆ）に示すように実部及び虚部を共に滑らかな構造として、実部がゼロとなる線分Ｌre0、虚部がゼロとなる線分Ｌim0を求める。図７（Ｄ）に示すように位相特異点の位置は実部がゼロとなる線分Ｌre0と虚部がゼロとなる線分Ｌim0の交わる点を位相特異点Ｐ0として設定することによりサブピクセルの分解能をもって位相特異点を決定することができる．
図８（Ａ）に示すように位相勾配の積分が±２πの整数倍になる性質を用いる方法では、１ピクセルの面積を持つ正方形の閉経路内のどこかに位相特異点が存在することはわかってもその位置を特定することができなかった。これに対して図８（Ｂ）に示すように複素解析信号の実部と虚部が共にゼロとなるという性質を用いて実部と虚部の内挿と実部と虚部のゼロの交差点から位相特異点を求めることにより１ピクセルより小さいサブピクセルの分解能をもって位相特異点の位置を決定できている。
【００６４】
〔位相特異点特定処理〕
図９は位相特異点特定処理の処理フローチャートを示す。
【００６５】
ステップＳ１−３の位相特異点特定処理は、位相特異点の位相構造から所定の特徴パラメータを取得し、特徴パラメータに基づいて位相特異点を特定する処理である。
【００６６】
処理部１３２は、ステップＳ３−１で位相特異点を特定する特徴パラメータとしてゼロクロス角を取得する。ゼロクロス角は、図７（Ｄ）、図８（Ｂ）に示す実部がゼロとなる線分Ｌre0と虚部がゼロとなる線分Ｌim0との交差する角度Δθである。
【００６７】
また、処理部１３２は、ステップＳ３−２で位相特異点を特定する特徴パラメータとして離心率ｅを取得し、ステップＳ３−３で渦度ωを取得する。
【００６８】
図１０、図１１、図１２は離心率ｅ、渦度ωを取得する動作を説明するための図を示す。図１０（Ａ）は擬似位相情報、図１０（Ｂ）は擬似位相情報の位相特異点近傍の振幅分布、図１０（Ｃ）は位相分布を示している。
【００６９】
図１０（Ｂ）に示すように振幅分布は位相特異点近傍に近づくにつれて暗くなっていて、これは値がゼロに近づいていることを示している。
【００７０】
また、図１１は特異点近傍の振幅分布の構造図を示している。図１１（Ａ）は振幅を等高線で表示したものを示しており、図１１（Ｂ）は三次元表示したものを示す。
【００７１】
図１１（Ａ）に示すように位相特異点近傍の等振幅の分布は略楕円形状をしている。一般に位相特異点近傍の振幅分布は略楕円形をしており、その形状は位相特異点に応じて異なっており、位相特異点を特定する特徴パラメータとして用いることができる。この楕円形状の特徴を二次曲線の特徴を表す値である離心率によって示すことによって、位相特異点を特定することが可能となる。
【００７２】
例えば、図１２に示すような楕円の場合、一般に数式で、
【００７３】
【数５】
で表せる。図１２に示す楕円の離心率ｅはａ：長軸半径，ｂ：短軸半径，ｆ：焦点距離とすると、離心率ｅは
【００７４】
【数６】
で表させる。
【００７５】
図１１（Ｂ）に示される振幅を三次元で表示した図を見ると位相特異点近傍の形状は楕円形の円錐のような構造をしていることが分かる。したがって、本実施例では、各位相特異点を特定する特徴パラメータの一つとして位相特異点近傍の渦度を求める。渦度とは微小部分の局所的な回転の強さと回転の方向を表す量である。図１１（Ａ）では振幅の等高線の間隔が狭いほど渦度は大きくなり、図１１（Ｂ）では三次元表示した位相特異点近傍の勾配（変化の度合い）が急になるほど渦度は大きくなる。すなわち、渦度は位相特異点近傍の振幅の大きさに依存する。
【００７６】
渦度ωは次式で表される。
【００７７】
【数７】
以上のように位相特異点は位置（ｘ、ｙ），ゼロクロス角Δθ，離心率e，渦度ωの特徴パラメータを持っていて、これらの特徴パラメータを用いて各位相特異点を一意に特定することができる。このとき、複数の特徴パラメータを用いることにより確実に位相特異点を特定することができる。
【００７８】
〔直線変位検出〕
次にステップＳ１−４の位相特異点変位検出処理について説明する。
【００７９】
位相変位検出処理では、変位前後の位相特異点を上記の特定方法により特定することにより変位量を検出する。
【００８０】
このとき、まず、変位前後で位相特異点をそれぞれ一意に特定するためのパラメータとして前述のように求めたゼロクロス角Δθ，離心率ｅ，渦度ωを用いて位相特異点を特定する。その方法は変位前の各パラメータをΔθ1、ｅ1、ω1、変位後の各パラメータをΔθ2、ｅ2、ω2とする。
【００８１】
上記パラメータ（Δθ1、ｅ1、ω1）、（Δθ2、ｅ2、ω2）のそれぞれの差（Δθ2−Δθ1、ｅ2−ｅ1、ω2−ω1）の絶対値がある値より小さい時に同一の位相特異点であるとみなす。各位相特異点をそれぞれ一意に特定できるため、変位前後で各位相特異点がどこに移動したのか分かる。このとき、変位前の位相特異点の位置（ｘ1、ｙ1）と変位後の位相特異点の位置（ｘ2、ｙ2）の差から変位量を求めることができる。
【００８２】
実際の変位量ΔＬは、Δｘ＝（ｘ2−ｘ1）、Δｙ＝（ｙ2−ｙ1）とすると、
【００８３】
【数８】
で表すことができる。
【００８４】
なお、図１３は、全位相特異点の位置の差を求めｘ方向，ｙ方向の変位をヒストグラムに表した図を示す。図１３（Ａ）はｘ方向の変位、図１３（Ｂ）はｙ方向の変位を示している。
【００８５】
図１３に示すように全位相特異点の位置がｘ方向に略平行に変位したことがわかる。
【００８６】
〔回転変位検出〕
また、本実施例では、各位相特異点を回転変位前後で特定できるため、回転の検出も可能である。
【００８７】
図１４は回転変位の検出動作を説明するための図を示す。
【００８８】
図１４（Ａ）に示すように回転前の２組の位相特異点Ｐ11、Ｐ12のｘｙ座標をそれぞれ
（ｘ1、ｙ1）、（ｘ2、ｙ2）
とし，回転後の対応する位相特異点Ｐ21、Ｐ22のｘｙ座標をそれぞれ
（ｘ'1、ｙ'1）、（ｘ'2、ｙ'2）
とすると、回転の中心
（ｘc、ｙc）
は図１４（Ａ）に示すように回転前後の位相特異点Ｐ11、Ｐ12と位相特異点Ｐ21、Ｐ22の垂直二等分線
【００８９】
【数９】
の交わる点となる。
【００９０】
一般的に回転の中心を求めるには２組の回転前後の位相特異点が分かれば十分であるが、実際は多数の位相特異点が存在するので回転前後の各位相特異点から得られる垂直二等分線が図１４（Ｂ）に示すに必ずしも１点では交わらない。このような場合には、回転の中心（ｘc、ｙc）と各垂直二等分線との距離の和が最小となるような回転の中心（ｘc、ｙc）を最小二乗法により求めるようにする。
【００９１】
以上により回転中心を決定できる。
【００９２】
次に図１４（Ｃ）に示すように回転前後の各位相特異点Ｐ11、Ｐ12と回転中心（ｘc、ｙc）からなる三角形に対して余弦定理を用いることによって、回転角θを求めることができる。
【００９３】
以上により、本実施例によれば、回転変位の検出を行うことが可能となる。
【００９４】
〔シミュレーション実験結果〕
実験で得た強度情報をコンピュータ内部で９０度だけ回転させて，回転させる前後の強度情報を用いてシミュレーション実験を行った。
【００９５】
図１５はシミュレーション実験を説明するための図を示す。なお、図１５は、画素数が１０２４×１０２４で、回転させる前後の強度情報からＬＧフィルタを用いて２つの解析信号を得たときの振幅分布と位相分布を表示したものである。図１５（Ａ）は回転変位前の擬似位相情報、図１５（Ｂ）は回転変位前の擬似位相情報の振幅分布、図１５（Ｃ）は回転変位前の擬似位相情報の位相分布、図１５（Ｄ）は回転変位後の擬似位相情報、図１５（Ｅ）は回転変位後の擬似位相情報の振幅分布、図１５（Ｆ）は回転変位後の擬似位相情報の位相分布を示している。
【００９６】
この回転前後の位相特異点の数は両方とも正の位相特異点が322個，負の位相区異点321個の合計643個で同数の位相特異点が存在していた。図１５（Ｂ）及び図１５（Ｅ）に示される振幅分布からも明らかなように回転前の三角で囲った２つのリング状の振幅が回転後には回転の中心（全画素の中心）に対して９０°回転している。
【００９７】
次に位相特異点を特定するための特徴パラメータであるRe=0とIm=0の直線の交わる角度：Δθ，離心率（eccentricity）：ｅ，渦度（vorticity）：ωを回転前後の全位相特異点に対して求める。回転前後で対応する位相特異点を探す条件は以下の条件を用いた。
【００９８】
条件１：Re=0とIm=0の直線の交わる角度Δθの差が±0.5度以下
条件２：回転前後の離心率ｅの差が±0.005以下
条件３：回転前後の渦度ωの差が±0.05以下
回転前の位相特異点に対応する回転後の位相特異点が１対１の対応になっているのかを調べた。
【００９９】
図１６は回転前の位相特異点に対応する回転後の位相特異点の数をヒストグラムで表した図を示す。
【０１００】
回転前後でほぼ全ての位相特異点が一意に特定できた。
シミュレーションで得られた測定結果の回転の中心は
（ｘc、ｙc）＝（511.5、511.5）
となった。これは設定した回転の中心
（ｘc、ｙc）＝（511.5、511.5）
と一致する。
【０１０１】
図１７は回転前の位相特異点に対応する回転後の位相特異点と回転の中心から求めたそれぞれの回転角θiをヒストグラムで表した図を示す。
【０１０２】
シミュレーション結果である回転角のヒストグラムは±0.001度の精度で全位相特異点が９０度回転していることを示している。シミュレーション実験では、測定対象１０２を９０度回転させており、図１７に示す測定結果も９０度回転していることを示しているので、この測定結果から本発明が回転に対しても有効に計測できることを確認できる。
【０１０３】
図１８は回転前後の対応する位相特異点の角度Δθの差をヒストグラムで表した図、図１９は回転前後の対応する位相特異点の離心率ｅの差をヒストグラムで表した図、図２０は回転前後の対応する位相特異点の渦度ωの差をヒストグラムで表した図を示す。
【０１０４】
図１８、図１９、図２０に示すようにゼロ付近に固まっており、誤差が少なく、位相特異点を確実に特定できていることがわかる。
【０１０５】
なお、上記実施例では、位相特異点を用いて回転変位を検出する処理について説明したが縮尺、移動、回転を簡単な処理により検出することも可能である。
【０１０６】
〔ポイントマッチング処理〕
ポイントマッチング処理は、任意の複数点の位置関係と、その点の空間構造の類似性を用いて変位前のパターンから変位後のパターンを検出する処理である。
【０１０７】
すなわち、位相構造が類似する点を候補点として、その候補点を含んだパターン内の点間の距離が最小となる点を変位後のパターンとして決定するというものである。
【０１０８】
以下にポイントマッチング処理の手順を説明する。
【０１０９】
まず、変位前後の２つの強度分布画像に対して各々前述の位相特異点を抽出するためのフィルタ演算を行うことにより位相特異点を抽出する。抽出した位相特異点近傍にある実部及び虚部を線形補間することによって、サブピクセルの精度で全位相特異点の位置を特定する。更に、位置が特定された各位相特異点に対して離心率ｅ、実軸と虚軸との角Δθ、渦度ωを求め、ラベリングを行う。
【０１１０】
図２１は、ポイントマッチング処理の動作説明図を示す。図２１（Ａ）は変位前の位相特異点の状態、図２１（Ｂ）は変位後の位相特異点の状態を示している。
【０１１１】
まず、図２１（Ａ）に示す位相特異点から黒丸で示す位相特異点を選択する。選択する位相特異点はその特徴パラメータである離心率ｅ、実軸と虚軸との角Δθ、渦度ωが他の多くの位相特異点とは明らかにことなるものを抽出することが望ましい。これによって、変位後の位相特異点を特定しやすくなる。
【０１１２】
図２１（Ａ）に示す変位前の画像から選択した位相特異点の特徴パラメータと位置と図２１（Ｂ）に示す変位後の画像から選択した位相特異点の特徴パラメータと位置とから図２１（Ｂ）に黒丸で示すように図２１（Ａ）に黒丸で示す位相特異点に対応する位相特異点を検出する。
【０１１３】
位相特異点の特徴パラメータ、離心率e、渦度ω、実軸と虚軸との角Δθは、変位前後で同じであるか、または、ほとんど変化しない。このため、変位前後の位相特異点を正確に特定することによって移動、回転、縮尺などを正確に測定することができる。
【０１１４】
図２２はポイントマッチング処理のフローチャートを示す。
【０１１５】
まず、処理部１３２は、ステップＳ４−１で、図３とともに説明した位相特異点取得処理によって変位前後の位相特異点を取得し、ラベリングする。
【０１１６】
次に、処理部１３２は、ステップＳ４−２で変位前の各位相特異点にほぼ同じ、あるいは、所定の範囲内に近似した特徴パラメータを持つ位相特異点を変位後の位相特異点から探索する。
【０１１７】
次に処理部１３２は、ステップＳ４−３で変位前の位相特異点から所定の２点を選択し、それを組み合わせて組み合わせＭとする。次に、処理部１３２は、ステップＳ４−４で変数ｋに「３」が代入される。
【０１１８】
次に、処理部１３２は、ステップＳ４−５で変数ｋが定数ｎより小さいか否かを判定する。なお、定数ｎは、探索する位相特異点の数である。
【０１１９】
次に処理部１３２は、ステップＳ４−６で選択された複数の位相特異点から次の探索対象となる位相特異点を抽出する。
【０１２０】
次に処理部１３２は、ステップＳ４−７で新たに抽出された位相特異点に対して特徴パラメータが近似する変位後の位相特異点を探索し、組み合わせＭとして追加する。
【０１２１】
次に処理部１３２は、ステップＳ４−８で組み合わせＭの位相特異点の特徴パラメータに対して下記の式（１）によって計算を行い、閾値以下のものを変位前の位相特異点と変位後の位相特異点との有力な組み合わせの候補として残す。
【０１２２】
【数１０】
式（１）においてＴｉとＰw(i)は変位前と変位後の点を表している。なお、ＴｉはＰw(i)と対応している。
Ｔ＝｛Ｔ１，Ｔ２，Ｔ３，・・・, Ｔｎ｝
は変位前の位相特異点の組み合わせを表している。そして、
Ｐ＝｛Ｐ１，Ｐ２，Ｐ３，・・・, Ｐｎ｝
は変位後の位相特異点の組み合わせを表している。
【０１２３】
また、｜Ｔi1−Ｔi2｜は、点Ｔi1と点Ｔi2との間の距離，αは縮尺のスケール(度合い)を表している。さらに、変数ｗはＴ→Ｐの対応関係を表している。
【０１２４】
なお、上記式（１）の関数をコスト関数（cost function）と呼ぶ。
【０１２５】
処理部１３２は、ステップＳ４−９で変数ｋに（ｋ＋１）を代入し、ステップＳ４−５に戻る。処理部１３２は、上記ステップＳ４−６〜Ｓ４−９をステップＳ４−５で変数ｋが定数ｎより大きくなるまで繰り返す。
【０１２６】
次に処理部１３２は、ステップＳ４−１０で残った組み合わせＭの位相特異点の特徴パラメータ間の距離が最小となるものをマッチングパターンとして抽出する。
【０１２７】
次に、処理部１３２は、ステップＳ４−１１で抽出した特徴パラメータ間の距離が最小となる位相特異点の変換パラメータ、例えば、平行移動(shift),回転(rotation)、縮尺(scale)を計算する。
【０１２８】
上記ポイントマッチング処理により、目標はパターン内の複数の点の位置関係において変位前後で最も良く一致するものを見つけ出すことができる。
【０１２９】
次に、上記マッチングアルゴリズムの処理動作を、図面を用いてさらに詳細に説明する。
【０１３０】
図２３、図２４、図２５はマッチングアルゴリズムの処理動作を説明するための図を示す。
【０１３１】
図２３（Ａ）は変位前の位相特異点の中から複数の位相特異点の組を選んだものであり、これをＴ空間とし、ｔi（ｉ＝１、．．．,ｎ）、Ｔi（Ｘi,Ｙi,ｅi,ωi,Δθi）で表す。図２３（Ｂ）は変位後の位相特異点でＰ空間とし、ｐj（ｊ＝１、．．．,ｋ）、Ｐj（Ｘj,Ｙj,ｅj,ωj,Δθj）で表される。
【０１３２】
まず、図２３（Ａ）に黒丸で示すように、変位前の位相特異点から任意の個数の位相特異点を選ぶ。
【０１３３】
変位前のパターン内のそれぞれの位相特異点ｔiの特徴パラメータである離心率ｅti、ゼロクロス角Δθti、渦度ωtiと変位後の位相特異点ｐjの特徴パラメータである離心率ｅpj、ゼロクロス角Δθpj、渦度ωpjとの差を閾値と比較する。それは、例えば、次式で表せる。
【０１３４】
【数１１】
なお、式（２）で２番目の分母の｜ωti＋ωpj｜は正規化するために付け加えたものである。実際のプログラム内ではこの形で使われている。また、ｅthは離心率の閾値、Δθthはゼロクロス角の閾値、ωthは渦度の閾値を示している。
【０１３５】
次にＴ空間の位相特異点の組とＰ空間の位相特異点の組とのマッチングを行う。
【０１３６】
図２４において、実線で示す丸印は正の位相特異点を示しており、破線で示した丸印は負の位相特異点を点線で示している。
【０１３７】
まず、図２４（Ａ）に示すＴ空間の２個の位相特異点の極性、正か負か符号と、その位置関係とに基づいて図２４（Ｂ）に破線の直線で示すようにＰ空間に存在しうる可能性のある２個の位相特異点の組を探索する。
【０１３８】
このマッチングでは移動(shift)，回転(rotation)，だけでなく縮尺(scale)も考慮しているので、Ｔ空間での初めの２つの組をＰ空間の中から２つを選ぶ組み合わせは、この例では正の位相特異点と負の位相特異点の組み合わせ全てであり、正が７個、負が６個なので４２通り存在する。
【０１３９】
なお、ここで、正の位相特異点と負の位相特異点というのは特徴パラメータの渦度(vorticity)の値が正のものを正の位相特異点、また、負のものを”負の位相特異点としている。
【０１４０】
この得られた位相特異点の組み合わせに対して，位相特異点の特徴パラメータの条件に合わない組み合わせを切り捨てる。このとき、各特徴パラメータは実験環境により変化しやすいことを考え、条件を緩くとることが好ましい。ここで言う条件とは、式（２）の閾値ｅth、Δθth、ωthのことである。
【０１４１】
次に、ステップＳ４−３で空間Ｔ内の２つの位相特異点ｔi1、ｔi2を選ぶ。この位相特異点ｔi1、ｔi2は少なくとも１つのマッチング候補点を持っている。位相特異点ｔi1はｎ1個のマッチング候補の点を、位相特異点ｔi2はｎ2個のマッチング候補点を持つ。したがって、変位後の空間Ｐ内の候補点の組み合わせはｎ1ｎ２通りあることがわかる。それをテンポラルマッチング（temporal matching）と呼ぶことにする。
【０１４２】
次に、図２５のようにＴ空間の位相特異点を１つ増やし、Ｐ空間の中から３点目として可能性のある点を探す。これにより３点目を含む組み合わせが得られる。なお、特徴パラメータを用いて明らかに間違った点は除外される。このような操作をＴ空間で選択される位相特異点がｎ個になるまで繰り返す。
【０１４３】
つまり、新しいマッチングポイントを加えつつ、上記テンポラルマッチングを最終的にすべての位相特異点のマッチング点が決定されるまで繰り返す。
【０１４４】
この操作を行うためには、新しく追加されていく点の地理的な位置関係が変位前と変位後で類似性があるかどうか調べる必要がある。
【０１４５】
式（１）では、テンポラルマッチングの類似性を調べると同時に距離の閾値は間違って検出されたマッチングポイントを除外する操作も行っていることになる。
【０１４６】
最終的にで得られた位相特異点の組み合わせＭの中から式（１）の評価式によりＴ空間の全位相特異点間の距離と対応するＰ空間の位相特異点間の空間的な距離及び特徴点の距離の差の総和が最小になるものを選び出してマッチングを行っている。
【０１４７】
さらに、図面を用いて説明を続ける。
【０１４８】
図２６、図２７はマッチング処理の動作説明図を示す。図２６（Ａ）〜（Ｃ）、図２７（Ａ）〜（Ｃ）は変位前の位相特異点の分布、すなわち、Ｔ空間における位相特異点の分布、図２６（Ｄ）〜（Ｆ）、図２７（Ｄ）〜（Ｆ）は変位後の位相特異点の分布、すなわち、Ｐ空間における位相特異点の分布を示す。
【０１４９】
まず、図２６（Ａ）、（Ｄ）に示すようにＴ空間及びＰ空間のそれぞれの位相特異点の特徴パラメータを計算し、注目ポイントを選択する。図２６（Ａ）では、ｔ1、ｔ2、ｔ4の組み合わせを選択した。
【０１５０】
次に、図２６（Ｂ）、（Ｅ）に示すようにＴ空間内のそれぞれの点の特徴パラメータと似たＰ空間の候補点Ｐ1〜Ｐ4を探索する。
【０１５１】
次に、図２６（Ｃ）、（Ｅ）に示すようにＴ空間内の２つの位相特異点ｔ1、ｔ2を選択して、Ｐ空間内の対応する２つの位相特異点ｐ1、ｐ2を選択し、組み合わせＭとする。
【０１５２】
次に、図２７（Ａ）、（Ｄ）に示すように既に選択した位相特異点以外で候補となる位相特異点ｔ3、ｔ4及びそれに対応するＰ空間内の位相特異点ｐ3、ｐ4を探索する。
【０１５３】
次に、図２７（Ｂ）、（Ｆ）に示すようにＴ空間の位相特異点ｔ1、ｔ2、ｔ4とＰ空間の位相特異点ｐ1、ｐ2、ｐ4の特徴パラメータを式（１）に代入して、式（１）の値を算出する。同様に、Ｔ空間の注目位相特異点ｔ1、ｔ2、ｔ4とＰ空間の位相特異点ｐ1、ｐ2、ｐ3の特徴パラメータを式（１）に代入して、式（１）の値を算出する。
【０１５４】
算出結果のうち閾値より小さい位相特異点の組を抽出する。例えば、ここでは、組み合わせＭ１（ｐ1、ｐ2、ｐ4）及びＭ2（ｐ1、ｐ2、ｐ3）が抽出されたとする。抽出された位相特異点の組Ｍ1、Ｍ2をマッチングした組み合わせとして抽出する。
【０１５５】
次に図２７（Ｃ）、（Ｇ）に示すようにマッチングした組Ｍ1、Ｍ2のうち最もコスト、すなわち、式（１）の値が低い組Ｍ1をマッチングした組として抽出する。
【０１５６】
以上は、説明を簡単にするために注目する位相特異点の数を３点とした場合について説明したが、３点以上の位相特異点に注目して同様な処理を行うことも可能である。
【０１５７】
次に、ステップＳ４−１１で算出される変換パラメータについて説明する。すなわち、位相特異点の移動が平行移動、回転、縮尺の３つのパラメータに分解できることについて説明する。
【０１５８】
一般に、流体力学の観点からどんな複雑な動きもそれを平行移動(transition)，回転(rotation)，縮尺(scale)で表すことができることが知られており、変位前後座標の回転関係式は次式で表すことができる。
【０１５９】
【数１２】
ここで、（ｘi、ｙi）は変位前の位相特異点の座標、（ｘi'、ｙi'）は変位後の位相特異点の座標、αは縮尺、θは回転角である。
【０１６０】
【数１３】
は変位前の位相特異点の重心、
【０１６１】
【数１４】
は変位後の位相特異点の重心を表している。
まず、はじめに、マッチングを表すのに適した関数について説明する。
【０１６２】
この関数はそれぞれ変位前後で対応している位相特異点の差を足し合わせた形で書ける。
【０１６３】
【数１５】
ここで、Ｃ＝αcosθ、Ｓ＝αsinθを表している。Ｅの値の最小値を見つけるためにはそれぞれのパラメータをまず計算しなければならない。例えば、
【０１６４】
【数１６】
などである。
【０１６５】
ここで、
【０１６６】
【数１７】
とおくと、
【０１６７】
【数１８】
から次式が得られる。
【０１６８】
【数１９】
これを解くと次のようになる。
【０１６９】
【数２０】
同様にして
【０１７０】
【数２１】
から次式が得られる。
【０１７１】
【数２２】
これを解くと次式が得られる。
【０１７２】
【数２３】
ＳとＣの定義より以下のパラメータが得られる。
【０１７３】
【数２４】
このように、位相特異点の移動は、平行移動Δｘ、Δｙ、回転θ、縮尺αの３つのパラメータで表せることがわかる。この３つのパラメータを変換パラメータと呼び、ステップＳ４−１１で算出される。
【０１７４】
〔実験・検証〕
次に、上記マッチング処理方法の実験・検証結果について説明する。
【０１７５】
図２８は実験・検証に用いられたＬＧフィルタの振幅特性図、図２９〜図３５は実験・検証結果を説明するための図を示す。図２９は実験画像を示す図であり、図２９（Ａ）は変位前の画像、図２９（Ｂ）は変位後の画像を示している。図３０は図２９に示す画像から抽出された位相特異点を示す画像であり、図２９（Ａ）は変位前の位相特異点の分布、図２９（Ｂ）は変位後の位相特異点の分布を示している。図３１は選択した位相特異点を拡大した図であり、図３１（Ａ）は変位前、図３１（Ｂ）は変位後の拡大図を示している。
【０１７６】
まず、今回の実験では、対象としてみどりふぐを使用し、高速ＣＣＤカメラを用いて秒間３０枚、連続撮影した画像を用いてマッチング処理を検証している。
【０１７７】
撮影した画像に対してＬＧフィルタ演算を行い、位相特異点を得る。このとき、ＬＧフィルタのフィルタ幅は、ω2＝３０としている。なお、ＬＧフィルタは、spiral phase filterにバンドパスフィルタの特性を加えたものであり、図２８に示すような振幅特性とした。なお、ＬＧフィルタの振幅幅ωを広くとれば、出てくる位相特異点の数は多くなり、逆に狭くすれば少なくなる。
【０１７８】
図２９に示す変位前後の画像に対して今回提案したマッチングアルゴリズムを適用した例について説明する。図２９に示す画像の変位前後の位相特異点の分布は図３０に示すようになる。図３０（Ａ）に示す変位前の位相特異点から黒丸で示される６つの位相特異点を選択してマッチングを行った。
【０１７９】
図３０（Ｂ）において、黒丸で示す位相特異点は変位前の黒丸で示す位相特異点に対応しており、実際にマッチングアルゴリズムにより検出されたパターンである。
【０１８０】
図３１は変位前に選んだ位相特異点の組み合わせと、ポイントマッチングアルゴリズムによって検出された変位後の位相特異点の組み合わせを拡大図示したものである。このように変位前，変位後のパターンを見比べると移動，回転，縮尺を組み合わせた変形が起こっているがマッチングアルゴリズムによって検証された。
【０１８１】
異なる２枚の画像に対して異なったパターンに着目して、本実施例のポイントマッチング処理を採用した場合の結果についても検証した。
【０１８２】
図３２は異なる２枚の画像に対して異なったパターンに着目して、本実施例のポイントマッチング処理を採用した場合の変位前後の位相特異点の分布を示す図であり、図３２（Ａ）は変位前の分布、図３２（Ｂ）は変位後の分布を示している。図３３は図３２に示した位相特異点の変位を示す図であり、図３３（Ｂ）は選択された位相特異点を拡大した図である。
【０１８３】
この検証例においても図３２（Ａ）に示すように変位前の４つの位相特異点が変位後、図３２（Ｂ）に示す位置に移動し、図３３に示すように注目した４つの位相特異点が変位していることが検証できた。
【０１８４】
さらに同じ作業をふぐの連続撮影画像１００枚に対して行った。
【０１８５】
図３４は連続した画像について実験・検証例を行った場合の画像を示す図である。図３４（Ａ）は最初の画像、図３４（Ｂ）は最後の画像を示している。図３５は図３４に示す画像から６つの位相特異点に注目し、その軌跡を追跡したものである。
【０１８６】
このように、本実施例のポイントマッチング処理を用いることにより、位相特異点を連続して追跡できる。
【０１８７】
さらに、任意の位相特異点の組み合わせに対して位相特異点が持つ位相構造と他の位相特異点との距離とによって対応する位相特異点の組を検出し、それらの位相特異点の組を追跡することによってランダムな動きを追跡できることがわかる。特徴パラメータが変化してしまう理由で計測不可能であった大きな変位を伴った物体の運動変位計測が可能となった。
【０１８８】
また、本実施例では、位相特異点の組のうち式（１）が最小となる組み合わせを検出することにより、マッチングを行ったが、位相特異点の点間総和距離の差を算出し、変位前のパターンと変位後のパターンとで点間総和距離の差が最小となる位相特異点を検出することによりマッチングを行うようにしてもよい。
【０１８９】
図３６、図３７は点間総和距離の差によりマッチングを行う動作を説明するための図を示す。
【０１９０】
図３６（Ａ）、図３７（Ａ）は変位前の位相特異点の分布、図３６（Ｂ）、図３７（Ｂ）は変位後の位相特異点の分布を示す。
【０１９１】
まず、変位前の画像から位相特異点ｐ1〜ｐ4を抽出するとともに、変位後の画像から位相特異点ｐ1'、ｐ2'、ｐ3'、ｐ4'を抽出し、変位前の位相特異点ｐ1〜ｐ4と変位後の位相特異点ｐ1'、ｐ2'、ｐ3'、ｐ4'との点間総和距離の差Ｓ1を求める。点間総和距離の差Ｓ1は、
Ｓ1＝（｜ｐ1ｐ3｜−｜ｐ1'ｐ3'｜）2＋（｜ｐ1ｐ2｜−｜ｐ1'ｐ2'｜）2
＋（｜ｐ2ｐ3｜−｜ｐ2'ｐ3'｜）2＋（｜ｐ3ｐ4｜−｜ｐ3'ｐ4'｜）2
＋（｜ｐ1ｐ4｜−｜ｐ1'ｐ4'｜）2＋（｜ｐ2ｐ4｜−｜ｐ2'ｐ4'｜）2
で求められる。
ここで、｜ｐiｐj｜は、点ｐiと点ｐjとの間の距離を表している。
【０１９２】
また、変位後の画像から位相特異点ｐ1'、ｐ2'、ｐ3'、ｐ4'とは異なる位相特異点ｐ1"、ｐ2"、ｐ3"、ｐ4"を抽出し、変位前の位相特異点ｐ1〜ｐ4と変位後の位相特異点ｐ1"、ｐ2"、ｐ3"、ｐ4"との点間総和距離の差Ｓ2を求める。点間総和距離の差Ｓ2は、
Ｓ2＝（｜ｐ1ｐ3｜−｜ｐ1"ｐ3"｜）2＋（｜ｐ1ｐ2｜−｜ｐ1"ｐ2"｜）2
＋（｜ｐ2ｐ3｜−｜ｐ2"ｐ3"｜）2＋（｜ｐ3ｐ4｜−｜ｐ3"ｐ4"｜）2
＋（｜ｐ1ｐ4｜−｜ｐ1"ｐ4"｜）2＋（｜ｐ2ｐ4｜−｜ｐ2"ｐ4"｜）2
で求められる。
【０１９３】
このとき、例えば、Ｓ2＜Ｓ1であれば、変位前の画像から位相特異点ｐ1〜ｐ4は、位相特異点ｐ1"、ｐ2"、ｐ3"、ｐ4"に対応すると判定する。
【０１９４】
なお、ここでは、説明を簡単にするために４点について説明したが、２点以上であれば、点間総和距離の差によりマッチングを行うことが可能であることは言うまでもない。
【０１９５】
本実施例のポイントマッチング処理は、複数の位相特異点の位置関係を用いるので縮尺、移動そして、回転により変位前後で完全にその位置関係が一致しなくてもポイントマッチングを行うことで一番近い組を選び出すことができるので、本発明はあらゆる変位を計測することができる。
【０１９６】
なお、本実施例では、特徴点として位置情報、離心率、渦度、実部がゼロとなる線分Ｌre0と虚部がゼロとなる線分Ｌim0との交叉角度であるゼロクロス角を特徴パラメータとして持つ位相特異点を例に説明したが、特徴パラメータはこれに限定されるものではなく、例えば、位置情報だけであってもよい。
【０１９７】
以上、本実施例によれば、干渉計なしで空間的なランダムパターンの強度情報を位相情報に対応づけ、位相情報の中に存在する全位相特異点の特徴を基に変位前後で各位相特異点がどこに移動したか特定することにより、計測物体の全体の変位計測に対しても、物体の局所領域における変位計測のいずれに対しても適用することができ、さらに場の回転の計測をも実現した変位検出装置を提供することできる。
【０１９８】
本国際出願は、２００６年２月１日に出願した日本国特許出願２００６−０２４８４６号、及び、２００６年１１月８日に出願した日本国特許出願２００６−３０３２３８号に基づく優先権を主張するものであり、日本国特許出願２００６−０２４８４６号、及び、日本国特許出願２００６−３０３２３８号の全内容を本国際出願に援用する。
【符号の説明】
【０１９９】
１０１ 変位検出装置、１０２ 測定対象
１１１ 画像検出部、１１２ 信号処理部
１２１ 光源、１２２ 撮像装置
１３１ フレームメモリ、１３２ 処理部、１３３ フーリエ変換部
１３４ ＲＯＭ
１３５ ハードディスクドライブ、１３６ ディスプレイ、１３７ 入力装置
Claims (51)
- 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを有し、
前記位相特異点特定手順で取得される特徴パラメータは、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度であることを特徴とする変位検出方法。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを有し、
前記位相特異点特定手順で取得される特徴パラメータは、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率であることを特徴とする変位検出方法。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを有し、
前記位相特異点特定手順で取得される特徴パラメータは、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積であることを特徴とする変位検出方法。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを有し、
前記位相特異点特定手順は、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度と、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率と、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積である渦度とを含むことを特徴とする変位検出方法。 - 前記変位検出手順は、前記位相特異点特定手順で特定された位相特異点の変位前後の位置を前記特徴パラメータに基づいて特定する変位前後位置検出手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置に基づいて前記位相特異点の変位を測定する変位測定手順とを有することを特徴とする請求項１乃至４のいずれか一項記載の変位検出方法。 - 前記変位測定手順は、前記変位前後位置検出手順で検出された変位前後の位相特異点の位置の座標の差分の二乗平方根を位相特異点の変位量として算出することを特徴とする請求項５記載の変位検出方法。
- 前記変位測定手順は、前記変位前後位置検出手順で検出された複数組の変位前後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる点を回転中心として求める回転中心取得手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置と前記回転中心取得手順で取得した前記回転中心とからなる三角形の前記回転中心を挟む辺の角度を回転角として求める回転角度取得手順とを有することを特徴とする請求項５記載の変位検出方法。 - 前記フィルタリング処理は、前記フーリエ変換手順でフーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振幅特性を持つラゲールガウスビームフィルタ特性による処理であることを特徴とする請求項１乃至４のいずれか一項記載の変位検出方法。
- 前記画像の強度信号の平均値を前記画像の各強度信号から減算した画像を前記画像とする前処理手順を有することを特徴とする請求項１乃至４のいずれか一項記載の変位検出方法。
- 前記フィルタリング処理は、前記フーリエ変換手順でフーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒルベルト変換による処理であることを特徴とする請求項１乃至４のいずれか一項記載の変位検出方法。
- 変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うマッチング手順を、更に、有することを特徴とする請求項１乃至４のいずれか一項記載の変位検出方法。
- 前記マッチング手順は、変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる特徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定することを特徴とする請求項１１記載の変位検出方法。
- 前記マッチング手順は、変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
- 取得した画像の強度信号をフーリエ変換するフーリエ変換手段と、
前記フーリエ変換手段でフーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行う手段と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手段と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手段と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手段と、
前記位相特異点特定手段で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手段とを有し、
前記位相特異点特定手段で取得される特徴パラメータは、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度であることを特徴とする変位検出装置。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手段と、
前記フーリエ変換手段でフーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行う手段と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手段と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手段と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手段と、
前記位相特異点特定手段で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手段とを有し、
前記位相特異点特定手段で取得される特徴パラメータは、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率であることを特徴とする変位検出装置。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手段と、
前記フーリエ変換手段でフーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行う手段と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手段と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手段と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手段と、
前記位相特異点特定手段で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手段とを有し、
前記位相特異点特定手段で取得される特徴パラメータは、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積であることを特徴とする変位検出装置。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手段と、
前記フーリエ変換手段でフーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理を行う手段と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手段と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手段と、
前記位相特異点の位相構造から所定の特徴パラメータを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手段と、
前記位相特異点特定手段で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手段とを有し、
前記位相特異点特定手段は、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度と、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率と、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積である渦度とを含むことを特徴とする変位検出装置。 - 前記変位検出手段は、前記位相特異点特定手段で特定された位相特異点の変位前後の位置を前記特徴パラメータに基づいて特定する変位前後位置検出手段と、
前記変位前後位置検出手段で検出された変位前後の位相特異点の位置に基づいて前記位相特異点の変位を測定する変位測定手段とを有することを特徴とする請求項１４乃至１７のいずれか一項記載の変位検出装置。 - 前記変位測定手段は、
前記変位前後位置検出手段で検出された変位前後の位相特異点の位置の座標の差の二乗平方根を位相特異点の変位量として算出することを特徴とする請求項１８記載の変位検出装置。 - 前記変位測定手段は、前記変位前後位置検出手段で検出された複数組の変位前後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる点を回転中心として求める回転中心取得手段と、
前記変位前後位置検出手段で検出された変位前後の位相特異点の位置と前記回転中心取得手段で取得した前記回転中心とからなる三角形の前記回転中心を挟む辺の角度を回転角として求める回転角度取得手段とを有することを特徴とする請求項１８記載の変位検出装置。 - 前記フィルタリング処理は、前記フーリエ変換手段でフーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振幅特性を持つラゲールガウスビームフィルタ特性による処理であることを特徴とする請求項１８記載の変位検出装置。
- 前記画像の強度信号の平均値を前記画像の各強度信号から減算した画像を前記画像とする前処理手段を有することを特徴とする請求項１４乃至１７のいずれか一項記載の変位検出装置。
- 前記フィルタリング処理は、前記フーリエ変換手段でフーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒルベルト変換による処理であることを特徴とする請求項１４乃至１７のいずれか一項記載の変位検出装置。
- 変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うマッチング手段を、更に、有することを特徴とする請求項１４乃至１７のいずれか一項記載の変位検出装置。
- 前記マッチング手段は、変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる特徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定することを特徴とする請求項２４記載の変位検出装置。
- 前記マッチング手段は、変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
- コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを実行させる変位検出プログラム。 - コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを実行させる変位検出プログラム。 - コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを実行させる変位検出プログラム。 - コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度と、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率と、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積である渦度とを取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出することにより変位を検出する変位検出手順とを実行させる変位検出プログラム。 - 前記変位検出手順は、前記位相特異点特定手順で特定された位相特異点の変位前後の位置を前記特徴パラメータに基づいて特定する変位前後位置検出手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置に基づいて前記位相特異点の変位を測定する変位測定手順とを有することを特徴とする請求項２７乃至３０のいずれか一項記載の変位検出プログラム。 - 前記変位測定手順は、前記変位前後位置検出手順で検出された変位前後の位相特異点の位置座標の差分の二乗平方根を位相特異点の変位量として算出することを特徴とする請求項３１記載の変位検出プログラム。
- 前記変位測定手順は、前記変位前後位置検出手順で検出された複数組の変位前後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる点を回転中心として求める回転中心取得手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置と前記回転中心取得手順で取得した前記回転中心とからなる三角形の前記回転中心を挟む辺の角度を回転角として求める回転角度取得手順とを有することを特徴とする請求項３１記載の変位検出プログラム。 - 前記フィルタリング処理は、前記フーリエ変換手順でフーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振幅特性を持つラゲールガウスビームフィルタ特性による処理であることを特徴とする請求項２７乃至３０のいずれか一項記載の変位検出プログラム。
- 前記画像の強度信号の平均値を前記画像の各強度信号から減算した画像を前記画像とする前処理手順を有することを特徴とする請求項２７乃至３０のいずれか一項記載の変位検出プログラム。
- 前記フィルタリング処理は、前記フーリエ変換手順でフーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒルベルト変換による処理であることを特徴とする請求項２７乃至３０のいずれか一項記載の変位検出プログラム。
- 変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うマッチング手順を、更に、有することを特徴とする請求項２７乃至３０のいずれか一項記載の変位検出プログラム。
- 前記マッチング手順は、変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる特徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定することを特徴とする請求項３７記載の変位検出プログラム。
- 前記マッチング手順は、変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
- 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とする位相特異点マッチング処理方法。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とする位相特異点マッチング処理方法。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とする位相特異点マッチング処理方法。 - 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、位置情報、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積である渦度、及び前記位相特異点近傍の前記複素解析信号の実部がゼロとなる線分と虚部がゼロとなる線分とが交わる角度を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とする位相特異点マッチング処理方法。 - 変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる位相特異点の組を変位前と変位後とで互いにマッチングする位相特異点として特定することを特徴とする請求項４０乃至４３のいずれか一項記載の位相特異点マッチング処理方法。
- 変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
- コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部と虚部の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とするコンピュータ読み取り可能な位相特異点マッチング処理プログラム。 - コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とするコンピュータ読み取り可能な位相特異点マッチング処理プログラム。 - コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とするコンピュータ読み取り可能な位相特異点マッチング処理プログラム。 - コンピュータに、
取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号のエルミート対称性を崩すようなフィルタリング処理を行う手順と、
前記フィルタリング処理の後にフーリエ逆変換することにより複素解析信号を得る手順と、
前記複素解析信号の実部と虚部の値が共にゼロとなる位相特異点を取得する位相特異点取得手順と、
前記位相特異点の位相構造から特徴パラメータとして、位置情報、前記位相特異点近傍の前記複素解析信号の強度分布の形状の離心率、前記位相特異点近傍の前記複素解析信号の実部の勾配と前記複素解析信号の虚部の勾配とのベクトル積である渦度、及び前記位相特異点近傍の前記複素解析信号の実部がゼロとなる線分と虚部がゼロとなる線分とが交わる角度を取得し、前記特徴パラメータに基づいて前記位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された特徴パラメータに基づいて変位後の位相特異点の位置を検出する手順とを有し、
変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの差異に基づいて変位前の位相特異点に対応する変位後の位相特異点を探しだし、その位相特異点の地理的な位置関係の類似性から変位前の位相特異点に対応する変位後の位相特異点のマッチングを行うことを特徴とするコンピュータ読み取り可能な位相特異点マッチング処理プログラム。 - 変位前の複数の位相特異点間の特徴パラメータの総和距離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小となる位相特異点の組を変位前と変位後とで互いにマッチングする位相特異点として特定することを特徴とする請求項４６乃至４９のいずれか一項記載の位相特異点マッチング処理プログラム。
- 変位前の所定の位相特異点の特徴パラメータをＴi1、Ｔi2、変位後の所定の位相特異点の特徴パラメータをＰw(i1)、Ｐw(i2)、αを縮尺としたとき、
Applications Claiming Priority (5)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2006024846 | 2006-02-01 | ||
| JP2006024846 | 2006-02-01 | ||
| JP2006303238 | 2006-11-08 | ||
| JP2006303238 | 2006-11-08 | ||
| PCT/JP2007/051095 WO2007088759A1 (ja) | 2006-02-01 | 2007-01-24 | 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラム |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPWO2007088759A1 JPWO2007088759A1 (ja) | 2009-06-25 |
| JP4385139B2 true JP4385139B2 (ja) | 2009-12-16 |
Family
ID=38327338
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2007556824A Active JP4385139B2 (ja) | 2006-02-01 | 2007-01-24 | 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラム |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US7813900B2 (ja) |
| JP (1) | JP4385139B2 (ja) |
| WO (1) | WO2007088759A1 (ja) |
Families Citing this family (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8116553B2 (en) * | 2007-10-03 | 2012-02-14 | Siemens Product Lifecycle Management Software Inc. | Rotation invariant 2D sketch descriptor |
| JP5224441B2 (ja) * | 2007-12-14 | 2013-07-03 | 独立行政法人産業技術総合研究所 | 高さを測定する方法及び高さ測定装置 |
| JP2010210571A (ja) * | 2009-03-12 | 2010-09-24 | Mitsutoyo Corp | 画像相関変位計、及び変位測定方法 |
| US8730332B2 (en) | 2010-09-29 | 2014-05-20 | Digitaloptics Corporation | Systems and methods for ergonomic measurement |
| US8917950B2 (en) * | 2011-01-18 | 2014-12-23 | Sony Corporation | Simplifying parametric loop filters |
| US8587666B2 (en) * | 2011-02-15 | 2013-11-19 | DigitalOptics Corporation Europe Limited | Object detection from image profiles within sequences of acquired digital images |
| US8705894B2 (en) | 2011-02-15 | 2014-04-22 | Digital Optics Corporation Europe Limited | Image rotation from local motion estimates |
| US9285206B1 (en) * | 2012-02-07 | 2016-03-15 | Pile Dynamics, Inc. | Measurement device for pile displacement and method for use of the same |
| JP6516134B2 (ja) * | 2013-06-29 | 2019-05-22 | 堀 健治 | 位相変換作用を持つフィルター、レンズ、結像光学系及び撮像システム |
| JP2016009043A (ja) * | 2014-06-24 | 2016-01-18 | ソニー株式会社 | イメージセンサ、演算方法、および電子装置 |
| KR102107145B1 (ko) * | 2018-04-25 | 2020-05-06 | 경북대학교 산학협력단 | 구조물 변형 측정 장치 및 방법, 기록 매체 |
| CN109883333A (zh) * | 2019-03-14 | 2019-06-14 | 武汉理工大学 | 一种基于图像特征识别技术的非接触位移应变测量方法 |
| CN112926356B (zh) * | 2019-12-05 | 2024-06-18 | 北京沃东天骏信息技术有限公司 | 一种目标跟踪方法和装置 |
| CN112614180B (zh) * | 2020-12-31 | 2022-05-31 | 中国科学院长春光学精密机械与物理研究所 | 一种位移测量方法 |
Family Cites Families (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5118935A (en) * | 1991-02-28 | 1992-06-02 | Eastman Kodak Company | Composite imaging mask |
| DE69631845T2 (de) | 1996-11-01 | 2005-01-05 | Yamatake Corp. | Musterextrahierungsgerät |
| US6577966B2 (en) * | 2000-06-21 | 2003-06-10 | Siemens Corporate Research, Inc. | Optimal ratio estimator for multisensor systems |
| JP2002171395A (ja) | 2000-11-29 | 2002-06-14 | Canon Inc | 画像処理装置及びその方法並びに記憶媒体 |
| US20020076116A1 (en) * | 2000-12-15 | 2002-06-20 | Xerox Corporation | Fast implementation of homomorphic filters for image enhancement |
| JP2002190020A (ja) | 2000-12-20 | 2002-07-05 | Monolith Co Ltd | 映像効果方法および装置 |
| JP3871309B2 (ja) * | 2001-01-31 | 2007-01-24 | フジノン株式会社 | 位相シフト縞解析方法およびこれを用いた装置 |
| FR2863080B1 (fr) * | 2003-11-27 | 2006-02-24 | Advestigo | Procede d'indexation et d'identification de documents multimedias |
| TW200604494A (en) * | 2004-04-22 | 2006-02-01 | Univ Electro Communications | Small displacement measuring method and instrument |
| JP2006038779A (ja) * | 2004-07-30 | 2006-02-09 | Hitachi High-Technologies Corp | パターン形状評価方法、評価装置、及び半導体装置の製造方法 |
| US7539354B2 (en) * | 2004-08-25 | 2009-05-26 | Canon Kabushiki Kaisha | Image database key generation method |
-
2007
- 2007-01-24 WO PCT/JP2007/051095 patent/WO2007088759A1/ja active Search and Examination
- 2007-01-24 US US12/162,370 patent/US7813900B2/en active Active
- 2007-01-24 JP JP2007556824A patent/JP4385139B2/ja active Active
Also Published As
| Publication number | Publication date |
|---|---|
| JPWO2007088759A1 (ja) | 2009-06-25 |
| WO2007088759A1 (ja) | 2007-08-09 |
| US7813900B2 (en) | 2010-10-12 |
| US20090037137A1 (en) | 2009-02-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4385139B2 (ja) | 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラム | |
| Shang et al. | Multi-point vibration measurement and mode magnification of civil structures using video-based motion processing | |
| Radkowski | Object tracking with a range camera for augmented reality assembly assistance | |
| Liu et al. | Structural motion estimation via Hilbert transform enhanced phase-based video processing | |
| Cicconet et al. | Finding mirror symmetry via registration and optimal symmetric pairwise assignment of curves: Algorithm and results | |
| Podsiadlo et al. | Characterization of surface topography of wear particles by SEM stereoscopy | |
| Luan et al. | Extracting full-field subpixel structural displacements from videos via deep learning | |
| Javed et al. | Vibration measurement of a rotating cylindrical structure using subpixel-based edge detection and edge tracking | |
| Réthoré et al. | Curve and boundaries measurement using B-splines and virtual images | |
| Wang et al. | Target-less approach of vibration measurement with virtual points constructed with cross ratios | |
| JP2006214893A (ja) | 対象物測定方法及びコンピュータシステムを用いて対象物の三次元形状を計測するためのコンピュータソフトウエアプログラム | |
| Khokhlov et al. | Multi-scale stereo-photogrammetry system for fractographic analysis using scanning electron microscopy | |
| Li et al. | Two-dimensional motion estimation using phase-based image processing with Riesz transform | |
| Lv et al. | A point tracking method of TDDM for vibration measurement and large-scale rotational motion tracking | |
| Coleman et al. | Gradient operators for feature extraction and characterisation in range images | |
| Cosco et al. | Towards phase-based defect detection: A feasibility study in vibrating panels | |
| Shang et al. | Multi-point vibration measurement for mode identification of bridge structures using video-based motion magnification | |
| Li et al. | A novel multispectral fusion defect detection framework with coarse-to-fine multispectral registration | |
| Prasad et al. | A real-time feature-based clustering approach for vibration-based SHM of large structures | |
| Hang et al. | Eulerian fast motion identification algorithm for deformation measurement of cable-stayed bridge | |
| JP2015206654A (ja) | 情報処理装置、情報処理方法、プログラム | |
| Miyazaki et al. | Dynamic projection mapping onto non-rigid objects with dot markers | |
| JP4820998B2 (ja) | 位相特異点検出方法、及び、位相特異点検出装置、並びに、プログラム | |
| KR101039183B1 (ko) | 보간이 필요없는 어파인변환 기반 하이브리드 피아이브이를 이용한 유동장 계측 시스템 | |
| Wang et al. | Full Period Three-dimensional (3-D) Reconstruction Method for a Low Cost Singlelayer Lidar. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20090507 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20090901 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20090907 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |