JP4775540B2 - 撮影画像における歪曲収差補正方法 - Google Patents

Description

本発明は、撮影画像における歪曲収差補正方法に関する。

近年、デジタルカメラなどの撮影機器の性能向上により、高画質の画像が比較的簡単に得られるようになっており、例えば、橋梁などの構造物をデジタルカメラで撮影し、その撮影画像から構造物の形状、大きさなどの各種寸法を、非接触にて且つ高精度に計測し得るシステムが開発されている。

しかし、デジタルカメラなどの光学レンズを用いたシステムでは、レンズの歪曲収差による計測誤差が含まれてしまう。
このため、計測対象物に予めターゲットマークが設けられたキャリブレーション板を取り付けておき、このキャリブレーション板を撮影することにより、歪曲収差を補正するようにした計測方法が提案されている（例えば、特許文献１参照）。
特開２００１−１３３２２５

上述したように、計測対象物に取り付けられたキャリブレーション板、例えば格子パターンを撮影することにより歪曲収差を補正する計測方法によると、カメラで撮影した際に、位置振れや回転振れが生じた場合、カメラの位置振れや回転振れによる誤差分を分離することができず、したがって歪曲収差を補正した場合でも、その計測結果に誤差が発生し、高精度な計測を行うことができないという問題があった。勿論、格子パターンが正しく配置されていない場合にも、同様の問題が生じる。

そこで、本発明は、カメラに振れが生じた場合、または格子パターンの位置がずれている場合でも、これらに起因して発生する誤差についても補正し得る撮影画像における歪曲収差補正方法を提供することを目的とする。

上記課題を解決するため、本発明の請求項１に係る撮影画像における歪曲収差補正方法は、
二次元の格子状パターンをレンズを介して撮影して得られた輝度値を示す撮影画像にフーリエ変換を施して周波数スペクトルを得るステップと、
フィルタリング処理により、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去するステップと、
上記１次調和波にフーリエ逆変換を施して位相値を求めるステップと、
上記ステップにて所定範囲に亘って求められた位相値の位相接続を行うステップと、
上記位相接続により得られた位相値（δ _ｘ ，δ _ｙ ）を下記（Ａ）式および（Ｂ）式に代入して格子状パターンの変形量（β_ｘ，β_ｙ）を求めるステップと、
上記格子状パターンの回転移動および並行移動を考慮した平面の式で表されるずれ量算出式に当該格子状パターンに生じる下記（Ｃ）式および（Ｄ）式に基づく収差モデル量（α_ｘ，α_ｙ）を加えてなる下記（Ｅ）式および（Ｆ）式で示される変形量算出式に上記変形量を代入し、最小二乗法を用いて当該変形量算出式の少なくとも収差モデル量を表す係数を決定するステップと、
上記決定された収差モデル量の係数を用いて歪曲収差を補正するステップと
を具備した補正方法である。
β _ｘ ＝−δ _ｘ ／２πω ・・・（Ａ）
β _ｙ ＝−δ _ｙ ／２πω ・・・（Ｂ）
但し、ωは変形前の格子の周波数である。
α _ｘ ＝（ｇ _１ ＋ｇ _３ ）ｘ ^２ ＋ｇ _４ ｘｙ＋ｇ _１ ｙ ^２ ＋ｋ _１ ｘ（ｘ ^２ ＋ｙ ^２ ）・・（Ｃ）
α _ｙ ＝ｇ _２ ｘ ^２ ＋ｇ _３ ｘｙ＋（ｇ _２ ＋ｇ _４ ）ｙ ^２ ＋＋ｋ _１ ｙ（ｘ ^２ ＋ｙ ^２ ）・・（Ｄ）
β _ｘ ＝ａ _１ ｘ＋ａ _２ ｙ＋ａ _３ ＋α _ｘ ・・・（Ｅ）
β _ｙ ＝ａ _４ ｘ＋ａ _５ ｙ＋ａ _６ ＋α _ｙ ・・・（Ｆ）

