JP6217227B2 - 校正装置、方法及びプログラム - Google Patents

Description

本発明は、三次元計測装置を校正する校正装置、方法及びプログラムに関する。

従来、所定のパターン光を投影した計測対象物を異なる視点から撮影することによって、当該計測対象物の三次元形状を計測する技術が知られている。

このような三次元計測は、投影したパターンと撮影されたパターンの対応付けを行った後、三角測量の原理により計測対象物までの距離を算出するものであるが、その計算の際、投影装置および撮影装置の光学系の内部パラメータ（焦点距離、光軸位置、レンズ歪み等）や外部パラメータ（投影装置および撮影装置の相対位置・姿勢）などが必要となる。

上述したパラメータを事前に取得することを校正（キャリブレーション）と呼ぶ。この校正のやり方としては、既知の形状を持つ専用の校正物体を使用して事前に校正を行う方法と、計測対象物自体を使用して計測時に同時的に校正を行う方法（自己校正法）がある。この自己校正法によれば、専用の校正物体を使用した事前の校正作業が不要となり、システムの経年変化にも柔軟に対応できるという利点があるため、近年、種々の検討がなされている。

このような自己校正法として、特許文献１は、計測対象物から複数の平面領域を抽出した後、各平面上の４点以上の対応点に基づいて投影画像と撮影画像の平面射影行列Ｈを導出し、導出された複数の平面射影行列Ｈ_１、Ｈ_２…Ｈ_ｎに基づいてパラメータを線形計画問題の解として求める方法を開示する。しかしながら、特許文献１の方法は、計測対象物が平面領域を持たない場合（例えば、計測対象物が曲面で構成されている場合）に適用することができない。

この点につき、特許文献２は、計測対象物上の投影点について、カメラとプロジェクターがエピポーラ拘束を満たす場合に、プロジェクターから当該投影点に向かう視線（第１の直線）とカメラから当該投影点に向かう視線（第２の直線）の三次元空間における距離が最小となる誤差関数を定義した上で、当該誤差関数を最小化するパラメータを非線形計画法により推定する方法を開示する。

しかしながら、特許文献２が開示する自己校正法は、エピポーラ拘束を前提とした誤差関数を使用するため、プロジェクターやカメラの光学系がピンホールモデルから乖離している場合（例えば、短焦点プロジェクターや全方位カメラを使用する場合）に適用することができないという問題がある。

本発明は、上記従来技術における課題に鑑みてなされたものであり、本発明は、非ピンホールモデルの光学系を有するシステムに適用することができる新規な自己校正法を提供することを目的とする。

本発明者は、非ピンホールモデルの光学系を有するシステムに適用することができる新規な自己校正法につき鋭意検討した結果、以下の構成に想到し、本発明に至ったのである。

すなわち、本発明によれば、投影装置を制御して計測対象物に校正用パターンを投影させる投影手段と、撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得する撮影手段と、前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出する特徴点抽出手段と、
抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定するパラメータ推定手段とを含み、前記パラメータ推定手段は、撮影画像平面上の前記特徴点と該特徴点に対応する投影光線を前記撮影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定する、校正装置が提供される。

上述したように、本発明によれば、非ピンホールモデルの光学系を有するシステムに適用することができる新規な自己校正法、当該方法を適用した校正装置およびプログラムが提供される。

本実施形態の校正装置の機能ブロック図。 ピンホールカメラ光学系と非ピンホールカメラ光学系を示す図。 本実施形態の校正装置が実行する処理のフローチャート。 校正用パターンとしてのチェッカーパターンを示す図。 本実施形態における評価関数を説明するための概念図。 本実施形態における評価関数を説明するための概念図。 本実施形態における評価関数を説明するための概念図。 本実施形態における評価関数を説明するための概念図。 本実施形態における評価関数を説明するための概念図。 本実施形態における評価関数を説明するための概念図。 本実施形態における評価関数を説明するための概念図。

以下、本発明を、実施形態をもって説明するが、本発明は後述する実施形態に限定されるものではない。なお、以下に参照する各図においては、共通する要素について同じ符号を用い、適宜、その説明を省略するものとする。

図１は、本発明の実施形態である校正装置１００の機能ブロック図を示す。本実施形態の校正装置１００は、投影装置２００および撮影装置３００の双方に通信可能に接続される情報処理装置であり、投影部１０、撮影部１２、記憶部１４、特徴点抽出部１６およびパラメータ推定部１８を含んで構成されている。ここで、撮影装置３００は、CCDやCMOSなどの撮像素子を備えたデジタルカメラとして参照することができ、投影装置２００は、プロジェクターとして参照することができる。

投影部１０は、記憶部１４から所定の校正用パターンを読み出し、投影装置２００を制御して計測対象物に当該校正用パターンを投影させるための機能部である。

撮影部１２は、撮影装置３００を制御して計測対象物を撮影し、取得した撮影画像を記憶部１４に格納するための機能部である。

特徴点抽出部１６は、撮影装置３００が撮影した撮影画像を記憶部１４から読み込み、当該撮影画像から計測対象物に投影されたパターンの特徴点を抽出するための機能部である。

