JP6482119B2

JP6482119B2 - 三次元計測法

Info

Publication number: JP6482119B2
Application number: JP2015047035A
Authority: JP
Inventors: 宣一吉川
Original assignee: Saitama University NUC
Current assignee: Saitama University NUC
Priority date: 2015-03-10
Filing date: 2015-03-10
Publication date: 2019-03-13
Anticipated expiration: 2035-03-10
Also published as: JP2016166806A

Description

本発明は、被測定物の形状を三次元計測する三次元計測法に関するものである。

従来、この種の三次元計測法としては、例えば、特許文献１に開示されたリアルタイム形状変形計測方法がある。この計測方法では、格子投影装置によって計測対象物体に格子パターンを投影する。投影する格子パターンは、１ピッチを明暗比１：２に分けた矩形波状の格子パターンを用い、その格子パターンが投影された物体をＣＣＤカメラで撮影する。格子パターンは、ＣＣＤカメラの１フレームの撮影時間である１／３０秒間に１／３ピッチで等速にシフトされ、ＣＣＤカメラは１／３０秒間に受けた光の強度を電荷として蓄え、それを輝度とした１つの格子画像に変換する。順次取得される各格子画像は、対応する第１〜第３フレームメモリに順に格納され、更新されていく。第１〜第３フレームメモリに格納された画像は、予め演算して求めた位相値に対する輝度値を格納する３次元位相算出テーブルによって即座に位相値に変換され、物体の位相分布画像としてリアルタイムに出力される。

特許第３５５４８１６号公報

格子パターン投影法で撮影された変形格子像の解析手法として、位相シフト法がある。位相シフト法は、投影する格子パターンの位相をシフトさせて３枚以上の変形格子像を撮影し、その変形格子像から位相情報を求めることで、計測対象物体の三次元情報を復元する手法である。この位相シフト法は、その撮影方式からステップ型と積分型とに分類することができる。

ステップ型は、計測対象物体に投影する格子パターンの位相をステップ移動させてシフトさせる手順、その位相シフトを一時停止させる手順、変形格子像を撮影する手順が順に行われ、その後に再び格子パターンの位相をシフトさせる手順以下の手順の流れが繰り返される方法である。一方、積分型は、計測対象物体に投影する格子パターンの位相を停止させずに連続的にシフトさせながら、変形格子像の撮影が行われる手法であり、連続フリンジスキャンと呼ばれることもある。

また、位相シフト法は、投影する格子パターンの位相シフト量の違いから、一定量位相シフト法と一般化位相シフト法とに分類することができる。一定量位相シフト法は、一般に用いられる位相シフト法で、例えばπ／２の一定量の位相を３ステップでシフトする３ステップ法が有名である。一方、一般化位相シフト法は、任意の位相シフト量で格子パターンの位相をシフトする方法である。

特許文献１に開示された上記従来の三次元計測法は、投影する格子パターンの撮影方式が１／３０秒間電荷を蓄える積分型で、位相シフト量がπ／３の一定量の、積分型・一定量位相シフト法である。投影する格子パターンの位相を一定量でシフトするこのような位相シフト法で、変形格子像を高速に撮影するためには、高速なデジタルカメラと位相シフタが必要になる。そして、格子投影装置の格子パターン表示，カメラの露光時間，シャッタータイミングなどを高精度に同期して制御する必要がある。これを実現するためには高性能な同期回路，デバイスと、十分な位相キャリブレーション作業が必要になる。

一方、任意の位相シフト量で格子パターンの位相をシフトする一般化位相シフト法は、位相シフト量の推定が必要になる。このため、従来、別の計測装置を使用して格子パターンの位相シフト量を計測したり、大量の変形格子像を用いて計算で位相シフト量を求めるなどの処理が必要になり、簡単に三次元計測が行えなかった。