また、請求項２に係る撮影画像における歪曲収差補正方法は、
二次元の格子状パターンをレンズを介して撮影して得られた輝度値を示す撮影画像にフーリエ変換を施して周波数スペクトルを得るステップと、
フィルタリング処理により、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去するステップと、
上記１次調和波を、収差が無い格子状パターンの周波数スペクトル分だけ原点側にシフトさせた後、フーリエ逆変換を施して位相値を求めるステップと、
上記ステップにて所定範囲に亘って求められた位相値の位相接続を行うステップと、
上記位相接続により得られた位相値（δ _ｘ ，δ _ｙ ）を下記（Ａ）式および（Ｂ）式に代入して格子状パターンの変形量（β_ｘ，β_ｙ）を求めるステップと、
上記格子状パターンの回転移動および並行移動を考慮した平面の式で表されるずれ量算出式に当該格子状パターンに生じる下記（Ｃ）式および（Ｄ）式に基づく収差モデル量（α_ｘ，α_ｙ）を加えてなる下記（Ｅ）式および（Ｆ）式で示される変形量算出式に上記変形量を代入し、最小二乗法を用いて当該変形量算出式の少なくとも収差モデル量を表す係数を決定するステップと、
上記決定された収差モデル量の係数を用いて歪曲収差を補正するステップと
を具備した補正方法である。
β _ｘ ＝−δ _ｘ ／２πω ・・・（Ａ）
β _ｙ ＝−δ _ｙ ／２πω ・・・（Ｂ）
但し、ωは変形前の格子の周波数である。
α _ｘ ＝（ｇ _１ ＋ｇ _３ ）ｘ ^２ ＋ｇ _４ ｘｙ＋ｇ _１ ｙ ^２ ＋ｋ _１ ｘ（ｘ ^２ ＋ｙ ^２ ）・・（Ｃ）
α _ｙ ＝ｇ _２ ｘ ^２ ＋ｇ _３ ｘｙ＋（ｇ _２ ＋ｇ _４ ）ｙ ^２ ＋＋ｋ _１ ｙ（ｘ ^２ ＋ｙ ^２ ）・・（Ｄ）
β _ｘ ＝ａ _１ ｘ＋ａ _２ ｙ＋ａ _３ ＋α _ｘ ・・・（Ｅ）
β _ｙ ＝ａ _４ ｘ＋ａ _５ ｙ＋ａ _６ ＋α _ｙ ・・・（Ｆ）

上記各歪曲収差補正方法によると、撮影画像による格子状パターンの位相値を求める際に、収差モデル量に格子状パターンの回転量および並行移動量を表す平面の式を加えた変形量算出式を用いるとともに、この変形量算出式の少なくとも収差モデル量の係数を、実際に求められた変形量を用いて最小二乗法により決定するようにしたので、格子状パターンを用いて歪曲収差を補正する際に、格子状パターンがずれている場合でも、そのずれ量を考慮した補正を行うことができるので、精度良く歪曲収差の補正を行うことができる。

［実施の形態］
以下、本発明の実施の形態に係る撮影画像における歪曲収差補正方法を、図面に基づき説明する。

本実施の形態においては、例えばデジタルカメラなどの撮影機器（以下、カメラという）を用いて、すなわちレンズを介して撮影対象物を撮影するとともに、このカメラで撮影された撮影対象物の撮影画像に基づき当該撮影対象物の形状、大きさなどの表面に関わる各種寸法を計測する際に、レンズの歪曲収差による変位分が誤差として含まれているため、この歪曲収差による変位分を除去するための歪曲収差補正方法について説明する。

本実施の形態に係るに歪曲収差補正方法は、基準となる二次元格子パターン（格子状パターンの一例であり、参照用パターンとも言える）を撮影し、この格子パターンの撮影画像における輝度値の信号パターン（信号周波数）にフーリエ変換を施して位相値を求め、この位相値に基づき変位量（変形量）を求め、この変位量を用いて当該格子パターンの歪曲収差を補正する際に、格子パターン自体が正しく配置されていないことによる例えば回転移動および並行移動によるずれ量についても補正し得る方法である。

以下の説明においては、歪曲収差により周辺部が変形した格子パターンの位相値をフーリエ変換法にて求める方法を説明した後、格子パターンの平面上でのずれ量を補正する方法について説明する。

まず、格子パターンの撮影画像における位相値をフーリエ変換を用いて求める方法について説明する。
なお、歪曲収差とは、例えば直線である物体が曲線となって撮影される現象を言い、収差の中心（多くの場合は、画像の中心であるが、そうでない場合もある）からの距離を用いてモデル化することができる。

例えば、図１に示すように、本来、Ａ点の位置に写るべき点が歪曲収差によりＢ点に写ったとする。
このＡ点とＢ点の距離（収差の量）は、半径方向（放射方向）歪曲収差α_ｒと円周方向歪曲収差α_θとに分解することができる。

