JPH09161067A

JPH09161067A - 画像マッチング方法

Info

Publication number: JPH09161067A
Application number: JP7345901A
Authority: JP
Inventors: Koichi Ejiri; 公一江尻; Hideaki Yamagata; 秀明山形
Original assignee: Ricoh Co Ltd
Current assignee: Ricoh Co Ltd
Priority date: 1995-12-08
Filing date: 1995-12-08
Publication date: 1997-06-20

Abstract

(57)【要約】【課題】２枚の画像の対応点を精度よく見つける。【解決手段】２枚の画像の走査範囲を設定し（１）、
その中に走査窓を設定する（２）。窓を移動させて特徴
パラメータを計算する（３）。特徴パラメータを平滑化
した（４）後、極小値を探索し（５）、極小座標値を算
出する（６）。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、２枚の画像の対応
点を見つける、画像マッチング方法に関する。

【０００２】

【従来の技術】同一の対象物を撮影した２枚の画像が与
えられたとき、２枚の画像から特定の対象物を計測した
り、あるいは２枚の画像を重ね合わせる必要が生じるこ
とがある。代表的な用途の一つは大型の図面を複写する
場合である。このとき、図面全体は一度に複写できない
から、部分的にオーバーラップさせながら読み取り、後
でオーバーラップ部分を基にして個々の小さい画像を重
ね合わせる。しかし、画素のレベルまで誤差なく正確に
重ねることは人手では困難である。

【０００３】

【発明が解決しようとする課題】従来、上記したような
必要が生じたとき、オーバーラップ部分で目立ちやすい
対象物を人が見つけ、これをコンピュータに指示して重
ね合わせるようにしていた。あるいは、これをコンピュ
ータで自動的に行う場合には、適当に定めた小領域を、
画像全体にわたって一画素ずつずらしながら最大マッチ
ング部分を探し、しかる後にこの部分を重ねていた。

【０００４】小領域に属する画素の数をＳ、探索する画
像領域の画素の数をＮとすると、比較に必要とする演算
回数は、Ｓ×Ｎとなる。小領域は１つでは足りず、最低
３つは必要である。従って、演算回数は少なくてもＳ×
Ｎ×３となる。なお、この種の関連する技術としては、
例えば特開平５−１８９０２号公報に記載された画像間
の位置合わせ方法がある。

【０００５】一方、白黒の２値画像に限られるが、次の
ような画像の重ね合わせ方法が知られている。まず、２
つの画像を重ね合わせる。重ねた画像中には、しばしば
同心円状のパターンが現れる。この同心円の中心が２つ
の画像に対応する部分（つまり同じパターン）である。

【０００６】この処理に要する演算回数は、重ね合わせ
のために２Ｎとなる。また、円状パターンの中心を見つ
けるには、小領域（画素数をＳとする）を定義し、その
小領域内部のランの数または周囲長などを計測する。小
領域の移動ウィンドウは単純なマッチングと異なり、通
常１００分の１くらいの割合で大まかに行う。従って、
演算回数は、２Ｎ＋Ｓ×Ｎ×０．０１≒２Ｎとなり、先
の方法より演算回数は少なくなる。ただし、Ｓ≒１００
〜１０００なので、２項目は無視した。

【０００７】しかし、後者の方法には、以下のような欠
点がある。すなわち、（１）もし、２つの画像の方向がたまたま揃っていると
（つまり回転していない）、同心円のパターンが現れな
い。（２）もし、２つの画像が白黒の２値画像でないなら
ば、同心円のパターンは現れない。

【０００８】本発明は、上記した事情を考慮してなされ
たもので、本発明の目的は、２枚の画像の対応点を精度
よく見つけることができる画像マッチング方法を提供す
ることにある。

【０００９】

【課題を解決するための手段】前記目的を達成するため
に、請求項１記載の発明では、２枚の画像の対応点を求
める画像マッチング方法において、該２枚の画像の内、
一方の画像を他方の画像に対して縮小または拡大させ、
しかる後に両画像を重ね合わせて、その相関をとること
によって前記対応点を求めることを特徴としている。