本発明はこのような課題を解決するためになされたもので、
移動する格子パターンを被測定物に投影し、被測定物の形状に応じて変形する格子パターンを撮影して３枚の変形格子像を取得し、取得した３枚の変形格子像間における輝度値について平均二乗誤差を計算し、輝度値の位相分布のランダム性が成り立つときには輝度値の振動成分の空間平均が統計的に０になることに基づいて振動成分を表す項が０にされて三角関数で表された平均二乗誤差から、各変形格子像間の相対的な位相シフト量の絶対値を逆三角関数を使って求め、求めた各位相シフト量の和で表した評価関数が巡回差分位相条件を満たす各位相シフト量の符号および評価関数の解を解探索アルゴリズムを用いて求め、求めた各位相シフト量の符号および評価関数の解から符号を含めた位相シフト量の値を求め、位相シフト量のこの値から変形格子像の位相値を求めて被測定物の形状を三次元計測する三次元計測法を構成した。

本構成によれば、撮影方式が積分型で、投影する格子パターンの位相シフト量が任意の一般化位相シフト法により、各変形格子像間の相対的な位相シフト量の値を簡単に算出することができる。したがって、算出した位相シフト量の値から変形格子像の位相値を求めることで、被測定物の形状を簡単に三次元計測することが可能になる。一般化位相シフト法は、一定量位相シフト法におけるような正確な位相キャリブレーションを必要とせず、また、撮影装置と格子投影装置との間で高精度な同期制御をとる必要がない。このため、一定量位相シフト法を用いて三次元計測する際に発生する様々な問題を生じることなく、計測システムを容易に構築して簡単な計測手順で三次元計測することが可能になる。

また、本発明は、解探索アルゴリズムが、値が０となる評価関数を求めて、求めた評価関数を構成する各位相シフト量の係数値から各位相シフト量の符号を求め、求めた評価関数から解を見つけることを特徴とする。

本構成によれば、簡単なアルゴリズムで各位相シフト量の符号および評価関数の解を見つけることができる。このため、変形格子像間の相対的な位相シフト量の値をより簡単に算出することができ、被測定物の形状をより簡単に三次元計測することが可能になる。

本発明によれば、撮影方式が積分型で、投影する格子パターンの位相シフト量が任意の一般化位相シフト法により、変形格子像間の相対的な位相シフト量の値を簡単に算出して、被測定物の形状を簡単に三次元計測することが可能になる。

本発明の一実施の形態による三次元計測法を実施するための計測システムを示す概略構成図である。（ａ）は、一実施の形態による三次元計測法を実施するための計測システムにおける格子投影法の概略図、（ｂ）は、第１の変形例による格子投影法の概略図、（ｃ）は、第２の変形例による格子投影法の概略図である。（ａ），（ｂ），（ｃ）は、それぞれ、一実施の形態による三次元計測法に用いられる位相シフトされた３枚の変形格子像を示す図である。一実施の形態による三次元計測法における解探索アルゴリズムで用いられる評価関数ｆ_ｎ（κ）を表すグラフである。（ａ）は、一実施の形態による三次元計測法で求められた位相θを用いて三次元復元を行うことで得られる画像、（ｂ）は、最終的に計測される被測定物の三次元形状を示す画像である。

次に、本発明による三次元計測法を実施するための形態について、説明する。

図１は、本発明の一実施の形態による三次元計測法を実施するための計測システムＳを示す概略構成図である。

三次元計測が行われる対象の被測定物１には、等速で移動する格子パターン２が格子投影装置３によって投影される。本実施の形態では、格子パターン２には、輝度が正弦波状に変化する正弦波格子パターンが縞パターンとして理論的に用いられる。しかし、実際には、輝度キャリブレーションしてグラデーションのある正弦波近似縞パターンを作ったり、白黒二値の格子パターンをぼかして正弦波近似縞パターンを作ることなどで、明るい帯と暗い帯とが交互に繰り返される縞パターンを正弦波近似格子パターンとして得る。格子投影装置３から出射される格子パターン２は、軸４を中心に回転するミラー５で反射し、被測定物１に投影されることで、等速ｖで移動させられる。被測定物１に投影される格子パターン２はデジタルカメラ等の撮影装置６によって撮影され、撮影装置６には、被測定物１の形状に応じて変形する格子パターン２が変形格子像として取得される。この際、被測定物１は、撮影装置６からｚ方向に距離Ｌ_０離れたｘ−ｙ平面を基準平面として撮影される。