そして、円周方向歪曲収差α_θは半径方向歪曲収差α_ｒに比べて十分に小さく、例えば円周方向歪曲収差α_θを無視した場合、半径方向歪曲収差α_ｒは、収差の中心からの距離ｒの関数として、下記（１）式にて表される。

ここで、ｋ_１，ｋ_２，ｋ_３，・・・は収差の量を表す係数である。多くの場合、高次項を無視できるので、（１）式は下記（２）式にて近似することができる。

すなわち、画像上の歪曲収差の量（幾何学的なずれ量）は収差の中心からの距離ｒの３乗に比例し、その係数ｋ_１を決定することで歪曲収差の量を知ることができ、これに基づき補正を行うことができる。

なお、半径方向歪曲収差α_ｒをｘ方向およびｙ方向に分解した場合、すなわち直交座標系における歪曲収差の量（収差モデル量）は、下記（３）式で表される。

ここで、θ＝arctan（ｘ，ｙ）である。

ところで、円周方向歪曲収差α_θについても考慮する必要がある場合、当該歪曲収差のモデルは、各係数ｇ_１，ｇ_２，ｇ_３，ｇ_４，ｋ_１を用いて下記（４）式にて表すことができる。この（４）式については公知である（例えば、「Weng, J., Cohen, P., and Herniou, M., Camera Calibration with Distortion Models and Accuracy Evaluation, IEEE Trans. Pattern Anal. Mach. Intell., 14(10), 965-980 (1992)」参照）。これらの各係数ｇ_１，ｇ_２，ｇ_３，ｇ_４，ｋ_１を決定することで歪曲収差の量を知ることができ、これに基づき補正を行うことができる。

以下、この円周方向歪曲収差α_θを考慮した場合について説明する。

すなわち、（４）式における各係数ｇ_１，ｇ_２，ｇ_３，ｇ_４，ｋ_１を決定する必要がある。

また、α_ｘ，α_ｙを用いると、収差により歪んだ点の座標ｘ′，ｙ′は、収差を含まない場合の座標ｘ，ｙに対し、下記（５）式にて表すことができる。

ところで、フーリエ変換を用いて、格子パターンにおける撮影画像の位相値を決定し得ることは一般的に知られている。

例えば、図２に示す変形していない格子パターンに歪曲収差が生じると、図３に示すような変形した画像として撮影される。
この撮影された輝度値の格子パターンにフーリエ変換を施して周波数スペクトルを得るとともに、この周波数スペクトルの１次調和波だけを残すようにフィルタリング処理（つまり、ノイズの除去処理）を行い、次にその１次調和波のスペクトルを、変形していない（歪曲収差が無い場合の）格子パターンにおける周波数スペクトル分だけ中心側（原点側）にシフトさせた後、フーリエ逆変換を行うことにより、格子パターンのｘ方向およびｙ方向での変位（変形）を表す位相分布ψ_ｘ，ψ_ｙを求める。

このようにして求められた位相分布については、格子パターンの一対の格子線間についてのものであり、これらを連続して表す必要がある。
したがって、これらの位相分布の位相接続が行われて、ｘ方向およびｙ方向での変形に係る位相値δ_ｘおよびδ_ｙが求められる。

上述した方法を、図面にて示すと、図４〜図６のようになる。
図４は歪曲収差の影響により変形して画像に記録された格子パターンにフーリエ変換を施した後、その周波数スペクトルの１次調和波成分（ｘ方向およびｙ方向それぞれ）を残し、他の成分を除去するフィルタリングを行った結果を示しており、図５（ａ）および（ｂ）は、各１次調和波成分をそれぞれ原点（図４の中心）側にシフトさせて、フーリエ逆変換を行った結果の位相分布を示している。図５の白い箇所は位相値πradを表し、黒い箇所は位相値−πradを表し、その間のグレーはその間の値を表している。

これらの位相分布の位相接続を行うと、図６（ａ）および（ｂ）に示すように、連続した位相分布を得ることができる。図６（ａ）および（ｂ）において、白は２πradを表し、黒は−２πradを表している。

この位相分布は格子パターンの変形状態、すなわち歪曲収差の分布状態を表しており、変形量（収差量）については、位相値から下記（６）式に基づき計算することができる。