【００１０】請求項２記載の発明では、前記相関は、所
定の小領域内部で、前記両画像の差分の和をとる演算で
あることを特徴としている。

【００１１】請求項３記載の発明では、前記重ね合わせ
た画像の部分領域において、２次元の周波数分布を求
め、該分布関数から所定の方向成分を抽出することを特
徴としている。

【００１２】請求項４記載の発明では、前記周波数分布
を求めるとき、元の画像の周波数成分を差し引いた残り
から、周波数成分を算出することを特徴としている。

【００１３】請求項５記載の発明では、前記周波数分布
を求めるとき、元の２つの画像のそれぞれの周波数成
分、両者の差分成分、および両者の加算成分を差し引い
た残りから、周波数成分を算出することを特徴としてい
る。

【００１４】請求項６記載の発明では、前記縮小率また
は拡大率を１に近い値に設定するとともに、回転角度を
０度に近い値に設定することを特徴としている。

【００１５】

【発明の実施の形態】以下、本発明の一実施例を図面を
用いて具体的に説明する。同一の対象物を撮影した２枚
の画像中から同一のパターンを見つける方法として、従
来、パターンマッチング法が使用されていて、多くの場
合、有効な方法である。しかし、より高速性を求める場
合や、ランダムパターンのように特徴のない画像を扱う
場合は、この方法は適さない。以下、説明する本発明の
方法は、従来困難であった、ランダムパターンに適用で
きるばかりでなく、さらに中間調のランダムパターンや
カラーのランダムパターンにも適用することができる点
が特徴である。また、本発明は、デジタル画像、すなわ
ち画素の単位で読み取られた画像に適用される。

【００１６】〈実施例１〉まず、対象となる画像をｆ
（ｘ，ｙ）、ｇ（ｘ’，ｙ’）で表す。以下の説明ではｆ（ｘ，ｙ）≒ｇ（ｘ’，ｙ’）（１）とする。ここで、（ｘ，ｙ）、（ｘ’，ｙ’）は、各画
素の座標値である。あるいはｇ＝Ａｆ

【００１７】

【数１】ここで扱う２つの画像は、歪が互いに小さいと仮定し
（式（１）の≒の意味）、ａ≒ｄ≒１ｂ＋ｃ＝０かつｂ≒０，ｃ≒０ｅ、ｆは任意であるとする。

【００１８】ここでは、画像のズームと回転、平行移動
のみを考慮することにすると（通常の光学系ではこれで
十分）、式（２）は、

【００１９】

【数２】となる。ここで、ｓは回転角度、Ｐ１はズーム倍率であ
る。従って、式（２）の下に表された条件は、｜ｓ｜≒０、Ｐ１≒１（４）と書き表せる。

【００２０】図１は、本発明の実施例１の機能ブロック
図である。走査範囲設定部１は、２つの画像ｆ（ｘ，
ｙ）、ｇ（ｘ’，ｙ’）を用意し、共通の座標系（ｉ，
ｊ）を定義し、走査範囲Ｒを設定する。ここで（ｉ，
ｊ）の座標系は、通常１枚目の画像の座標系を使う。ま
た、２つの画像ｆとｇは、同一の対象物を含み、かつ共
通座標系上で同一対象物が重なる部分を持つものとす
る。

【００２１】走査窓設定部２は、画像ｆ（ｉ，ｊ）、ｇ
（ｉ，ｊ）に共通する小さな走査窓Ｗ（ｋ，ｌ）を設定
する。すなわち、Ｗ（ｋ，ｌ）は、範囲Ｒ内部を走査す
るものとし、その位置関係は図２に示される。

【００２２】特徴パラメータ計算部３は、窓ＷをＲ内部
で走査しながら、例えば以下のような特徴パラメータを
計算する。走査ステップは、Ｗの一辺の半分程度の距離
とする（Ｗの一辺が１６画素ならば、走査ステップは８
程度となる）。

