JP2017123532A - 画像処理装置、撮像装置、画像処理方法、画像処理プログラムおよび記録媒体 - Google Patents

Abstract

【課題】撮像光学系と撮像素子の任意の組み合わせにおける画像内のシフトバリアントな劣化特性を高速かつ良好に復元可能な画像処理装置、撮像装置、画像処理方法、画像処理プログラムおよび記録媒体を提供すること。
【解決手段】撮影画像を取得する第１の取得手段６１と、撮像光学系および撮像素子の情報を取得する第２の取得手段６２と、撮像素子の画素ピッチと第１のタップ数に基づいて、撮像光学系に関する第１の点像強度分布関数を生成する第１の生成手段６６と、第１の点像強度分布関数に関する分布指標値を取得する第３の取得手段６７と、第１点像強度分布関数と、分布指標値に基づいて算出される第１のタップ数以下の第２のタップ数に基づいて第２の点像強度分布関数を生成する第２の生成手段６８と、第２の点像強度分布関数に基づいて撮影画像を復元する復元手段６９と、を有する。
【選択図】図６

Description

本発明は、画像処理装置、撮像装置、画像処理方法、画像処理プログラムおよび記録媒体に関する。

撮像光学系により撮影された被写体は、撮像光学系で発生する回折や収差等の影響により、１点から発生した光が１点に収束せず、微小な広がりを持つ。このような微小な広がりを持つ分布をＰＳＦ（点像強度分布関数）という。このとき、撮影画像はＰＳＦが畳み込まれて形成されるため、解像度が劣化する。近年、撮影画像を電子データとして保持することが一般的になり、画像処理によって光学系による画像劣化を補正する技術（画像回復処理、画像復元処理）が提案されている。

以下、一般的な画像回復処理について説明する。実空間（ｘ，ｙ）上で、撮像光学系による劣化を受ける前の画像をｆ（ｘ，ｙ）、ＰＳＦをｈ（ｘ，ｙ）、劣化した画像をｇ（ｘ，ｙ）とすると、これらは以下の式（１）で表される。

ｇ（ｘ，ｙ）＝∫∫ｆ（Ｘ，Ｙ）＊ｈ（ｘ−Ｘ，ｙ−Ｙ）ｄＸｄＹ （１）
式（１）にフーリエ変換を施し、実空間（ｘ，ｙ）から周波数空間（ｕ，ｖ）への変換を行うと、式（１）は以下の式（２）で表される。

Ｇ（ｕ，ｖ）＝Ｆ（ｕ，ｖ）＊Ｈ（ｕ，ｖ） （２）
ここで、Ｆ（ｕ，ｖ）はｆ（ｘ，ｙ）のフーリエ変換、Ｇ（ｕ，ｖ）はｇ（ｘ，ｙ）のフーリエ変換、Ｈ（ｕ，ｖ）はｈ（ｘ，ｙ）のフーリエ変換である。また、式（２）を変形すると、以下の式（３）で表される。

Ｆ（ｕ，ｖ）＝Ｇ（ｕ，ｖ）／Ｈ（ｕ，ｖ） （３）
式（３）は、周波数空間上で、劣化画像ｇ（ｘ，ｙ）のフーリエ変換Ｇ（ｕ，ｖ）をＰＳＦであるｈ（ｘ，ｙ）のフーリエ変換Ｈ（ｕ，ｖ）で割ることにより、劣化を受ける前の画像ｆ（ｘ，ｙ）のフーリエ変換Ｆ（ｕ，ｖ）を得ることができることを意味する。したがって、Ｆ（ｕ，ｖ）にフーリエ逆変換を施せば、劣化を受ける前の画像ｆ（ｘ，ｙ）を取得することができる。

しかしながら、実際にこのような処理により劣化を受ける前の画像を取得する場合、特にＨ（ｕ，ｖ）上に０近傍の値がある場合、１／Ｈ（ｕ，ｖ）が非常に大きな値となってしまう。この問題を解決するために以下の式（４）で表されるウィナーフィルタを用いることが知られている。

Ｒ（ｕ，ｖ）＝１／Ｈ（ｕ，ｖ）＊｜Ｈ（ｕ，ｖ）｜＾２／（｜Ｈ（ｕ，ｖ）｜＾２＋Γ） （４）
ここで、Ｈ（ｕ，ｖ）は光学伝達関数（ＯＴＦ）であり、Γは処理上の０除算を回避するための定数項である。式（４）を、劣化画像ｇ（ｘ，ｙ）のフーリエ変換Ｇ（ｕ，ｖ）に乗算すれば、光学系やセンサ等による劣化を補正することができる。すなわち、ウィナーフィルタを用いた画像復元処理により復元された画像の周波数空間情報Ｆ＾（ｕ，ｖ）は、以下の式（５）で表される。

Ｆ＾（ｕ，ｖ）＝Ｇ（ｕ，ｖ）＊Ｗ（ｕ，ｖ） （５）
式（５）が成立するのは、ＰＳＦがシフトインバリアントな領域に限定される。ＰＳＦがシフトインバリアントな領域とは、簡単に言えばＰＳＦが変化しない領域のことである。デジタルカメラに用いられる非常に高画角なレンズやズームレンズ等では、シフトインバリアントな領域は非常に狭い。これを回避するために実空間で領域ごとに、以下の式（６）で表されるコンボリューション演算で処理する手法などがある。