ここで、δ_ｘおよびδ_ｙは位相接続後のｘ方向およびｙ方向における位相値であり、ωは変形前の格子の周波数である。

ところで、格子パターンをカメラの水平姿勢および垂直姿勢（鉛直姿勢でもある）に対して傾けずに設置することは実際には困難であるとともに、１次調和波成分を正確に原点側にシフトさせることが難しく、したがって得られた変形量β_ｘおよびβ_ｙについては、正確な収差を表しているとは言えず、格子パターンのずれ量（回転量および並行移動量）も含んでいる。

ところで、この格子パターンのずれ量、すなわち格子パターンの傾き（回転量）および周波数スペクトルのシフト（並行移動量）の誤差により生じる位相分布については、平面を表す式でもって近似し得るという知見を得た。

すなわち、図７の（ａ）に示す基準となる２次元の格子パターンが、（ｂ）に示すように、僅かに反時計回りに傾いていた場合、フーリエ変換により得られた格子パターンの変位を表す位相値には、格子パターンの回転量が含まれる。

回転する前の面内の点、すなわち図７（ａ）上の点の座標を（ｘ，ｙ）、回転後の座標、すなわち図７（ｂ）上の（ｘ，ｙ）に対応する点の座標を（ｘ_ｒ，ｙ_ｒ）とする｛点（ｘ，ｙ）が剛体回転運動により（ｘ_ｒ，ｙ_ｒ）に移動したとする｝と、点（ｘ，ｙ）と点（ｘ_ｒ，ｙ_ｒ）との関係は回転行列を用いて、下記（７）式のように表すことができる。

ここで、θは剛体回転角度である。この場合の各点のｘ方向変位は、ｕ_ｘ＝ｘ_２−ｘ_１、ｙ方向変位は、ｕ_ｙ＝ｙ_２−ｙ_１であるので、下記（８）式が得られる。

θは回転角度であるので、撮影した特定の基準となる格子パターンの画像に対して定数となる。

ところで、周波数スペクトルの１次調和波成分を原点側にシフトする際のシフト量に誤差が含まれていた場合には、基準となる変形前の格子パターンの周波数が実際とは異なることになる。

それは、変形前の格子パターンの間隔（ピッチ）が異なることを意味している。例えば、実際には、図７（ａ）の格子パターンを用いているにも拘わらず、ｘ方向の周波数スペクトルの１次調和波成分を正確に原点側にシフトできなかった場合には、図８に示すように、ｘ方向に伸びた（または縮んだ）格子パターンを用いていることになってしまう。この場合、ｘ方向単軸の変形と等価な位相分布が得られる。

すなわち、ｘ方向の変位は、下記（９）式にて表される。

ここで、ε_ｘはｘ方向の歪であり、この場合には画像上で一様な値となる。ｙ方向についても同様（ｕ_ｙ＝ε_ｙｙ）である。

したがって、剛体回転変位｛（８）式｝と１次調和波成分のシフト誤差｛（９）式｝とを考慮すると、ｘ方向およびｙ方向ともに、下記（１０）式のように平面の式にて表すことができる。

ここで、定数ｃは理論的には必ずしも必要ではないが、計測データを取り扱うということで組み込んでいる。例えば、フーリエ変換による位相解析の場合の座標の原点と、その後のデータ解析（最小二乗法など）の原点（収差の中心）とが異なる場合に、定数項ｃが必要となるからである。

したがって、格子パターンの収差量（α_ｘ，α_ｙ）とずれ量（ｕ_ｘ，ｕ_ｙ）を考慮した変形量β_ｘおよびβ_ｙは下記（１１）式および（１２）式にて近似することがきる。

なお、半径方向成分（ｒ方向）を考えた場合、その変形量β_ｒについては、下記（１３）式にて近似することができるが、ここでは、β_ｘおよびβ_ｙを用いて説明する。

そして、（６）式により得られた値の分布から、（１１）〜（１２）式を用いて、各係数を決定することができる。

このとき、（６）式により得られるβ_ｘ，β_ｙの値は収差を含む点ｘ′，ｙ′の関数として得られる。
一方、（４）式および（１１）〜（１２）式の収差のモデルは当該収差を含まない点（ｘ，ｙ）の関数として与えられている。そのため、各係数を直接決定することはできない。そこで、最小二乗法と繰返し演算とを組み合わせて、各係数を決定する。