【００２３】特徴パラメータＰｗは、Ｐｗ＝Σ｜ｆ（ｉ，ｊ）−ｇ（ｉ，ｊ）｜／ＤｓｉｆΣ（）＞０Ｄｓ＝Σ（｜ｄｆ／ｄｉ｜＋｜ｄｆ／ｄｊ｜＋｜ｄｇ／
ｄｉ｜＋｜ｄｇ／ｄｊ｜ｉｆΣ（）＝０Ｄｓ＝十分大きな値ここで、加算領域は、窓Ｗの内部とする。カウントした
Ｐｗの値を（ｉ，ｊ）平面上にプロットする。図３は、
特徴パラメータの値をプロットした図である。また、
（ｉ、ｊ）座標系を画像ｆに合わせると、一般に画像ｇ
の画素位置は整数値にならないが、ここの計算では最も
近接する整数画素値を利用するものとする。

【００２４】平滑化処理部４は、特徴パラメータ計算部
３でプロットしたＰｗの値を平滑化する。これを〈Ｐ
ｗ〉（ｉ，ｊ）とおく。

【００２５】極小値探索部５は、〈Ｐｗ〉の最も鋭い極
小値を探す。具体的には、〈Ｐｗ〉（ｉ，ｊ）の２次微
分の絶対値が大きく周囲の平均値より小さな極小値を持
つ座標値（Ｉｍｉｎ，Ｊｍｉｎ）を求める。ここで
（Ｉ，Ｊ）はＰｗ値の存在する座標値のみを指すものと
する。すなわち、図３の例では、Ｐｗ（１，１）＝２
１、Ｐｗ（２，１）＝１８であり、Ｉｍｉｎ＝２４、Ｊ
ｍｉｎ＝１６である。

【００２６】上記した特徴パラメータ計算部３、平滑化
処理部４、極小値探索部５での演算は、図４に示すよう
に重なり合った、あるいは互いに比較された画像上の
「渦」の中心点、つまり両画像の対応点を求めることに
相当する。

【００２７】極小座標値算出部６は、より正確に極小点
を求めるために、極小点の周りの値から新たな極小点を
算出する。極小値探索部５で求めた極小点の周りの頻度
分布を２次関数で近似すると、極小値を持つ座標を（Ｘ
ｍｉｎ，Ｙｍｉｎ）で表し、その分布関数のｘ軸に沿っ
た値をＨ１（ｘ１，ｙ）、Ｈ２（ｘ２，ｙ）、Ｈ３（ｘ
３，ｙ）と表し、ｙ軸に沿った値をＨ１（ｘ，ｙ１）、
Ｈ２（ｘ，ｙ２）、Ｈ３（ｘ，ｙ３）と表すとき、Ｘｍｉｎ＝｛（Ｈ１−Ｈ２）（ｘ２−ｘ３）（ｘ２＋ｘ
３）−（Ｈ２−Ｈ３）（ｘ１−ｘ２）（ｘ１＋ｘ２）｝
／［２｛（Ｈ１−Ｈ２）（ｘ２−ｘ３）−（Ｈ２−Ｈ
３）（ｘ１−ｘ２）｝］Ｙｍｉｎ＝｛（Ｈ１−Ｈ２）（ｙ２−ｙ３）（ｙ２＋ｙ
３）−（Ｈ２−Ｈ３）（ｙ１−ｙ２）（ｙ１＋ｙ２）｝
／［２｛（Ｈ１−Ｈ２）（ｙ２−ｙ３）−（Ｈ２−Ｈ
３）（ｙ１−ｙ２）｝］によって極小座標値が求まる。

【００２８】求めた座標値をそれぞれの座標系に戻して
（Ｘｍｉｎ１，Ｙｍｉｎ１），（Ｘｍｉｎ１’，Ｙｍｉ
ｎ１’）とおく。ここで、前者は元の座標値、後者は参
照画像（図２の例では画像ｇ）の座標値を示す。この２
点を重ねれば画像が重なることになる。

【００２９】なお、式（３）のＰ１、ｓ、ｅまたはｆを
変化させ、上記した平滑化処理部４、極小値探索部５、
極小座標点算出部６での演算を実行する。このとき画像
は以前と異なった位置関係にあるから、例えば図５のよ
うになる（一方の画像を少し移動したもので、図４に比
べて渦の中心が大きく移動している）。