ｆ＾（ｘ，ｙ）＝∫∫ｇ（Ｘ，Ｙ）＊ｒ（ｘ−Ｘ，ｙ−Ｙ）ｄＸｄＹ （６）
ここで、ｆ＾（ｘ，ｙ）、ｒ（ｘ，ｙ）はそれぞれ、実空間上の復元された画像、実空間上のウィナーフィルタである。

また、画像復元処理の別の手法として、非特許文献１のＬｕｃｙ−Ｒｉｃｈａｒｄｓｏｎ法（以下、ＬＲ法という）等が知られている。ＬＲ法は、ウィナーフィルタとは異なり複数回繰り返し処理を行うことで画像を復元する手法であり、繰り返し処理による計算を実行する必要があるため、ウィナーフィルタよりも計算コストが高い。ＬＲ法は、実空間において以下の式（７）で表される。

ここで、ｆ_ｎは更新後の画像、ｆ_ｎ−１は更新前の画像である。ｈは、光学系の劣化特性であり、例としてＰＳＦがあげられる。また、ｇは、撮影された画像である。

ＬＲ法は、画像ｆ_ｎ−１を逐次更新し、例えば、ｆ_ｎとｆ_ｎ−１の差分値等の所定の評価基準に基づいて処理を打ち切り、ｆ_ｎを出力する。ＬＲ法の開始時では、ｆ_０＝ｇなどとして処理を開始すればよい。

これらの処理に要する時間は、画像の画素数と畳み込まれるフィルタのタップ数が支配的である。これは式（６）においては画像ｇ（ｘ，ｙ）、フィルタｒ（ｘ，ｙ）であり、式（７）においては画像ｇ（ｘ，ｙ）またはｆ_ｎ−１（ｘ，ｙ）、フィルタｈ（ｘ，ｙ）である。フィルタｈ（ｘ，ｙ）は、ＰＳＦそのものである。

特開２０１３―２５４７３号公報

"ＪＯＵＲＮＡＬ ＯＦ ＴＨＥ ＯＰＴＩＣＡＬ ＳＯＣＩＥＴＹ ＯＦ ＡＭＥＲＩＣＡ"，（Ａｍｅｒｉｃａ），Ｔｈｅ Ｏｐｔｉｃａｌ Ｓｏｃｉｅｔｙ，ＪＡＮＵＡＲＹ １９７２，ＶＯＬＵＭＥ ６２，ＮＵＭＢＥＲ １，ｐ．５５−５９

一般的なデジタルカメラで撮影された画像は、画像の中心から周辺部になるにつれて撮像光学系の収差が大きくなり、結果としてＰＳＦのスポットサイズが大きくなる。また、レンズ交換式カメラのように、センサピッチが異なる様々な撮像素子と撮像光学系の組み合わせは無数にある。したがって、様々な撮像素子と撮像光学系の組み合わせにおいて、フィルタのタップ数を最適に設定して画像復元処理を高速化することが大きな課題となる。特許文献１には、レンズ特性により画像劣化が画面内位置で異なる場合にも、各位置の劣化状態に応じて最適に復元処理を行うとともに、演算量を低減させる手法が開示されている。しかしながら、特許文献１の手法では、測定したＰＳＦをガウス分布の縦横の分散値に変換してしまうため、大きなコマ収差のように、ガウス分布では高精度な近似が不可能な劣化特性に対する補正効果が大きく低下する。また、ＰＳＦの縦横の分散値を用いた手法は、座標の単位が変わってしまうため、撮像光学系と撮像素子の組み合わせが無数にあるレンズ交換式のカメラには適用できない。すなわち、撮像光学系と撮像素子の様々な組み合わせにおける画像復元処理において、演算時間を低減するために最適なフィルタタップサイズを設定し、撮影画像内のシフトバリアントな劣化特性を良好に復元する手法はこれまで開示されていない。

このような課題に鑑みて、撮像光学系と撮像素子の任意の組み合わせにおける画像内のシフトバリアントな劣化特性を高速かつ良好に復元可能な画像処理装置、撮像装置、画像処理方法、画像処理プログラムおよび記録媒体を提供することを目的とする。

本発明の一側面としての画像処理装置は、撮影画像を取得する第１の取得手段と、撮像光学系および撮像素子の情報を取得する第２の取得手段と、前記撮像素子の画素ピッチと第１のタップ数に基づいて、前記撮像光学系に関する第１の点像強度分布関数を生成する第１の生成手段と、前記第１の点像強度分布関数に関する分布指標値を取得する第３の取得手段と、前記第１点像強度分布関数と、前記分布指標値に基づいて算出される前記第１のタップ数以下の第２のタップ数に基づいて第２の点像強度分布関数を生成する第２の生成手段と、前記第２の点像強度分布関数に基づいて前記撮影画像を復元する復元手段と、を有することを特徴とする。

本発明の他の側面としての画像処理方法は、撮影画像を取得するステップと、撮像光学系および撮像素子の情報を取得するステップと、前記撮像素子の画素ピッチと第１のタップ数に基づいて、前記撮像光学系に関する第１の点像強度分布関数を生成するステップと、前記第１の点像強度分布関数に関する分布指標値を取得するステップと、前記分布指標値に基づいて前記第１のタップ数以下の第２のタップ数を算出するステップと、前記第１点像強度分布関数と前記第２のタップ数に基づいて第２の点像強度分布関数を生成するステップと、前記第２の点像強度分布関数に基づいて前記撮影画像を復元するステップと、を有することを特徴とする。