まず、収差を含む座標値（ｘ′，ｙ′）を（ｘ，ｙ）に初期値として代入し、最小二乗法により各係数の近似値を決定する。すなわち、円周方向歪曲収差を考慮し、ｘ方向およびｙ方向の格子の位相分布を利用すると、下記（１４）式にて計算することができる。勿論、最小二乗法を用いる場合には、複数の画素について計算が行われる。

ここで、

以上の計算により得られた各係数の近似値を用い、収差を含まない座標（ｘ，ｙ）の近似値を計算する。その座標値（ｘ，ｙ）を用い、再び、同様の最小二乗法により各係数を計算する。これを各係数の値が収束するまで繰り返すことにより、各係数を決定することができる。

多くの場合、収差の中心は画像の中心でよいと考えられるが、収差の中心が分からない場合には、収差の中心の座標を未知数として非線形最小二乗法を用いればよい。このとき、ｒおよびθと収差の中心座標（ｘ_０，ｙ_０）との関係は下記（１８）式で表される。

以上の方法により、収差モデル量すなわち収差量を決定した後、その値を用いて撮影画像が補正される。撮影画像の歪みを補正する場合、例えば画像の輝度値補間法などが用いられる。

ここで、デジタル画像相関法を用いた変位測定結果に対して、本発明のフーリエ変換を用いた収差補正方法を適用し、その有効性を検証した結果について説明しておく。
図９は平板をｘ方向に移動（剛体並進変位）させた場合の変位をデジタル画像相関法により求めた等変位線図（変位分布）を示している。剛体並進変位であるので、変位は場所により異なる事は無く一様になるはずである。しかしながら、歪曲収差の影響により図のように変位に分布が現れている。これらの図のレジェンドの単位は画素（ピクセル）である。図９（ａ）を見ると、画像の中心での変位が９１．５画素であるのに対し、画像の周囲では９６．５画素と５画素もの差が生じている。なお、図９（ｂ）はｙ方向の変位を表している。

図１０は本発明に係る収差補正方法を用いて補正した変位分布（ａ）と補正前の変位分布（ｂ）との比較を示す図で、収差の除去が行われていることが分かる。この例では、収差のモデルとして（４）式を用いたが、（３）式を用いた場合においても同様の結果が得られた。

図１１は平板に単軸引張負荷を与えた際の板の中心付近における変位分布を表しており、（ａ）はｘ方向変位を、（ｂ）はｙ方向変位を示している。この図の収差補正前の分布では丸みを帯びた等変位線が現れている。一方、収差補正後の図１２ではほぼ直線の等変位線が現れており｛（ａ）はｘ方向変位で、（ｂ）はｙ方向変位を示している｝、妥当な結果が得られていると言える。なお、この図で等変位線が斜めに現れているのは、引張負荷を与えた試験片が負荷後に回転（剛体回転）を生じているためである。

ここで、上述した歪曲収差補正方法をステップ様式にて記載しておく。
すなわち、この歪曲収差補正方法は、二次元の格子状パターンをレンズを介して撮影して得られた輝度値を示す撮影画像にフーリエ変換を施して周波数スペクトルを得るステップと、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去するステップと、上記１次調和波を、収差が無い格子状パターンの周波数スペクトル分だけ原点側にシフトさせた後、フーリエ逆変換を施して位相値を求めるステップと、上記ステップにて所定範囲に亘って求められた位相値の位相接続を行うステップと、上記位相接続により得られた位相値に基づき格子状パターンの変形量を求めるステップと、上記格子状パターンの回転移動および並行移動を考慮した平面の式で表されるずれ量算出式に当該格子状パターンに生じる収差モデル量を加えてなる変形量算出式に上記変形量を代入し、最小二乗法を用いて当該変形量算出式の少なくとも収差モデル量を表す係数を決定するステップと、上記決定された収差モデル量の係数を用いて歪曲収差の補正を行うステップとが具備されている。

