(57)【要約】【課題】高精度にナイフエッジの傾斜角を推定する傾斜角推定装置を提供する。【解決手段】傾斜角推定装置１は、ナイフエッジ画像を２次元離散フーリエ変換するフーリエ変換手段１２１と、空間周波数領域のデータを用いて一方の空間周波数成分ｆｘ毎に、最大振幅を与える他方の空間周波数成分ｆｙを検出する空間周波数成分検出手段１２２と、検出された（ｆｘ，Ｆｙ（ｆｘ））の組について、直線の方程式である式（１）にフィットする角度θを算出して、ナイフエッジ画像の延伸方向が、他方の空間周波数の方向と成す角度であると推定する角度推定手段１３と、を備える。ここで、式（１）は、０＝ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｘ）・ｓｉｎθであり、Ｆｙ（ｆｘ）は前記一方の空間周波数成分ｆｘにおいて最大振幅を与える前記他方の空間周波数成分ｆｙを示す。【選択図】図１
2014203162.tif 000002
【特許請求の範囲】
【請求項１】
ナイフエッジ画像について、前記ナイフエッジの傾斜角を推定する傾斜角推定装置であって、
前記ナイフエッジ画像を２次元離散フーリエ変換することにより空間周波数領域のデータに変換するフーリエ変換手段と、
前記空間周波数領域のデータについて、一方の空間周波数成分ｆｘ毎に、最大振幅を与える他方の空間周波数成分ｆｙを検出する空間周波数成分検出手段と、
前記空間周波数成分検出手段により検出された（ｆｘ，ｆｙ）の組について、直線の方程式である式（１）にフィットする角度θを算出して、前記ナイフエッジの延伸方向が前記他方の空間周波数の方向と成す角度であると推定する角度推定手段と、
を備えることを特徴とする傾斜角推定装置。
０＝ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｘ）・ｓｉｎθ  …式（１）
但し、式（１）において、Ｆｙ（ｆｘ）は前記一方の空間周波数成分ｆｘにおいて最大振幅を与える前記他方の空間周波数成分ｆｙを示し、θは０度以上１８０度未満の、式（１）で表わされる直線が前記一方の空間周波数軸と成す角度である。
【請求項２】
前記ナイフエッジ画像の点対称画像を生成し、前記ナイフエッジ画像の一の辺と当該一の辺に対応する前記点対称画像の辺とが接するように、前記ナイフエッジ画像と前記点対称画像とを合成した合成ナイフエッジ画像を生成する画像調製手段を更に有し、
前記フーリエ変換手段は、前記ナイフエッジ画像に代えて、前記合成ナイフエッジ画像を２次元離散フーリエ変換することにより空間周波数領域のデータに変換することを特徴とする請求項１に記載の傾斜角推定装置。
【請求項３】
前記角度推定手段は、式（２）で示される評価式の値Ｊが最小となるθを、前記ナイフエッジの傾斜角であると推定することを特徴とする請求項１又は請求項２に記載の傾斜角推定装置。
Ｊ＝Σ｜ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｙ）・ｓｉｎθ｜ｐ  …式（２）
但し、式（２）において、Σは複数の空間周波数成分ｆｘについて右式の和を算出することを示し、ｐは正の実数を示す。
【請求項４】
前記式（２）において、ｐが２未満であることを特徴とする請求項３に記載の傾斜角推定装置。
【請求項５】
撮像手段により撮像されたナイフエッジチャートの画像から、ナイフエッジ画像を抽出する画像抽出手段と、
前記ナイフエッジ画像について、ナイフエッジの傾斜角を推定する傾斜角推定手段と、
前記推定された傾斜角の方向に各画素を投影し、当該ナイフエッジの傾斜に直交する次元に関する画像値の分布であるエッジ広がり関数を算出するエッジ広がり関数算出手段と、
前記エッジ広がり関数を微分し、更に離散フーリエ変換することにより、ＭＴＦ（Modulation Transfer Function）を算出するＭＴＦ算出手段と、を備え、
前記傾斜角推定手段として、請求項１乃至請求項４の何れか一項に記載の傾斜角推定装置を用いることを特徴とするＭＴＦ測定装置。
【請求項６】
コンピュータを、請求項１乃至請求項４の何れか一項に記載の傾斜角推定装置として機能させるための傾斜角推定プログラム。
【請求項７】
コンピュータを、請求項５に記載のＭＴＦ測定装置として機能させるためのＭＴＦ測定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、図形の傾斜角を推定する傾斜角推定装置、傾斜角推定装置により推定された傾斜角を利用して撮像手段についてのＭＴＦを測定するＭＴＦ測定装置及びこれらの装置を実現するプログラムに関する。
【背景技術】
【０００２】
スキャナやデジタルカメラ、撮像レンズなどの撮像系の品質評価を行うための代表的な指標として、空間周波数に対する振幅伝達特性を表わすＭＴＦ（Modulation Transfer Function）が用いられてきた。ＭＴＦ測定方法として、水平又は垂直から傾斜した辺を有する図形が含まれるＭＴＦ測定用チャートを撮像したナイフエッジ画像を用いる方法であるSlanted-edge法（以下、適宜に「エッジ法」という）がある（例えば、非特許文献１、特許文献１参照）。
このエッジ法では、ＭＴＦ算出の過程でナイフエッジ画像の傾斜角が用いられる。
【０００３】
また、一般の画像について、直線等の図形を自動的に検出する方法として、ハフ変換が知られている（例えば、非特許文献２参照）。具体的には、撮像した画像に対して図形の輪郭抽出処理を行い、抽出した輪郭を構成する画素についてハフ変換を適用して直線を検出するものである。ここで、ハフ変換による直線の検出結果として、その直線の傾斜角（勾配）を得ることができる。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特許第４０９２８５３号公報
【非特許文献】
【０００５】
【非特許文献１】ISO 12233:2000,”Photography - Electronic still-picture cameras - Resolution measurements.”
【非特許文献２】Duda,R.O. and Hart,P.E,”Use of the Hough Transformation to Detect Lines and Curves in Pictures”,Comm.ACM,Vol.15,pp.11-15,Jan,1972
【発明の概要】
【発明が解決しようとする課題】
【０００６】
従来、ハフ変換を用いた手法や、線形回帰法などによりフィッティングする手法による傾斜角の推定においては、前記したように処理の初期段階として図形の輪郭抽出処理を行っている。この処理は画像中のノイズ要因に強く影響されるため、従来の手法ではナイフエッジ画像の傾斜角を高精度に推定できないことがあった。
【０００７】
また、エッジ法によるＭＴＦ測定では、その測定精度にエッジ画像の傾斜角の推定精度が大きく影響する。ここで、図１０及び図１１を参照して、その例について説明する。
なお、図１０は、従来のエッジ法によるＭＴＦ測定において、ナイフエッジ傾斜角の推定値の偏差がＭＴＦ測定値に与える影響をシミュレーションした結果を示す図であり、（ａ）は真の傾斜角が１度の場合、（ｂ）は真の傾斜角が５度の場合を示す。なお、図１０において、横軸は空間周波数、縦軸はＭＴＦ値を示し、空間周波数＝０．５がナイキスト周波数である。
【０００８】
図１０(ａ)及び図１０（ｂ）において、太い実線で示したＭＴＦ特性が、傾斜角の偏差（すなわち、誤差）が「０」の場合、すなわち真のＭＴＦ特性を示すものである。これに対して、誤差として±０．１度（破線及び一点鎖線）及び±０．５度（細い破線及び二点鎖線）だけ偏差した場合のＭＴＦ特性は、真のＭＴＦ特性から大きく外れており、僅か±０．１度の誤差であっても、ＭＴＦ測定値に大きく影響することが分かる。
【０００９】
また、図１１は、従来のエッジ法によるＭＴＦ測定において、ナイフエッジ傾斜角の推定誤差がナイキスト周波数におけるＭＴＦ測定値に与える影響をシミュレーションした結果を示す図である。図１１において、横軸が推定誤差である偏差を示し、縦軸がＭＴＦ値を示す。また、真の傾斜角が５度の場合を実線で、１度の場合を破線で示している。
図１１に示すように、偏差が「０」の場合はＭＴＦ値が約「０．９」であり、偏差の絶対値の増加とともに、急激にＭＴＦ値が低下することがわかる。なお、偏差「０」のときのＭＴＦ値が「１」よりも小さいのは、アパーチャ開口が有限の値を有しているとしているためである。
このように、エッジ法によりＭＴＦを精度よく測定するには、高精度にナイフエッジの傾斜角を推定することが必要である。
【００１０】
本発明は、かかる問題に鑑みて創案されたものであり、高精度にナイフエッジの傾斜角を推定する傾斜角推定装置及びそのプログラムを提供することを課題とする。
また、本発明は、前記した傾斜角推定装置によって推定されたナイフエッジの傾斜角を用いて、Slanted-edge法によりＭＴＦ測定を行うＭＴＦ測定装置及びそのプログラムを提供することを他の課題とする。
【課題を解決するための手段】
【００１１】
前記した課題を解決するために、本発明に係る傾斜角推定装置は、ナイフエッジ画像について、前記ナイフエッジの傾斜角を推定する傾斜角推定装置であって、フーリエ変換手段と、空間周波数成分検出手段と、角度推定手段とを備えることとした。
【００１２】
かかる構成によれば、傾斜角推定装置は、フーリエ変換手段によって、前記ナイフエッジ画像を２次元離散フーリエ変換することにより空間周波数領域のデータに変換する。また、傾斜角推定装置は、空間周波数成分検出手段によって、前記空間周波数領域のデータについて、一方の空間周波数成分ｆｘ毎に、最大振幅を与える他方の空間周波数成分ｆｙを検出する。そして、傾斜角推定装置は、角度推定手段によって、前記空間周波数成分検出手段により検出された（ｆｘ，ｆｙ）の組について、直線の方程式である式（１）にフィットする角度θを算出して、前記ナイフエッジの延伸方向が前記他方の空間周波数の方向と成す角度であると推定する。すなわち、前記ナイフエッジの傾斜角が、前記他方の空間周波数の方向であるｙ軸方向と成す角度として推定される。
０＝ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｘ）・ｓｉｎθ  …式（１）
但し、式（１）において、Ｆｙ（ｆｘ）は前記一方の空間周波数成分ｆｘにおいて最大振幅を与える前記他方の空間周波数成分ｆｙを示し、θは０度以上１８０度未満の、式（１）で表わされる直線が前記一方の空間周波数軸と成す角度である。
【００１３】
また、本発明に係る傾斜角推定装置は、更に画像調製手段を有するようにしてもよい。
かかる構成によれば、傾斜角推定装置は、画像調製手段によって、前記ナイフエッジ画像の点対称画像を生成し、前記ナイフエッジ画像の一の辺と当該一の辺に対応する前記点対称画像の辺とが接するように、前記ナイフエッジ画像と前記点対称画像とを合成した合成ナイフエッジ画像を生成する。そして、傾斜角推定装置は、前記フーリエ変換手段によって、前記ナイフエッジ画像に代えて、前記合成ナイフエッジ画像を２次元離散フーリエ変換することにより空間周波数領域のデータに変換する。
これによって、２次元離散フーリエ変換した際に、ナイフエッジ画像の両端におけるエッジ成分の混入によるアーティファクトを避けることができる。
【００１４】
また、本発明に係る傾斜角推定装置は、前記角度推定手段は、式（２）で示される評価式の値Ｊが最小となるθを、前記ナイフエッジの傾斜角であると推定することが好ましい。
Ｊ＝Σ｜ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｙ）・ｓｉｎθ｜ｐ  …式（２）
但し、式（２）において、Σは複数の空間周波数成分ｆｘについて右式の和を算出することを示し、ｐは正の実数を示す。
かかる構成によれば、傾斜角推定装置は、角度推定手段によって、精度よく傾斜角を推定することができる。
【００１５】
また、本発明に係る傾斜角推定装置は、前記式（２）において、ｐが２未満であることが好ましい。
かかる構成によれば、傾斜角推定装置は、角度推定手段によって、ノイズレベルの変化に対して頑健に傾斜角を推定することができる。
【００１６】
また、本発明に係るＭＴＦ測定装置は、画像抽出手段と、傾斜角推定手段と、エッジ拡がり関数算出手段と、ＭＴＦ算出手段とを備え、前記傾斜角推定手段として、本発明に係る傾斜角推定装置を用いて構成した。
【００１７】
かかる構成によれば、ＭＴＦ測定装置は、画像抽出手段によって、撮像手段により撮像されたナイフエッジチャートの画像から、ナイフエッジ画像を抽出する。また、ＭＴＦ測定装置は、傾斜角推定手段によって、前記ナイフエッジ画像について、ナイフエッジの傾斜角を推定する。また、ＭＴＦ測定装置は、エッジ拡がり関数算出手段によって、前記推定された傾斜角の方向に各画素を投影し、当該ナイフエッジの傾斜に直交する次元に関する画像値の分布であるエッジ広がり関数を算出する。そして、ＭＴＦ測定装置は、ＭＴＦ算出手段によって、前記エッジ広がり関数を微分し、更に離散フーリエ変換することにより、ＭＴＦを算出する。
これによって、ＭＴＦ測定装置は、傾斜角推定手段によって高精度に推定されたナイフエッジの傾斜角を用いて、ＭＴＦを高精度に算出することができる。
【００１８】
また、本発明に係る傾斜角推定装置及びＭＴＦ測定装置は、ＣＰＵ（Central Processing Unit）、記憶手段（例えば、メモリ、ハードディスク）等のハードウェア資源を備えるコンピュータを、前記した各手段として協調動作させるための、それぞれ傾斜角推定プログラム及びＭＴＦ測定プログラムによって実現することもできる。このプログラムは、通信回線を介して配布してもよく、光ディスクや磁気ディスク、フラッシュメモリ等の記録媒体に書き込んで配布してもよい。
【発明の効果】
【００１９】
本発明によれば、ナイフエッジ画像について空間周波数領域におけるフィッティングを行うことにより傾斜角を推定するため、ノイズが重畳された画像であっても頑健で高精度に傾斜角を推定することができる。
また、本発明によれば、頑健で高精度に推定されたナイフエッジの傾斜角を用いて、Slanted-edge法によるＭＴＦ測定を行うため、高精度にＭＴＦを算出することができる。
【図面の簡単な説明】
【００２０】
【図１】本発明の第１実施形態に係る傾斜角推定装置の構成を示すブロック図である。
【図２】本発明の第１実施形態に係る傾斜角推定装置において用いられる画像の例を示し、（ａ）は撮像画像からナイフエッジ画像を抽出し傾斜角推定用画像を調整する様子を示し、（ｂ）は調製したナイフエッジ画像を示し、（ｃ）は（ｂ）のナイフエッジ画像のパワースペクトルを示し、（ｄ）は検出点を示す。
【図３】本発明の第１実施形態に係る傾斜角推定装置の動作を示すフローチャートである。
【図４】ノイズを重畳した種々の傾斜角のナイフエッジ画像について、発明法及び従来法により傾斜角を推定した結果を示す図である。
【図５】ノイズレベルを変えて重畳したナイフエッジ画像について、発明法及び従来法により傾斜角を推定した結果を示す図であり、（ａ）は傾斜角が１度の場合、（ｂ）は傾斜角が５度の場合である。
【図６】本発明の第１実施形態に係る傾斜角推定装置において、傾斜角の推定値を評価するための評価式のｐ値と推定誤差との関係を示す図である。
【図７】本発明の第２実施形態に係るＭＴＦ測定装置の構成を示すブロック図である。
【図８】本発明の第２実施形態に係るＭＴＦ測定装置におけるエッジ拡がり関数の算出処理を説明するための図である。
【図９】本発明の第２実施形態に係るＭＴＦ測定装置の動作を示すフローチャートである。
【図１０】従来のSlanted-edge法によるＭＴＦ測定法において、ナイフエッジ傾斜角の推定誤差がＭＴＦ測定値に与える影響をシミュレーションした結果を示す図であり、（ａ）は真の傾斜角が１度の場合、（ｂ）は真の傾斜角が５度の場合を示す。
【図１１】従来のSlanted-edge法によるＭＴＦ測定法において、ナイフエッジ傾斜角の推定誤差がナイキスト周波数におけるＭＴＦ測定値に与える影響をシミュレーションした結果を示す図である。
【発明を実施するための形態】
【００２１】
以下、本発明の実施形態について、適宜に図面を参照して詳細に説明する。
＜第１実施形態＞
［傾斜角推定装置の構成］
まず、第１実施形態に係る傾斜角推定装置の構成について、図１及び図２を参照して説明する。
第１実施形態に係る傾斜角推定装置１は、外部の撮像手段などから入力したナイフエッジ画像を２次元離散フーリエ変換処理によって空間周波数領域のデータに変換し、空間周波数領域において頑健なフィッティングを行うことにより、当該ナイフエッジの傾斜角を推定するものである。空間周波数領域でのフィッティングについては、各部の構成とともに順次説明する。
図１に示すように、本実施形態に係る傾斜角推定装置１は、画像調製手段１１と、周波数スペクトル算出手段１２と、角度推定手段１３とを備えて構成されている。
【００２２】
画像調製手段１１は、ナイフエッジ画像Ｇ２を入力し、傾斜角推定用の画像Ｇ４を調製するものである。このために、画像調製手段１１は、画像変換手段１１１と、画像合成手段１１２と、画像記憶手段１１３とを備えている。
【００２３】
画像変換手段１１１は、図２（ａ）に示すように、外部のカメラなどの撮像手段で撮像された画像Ｇ１から抽出されたナイフエッジ画像Ｇ２を入力し、入力したナイフエッジ画像Ｇ２を、ナイフエッジの延伸方向に沿った一方の辺の中点（図２（ａ）及び図２（ｂ）においてナイフエッジ画像Ｇ２の右辺に×印で示した点）を点対称移動の中心点として、点対称移動変換した画像（点対称画像）Ｇ３を生成する。画像変換手段１１１は、生成した画像Ｇ３を画像合成手段１１２に出力する。
なお、本実施形態では、点対称移動の中心点をナイフエッジ画像Ｇ２の右辺の中点として画像Ｇ３を生成するようにしたが、これに限定されるものではなく、任意の点を点対称移動の中心点として画像Ｇ３を生成することができる。
【００２４】
画像合成手段１１２は、外部からナイフエッジ画像Ｇ２を、画像変換手段１１１から画像Ｇ３を、それぞれ入力し、画像Ｇ２の右辺と画像Ｇ３の左辺とが接するように画像Ｇ２と画像Ｇ３とを合成することで、図２（ａ）及び図２（ｂ）に×印で示した画像中心に対して点対称な画像（合成ナイフエッジ画像）Ｇ４を生成する。画像合成手段１１２は、生成した画像Ｇ４を画像記憶手段１１３に記憶する。
なお、本実施形態では、ナイフエッジ画像Ｇ２と、その点対称画像である画像Ｇ３とを、ナイフエッジ画像Ｇ２の右辺と対応する画像Ｇ３の辺とが接するように合成したが、これに限定されるものではなく、ナイフエッジ画像Ｇ２の他の辺と当該他の辺に対応する画像Ｇ３の辺とが接するように合成してもよい。
【００２５】
ここで、一方向に延伸するナイフエッジのみを有するナイフエッジ画像Ｇ２と、点対称移動変換した画像Ｇ３とを合成することにより、ナイフエッジ画像を２次元離散フーリエ変換した際に、画像の両端における垂直方向又は水平方向のエッジ成分の混入によるアーティファクトを避け、傾斜角の推定精度を高めることができる。
なお、図２（ａ）に示した撮像画像Ｇ１が、平行なナイフエッジが連続するバーコードチャートを撮像した画像である場合は、図２（ｂ）に示すような画像領域を抽出して用いることができる。このため、図２（ｂ）に示すような画像をナイフエッジ画像Ｇ２として入力する場合は、点対称移動した画像との合成処理は不要である。
【００２６】
画像記憶手段１１３は、画像合成手段１１２から入力した画像Ｇ４を記憶するものである。画像記憶手段１１３に記憶した画像Ｇ４は、周波数スペクトル算出手段１２によって参照される。
【００２７】
周波数スペクトル算出手段１２は、画像調製手段１１の画像記憶手段１１３から点対称にナイフエッジが合成された画像Ｇ４を入力し、空間周波数領域のデータに変換してナイフエッジ画像の構成に寄与する空間周波数成分を検出し、検出した空間周波数成分を角度推定手段１３に出力するものである。
このために、周波数スペクトル算出手段１２は、フーリエ変換手段１２１と、空間周波数成分検出手段１２２とを備えている。
【００２８】
フーリエ変換手段１２１は、画像記憶手段１１３から画像Ｇ４を入力し、入力した画像Ｇ４を２次元離散フーリエ変換処理により空間周波数領域のデータに変換する。図２（ｃ）において、Ｇ５は、この２次元離散フーリエ変換のパワースペクトルを、パワーの大きさを画像の明暗（明るい方がパワーが大きい）で示したものである。なお、図２（ｃ）において、直交する２本の空間周波数軸ｆｘ，ｆｙの交点（原点）は、直流成分を示している。
フーリエ変換手段１２１は、画像Ｇ４をフーリエ変換したデータを空間周波数成分検出手段１２２に出力する。
【００２９】
空間周波数成分検出手段１２２は、フーリエ変換手段１２１から画像Ｇ４をフーリエ変換により空間周波数領域のデータに変換されたナイフエッジ画像を入力し、一方の空間周波数成分ｆｘ毎に、振幅が最大となる他方の空間周波数成分ｆｙを検出する。
ここで、空間周波数成分ｆｘにおいて、最大振幅を与えるｆｙを示す関数として、Ｆｙ（ｆｘ）を定義する。そして、空間周波数成分ｆｘの成分数をＮとすると、Ｎ個の点（ｆｘ，Ｆｙ（ｆｘ））が検出される。空間周波数成分検出手段１２２は、検出したＮ個の点（ｆｘ、Ｆｙ（ｆｘ））を、角度推定手段１３に出力する。
なお、図２（ｄ）において、これらの検出点を「●」で示している。
【００３０】
図２（ｂ）に示すように、画像Ｇ４において、実空間領域における座標軸であるｙ軸とナイフエッジの延伸方向とが成す角度（ナイフエッジ傾斜角）をθとする。
このように、画像の一の辺の端部から対向する辺の端部まで延伸する１つのナイフエッジのみを含む画像を空間周波数領域のデータに変換すると、図２（ｃ）に示すように、ナイフエッジの延伸方向に垂直な方向についての空間周波数成分が大きな振幅を有することになる。従って、前記したＮ個の検出点（ｆｘ，Ｆｙ（ｆｘ））は、図２（ｄ）に示すように、空間周波数領域における座標軸であるｆｘ軸と成す角度がθの直線上にプロットされる。
なお、図２（ｄ）において、当該直線から外れた位置にプロットされている検出点は、ノイズなどの影響を受けている空間周波数成分である。
【００３１】
ここで、（ｆｘ，Ｆｙ（ｆｘ））がプロットされる直線は、式（１）で表される。
ρ＝ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｘ）・ｓｉｎθ  …式（１）
但し、式（１）において、ρは空間周波数領域における原点から式（１）で表される直線との距離を示し、θは０度以上１８０度未満の、式（１）で表わされる直線が前記一方の空間周波数軸と成す角度である。また、エッジが直線状に延伸するナイフエッジ画像の場合は、ρ＝０となり、式（１）は原点を通る直線となる。
【００３２】
なお、本実施形態では、ｘ軸方向の空間周波数成分ｆｘ毎に、最大振幅を与えるｙ軸方向の空間周波数成分ｆｙを検出するようにしたが、ｆｘとｆｙとを入れ替えて、ｙ軸方向の空間周波数成分ｆｙ毎に、最大振幅を与えるｆｘを検出するようにしてもよい。
また、水平方向に近い傾斜角を有するナイフエッジ画像Ｇ２を用いる場合は、画像Ｇ２を時計回り又は反時計回りに９０度回転した画像に変換し、更に点対称移動変換した画像Ｇ３と合成して用いるようにしてもよい。これによって、以降の処理を垂直方向に近い傾斜角を有するナイフエッジ画像Ｇ２を用いる場合と同じ手順で行うことができる。
【００３３】
角度推定手段１３は、周波数スペクトル算出手段１２の空間周波数成分検出手段１２２からＮ個の検出点（ｆｘ、Ｆｙ（ｆｘ））を入力し、入力したこれらの検出点の集合が最もフィットする直線を算出して、当該算出した直線の勾配をナイフエッジ傾斜角、すなわちナイフエッジの延伸方向がｙ軸方向と成す角度であると推定するものである。
フィットする直線の算出方法は特に限定されるものではないが、以下に具体例を挙げて説明する。本例では、図１に示すように、角度推定手段１３は、評価値算出手段１３１と傾斜角推定値更新手段１３２と傾斜角推定値記憶手段１３３とを備えている。
【００３４】
評価値算出手段１３１は、周波数スペクトル算出手段１２からＮ個の検出点（ｆｘ、Ｆｙ（ｆｘ））を入力し、入力したこれらの検出点について、式（２）で示される評価式の値である評価値Ｊを算出する。なお、評価値算出手段１３１は、傾斜角推定値記憶手段１３３に記憶されている傾斜角の推定値を式（２）におけるθの値として用いる。
Ｊ＝Σ｜ｆｘ・ｃｏｓθ＋Ｆｙ（ｆｘ）・ｓｉｎθ｜ｐ  …式（２）
但し、式（２）において、Σは複数の空間周波数成分ｆｘについての和を算出することを示し、ｐは正の実数を示す。
【００３５】
すなわち、評価値算出手段１３１は、検出点（ｆｘ，Ｆｙ（ｆｘ））の集合が、式（１）にフィットする度合いを評価するための評価値を式（２）によって算出するものである。
評価値算出手段１３１は、算出した評価値Ｊを傾斜角推定値更新手段１３２に出力する。
【００３６】
なお、式（２）を用いて評価値Ｊを算出する際に、すべての空間周波数成分ｆｘについての検出点（ｆｘ，Ｆｙ（ｆｘ））を用いるようにしてもよいが、所定の空間周波数以下の空間周波数成分ｆｘについての検出点を用いるようにしてもよい。ノイズの影響を受けやすい高周波成分を評価値Ｊの計算から除外することで、安定した精度で、より頑健に傾斜角を推定することができる。
【００３７】
式（２）において、ｐ値は正の実数であれば限定されるものではないが、例えば、ｐ＝２の場合は、式（２）は最小二乗距離を基準とするフィッティングを行うこととなる。また、ｐ＜２、例えば、ｐ＝１、ｐ＝１／２とすることにより、検出点（ｆｘ、Ｆｙ（ｆｘ））の中のはずれ値の寄与を抑えた頑健なフィッティングを行うことができる。また、ナイフエッジ画像Ｇ２に重畳されているノイズの性質を予め調べ、ノイズの性質に応じて、適宜に最適なｐ値を選択するようにしてもよい。
なお、ｐ値と傾斜角推定精度との関係の詳細については、後記する。
【００３８】
傾斜角推定値更新手段１３２は、評価値算出手段１３１から評価値Ｊを入力し、評価値Ｊがθに対して最小値となっているかどうかを判定する。評価値Ｊがθに対して最小値となっていない場合は、傾斜角推定値θを、評価値Ｊを最小化させるように逐次更新し、更新した傾斜角推定値θを傾斜角推定値記憶手段１３３に記憶する。
なお、傾斜角推定値の更新値の算出は、例えば、シンプレックス法を用いることができる。また、必要な精度よりも小さな差分量（例えば、０．００１度）を加減することで更新するようにしてもよい。また、評価値が最小値となっているかどうかの判定は、前回の評価値Ｊとの変化量が所定の閾値より小さい場合に最小値になっていると判定するようにすることができる。
また、傾斜角推定値更新手段１３２が、評価値Ｊが最小値となっていると判定した場合は、傾斜角推定装置１は、そのときに傾斜角推定値記憶手段１３３に記憶されている傾斜角推定値をナイフエッジ画像の傾斜角の推定値として出力する。
【００３９】
なお、本例では、評価値Ｊを算出すための式（２）において、右辺の絶対値記号内の式である式（１）は、理想的には「０」となる。従って、そのｐ乗の総和である評価値Ｊも、理想的には「０」となる。そこで、評価値Ｊがθに対して最小かどうかの判定方法として、評価値Ｊが十分に小さいかどうかで判定することもできる。十分に小さいかどうかは、前記したように前回の評価値Ｊとの変化量が所定の閾値よりも小さい場合に十分に小さくなっている（すなわち、十分に小さい）と判定することもできるが、変化量ではなく当該評価値Ｊが一定の値よりも小さいかどうかで判定することもできる。
【００４０】
傾斜角推定値記憶手段１３３は、傾斜角推定値更新手段１３２から、傾斜角推定値を入力し、記憶している傾斜角推定値と置き換えて記憶する。
また、傾斜角推定値記憶手段１３３に記憶されている傾斜角推定値は、評価値算出手段１３１によって参照される。なお、傾斜角推定値が傾斜角推定値更新手段１３２によって最初に更新される前の初期値として、キーボードなどの不図示の入力手段を介して、ナイフエッジ画像の傾斜角のおよその値を入力して設定するようにしてもよい。
【００４１】
なお、傾斜角推定装置１によって出力されるナイフエッジ画像の傾斜角の推定値は、例えば、後記するＭＴＦ測定装置１００（図６参照）において、エッジ法を用いたＭＴＦ測定の際に、ナイフエッジの傾斜角として用いられる。また、傾斜角の推定値は、ＭＴＦ測定での利用に限定されるものではなく、例えば画像認識処理において、種々の図形が含まれる画像の個々のエッジについて傾斜角を推定することにより、画像中の図形の形状を高精度に認識することもできる。
また、垂直や水平など、傾斜角が既知の図形などの指標を撮像した画像について、本発明の手法を適用して画像座標に対する当該指標の傾斜角を推定し、既知の傾斜角からのズレを算出することにより、撮像系の傾きを検出することができる。そして、検出した撮像系の傾きを用いて当該撮像系の傾きを補正することもできる。
【００４２】
［傾斜角推定装置の動作］
次に、図３を参照（適宜図１及び図２参照）して、本実施形態に係る傾斜角推定装置１の動作について説明する。
図３に示すように、傾斜角推定装置１は、まず、画像調製手段１１によって、外部から入力したナイフエッジ画像Ｇ２から、傾斜角推定用の画像Ｇ４を調製する（ステップＳ１１）。このステップＳ１１では、より詳細には、傾斜角推定装置１は、まず、画像変換手段１１１によって、ナイフエッジ画像Ｇ２を点対称移動変換した画像Ｇ３を生成する。次に、傾斜角推定装置１は、画像合成手段１１２によって、画像Ｇ２と画像Ｇ３とを合成し、画像の中心が点対称の中心となる画像Ｇ４を生成し、画像記憶手段１１３に記憶する。
【００４３】
次に、傾斜角推定装置１は、周波数スペクトル算出手段１２のフーリエ変換手段１２１によって、画像Ｇ４を２次元離散フーリエ変換することにより、空間周波数領域のデータに変換する（ステップＳ１２）。
そして、傾斜角推定装置１は、周波数スペクトル算出手段１２の空間周波数成分検出手段１２２によって、空間周波数領域のデータに変換されたナイフエッジ画像の一方の空間周波数成分ｆｘ毎に、Ｆｙ（ｆｘ）を、すなわち最大振幅を与える他方の空間周波数成分ｆｙを検出する（ステップＳ１３）。
【００４４】
次に、傾斜角推定装置１は、角度推定手段１３の評価値算出手段１３１によって、式（２）を用いて評価値Ｊを算出する（ステップＳ１４）。なお、このとき、式（２）における傾斜角推定値θとして、傾斜角推定値記憶手段１３３に記憶されている推定値を用いる。
【００４５】
次に、傾斜角推定装置１は、角度推定手段１３の傾斜角推定値更新手段１３２によって、ステップＳ１４で算出された評価値Ｊが十分に小さいかどうかを判定する（ステップＳ１５）。
十分に小さい場合は（ステップＳ１５でＹｅｓ）、傾斜角推定装置１は処理を終了する。このとき、傾斜角推定値記憶手段１３３に記憶されている推定値が、ナイフエッジ画像の傾斜角の最終的な推定値となる。
【００４６】
一方、十分に小さくなっていない場合は（ステップＳ１５でＮｏ）、傾斜角推定装置１は、傾斜角推定値更新手段１３２によって、評価値Ｊを最小化させるように傾斜角の推定値を変化させた値に更新し（ステップＳ１６）、更新した推定値を傾斜角推定値記憶手段１３３に記憶し（ステップＳ１７）、ステップＳ１４に戻る。
以上の手順により、傾斜角推定装置１は、ナイフエッジの傾斜角を推定する。
【００４７】
［シミュレーション］
次に、本発明による傾斜角推定方法（以下、適宜「発明法」という）と、実空間のナイフエッジ画像から輪郭抽出した画素に線形回帰法を用いて直線の勾配を算出する従来法とにより、ノイズが重畳されたナイフエッジ画像について、ナイフエッジの傾斜角の推定を行ったシミュレーション結果について説明する。
【００４８】
（第１実験）
まず、第１実験として、コンピュータグラフィックスにより真の傾斜角が１度から２０度までのノイズのないナイフエッジ画像を生成し、更に−１０ｄＢのガウシアンホワイトノイズを重畳した画像を用いて、傾斜角推定シミュレーションを行った。なお、発明法において、式（２）におけるｐ値は「０．５」とした。
【００４９】
図４に第１実験のシミュレーション結果を示す。図４において、横軸は真の傾斜角を示し、縦軸は傾斜角推定値を示す。また、発明法による結果を「○」印で示し、従来法による結果を「×」印で示し、更に推定の信頼度を誤差棒で示している。
図４に示すように、画像にノイズが重畳されている場合であっても、従来法に比べて、発明法では高い精度で傾斜角を推定できることが分かる。また、従来法による推定値が、真の傾斜角からの乖離が大きい５度以上の傾斜角についても、発明法では高い精度、かつ高い信頼度で推定できることが分かる。
【００５０】
（第２実験）
次に、第２実験として、コンピュータグラフィックスにより真の傾斜角が１度及び５度のノイズのないナイフエッジ画像を生成し、更に「−１３．０１０３ｄＢ」、「−１０ｄＢ」及び「−６．０２０６ｄＢ」のガウシアンホワイトノイズを重畳した画像を用いて、傾斜角推定シミュレーションを行い、ノイズレベルの傾斜角推定値の誤差への影響について評価した。なお、発明法において、式（２）におけるｐ値は「０．５」とした。
【００５１】
図５に第２実験のシミュレーション結果を示す。図５において、上段に傾斜角１度の場合を、下段に傾斜角５度の場合を示し、左側から順にノイズレベルが「−１３．０１０３ｄＢ」、「−１０ｄＢ」及び「−６．０２０６ｄＢ」の場合を示す。また、縦軸は、傾斜角推定値の誤差の絶対値を示す。
図５に示したように、発明法は従来法に比べて、ノイズレベルの変化に対しても頑健に傾斜角を推定できることが分かる。
【００５２】
（第３実験）
次に、第３実験として、コンピュータグラフィックスにより真の傾斜角が５度のノイズのないナイフエッジ画像を生成し、更に「−１３．０ｄＢ」、「−１０．０ｄＢ」及び「−６．９９ｄＢ」のガウシアンホワイトノイズを重畳した画像を用い、発明法におけるｐ値を０．１から５まで変化させて傾斜角推定シミュレーションを行い、ｐ値の傾斜角推定値の誤差への影響について評価した。
【００５３】
図６に第３実験のシミュレーション結果を示す。図６において、横軸はｐ値を示し、縦軸は推定誤差を示す。また、ノイズレベルが「−１３．０ｄＢ」の場合を「○」で、「−１０．０ｄＢ」の場合を「□」で、「−６．９９ｄＢ」の場合を「△」で示した。また、この順に、対応するノイズレベルの画像についての従来法による結果を「●」、「■」及び「▲」で図６の右端に示した。
図６に示すように、ｐ値が２以下の場合に、従来法よりも精度よく傾斜角を推定できることが分かる。特に、ｐ値を２未満、更に好ましくは１以下とすることにより、ノイズレベルの変化に対しても頑健に、かつ高精度に傾斜角を推定できることが分かる。
【００５４】
＜第２実施形態＞
次に、本発明の第２実施形態として、第１実施形態に係る傾斜角推定装置１を用いたＭＴＦ測定装置について説明する。
まず、図７を参照（適宜図１及び図２参照）して、第２実施形態に係るＭＴＦ測定装置の構成について説明する。
【００５５】
図７に示すように、第２実施形態に係るＭＴＦ測定装置１００は、画像記憶手段１０１と、ＲＯＩ取得手段１０２と、画像抽出手段１０３と、傾斜角推定手段１０４と、ＥＳＦ算出手段１０５と、ＭＴＦ算出手段１０６とを備えて構成されている。
【００５６】
また、ＭＴＦ測定装置１００は、ＭＴＦ測定用チャートを撮像するカメラやスキャナなどの撮像手段１１０からナイフエッジ画像が含まれる画像を入力するとともに、キーボード、マウスなどの入力装置との入力インターフェイス１２０を介して、入力画像からナイフエッジ画像を抽出する領域を示す情報を入力する。また、ＭＴＦ測定装置１００は、算出したＭＴＦのデータを、出力インターフェイス１３０を介して、外部の表示装置や記憶装置などに出力する。
【００５７】
なお、本実施形態に係るＭＴＦ測定装置１００は、第１実施形態に係る傾斜角推定装置１を傾斜角推定手段１０４として用いることでSlanted-edge法におけるナイフエッジ画像の傾斜角を高精度に推定し、この傾斜角の推定値を用いてSlanted-edge法によりＭＴＦを測定するものである。
【００５８】
画像記憶手段１０１は、カメラやスキャナなどの撮像手段１１０によって撮像されたＭＴＦ測定用のチャートの画像Ｇ１（図２（ａ）参照）を記憶するものである。画像記憶手段１０１に記憶された画像Ｇ１は、ＲＯＩ取得手段１０２及び画像抽出手段１０３によって参照される。
ＭＴＦ測定用チャートとしては、例えば、非特許文献１に記載された「Resolution test chart」を用いることができるが、これに限定されるものではない。例えば、図２（ａ）に示すように、白色の背景Ｇ１ｗに傾斜して配置された黒べたの矩形図形Ｇ１ｋを有するチャートを用いることができる。
【００５９】
なお、水平方向（ｘ軸方向）についてのＭＴＦを測定する場合は、ナイフエッジの延伸方向が垂直よりも僅かに傾斜（例えば、１度〜５度程度）するようにチャートを配置することが好ましい。また、垂直方向（ｙ軸方向）についてのＭＴＦを測定する場合は、ナイフエッジの延伸方向が水平よりも僅かに傾斜するようにチャートを配置することが好ましい。更にまた、他の方向のＭＴＦを測定する場合は、当該他の方向に垂直な方向よりも僅かにナイフエッジの延伸方向が傾斜するようにチャートを配置することが好ましい。
【００６０】
ＲＯＩ（Region of interest）取得手段１０２は、外部から入力した画像Ｇ１から、ＭＴＦ測定のために用いるナイフエッジ画像の領域を示す情報を取得する手段である。ここで、ＲＯＩ取得手段１０２は、オペレータが指定したナイフエッジ画像の領域を示す情報（例えば、画像Ｇ１における矩形領域の右上角及び左下角の画素の座標値など）を、入力インターフェイス１２０を介して取得する。
ＲＯＩ取得手段１０２は、取得したナイフエッジ画像の領域を示す情報を、画像抽出手段１０３に出力する。
【００６１】
画像領域の指定は、例えば、画像Ｇ１を表示する画像表示装置（不図示）とマウスなどの入力デバイス（不図示）を用いる方法を用いることができる。具体的には、図２（ａ）に示すように撮像された画像Ｇ１を画像表示装置に表示しつつ、マウスを用いて移動可能な矩形枠（図２（ａ）に破線で示した矩形枠を参照）を重畳表示し、所望の位置に矩形枠が表示されたときにマウスのボタンをクリックする。ＲＯＩ取得手段１０２は、このときに矩形枠が表示された領域に対応する画像Ｇ１の領域情報を、入力インターフェイス１２０を介して取得するようにすることができる。
【００６２】
画像抽出手段１０３は、画像記憶手段１０１から撮像手段１１０により撮像された画像Ｇ１を、ＲＯＩ取得手段１０２からナイフエッジ画像の領域を示す情報をそれぞれ入力し、画像４からナイフエッジ画像の領域を示す情報で指示された領域の画像をナイフエッジ画像Ｇ２として抽出する。画像抽出手段１０３は、抽出したナイフエッジ画像Ｇ２を、傾斜角推定手段１０４及びＥＳＦ算出手段１０５に出力する。
【００６３】
なお、ナイフエッジ画像Ｇ２として抽出する領域は、図２（ａ）及び図２（ｂ）に示すように、ナイフエッジの延伸方向を示す直線が、矩形領域の対向する２辺と交差するように選択する。図２（ａ）及び図２（ｂ）に示した例では、ナイフエッジの延伸方向を示す直線が、画像Ｇ２を示す矩形領域の上下の辺と交差している。すなわち、画像Ｇ２中には、１方向のナイフエッジのみが含まれる。１方向のナイフエッジのみが含まれる画像Ｇ２とすることにより、当該ナイフエッジの傾斜角を精度よく推定することができる。
なお、図２（ａ）及び図２（ｂ）に示した例は、水平方向（ｘ軸方向）のＭＴＦを測定するための画像である。
【００６４】
傾斜角推定手段１０４は、図１に示した第１実施形態に係る傾斜角推定装置１を用いるものである。そして、傾斜角推定手段１０４は、画像抽出手段１０３からナイフエッジ画像Ｇ２を入力し、当該ナイフエッジ画像Ｇ２の傾斜角を推定し、傾斜角の推定値をＥＳＦ算出手段１０５に出力する。
傾斜角推定手段１０４の構成は、図１に示した傾斜角推定装置１と同様であるから、詳細な説明は省略する。
【００６５】
ＥＳＦ（Edge Spread Function；エッジ拡がり関数）算出手段（エッジ拡がり関数算出手段）１０５は、画像抽出手段１０３からナイフエッジ画像Ｇ２を、傾斜角推定手段１０４からナイフエッジ画像Ｇ２におけるナイフエッジの傾斜角の推定値をそれぞれ入力する。ＥＳＦ算出手段１０５は、ナイフエッジ画像Ｇ２におけるナイフエッジの延伸方向から角度θだけ傾斜した方向に垂直な方向、すなわちｘ軸方向についてのプロファイルであるＥＳＦを、本来の画素ピッチよりも細かいサブ画素ピッチで算出するものである。
ＥＳＦ算出手段１０５は、算出したＥＳＦをＭＴＦ算出手段１０６に出力する。
【００６６】
ここで、サブ画素ピッチでＥＳＦを算出する方法について、図８を参照して説明する。
図８の下段の図において、エッジ位置より右側のハッチングを施した領域がテストチャートの黒べた図形を撮像した領域に該当し、ハッチングを施していない左側の領域が白い背景部を撮像した領域に該当する。また、ナイフエッジ延伸方向は、垂直方向を示すｙ軸に対して角度θだけ傾斜している。
【００６７】
ここで、ナイフエッジは直線であるから、ｙ軸方向に配列する各画素のナイフエッジ延伸方向に垂直な軸への投影像は、本来の画素ピッチの（ｓｉｎθ）倍のピッチのサブ画素で構成される１次元画像となる。
図８の上段の図は、縦軸を画素値として、サブ画素ピッチで構成された１次元画像であるＥＳＦを示している。
【００６８】
なお、ｘ軸はナイフエッジ延伸方向に垂直な軸から角度θだけ傾斜しているから、この１次元画像のｘ軸への投影像を考える。そうすると、図８の上段に示したＥＳＦは、本来の画素ピッチの（ｔａｎθ）倍のピッチのサブ画素でサンプリングされたｘ軸方向のＥＳＦであるとみることができる。
このように、傾斜したナイフエッジ画像を用いた超解像処理により、本来の画素ピッチにより規定されるサンプリング周波数を超える解像度のＥＳＦを算出することができる。
これによって、頑健かつ高精度にＭＴＦを測定することができる。
【００６９】
図７に戻って、ＭＴＦ測定装置１００の構成について説明を続ける。
ＭＴＦ算出手段１０６は、ＥＳＦ算出手段１０５からＥＳＦを入力し、ナイフエッジ図形の配置方向に応じた方向（本例の場合は、ｘ軸方向）についてのＭＴＦを算出する。
ＥＳＦを用いたＭＴＦの算出は、微分とフーリエ変換を用いた公知の手法により行うことができる。
【００７０】
具体的には、まず、１次元画像であるＥＳＦを微分することで、ＬＳＦ（Line Spread Function；線広がり関数）を算出する。次に、ＬＳＦを１次元離散フーリエ変換し、そのパワースペクトルを算出することにより、ｘ軸方向の各空間周波数成分についてのＭＴＦを算出することができる。
【００７１】
なお、本実施形態では、超解像処理を行うことでＥＳＦを算出し、超解像のＥＳＦを微分して超解像のＬＳＦを算出するようにしたがこれに限定されるものではない。非特許文献１に記載された手法のように、本来の画素ピッチで複数のＥＳＦを算出し、算出した複数のＥＳＦ毎に微分して複数のＬＳＦを算出し、超解像処理により、複数のＬＳＦから超解像のＬＳＦを算出するようにしてもよい。この場合も、ＬＳＦを用いたＭＴＦの算出は前記した手順と同様であるから説明は省略する。
【００７２】
［ＭＴＦ測定装置の動作］
次に、図９を参照（適宜図２及び図７参照）して、本実施形態に係るＭＴＦ測定装置１００の動作について説明する。
まず、ＭＴＦ測定装置１００は、画像記憶手段１０１によって、撮像手段１１０から入力したＭＴＦ測定用チャートを撮像した画像Ｇ１を記憶する（ステップＳ２１）。
次に、ＭＴＦ測定装置１００は、ＲＯＩ取得手段１０２によって、ステップＳ２１で記憶した画像Ｇ１から、ＭＴＦの測定に用いるナイフエッジ画像領域を示す情報を取得する（ステップＳ２２）。
次に、ＭＴＦ測定装置１００は、画像抽出手段１０３によって、ステップＳ２２で取得したナイフエッジ画像領域を示す情報を用いて、ステップＳ２１で記憶した画像Ｇ１から当該情報で示された領域の画像を、ナイフエッジ画像Ｇ２として抽出する（ステップＳ２３）。
【００７３】
次に、ＭＴＦ測定装置１００は、傾斜角推定手段１０４によって、ステップＳ２３で抽出したナイフエッジ画像Ｇ２を用いて、ナイフエッジの傾斜角を推定する（ステップＳ２４）。なお、ナイフエッジの傾斜角の推定処理は、図３に示した第１実施形態に係る傾斜角推定装置１による傾斜角の推定と同様であるから、詳細な説明は省略する。
次に、ＭＴＦ測定装置１００は、ＥＳＦ算出手段１０５によって、ステップＳ２４で推定したナイフエッジの傾斜角推定値を用いて、ステップＳ２３で抽出したナイフエッジ画像Ｇ２について超解像処理によりＥＳＦを算出する（ステップＳ２５）。
【００７４】
そして、ＭＴＦ測定装置１００は、ＭＴＦ算出手段１０６によって、ステップＳ２５で算出したＥＳＦを微分してＬＳＦを算出し、更にこのＬＳＦを１次元離散フーリエ変換してパワースペクトルを算出することにより、空間周波数成分毎のＭＴＦを算出する（ステップＳ２６）。
以上の手順により、ＭＴＦ測定装置１００は、ＭＴＦを高精度に測定することができる。
【００７５】
以上、本発明の実施形態に係る傾斜角推定装置及びＭＴＦ測定装置について、発明を実施するための形態により具体的に説明したが、本発明の趣旨はこれらの記載に限定されるものではなく、特許請求の範囲の記載に基づいて広く解釈されなければならない。また、これらの記載に基づいて種々変更、改変等したものも本発明の趣旨に含まれることはいうまでもない。
【符号の説明】
【００７６】
１      傾斜角推定装置
１１    画像調製手段
１１１  画像変換手段
１１２  画像合成手段
１１３  画像記憶手段
１２    周波数スペクトル算出手段
１２１  フーリエ変換手段
１２２  空間周波数成分検出手段
１３    角度推定手段
１３１  評価値算出手段
１３２  傾斜角推定値更新手段
１３３  傾斜角推定値記憶手段
１００  ＭＴＦ測定装置
１０１  画像記憶手段
１０２  ＲＯＩ取得手段
１０３  画像抽出手段
１０４  傾斜角推定手段
１０５  ＥＳＦ算出手段（エッジ拡がり関数算出手段）
１０６  ＭＴＦ算出手段
１１０  撮像手段
１２０  入力インターフェイス
１３０  出力インターフェイス
【図１】
2014203162.tif 000003
【図３】
2014203162.tif 000004
【図４】
2014203162.tif 000005
【図５】
2014203162.tif 000006
【図６】
2014203162.tif 000007
【図７】
2014203162.tif 000008
【図８】
2014203162.tif 000009
【図９】
2014203162.tif 000010
【図１０】
2014203162.tif 000011
【図１１】
2014203162.tif 000012
【図２】
2014203162.tif 000013