本発明によれば、撮像光学系と撮像素子の任意の組み合わせにおける画像内のシフトバリアントな劣化特性を高速かつ良好に復元可能な画像処理方法、画像処理装置、撮像装置、画像処理プログラムおよび記録媒体を提供することできる。

実施例１の画像復元処理を示すフローチャートである。 分布指標値を示す図である。 第２のタップ数に対応する第２のＰＳＦを示す図である。 実施例２と実施例３の画像復元処理の説明図である。 画像処理プログラムをインストールして実行する装置の外観図である（実施例４）。 実施例５の画像処理装置の構成図である。 実施例６の撮像装置の構成を示すブロック図である。

以下、本発明の実施例について、図面を参照しながら詳細に説明する。各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。

本実施例の画像処理方法について、図１を参照して説明する。図１は、本実施例の画像復元処理を示すフローチャートである。本実施例の画像復元処理は、ソフトウエア、ハードウエア上で動作するコンピュータプログラムとしての画像処理プログラムにしたがって実行される。画像処理プログラムは、例えば、コンピュータに読み取り可能な記録媒体に記録されていてもよい。

画像復元処理が開始されると、ステップＳ１０１では、画像復元処理が施される撮影画像が取得される。

ステップＳ１０２では、撮影時の撮影条件（撮像光学系および撮像素子の情報）が取得される。撮影情報は、Ｅｘｉｆ情報等の撮影画像に付随する情報から取得してもよいし、撮像装置等から直接取得してもよい。撮像光学系の情報とは、撮影時のレンズを特定可能なレンズ識別番号（以下、レンズＩＤという）や撮影時の焦点距離、Ｆ値、または被写体距離のことである。レンズＩＤ、焦点距離、Ｆ値、および被写体距離を組み合わることで、撮像光学系の光学伝達関数を特定することができる。また、撮像素子の情報は、画素ピッチ、カラーフィルタ配列、または光学ローパスフィルタに関する情報のことである。

撮像光学系の光学伝達関数とは、折り返し信号が発生しない程度のサンプリングピッチで離散化されたデータである。光学伝達関数は、例えば、レンズの設計データ等からコンピュータ上のシミュレーションで算出されることが好ましいが、画素ピッチの小さい撮像素子を用いて測定してもよい。このときの画素ピッチは、折り返し信号が発生しない程度の画素ピッチであることが好ましい。コンピュータシミュレーションでは、想定する光源の分光特性や撮像素子のカラーフィルタの分光特性を用いて、Ｒ，Ｇ，Ｂに対応する光学伝達関数を生成しておくことが好ましい。撮像光学系の光学伝達関数は、焦点距離、Ｆ値、被写体距離を決定した場合に撮像素子上の位置と主にＲ，Ｇ，Ｂのカラーフィルタの分光特性の組み合わせの数だけ必要となる。例えば、撮像素子上の代表点１００点とＲ，Ｇ，Ｂのカラーフィルタの分光特性の組み合わせの数は３００となる。したがって、撮像光学系の光学伝達関数を特定するとは、この３００点を特定することを意味している。

次に、このように特定した光学伝達関数を取得する。このとき、使用される撮像素子のサイズによって、光学伝達関数の取得範囲を変えることが望ましい。この取得範囲は、撮像光学系によって決定される範囲と、撮像素子が光学像を取得できる範囲の比率等に基づいて決定することができる。撮像光学系によって決定される範囲とは、設計時に収差補正が行われている範囲である。

取得した光学伝達関数は、そのままでは画像復元処理に使用できないので、それぞれ修正される。本実施例では、様々な撮像素子に対して、その素子が持つ帯域や画素の形状などの周波数特性を考慮して、光学伝達関数を修正する。撮像素子の特性として重要なものとして、撮像素子の素子サイズと画素ピッチ、カラーフィルタ配列情報等がある。撮像素子の素子サイズは光学空中像を切り出すイメージサークルを指定する情報であり、画素ピッチは光学伝達関数の折り返し信号に関する情報である。また、光学ローパスフィルタに関する特性を考慮してもよい。光学ローパスフィルタに関する特性とは、光学ローパスフィルタの有無や、光学ロータパスフィルタにより点像強度分布の分離幅、分離方向、分離方法といった光学伝達関数が受ける特性に関する情報のことである。また、単板センサの場合、Ｒ，Ｇ，Ｂの配列並びもベイヤー配列のサンプリング間隔における折り返し信号についての情報として使用できる。また、ベイヤー配列になっている撮像素子の場合、Ｒ，Ｇ，Ｂについてそれぞれのサンプリングピッチが異なる場合や特殊な配置になっている場合もがあり、それらを考慮に入れることが可能である。これらは撮影条件の情報として取得される。

ここで、撮像素子の特性を考慮して修正された光学伝達関数Ｈ_Ｓ（ｕ，ｖ）は、以下の式（８）で表される。

ここで、Ｈ（ｕ，ｖ）は撮像光学系の光学伝達関数、ＯＬＰＦ（ｕ，ｖ）は光学ローパスフィルタの周波数特性、ＰＩＸ（ｕ，ｖ）は撮像素子の画素形状による周波数特性、[ ]_Ｓは撮像素子の帯域内の周波数特性の取得である。また、ｕ，ｖはそれぞれ、Ｕ，Ｖタップに離散化された状態である。