【００３０】例えば、もしＰ１≠１ならば、Ｅｏだけ画
像ｆを平行移動したとすると、ｅ’＝ｅ＋Ｅｏこのとき、極小値探索部５で求めた座標値を（Ｘｍｉｎ
２，Ｙｍｉｎ２），（Ｘｍｉｎ２’，Ｙｍｉｎ２’）と
おくとき、Ｘｍｉｎ１’→Ｘｍｉｎ１Ｙｍｉｎ１’→Ｙｍｉｎ１Ｘｍｉｎ２’→Ｘｍｉｎ２Ｙｍｉｎ２’→Ｙｍｉｎ２への変換を行えばよい。すなわち、式（３）において、
Ｐ１，ｓ，ｅ，ｆを変数として連立方程式を解き、Ｐ
１，ｓ，ｅ，ｆを求める。

【００３１】上記した方法による演算回数は、従来の最
短処理法である「同心円を使う方法」に比べて若干多
い。つまり、アフィン変換が加わるので約３Ｎとなる。

【００３２】次に実際の画像を例にして説明する。図４
は、２つのランダムなパターンの画像の内、一方の画像
を２％縮小し、約６度回転して互いに重ね合わせたもの
である。渦巻き状のパターンが形成され、渦巻き状のパ
ターンの中心が両画像の対応する点である。その中心で
の白黒パターンは、白の割合が大きい、言い換えれば２
次微分の絶対値が小さな値を持っている。図５は、図４
の重なり合ったパターンの一方を微少量移動したもの
で、渦巻パターンが大きく移動していることが分かる。

【００３３】渦巻の中心の移動量をＳｘ、Ｓｙで表し、
２つの画像の相対的な平行移動量（回転を含んでいても
よい）をΔｘ、Δｙとおくとき、両画像間の縮率の違い
Ｒｘ、Ｒｙは、Ｒｘ＝Ｓｘ／（Ｓｘ−Δｘ）Ｒｙ＝Ｓｙ／（Ｓｙ−Δｙ）となる。

【００３４】図６は、一般の文書画像に適用した例であ
り、渦巻の中心部分に、白／黒の変化の小さい領域が容
易に見える。また、図７は、２枚の画像の回転による変
位が加わらない例である。従来、図７のような縮率が同
じで、方位も同じ画像を重ねると、系統的なパターンが
観察できなかった。従って、両画像の対応する場所を見
つけることができないという問題があった。また、図８
は、従来知られている、同一のパターンを微少量回転し
て重ね合わせたもので、同心円のパターンが見えるが、
対応点が見えない。

【００３５】これに対して、本発明では、上記したよう
に縮率が異なれば明瞭な渦巻状パターンや放射状のパタ
ーンが観察され、その中心では２次微分値の絶対値が小
さくなっていて、両画像の対応点を精度よく求めること
ができる。

【００３６】〈実施例２〉上記した実施例１では、２つ
の画像の差分の総和Ｐｗを利用して、対応する領域を見
つけた。本実施例２では、これと異なるアルゴリズムで
対応するパターンを見つける。

【００３７】式（３）で定義されるように、２つの画像
は互いに変換された関係にある。また、画像ｆと画像ｇ
を重ねたものは、図４の例に示すように特有の渦巻パタ
ーンを生じる。ただし、この方法は文書画像のような高
い周波数成分を有する白黒画像についてのみ有効であ
る。

【００３８】重ね合わせた２つの白黒画像をｈ（ｘ，
ｙ）で表すと、ｈ（ｘ，ｙ）＝（１−ｇ（ｘ，ｙ））（１−ｆ（ｘ，
ｙ））となる。ここで、画像は濃度値をパラメータにして表さ
れている。

【００３９】画像ｆとｇは、式（３）の関係が成り立つ
とする。まず、画像ｈの渦巻の接線方向ベクトルを検出
する。このためには、画像ｈ内部に小領域Ｔを定義し、
その内部で２次元フーリエ変換を行う。小領域Ｔは実施
例１の小領域Ｒより一般に大きく、その一辺のサイズは
１００画素程度である。