縞画像のアニメーションをプロジェクタを用いてＰＣ（パーソナルコンピュータ）の制御で被測定物１に投影すると、格子パターン２が等速ｖで移動しなかったり、撮影装置６との同期制御が必要になる。このため、本実施の形態では、図２（ａ）に示すように、格子投影装置３の内蔵光源３ａから出射される光を縞パターンが印刷されたフィルム３ｂに通し、軸４を中心に回転するミラー５で反射させることで、格子パターン２を等速ｖで移動させ続ける。

なお、格子パターン２を投影するデバイスは、上記のような写真フィルムを用いたスライド式のプロジェクタに限られることはなく、透過型や反射型の液晶パネルをＰＣで制御して投影する方式のプロジェクタ等であってもよく、固定パターンを投影することができるものであればよい。また、同図（ｂ）に示すように、格子投影装置３の内部で、円筒状に形成したフィルム３ｂを光源３ａを中心に等速で回転させることによっても、格子パターン２を等速ｖで移動させ続けることができる。また、同図（ｃ）に示すように、格子投影装置３の内部でフィルム３ｂを等速ｖで平行移動させることで、格子パターン２を等速ｖで移動させることもできる。図２に例示するこれらの構成によれば、格子投影装置３は、ＰＣ制御のプロジェクタである必要は無く、単なる縞パターンを投影する機能だけを備えた安価なものでよくなる。

撮影装置６によって撮影される被測定物１の変形格子像は位相シフト解析装置７へ出力され、以下のようにして、変形格子像の位相値、すなわち、被測定物１のｚ方向における高さが算出されて、被測定物１の三次元形状が計測される。位相シフト解析装置７は、並列処理に優れるＧＰＵ（グラフィック・プロセッシング・ユニット）を備えて構成される場合があり、高速動作するマイクロコンピュータによって以下の計算処理がリアルタイムに行われることで、撮影方式が積分型で、投影する格子パターン２の位相シフト量が任意の一般化位相シフト法により、被測定物１の三次元計測が高速にかつ高解像度で行われる。

位相シフト解析装置７は撮影装置６から連続する３枚の変形格子像、例えば、図３（ａ），（ｂ），（ｃ）に示すお面の３枚の変形格子像ａ，ｂ，ｃを入力する。これら各変形格子像ａ，ｂ，ｃにおける縞の輝度分布は余弦波形と見なすことができ、縞の輝度値Ｉφ（ｘ，ｙ）は空間（ｘ，ｙ）上に余弦波状に分布する。このとき、その波の、被測定物１により作られる位相をθ_０（ｘ，ｙ）とすると、輝度値Ｉφ（ｘ，ｙ）は次の（１）式のように表される。

ここで、Ａは背景輝度を表すバイアス項、Ｂは輝度の振幅を表すコントラスト項、ｆ_０は縞の周波数、Ｔは縞の周期、φは縞の位相シフト量である。θ_Ｒを次の（２）式

で定義し、θ＝θ_０＋θ_Ｒとおいて、座標系（ｘ，ｙ）を省略すると、（１）式に表される輝度値Ｉφ（ｘ，ｙ）は次の（３）式のように簡単化することができる。

変形格子像の位相φが連続的に変化している状態においては、撮影装置６の露光時間ｔの間にシフトする変形格子像の位相変化量をΔt とすると、位相シフト解析装置７に取得されて記録される変形格子像ａ，ｂ，ｃの輝度値Ｉi（ｘ，ｙ）はその強度を積分することで求めることができ、次の（４）式のように表される。

ここで、φi は変形格子像ｉの位相シフト量で、また、関数sinc(x)は次の（５）式のように表される。

（４）式は、元の（３）式の振動項にsinc関数の形の定数が掛算されたものである。ここで、Δt＝０ならばステップ型、Δt＞０ならば積分型の位相シフトになる。つまり、ステップ型と積分型の位相シフトは、次の（６）式に示す、実質的に同じ変形格子像の式で表すことができる。

上記の（４）式から、変形格子像ａ，ｂ，ｃの輝度値Ｉ_０，Ｉ_１，Ｉ_２は、変形格子像ａ，ｂ，ｃの位相シフト量をそれぞれφ_０，φ_１，φ_２とすると、次の（７）式のように表される。