撮像素子のナイキスト周波数をそれぞれＵｎｙｑ，Ｖｎｙｑとすると、実空間サイズは水平方向に対してＵ／（２＊Ｕｎｙｑ），垂直方向に対してＶ／（２＊Ｖｎｙｑ）となる。ナイキスト周波数Ｕｎｙｑ，Ｖｎｙｑの単位は、慣例的に本／ｍｍが使用されることが多いが、特に限定はされない。ナイキスト周波数Ｕｎｙｑ，Ｖｎｙｑは、使用されうる撮像素子群の情報からあらかじめ見積もっておくことができる。

一般的な例として、撮像素子のセンサピッチ１．５μｍ〜７μｍ程度の範囲で使用する場合、ナイキスト周波数Ｕｎｙｑは７１〜３３３本／ｍｍとなる。ナイキスト周波数Ｖｎｙｑも同様である。例えば、フィルタタップ数Ｕを６４タップとすると、実空間サイズはＵ／（２＊Ｕｎｙｑ）であるから６４／（２＊Ｕｎｙｑ）となる。したがって、センサピッチ１．５μｍの撮像素子においては６４タップで９６μｍを表すことができる。また、センサピッチ７μｍの撮像素子においては６４タップで４４８μｍを表すことができる。撮像光学系によるＰＳＦの主要な広がりが１５０ｕｍ程度の場合、センサピッチ７μｍの撮像素子では問題なくフィルタ領域に納めることができるが、センサピッチ１．５μｍの撮像素子では領域が不足する。一方、ＰＳＦの主要な広がりが６０ｕｍ程度の場合、センサピッチ１．５μｍの撮像素子にでは４０タップ程度の使用で適切であるが、センサピッチ７μｍの場合６タップ程度しか有効に使用されておらず、メモリや計算コスト等の処理上の無駄がある。

そこで、様々な撮像光学系と撮像素子の組み合わせから、最もＰＳＦが大きくなる撮像光学系に対して、最も狭いセンサピッチを持つ撮像素子においてＰＳＦを収めることができる十分な第１のタップ数をあらかじめ算出しておく。

ステップＳ１０３では、上述のような手順で修正された光学伝達関数に対して、それぞれ周波数変換を行って実空間上のＰＳＦデータに変換することで、撮像素子の特性を含む第１のＰＳＦが生成される。

ステップＳ１０４では、第１のＰＳＦに基づいて、色成分と像高ごとに分布指標値が取得される。分布指標値とは、第１のＰＳＦの第１のタップ数における空間的な広がりを表す指標値であり、本実施例では、色成分、像高ごとにおける第１のＰＳＦの値の最大値とする。また、第１のＰＳＦの値の合計値は、色成分ごとに１に正規化されている。本実施例では、分布指標値は、第１のＰＳＦの値の最大値に基づいて取得されたが、所定の閾値より大きい第１のＰＳＦの値に基づいて取得してもよい。また、本実施例では、分布指標値は演算処理により算出されるが、あらかじめ算出しておき、係数テーブルのような形式で保持しておいてもよい。

ステップＳ１０５では、第１のタップ数以下になるように、色成分と像高ごとの分布指標値に基づいて、第２のタップ数が算出される。分布指標値、第１のタップ数、および第２のタップ数をそれぞれＤＶａｌ１、Ｔ１、Ｔ２、像高と色成分をそれぞれＨｅｉｇｈｔ、Ｃｏｌｏｒとすると、分布指標値ＤＶａｌ１、第２のタップ数Ｔ２はそれぞれ、式（９），（１０）で表される。なお、Ｍａｘ＿ｘ（ｆ（ｘ））は、変数ｘに関するｆ（ｘ）の最大値である。
ＤＶａｌ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ）＝Ｍａｘ＿ｘｙ（ＰＳＦ（ｘ，ｙ）） （９）
Ｔ２（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ）＝Ｔ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ）＊（１．０−ＤＶａｌ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ）） （１０）
本実施例では、色成分と像高ごとに第２のタップ数が算出されるが、像高は共通化してもよい。例えば、式（１１）や式（１２）を用いて、第２のタップ数を算出してもよい。なお、Ｍｉｎ＿ｘ（ｆ（ｘ））、Ａｖｅ＿ｘ（ｆ（ｘ））はそれぞれ、変数ｘに関するｆ（ｘ）の最小値、平均値である。

Ｔ２（Ｃｏｌｏｒ）＝Ｍｉｎ＿Ｈｅｉｇｈｔ（Ｔ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ）＊（１．０−ＤＶａｌ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ））） （１１）
Ｔ２（Ｃｏｌｏｒ）＝Ａｖｅ＿Ｈｅｉｇｈｔ（Ｔ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ）＊（１．０−ＤＶａｌ１（Ｈｅｉｇｈｔ，Ｃｏｌｏｒ））） （１２）
図２は式（９）から算出される特定の色成分（Ｒ）の領域ごとの分布指標値、図３は図２の破線で囲まれる分布指標値の領域の第２のタップ数に対応する第２のＰＳＦを示す図である。図３（ａ）−図３（ｃ）はそれぞれ、式（１０），（１１），（１２）から算出される。図３（ａ）−図３（ｃ）の実線で囲まれる領域は、第２のタップ数に対応する代表的な第２のＰＳＦを示しており、全てのＰＳＦを示してはいない。

