(57)【特許請求の範囲】
【請求項１】
ぼけ画像を取得する取得手段と、
前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成する生成手段と、を有し、
前記生成手段は、
前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、
前記生成手段は、前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する、ことを特徴とする画像処理装置。
【請求項２】
前記生成手段は、前記複数の演算式の種類または該複数の演算式のパラメータを変更して前記複数の補正信号に関する情報を生成することを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記生成手段は、前記信号に関する情報として、前記ぼけ推定領域における輝度分布または該輝度分布の微分値に関する情報を補正することを特徴とする請求項１または２に記載の画像処理装置。
【請求項４】
前記生成手段は、前記補正処理において前記複数の補正信号に関する情報を生成した場合、前記推定処理において、少なくとも２つ以上の補正信号に関する情報に関して共通のぼけを推定し、
前記推定処理において推定されるぼけの総数は、前記複数の補正信号に関する情報の総数よりも少ないことを特徴とする請求項１乃至３のいずれか１項に記載の画像処理装置。
【請求項５】
前記生成手段は、前記補正処理において前記複数の補正信号に関する情報を生成した場合、前記推定処理において、該複数の補正信号に関する情報の全てに関して共通の１つのぼけを推定することを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記生成手段は、同一のぼけを用いて前記信号に関する情報を補正し、前記複数の補正信号に関する情報を生成することを特徴とする請求項１乃至５のいずれか１項に記載の画像処理装置。
【請求項７】
前記生成手段は、前記反復演算処理の全ての補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成することを特徴とする請求項１乃至６のいずれか１項に記載の画像処理装置。
【請求項８】
前記生成手段は、前記補正処理において前記複数の補正信号に関する情報を生成した場合、前記推定処理において、該複数の補正信号に関する情報のそれぞれに重み付けを施して前記ぼけを推定することを特徴とする請求項１乃至７のいずれか１項に記載の画像処理装置。
【請求項９】
前記生成手段は、
前記反復演算処理の前記推定処理として、第１の推定処理および該第１の推定処理よりも後に第２の推定処理を行い、
前記第１の推定処理において推定された第１のぼけと、前記第２の推定処理において推定された第２のぼけとの差異を判定し、
前記差異が所定のしきい値よりも大きい場合、前記複数の補正信号に関する情報の数を増加して、再度、前記推定処理を行うことを特徴とする請求項１乃至８のいずれか１項に記載の画像処理装置。
【請求項１０】
前記生成手段は、前記ぼけ推定領域における信号変化の量に関する情報に基づいて、前記複数の補正信号に関する情報の数を決定することを特徴とする請求項１乃至９のいずれか１項に記載の画像処理装置。
【請求項１１】
前記生成手段は、
前記ぼけ推定領域の解像度を低下させて低解像ぼけ推定領域を生成し、
前記低解像ぼけ推定領域において、前記信号に関する情報に含まれるぼけを補正して前記補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理とを反復する反復演算処理を行って前記推定ぼけを生成し、
前記反復演算処理の過程において、前記低解像ぼけ推定領域の解像度を前記ぼけ推定領域の前記解像度に近づける、ことを特徴とする請求項１乃至１０のいずれか１項に記載の画像処理装置。
【請求項１２】
前記生成手段は、前記ぼけ推定領域のデノイジングを行い、
前記ぼけ推定領域に含まれるノイズ量を取得し、
前記ぼけ推定領域を周波数分解して周波数分解ぼけ推定領域を生成し、
前記ノイズ量に基づいて、前記周波数分解ぼけ推定領域に対するデノイジング処理を行い、
デノイジング処理後の周波数分解ぼけ推定領域を再合成する、ことを特徴とする請求項１乃至１１のいずれか１項に記載の画像処理装置。
【請求項１３】
前記推定ぼけを用いて前記ぼけ画像の少なくとも一部を補正する補正手段を更に有することを特徴とする請求項１乃至１２のいずれか１項に記載の画像処理装置。
【請求項１４】
光学系を介して形成された光学像を光電変換して画像信号を出力する撮像素子と、
ぼけ画像を取得する取得手段と、
前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成する生成手段と、を有し、
前記生成手段は、
前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、
前記生成手段は、前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する、ことを特徴とする撮像装置。
【請求項１５】
ぼけ画像を取得するステップと、
前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成するステップと、を有し、
前記推定ぼけを生成するステップにおいて、
前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、
前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する、ことを特徴とする画像処理方法。
【請求項１６】
ぼけ画像を取得するステップと、
前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成するステップと、をコンピュータに実行させるように構成された画像処理プログラムであって、
前記推定ぼけを生成するステップにおいて、
前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、
前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する、ことを特徴とする画像処理プログラム。
【請求項１７】
請求項１６に記載の画像処理プログラムを記憶していることを特徴とする記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、単一のぼけ画像から、ぼけ画像に作用している劣化成分であるぼけを推定する画像処理方法に関する。
【背景技術】
【０００２】
近年、表示装置の高精細化に伴い、撮影画像の高画質化が望まれている。撮影画像においては、撮影に用いられる光学系の収差や回折、または、撮影時の手ぶれなどの劣化要因により、被写体空間の情報が失われている。このため従来から、これらの要因に基づく撮影画像の劣化を補正し、より高画質な画像を得る手法が提案されている。そのような手法として、例えば、ＷｉｅｎｅｒフィルタやＲｉｃｈａｒｄｓｏｎ−Ｌｕｃｙ法などが挙げられる。しかし、これらの手法では、画像に作用している劣化成分（ぼけ）が既知でない場合、高い補正効果を得ることができない。
【０００３】
一方、従来から、手ぶれにより劣化した一枚の画像から手ぶれ成分を推定する手法が提案されている。特許文献１には、手ぶれにより劣化した一枚の画像から、既知である自然画像の強度勾配分布に関する統計情報を用いることにより、手ぶれ成分を推定する手法が開示されている。ここで自然画像とは、近代の人間が生活する上で自然と目にするような画像という意味である。このため自然画像は、樹木や動物などが写った画像に限定されるものではなく、人間、建築物、電子機器などが写った画像も含まれる。自然画像の性質として、信号の強度勾配に関するヒストグラム（強度勾配ヒストグラム）は、勾配の強さに応じて裾の重い分布（Ｈｅａｖｙ  Ｔａｉｌｅｄ  Ｄｉｓｔｒｉｂｕｔｉｏｎ）に従うことが知られている。特許文献１には、手ぶれが補正された画像の強度勾配ヒストグラムが、裾の重い分布に従うように制約をかけることにより、手ぶれ画像のみから手ぶれ補正画像を推定する手法が開示されている。そして、その手ぶれ補正画像と手ぶれ画像との比較結果に基づいて、ぶれ成分が推定される。このときのぶれ成分の推定精度は、手ぶれ画像内にエッジが多く含まれているほど上昇する。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】米国特許第７６１６８２６号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、ぼけ画像内のエッジ量は、写っている被写体の構造に依存する。このため、常に推定に充分なエッジ量が画像内に存在しているとは限らない。特に、画像内の領域に応じてぼけの形状が変化する（Ｓｈｉｆｔ−ｖａｒｉａｎｔと呼ばれる）場合、ぼけ画像内の特定の部分領域のみを用いてぼけの推定を行う必要があるため、エッジ量が不足しやすい。エッジやテクスチャなどの信号変化が少ない領域においては、ぼけの推定精度が極端に低下する。このような領域に対して特許文献１の推定手法を用いると、推定された手ぶれ成分には大きな誤差が生じる。これは、手ぶれ以外の回折、収差、デフォーカス、または、擾乱（これらの劣化要因を総称して、「ぼけ」と呼ぶ）においても同様である。
【０００６】
そこで本発明は、ぼけ画像における信号変化が少ない領域からでも、高精度にぼけを推定可能な画像処理装置、撮像装置、画像処理方法、画像処理プログラム、および、記憶媒体を提供する。
【課題を解決するための手段】
【０００７】
本発明の一側面としての画像処理装置は、ぼけ画像を取得する取得手段と、前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成する生成手段とを有し、前記生成手段は、前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、前記生成手段は、前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する。
【０００８】
本発明の他の側面としての撮像装置は、光学系を介して形成された光学像を光電変換して画像信号を出力する撮像素子と、ぼけ画像を取得する取得手段と、前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成する生成手段とを有し、前記生成手段は、前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、前記生成手段は、前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する。
【０００９】
本発明の他の側面としての画像処理方法は、ぼけ画像を取得するステップと、前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成するステップとを有し、前記推定ぼけを生成するステップにおいて、前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する。
【００１０】
本発明の他の側面としてのプログラムは、ぼけ画像を取得するステップと、前記ぼけ画像における少なくとも一部のぼけ推定領域を取得し、該ぼけ推定領域から推定ぼけを生成するステップと、をコンピュータに実行させるように構成された画像処理プログラムであって、前記推定ぼけを生成するステップにおいて、前記ぼけ推定領域における信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、該信号に関する情報と該補正信号に関する情報とに基づいてぼけを推定する推定処理と、を反復する反復演算処理を行って前記推定ぼけを生成し、前記反復演算処理における少なくとも一回の補正処理において、前記補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する。
【００１１】
本発明の他の側面としての記憶媒体は、前記プログラムを記憶している。
【００１２】
本発明の他の目的及び特徴は、以下の実施例において説明される。
【発明の効果】
【００１３】
本発明によれば、ぼけ画像における信号変化が少ない領域からでも、高精度にぼけを推定可能な画像処理装置、撮像装置、画像処理方法、画像処理プログラム、および、記憶媒体を提供することができる。
【図面の簡単な説明】
【００１４】
【図１】実施例１、３における推定ぼけの生成方法を示すフローチャートである。
【図２】実施例１における画像処理システムのブロック図である。
【図３】実施例１、２における画像処理システムの外観図である。
【図４】実施例１、３における画像処理方法を示すフローチャートである。
【図５】実施例１〜３におけるぼけ推定の一例を示す図である。
【図６】実施例１〜３におけるぼけ推定の一例を示す図である。
【図７】実施例２における画像処理システムのブロック図である。
【図８】実施例２における画像処理方法を示すフローチャートである。
【図９】実施例３における撮像システムのブロック図である。
【図１０】実施例３における撮像システムの外観図である。
【発明を実施するための形態】
【００１５】
以下、本発明の実施例について、図面を参照しながら詳細に説明する。各図において、同一の部材については同一の参照符号を付し、重複する説明は省略する。
【００１６】
本実施形態において、被写体空間の情報を劣化させる要因を「ぼけ」と総称するが、ぼけの要因としては、回折、収差、デフォーカス、ぶれ、および、擾乱が含まれる。本実施形態を説明する前に、これらのぼけの要因について詳述する。
【００１７】
回折とは、画像を撮影した撮像装置の光学系で発生する、回折に起因する劣化である。これは、光学系の開口径が有限であるために生じる。
【００１８】
収差とは、光学系で発生する理想波面からのずれに起因する劣化である。光学系の設計値により発生する収差を設計収差、光学系の製造誤差や環境変化により発生する収差を誤差収差と呼ぶ。環境変化とは、温度、湿度、気圧などの変化を示し、これらの変化に応じて光学系の性能も変化する。単に収差という場合、設計収差および誤差収差の両方が含まれる。
【００１９】
デフォーカスとは、光学系の焦点と被写体とが合致しないことに起因する劣化である。デフォーカスが画像の全体に及ぶものをピントずれ、デフォーカスが画像の一部のみであるものをデフォーカスぼけと呼ぶ。具体的は、デフォーカスぼけとは、画像内に異なる距離の主要被写体と背景とが存在し、主要被写体に合焦していた場合において、背景の情報を劣化させている成分をいう。この劣化の大きさは、背景がどれだけ主要被写体から奥行き方向に離れているかに依存して変化する。デフォーカスぼけは、例えば、画像内の主要被写体を際立たせる表現方法として用いられる。単にデフォーカスという場合、ピントずれとデフォーカスぼけとの両方が含まれる。
【００２０】
ぶれとは、撮影の際、露光中に被写体と撮像装置との間の相対関係（位置および角度）が変化することにより生じる劣化である。ぶれのうち、画像の全体を劣化させるものを手ぶれ、画像の一部が劣化するものを被写体ぶれと呼ぶ。単にぶれという場合、手ぶれおよび被写体ぶれの両方が含まれる。
【００２１】
擾乱とは、撮影中に被写体と撮像装置との間に存在する物質が揺らぐことにより生じる劣化である。例えば、大気の揺らぎや水中撮影における水の揺らぎが挙げられる。短秒露光中に擾乱が発生している場合、ぼけ画像では被写体空間の直線が揺らめいた曲線になる。この曲がりを補正するため、擾乱成分の異なる複数のフレーム（または連写）画像を合成することがある。しかし、このような合成処理では、エッジの曲がりを補正することができても、周波数成分の劣化は残存する。ここで取り扱うぼけ画像は、１回の露光中に擾乱に伴う周波数成分の劣化が起きた画像（長秒露光画像）だけでなく、前述のような複数フレーム（または連写）の画像を合成した画像を含む。
【００２２】
本実施形態において、点像強度分布はＰＳＦ（Ｐｏｉｎｔ  Ｓｐｒｅａｄ  Ｆｕｎｃｔｉｏｎ）と記す。また本実施形態において、１次元（モノクロ）画像に対する処理について説明するが、多次元（例えばＲＧＢ）画像にも同様に適用可能であり、ＲＧＢの各チャンネルで処理を行えばよい。チャンネル間における差がない、または無視できる程度のぼけでは、チャンネル数を減らして（ＲＧＢをモノクロ化するなど）、処理しても構わない。各チャンネルが異なる波長を取得している場合、チャンネル間における差がないぼけとしては、ぶれが挙げられる。
【００２３】
一方、収差、回折、デフォーカス、および、擾乱は、波長に依存してぼけが変化する。ここでデフォーカスは、ピント位置から大きく離れていても、軸上色収差の影響で波長によってぼけの広がりが変わる。しかし、これらの波長依存性を有するぼけに対しても、撮影に用いる撮像素子のサンプリング周波数に対して、チャンネル間の性能差が充分に小さければ、チャンネル間の差が無視できると考えてよい。ただし、これらチャンネル間における差が無視できるぼけを推定する場合、前述のチャンネル数を減らす処理を行うよりも、複数のチャンネルのままで１つのぼけを推定することが望ましい。ぼけの推定は、画像の信号勾配に関する情報を用いて行うため、その情報が多いほど推定精度が向上する。すなわち、チャンネル数を減らさずに複数の画像で推定したほうが、信号勾配の情報が増加する（ただし、各チャンネルの画像がそれぞれの比例倍に一致する場合、信号勾配の情報は増加しない）ため、ぼけをより高精度に推定することができる。また本実施形態では、チャンネルが異なる波長を撮影する例を挙げているが、それ以外のパラメータ（偏光など）でも同様に考えればよい。
【００２４】
ここで、具体的な実施例を説明する前に、本実施形態の簡単な概要について説明する。ぼけ画像からぼけを推定するためのぼけ推定領域を取得し、ぼけ推定領域内における信号勾配と、信号勾配のぼけが補正された補正信号勾配とを比較することにより、推定ぼけを生成する。このとき、ぼけ推定領域内における信号変化の量が少ないと、推定ぼけの推定精度が低下してしまう。そこで、異なる補正処理または補正パラメータを用いて、信号勾配から複数の補正信号勾配を生成する。これにより、ぼけ推定領域内における信号変化の量を増加させることができるため、推定ぼけの推定精度を向上させることができる。
【実施例１】
【００２５】
まず、図２および図３を参照して、本発明の実施例１における画像処理システムについて説明する。図２は、本実施例における画像処理システム１００のブロック図である。図３は、画像処理システム１００の外観図である。
【００２６】
画像処理システム１００は、撮像装置１０１、記録媒体１０２、表示装置１０３、出力装置１０４、および、画像処理装置１０５を備えて構成される。画像処理装置１０５は、通信部１０６、記憶部１０７、および、ぼけ補正部１０８（画像処理部）を有する。ぼけ補正部１０８は、取得部１０８１（取得手段）、生成部１０８２（生成手段）、および、補正部１０８３（補正手段）を有する。ぼけ補正部１０８の各部は、後述のように本実施例の画像処理方法を実行する。なお、撮像装置１０１の内部に画像処理装置１０５（ぼけ補正部１０８）を備えるように構成してもよい。
【００２７】
撮像装置１０１は、光学系１０１１（撮像光学系）および撮像素子１０１２を有する。光学系１０１１は、被写体空間からの光線を撮像素子１０１２上に結像させる。撮像素子１０１２は、複数の画素を有し、光学系１０１１を介して形成された光学像（被写体像）を光電変換して画像信号を出力する。撮像装置１０１は、撮像素子１０１２から出力された画像信号に基づいて撮影画像（ぼけ画像）を生成する。撮像装置１０１により得られたぼけ画像は、通信部１０６を介して画像処理装置１０５へ出力される。ぼけ画像に関しては、前述のような複数の種類のぼけの少なくとも一つが作用することにより、被写体空間の情報が劣化している。
【００２８】
記憶部１０７は、画像処理装置１０５に入力されたぼけ画像、および、ぼけ画像を撮影した際の撮影条件に関する情報を記憶する。ここで撮影条件とは、撮影時における撮像装置１０１の焦点距離、絞り、シャッタースピード、ＩＳＯ感度などである。ぼけ補正部１０８は、ぼけ画像に基づいて特定のぼけ成分を推定して補正し、ぼけ補正画像を生成する。ぼけ補正画像は、通信部１０６を介して、表示装置１０３、記録媒体１０２、および、出力装置１０４のいずれか一つまたは複数に出力される。表示装置１０３は、例えば液晶ディスプレイやプロジェクタなどである。ユーザは、表示装置１０３を介して、処理途中の画像を確認しながら作業を行うことができる。記録媒体１０２は、半導体メモリ、ハードディスク、ネットワーク上のサーバなどである。出力装置１０４は、プリンタなどである。画像処理装置１０５は、必要に応じて、現像処理やその他の画像処理を行う機能を有してもよい。
【００２９】
本実施例の画像処理方法を実現するため、ソフトウエア（画像処理プログラム）を、ネットワークまたはＣＤ−ＲＯＭなどの記憶媒体を介して、画像処理装置１０５に供給することもできる。この際、画像処理プログラムは、画像処理装置１０５のコンピュータ（又は、ＣＰＵやＭＰＵなど）によって読み出され、ぼけ補正部１０８の機能を実行する。
【００３０】
ぼけ補正部１０８は、図１および図４のフローチャートで示される画像処理を実行する。その画像処理に関する説明を行う前に、ぼけ推定領域（ぼけ画像の少なくとも一部）における信号勾配のぼけを補正する方法（補正処理）、および、ぼけを推定する方法（推定処理）に関して説明する。なお本実施例において、信号勾配とは、ぼけ推定領域における信号に関する情報であり、画像そのもの（画像の輝度分布）に関する情報またはそのｎ次微分（ｎは０を含む自然数）に関する情報（微分値に関する情報）を含む。
【００３１】
本実施例において、ぼけ（劣化成分）はＰＳＦの形で推定されるが、これに限定されるものではない。例えば、ＰＳＦをフーリエ変換したＯＴＦ（Ｏｐｔｉｃａｌ  Ｔｒａｎｓｆｅｒ  Ｆｕｎｃｔｉｏｎ）の形で推定しても構わない。まず、適当なＰＳＦの初期値（Ｇａｕｓｓ分布や１次元のラインなど）を与え、そのＰＳＦでぼけ推定領域における信号勾配のぼけを補正し、補正信号勾配を生成する。続いて、信号勾配と補正信号勾配とに基づいてＰＳＦを推定する。推定されたＰＳＦを用いて、再度ぼけ推定領域の信号勾配を補正して新たな補正信号勾配を生成し、ＰＳＦを推定する。これを反復すること（ループ処理）により、ぼけ画像（ぼけ推定領域）のみからぼけを推定することができる。
【００３２】
次に、信号勾配のぼけ補正方法（補正処理）について、具体的に説明する。ぼけ推定領域とぼけの関係は、以下の式（１）のように表される。
【００３３】
【数１】
0006541454.tif 000002
【００３４】
式（１）において、ｂ-ｉはｉ番目のループ処理におけるぼけ推定領域の信号分布、ｋｉはぼけ、ａｉはぼけｋｉによる劣化がない信号分布、ｎｉはノイズである。「＊」はコンボリューション演算を表す。
【００３５】
以下の式（２）のように表される最適化問題を解くことにより、ｉ番目のループ処理におけるａｉの推定値ｄ-ｉ（補正信号分布に相当）を推定する。
【００３６】
【数２】
0006541454.tif 000003
【００３７】
式（２）において、Ｌは損失関数、Φはｄｉに対する正則化項であり、それぞれの具体的な例については後述する。損失関数Ｌは、解をモデル（ここでは式（１）を指す）へフィッティングさせる効果を持つ。正則化項Φは、解を尤もらしい値へ収束させる効果を有する。正則化項には、事前知識と呼ばれる解（ａｉ）が持つべき性質を用いる。また、正則化項は損失関数のみを考えた際に起こる過剰なフィッティング（ノイズｎｉの影響をｄｉへ反映させてしまうこと）を防ぐ役割を有する。なお、ぼけ推定領域を複数のチャンネル（例えばＲＧＢ）から取得している場合、各チャンネルに対して式（２）を解く。
【００３８】
本実施例では、式（２）を解くことにより、補正信号勾配を求めるが、その他の手法を用いてもよい。例えば、Ｗｉｅｎｅｒフィルタなどの逆フィルタを用いた手法、または、Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙ法などの超解像手法を用いることができる。さらに、逆フィルタや超解像手法と合わせて、ショックフィルタなどの先鋭化フィルタを適用して補正信号勾配を生成してもよい。さらに、先鋭化フィルタを併用する際には、バイラテラルフィルタやガイデッドフィルタなどのエッジの解像感を保つ平滑化フィルタを用いて、ノイズやリンギングを抑制することが好ましい。複数の補正信号勾配を生成する際には、それぞれ異なる手法（例えば、式（２）の最適化と逆フィルタ）を用いて生成してもよい。
【００３９】
次に、式（２）中の損失関数Ｌおよび正則化項Φの具体例について説明する。損失関数Ｌとしては、以下の式（３）のように示される関数が考えられる。
【００４０】
【数３】
0006541454.tif 000004
【００４１】
式（３）において、以下の式（４）のように表される記号はｐ次平均ノルムを表し、ｐ＝２の場合にはユークリッドノルムを示す。
【００４２】
【数４】
0006541454.tif 000005
【００４３】
正則化項Φの例としては、以下の式（５）のように示される１次平均ノルムがある。
【００４４】
【数５】
0006541454.tif 000006
【００４５】
式（５）において、λは正則化項Φの重みを表すパラメータ、Ψは画像に対する基底変換を表す関数であり、その例としてウェーブレット変換や離散コサイン変換などがある。式（５）の正則化項Φは、画像がウェーブレット変換や離散コサイン変換などの基底変換を施されることで、信号成分がスパース（疎）になる、すなわち、より少ない数の信号で表すことができるという性質に基づく。これに関しては、例えば、「Ｒｉｃｈａｒｄ  Ｇ．Ｂａｒａｎｉｕｋ，‘‘Ｃｏｍｐｒｅｓｓｉｖｅ  Ｓｅｎｓｉｎｇ’’，ＩＥＥＥ  ＳＩＧＮＡＬ  ＰＲＯＣＥＳＳＩＮＧ  ＭＡＧＡＺＩＮＥ  ［１１８］  ＪＵＬＹ  ２００７」に説明されている。また、その他の正則化項の例として、Ｔｉｋｈｏｎｏｖ正則化項やＴＶ（Ｔｏｔａｌ  Ｖａｒｉａｔｉｏｎ）ノルム正則化項などを用いてもよい。
【００４６】
最適化問題である式（２）で表される推定式を解くには、繰り返し演算を用いた手法を利用し、例えばＴｉｋｈｏｎｏｖ正則化項を採用した場合は、共役勾配法などを使用すればよい。また、式（５）またはＴＶノルム正則化項を採用する場合、ＴｗＩＳＴ（Ｔｗｏ−ｓｔｅｐ  Ｉｔｅｒａｔｉｖｅ  Ｓｈｒｉｎｋａｇｅ／Ｔｈｒｅｓｈｏｌｄｉｎｇ）などを用いるとよい。ＴｗＩＳＴに関しては、「Ｊ．Ｍ．Ｂｉｏｕｃａｓ−Ｄｉａｓ，ｅｔ  ａｌ．，“Ａ  ｎｅｗ  ＴｗＩＳＴ：ｔｗｏ−ｓｔｅｐ  ｉｔｅｒａｔｉｖｅ  ｓｈｒｉｎｋａｇｅ／ｔｈｒｅｓｈｏｌｄｉｎｇ  ａｌｇｏｒｉｔｈｍｓ  ｆｏｒ  ｉｍａｇｅ  ｒｅｓｔｏｒａｔｉｏｎ”，ＩＥＥＥ  Ｔｒａｎｓ．ｏｎ  Ｉｍａｇｅ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．１６，Ｄｅｃ．２００７」に説明されている。
【００４７】
また、これらの繰り返し演算を行う際、繰り返しの度に正則化の重みなどのパラメータを更新してもよい。なお、式（１）、（２）は、画像（信号分布）に対して記載されているが、画像の微分に関しても同様に成り立つ。このため、画像に代えて、画像の微分（両者を含めて信号勾配と表現する）に対してぼけ補正を行ってもよい。
【００４８】
次に、ぼけの推定方法（推定処理）について具体的に説明する。ＰＳＦに関しても、式（２）と同様に、以下の式（６）を用いて推定することができる。
【００４９】
【数６】
0006541454.tif 000007
【００５０】
なお、ぼけ推定領域を複数のチャンネルから取得している（チャンネルによるぼけの変化がない、または無視できる）場合、式（６）は以下の式（６ａ）に示されるように変形される。
【００５１】
【数６ａ】
0006541454.tif 000008
【００５２】
式（６ａ）において、Ｈはぼけ推定領域に含まれるチャンネルの総数、ｄｉ，ｈはｈ番目のチャンネルにおけるｄｉ、ｖｈは重みを示す。信号勾配の補正の場合、式（２）を各チャンネルに対して解けばよいが、ぼけ推定では式（６ａ）のように全チャンネルで纏められた形になる。これは、推定する対象が、式（２）では各チャンネルで異なるのに対し、式（６ａ）では共通している（同一のＰＳＦである）ためである。
【００５３】
式（６）、（６ａ）の損失関数Ｌとして、以下の式（７）が考えられる。
【００５４】
【数７】
0006541454.tif 000009
【００５５】
式（７）において、∂ｊは微分演算子を表す。∂０（ｊ＝０）は恒等演算子であり、∂ｘ、∂ｙはそれぞれ、画像の横方向と縦方向の微分を表す。さらに高次の微分は、例えば∂ｘｘや∂ｘｙｙのように表される。なお、ぼけ推定領域の信号勾配とは、これら全てのｊ（ｊ＝０，ｘ，ｙ，ｘｘ，ｘｙ，ｙｙ，ｙｘ，ｘｘｘ，……）を含むが、本実施例ではｊ＝ｘ，ｙのみを考える。ｕｊは重みである。
【００５６】
式（６）の正則化項Φとしては、推定するＰＳＦの性質に沿ったものを用いることができる。例えば、ＰＳＦがデフォーカスのように潰れた形になる場合、以下の式（８）のように表されるＴＶノルム正則化を用いることが好ましい。
【００５７】
【数８】
0006541454.tif 000010
【００５８】
式（８）において、ζは正則化項の重みである。
【００５９】
また、ＰＳＦがぶれのようにライン状の形になる場合、以下の式（９）のように表される１次平均ノルム正則化を用いることが好ましい。
【００６０】
【数９】
0006541454.tif 000011
【００６１】
また、ＰＳＦの形に依存せず、ある程度の精度が確保され、かつ、高速に解ける正則化項として、Ｔｉｋｈｏｎｏｖ正則化を用いてもよい。
【００６２】
次に、本実施例の画像処理方法に関して説明する。本実施例では、補正信号勾配とぼけの推定を反復するループ処理の中で、少なくとも１回以上、複数の補正処理または補正パラメータを用いて複数の補正信号勾配を生成する。この効果について、以下に詳述する。まず、生成された補正信号勾配が単数の場合について説明し、その後、補正信号勾配を複数生成することによる効果について述べる。
【００６３】
式（１）から分かるように、ぼけ推定の際の推定対象であるぼけｋｉは、１次方程式で表現されている。この１次方程式は、ぼけ推定領域ｂｉの画素数だけ存在するため、ｋｉの推定は、ｂｉの画素数個ある連立方程式を解く問題に相当する。このため、未知数の数よりも独立な１次方程式の数が多ければ、解は一意的に決定されることになる。しかし、一般的に、未知数の数より独立な１次方程式の数が多くなることはない。これは、ぼけｋｉだけでなく、ノイズｎｉも未知であることによる（ぼけを求める際には、ａｉに推定したｄｉの値を代入するため、ａｉは既知となる）。すなわち、未知数の数は、ぼけｋｉの成分数と、ノイズｎｉの成分数（ｂｉの画素数と等しい）との和に相当し、常に１次方程式の数（ｂｉの画素数）より多い。このため、ぼけｋｉを求めるには、推定処理が必要となる。
【００６４】
ここで、ノイズｎｉの影響を抑え、高精度にぼけｋｉを求める方法について考える。まず、式（１）を、未知数ｋｉ（数はｋｉの成分数）の連立１次方程式（数はｂｉの画素数）それぞれに、ノイズが含まれていると考える。ただし、ｂ-ｉの画素数がｋｉの成分数より多いとする。仮に、ノイズが含まれていない場合、連立１次方程式の数が未知数の数よりも多いため、１次方程式の中には必ず従属関係となるものが存在する。すなわち、全く同じ方程式が複数存在することになる。しかし、実際にはノイズが存在するため、それらの方程式はノイズの分だけずれを生じている。一方、ノイズは、原信号が同じで異なるノイズの値を有する信号を平均化することにより、抑制可能である。このため、ノイズの出方が異なる従属な１次方程式が多ければ、平均化の効果によりノイズの影響を抑えたぼけ推定が可能になる。この平均化によるノイズ抑制効果が充分に作用している場合、式（６）、（６ａ）の正則化項Φを０にしてもよい。
【００６５】
また、ぼけｋｉを高精度に推定するには、更なる条件が存在する。それは、ノイズを除いた連立１次方程式中で、独立な方程式の数ができるだけ多いことである。ぼけｋｉの成分数（つまり未知数の数）よりも独立な方程式の数が少なすぎると、推定は充分な精度を得られない。連立１次方程式の係数は、補正信号勾配ｄｉであるから、独立な方程式の数が多くなるためには、補正信号勾配ｄｉが様々な組み合わせを持っている必要がある。換言すると、ぼけ推定領域内において、様々な方向または強度の信号変化（エッジ、テクスチャ、または、グラデーション）が存在している程、推定精度が向上する。例えば、ぼけ推定領域ｂｉおよび補正信号勾配ｄｉの信号値が一定（画像で言うと、完全な平坦部）であると仮定すると、式（１）の連立１次方程式は（ノイズを除いて）全て同じとなるため、ぼけｋｉを推定することができない。
【００６６】
以上の議論を纏めると、ぼけ推定領域内において様々な方向や強度を持つ信号変化が多く存在し、かつ、同じ信号変化でありながらノイズの出方が異なるものが複数存在すると、ぼけの推定精度が向上する。これを実現するため、本実施例では、補正信号勾配の生成の際において、１つのぼけ推定領域の信号勾配ｂｉから、複数の補正信号勾配ｄｉを生成する。複数の補正信号勾配（複数の補正信号に関する情報）の生成には、互いに異なる複数の演算式、すなわち、互いに異なる複数の補正処理またはパラメータ（補正パラメータ）を用いる。
【００６７】
補正処理とは、ぼけ補正の手法（式（２）による最適化、Ｗｉｎｅｒフィルタなどの逆フィルタ、Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙ法などの超解像手法）や、式（２）の最適化で用いる正則化項（式（５）の１次平均ノルム正則化やＴＶノルム正則化など）を指す。また、ぼけ補正後のフィルタリング（先鋭化フィルタおよび平滑化フィルタ）の有無も補正処理に含まれる。補正パラメータとは、式（２）の正則化項の重み（例えば、式（５）のλ）や、Ｒｉｃｈａｒｄｓｏｎ−Ｌｕｃｙ法などの超解像手法における高解像度化パラメータ、または、フィルタリングにおける先鋭化や平滑化の強さを表すパラメータである。
【００６８】
このような補正処理または補正パラメータを変更して複数の補正信号勾配を生成すると、それらは補正のかかり方が互いに異なるため、信号変化の振る舞い、および、ノイズの出方に差異が生じる。すなわち、複数の補正信号勾配を生成することにより、式（１）の連立方程式の数が増加する。例えば、２つの補正信号勾配を生成すれば、式（１）の方程式の数は２倍になる。これは、１つの補正信号勾配とぼけた信号勾配でぼけ推定領域の画素数個だけ連立方程式が立ち、もう一方の補正信号勾配とぼけた信号勾配で同じ数だけ連立方程式が立つためである。すなわち、１つのみの補正信号勾配を生成した場合では、連立１次方程式の数はぼけ推定領域の画素数個になるが、複数の補正信号勾配を生成した場合では、その生成した数だけ連立一次方程式の数が比例倍されることとなる。このため、複数の補正信号勾配を生成することにより、信号変化が少ないぼけ推定領域からでも、高精度にぼけを推定することができる。
【００６９】
複数の補正信号勾配が生成された場合、式（６）、（６ａ）はそれぞれ、以下の式（１０）、式（１０ａ）のようになる。
【００７０】
【数１０】
0006541454.tif 000012
【００７１】
【数１０ａ】
0006541454.tif 000013
【００７２】
式（１０）、（１０ａ）において、Ｇは共通のぼけの推定に用いられる補正信号勾配の総数、ｔｇは重みである。
【００７３】
次に、図４を参照して、ぼけ補正部１０８で行われる画像処理に関して説明する。図４は、本実施例の画像処理方法を示すフローチャートである。図４の各ステップは、ぼけ補正部１０８の取得部１０８１、生成部１０８２、および、補正部１０８３により実行される。
【００７４】
まずステップＳ１０１において、ぼけ補正部１０８の取得部１０８１は、ぼけ画像（撮影画像）を取得する。続いてステップＳ１０２において、ぼけ補正部１０８の生成部１０８２は、ぼけ画像内からぼけ推定領域を取得する。生成部１０８２は、ぼけ推定領域内において同じぼけが作用していると仮定して、ぼけを推定する。取得部１０８１は、ぼけ推定領域として、ユーザにより指定された領域を取得するか、または、ぼけ推定領域を自動的に取得する。ぼけ画像の全体においてぼけが一定の場合（Ｓｈｉｆｔ−ｉｎｖａｒｉａｎｔという）、ぼけ画像の全体をぼけ推定領域としてもよい。Ｓｈｉｆｔ−ｉｎｖａｒｉａｎｔの例として、露光中に撮像装置１０１の全体がシフトする手ぶれや、撮像装置１０１内の光学系１０１１（ヴィネッティングが充分に小さい場合）による回折などがある。また、ぼけ画像の位置に応じてぼけが変化する場合（Ｓｈｉｆｔ−ｖａｒｉａｎｔ）、ぼけ画像を複数の部分領域に分割し、分割した複数の部分領域のうちの１つの部分領域をぼけ推定領域とすることが好ましい。
【００７５】
続いてステップＳ１０３において、生成部１０８２は、ぼけ推定領域から推定ぼけを生成する。ステップＳ１０３の詳細は、図１のフローチャートに示されており、その説明に関しては後述する。続いてステップＳ１０４において、生成部１０８２は、ぼけ画像における所定の領域の全て（例えば、ぼけ画像の全体）に対して、推定ぼけの生成が完了したか否かを判定する。推定ぼけの生成が完了している場合、ステップＳ１０５へ進む。一方、推定ぼけの生成が完了していない場合、ステップＳ１０２へ戻る。そして生成部１０８２は、ぼけ画像における所定の領域内から、まだ推定ぼけが生成されていない部分領域を、新たなぼけ推定領域として取得する。
【００７６】
ステップＳ１０５において、生成部１０８２は、生成された推定ぼけを出力する。ぼけ画像における所定の領域内に複数のぼけ推定領域が含まれる場合（すなわち、Ｓｈｉｆｔ−ｖａｒｉａｎｔの場合）、複数の推定ぼけが存在する。このため生成部１０８２は、複数の推定ぼけを出力する。本実施例において、推定ぼけはＰＳＦに関するデータである。ただし、本実施例はこれに限定されるものではなく、ＯＴＦ、または、ＰＳＦもしくはＯＴＦを所定の基底でフィッティングした係数データ、または、ＰＳＦもしくはＯＴＦを画像データに変換した画像などの形で出力してもよい。
【００７７】
続いてステップＳ１０６において、ぼけ補正部１０８の補正部１０８３は、ステップＳ０１５にて生成部１０８２から出力された推定ぼけに基づいて、ぼけ画像に含まれるぼけ（劣化成分）を補正する。この補正には、例えば、Ｗｉｅｎｅｒフィルタなどの逆フィルタを用いた手法やＲｉｃｈａｒｄｓｏｎ−Ｌｕｃｙ法などの超解像手法を用いることができる。複数の推定ぼけが存在する場合、それぞれの推定ぼけが得られたぼけ推定領域に対して、対応する推定ぼけを用いて補正を行う。また、全ての推定ぼけが出力される前に、ステップＳ１０３にて生成された推定ぼけを用いて、逐次ぼけ推定領域のぼけを補正し、最後にそれらを合わせることにより、ぼけ補正画像を取得してもよい。
【００７８】
次に、図１を参照して、推定ぼけの生成処理（ステップＳ１０３）について詳述する。図１は、推定ぼけの生成方法を示すフローチャートである。図１の各ステップは、ぼけ補正部１０８の生成部１０８２により実行される。
【００７９】
まず、ステップＳ２０１において、生成部１０８２は、ぼけ推定領域のデノイジングを行う。ぼけ推定領域のデノイジングは、ぼけ推定領域にノイズが存在していることによる、ぼけの推定精度の劣化を低減するために行われる。なお、ステップＳ２０１に代えて、図４のステップＳ１０２の前に、ぼけ画像の全体をデノイジングするステップを挿入してもよい。デノイジング手法としては、バイラテラルフィルタやＮＬＭ（Ｎｏｎ  Ｌｏｃａｌ  Ｍｅａｎｓ）フィルタなどを用いる手法がある。
【００８０】
好ましくは、生成部１０８２は、以下の手法でぼけ推定領域をデノイジングする。まず、生成部１０８２は、ぼけ推定領域を周波数分解して、周波数分解ぼけ推定領域を生成する。そして生成部１０８２は、ぼけ推定領域のノイズ量に基づいて、周波数分解ぼけ推定領域をデノイジングする。次に、生成部１０８２は、周波数分解ぼけ推定領域を再合成することにより、ノイズ低減されたぼけ推定領域を取得する。一般に、デノイジング処理は、画像からノイズを低減するとともに、画像をぼかすという問題がある。ぼけ推定領域がデノイジングによってぼけると、後の推定処理を行う際に、画像を元来劣化させているぼけと、デノイジングによるぼけが混在したＰＳＦを推定してしまうことになる。このため、画像に与えるぼけが小さいデノイジング手法を用いることが好ましい。このような手法として、画像の周波数分解を利用したデノイジング処理を適用する。ここでは、周波数分解としてウェーブレット変換を用いた例を説明する。なお、この詳細については、「Ｄｏｎｏｈｏ  Ｄ．Ｌ．，‘‘Ｄｅ−ｎｏｉｓｉｎｇ  ｂｙ  ｓｏｆｔ−ｔｈｒｅｓｈｏｌｄｉｎｇ’’，ＩＥＥＥ  Ｔｒａｎｓ．  ｏｎ  Ｉｎｆ．  Ｔｈｅｏｒｙ，４１，３，ｐｐ．６１３−６２７」に記載されている。
【００８１】
ウェーブレット変換とは、局在する小さな波（ウェーブレット）を用いて画像の位置ごとに周波数解析を行い、信号を高周波成分と低周波成分に分解していく変換である。画像のウェーブレット変換では、画像の水平方向に対してウェーブレット変換を行って低周波成分と高周波成分に分解し、さらに分解して得られた低周波成分と高周波成分の垂直方向に対してウェーブレット変換を行う。ウェーブレット変換により、画像は４分割され、周波数帯域が互いに異なる４つのサブバンド画像に周波数分解される。このとき、左上の低周波帯域成分（スケーリング係数）のサブバンド画像をＬＬ１とし、右下の高周波帯域成分（ウェーブレット係数）のサブバンド画像をＨＨ１とする。また、右上（ＨＬ１）と左下（ＬＨ１）のサブバンド画像はそれぞれ、水平方向に高周波帯域成分をとって垂直方向に低周波帯域成分を取り出したものと、水平方向に低周波帯域成分をとって垂直方向に高周波帯域成分を取り出したものである。
【００８２】
さらに、サブバンド画像ＬＬ１をウェーブレット変換すると、画像サイズを半分にして、サブバンド画像ＬＬ２、ＨＬ２、ＬＨ２、ＨＨ２に分解することができ、分解して得られたサブバンド画像ＬＬに対して変換レベル回数だけ分解することができる。
【００８３】
ウェーブレット変換を用いてノイズ低減処理を行う方法として、スレッショルディングが知られている。これは、設定したしきい値より小さい量の成分をノイズとみなし、そのノイズを低減する方法である。ウェーブレット空間上でのしきい値処理は、サブバンド画像ＬＬ以外のサブバンド画像に対して行われ、以下の式（１１）に示されるように、しきい値以下の絶対値を有するウェーブレット係数ｗｓｕｂｂａｎｄ（ｘ，ｙ）を０に置き換えてデノイジングを行う。
【００８４】
【数１１】
0006541454.tif 000014
【００８５】
式（１１）において、ｘ、ｙはそれぞれ、画像の縦と横の座標であり、ρｓｕｂｂａｎｄは重みパラメータ、σはノイズの標準偏差である。ぼけ推定領域に含まれるノイズ量σは、ぼけ推定領域から測定または推定することにより求められる。ノイズが実空間および周波数空間において一様なホワイトガウシアンノイズである場合、ぼけ推定領域のノイズを、以下の式（１２）に示されるようなＭＡＤ（Ｍｅｄｉａｎ  Ａｂｓｏｌｕｔｅ  Ｄｅｖｉａｔｉｏｎ）から推定する手法が知られている。
【００８６】
【数１２】
0006541454.tif 000015
【００８７】
ＭＡＤは、ぼけ推定領域をウェーブレット変換して得られたサブバンド画像ＨＨ１におけるウェーブレット係数ｗＨＨ１のメディアン（中央値）を用いて求められる。標準偏差とＭＡＤとは、以下の式（１３）に示される関係であることから、ノイズ成分の標準偏差を推定することができる。
【００８８】
【数１３】
0006541454.tif 000016
【００８９】
なお、式（１２）、（１３）に代えて、撮影条件（撮影の際におけるＩＳＯ感度）に基づいてノイズ量σを取得してもよい。
【００９０】
続いてステップＳ２０２において、生成部１０８２は、ぼけ推定領域の解像度を低下させ、低解像ぼけ推定領域を生成する。ぼけ推定領域の解像度を低下させることにより、推定するぼけの解像度も同様に低下する。その結果、後述するぼけ推定の収束性が向上し、かつ、推定結果が最適解とは異なる局所解へ落ち込む可能性を低減することができる。ぼけ推定領域の解像度を下げる割合は、ダウンサンプリングパラメータに応じて決定される。ステップＳ２０２は、ループ処理（反復演算）によって複数回実行されるが、初回の実行時には、所定のダウンサンプリングパラメータが用いられる。２回目以降の実行の際には、後述のステップＳ２０９にて設定されたダウンサンプリングパラメータを用いて、低解像ぼけ推定領域を生成する。ループを反復する度に、解像度の低下量は小さくなり、低解像ぼけ推定領域の解像度は、ぼけ推定領域に近づく。すなわち、最初は低解像のぼけを推定し、その推定結果を新たな初期値として、解像度を少しずつ上げながら推定を繰り返す。これにより、局所解を回避して最適なぼけを推定することができる。なお、低解像ぼけ推定領域の解像度は、ぼけ推定領域の解像度以下であり、両者の解像度が一致することもある。
【００９１】
続いてステップＳ２０３において、生成部１０８２は、低解像ぼけ推定領域の信号変化量情報を取得し、信号変化量情報に基づいて後述のステップＳ２０４にて生成する補正信号勾配の数を決定する。信号変化量情報とは、低解像ぼけ推定領域内における信号変化（信号変化の量）、換言すると、低解像ぼけ推定領域内に存在するエッジ、テクスチャ、または、グラデーションの量を示す情報である。信号変化量情報は、ぼけている信号勾配から求めてもよいし、その信号勾配を補正して得られた補正信号勾配から求めてもよい。後者の場合、予め補正信号勾配（単数または複数の補正信号勾配）を生成し、その補正信号勾配から得られた信号変化量情報に基づいて、ステップＳ２０４にて生成する補正信号勾配の数を決定する。前述のように、補正信号勾配に信号変化が少ないほど、式（１）の連立１次方程式の線型独立性が低下する。このため生成部１０８２は、信号変化量情報に基づいて、生成する補正信号勾配の数を決定する。本実施例において、生成部１０８２は、信号変化の量が少ないほど、補正信号勾配の数が多くなるように決定する。
【００９２】
ここで、信号変化量情報の例を説明する。ただし本実施例は、これらの例に限定されるものではない。１つ目の例は、低解像ぼけ推定領域内の信号勾配（ぼけている）を用いる方法である。低解像ぼけ推定領域内の信号勾配を微分し、その二乗（絶対値でもよい）の総和を信号変化量情報とする。微分は、２次元的な微分（例えば、縦微分と横微分）である。画像を微分することにより各画素における信号変化量が求まり、その二乗和を取ることでぼけ推定領域内においてどの程度の信号変化が生じているが分かる。
【００９３】
好ましくは、微分値の二乗または絶対値を量子化し、その後にそれらの総和を取るようにする。量子化を行うことなく総和をとる場合、信号変化の領域が多く存在するのか、または、非常に大きな信号変化が僅かな領域に存在するのかを区別することができない。前者の場合、エッジやテクスチャなどの方向や強度にバリエーションが多く存在する。一方、後者の場合、そのようなバリエーションが少ない。エッジやテクスチャなどのバリエーションが少ないと、式（１）により表される連立１次方程式の中で独立な方程式の数が小さくなるため、ぼけの推定精度は低下する。このため、ぼけ推定領域内で信号変化が複数の領域に存在するのか、または、非常に大きな信号変化が僅かな領域に存在するのかを区別可能であることが好ましい。量子化を行うことにより、この区別を実現することができる。
【００９４】
例えば、微分量の二乗（または絶対値）が、所定のしきい値よりも大きい場合には１、所定のしきい値よりも小さい場合（平坦部の場合）には０というように、二値化してから総和を取る。これにより、エッジやテクスチャが多いのか否かを区別することができる。ここでは量子化として二値化による手法を用いているが、それより多いビット数で量子化しても構わない。以上の例は、ぼけている信号勾配の代わりに、補正信号勾配を生成してそれを用いる場合にも適用可能である。
【００９５】
または、予め１つの補正信号勾配を生成し、その局所領域の信号を行ベクトルで表し、異なる局所領域の行ベクトルを配列した行列を求め、その行列の階数を信号勾配情報としてもよい。これは、式（１）を行列の形に書き直した以下の式（１４）における行列Ｄｉの階数を示す。
【００９６】
【数１４】
0006541454.tif 000017
【００９７】
式（１４）において、ｂｉは低解像ぼけ推定領域の各信号値を成分とした縦ベクトル、ｋｉはぼけの各成分からなる縦ベクトル、ｎｉはノイズの各成分からなる縦ベクトル、Ｄｉは補正信号勾配の各信号値を成分とする行列である。ここで、予め生成した補正信号勾配は１つである。このため、行列Ｄｉは低解像ぼけ推定領域の画素数だけ行数を持ち、その列数はぼけの成分数となる。行列Ｄｉの階数は、連立１次方程式の線型独立性を如実に表しているため、これを信号変化量情報とすることが好ましい。前述の２つの例では、同じエッジやテクスチャが複数存在しているかが分からないが、行列Ｄｉの階数ではそれらも評価することができる。
【００９８】
続いてステップＳ２０４において、生成部１０８２は、低解像ぼけ推定領域における信号勾配のぼけを補正し、補正信号勾配を生成する（ぼけ推定領域における信号に関する情報の補正処理）。補正信号勾配の生成には、前述のように、式（２）を用いた最適化やその他の超解像、逆フィルタなどを用いる。ぼけ補正に用いるＰＳＦは、前のループにおけるステップＳ２０５にて推定した結果（ぼけ）を用いる。ループの初回では、適当な形状のＰＳＦ（Ｇａｕｓｓ分布や１次元のラインなど）を用いる。生成する補正信号勾配の数は、ステップＳ２０３にて決定した数であり、本実施例では、ループのうち少なくとも１回は複数の補正信号勾配を生成し、ぼけ推定の高精度化を図る。好ましくは、ぼけの推定精度を向上させるため、全てのループ（反復演算）に対して複数の補正信号勾配を生成する。複数の補正信号勾配を生成する際には、共通のＰＳＦ（ぼけ）を用い、互いに異なる補正処理または補正パラメータに基づいて補正信号勾配を生成する。また、複数の補正信号勾配を生成する際に、ＰＳＦを微小に変化させてもよい。ただし、ＰＳＦを変化させると後のぼけ推定に誤差が生じる。このため、生成される複数の補正信号勾配は全て、共通のＰＳＦで生成されることが好ましい。
【００９９】
続いてステップＳ２０５において、生成部１０８２は、低解像ぼけ推定領域の信号勾配と、補正信号勾配とに基づいて、ぼけを推定する（ぼけの推定処理）。ぼけの推定には、補正信号勾配の数に応じて、式（６）、式（６ａ）、または、式（１０）、（１０ａ）の最適化などを用いる。ステップＳ２０４にて複数の補正信号勾配が生成されている場合、少なくとも２つ以上の補正信号勾配に対して共通のぼけを推定する。共通のぼけの推定には、式（１０）または式（１０ａ）を用いるが、式（１０）、（１０ａ）中のＧは、ステップＳ２０３にて決定した補正信号勾配の数と一致してもよいし、それ未満の値でもよい。これについて、図５および図６を参照して説明する。
【０１００】
図５および図６は、本実施例におけるぼけ推定の一例を示す図であり、ぼけ推定領域の信号勾配２０１から、６つの補正信号勾配２０２ａ〜２０２ｆを生成し、それらの補正信号勾配２０２ａ〜２０２ｆに基づいてぼけ２０３を推定する様子を示している。ただし、補正信号勾配の生成数は、これに限定されるものではない。図５においては、６つの補正信号勾配２０２ａ〜２０２ｆの全てに対して、共通のぼけ２０３を推定している。この際、式（１０）または式（１０ａ）のＧは６となる。一方、図６においては、６つの補正信号勾配を３つのグループに分けて、それぞれのグループに対して共通のぼけ２０３ａ〜２０３ｃを推定し、そこからぼけ２０３を求めている。ぼけ２０３ａは、２つの補正信号勾配２０２ａ、２０２ｂから推定され、この際、式（１０）または式（１０ａ）のＧは２となる。ぼけ２０３ｂは、３つの補正信号勾配２０２ｃ〜２０２ｅから推定され、この際、Ｇは３となる。ぼけ２０３ｃは、１つの補正信号勾配２０２ｆから推定され、この際、推定には式（６）または式（６ａ）が用いられる。最終的なぼけ２０３は、ぼけ２０３ａ〜２０３ｃの加重平均などで求めることができる。図５の場合、６つの補正信号勾配に対して、推定されたぼけの数は１つである。一方、図６の場合、６つの補正信号勾配に対して、推定されたぼけの数は３つである（図６のぼけ２０３は、推定されたぼけ２０３ａ〜２０３ｃから求められたぼけであるため、数に入れない）。
【０１０１】
このように、生成された補正信号勾配が複数存在する場合、推定されるぼけの数は、補正信号勾配の数よりも少なくなる。本実施例は、ぼけの推定において、図５および図６のいずれの場合も含むが、一度に考える連立方程式の数が多くなる図５の場合のほうが、ぼけ推定の精度が向上するため、より好ましい。
【０１０２】
また、ステップＳ２０４にて互いに異なる補正処理または補正パラメータを用いて複数の補正信号勾配が生成されている場合、それぞれの補正信号勾配を重み付けしてぼけを推定することが好ましい。これは、式（１０）または式（１０ａ）において、重みｔｇを一定の値でなく、変化させることに相当する。これは、補正処理や補正パラメータに応じて、補正信号勾配の補正精度またはリンギングなどの弊害の発生しやすさが変化するためである。補正精度が高く、弊害の出にくい補正処理または補正パラメータの補正信号勾配に対して、重みをより大きくすることが好ましい。
【０１０３】
続いてステップＳ２０６において、生成部１０８２は、ステップＳ２０５（第２の推定処理）にて推定されたぼけ（第２のぼけ）と、前ループのステップＳ２０５（第１の推定処理）にて推定されたぼけ（第１のぼけ）とを比較する。具体的には、生成部１０８２は、第１のぼけと第２のぼけとの差異が所定値（所定のしきい値）以下であるか否かを判定する。この差異が所定値以下である場合、ステップＳ２０８へ進む。一方、この差異が所定値よりも大きい場合、ステップＳ２０７へ進む。ステップＳ２０２にて説明したとおり、本実施例においては、ぼけ推定領域内の解像度を変更しながら、補正信号勾配の生成（補正処理）とぼけの推定（推定処理）をループさせる。ループ初期の低解像ぼけ推定領域は、画素数が少ないため、ぼけの推定結果は収束しやすい。一方、ループを繰り返して低解像ぼけ推定領域の画素数が多くなると、結果は収束しにくくなる。さらにエッジなどの信号変化が不足していると、推定したぼけは正しい値から大きく外れてしまう。これを回避するため、生成部１０８２は、現在のループで推定したぼけと前のループのぼけの差異を評価する。ここでは、現ループと前ループにおけるぼけのいずれか一方をバイリニア補間により、他方と同じ成分数にリサイズし、各成分の二乗和を取って差異を評価する。ただし本実施例は、これに限定されるものではなく、他の手法を用いてもよい。両方のぼけの差異が所定値よりも大きい場合、現ループにおいてぼけ推定に必要な信号変化の情報が不足して推定精度が低下していると考えられる。このため、補正信号勾配の数を増やし、ぼけ推定をやり直す。ループ初回においては、前のループの推定ぼけが存在しないため、このステップを実行することなくステップＳ２０８へ進む。
【０１０４】
ステップＳ２０７において、生成部１０８２は、生成する補正信号勾配の数を増やす。この理由は、前述のとおりである。そしてステップＳ２０４に戻り、生成部１０８２は、新たに複数の補正信号勾配を生成する。この際、生成部１０８２は、既に生成されている補正信号勾配を新たに生成することはせず、そのまま流用してもよい。
【０１０５】
ステップＳ２０８において、生成部１０８２は、反復（ループ）演算が完了したか否か、すなわちぼけ推定領域と低解像ぼけ推定領域の解像度差が所定値よりも小さいか否か、を判定する。この判定は、ぼけ推定領域の解像度と低解像ぼけ推定領域の解像度とを比較することにより行われる。両者の解像度差が所定値よりも小さい場合、反復演算を終了する。そして、ステップＳ２０５にて推定されたぼけを最終的な推定ぼけとして、ステップＳ２１０へ進む。両者の解像度差が所定値よりも小さいか否かは、両者の解像度差の絶対値が所定値よりも小さいか否かに限定されるものではなく、両者の解像度の比が所定値よりも１に近いか否かなどで判定してもよい。所定の条件を満たさない場合（例えば、両者の解像度差の絶対値が所定値よりも大きい場合）、推定されたぼけの解像度がまだ不充分であるため、ステップＳ２０９へ進み、反復演算を行う。
【０１０６】
ステップＳ２０９において、生成部１０８２は、ステップＳ２０２にて用いられるダウンサンプリングパラメータを設定する。ステップＳ２０２〜Ｓ２０７を反復する過程で解像度を上げていくため、ここでは前のループよりもダウンサンプリングの度合いを弱めるように（解像度低下量を下げるように）パラメータを設定する。また、反復演算では、前ループのステップＳ２０５にて推定されたぼけを新たな初期解として用いるが、その際、ぼけの解像度を上げる必要がある。解像度の向上には、バイリニア補間やバイキュービック補間を用いることが好ましい。
【０１０７】
ステップＳ２１０において、生成部１０８２は、ステップＳ２０５にて推定されたぼけ（推定ぼけ）のデノイジングを行う。式（１）に示されるように、ぼけ画像（ぼけ推定領域）にはノイズが存在するため、その影響を受けて推定ＰＳＦにもノイズが発生する。このため、生成部１０８２はデノイジング処理を行う。デノイジング処理を行う際は、スレッショルディング処理、孤立点を除去するオープニング処理、または、平滑化フィルタや前述の周波数分解を用いたデノイジング処理などを用いることが好ましい。
【０１０８】
本実施例によれば、ぼけ画像における信号変化が少ない領域からでも、高精度にぼけを推定可能な画像処理システムを提供することができる。
【実施例２】
【０１０９】
次に、図７を参照して、本発明の実施例２における画像処理システムについて説明する。図７は、本実施例における画像処理システム３００のブロック図である。
【０１１０】
本実施例の画像処理システム３００は、ぼけ補正部１０８を有する画像処理装置１０５に代えて、推定ぼけ生成部３０８を有する画像処理装置３０５を備えている点で、実施例１の画像処理システム１００とは異なる。本実施例の画像処理システム３００は、実施例１の画像処理システム１００よりも簡易な処理により、ぼけ推定を行うことが可能である。推定ぼけ生成部３０８は、撮像装置１０１にて撮像されたぼけ画像からぼけ成分を推定して出力する。画像処理システム３００のそれ以外の部位は、実施例１の画像処理システム１００と同様であるため、それらの説明は省略する。
【０１１１】
次に、図８を参照して、推定ぼけ生成部３０８で行われる画像処理に関して説明する。図８は、本実施例の画像処理方法を示すフローチャートである。図８の各ステップは、推定ぼけ生成部３０８の取得部３０８１（取得手段）および生成部３０８２（生成手段）により実行される。
【０１１２】
図８のステップＳ３０１は推定ぼけ生成部３０８の取得部３０８１により実行され、ステップＳ３０２は推定ぼけ生成部３０８の生成部３０８２により実行される。なお、ステップＳ３０１、Ｓ３０２は、図４を参照して説明した実施例１のステップＳ１０１、Ｓ１０２とそれぞれ同様である。
【０１１３】
続いてステップＳ３０３において、生成部３０８２は、ぼけ推定領域の信号勾配を補正し、複数の補正信号勾配を生成する。複数の補正信号勾配の生成には、互いに異なる補正処理または補正パラメータが用いられる。本実施例において、生成部３０８２は、常に所定の数だけ補正信号勾配を生成する。ステップＳ３０３における補正処理は、実施例１のステップＳ２０４にて説明した補正処理と同様である。続いてステップＳ３０４において、生成部３０８２は、ぼけ推定領域の信号勾配と複数の補正信号勾配とに基づいて、ぼけを推定する。ステップＳ３０４における推定処理は、実施例１のステップＳ２０５にて説明した推定処理と同様である。
【０１１４】
続いてステップＳ３０５において、生成部３０８２は、ステップＳ３０４にて推定したぼけ（推定結果）が収束したか否かを判定する。ぼけが収束した場合、ステップＳ３０６へ進む。一方、ぼけが収束していない場合、ステップＳ３０３へ戻り、生成部３０８２は、ステップＳ３０４にて推定したぼけを用いて、ステップＳ３０３にて複数の補正信号勾配を新たに生成する。推定したぼけが収束したか否かの判定は、例えば、補正信号勾配を推定したぼけで劣化させたものと、ぼけ推定領域の信号勾配との差または比を求め、所定値と比べればよい。または、ステップＳ３０４にて共役勾配法などの繰り返し演算でぼけを推定している場合、この繰り返し演算によるぼけの更新量が所定値よりも小さいか否かにより判定してもよい。推定したぼけが収束している場合、生成部３０８２は、推定されたぼけをぼけ推定領域における最終的な推定ぼけとする。
【０１１５】
続いてステップＳ３０６において、生成部３０８２は、ぼけ画像における所定の領域の全て（例えば、ぼけ画像の全体）に対して、推定ぼけの生成が完了したか否かを判定する。ステップＳ３０６は、実施例１のステップＳ１０４と同様である。推定ぼけの生成が完了している場合、ステップＳ３０７へ進む。一方、推定ぼけの生成が完了していない場合、ステップＳ３０２へ戻る。
【０１１６】
ステップＳ３０７において、生成部３０８２は、生成された推定ぼけを出力する。ぼけ画像における所定の領域内に複数のぼけ推定領域が含まれる場合（すなわち、Ｓｈｉｆｔ−ｖａｒｉａｎｔの場合）、複数の推定ぼけが存在する。このため生成部３０８２は、複数の推定ぼけを出力する。なお、複数の推定ぼけが存在する場合、それらを求める途中で逐次出力しても構わない。出力された推定ぼけは、ぼけ画像の補正、撮影した光学系の光学性能測定、または、撮影の際の手ぶれの解析などに用いることができる。
【０１１７】
本実施例によれば、ぼけ画像における信号変化が少ない領域からでも、高精度にぼけを推定可能な画像処理システムを提供することができる。
【実施例３】
【０１１８】
次に、図９および図１０を参照して、本発明の実施例３における撮像システムについて説明する。図９は、本実施例における撮像システム４００のブロック図である。図１０は、撮像システム４００の外観図である。
【０１１９】
撮像システム４００は、撮像装置４０１、ネットワーク４０２、および、サーバ４０３（画像処理装置）により構成される。撮像装置４０１とサーバ４０３とは無線で接続されており、撮像装置４０１からの画像はサーバ４０３へ転送され、サーバ４０３がぼけの推定および補正を行う。
【０１２０】
サーバ４０３は、通信部４０４、記憶部４０５、および、ぼけ補正部４０６（画像処理部）を有する。サーバ４０３の通信部４０４は、ネットワーク４０２を介して撮像装置４０１と接続されている。本実施例において、撮像装置４０１とサーバ４０３とは無線で接続されているが、これに限定されるものではなく、有線で接続されてもよい。サーバ４０３の通信部４０４は、撮像装置４０１からぼけ画像を受信するように構成されている。撮像装置４０１により撮影が行われると、ぼけ画像（入力画像または撮影画像）は、自動的または手動でサーバ４０３に入力され、記憶部４０５およびぼけ補正部４０６に送られる。記憶部４０５は、ぼけ画像およびそのぼけ画像を撮影した撮影条件に関する情報を記憶する。ぼけ補正部４０６は、ぼけ画像に基づいてぼけ（推定ぼけ）を推定する。そしてぼけ補正部４０６は、推定ぼけに基づいて、ぼけ補正画像を生成する。ぼけ補正画像は、記憶部４０５に記憶されるか、または、通信部４０４を介して撮像装置４０１へ送られる。
【０１２１】
本実施例の画像処理方法を実現するため、ソフトウエア（画像処理プログラム）を、ネットワーク、またはＣＤ−ＲＯＭなどの記憶媒体を介して、サーバ４０３に供給することもできる。この際、画像処理プログラムは、サーバ４０３のコンピュータ（またはＣＰＵやＭＰＵなど）により読み出されて、サーバ４０３の機能を実行する。
【０１２２】
なお、ぼけ補正部４０６で行われる処理は、図１および図４を参照して説明した実施例１の画像処理方法と同様であるため、その説明は省略する。本実施例によれば、ぼけ画像における信号変化の量が少ない領域からでも、高精度にぼけを推定可能な撮像システムを提供することができる。
【０１２３】
このように各実施例において、画像処理装置（画像処理装置１０５、３０５、サーバ４０３）は、取得手段（取得部１０８１、３０８１）および生成手段（生成部１０８２、３０８２）を有する。取得手段は、ぼけ画像（撮影画像または入力画像）を取得する。生成手段は、ぼけ画像における少なくとも一部のぼけ推定領域を取得し、ぼけ推定領域から推定ぼけを生成する。生成手段は、ぼけ推定領域における信号に関する情報に含まれるぼけ（信号勾配）を補正して補正信号に関する情報（補正信号勾配）を生成する補正処理（Ｓ２０４）を行う。また生成手段は、信号に関する情報と補正信号に関する情報とに基づいてぼけを推定する推定処理（Ｓ２０５）を行う。また生成手段は、補正処理および推定処理を反復する反復演算処理（Ｓ２０４〜Ｓ２０７）を行って推定ぼけを生成する。また生成手段は、反復演算処理における少なくとも一回の補正処理において、補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する。
【０１２４】
好ましくは、生成手段は、複数の演算式の種類または複数の演算式のパラメータ（係数）を変更して複数の補正信号に関する情報を生成する。また好ましくは、生成手段は、信号に関する情報として、ぼけ推定領域における輝度分布または輝度分布の微分値に関する情報（信号勾配）を補正する。また好ましくは、生成手段は、補正処理において複数の補正信号に関する情報を生成した場合、推定処理において、少なくとも２つ以上の補正信号に関する情報に関して共通のぼけを推定する。このとき、推定処理において推定されるぼけの総数は、複数の補正信号に関する情報の総数よりも少ない。より好ましくは、生成手段は、補正処理において複数の補正信号に関する情報を生成した場合、推定処理において、複数の補正信号に関する情報の全てに関して共通の１つのぼけを推定する。
【０１２５】
好ましくは、生成手段は、同一のぼけを用いて信号に関する情報を補正し、複数の補正信号に関する情報を生成する。また好ましくは、生成手段は、反復演算処理の全ての補正処理において、補正信号に関する情報として、互いに異なる複数の演算式を用いて複数の補正信号に関する情報を生成する。また好ましくは、生成手段は、補正処理において複数の補正信号に関する情報を生成した場合、推定処理において、複数の補正信号に関する情報のそれぞれに重み付けを施してぼけを推定する（Ｓ２０４）。また好ましくは、生成手段は、反復演算処理の推定処理として、第１の推定処理および第１の推定処理よりも後に第２の推定処理を行う。また生成手段は、第１の推定処理において推定された第１のぼけと、第２の推定処理において推定された第２のぼけとの差異を判定する。そして生成手段は、差異が所定のしきい値よりも大きい場合、複数の補正信号に関する情報の数を増加して、再度、推定処理を行う（Ｓ２０６）。
【０１２６】
好ましくは、生成手段は、ぼけ推定領域における信号変化の量に関する情報（信号変化量情報）に基づいて、複数の補正信号に関する情報の数を決定する（Ｓ２０３）。好ましくは、生成手段は、ぼけ推定領域の解像度を低下させて低解像ぼけ推定領域を生成する。また生成手段は、低解像ぼけ推定領域において、信号に関する情報に含まれるぼけを補正して補正信号に関する情報を生成する補正処理と、信号に関する情報と補正信号に関する情報とに基づいてぼけを推定する推定処理とを反復する反復演算処理を行う。これにより、生成手段は推定ぼけを生成する。また生成手段は、反復演算処理の過程において、低解像ぼけ推定領域の解像度をぼけ推定領域の解像度に近づける（ぼけ推定領域の解像度低下量を下げる：Ｓ２０９）。
【０１２７】
好ましくは、生成手段は、ぼけ推定領域のデノイジングを行い、ぼけ推定領域に含まれるノイズ量を取得し、ぼけ推定領域を周波数分解して周波数分解ぼけ推定領域を生成する。そして生成手段は、ノイズ量に基づいて、周波数分解ぼけ推定領域に対するデノイジング処理を行い、デノイジング処理後の周波数分解ぼけ推定領域を再合成する（Ｓ２０１）。また好ましくは、画像処理装置は、推定ぼけを用いてぼけ画像の少なくとも一部を補正する補正手段（補正部１０８３）を更に有する。
【０１２８】
（その他の実施例）
本発明は、上述の実施例の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサーがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。
【０１２９】
各実施例によれば、ぼけ画像における信号変化が少ない領域からでも、高精度にぼけを推定可能な画像処理装置、撮像装置、画像処理方法、画像処理プログラム、および、記憶媒体を提供することができる。
【０１３０】
以上、本発明の好ましい実施例について説明したが、本発明はこれらの実施例に限定されず、その要旨の範囲内で種々の変形及び変更が可能である。
【符号の説明】
【０１３１】
１０５、３０５  画像処理装置
１０８１、３０８１  取得部（取得手段）
１０８２、３０８２  生成部（生成手段）
【図１】
0006541454.tif 000018
【図２】
0006541454.tif 000019
【図３】
0006541454.tif 000020
【図４】
0006541454.tif 000021
【図５】
0006541454.tif 000022
【図６】
0006541454.tif 000023
【図７】
0006541454.tif 000024
【図８】
0006541454.tif 000025
【図９】
0006541454.tif 000026
【図１０】
0006541454.tif 000027