パラメータ推定部１８は、特徴点抽出部１６が抽出した複数の特徴点に基づいて、三次元計測に必要なパラメータを推定するための機能部である。ここで、三次元計測に必要なパラメータとしては、投影装置２００の内部パラメータＬ、撮影装置３００の内部パラメータＫ、投影座標系と撮影座標系の間の回転行列Ｒおよび並進ベクトルｔの他、投影装置２００や撮影装置３００のレンズの歪みなどを挙げることができる。

図２（ａ）は、全ての光線が１つの視点を通過するピンホールカメラモデルの光学系（以下、ピンホール光学系という）を示し、図２（ｂ）は、全ての光線が１点で交差しないピンホールカメラモデルから乖離した光学系（以下、非ピンホール光学系という）を示す。本実施形態の自己校正法は、使用する光学装置（投影装置２００および撮影装置３００）がピンホール光学系および非ピンホール光学系のいずれを採用する場合であっても適用することができる。

以上、本実施形態の校正装置１００を構成する各機能部について概説したが、続いて、上述した各機能部が協働して実行する処理の内容を、図３に示すフローチャートに基づいて説明する。なお、以下の説明においては、適宜、図１を参照するものとする。

校正に先だって、まず、投影装置２００と撮影装置３００の相対位置・姿勢が固定される。その後、ステップ１０１では、投影部１０が記憶部１４から校正用パターンを読み出し、投影装置２００を制御して当該校正用パターンを計測対象物に投影する。校正用パターンは、特徴点を抽出し得るものであればどのようなパターンであってもよいが、以下においては、図４に示すチェッカーパターン５０を校正用パターンとした場合を例にとって説明を行う。

続くステップ１０２では、撮影部１２が撮影装置３００を制御して計測対象物を撮影する。このとき、撮影部１２は、チェッカーパターン５０が投影された計測対象物の撮影画像を取得し、当該撮影画像を記憶部１４に保存する。

続くステップ１０３では、特徴点抽出部１６が記憶部１４に保存された撮影画像を読み込んで画像解析を行い、撮影画像から特徴点を抽出する。具体的には、撮影画像から、計測対象物に投影されたチェッカーパターン５０の正方形部分（白または黒）のコーナー（特徴点）をハリスのコーナー検出法などを用いて検出し、当該特徴点の撮影画像平面上の位置（二次元座標）を取得する。

続くステップ１０４では、ステップ１０３で取得した複数の特徴点に基づいて、パラメータ推定部１８が上述した三次元計測に必要なパラメータ（Ｌ、Ｋ、Ｒ、ｔなど）のうちの少なくとも１つを推定する。以下、ステップ１０４で実行される推定処理の内容を説明する。

ステップ１０４において、パラメータ推定部１８は、図５に示すように、先のステップ１０３で抽出した撮影画像平面上の特徴点に対応する計測対象物上の投影点を通る投影光線Ｍを撮影画像平面に射影した直線Ｎの式を求める。その結果、ステップ１０４では、ステップ１０３で抽出した複数の特徴点１、２、３、…ｎに対して、それぞれ、直線Ｎ_１、Ｎ_２、Ｎ_３…Ｎ_ｎが求められる。

このとき、撮影画像平面では、特徴点の位置抽出誤差に起因して、図６に示すように、各特徴点とこれに対応する各直線Ｎとの間に距離が生じる。この点に関し、パラメータ推定部１８は、直線Ｎと特徴点の間のユークリッド距離（d₁、d₂、d₃、…d_n）を評価するための評価関数を有しており、当該評価関数の出力値を最小化するパラメータを非線形計画法により推定する。

つまり、多くの応用では、撮影画像平面上の特徴点の位置抽出誤差が支配的な誤差となるところ、本実施形態においては、直線Ｎと特徴点の間のユークリッド距離を直接最小化するようにパラメータを推定するので、高精度な校正が可能になる。

なお、全てのパラメータを未知数とすると不定性があり、そのままでは値が一意に決まらないため、推定するパラメータを絞るか、あるいは、評価関数に正則化項をつけるなどして最適化を行う。最適化には最急降下法、ニュートン法、共役勾配法といった非線形計画法のアルゴリズムを用いることができる。

以上、ステップ１０４で実行される推定処理の概要を説明したが、ステップ１０４で実行される処理は、ピンホール光学系の光学装置（投影装置２００および撮影装置３００）を使用する場合と非ピンホール光学系の光学装置（投影装置２００および撮影装置３００）を使用する場合とで、その内容が異なる。この点につき、以下、順を追って説明する。なお、以下の説明においては、図３および図５を適宜参照するものとする。

＜２つの光学装置が両方ともピンホール光学系である場合＞
ここでは、まず、２つの光学装置（投影装置２００および撮影装置３００）の光学系が両方ともピンホール光学系である場合に、パラメータ推定部１８が実行する処理について説明する。

最初に、パラメータ推定部１８が使用する評価関数をその導出過程に沿って説明する。

ステップ１０３で抽出された撮影画像平面上のi番目の特徴点に対応する投影画像平面上のi番目の対応点の位置座標を(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)とすると、i番目の対応点に対応する投影光線Ｍ⁽ⁱ⁾の投影座標系における方向ベクトルe_p ⁽ⁱ⁾は、投影装置２００の焦点距離f_pを用いて下記式（１）で表される。

ここで、投影装置２００の内部パラメータ行列Ｌを下記式（２）で表す場合、方向ベクトルe_p ⁽ⁱ⁾は、下記式（３）で表される。なお、下記式（３）において×は外積を表す。