ステップＳ１０６では、第２のＰＳＦが、第１のＰＳＦの中央部から、第２のタップ数分、値の合計値が１になるように正規化して取得される。

ステップＳ１０７では、第２のＰＳＦを用いて画像復元処理が実行され、復元画像が生成される。画像復元処理では、第２のＰＳＦをそのまま用いて処理を行ってもよいし、周波数空間に変換して処理を行ってもよい。

本実施例では、図４を参照して、ＰＳＦの特性が画像中の位置によって大きく変化する場合でも、良好に画像復元処理を実行可能な画像処理手法について説明する。図４は、本実施例の撮像素子の特性を考慮したＰＳＦを用いた画像復元処理の説明図である。

ＰＳＦが位置ごとに大きく変化する場合、位置ごとにそれぞれの正確なＰＳＦを用いて処理することが理想的である。しかしながら、全画素数分の正確なＰＳＦをデータとしてメモリに保持することは現実的ではない。そこで、本実施例では、あらかじめ設定された位置において、その各々の位置に対応する撮像光学系の光学伝達関数を撮像素子の特性を考慮して修正し、さらに周波数変換により第１のＰＳＦに変換して用いる。そして、第１のＰＳＦが第１のタップ数に基づいて生成され、第１のＰＳＦに基づいて分布指標値が算出される。その後、分布指標値から第２のタップ数が算出され、第２のタップ数に基づいて第２のＰＳＦが生成される。これらの処理は、実施例１と同様であるので、詳細な説明は省略する。以下の説明では、簡単のために第２のＰＳＦをＰＳＦと記載する。

図４（ａ）に示されるように、画像のＲ，Ｇ，Ｂのそれぞれ縦、横のデータ数をＶ，Ｈとする。Ｒ，Ｇ，Ｂ画素をベイヤーとして処理する場合は、Ｖ，Ｈがそれぞれ異なる数でもよい。画像の縦横にはそれぞれＭ，Ｎ点のＰＳＦを対応付け、それぞれの配置間隔をＡ，Ｂとし（Ｍ−１）×（Ｎ−１）個の複数の長方形領域を考える。ここで、Ａ＝Ｖ／（Ｍ−１），Ｂ＝Ｈ／（Ｎ−１）であることが計算上望ましい。このようにすることで、処理対象の着目画素（ｉ，ｊ）ごとに、ｉ／Ａ，ｊ／Ｂの整数除算をすることで、条件分岐等の処理を用いることなく上記複数の長方形領域の何番目に属するかを直ちに算出できる。本実施例では、ｍ＝ｉ／Ａ，ｎ＝ｊ／Ｂより、図４に示されるように、ＰＳＦ［ｍ，ｎ］，ＰＳＦ［ｍ，ｎ＋１］，ＰＳＦ［ｍ＋１，ｎ］，ＰＳＦ［ｍ＋１，ｎ＋１］の４つのＰＳＦが選択され、この長方形領域を領域＜ｍ,ｎ＞と表す。図４（ｂ）は、図４（ａ）の領域＜ｍ,ｎ＞を抜き出した図である。

次に、図４（ｃ）に示される着目画素（ｉ，ｊ）に対する画像復元処理について説明する。着目画素（ｉ，ｊ）に対応するｒｉｊは、ｒ［ｍ，ｎ］，ｒ［ｍ，ｎ＋１］，ｒ［ｍ＋１，ｎ］，ｒ［ｍ＋１，ｎ＋１］の４つのＰＳＦであり、以下の式（１３）で表される。

ここで、Ｗ（ｍ，ｎ）（＝（Ａ−ａ）＊（Ｂ−ｂ）／（Ａ＊Ｂ））はＰＳＦ［ｍ，ｎ］の重み、Ｗ（ｍ，ｎ＋１）（＝（Ａ−ａ）＊ｂ／（Ａ＊Ｂ））はＰＳＦ［ｍ，ｎ＋１］の重みである。また、Ｗ（ｍ＋１，ｎ）（＝ａ＊（Ｂ−ｂ）／（Ａ＊Ｂ））はＰＳＦ［ｍ＋１，ｎ］の重み、Ｗ（ｍ＋１，ｎ＋１）（＝ａ＊ｂ／（Ａ＊Ｂ））はＰＳＦ［ｍ＋１，ｎ＋１］の重みである。ただし、ａ＝ｍｏｄ（ｉ，Ａ），ｂ＝ｍｏｄ（ｊ，Ｂ）であり、ｍｏｄ（ｘ，ｙ）はｘのｙによる剰余である。

本実施例では、式（１５）に基づいて式（７）が以下の式（１４）で表される。

式（１４）による処理は、それぞれ画像のＲ，Ｇ，Ｂごとに所定の繰り返し回数だけ実行する。この繰り返し回数はあらかじめ、画像復元処理の実行前に設定しておくことが望ましい。特に、本実施例では式（１３）を少ないＭとＮにおいて高精度に行うため、あらかじめそれぞれのＰＳＦの重心位置を一致させる処理を施しておくことが望ましい。また、重心位置だけでなく、最大値を一致させてもよいし、両方を評価する評価値を満たすＰＳＦの補正を行ってもよい。