これらの式からθを求めると、次の（８）式となる。

ただし、ΔＩpq＝Ｉq−Ｉp、Δφpq＝φq−φp（ｐ，ｑ＝０，１，２）である。この（８）式の導出は次のように行われる。ΔＩ_０１，ΔＩ_１２，ΔＩ_２０は、各輝度値Ｉ_０，Ｉ_１，Ｉ_２の差分をとり、複素形式で表すと、次の（９）式のようになる。

ただし、Ｃは次の（１０）式に表され、Ｃ^＊はＣの共役複素数である。

次に、（９）式のΔＩ_０１とΔＩ_２０の式からなる以下の（１１）式に表される連立方程式を考え、

この連立方程式をＣについて解くと、Ｃは次の（１２）式に表される。

ただし、Ｄは次の（１３）式に表される。

（１２）式において、初期位相φ_０をゼロとしても一般性は失われないので、Ｃは次の（１４）式に表すことができる。

この（１４）式に（１３）式に表されるＤを代入して整理すると、Ｃは次の（１５）式のように表される。

位相θは、（１５）式に表されるＣの位相から次の（１６）式のように求めることができる。

この結果、変形格子像ａ，ｂ，ｃのｚ方向の高さに相当する波の位相θは、（８）式に表されることになる。また、変形格子像ａ，ｂ，ｃの波の振幅（コントラスト）、つまり、縞の輝度値Ｉiの大きさは、（１５）式に表されるＣの絶対値を次の（１７）式で算出することで、求まる。

縞のコントラストが低い画像領域の位相θは信頼性の低い性質がある。このため、振幅が低い位相θの情報を除去することで、信頼性の高い三次元計測を行うことができる。

（８）式に表される位相θの値を求めるには、変形格子像ａ，ｂ，ｃ間の相対的な位相シフト量Δφpqを知る必要がある。本実施の形態では、次のようにして位相シフト量Δφpqの値を統計手法を用いて推定する。まず、（７）式に示される縞の輝度値Ｉ_０，Ｉ_１，Ｉ_２から、取得した２枚の変形格子像ａ，ｂ、ｂ，ｃ、ｃ，ａ間における輝度値Ｉｉの差ΔＩpqを次の（１８）式により求める。

次に、ΔＩpqの絶対値の二乗の空間平均を平均二乗誤差として次の（１９）式のように求める。

（１９）式における記号<>は、１フレームの画像全体について平均をとる演算子である。輝度値Ｉiの位相分布のランダム性が成り立つときには、空間に分布する輝度値Ｉiの振動成分の空間平均は統計的に０になる。本実施形態のように、縞状の格子パターン２を投影して撮影する縞パターン投影方式では、輝度値Ｉiの位相分布のランダム性が十分に成り立つことが分かっている。このことに基づき、（１９）式における輝度値Ｉiの振動成分を表す余弦関数の項<cos(２θ＋φp＋φq)>を０にすると、（１９）式は次の（２０）式で表される。

各変形格子像ａ，ｂ、ｂ，ｃ、ｃ，ａ間における格子パターン２の相対的な位相シフト量Δφpqの絶対値は、（２０）式に示す三角関数sinで表された平均二乗誤差から、次の（２１）式のように、逆三角関数arccosを使って求めることができる。

ただし、Ｅpqおよびκはそれぞれ次の（２２）式に表されるものとする。

ここで、逆三角関数arccosによって表される位相シフト量Δφpqは、０≦Δφpq＜πの位相シフト範囲に制限され、その符号はこれだけでは決定できない。このため、位相シフト範囲を−π≦Δφpq＜πに拡張するため、位相シフト量Δφpqの総和が２ｍπ（ｍは整数）となる次の（２３）式に表される巡回差分位相条件を使って、位相シフト量Δφpqの符号推定を行う。

ここで、位相シフト量Δφpqの符号を考慮した、次の（２４）式に表される評価関数ｆ_ｎ（κ）を考える。評価関数ｆ_ｎ（κ）は、変形格子像ａ，ｂ、ｂ，ｃ、ｃ，ａ間について（２１）式で求めた各位相シフト量Δφpqの和で表される。