次に、撮影座標系における投影光線Ｍ⁽ⁱ⁾は、方向ベクトルe_p ⁽ⁱ⁾、回転行列Ｒ、並進ベクトルｔおよび媒介変数sを用いて下記式（４）で表される。

次に、上記式（４）で表される直線を撮影画像平面に投影した直線Ｎ⁽ⁱ⁾は、撮影画像平面の二次元座標を(x_c, y_c)とすると、方向ベクトルe_p ⁽ⁱ⁾、回転行列Ｒ、並進ベクトルｔ、撮影装置３００の内部パラメータ行列Ｋおよび媒介変数sを用いて下記式（５）で表される。

上記式（５）において媒介変数sの値を変動させたときの(x_c, y_c)の軌跡が直線Ｎ⁽ⁱ⁾となる。上記式（５）は、左辺が右辺のスカラー倍であることを表現しており、上記式（５）からスカラー倍の不定性と媒介変数sを消去すると、上記式（５）は下記式（６）に集約される。

ここで、上記式（５）を上記式（６）に集約する考え方について、図７に基づいて説明する。撮影装置３００の撮像中心Ｏを撮影座標系の原点とした場合、z=1の平面が撮影画像平面となる。そして、上記式(５)の右辺はs=0のとき“Kt”となり、s=1のとき“K(Re_p ⁽ⁱ⁾+t）”となり、これらの２点ならびに撮像中心Ｏの三点を通る平面と撮影画像平面(z=1)が交わる線が求める直線Ｎ⁽ⁱ⁾となる。

ここで、一般に、平面の方程式は、当該平面の法線ベクトルを“ｎ”、当該平面上の座標を“ｘ”とした場合、“ｎ^Ｔｘ＝０”で表現することができる（すなわち、平面は、当該平面の法線方向に直交する点の集合として表される）。

ここで、“Kt”、“K(Re_p ⁽ⁱ⁾+t）”および撮像中心Ｏの三点を通る平面の法線ベクトルｎは、撮像中心Ｏから“Kt”に向かうベクトルと、撮像中心Ｏから“K(Re_p ⁽ⁱ⁾+t）”に向かうベクトルの外積で表現することができるので、当該法線ベクトルｎは下記式（７）で表される。

ここで、撮影画像平面を点の集合(x_c, y_c, 1)とすれば、求めるべき直線Ｎ⁽ⁱ⁾の方程式として、集合(x_c, y_c, 1)と上記式（７）で表す法線ベクトルｎとの内積が０となる式として上記式（６）が導出されるのである。

次に、上記式（７）で表す法線ベクトルｎをc⁽ⁱ⁾と置くと、下記式（８）が成立する。

この場合、撮影画像平面上のi番目の特徴点の位置座標(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)と撮影画像平面上の直線Ｎ⁽ⁱ⁾のユークリッド平方距離d_i ²は、下記式（９）で表せる。

そして、本実施形態においては、特徴点抽出部１６が抽出したＮ個の特徴点について求めたユークリッド平方距離d_i ²の和を評価関数Ｊとして定義する。下記式（１０）は評価関数Ｊを表す。

撮影画像平面上の特徴点の位置抽出誤差がガウス分布に従う場合、上記式（１０）に示す評価関数Ｊの最小化によるパラメータ推定は最尤推定となる。また、ガウス分布以外の場合でも、位置抽出誤差の分布が既知であれば、評価関数に分布の形状を反映させることで最尤推定が可能となる。

以上、２つの光学装置が両方ともピンホール光学系である場合に使用する評価関数について説明したが、次に、上記式（１０）に示す評価関数Ｊの出力値を最小化するパラメータを非線形計画法により推定する手順を例示的に説明する。なお、２つの光学装置が両方ともピンホール光学系である場合、スケールが決まらないため、以下のケースにおいては、いずれも、ｔのノルム（スケール）として実測値または設計値を与えるものとする。

（ケース１：Ｒ、ｔを推定する場合）
この場合、Ｌ、Ｋは、従来のピンホールカメラの校正法を利用して予め求めておく。そして、ｔ（三次元ベクトル）のうち、ノルムを除く残りの自由変数をt1、t2とする。具体的には、ｔの第一要素と第二要素を取り出してそれぞれをt1、t2したり、ｔのx軸、y軸に対する偏角をそれぞれt1、t2としたりすることができる。

一方、Ｒはx, y, z軸回りの回転角θx、θy、θzの３変数で表現される。その結果、本ケースでは、t1、t2、θx、θy、θzの５変数が最適化の対象となり、この５変数について、上記式（１０）に示す評価関数Ｊの出力値を最小化する値を非線形計画法により求める。例えば、最急降下法を使用する場合、以下のステップ１〜４を実施してパラメータを推定する。

ステップ１：t1、t2、θx、θy、θzに初期値を設定する。
ステップ２：t1、t2、θx、θy、θzに対する評価関数Ｊの勾配を算出する。
ステップ３：勾配の方向に直線探索を行うことで評価関数Ｊの最小点を求め、t1、t2、θx、θy、θzを更新する。
ステップ４：評価関数Ｊの改善量が閾値以下になるまで、ステップ２→ステップ３を繰り返し、改善量が閾値以下になった時点で解が収束したと判定して、その時点の値を推定値として得る。

（ケース２：Ｒ、ｔに加えて、Ｌ、Ｋの焦点距離を推定する場合）
Ｌ、Ｋのうち焦点距離以外の値は事前に与えておく。その上で、撮影装置３００の焦点距離をf_c、投影装置２００の焦点距離をf_pとし、t1、t2、θx、θy、θz、f_c、f_pの７変数について、ケース１と同様の手順で推定値を得る。なお、f_cおよびf_pのうちの少なくとも一方を他の任意のパラメータ（例えば、光軸位置）に置き換えて推定することもできる。

（ケース３：Ｒ、ｔ、Ｌ、Ｋのすべてを推定する場合）
すべてのパラメータを未知とする場合には、評価関数Ｊに正則化項を付与して最適化を行う。ここで、正則化項としては、事前値からのずれ量を設定することができ、例えば、パラメータ（Ｒ、ｔ、Ｌ、Ｋ）の事前値が（R0、t0、L0、K0）であった場合、下記式（１１）に示す正則化項を評価関数Ｊに与えることができる。なお、下記式（１１）において、事前値（L0、K0）は、従来のピンホールカメラの校正法を利用して予め求めておき、事前値（R0、t0）は、撮影装置３００と投影装置２００の相対位置・姿勢の実測値または設計値を与える。

上記式（１１）に示す正則化項を評価関数Ｊに加算したものを新しい評価関数Ｊ’とし、評価関数Ｊ’を使用して上述したのと同様の手順でＲ、ｔ、Ｌ、Ｋの推定値を得る。なお、回転行列Ｒの正則化項は“||R-R0||^2”の代わりに、θx、θy、θzの事前値θ0x、θ0y、θ0zを使って、“(θx-θ0x)^2 + (θy-θ0y)^2 + (θz-θ0z)^2”とすることもできる。

以上、２つの光学装置（投影装置２００および撮影装置３００）の光学系が両方ともピンホール光学系である場合について説明してきたが、次に、２つの光学装置（撮影装置３００および投影装置２００）の一方の光学系がピンホール光学系であり、他方の光学系が非ピンホール光学系である場合にパラメータ推定部１８が実行する処理について説明する。

＜２つの光学装置の一方の光学系がピンホール光学系であり、他方の光学系が非ピンホール光学系である場合＞
ここでは、撮影装置３００の光学系がピンホール光学系であり、投影装置２００の光学系が非ピンホール光学系である場合について説明する。投影装置２００が非ピンホール光学系である場合、投影装置２００は、投影画像平面上の各画素に投影光線の方程式を紐付けてなる光線モデルを具備している。

ここで、投影光線は、投影座標系における投影光線の方向ベクトルe_ｐと当該投影光線の通過位置p_pで表される直線の方程式群としてモデル化される。例えば、本実施形態における光線モデルは、投影画像平面上の画素位置（x,y）と方向ベクトルe_ｐおよび通過位置p_pを対応付けたルックアップテーブルとして構築することができ、投影画像平面上の画素位置（x,y）を変数とする二変数多項式として表現することもできる。

上述した光線モデルの存在を前提として、投影装置２００が非ピンホール光学系である場合にパラメータ推定部１８が使用する評価関数をその導出過程に沿って説明する。

ステップ１０３で抽出された撮影画像平面上のi番目の特徴点に対応する投影画像平面上のi番目の対応点の位置座標を(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)とした場合、図８に示すように、上述した光線モデルに(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)を入力することによって、i番目の対応点に対応する投影光線Ｍ⁽ⁱ⁾の投影座標系における方向ベクトルe_p ⁽ⁱ⁾および投影光線Ｍ⁽ⁱ⁾の通過位置p_p ⁽ⁱ⁾が得られる。

そして、撮影座標系における投影光線Ｍ⁽ⁱ⁾は、回転行列Ｒ、並進ベクトルｔおよび媒介変数sを用いて下記式（１２）で表される。

次に、上記式（１２）で表される直線を撮影画像平面に投影した直線Ｎ⁽ⁱ⁾は、撮影画像平面の二次元座標を(x_c, y_c)とすると方向ベクトルe_p ⁽ⁱ⁾、通過位置p_p ⁽ⁱ⁾、回転行列Ｒ、並進ベクトルｔ、撮影装置３００の内部パラメータ行列Ｋおよび媒介変数sを用いて下記式（１３）で表される。

次に、先に説明したのと同様に直線Ｎ⁽ⁱ⁾を含む平面の法線ベクトルをc⁽ⁱ⁾と置くと、下記式（１４）が成立する。

この場合、撮影画像平面上のi番目の特徴点の位置座標(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)と撮影画像平面上の直線Ｎ⁽ⁱ⁾とのユークリッド平方距離d_i ²は、下記式（１５）で表せる。

本実施形態においては、特徴点抽出部１６が抽出したＮ個の特徴点について求めたユークリッド平方距離d_i ²の和を評価関数Ｊとして定義する。下記式（１６）は評価関数Ｊを表す。

続いて、上記式（１６）に示す評価関数Ｊの出力値を最小化するパラメータを非線形計画法により推定する手順を例示的に説明する。

（ケース１：Ｒ、ｔを推定する場合）
この場合、Ｋが既知であり、且つ、Ｌ（e_pおよびp_p）が予め用意された光線モデルにより与えられることを前提として、Ｒ、ｔを推定する。本ケースでは、投影装置２００の光学系の非ピンホール性を前提としているため、ｔのノルム（スケール）を推定対象とすることができる（以下のケースにおいて同様）。したがって、本ケースでは、ｔの自由変数をt1、t2、t3とし、Ｒの自由変数をx, y, z軸回りの回転角θx、θy、θzとした６変数について、上記式（１６）に示す評価関数Ｊの出力値を最小化する値を非線形計画法により求める。

（ケース２：Ｒ、ｔ、Ｌ、Ｋのすべてを推定する場合）
Ｋのうち撮影装置３００の焦点距離f_c以外の値は事前に与えておき、f_c t1、t2、t3、θx、θy、 θzの７変数を推定対象とする。

加えて、Ｌについては、方向ベクトルe_pおよび通過位置p_pを多項式で近似した上で、その係数を推定対象としたり、方向ベクトルe_pおよび通過位置p_pを少数の基底の線形和で表現した上で、その結合係数を推定対象としたりすることができる。

ここでは、まず、方向ベクトルe_pおよび通過位置p_pを投影画像平面上の対応点の位置（x,y）を変数とする二変数多項式で近似する場合について説明する。下記式（１７）は、方向ベクトルe_pに対応する多項式を、下記式（１８）は、通過位置p_pに対応する多項式をそれぞれ例示的に示す。

この場合、上述した７変数（f_c t1、t2、t3、θx、θy、 θz）、ならびに、上記式（１７）の係数（行列a⁽⁰⁾,a⁽¹⁾,a⁽²⁾）、ならびに、上記式（１８）の係数（行列b⁽⁰⁾,b⁽¹⁾,b⁽²⁾）を推定対象として、ケース１と同様の手順で非線形計画法により推定値を求める。

次に、方向ベクトルe_pおよび通過位置p_pを、投影画像平面上の対応点の位置（x,y）を変数とする基底の線形和で表現する場合について説明する。例えば、投影装置２００がズーム機構を備えている場合、最小ズーム時および最大ズーム時のそれぞれにおけるＬ（e_pおよびp_p）を事前に求めておく。その上で、方向ベクトルe_pの値がズームに応じて線形に変化するという前提のもと、方向ベクトルe_pを下記式（１９）で表現する。

上記式（１９）において、e₀(x,y)は最小ズーム時における投影光線の方向ベクトルを示し、e₁(x,y)は最大ズーム時における投影光線の方向ベクトルを示し、wは線形和の重み係数を示す。

同様に、通過位置p_pの値がズームに応じて線形に変化するという前提のもと、通過位置p_pを下記式（２０）で表現する。

上記式（２０）において、p₀(x,y)は最小ズーム時における投影光線の通過位置を示し、p₁(x,y)は最大ズーム時における投影光線の通過位置を示し、wは線形和の重み係数を示す。

この場合、上述した７変数（f_c t1、t2、t3、θx、θy、 θz）、ならびに、上記式（１９）および上記式（２０）において共有される重みwを推定対象として、ケース１と同様の手順で非線形計画法により推定値を求める。

次に、投影装置の個体変動を求める場合について説明する。この場合、複数の個体について、Ｌ（e_pおよびp_p）の平均値、ならびに、各個体の変動ベクトルを事前に求めておく。その上で、方向ベクトルe_pを下記式（２１）で表現する。

上記式（２１）において、e₀(x,y)は複数の個体の投影光線の方向ベクトルの平均値を示し、e₁(x,y)は投影光線の方向ベクトル変動を示し、wは線形和の重み係数を示す。

同様に、通過位置p_pを下記式（２２）で表現することができる。

上記式（２２）において、p₀(x,y)は複数の個体の投影光線の平均通過位置を示し、p₁(x,y)は通過位置変動を示し、wは線形和の重み係数を示す。

この場合、上述した７変数（f_c t1、t2、t3、θx、θy、 θz）、ならびに、上記式（２１）および上記式（２２）において共有される重みwを推定対象として、ケース１と同様の手順で非線形計画法により推定値を求める。

以上、Ｌについて、e_pおよびp_pを多項式で近似してその係数を推定対象とする場合、ならびに、e_pおよびp_pをより少数の基底の線形和で表現しその重み係数を推定対象とする場合について説明してきたが、より高次の多項式近似を行ったり、線形和の基底の数を増やしたりする場合は、変数が増加するので、上述したのと同様に、評価関数Ｊに対して正則化項を付与し、各変数に対して用意した事前値とのずれ量を評価関数に加算して最適化を行う。

以上、撮影装置３００がピンホール光学系であり、投影装置２００が非ピンホール光学系である場合について説明してきたが、次に、２つの光学装置（投影装置２００および撮影装置３００）の光学系が両方とも非ピンホール光学系である場合にパラメータ推定部１８が実行する処理について説明する。

＜２つの光学装置の光学系が両方とも非ピンホール光学系である場合＞
既に、上述した通り、非ピンホール光学系の投影装置２００から出射する投影光線の方程式群は、投影光線の投影座標系における方向ベクトルe_ｐと当該投影光線の通過位置p_pによってモデル化される。

同様に、非ピンホール光学系の撮影装置３００に入射する撮影光線の方程式群は、撮影座標系における撮影光線の方向ベクトルe_ｃと当該撮影光線の通過位置p_cによってモデル化される。

上述した光線モデルの存在を前提として、２つの光学装置の光学系が両方とも非ピンホール光学系である場合にパラメータ推定部１８が使用する評価関数をその導出過程に沿って説明する。

ステップ１０３で抽出された撮影画像平面上のi番目の特徴点に対応する投影画像平面上のi番目の対応点の位置座標を(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)とした場合、上述した光線モデルに(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)を入力することによって、図８に示すように、i番目の対応点に対応する投影光線Ｍ⁽ⁱ⁾の投影座標系における方向ベクトルe_p ⁽ⁱ⁾および投影光線Ｍ⁽ⁱ⁾の通過位置p_p ⁽ⁱ⁾が得られる。

そして、撮影座標系における投影光線Ｍ⁽ⁱ⁾は、方向ベクトルe_p ⁽ⁱ⁾、通過位置p_p ⁽ⁱ⁾回転行列Ｒ、並進ベクトルｔおよび媒介変数sを用いて下記式（２３）で表される。

ここで、上記式（２３）で表される投影光線Ｍ⁽ⁱ⁾に対応する撮影画像平面上の像を考える。撮影装置３００がピンホール系の場合は、先に説明したように、上記式（２３）を射影変換するだけで、投影光線Ｍ⁽ⁱ⁾に対応する撮影画像平面上の像（直線）を得ることができた。しかしながら、撮影装置３００が非ピンホール系である場合、投影光線Ｍ⁽ⁱ⁾は撮影画像平面上に曲線として投影されるので、投影光線Ｍ⁽ⁱ⁾に対応する撮影画像平面上の射影像を上記式（２３）の射影変換によっては求めることができない。

この点につき、本実施形態においては、図９に示すように、撮影座標系における任意の点Ｘ（三次元ベクトル）を通る撮影光線が入射する撮影画像平面上の撮影点（二次元座標）を導出する関数をｍ（Ｘ）と置く（この関数については後述する）。この場合、投影光線Ｍ⁽ⁱ⁾に対応する撮影画像平面上の曲線Ｎ⁽ⁱ⁾上の点の二次元座標は、下記式（２４）で表される。

ここで、ステップ１０３で抽出された撮影画像平面上のi番目の特徴点の二次元座標を(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)とすると、特徴点(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)と曲線Ｎ⁽ⁱ⁾上の任意の点のユークリッド平方距離は、媒介変数sを用いて下記式（２５）で表される。

そして、本実施形態においては、図１０に示すように、撮影画像平面上の曲線Ｎ⁽ⁱ⁾と特徴点(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)のユークリッド平方距離を曲線Ｎ⁽ⁱ⁾上の点と特徴点(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)のユークリッド平方距離の最小値と定義すると、撮影画像平面上の曲線Ｎ⁽ⁱ⁾と特徴点(x_c ⁽ⁱ⁾, y_c ⁽ⁱ⁾)のユークリッド平方距離d_i ²は、下記式（２６）で表せる。

本実施形態においては、特徴点抽出部１６が抽出したＮ個の特徴点について求めたユークリッド平方距離d_i ²の和を評価関数Ｊとして定義する。下記式（２７）は評価関数Ｊを表す。

そして、本実施形態においても、上記式（２７）に示す評価関数Ｊの出力値を最小化するパラメータを非線形計画法により推定する。なお、パラメータの推定手順の詳細は、投影装置２００が非ピンホール光学系である場合について上述したのと同様であるので、これ以上の説明を省略する。

ここで、本実施形態においては、撮影座標系における任意の点Ｘ（三次元ベクトル）を通る撮影光線が撮像される撮影画像平面上の撮影点（二次元座標）を導出する関数ｍ（Ｘ）を以下の方法で実現することができる。

（最近傍法）
撮影装置３００に係る撮影光線の方程式群の中から、点Ｘに最も近い撮影光線を抽出し、当該撮影光線に対応付けられた撮影画像平面上の二次元座標(x, y)を出力するｍ（Ｘ）を定義する。この点につき、図１１に基づいて説明する。

撮影画像平面上の画素位置(x, y)に撮像される撮影光線の撮影座標系における方向ベクトルをe_c(x, y)、その通過位置をp_c(x, y)とすると、点Ｘと当該撮影光線のユークリッド平方距離d_i ²（Ｘ，ｘ、ｙ）は下記式（２８）で表される。

この場合、ｍ（Ｘ）は、下記式（２９）に示すように、このユークリッド平方距離d_i ²（Ｘ，ｘ、ｙ）を最小化する(x, y)を出力する関数として定式化される。

ここで、d_i ²（Ｘ，ｘ、ｙ）を最小化する(x, y)の導出は、(x, y)が離散的である場合（x, yが格子点となっているなど）は、離散的な(x, y)の組すべてについて距離を算出して最小値を見つける方法（全探索）の他、山登り法やシミュレーティッドアニーリングといったより計算コストの小さい離散最適化手法を適用する方法がある。一方、(x, y)が連続値を取る場合は、(x, y)に対するd_i ²（Ｘ，ｘ、ｙ）の勾配を計算して、最急降下法、ニュートン法といった方法でd_i ²（Ｘ，ｘ、ｙ）の最小化を行う。

（k-最近傍法）
(x, y)が離散的である場合、上述した最近傍法では撮影光線が空間中に密にないケースで精度が低くなる。そこで、点Ｘに近いものから順番にk個の撮影光線を抽出し、当該k個の撮影光線に対応する撮影画像平面上の画素位置の平均座標を出力するｍ（Ｘ）を定義する。この場合、平均座標を求めるにあたって、点Ｘに近い撮影光線ほど大きい重みを付与するようにしてもよい。

以上、説明したように、本実施形態の自己校正法は、非ピンホールモデルの光学系を有するシステムに適用することができる。また、本実施形態の自己校正法は、撮影画像平面上の特徴点の位置抽出誤差が支配的であることに着目して、直線Ｎと特徴点の間のユークリッド距離を直接最小化するようにパラメータを推定するので、高精度な校正が可能になる。

さらに、非ピンホール光学系のために用意された上記アルゴリズムによれば、ｔのノルム（スケール）を推定対象とすることができる。また、現実には、ピンホールモデルに完全に合致する光学装置は存在しないので、２つの光学装置の両方がピンホール光学系を採用するケースにおいても、非ピンホール光学系のために用意された上記アルゴリズムを適用することは排除されず、その場合、ｔのノルム（スケール）を推定対象とすることができる。

以上、本発明について実施形態をもって説明してきたが、本発明は上述した実施形態に限定されるものではない。

まず、第１に、２つの光学装置の光学系が両方ともピンホール光学系である場合、上述した手順に代えて、撮影画像平面上のi番目の特徴点に対応する投影画像平面上のi番目の対応点の位置座標(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)と、当該i番目の対応点に対応する撮影光線を投影画像平面に射影してなる直線とのユークリッド平方距離d_i ²の和を評価関数とし、当該評価関数の出力値を最小化するパラメータを非線形計画法により推定することもできる。

第２に、２つの光学装置の光学系が両方とも非ピンホール光学系である場合、上述した手順に代えて、投影座標系における任意の点Ｘ（三次元ベクトル）を通る投影光線が出射する投影画像平面上の対応点（二次元座標）を導出する関数をｍ（Ｘ）と置き、当該ｍ（Ｘ）を使用して、撮影画像平面上のi番目の特徴点に対応するｉ番目の撮影光線を投影画像平面に射影してなる直線を求め、投影画像平面上のi番目の対応点の位置座標(x_p ⁽ⁱ⁾, y_p ⁽ⁱ⁾)と当該直線とのユークリッド平方距離d_i ²の和を評価関数とし、当該評価関数の出力値を最小化するパラメータを非線形計画法により推定することもできる。

つまり、本発明によれば、
（１）撮影装置および投影装置が両方ともピンホール光学系
（２）撮影装置がピンホール光学系で投影装置が非ピンホール光学系
（３）撮影装置が非ピンホール光学系で投影装置がピンホール光学系
（４）撮影装置および投影装置が両方とも非ピンホール光学系
の４通りの組み合わせすべてについて自己校正を行うことができる。

さらに、上述した実施形態においては、撮影装置と投影装置を1つずつ使用する態様について説明したが、２以上の撮影装置ならびに２以上の投影装置を使用する場合にも本発明を適用することができる。

ここでは、N_p個の投影装置とN_c個の撮影装置を使用する場合を考える。ここで、i番目の投影装置とj番目の撮影装置の組み合わせに係る評価関数J_ijを上述した手順で定義した上で、さらに評価関数J_ijを用いて、下記式（３０）に示す評価関数Ｊを新たに定義する。

そして、上記式（３０）が示す評価関数Ｊの出力値を最小化するパラメータを非線形計画法により推定する。

さらに、上述した実施形態においては、投影装置および撮影装置のレンズの歪みについて言及しなかったが、放射方向歪みモデル、接線方向歪みモデルなどを用いてレンズの歪みをパラメタライズすることで、上述したのと同様の方法でレンズ歪みを含めたパラメータを推定することもできる。

さらに、上述した実施形態においては、校正用パターンとしてチェッカーパターンを例示したが、この他にも、円パターンを用いて円の中心点を特徴点として抽出したり、正弦波状のパターンを複数個組み合わせて用い、正弦波の位相によって特徴点の位置を識別するようにしたりすることもできる（位相シフト法）。

その他、当業者が推考しうる実施態様の範囲内において、本発明の作用・効果を奏する限り、本発明の範囲に含まれるものである。

また、上述した実施形態の各機能は、Ｃ、Ｃ＋＋、Ｃ＃、Ｊａｖａ（登録商標）などで記述された装置実行可能なプログラムにより実現でき、本実施形態のプログラムは、ハードディスク装置、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ、フレキシブルディスク、ＥＥＰＲＯＭ、ＥＰＲＯＭなどの装置可読な記録媒体に格納して頒布することができ、また他装置が可能な形式でネットワークを介して伝送することができる。

１０…投影部
１２…撮影部
１４…記憶部
１６…特徴点抽出部
１８…パラメータ推定部
５０…チェッカーパターン
１００…校正装置
２００…投影装置
３００…撮影装置

特開２００５−３２６２４７号公報 特許第４２３０５２５号公報

Claims (13)

1. 投影装置を制御して計測対象物に校正用パターンを投影させる投影手段と、
   撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得する撮影手段と、
   前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出する特徴点抽出手段と、
   抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定する手段であって、撮影画像平面上の前記特徴点と該特徴点に対応する投影光線を撮影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するパラメータ推定手段と
   を含み、
   前記投影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、投影画像平面上の二次元座標を変数とする多項式として表現されており、
   前記パラメータ推定手段は、前記多項式の係数を前記パラメータとして推定する、
   校正装置。
2. 投影装置を制御して計測対象物に校正用パターンを投影させる投影手段と、
   撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得する撮影手段と、
   前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出する特徴点抽出手段と、
   抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定する手段であって、撮影画像平面上の前記特徴点と該特徴点に対応する投影光線を撮影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するパラメータ推定手段と
   を含み、
   前記投影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、投影画像平面上の二次元座標を変数とする基底の線形和として表現されており、
   前記パラメータ推定手段は、前記線形和の重み係数を前記パラメータとして推定する、
   校正装置。
3. 前記撮影手段の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする多項式として表現されており、
   前記パラメータ推定手段は、前記多項式の係数を前記パラメータとして推定する、
   請求項１または２に記載の校正装置。
4. 前記撮影手段の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする基底の線形和として表現されており、
   前記パラメータ推定手段は、前記線形和の重み係数を前記パラメータとして推定する、
   請求項１または２に記載の校正装置。
5. 投影装置と撮影装置を含む三次元計測装置を校正する方法であって、
   投影装置を制御して計測対象物に校正用パターンを投影するステップと、
   撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得するステップと、
   前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出するステップと、
   抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定するステップであって、撮影画像平面上の前記特徴点と該特徴点に対応する投影光線を撮影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するステップと
   を含み、
   前記投影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、投影画像平面上の二次元座標を変数とする多項式として表現されており、
   前記パラメータを推定するステップは、前記多項式の係数を前記パラメータとして推定するステップを含む、
   方法。
6. 投影装置と撮影装置を含む三次元計測装置を校正する方法であって、
   投影装置を制御して計測対象物に校正用パターンを投影するステップと、
   撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得するステップと、
   前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出するステップと、
   抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定するステップであって、撮影画像平面上の前記特徴点と該特徴点に対応する投影光線を撮影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するステップと
   を含み、
   前記投影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、投影画像平面上の二次元座標を変数とする基底の線形和として表現されており、
   前記パラメータを推定するステップは、前記線形和の重み係数を前記パラメータとして推定するステップを含む、
   方法。
7. 前記撮影装置の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする多項式として表現されており、
   前記パラメータを推定するステップは、前記多項式の係数を前記パラメータとして推定するステップを含む、
   請求項５または６に記載の方法。
8. 前記撮影装置の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする基底の線形和として表現されており、
   前記パラメータを推定するステップは、前記線形和の重み係数を前記パラメータとして推定するステップを含む、
   請求項５または６に記載の方法。
9. 請求項５〜８のいずれか一項に記載の方法の各ステップをコンピュータに実行させるためのコンピュータ実行可能なプログラム。
10. 投影装置を制御して計測対象物に校正用パターンを投影させる投影手段と、
    撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得する撮影手段と、
    前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出する特徴点抽出手段と、
    抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定する手段であって、撮影画像平面上の前記特徴点に対応する投影画像平面上の対応点と、該対応点に対応する撮影光線を投影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するパラメータ推定手段と
    を含み、
    前記撮影手段の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする多項式として表現されており、
    前記パラメータ推定手段は、前記多項式の係数を前記パラメータとして推定する、校正装置。
11. 投影装置を制御して計測対象物に校正用パターンを投影させる投影手段と、
    撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得する撮影手段と、
    前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出する特徴点抽出手段と、
    抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定する手段であって、撮影画像平面上の前記特徴点に対応する投影画像平面上の対応点と、該対応点に対応する撮影光線を投影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するパラメータ推定手段と
    を含み、
    前記撮影手段の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする基底の線形和として表現されており、
    前記パラメータ推定手段は、前記線形和の重み係数を前記パラメータとして推定する、
    校正装置。
12. 投影装置と撮影装置を含む三次元計測装置を校正する方法であって、
    投影装置を制御して計測対象物に校正用パターンを投影するステップと、
    撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得するステップと、
    前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出するステップと、
    抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定するステップであって、撮影画像平面上の前記特徴点に対応する投影画像平面上の対応点と、該対応点に対応する撮影光線を投影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するステップと
    を含み、
    前記撮影装置の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする多項式として表現されており、
    前記パラメータを推定するステップは、前記多項式の係数を前記パラメータとして推定するステップを含む、
    方法。
13. 投影装置と撮影装置を含む三次元計測装置を校正する方法であって、
    投影装置を制御して計測対象物に校正用パターンを投影するステップと、
    撮影装置を制御して前記校正用パターンが投影された前記計測対象物の撮影画像を取得するステップと、
    前記撮影画像から前記計測対象物に投影された前記校正用パターンの特徴点を抽出するステップと、
    抽出された複数の前記特徴点に基づいて三次元計測に必要なパラメータを推定するステップであって、撮影画像平面上の前記特徴点に対応する投影画像平面上の対応点と、該対応点に対応する撮影光線を投影画像平面に射影した線との距離を最小化する前記パラメータを非線形計画法により推定するステップと
    を含み、
    前記撮影装置の撮影光線は、方向ベクトルと通過位置で表される直線の方程式群としてモデル化され、該方向ベクトルおよび該通過位置のそれぞれは、撮影画像平面上の二次元座標を変数とする基底の線形和として表現されており、
    前記パラメータを推定するステップは、前記線形和の重み係数を前記パラメータとして推定するステップを含む、
    方法。