本実施例では、ＰＳＦの特性が画像中の位置によって大きく変化する場合でも、良好に画像復元処理を実行可能な画像処理手法について説明する。本実施例では、実施例２と異なり、ＰＳＦをそのまま用いるのではなく、回復フィルタを設計して使用する。本実施例では回復フィルタとしてウィナーフィルタに準じたものを使用するが、本発明はウィナーフィルタに限定されない。また、本実施例では、図４（ａ）−図４（ｃ）のＰＳＦを実空間上のウィナーフィルタｒとして置き換えて説明する。

ＰＳＦが位置ごとに大きく変化する場合、位置ごとにそれぞれの正確なＰＳＦを用いて処理することが理想的である。しかしながら、全画素数分の正確なＰＳＦをデータとしてメモリに保持することは現実的ではない。そこで、本実施例では、あらかじめ設定された位置において、その各々の位置に対応する撮像光学系の光学伝達関数を撮像素子の特性を考慮して修正し、さらに周波数変換により第１のＰＳＦに変換して用いる。そして、第１のＰＳＦが第１のタップ数に基づいて生成され、第１のＰＳＦに基づいて分布指標値が算出される。その後、分布指標値から第２のタップ数が算出され、第２のタップ数に基づいて第２のＰＳＦが生成される。これらの処理は、実施例１と同様であるので、詳細な説明は省略する。以下の説明では、簡単のために第２のＰＳＦをＰＳＦと記載する。

本実施例では、ＰＳＦにフーリエ変換を施しＯＴＦを生成する。そして、生成したＯＴＦと式（４）を用いて生成された周波数空間上のウィナーフィルタＲ（ｕ，ｖ）にフーリエ逆変換を施すことで実空間上のウィナーフィルタｒを生成する。

本実施例では、図４（ａ）に示されるように、画像のＲ，Ｇ，Ｂのそれぞれ縦、横のデータ数をＶ，Ｈとする。Ｒ，Ｇ，Ｂをベイヤーとして処理する場合は、Ｖ，Ｈがそれぞれ異なる数でもよい。画像の縦横にはそれぞれＭ，Ｎ点の実空間上のウィナーフィルタｒを対応付け、それぞれの配置間隔をＡ，Ｂとし（Ｍ−１）×（Ｎ−１）個の複数の長方形領域を考える。ここで、Ａ＝Ｖ／（Ｍ−１），Ｂ＝Ｈ／（Ｎ−１）であることが計算上望ましい。このようにすることで、処理対象の着目画素（ｉ，ｊ）ごとに、ｉ／Ａ，ｊ／Ｂの整数除算をすることで、条件分岐等の処理を用いることなく複数の長方形領域の何番目に属するかを直ちに算出できる。本実施例では、ｍ＝ｉ／Ａ，ｎ＝ｊ／Ｂより、図４に示されるように、ｒ［ｍ，ｎ］，ｒ［ｍ，ｎ＋１］，ｒ［ｍ＋１，ｎ］，ｒ［ｍ＋１，ｎ＋１］の４つの実空間上のウィナーフィルタが選択され、この長方形領域を領域＜ｍ,ｎ＞と表す。図４（ｂ）は、図４（ａ）の領域＜ｍ,ｎ＞を抜き出した図である。

次に、図４（ｃ）に示される着目画素（ｉ，ｊ）に対する画像復元処理について説明する。着目画素（ｉ，ｊ）に対応するｒｉｊは、ｒ［ｍ，ｎ］，ｒ［ｍ，ｎ＋１］，ｒ［ｍ＋１，ｎ］，ｒ［ｍ＋１，ｎ＋１］の４つの実空間上のウィナーフィルタｒであり、以下の式（１５）で表される。

ここで、Ｗ（ｍ，ｎ）（＝（Ａ−ａ）＊（Ｂ−ｂ）／（Ａ＊Ｂ））はｒ［ｍ，ｎ］の重み、Ｗ（ｍ，ｎ＋１）（＝（Ａ−ａ）＊ｂ／（Ａ＊Ｂ））はｒ［ｍ，ｎ＋１］の重みである。また、Ｗ（ｍ＋１，ｎ）（＝ａ＊（Ｂ−ｂ）／（Ａ＊Ｂ））はｒ［ｍ＋１，ｎ］の重み、Ｗ（ｍ＋１，ｎ＋１）（＝ａ＊ｂ／（Ａ＊Ｂ））はｒ［ｍ＋１，ｎ＋１］の重みである。ただし、ａ＝ｍｏｄ（ｉ，Ａ），ｂ＝ｍｏｄ（ｊ，Ｂ）であり、ｍｏｄ（ｘ，ｙ）はｘのｙによる剰余である。

本実施例では、式（１５）に基づいて式（６）が以下の式（１６）で表される。

式（１６）による処理は、それぞれ画像のＲ，Ｇ，Ｂごとに実行することが望ましい。特に本実施例では式（１５）を少ないＭとＮにおいて高精度に行うため、あらかじめそれぞれのＰＳＦの重心位置を一致させる処理を施しておくことが望ましい。また、重心位置だけでなく、最大値を一致させてもよいし、両方を評価する評価値を満たすＰＳＦの補正を行ってもよい。

本実施例では、本発明の画像処理方法の各工程をコンピュータに実行させる画像処理プログラムを用いて画像復元処理を実行する方法について、図５を参照して説明する。図５は、本発明の画像処理プログラムをインストールして実行する装置の一例であるＰＣ本体５０の外観図である。