【００４０】Ｔ内部のフーリエ変換結果を２次元表示す
れば、その中のピーク値の方向が渦の接線方向を表す。
小領域Ｔを動かしながら、各位置における渦接線がこの
ようにして得られる。図９は、画像の渦巻の接線方向ベ
クトルを示す。

【００４１】領域Ｔの中心座標値（Ｘｉ，Ｙｉ）は変換
式（３）によって（ｘ’ｉ，ｙ’ｉ）に移動しているは
ずであるから、変形量が小さいときは、すなわち（４）
の条件が成立するときは、図９の各位置での渦巻の勾配
Ｄｉは、（Ｙｉ−ｙ’ｉ）／（Ｘｉ−ｘ’ｉ）＝Ｄｉ（５）の関係を満たす。

【００４２】さて、回転角度ｓは小さいと仮定し、(ｘ
ｏ，ｙｏ）を回転中心として式（３)を書き下すと、ｘ’ｉ＝Ｐ１（ｘｉ−ｘｏ）−Ｐ１・ｓ（ｙｉ−ｙｏ）（６）ｙ’ｉ＝Ｐ１・ｓ（ｘｉ−ｘｏ）＋Ｐ１（ｙｉ−ｙｏ）（７）式（５）、（６）、（７）より、Ｐ１（ｓ−Ｄｉ）ｘｏ＋Ｐ１（Ｄｉ・ｓ＋１）ｙｏ＋（Ｄｉ・Ｐ１−Ｄｉ− Ｐ１・ｓ）Ｘｉ−（Ｄｉ・ｓ＋Ｐ１−１）Ｙｉ＝０（８）従って、多数の小領域Ｔについて、観測値（Ｘｉ，Ｙ
ｉ）、Ｄｉを求め、式（８）に代入して連立方程式が得
られる。これから、未知変数Ｐ１，ｓ，ｘｏ，ｙｏにつ
いて解く。その結果、縮率Ｐ１、回転角ｓ、渦パターン
の中心座標、すなわち対応パターンの座標値（ｘｏ，ｙ
ｏ）が求まる。なお、観測値の精度が高くないから、式
（８）の解はできるだけ多くの実測値で求め、その結果
を平均化するか、中央値のようなノイズに強いサンプリ
ング方法を採用する必要がある。式（８）の解法は多数
知られているので、その説明を省略する。

【００４３】この実施例２の方法による演算回数は、元
の画像と縮小画像（アフィン変換）重ね合わせに３Ｎ、
渦の接線方位算出に０．００１ＮＳ程度、方程式（３）
の解を求めるのをほとんど無視すると、従来法より若干
増加する。

【００４４】〈実施例３〉画像を濃度関数ｆ（ｘ），ｇ
（ｙ）で表現するとき、２つの画像を重ね合わせたとき
に見える画像は（１−ｆ（ｘ））（１−ｇ（ｙ））で表
わされる。これを分解すると、１−（ｆ（ｘ）＋ｇ（ｙ））＋ｆ（ｘ）ｇ（ｙ）となる。もし、ｆ（ｘ），ｇ（ｙ）が三角関数ならば、
上記式は周波数成分、Ａ＋Ｂ，Ａ−Ｂ，Ａ＋Ｂを持つこ
とはよく知られている。従って、以上の５成分（Ａ，
Ｂ，Ａ＋Ｂ，Ａ−Ｂ，Ａ＋Ｂ）を予め除去しておけば、
より鋭敏な周波数検出が可能である。

【００４５】２つの対応する画像を重ね合わせるか、あ
るいは差分をとり、その２次元周波数成分を求める。こ
のとき求められた周波数成分から、元の画像の周波数成
分を引き去った残りから、局在するピークを見つけて周
波数を抽出すると、より鋭敏な周波数を抽出することが
できる。