ただし、係数ｃpq＝＋１または−１である。逆三角関数arccosによって（２１）式のように表される位相シフト量Δφpqは正の値であるので、係数ｃpqの値は位相シフト量Δφpqの符号に対応している。ｆ_ｎ（κ）＝２ｍπの制約条件式を成立させる評価関数ｆ_ｎ（κ）は巡回差分位相条件を満たし、その係数ｃpqの値は位相シフト量Δφpqの正しい符号である。

巡回差分位相条件を満たす各位相シフト量Δφpqの符号および評価関数ｆ_ｎ（κ）の解は、一般的な解探索アルゴリズムを用いて求めることができる。本実施の形態に用いる解探索アルゴリズムでは、評価関数ｆ_ｎ（κ）の値が０となる以下に示すゼロクロス条件で評価関数ｆ_ｎ（κ）を選択して、その評価関数ｆ_ｎ（κ）の各位相シフト量Δφpqに係る係数ｃpqの値、およびその評価関数ｆ_ｎ（κ）の解κ_０を見つける。

評価関数ｆ_ｎ（κ）は、３つの位相シフト量Δφpqの和で表され、各係数ｃpqの取り得る値の組み合わせによって次の８つの式ｆ_１（κ）〜ｆ_８（κ）で表される。

ｆ_５（κ）〜ｆ_８（κ）はｆ_１（κ）〜ｆ_４（κ）の符号反転なので、位相の回転方向が逆になっていることを示している。つまり、実質的にはｆ_１（κ）〜ｆ_４（κ）までの４通りを考えればよいことになる。これらの各式ｆ_１（κ）〜ｆ_８（κ）は、図４に示すグラフに表すことができる。同グラフの横軸はκの値、縦軸は評価関数ｆ_ｎ（κ）の値である。

同グラフから、ｆ_ｎ（κ）＝０のラインにクロスする評価関数ｆ_ｎ（κ）はｆ_１（κ）であることが分かる。したがって、この例では、各係数ｃpq（ｃ_０１，ｃ_１２，ｃ_２０）の値は、（２５）式のｆ_１（κ）の式からそれぞれ＋１となる。また、ｆ_ｎ（κ）＝２ｍπの制約条件を満たすκの解κ_０は、ｆ_１（κ）がｆ_ｎ（κ）＝０のラインにクロスするκの値になる。

最適な評価関数ｆ_ｎ（κ）は、詳細には次のように求めることができる。（２１）式に示す逆三角関数arccosより、次の（２６）式に示す不等式が成り立つ。

したがって、κの範囲は次の（２７）式に示す範囲となる。

ただし、Ｅmaxは次の（２８）式の大括弧｛｝中の最も大きなＥpqとする。

このκの範囲で、巡回差分位相条件を満たす評価関数ｆ_ｎ（κ）を調べる。ここで、αを０＜α＜κ_０＜２／Ｅmaxの関係を満たす値とする（実際の計算ではαは０に近い値とすればよい）。また、ｆ_１（κ）の場合は評価関数ｆ_ｎ（κ）＝ｆ_１（κ）−２π、ｆ_２（κ），ｆ_３（κ），ｆ_４（κ）の場合は評価関数ｆ_ｎ（κ）＝ｆ_ｎ（κ）、ｎ＝２，３，４の式を用いる。評価関数ｆ_ｎ（κ）がｆ_ｎ（κ）＝０の軸と交差する場合、０の前後で評価関数ｆ_ｎ（κ）の符号が変わるので、ｆ_ｎ（α）・ｆ_ｎ（２／Ｅmax）＜０となる関数が最適な評価関数ｆ_ｎ（κ）となる。最適な評価関数ｆ_ｎ（κ）が決定したら、その解κ_０は２分法のような単純なアルゴリズムによって求めることができる。また、各位相シフト量Δφpqの係数ｃpqの値は（２５）式のｆ_１（κ）〜ｆ_８（κ）の式を参照することで求めることができる。

符号を含めた位相シフト量Δφpqの値は、このように求めた各位相シフト量Δφpqの係数ｃpqの値および評価関数ｆ_ｎ（κ）の解κ_０から、次の（２９）式によって求めることができる。