ＰＣ本体５０には、本発明の画像処理プログラムがインストールされる。ＰＣ本体上でインストールが完了しているソフトを起動することで、本発明の画像処理方法をソフト上で実行させることができる。なお、ＰＣ本体５０は、画像処理プログラムが記録されている記録媒体を用いて、本発明の画像処理方法を実行してもよい。また、ＰＣ本体５０の所定のフォルダには、撮像光学系の光学伝達関数が格納されている。ここでの光学伝達関数は、撮像光学系の様々な撮影条件ごとの全ての光学伝達関数のことである。画像処理プログラムおよび光学伝達関数は、可搬メディアからインストールしてもよいし、インターネットを通じてダウンロードしてもよい。

ソフト上では、第２のタップ数を決定する方法を、ユーザーが式（１０），（１１），（１２）の方式から選択できるようにしてもよい。また、画像復元処理が繰り返し演算を必要とする場合は、処理の強さを決める繰り返し回数や、画像の所定の対象領域にのみ処理をかけるなど様々な設定を行えるようにしてもよい。ユーザーは、モニター５１上で処理後の画像を確認しつつ、上記パラメータが変更できることがより好ましい。処理の内容については実施例１〜３に詳述したので、説明は省略する。

本発明の画像処理方法を実行する画像処理装置について、図６を参照して説明する。図６（ａ）は本実施例の画像処理装置６０Ａの構成図、図６（ｂ）は画像処理装置６０Ｂの構成図である。まず、図６（ａ）を用いて画像処理装置６０Ａについて説明する。

撮影画像取得部（第１の取得手段）６１は、図６撮影条件に関する情報が付加された撮影画像を取得し、取得した撮影画像を撮影条件取得部６２と画像復元部６９に送信する。

撮影条件取得部（第２の取得手段）６２は、撮像光学系に関する情報と、撮像素子に関する情報を撮影条件として取得する。撮像光学系に関する情報とは、例えば、焦点距離、Ｆ値、被写体距離であり、その他にも撮像光学系のイメージサークル径等である。撮像光学系のイメージサークルとは、撮像光学系の設計時に、収差補正上あらかじめ定めた単位円の半径または直径である。撮像素子に関する情報とは、例えば、撮像素子のサイズ、素子のピッチ、画素の形状や画素のサイズである。また、撮像素子側に光学ローパスフィルタ等の種々の素子が装着されている場合があるが、そのような場合には種々の素子の特性に関する情報も撮像素子に関する情報に含めてよい。

光学伝達関数記憶部６３はあらかじめ撮像光学系の光学伝達関数を記憶し、光学伝達関数取得部６４は撮影条件取得部６２で取得した撮像光学系の情報とあらかじめ設定されている第１のタップ数に基づいて複数の光学伝達関数を取得する。光学伝達関数修正部６５は、撮像素子の情報に基づいて、光学伝達関数取得部６４が取得する光学伝達関数に対して修正を加える。ここで、修正とは、光学ローパスフィルタや画素の特性を付与することや、撮像素子のピッチに合わせた周波数帯域までの光学伝達関数に対する変換処理である。

第１点像強度分布関数生成部（第１の生成手段）６６は、光学伝達関数修正部６５が修正した光学伝達関数を全てフーリエ逆変換し、撮像光学系と撮像素子の特性に基づく第１の点像強度分布関数を生成する。分布指標値取得部（第３の取得手段）６７は、第１点像強度分布関数６６に基づいて、各色成分、像高ごとに分布指標値を取得する。

第２点像強度分布関数生成部（第２の生成手段）６８は、分布指標値取得部６７により取得される分布指標値に基づいて第２のタップ数を算出し、算出した第２のタップ数と第１の点像強度分布関数に基づいて第２の点像強度分布関数を生成する。

画像復元部（復元手段）６９は、第２の点像強度分布関数を用いて、撮影画像取得部６１が取得する撮影画像に対して画像復元処理を実行する。

画像処理装置６０Ｂでは、撮影画像取得部６１は撮影条件を付加されていない撮影画像を取得するため、撮影画像を画像復元部６９にだけ送信し、撮影条件取得部６２は撮像装置等から直接撮影条件を取得する構成が画像処理装置６０Ａと異なる。他の構成については画像処理装置６０Ａと同様であるので、詳細な説明は省略する。

本発明の画像処理方法を実行する画像処理装置を内蔵する撮像装置について、図７を参照して説明する。図７は、本実施例の撮像装置７０の構成を示すブロック図である。

撮像装置７０は、撮像光学系７０１、撮像素子７０２、Ａ／Ｄ変換部７０３、画像処理部（画像処理装置）７０４、表示部７０５、撮像光学系制御部７０６、状態検出部７０７、記憶部７０８、画像記録媒体７０９およびシステムコントローラ７１０を備える。

撮像光学系７０１は、絞り７０１ａとフォーカスレンズ７０１ｂを備え、被写体像を取得する。撮像素子７０２は、撮像光学系７０１により取得される被写体像を光電変換して撮影画像を生成する。このとき、撮像光学系制御部７０６や状態検出部７０７が、例えば、撮影時の焦点距離、Ｆ値、撮影距離情報といった撮影条件を取得する。記憶部７０８は、撮像光学系の光学伝達関数をデータ化し、撮影条件（例えば、焦点距離、Ｆ値、撮影距離、像高の組み合わせ）ごとに保持する。