また、上記収差モデル量として半径方向の歪曲収差、若しくは収差モデル量として半径方向および円周方向の歪曲収差を考慮したものである。
さらに、詳しくは、二次元の格子状パターンをレンズを介して撮影して得られた輝度値を示す撮影画像にフーリエ変換を施して周波数スペクトルを得るステップと、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去するステップと、上記１次調和波を、収差が無い格子状パターンの周波数スペクトル分だけ原点側にシフトさせた後、フーリエ逆変換を施して位相値を求めるステップと、上記ステップにて所定範囲に亘って求められた位相値の位相接続を行うステップと、上記位相接続により得られた位相値に基づき格子状パターンの変形量（β_ｘ，β_ｙ）を求めるステップと、上記格子状パターンの回転移動および並行移動を考慮した平面の式で表されるずれ量算出式に当該格子状パターンに生じる収差モデル量（α_ｘ，α_ｙ）を加えてなる下記（Ａ）式および（Ｂ）式で示される変形量算出式に上記変形量を代入し、最小二乗法を用いて当該変形量算出式の少なくとも収差モデル量を表す係数を決定するステップと、上記決定された収差モデル量の係数を用いて歪曲収差の補正を行うステップとが具備されている。

β_ｘ＝ａ_１ｘ＋ａ_２ｙ＋ａ_３＋α_ｘ ・・・（Ａ）
β_ｙ＝ａ_４ｘ＋ａ_５ｙ＋ａ_６＋α_ｙ ・・・（Ｂ）
このように、格子パターンの撮影画像での位相値を求める際に、収差モデル量に格子パターンの回転量および並行移動量を表す平面の式を加えた変形量算出式を用いるとともに、この変形量算出式を表す平面の方程式の各係数を、実際に得られた変形量を用いて最小二乗法により決定するようにしたので、格子パターンを用いて歪曲収差を補正する際に、格子パターンがずれている場合でも、そのずれ量を考慮した補正を行うことができるので、精度良く歪曲収差の補正を行うことができる。

ところで、上記実施の形態においては、格子パターンの撮影画像にフーリエ変換を施して周波数スペクトルを得るとともに、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去し、そしてこの１次調和波を、収差が無い格子状パターンの周波数スペクトル分だけ原点側にシフトさせた後、フーリエ逆変換を施して位相値を求めるように説明したが、変形量算出式｛（１１）式および（１２）式｝そのものに、格子パターンの回転量および並行移動量が含まれているため、１次調和波を原点側にシフトさせなくてもよい。

また、上記実施の形態においては、円周方向歪曲収差を考慮した場合について説明したが、半径方向歪曲収差だけを考慮することもできる。この場合、（４）式の替わりに、（３）式が用いられる。

さらに、上記実施の形態においては、格子状パターンとして、互いに平行に設けられた格子線が、互いに直交するもの以外に、互いに９０度以外の角度で傾斜するものであってもよい。

本発明の実施の形態に係る歪曲収差補正方法を説明するための歪曲収差を示す図である。 同歪曲収差補正方法における基準となる格子パターンの図である。 同変形後の格子パターンの図である。 同変形後の格子パターンにおける周波数スペクトルを示す。 同変形後の格子パターンにおける位相分布を示し、（ａ）はｘ方向位相分布、（ｂ）はｙ方向位相分布である。 同変形後の格子パターンにおける位相接続後の位相分布を示し、（ａ）はｘ方向位相分布、（ｂ）はｙ方向位相分布である。 同歪曲収差補正方法における格子パターンを示す図で、（ａ）は基準となるものを示し、（ｂ）は回転したものを示す。 同歪曲収差補正方法における格子パターンで、ｘ方向に変形したものを示す。 板体を剛体並進変位させた場合の変位分布を示す図で、（ａ）はｘ方向変位、（ｂ）はｙ方向変位を示す。 同歪曲収差補正方法による補正後の変位分布と補正前の変位分布とを比較する図で、（ａ）はｘ方向変位、（ｂ）はｙ方向変位を示す。 平板の単軸引張における収差補正前の変位分布を示し、（ａ）はｘ方向変位、（ｂ）はｙ方向変位を示す。 同平板の単軸引張における収差補正後の変位分布を示し、（ａ）はｘ方向変位、（ｂ）はｙ方向変位を示す。

Claims (2)