前述したように、この符号を含めた位相シフト量Δφpqの値を（８）式に代入することで、変形格子像ａ，ｂ，ｃのｚ方向の高さに相当する波の位相θが求められる。この位相θから一般的な方法を用いて三次元復元を行うことで図５（ａ）に示す画像が得られ、最終的に同図（ｂ）に示す被測定物１の三次元形状が計測される。

このような本実施の形態による三次元計測法によれば、撮影方式が積分型で、投影する格子パターン２の位相シフト量Δφpqが任意の一般化位相シフト法により、各変形格子像ａ，ｂ，ｃ間の相対的な位相シフト量Δφpqの値を上記のように簡単に算出することができる。したがって、算出した位相シフト量Δφpqの値から（８）式を用いて変形格子像ａ，ｂ，ｃの位相θの値を求めることで、被測定物１の形状を簡単に三次元計測することが可能になる。一般化位相シフト法は、一定量位相シフト法におけるような正確な位相キャリブレーションを必要とせず、また、撮影装置６と格子投影装置３との間で高精度な同期制御をとる必要がない。さらに、位相シフトの安定性が良い。このため、一定量位相シフト法を用いて三次元計測する際に発生する様々な問題を生じることなく、ＰＣシステム等で計測システムＳを容易に構築して、簡単な計測手順で短時間に三次元計測することが可能になる。

また、本実施の形態による三次元計測法では、値が０となる評価関数ｆ_ｎ（κ）を求めて、求めた評価関数ｆ_ｎ（κ）を構成する各位相シフト量Δφpqの係数ｃpqの値から各位相シフト量Δφpqの符号を求め、また、求めた評価関数ｆ_ｎ（κ）から解κ_０を見つける解探索アルゴリズムが採用されている。本構成によれば、簡単なアルゴリズムで各位相シフト量Δφpqの符号および評価関数ｆ_ｎ（κ）の解κ_０を見つけることができる。このため、変形格子像ａ，ｂ，ｃ間における格子パターン２の相対的な位相シフト量Δφpqの値をより簡単に算出することができ、被測定物１の形状をより簡単に三次元計測することが可能になる。

上記の本発明による三次元計測法は、デジタル３Ｄ（三次元）アーカイブのための３Ｄスキャナや、腹腔鏡下での手術における三次元計測、義肢・装具製作，義歯製作における三次元計測、ロボットビジョンや人間の運動解析における三次元計測、ヒューマンインターフェース、生体認証等のセキュリティ技術、リバースエンジニアリング、３ＤＣＧ（コンピュータ・グラフィック）や３Ｄプリンタのための３Ｄデータ取得技術等に利用することが可能である。

Ｓ…計測システム
１…被測定物
２…格子パターン
３…格子投影装置
３ａ…光源
３ｂ…フィルム
４…軸
５…ミラー
６…撮影装置
７…位相シフト解析装置

Claims

移動する格子パターンを被測定物に投影し、被測定物の形状に応じて変形する前記格子パターンを撮影して３枚の変形格子像を取得し、取得した３枚の前記変形格子像間における輝度値について平均二乗誤差を計算し、前記輝度値の位相分布のランダム性が成り立つときには前記輝度値の振動成分の空間平均が統計的に０になることに基づいて前記振動成分を表す項が０にされて三角関数で表された前記平均二乗誤差から、各前記変形格子像間の相対的な位相シフト量の絶対値を逆三角関数を使って求め、求めた各前記位相シフト量の和で表した評価関数が巡回差分位相条件を満たす各前記位相シフト量の符号および前記評価関数の解を解探索アルゴリズムを用いて求め、求めた各前記位相シフト量の符号および前記評価関数の解から符号を含めた前記位相シフト量の値を求め、前記位相シフト量のこの値から前記変形格子像の位相値を求めて被測定物の形状を三次元計測する三次元計測法。
前記解探索アルゴリズムは、値が０となる前記評価関数を求めて、求めた前記評価関数を構成する各前記位相シフト量の係数値から各前記位相シフト量の符号を求め、求めた前記評価関数から前記解を見つけることを特徴とする請求項１に記載の三次元計測法。