画像処理部７０４は、取得した撮影条件とあらかじめ設定された第１のタップ数に基づいて撮像光学系の光学伝達関数を取得し、撮像素子７０２の情報を用いて撮像光学系の光学伝達関数を修正する。画像処理部７０４の第１点像強度分布関数生成部は、修正された光学伝達関数を全てフーリエ逆変換し、撮像光学系と撮像素子の特性に基づく第１の点像強度分布関数を生成する。画像処理部７０４の分布指標値取得部は、第１の点像強度分布関数に基づいて、各色成分、像高ごとに分布指標値を取得する。画像処理部７０４の第２点像強度分布関数生成部は、取得された分布指標値に基づいて第２のタップ数を算出し、算出した第２のタップ数と第１の点像強度分布関数に基づいて第２の点像強度分布関数を生成する。画像処理部７０４の画像復元部は、記憶部７０８から画像データを読み出し、第２の点像強度分布関数を用いて読み出した画像データに対して画像復元処理を実行する。

本実施例の画像復元処理は、画像の所定の対象領域にのみに実行するなど様々な設定を行えるようにしてもよい。ユーザーは、表示部７０５上で処理後の画像を確認しつつ、各種のパラメータが変更できることがより好ましい。

以上、本発明の好ましい実施形態について説明したが、本発明はこれらの実施形態に限定されず、その要旨の範囲内で種々の変形及び変更が可能である。

６０Ａ，６０Ｂ 画像処理装置
６１ 撮影画像取得部（第１の取得手段）
６２ 撮影条件取得部（第２の取得手段）
６６ 第１点像強度分布関数生成部（第１の生成手段）
６７ 分布指標値取得部（第３の取得手段）
６８ 第２点像強度分布関数生成部（第２の生成手段）
６９ 画像復元部（復元手段）

Claims (13)

1. 撮影画像を取得する第１の取得手段と、
   撮像光学系および撮像素子の情報を取得する第２の取得手段と、
   前記撮像素子の画素ピッチと第１のタップ数に基づいて、前記撮像光学系に関する第１の点像強度分布関数を生成する第１の生成手段と、
   前記第１の点像強度分布関数に関する分布指標値を取得する第３の取得手段と、
   前記第１点像強度分布関数と、前記分布指標値に基づいて算出される前記第１のタップ数以下の第２のタップ数に基づいて第２の点像強度分布関数を生成する第２の生成手段と、
   前記第２の点像強度分布関数に基づいて前記撮影画像を復元する復元手段と、を有することを特徴とする画像処理装置。
2. 前記撮像光学系の情報は、前記撮像光学系の光学伝達関数を含むことを特徴とする請求項１に記載の画像処理装置。
3. 前記分布指標値は、前記第１の点像強度分布関数の前記第１のタップ数における空間的な広がりを表す指標値であることを特徴とする請求項１または２に記載の画像処理装置。
4. 前記分布指標値は、前記第１の点像強度分布関数の前記第１のタップ数における最大値に基づいて算出されることを特徴とする請求項３に記載の画像処理装置。
5. 前記第２のタップ数は、前記撮影画像の領域ごとに取得される分布指標値に基づいて算出されることを特徴とする請求項１から４のいずれか１項に記載の画像処理装置。
6. 前記第２のタップ数は、前記撮影画像の領域ごとに取得される分布指標値の平均値に基づいて算出されることを特徴とする請求項１から４のいずれか１項に記載の画像処理装置。
7. 前記第２のタップ数は、前記撮影画像の領域ごとに取得される分布指標値の最小値に基づいて算出されることを特徴とする請求項１から４のいずれか１項に記載の画像処理装置。
8. 前記撮像光学系の情報は、前記撮像光学系の焦点距離、Ｆ値、および被写体距離の少なくとも１つを特定可能な情報を含むことを特徴とする請求項１から７のいずれか１項に記載の画像処理装置。
9. 前記撮像光学系および前記撮像素子の情報は、前記撮像光学系および前記撮像素子を特定可能な情報を含むことを特徴とする請求項１から８のいずれか１項に記載の画像処理装置。
10. 請求項１から９のいずれか１項に記載の画像処理装置と、
    被写体像を光電変換して撮影画像を取得する撮像素子と、を有することを特徴とする撮像装置。
11. 撮影画像を取得するステップと、
    撮像光学系および撮像素子の情報を取得するステップと、
    前記撮像素子の画素ピッチと第１のタップ数に基づいて、前記撮像光学系に関する第１の点像強度分布関数を生成するステップと、
    前記第１の点像強度分布関数に関する分布指標値を取得するステップと、
    前記分布指標値に基づいて前記第１のタップ数以下の第２のタップ数を算出するステップと、
    前記第１点像強度分布関数と前記第２のタップ数に基づいて前記第２の点像強度分布関数を生成するステップと、
    前記第２の点像強度分布関数に基づいて前記撮影画像を復元するステップと、を有することを特徴とする画像処理方法。
12. 請求項１１に記載の画像処理方法の各工程をコンピュータに実行させるための画像処理プログラム。
13. 請求項１２に記載の画像処理プログラムを記録するコンピュータで読み取り可能な記録媒体。