(57)【要約】【課題】  ビデオカメラ等を使って手動操作により撮影した画像から外界の物体の空間情報を獲得、復元するときに、それぞれ焦点距離が変化するため、従来の因子分解法を利用することができず、ピント調整のために変化する焦点距離を同時に復元しなければならないという問題を解決する。【解決手段】  カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測部４２、画像座標値から行列データを作成し、カメラ手段のカメラ視点運動、及び前対象物の３次元形状を復元する因子分解処理部４３、復元したカメラ視点運動及び３次元形状から再投影座標値を求め、この再投影座標値、及び観測した画像座標値からカメラ手段のレンズ歪パラメータを推定するパラメータ推定部４４を備えるレンズ歪推定装置。【選択図】  図４
2006195790.tif 000002
【特許請求の範囲】
【請求項１】
カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測手段と、
前記画像座標値から行列データを作成し、当該行列データから前記カメラ手段のカメラ視点運動、及び前記対象物の３次元形状を復元する因子分解処理手段と、
前記復元したカメラ視点運動及び３次元形状から再投影座標値を求め、当該再投影座標値、及び前記観測した画像座標値から前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定手段と、
前記行列データから評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、収束していない場合には現時点で有しているレンズ歪パラメータを用いて前記行列データを更新して前記因子分解処理手段からの処理を反復させる反復処理手段と、を備えることを特徴とするレンズ歪推定装置。
【請求項２】
カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測手段と、
前記画像座標値にレンズ歪係数及び主点を用いて座標変換した座標値を要素とする漸近行列を生成し、当該漸近行列を特異値分解し、雑音除去を行って運動情報を表す行列データ、及び３次元情報を表す行列データを得て、前記運動情報の行列データの行列要素において、運動を規定するために設定した条件を満足する変換行列を求め、運動情報を表す行列データにこの変換行列を作用させてカメラ視点運動を復元し、前記３次元情報を表す行列データにこの変換行列の逆行列を作用させて物体形状を構成する３次元情報を復元する因子分解処理手段と、
前記復元したカメラ視点運動及び３次元形状から再投影座標値を求め、当該再投影座標値を関係付ける未校正投影モデルにおいてレンズ歪係数及び主点を未知数とした評価関数を設定し、当該評価関数を最適にするレンズ歪係数及び主点を求めることにより前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定手段と、
前記漸近行列を行列分解したときの主成分からカメラ視点運動への漸近値を表す評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、収束していない場合には現時点で有しているレンズ歪パラメータを用いて前記行列データを更新して前記因子分解処理手段からの処理を反復させる反復処理手段と、を備えることを特徴とするレンズ歪推定装置。
【請求項３】
カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測手段と、
前記画像座標値から行列データを作成し、当該行列データから前記カメラ手段のカメラ視点運動、及び前記対象物の３次元形状を復元する因子分解処理手段と、
前記復元したカメラ視点運動及び３次元形状から再投影座標値を求め、当該再投影座標値、及び前記観測した画像座標値から前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定手段と、
前記再投影座標値、及び前記画像座標値から焦点距離を推定する焦点距離推定手段と、
前記行列データから評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、当該判断に基づいて前記パラメータ推定手段、又は前記焦点距離推定手段のいずいれに処理を行わせるかを決定する反復処理手段と、を備えることを特徴とするレンズ歪推定装置。
【請求項４】
カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測手段と、
前記画像座標値にレンズ歪係数及び主点を用いて座標変換した座標値を要素とする漸近行列を生成し、当該漸近行列を特異値分解し、雑音除去を行って運動情報を表す行列データ、及び３次元情報を表す行列データを得て、前記運動情報の行列データの行列要素において、運動を規定するために設定した条件を満足する変換行列を求め、運動情報を表す行列データにこの変換行列を作用させてカメラ視点運動を復元し、前記３次元情報を表す行列データにこの変換行列の逆行列を作用させて物体形状を構成する３次元情報を復元する因子分解処理手段と、
前記復元したカメラ視点運動及び３次元形状から再投影座標値を求め、当該再投影座標値を関係付ける未校正投影モデルにおいてレンズ歪係数及び主点を未知数とした評価関数を設定し、当該評価関数を最適にするレンズ歪係数及び主点を求めることにより前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定手段と、
前記再投影座標値、並びに前記画像座標値に前記レンズ歪係数及び主点を用いて変換した座標値から構成される値を要素とする行列データから焦点距離を推定する焦点距離推定手段と、
前記漸近行列を行列分解したときの主成分からカメラ視点運動への漸近値を表す評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、当該判断に基づいて前記パラメータ推定手段、又は前記焦点距離推定手段のいずいれに処理を行わせるかを決定する反復処理手段と、を備えることを特徴とするレンズ歪推定装置。
【請求項５】
特徴点観測手段が、カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測ステップと、
因子分解処理手段が、前記画像座標値から行列データを作成し、当該行列データから前記カメラ手段のカメラ視点運動、及び前記対象物の３次元形状を復元する因子分解処理ステップと、
パラメータ推定手段が、前記復元したカメラ視点運動及び３次元形状から再投影座標値を求め、当該再投影座標値、及び前記観測した画像座標値から前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定ステップと、
反復処理手段が、前記行列データから評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、収束していない場合には現時点で有しているレンズ歪パラメータを用いて前記行列データを更新して前記因子分解処理ステップからの処理を反復させる反復処理ステップと、を備えることを特徴とするレンズ歪推定方法。
【請求項６】
特徴点観測手段が、カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測ステップと、
因子分解処理手段が、前記画像座標値にレンズ歪係数及び主点を用いて座標変換した座標値を要素とする漸近行列を生成し、当該漸近行列を特異値分解し、雑音除去を行って運動情報を表す行列データ、及び３次元情報を表す行列データを得て、前記運動情報の行列データの行列要素において、運動を規定するために設定した条件を満足する変換行列を求め、運動情報を表す行列データにこの変換行列を作用させてカメラ視点運動を復元し、前記３次元情報を表す行列データにこの変換行列の逆行列を作用させて物体形状を構成する３次元情報を復元する因子分解処理ステップと、
パラメータ推定手段が、前記復元したカメラ視点運動及び３次元形状から再投影座標値を求め、当該再投影座標値を関係付ける未校正投影モデルにおいてレンズ歪係数及び主点を未知数とした評価関数を設定し、当該評価関数を最適にするレンズ歪係数及び主点を求めることにより前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定ステップと、
反復処理手段が、前記漸近行列を行列分解したときの主成分からカメラ視点運動への漸近値を表す評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、収束していない場合には現時点で有しているレンズ歪パラメータを用いて前記行列データを更新して前記因子分解処理ステップからの処理を反復させる反復処理ステップと、を備えることを特徴とするレンズ歪推定方法。
【請求項７】
特徴点観測手段が、カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測ステップと、
因子分解処理手段が、前記画像座標値から行列データを作成し、当該行列データから前記カメラ手段のカメラ視点運動、及び前記対象物の３次元形状を復元する因子分解処理ステップと、
パラメータ推定手段が、前記復元したカメラ視点運動、及び３次元形状から再投影座標値を求め、当該再投影座標値、及び前記観測した画像座標値から前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定ステップと、
焦点距離推定手段が、前記再投影座標値、及び前記画像座標値から焦点距離を推定する焦点距離推定ステップと、
反復処理手段が、前記行列データから評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、当該判断に基づいて前記パラメータ推定ステップ、又は前記焦点距離推定ステップのいずいれの処理を行うかを決定する反復処理ステップと、を備えることを特徴とするレンズ歪推定方法。
【請求項８】
特徴点観測手段が、カメラ手段によって取得した対象物の時系列画像に設定した画像座標系において、各画像の特徴点の時間的変化量を示す画像座標値を観測する特徴点観測ステップと、
因子分解処理手段が、前記画像座標値にレンズ歪係数及び主点を用いて座標変換した座標値を要素とする漸近行列を生成し、当該漸近行列を特異値分解し、雑音除去を行って運動情報を表す行列データ、及び３次元情報を表す行列データを得て、前記運動情報の行列データの行列要素において、運動を規定するために設定した条件を満足する変換行列を求め、運動情報を表す行列データにこの変換行列を作用させてカメラ視点運動を復元し、前記３次元情報を表す行列データにこの変換行列の逆行列を作用させて物体形状を構成する３次元情報を復元する因子分解処理ステップと、
パラメータ推定手段が、前記復元したカメラ視点運動、及び３次元形状から再投影座標値を求め、当該再投影座標値を関係付ける未校正投影モデルにおいてレンズ歪係数及び主点を未知数とした評価関数を設定し、当該評価関数を最適にするレンズ歪係数及び主点を求めることにより前記カメラ手段のレンズ歪パラメータを推定するパラメータ推定ステップと、
焦点距離推定手段が、前記再投影座標値、並びに前記画像座標値に前記レンズ歪係数及び主点を用いて変換した座標値から構成される値を要素とする行列データから焦点距離を推定する焦点距離推定ステップと、
反復処理手段が、前記漸近行列を行列分解したときの主成分からカメラ視点運動への漸近値を表す評価値を算出して、当該評価値が規定した許容値以下に収束しているかを判断し、当該判断に基づいて前記パラメータ推定ステップ、又は前記焦点距離推定ステップのいずいれの処理を行うかを決定する反復処理ステップと、を備えることを特徴とするレンズ歪推定方法。
【請求項９】
上記請求項１〜８のいずれかに記載のレンズ歪推定装置、又はレンズ歪推定方法を、コンピュータプログラムで記載してそれを実行可能にしたことを特徴とするレンズ歪推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、カメラ等の画像入力装置を使って取得した車載画像、室内画像、歩行撮影画像、船上からの海上画像、空撮画像などの時系列画像全般に利用可能で、取得した時系列画像から、カメラ視点の並進運動と光軸周りの回転（ヨー回転）、及び、時系列画像に映っている外界の３次元形状、すなわち被写体（物体）の外観形状を構成する３次元情報を復元し、同時にカメラ内部のレンズ歪パラメータ、及び焦点距離を推定する技術に関するものである。
【背景技術】
【０００２】
コンピュータビジョン分野では、時系列画像から対象物の形状を計測又は獲得する手法として、ステレオ計測やエピポーラ面解析を用いた３次元解析手法がある。この手法によれば、物体が撮影されている複数の時系列画像から、空間形状又は空間構造に関する３次元位置情報、及びカメラ視点に関する運動を復元することができる。しかし、カメラを動かしながら撮影した時系列画像にはランダム性の雑音が混入し、カメラ運動や物体形状を高精度に復元することが困難である。これに対して、因子分解法はランダム性の雑音が混入するようなビデオ画像などの時系列画像において、ユークリッド空間でのカメラ運動と物体形状をロバスト、かつ、高精度に復元することができる。しかし、この手法では、カメラの内部パラメータ（レンズ歪、焦点距離など）を既知としており、事前にカメラキャリブレーションによりカメラ内部パラメータを校正しておく必要がある。
【０００３】
ディジタル・スチルカメラやハンディカメラなどの市販カメラの光学系は複数のレンズで構成されており、理想的なピンホールカメラモデルから離れている。例えば、図１２に示すように、矩形形状１２２をカメラ１２１で撮影すると得られる画像１２３においてその矩形形状は樽型に歪む。この歪は主点を中心とし、主点から離れるほどその歪度合いは大きくなる。このようなレンズ歪を求めるには、図１３に示すように、三脚１３１上にカメラ１３２を固定し、市松模様などの３次元形状のターゲット１３３を奥行き方向に正確に移動させて、各奥行きでの３次元パターンを複数枚撮影し、得られた画像と既知の３次元座標値との対応関係からレンズ歪を求める手法が主流である。これに対して、図１４では、プリンタで印刷した程度の簡易な２次元のパターンを複数枚撮影した画像から、画像と既知の２次元座標値との対応関係からレンズ歪を算出する手法も存在する(例えば、非特許文献１参照。)。これらは精度よくレンズ歪を求めることができるが、いずれも既知の空間中の３次元座標値を与える必要があり、特に、図１３の場合はキャリブレーションするための高価で専用の器具が必要となる。
【０００４】
一方、市販カメラにはオートフォーカス機能が搭載されており、被写体との距離に合わせてピント調節を行い、ピントの合った画像の取得が可能である。レンズ歪を計測するときにはオートフォーカスで撮影された画像は各画像において焦点距離が微妙に調整される。この場合、レンズ歪と合わせて、このような焦点距離までも同時に推定する必要がある。
【非特許文献１】Ｚ．Ｚｈａｎｇ；“Ａ  Ｆｌｅｘｉｂｌｅ  Ｎｅｗ  Ｔｅｃｈｎｉｑｕｅ  ｆｏｒ  Ｃａｍｅｒａ  Ｃａｌｉｂｒａｔｉｏｎ”，ＩＥＥＥ  Ｔｒａｎｓ．Ｐａｔｔｅｒｎ  Ａｎａｌ．＆  Ｍａｃｈ．Ｉｎｔｅｌｌ，Ｖｏｌ．２２，Ｎｏ．１１，ｐｐ．１３３０−１３３４，２０００．
【発明の開示】
【発明が解決しようとする課題】
【０００５】
ディジタル・スチルカメラやハンディタイプのビデオカメラを使って手動操作により撮影した画像において、被写体との距離に応じてフォーカスを調整するというオートフォーカス機能により、画像ごとにピントが調整された画像が得られる。このようなオートフォーカスでは時系列画像においては、それぞれ焦点距離が変化するため、外界の物体の空間情報を獲得・復元するとき従来の因子分解法を利用することができず、ピント調整のために変化する焦点距離を同時に復元しなければならないという問題がある。
【０００６】
さらに、市販カメラにはレンズ歪があり、カメラで取得した画像では矩形形状が樽型に歪む。このようなレンズ歪があると、画像計測では正確な計測が困難であり、コンピュータビジョンにおいても３次元形状復元に影響を及ぼす。このために、事前にレンズ歪パラメータを求めておき、処理の対象とする画像に対してレンズ歪の補正を行い、レンズ歪がない画像にしておく必要がある。又は、画像座標値をレンズ歪に対して補正する幾何的変換を行うことでも同じことができるようになる。しかし、このようなカメラキャリブレーション作業は特殊な装置や既知の情報を必要としており、ユーザが簡易に操作してレンズ歪のない画像を取得することが容易ではなかった。
【０００７】
かかる事情に鑑みて本発明はなされたものであり、その目的は、上記課題を解決したレンズ歪パラメータ推定技術を提供することにある。
【課題を解決するための手段】
【０００８】
本発明では、レンズ歪及び主点であるレンズ歪パラメータ、及び焦点距離を推定するための投影モデルを
【０００９】
【数１】
2006195790.tif 000003
【００１０】
【数２】
2006195790.tif 000004
【００１１】
【数３】
2006195790.tif 000005
【００１２】
【数４】
2006195790.tif 000006
【００１３】
と仮定し、この投影モデルに基づいてカメラ運動、３次元形状（Ｘj，Ｙj，Ｚj）、レンズ歪係数ｋl、主点（Ｃｘ，Ｃｙ）、及び焦点距離ｆiを復元する。数式１〜数式４では、第ｉフレームの画像において観測される第ｊ番目の特徴点の画像座標値を（ｘij，ｙij）、そのときの焦点距離をｆi、ＸＹ並進運動を（Ｔｘi，Ｔｙi）、光軸周りの回転をθiとする。また、画像座標値（ｘij，ｙij）は主点（ｐｒｉｎｃｉｐｌｅ  ｐｏｉｎｔ）を中心とした画像面において観測された座標値であり、画像面の横縦のアスペクト比は１：１とする。
【００１４】
上記を全フレーム（ｉ＝１，２，・・・，Ｆ）、全特徴点（ｊ＝１，２，・・・，Ｐ）に対して行列表記すると、
【００１５】
【数５】
2006195790.tif 000007
【００１６】
となる。焦点距離ｆiはオートフォーカス機能等により変化してフレームごとに異なっている。また、画像から得られる観測座標値は（ｘij，ｙij）であり、数式５を満たすように平面運動と３次元形状を復元するにはレンズ歪パラメータ及び焦点距離を求めて、数式５左辺の行列要素を準備しなければならない。
【００１７】
本発明では、数式５の左辺の行列要素を段階的に得ながら、数式５の右辺に示す行列に分解できるように、反復的に平面運動投影型の因子分解法を使う(例えば、特開２００３−２７１９２５号参照。)。また、数式５の左辺の行列要素を得るには焦点距離とレンズ歪パラメータが必要であるが、焦点距離は時間的に可変であり、レンズ歪パラメータは時間的に一定であるために、同時に推定することは容易ではない。そこで、本発明では、焦点距離を推定する処理、及びレンズ歪パラメータを推定する処理に分けて、この問題に対処する。
【００１８】
すなわち、焦点距離、及び平面運動・３次元形状が与えられたとすると、数式５は以下のように考えることができる。
【００１９】
【数６】
2006195790.tif 000008
【００２０】
【数７】
2006195790.tif 000009
【００２１】
数式７は平面運動投影型の因子分解法(例えば、特開２００３−２７１９２５号参照。)で得られる平面運動及び３次元形状で算出される再投影座標値（ｘ’ij，ｙ’ij）であり、この座標値及び観測座標値（ｘij，ｙij）が与えられたときに、数式６の等式を満たすように、レンズ歪パラメータを求める問題に帰着する。この場合、全フレーム（ｉ＝１，２，・・・，Ｆ）、全特徴点（ｊ＝１，２，・・・，Ｐ）で加算して、
【００２２】
【数８】
2006195790.tif 000010
【００２３】
を評価関数として考え、数式８を最小化にするレンズ歪パラメータを求める。パラメータ最適化には、最急降下法、及びＰｏｗｅｌｌ法やＴａｙｌｏｒの微分補正法などの非線形最小二乗法、さらには、ニューラルネットワークを利用することが可能である。
【００２４】
一方、レンズ歪パラメータ及び平面運動・３次元形状が与えられたとすると、数式５は
【００２５】
【数９】
2006195790.tif 000011
【００２６】
として考えることができる。
【００２７】
上記を全フレーム（ｉ＝１，２，・・・，Ｆ）、全特徴点（ｊ＝１，２，・・・，Ｐ）に対して行列表記すると、
【００２８】
【数１０】
2006195790.tif 000012
【００２９】
となる。数式１０から、焦点距離は
【００３０】
【数１１】
2006195790.tif 000013
【００３１】
で得られるが、数式１０の左辺はランクが１であるため、数式１１の計算の前に、数式１０の右辺に示すｗijを要素とする行列を特異値分解し、ランク２以降の成分を除去しておく。これにより、雑音に対してロバストな焦点距離の推定が可能となる。
【００３２】
このように、焦点距離を推定する処理及びレンズ歪パラメータを推定する処理を反復的に実行する中で、数式５の左辺の行列要素をもつ行列（漸近行列）を得て、平面運動投影に基づく因子分解法を使用して、数式５の右辺の行列分解を行う。この反復において、平面運動への漸近を示す評価値として、漸近行列の要素と平面運動での再投影座標値との間の誤差、又は、漸近行列を特異値分解したときの特異値行列の特定の成分を使用することにより、漸近値の度合いで焦点距離を推定するか、レンズ歪パラメータを推定するかを決定し、漸近値を０へ近づける。この漸近が完了することで、焦点距離、及びレンズ歪パラメータを推定し、さらに、カメラ運動及び３次元形状を復元することができる。
【００３３】
上記説明は、焦点距離及びレンズ歪パラメータ双方を推定する方法であるが、は焦点距離推定を行わずにレンズ歪パラメータを求めることも可能であり、この場合には上記において焦点距離ｆiを無視した場合の投影モデルを使って上記と同様に扱うことができる。つまり、上記においてレンズ歪パラメータだけを反復的に推定する。
【発明の効果】
【００３４】
請求項１〜９に記載の発明によれば、カメラ手段を使って取得した時系列画像全般（移動手段を利用して撮影した車載画像、海上画像、空撮画像、屋内画像など）から、オートフォーカス機能が働いた状況でもレンズ歪パラメータを推定することが可能であり、さらに請求項３、４、７、及び８に記載の発明によれば焦点距離も推定することができる。
【００３５】
また、請求項１〜９に記載の発明はカメラが平面運動する状況に適用可能である。
【００３６】
また、請求項１〜９に記載の発明ではカメラ運動及び次元形状、及びカメラ内部パラメータをより高精度に復元、推定することが可能であり、これらの処理に使用する計算は、大半が線形演算で構成されるため、コンピュータ言語での実装が容易である。
【００３７】
また、請求項１〜９に記載の発明ではレンズ歪パラメータを推定できるので、取得した画像に幾何的補正を行い、レンズ歪を低減した透視投影像を取得できる。すなわち、レンズ歪パラメータを用いて各画素の座標値を変換し、その座標値に画素を埋め込む（画像を幾何的に補正する）ことにより、レンズ歪を低減した透視投影での画像を取得することが可能となる。
【発明を実施するための最良の形態】
【００３８】
(実施形態１)
図１はカメラを平行移動させて時系列画像を取得するシステムであり、回転式ハンドル１１の操作により平面レール１２上のワイヤ１３に固定したカメラ１４を平行移動させながら撮像を行うものである。
【００３９】
図２はカメラを回転移動させて時系列画像を取得するシステムであり、平面板２１のカメラを移動させる円２２に光軸と平行な支柱２３により固定されたカメラ２４を半径Ｒの円周上を平面運動させながら撮像を行うものである。
【００４０】
図３は図１、図２に記載のシステムにより時系列画像を取得したときのカメラと対象物との位置関係を表す図である。運動の中心は視点としており、視点を原点としたカメラ座標系ＸＹＺ，原点Ｏとした世界座標系ＸｗＹｗＺｗを設定する。この座標系において、カメラ運動として、Ｔi（Ｔｘi，Ｔｙi）の並進運動と光軸（Ｚ軸）周りのヨー回転θiとする。対象物の点Ｐj（Ｘj，Ｙj，Ｚj）はカメラにより画像面において投影中心（主点）を原点とした画像面へ投影され、画像座標値（ｘij，ｙij）として観測されるものとする。なお、説明の都合上、初期フレームでの視点とＯは一致しているものとし、光軸はＺｗ軸と平行関係にあり、θiはＸとＸｗ軸との成す角とするが、このことは一般性を損なうものではない。
【００４１】
図４は実施形態１におけるパラメータ推定装置の基本構成図であり、実施形態１を図４を用いて説明する。実施形態１のパラメータ推定装置は、時系列画像を蓄積した時系列画像データベース４１、各画像における特徴点の座標値を観測する特徴点観測部４２、画像座標値からなる行列データからカメラ運動と３次元形状を復元する因子分解処理部４３、レンズ歪パラメータを推定するパラメータ推定部４４、カメラ運動及び３次元形状、並びにレンズ歪パラメータが得られるまで処理を繰り返すための反復処理部４５から構成される。
【００４２】
時系列画像データベース４１には、ハードディスク、ＲＡＩＤ装置、ＣＤ−ＲＯＭなどの記録媒体を利用することが可能であり、ネットワークを介したリモートなデータ資源を利用することも可能である。また、図５はリアルタイムで処理を行う場合のパラメータ推定装置の基本構成図である。図５に示すように、パラメータ推定装置には必ずしも各データベース部などの記憶装置を必要としない。
【００４３】
このような構成のパラメータ推定装置において、まず、特徴点観測部４２が、対象物を撮影した時系列画像として時系列画像データベース４１からフレーム数Ｆの画像系列を取り出す。この取り出した画像系列において特徴点追跡を行う。特徴点は従来から用いられているような以下の手順により抽出する。
【００４４】
図６に示すように画像１の領域１において、
（１）各画素に対する２×２のヘッセ行列を求める。
（２）各点の３×３近傍領域において極大点かどうか判定し、極大点以外の点を削除する（ｎｏｎ−ｍａｘｉｍａ  ｓｕｐｐｒｅｓｓｉｏｎ）。
（３）得られた各点のヘッセ行列の固有値σl，σs（σs≦σl）を求め、σsがσp以上となる点を抽出する。
（４）抽出した点のσsの大きさの順にソートし、上位の点から順番にその点（ｐl）より上位の点（ｐh）がσd画素以内の距離に存在するかどうかを判定し、もし、存在する場合は下位の点ｐlを削除する。さらに、抽出した特徴点（ｊ＝１，２，・・・，Ｐ）をＫＬＴ法（Ｋａｎａｄｅ−Ｌｕｃａｓ−Ｔｏｍａｓｉ）により画像ｉ（ｉ＝１，２，・・・，Ｆ）にわたって追跡し、画像座標値（ｘij，ｙij）を観測する。
【００４５】
このようにして得られた特徴の画像座標値を
【００４６】
【数１２】
2006195790.tif 000014
【００４７】
に示す配列に並べた２Ｆ×Ｐの行列データ（行列データ[Ａ]）を用意する。
【００４８】
次に、因子分解処理部４３、パラメータ推定部４４、及び反復処理部４５の処理内容について、図７の処理フロー図にそって説明する。
【００４９】
因子分解処理部４３では、まず、行列［Ａ］の入力を行う(Ｓ７０１)。
【００５０】
漸近行列の生成(Ｓ７０２)では、レンズ歪係数ｋl、主点（Ｃｘ，Ｃｙ）、を
【００５１】
【数１３】
2006195790.tif 000015
【００５２】
のように初期化し、保持しておいた数式１２の２Ｆ×Ｐの行列データ［Ａ］の行列要素（ｘij，ｙij）から、数式１〜数式４にしたがってκij，（ｕij，ｖij）を計算する。
【００５３】
求めた要素を
【００５４】
【数１４】
2006195790.tif 000016
【００５５】
に示す配列に並べた２Ｆ×Ｐの行列を漸近行列［Ｂ］として保持しておく。
【００５６】
この漸近行列［Ｂ］を入力データとして特異値分解により
【００５７】
【数１５】
2006195790.tif 000017
【００５８】
に示す３つの行列［Ｕ］，［Ｗ］，［Ｖ］に行列分解する(Ｓ７０３)。
【００５９】
ここで、［Ｕ］は２Ｆ×Ｐサイズの行列、［Ｗ］はＰ×Ｐサイズの対角行列、［Ｖ］はＰ×Ｐサイズの行列である。
【００６０】
雑音除去(Ｓ７０４)では、
【００６１】
【数１６】
2006195790.tif 000018
【００６２】
の第二項に示すように、ランク４以降の各行列の成分を削除する。この削除のときは、行列［Ｕ］を取り出し、この行列の要素において第４から第Ｐ列目までを削除し、残りの成分からなる行列を保持し、行列［Ｗ］を取り出し、この行列の要素において第４から第Ｐ行目並びに第４から第Ｐ列目までを削除し、残りの成分からなる行列を保持し、行列［Ｖ］を取り出し、この行列の要素において第４から第Ｐ行目までを削除し、残りの成分からなる行列をそれぞれ保持する。
【００６３】
【数１７】
2006195790.tif 000019
【００６４】
ただし、特異値行列の（４，４）成分Ｗ44は保持しておく。
【００６５】
次に、第４から第Ｐ行目、及び第４から第Ｐ列目までを削除した行列［Ｗ］の対角要素の平方根をとった行列から、
【００６６】
【数１８】
2006195790.tif 000020
【００６７】
【数１９】
2006195790.tif 000021
【００６８】
に示す行列［Ｕ’］と行列［Ｖ’］を得る。
【００６９】
変換行列算出(Ｓ７０５)では、保持してある行列［Ｕ’］を取り出し、
【００７０】
【数２０】
2006195790.tif 000022
【００７１】
【数２１】
2006195790.tif 000023
【００７２】
【数２２】
2006195790.tif 000024
【００７３】
に示す連立する条件式における、下記数式２３で表されている対称行列［Ｃ］
【００７４】
【数２３】
2006195790.tif 000025
【００７５】
の各要素にかかる係数を計算する。これらの係数計算は、行列演算により容易に得られるものであり、この条件式を全フレームに対して計算する。
【００７６】
次に、全フレームの数式２０〜数式２２に示す連立する条件式に対して、最小二乗法などの数値計算を利用して３×３サイズの行列［Ｃ］の各要素を決定する。
【００７７】
求めた行列［Ｃ］を
【００７８】
【数２４】
2006195790.tif 000026
【００７９】
に示すように固有値分解する。ここで、固有値行列の平方根と固有値行列から、
【００８０】
【数２５】
2006195790.tif 000027
【００８１】
の行列［Ｃ’］を生成し、この行列要素を成分にもつ行列［Ｑ］を
【００８２】
【数２６】
2006195790.tif 000028
【００８３】
に従って算出する。
【００８４】
平面運動復元(Ｓ７０６)では、求めた行列［Ｑ］と、保持しておいた行列［Ｕ’］とから、
【００８５】
【数２７】
2006195790.tif 000029
【００８６】
の行列演算により行列［Ｍ’］を計算する。行列［Ｍ’］から各フレーム（第ｉフレーム）の行列要素（ｍix，ｎix）または（ｍiy，ｎiy）を取り出し、
【００８７】
【数２８】
2006195790.tif 000030
【００８８】
を使って、ヨー回転θiを復元する。
【００８９】
次に、行列［Ｍ’］から各フレーム（第ｉフレーム）の行列要素（Ｔiu，Ｔiv）を取り出す。この（Ｔiu，Ｔiv）から、
【００９０】
【数２９】
2006195790.tif 000031
【００９１】
を使って第ｉフレームにおけるユークリッド空間でのＸＹ並進運動（Ｔｘi，Ｔｙi）を計算する。
【００９２】
３次元情報復元(Ｓ７０７)では、先に保持しておいた行列［Ｖ’］と、変換行列算出で得られた行列［Ｑ］から、
【００９３】
【数３０】
2006195790.tif 000032
【００９４】
に示す行列演算を行い、行列［Ｓ’］を求める。
【００９５】
次に、行列［Ｓ’］の要素に対して、
【００９６】
【数３１】
2006195790.tif 000033
【００９７】
に示す変換を行い、これを要素とする行列を［Ｐ］とする。行列［Ｐ］の列ベクトルは、それぞれ第ｊ番目の特徴点のユークリッド空間での３次元座標値（Ｘj，Ｙj，Ｚj）になっている。
【００９８】
パラメータ推定部４４が行うパラメータ推定(Ｓ７０８)では、行列データの入力から観測座標値が与えられる。
【００９９】
一方、平面運動と３次元情報を復元した後で、数式７に示す再投影座標値を求めておく。
【０１００】
ここで、図８にパラメータ推定(Ｓ７０８)での処理フローを示す。図８において、まず、レンズ歪パラメータを初期化する(Ｓ８０１)。
【０１０１】
次に、パラメータ当てはめに必要な観測座標値（ｘij，ｙij）を行列［Ａ］入力から取り出し、求めておいた数式７の再投影座標値（ｘ’ij，ｙ’ij）を準備する(Ｓ８０２)。
【０１０２】
微分係数の計算は
【０１０３】
【数３２】
2006195790.tif 000034
【０１０４】
【数３３】
2006195790.tif 000035
【０１０５】
【数３４】
2006195790.tif 000036
【０１０６】
【数３５】
2006195790.tif 000037
【０１０７】
【数３６】
2006195790.tif 000038
【０１０８】
【数３７】
2006195790.tif 000039
【０１０９】
【数３８】
2006195790.tif 000040
【０１１０】
【数３９】
2006195790.tif 000041
【０１１１】
【数４０】
2006195790.tif 000042
【０１１２】
に示す値ａ11，ａ12，ａ13，ａ22，ａ23，ａ33，ｂ1，ｂ2，ｂ3を求める。
【０１１３】
続いて、補正量の計算(Ｓ８０３)では、数式３２〜数式４０の値を要素とする行列から
【０１１４】
【数４１】
2006195790.tif 000043
【０１１５】
に示す計算を実行し、補正量を△ｋ1，△Ｃｘ，△Ｃｙを得る。
【０１１６】
この補正量を
【０１１７】
【数４２】
2006195790.tif 000044
【０１１８】
に示すように加算したものを新たなレンズ歪パラメータとする(Ｓ８０４)。
【０１１９】
このとき、各補正量が許容値以下になる場合はそのまま処理を終了し、そのときのレンズ歪パラメータを保有し、図７でのパラメータ更新(Ｓ７０８)では、この保有したパラメータで更新する。一方、補正量が許容値以下とならない場合は、補正量を数式４２に示すように加算したものを新たなレンズ歪パラメータとして、数式３２〜数式４０に示す値ａ11，ａ12，ａ13，ａ22，ａ23，ａ33，ｂ1，ｂ2，ｂ3を求め直し、数式４１に示す計算を実行し、次の補正量を△ｋl，△Ｃｘ，△Ｃｙを得る。補正量を数式４２に示すように加算したものを新たなレンズ歪パラメータとし、許容値以下かどうかを判定する(Ｓ８０５)。この反復は許容値を満たすまで続ける。
【０１２０】
これにより、数式８の評価関数において、観測座標値（ｘij，ｙij）と再投影座標値（ｘ’ij，ｙ’ij）が与えられた場合の最適なレンズ歪パラメータが得られる。
【０１２１】
なお、実施形態１では放射型のレンズ歪を対象としているが、その他のレンズ歪モデルにも応用ができ、解析的に微分係数と補正量が得られれば、上記と同様に他のレンズ歪モデルにおいても本実施形態を利用することができる。
【０１２２】
反復処理部４５では、保持しておいたＷ44と１つ前の反復でのＷ44との絶対差分ΔＷを算出する(Ｓ７０９)。
【０１２３】
この差ΔＷが、所定の許容値以下に収束しているかどうかを判断し(Ｓ７１０)、収束している場合は処理を終了する。
【０１２４】
収束していない場合は、現時点で保有しているレンズ歪パラメータを使って数式１〜数式４に示す座標値（Ｕij，Ｖij）と係数κijを更新する(Ｓ７１１)。この係数更新に伴い、漸近行列の生成では数式１４により漸近行列［Ｂ］を更新し、特異値分解の処理に戻る。この繰り返しはΔＷが所定の許容値以下に収束するまで実行する。
【０１２５】
以上、実施形態１により、時系列画像の特徴点の画像座標値から、カメラ視点の運動：光軸周りの回転及びＸＹ軸方向の並進運動、並びに、物体形状を構成する３次元情報を復元し、同時に、レンズ歪パラメータを推定することができる。
【０１２６】
(実施形態２)
図９は実施形態２におけるパラメータ推定装置の基本構成図であり、実施形態２を図９を用いて説明する。実施形態２のパラメータ推定装置は、時系列画像を蓄積した時系列画像データベース９１、各画像における特徴点の座標値を観測する特徴点観測部９２、画像座標値からなる行列データからカメラ運動と３次元形状を復元する因子分解処理部９３、レンズ歪パラメータを推定するパラメータ推定部９４、焦点距離を推定する焦点距離推定部９５、カメラ運動と３次元形状、及びレンズ歪パラメータと焦点距離が得られるまで処理を繰り返すための反復処理部９６から構成される。
【０１２７】
時系列画像データベース９１には、ハードディスク、ＲＡＩＤ装置、ＣＤ−ＲＯＭなどの記録媒体を利用することが可能であり、ネットワークを介したリモートなデータ資源を利用する形態でもどちらでも構わない。また、図１０はリアルタイムで処理する場合の処理構成図である。図１０に示すように、必ずしも各データベース部などの記憶装置を必要としない。なお、本実施形態２は、実施形態１に対して焦点距離推定部９５の処理が異なっているため、以下では、この部分と関連する処理についてのみ記載する。
【０１２８】
ここで、因子分解処理部９３、パラメータ推定部９４、焦点距離推定部９５、及び反復処理部９６の処理内容について、図１１の処理フロー図にそって説明する。
【０１２９】
まず、因子分解部９３が、レンズ歪係数ｋl、主点（Ｃｘ，Ｃｙ）、及び焦点距離ｆiを
【０１３０】
【数４３】
2006195790.tif 000045
【０１３１】
のように初期化し(Ｓ１１０１)、漸近行列の生成において保持しておいた数式１２の２Ｆ×Ｐの行列データ[Ａ]の行列要素（ｘij，ｙij）から、数式１〜数式４にしたがってκij，（ｕij，ｖij）を計算する。このようにして求めた要素を
【０１３２】
【数４４】
2006195790.tif 000046
【０１３３】
に示すデータを求めておき、これを漸近行列［Ｂ］として保持しておく(Ｓ１１０２)。
【０１３４】
図１１での特異値分解(Ｓ１１０３)から、雑音除去(Ｓ１１０４)、変換行列算出(Ｓ１１０５)、平面運動復元(Ｓ１１０６)、３次元情報復元(Ｓ１１０７)の処理は実施実施形態１で説明しているため、ここでは処理モードの扱いについて説明する。
【０１３５】
初期化のとき、処理モードを「パラメータ推定」にしておく。これにより、最初の計算ではパラメータ推定が選択されて(Ｓ１１０８)、レンズ歪パラメータが推定される(Ｓ１１０９)。
【０１３６】
この処理は実施形態１で説明しているため処理内容を省略する。
【０１３７】
実施形態２では△Ｗが所定の許容値以下かどうかを判断し(Ｓ１１１２)、許容値以下となった場合、処理モードを「焦点距離推定」に切り替える(Ｓ１１１３)。これにより、レンズ歪パラメータを求める処理が完了すると、今度は焦点距離推定(Ｓ１１１０)の処理を開始する。
【０１３８】
焦点距離推定(Ｓ１１１０)では、まず、係数κijと再投影座標値（ｘ’ij，ｙ’ij）と（ｕij，ｖij）から数式９に示す値ｗijを求める。
【０１３９】
次に、数式１１に従い各フレームでの焦点距離ｆiを得る。この計算はｗij（ｊ＝１，２，・・・，Ｐ）の各ｉフレームに対する全特徴点による平均値となっている。
【０１４０】
図９の反復処理部９６では、保持しておいたＷ44と１つ前の反復でのＷ44との絶対差分ΔＷを算出し(Ｓ１１１１)、この差ΔＷが所定の許容値以下に収束しているか判断し(Ｓ１１１２)、収束していない場合は、先の計算で求めた焦点距離と現時点で保有しているレンズ歪パラメータを使って数式１〜数式４に示す座標値（ｕij，ｖij）と係数κijを更新する(Ｓ１１１３)。
【０１４１】
この係数更新に伴い、漸近行列の生成では数式４４により漸近行列［Ｂ］を更新する。この新たな漸近行列［Ｂ］を入力データとして特異値分解の処理に戻る。
【０１４２】
一方、差ΔＷが所定の許容値以下に収束した場合は、焦点距離推定の処理を停止し、処理モードを「パラメータ推定」に切り替える(Ｓ１１１５)。つまり、実施形態１で説明したパラメータ推定の処理を行い、その処理モードでの絶対差分ΔＷが収束するまでパラメータ推定を繰り返す。
【０１４３】
ただし、「パラメータ推定」と「焦点距離推定」との間の切り替えの際、切り替え回数を判断し(Ｓ１１１４)、切り替え回数がＮ回を超えた時点で全ての処理を停止し、その時点でのレンズ歪パラメータと焦点距離を出力する。
【０１４４】
以上、実施形態２により、時系列画像の特徴点の画像座標値から、カメラ視点の運動：光軸周りの回転とＸＹ軸方向の並進運動、及び物体形状を構成する３次元情報を復元し、同時に、レンズ歪パラメータ、焦点距離を推定することができる。
【０１４５】
なお、パラメータ推定装置は、例えば、パラメータ推定装置を構成するコンピュータ装置が有するＣＰＵによって実現され、特徴点観測機能、因子分解処理機能、パラメータ推定機能、焦点距離推定機能、反復処理機能などをアプリケーションプログラムとして搭載することができる。
【０１４６】
また、コンピュータ装置には上記説明した特徴点観測、因子分解処理、パラメータ推定、焦点距離推定、反復処理など行った処理結果や計算結果等のデータを内部メモリや外部記憶装置等に書き込み・読み出しができるようにしてもよい。
【０１４７】
また、本発明は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記録媒体を、システム又は装置に供給し、そのシステム又は装置のＣＰＵ(ＭＰＵ)が記録媒体に格納されたプログラムコードを読み出し実行することによっても実現できる。その場合、記憶媒体から読み出されたプログラムコード自体が上述した実施形態の機能を実現することになり、そのプログラムコードを記憶した記憶媒体、例えば、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＭＯ、ＨＤＤ等は本発明を構成する。
【図面の簡単な説明】
【０１４８】
【図１】カメラが平面レールの上で移動する場合を示す図。
【図２】カメラが平面板の上で移動する場合を示す図。
【図３】カメラ視点と３次元位置、及び、本発明の説明で使用する座標系を示す図。
【図４】実施形態１におけるパラメータ推定装置の構成図。
【図５】実施形態１におけるパラメータ推定装置の構成図。
【図６】特徴点抽出の説明図。
【図７】実施形態１における因子分解処理部、パラメータ推定部、及び反復処理部の処理フロー図。
【図８】パラメータ推定の処理フロー図。
【図９】実施形態２におけるパラメータ推定装置の構成図。
【図１０】実施形態２におけるパラメータ推定装置の構成図。
【図１１】実施形態２における因子分解処理部、パラメータ推定部、焦点距離推定部、及び反復処理部の処理フロー図。
【図１２】レンズ歪の説明図。
【図１３】３次元ターゲットを被写体としたレンズ歪測定の説明図。
【図１４】２次元ターゲットを被写体としたレンズ歪測定の説明図。
【符号の説明】
【０１４９】
１１…回転式ハンドル
１２…平面レール
１３…ワイヤ
１４…カメラ
２１…平面板
２２…カメラを移動させる円
２３…支柱
２４…カメラ
４１…時系列画像データベース
４２…特徴点観測部
４３…因子分解処理部
４４…パラメータ推定部
４５…反復処理部
５１…画像入力部
５２…特徴点観測部
５３…因子分解処理部
５４…パラメータ推定部
５５…反復処理部
９１…時系列画像データベース
９２…特徴点観測部
９３…因子分解処理部
９４…パラメータ推定部
９５…焦点距離推定部
９６…反復処理部
１０１…画像入力部
１０２…特徴点観測部
１０３…因子分解処理部
１０４…パラメータ推定部
１０５…焦点距離推定部
１０６…反復処理部
１２１…カメラ
１２２…被写体形状
１２３…投影像
１３１…三脚
１３２…カメラ
１３３…３次元ターゲット
【図１】
2006195790.tif 000047
【図２】
2006195790.tif 000048
【図３】
2006195790.tif 000049
【図４】
2006195790.tif 000050
【図５】
2006195790.tif 000051
【図６】
2006195790.tif 000052
【図７】
2006195790.tif 000053
【図８】
2006195790.tif 000054
【図９】
2006195790.tif 000055
【図１０】
2006195790.tif 000056
【図１１】
2006195790.tif 000057
【図１２】
2006195790.tif 000058
【図１３】
2006195790.tif 000059
【図１４】
2006195790.tif 000060