【００４６】図１０は、元の画像１、２を示し、画像１
と画像２は互いに少し回転している。図１１は、実施例
３の処理によって作成された画像を示す。画像３は、画
像１と画像２を重ね合わせ、差分をとった後にフーリエ
変換したものから、元の画像のフーリエ成分を差し引い
た画像である。画像４は、画像１と画像２を重ねて、フ
ーリエ変換した画像である。画像３では、「黒い筋」
（これはフーリエ成分が極小値をとる場所）が、画像１
と画像２の対応する点に向かって揃っていることが分か
る。これから、両画像の対応する点が容易に求められ
る。

【００４７】演算の対象となる画素の数の参照回数は、
小領域内部の画素数をｎとすると、ＦＦＴを使った場
合、ｎｌｏｇ（ｎ）、２次ＦＦＴの場合は２ｎｌｎ
（ｎ）、小領域の数はＮ／ｎ存在するから、２Ｎｌｎ
（ｎ）と表わせる。

【００４８】〈実施例４〉画像を濃度関数ｆ（ｘ），ｇ
（ｙ）で表現するとき、２つの画像が白黒の２値画像で
ある場合は重ねたときに見える画像には、しばしば２つ
の画像の移動量に対応する周波数成分が弱いことがあ
る。これは、黒部分で輝度信号が弱くなり、重ねた画像
中から信号の一部が消失するためである。

【００４９】これを防ぐために、｜ｆ（ｘ）−ｇ（ｙ）｜に対して２次元周波数成分を求めると効果的である。こ
のとき求められた周波数成分から、元の画像の周波数成
分を引き去った残りから、局在するピークを見つけて周
波数を抽出すると、より鋭敏な周波数成分を抽出するこ
とができる。

【００５０】

【発明の効果】以上、説明したように、本発明によれ
ば、同一の対象物を撮影した２枚の画像中から同一のパ
ターンを精度よく求めることができる。また、ランダム
パターン、さらに中間調のランダムパターンやカラーの
ランダムパターンのように特徴のない画像にも適用する
ことができる。

【図面の簡単な説明】

【図１】本発明の実施例１の機能ブロック図である。

【図２】走査窓を用いて画像を走査する図である。

【図３】特徴パラメータの値をプロットした図である。

【図４】２枚のランダムパターンの画像を重ねた例を示
す。

【図５】図４の一方の画像を微少位置だけ移動したとき
の例を示す。

【図６】一般の文書画像を重ねた例を示す。

【図７】２枚の画像の回転による変位が加わらない画像
を示す。

【図８】同一のパターンを微少量回転して重ね合わせた
画像を示す。

【図９】画像の渦巻の接線方向ベクトルを示す。

【図１０】元の画像１、２を示す。

【図１１】実施例３の処理によって作成された画像３、
４を示す。

【符号の説明】

１走査範囲設定部２走査窓設定部３特徴パラメータ計算部４平滑化処理部５極小値探索部６極小座標値算出部

Claims

【特許請求の範囲】

【請求項１】２枚の画像の対応点を求める画像マッチ
ング方法において、該２枚の画像の内、一方の画像を他
方の画像に対して縮小または拡大させ、しかる後に両画
像を重ね合わせて、その相関をとることによって前記対
応点を求めることを特徴とする画像マッチング方法。
【請求項２】前記相関は、所定の小領域内部で、前記
両画像の差分の和をとる演算であることを特徴とする請
求項１記載の画像マッチング方法。
【請求項３】前記重ね合わせた画像の部分領域におい
て、２次元の周波数分布を求め、該分布関数から所定の
方向成分を抽出することを特徴とする請求項１記載の画
像マッチング方法。
【請求項４】前記周波数分布を求めるとき、元の画像
の周波数成分を差し引いた残りから、周波数成分を算出
することを特徴とする請求項３記載の画像マッチング方
法。
【請求項５】前記周波数分布を求めるとき、元の２つ
の画像のそれぞれの周波数成分、両者の差分成分、およ
び両者の加算成分を差し引いた残りから、周波数成分を
算出することを特徴とする請求項３記載の画像マッチン
グ方法。
【請求項６】前記縮小率または拡大率を１に近い値に
設定するとともに、回転角度を０度に近い値に設定する
ことを特徴とする請求項１記載の画像マッチング方法。