1. 二次元の格子状パターンをレンズを介して撮影して得られた輝度値を示す撮影画像にフーリエ変換を施して周波数スペクトルを得るステップと、
   フィルタリング処理により、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去するステップと、
   上記１次調和波にフーリエ逆変換を施して位相値を求めるステップと、
   上記ステップにて所定範囲に亘って求められた位相値の位相接続を行うステップと、
   上記位相接続により得られた位相値（δ _ｘ ，δ _ｙ ）を下記（Ａ）式および（Ｂ）式に代入して格子状パターンの変形量（β_ｘ，β_ｙ）を求めるステップと、
   上記格子状パターンの回転移動および並行移動を考慮した平面の式で表されるずれ量算出式に当該格子状パターンに生じる下記（Ｃ）式および（Ｄ）式に基づく収差モデル量（α_ｘ，α_ｙ）を加えてなる下記（Ｅ）式および（Ｆ）式で示される変形量算出式に上記変形量を代入し、最小二乗法を用いて当該変形量算出式の少なくとも収差モデル量を表す係数を決定するステップと、
   上記決定された収差モデル量の係数を用いて歪曲収差を補正するステップと
   を具備したことを特徴とする撮影画像における歪曲収差補正方法。
   β _ｘ ＝−δ _ｘ ／２πω ・・・（Ａ）
   β _ｙ ＝−δ _ｙ ／２πω ・・・（Ｂ）
   但し、ωは変形前の格子の周波数である。
   α _ｘ ＝（ｇ _１ ＋ｇ _３ ）ｘ ^２ ＋ｇ _４ ｘｙ＋ｇ _１ ｙ ^２ ＋ｋ _１ ｘ（ｘ ^２ ＋ｙ ^２ ）・・（Ｃ）
   α _ｙ ＝ｇ _２ ｘ ^２ ＋ｇ _３ ｘｙ＋（ｇ _２ ＋ｇ _４ ）ｙ ^２ ＋＋ｋ _１ ｙ（ｘ ^２ ＋ｙ ^２ ）・・（Ｄ）
   β _ｘ ＝ａ _１ ｘ＋ａ _２ ｙ＋ａ _３ ＋α _ｘ ・・・（Ｅ）
   β _ｙ ＝ａ _４ ｘ＋ａ _５ ｙ＋ａ _６ ＋α _ｙ ・・・（Ｆ）
2. 二次元の格子状パターンをレンズを介して撮影して得られた輝度値を示す撮影画像にフーリエ変換を施して周波数スペクトルを得るステップと、
   フィルタリング処理により、この周波数スペクトルの内、１次調和波だけを残してノイズ分を除去するステップと、
   上記１次調和波を、収差が無い格子状パターンの周波数スペクトル分だけ原点側にシフトさせた後、フーリエ逆変換を施して位相値を求めるステップと、
   上記ステップにて所定範囲に亘って求められた位相値の位相接続を行うステップと、
   上記位相接続により得られた位相値（δ _ｘ ，δ _ｙ ）を下記（Ａ）式および（Ｂ）式に代入して格子状パターンの変形量（β_ｘ，β_ｙ）を求めるステップと、
   上記格子状パターンの回転移動および並行移動を考慮した平面の式で表されるずれ量算出式に当該格子状パターンに生じる下記（Ｃ）式および（Ｄ）式に基づく収差モデル量（α_ｘ，α_ｙ）を加えてなる下記（Ｅ）式および（Ｆ）式で示される変形量算出式に上記変形量を代入し、最小二乗法を用いて当該変形量算出式の少なくとも収差モデル量を表す係数を決定するステップと、
   上記決定された収差モデル量の係数を用いて歪曲収差を補正するステップと
   を具備したことを特徴とする撮影画像における歪曲収差補正方法。
   β _ｘ ＝−δ _ｘ ／２πω ・・・（Ａ）
   β _ｙ ＝−δ _ｙ ／２πω ・・・（Ｂ）
   但し、ωは変形前の格子の周波数である。
   α _ｘ ＝（ｇ _１ ＋ｇ _３ ）ｘ ^２ ＋ｇ _４ ｘｙ＋ｇ _１ ｙ ^２ ＋ｋ _１ ｘ（ｘ ^２ ＋ｙ ^２ ）・・（Ｃ）
   α _ｙ ＝ｇ _２ ｘ ^２ ＋ｇ _３ ｘｙ＋（ｇ _２ ＋ｇ _４ ）ｙ ^２ ＋＋ｋ _１ ｙ（ｘ ^２ ＋ｙ ^２ ）・・（Ｄ）
   β _ｘ ＝ａ _１ ｘ＋ａ _２ ｙ＋ａ _３ ＋α _ｘ ・・・（Ｅ）
   β _ｙ ＝ａ _４ ｘ＋ａ _５ ｙ＋ａ _６ ＋α _ｙ ・・・（Ｆ）

Images