JP4724826B2

JP4724826B2 - 超解像処理のパラメータ設計方法

Info

Publication number: JP4724826B2
Application number: JP2005008593A
Authority: JP
Inventors: 正行田中; 正敏奥富
Original assignee: Tokyo Institute of Technology NUC
Current assignee: Tokyo Institute of Technology NUC
Priority date: 2005-01-17
Filing date: 2005-01-17
Publication date: 2011-07-13
Anticipated expiration: 2025-01-17
Also published as: JP2006195856A

Description

本発明は、複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法に関し、特に、再構成型超解像処理において、高解像度画像を安定に推定できるパラメータを設計できる超解像処理のパラメータ設計方法に関する。

超解像処理とは、位置ずれのある複数の低解像度画像（観測画像）より一つの高解像度画像を推定する処理である。高解像度画像から低解像度画像への画像変換は、点広がり関数(ＰＳＦ)と呼ばれる積分核との積分変換として、定式化される。超解像処理は、この積分変換の逆変換に相当する処理である。

これまでにも、例えば、非特許文献１〜非特許文献６に開示されるように、数多くの超解像処理方法が提案されている。その中でも、再構成型超解像処理と分類される方法が広く利用されている。例えば、非特許文献２に開示されているＭＬ（maximum likelihood）法、非特許文献３に開示されているＭＡＰ（maximum ａ posteriori）法や非特許文献４に開示されているＰＯＣＳ（projection onto convex sets）法が、再構成型超解像処理方法として有名である。

また、多くの超解像処理に関する研究は、グレイスケール画像を対象としているが、ＭＡＰ法を応用してBayer配列の低解像度画像からカラーデモザイキング処理と超解像処理を同時に行う方法も提案されている（非特許文献５及び非特許文献６参照）。さらに、再構成型超解像処理は、計算コストが大きいことも知られており、効率的な計算方法に関する報告もされている（非特許文献７及び非特許文献８参照）。
サングシー．ピー．（Sung C.P.）・ミンケイ．ピー．（Min K.P.）共著、「スーパーレゾルーションイメージリコンストラクション:アーテクニカルオーバービュー（Super-ResolutionImage Reconstruction: A Technical Overview）」,ＩＥＥＥシグナルプロク．マガジン（IEEE SignalProc. Magazine）,第26巻,第3号, p.21-36，2003年ビー．シー．トム（B.C.Tom）・エイ．ケイ．カトサゲロス（A.K.Katsaggelos）共著、「リコンストラクションオフアーハイレゾルーションイメージバイシマルテイニアスレジストレイションレストレイションアンドインターポレイションオフローレゾルーションイメージズ（Reconstructionof a high-resolution image by simultaneous registration, restoration, andinterpolation of low-resolution images）」,プロク．ＩＥＥＥイント．コンフ．イメージプロセシング（Proc.IEEE Int. Conf. Image Processing）,第2巻,p.539-542,1995年アール.アール.シュルツ（R.R.Schulz）・アール.エル.スティブンスン（R.L.Stevenson）共著、「エクストラクションオフハイレゾルーションフレームズフロムビデオシーケンス（Extractionof high-resolution frames from video sequences）」,ＩＥＥＥトランス. イメージプロセシング（IEEETrans. Image Processing）,第5巻,p.996-1011,1996年エイチ.スターク（H.Stark）・ピー.オスコウイ（P.Oskoui）共著、「ハイレゾルーションイメージリカバリーフロムイメージプレーンアレーズユーシングコンベックスプロジェクションズ（Highresolution image recovery from image-plane arrays, using convex projections）」,ジェー. オプト. ソク. エム. エイ（J.Opt. Soc. Am. A）,第6巻,p.1715-1726,1989年後藤・奥富共著、「単板カラー撮像素子のＲＡＷデータを利用した高精細画像復元」、情報処理学会論文誌：コンピュータビジョンとイメージメディア,第45巻,第SIG8(CVIM９)号,p.15-25,2004年ティー.ゴトウ（T.Goto）・エム.オクトミ（M.Okutomi）共著、「ダイレクトスーパーレゾリューションアンドレジストレーションユーシングローＣＦＡイメージズ（DirectSuper-Resolution and Registration Using Raw CFA Images）」,プロク．オフＩＥＥＥコンピュータソサイエティカンフェレンスオンＣＶＰＲ（Proc.of IEEE Computer Society Conference on CVPR）,第II巻,p.600-607,2004年田中・奥富共著、「再構成型超解像処理の高速化アルゴリズム」、情報処理学会研究報告2004-CVIM-146,第2004巻,第113号,p.97-104,2004年エヌ.グエン（N.Nguyen）・ピー.ミランファー（P.Milanfar）共著、「アコンピューティションナリエフィシェントスーパーレゾリューションイメージレコンストラクションアルゴリズム（AComputationally Efficient Super-Resolution Image Reconstruction Algorithm）」,ＩＥＥＥトランス. イメージプロセシング（IEEETrans. Image Processing）,第10巻,第4号,p.573-583,2001年伊里正夫著、「線形代数Ｉ,II」、岩波書店、1994年エム.エラッド（M.Elad）・エイ.フォイヤ（A.Feuer）共著、「レストレイションオフアシングルスーパーレゾリューションイメージフロムセブラルブラーレドノイジーアンドアンダーサンプレドメジャードイメージズ（Restorationof a Single Super-Resolution Image from Several Blurred, Noisy, andUndersampled Measured Images）」,ＩＥＥＥトランス. イメージプロセシング（IEEE Trans. Image Processing）,第6巻,第12号,p.1646-1658,1997年エス.ベイカ（S.Baker）・ティー.カナデ（T.Kanade）共著、「リミッズオンスーパーレゾリューションアンドハウトゥーブレイクゼム（Limitson Super-Resolution and How to Break Them）」,ＩＥＥＥトランス. パターンアナリシスアンドマシンインテリジェンス（IEEETrans. Pattern Analysis and Machine intelligence）,第24巻,第9号,p.1167-1183,2002年ゼット.リン（Z.Lin）・エイチ.ワイ.シュム（H.Y.Shum）共著、「ファンダメンタルリミッズオフレコンストラクションベーセドスーパーレゾリューションアルゴリズムズアンダーローカルトランスレイション（Fundamentallimits of Reconstruction-Based Super-Resolution Algorithms under localTranslation）」,ＩＥＥＥトランス. パターンアナリシスアンドマシンインテリジェンス（IEEE Trans. PatternAnalysis and Machine intelligence）,第26巻,第1号,p.83-97,2004年ディー.ケープル（D.Capel）著、「イメージモザイキングアンドスーパーレゾリューション（Image Mosaicing andSuper-resolution）」,スプリンガー出版（Springer Verlag）,2004年ディー.ロビンスン（D.Robinson）・ピー.ミランファー（P.Milanfar）共著、「ファンダメンタルパフォーマンスリミッズインイメージレジストレーション（FundamentalPerformance Limits in Image Registration）」,ＩＥＥＥトランス. イメージプロセシング（IEEE Trans.Image Processing）,第13巻,第9号,p.1185-1199,2004年清水・矢野・奥富共著、「２次元サブピクセル同時推定を拡張した画像変形Ｎパラメータ同時推定」、情報処理学会研究報告,2004-CVIM-143,第2004巻,第26号,p.81-88,2004年

ところで、実際に、再構成型超解像処理を行うためには、高解像度画像の低解像度画像に対する解像度の倍率Ｍと、ＰＳＦを事前に設定する必要がある。超解像処理のパラメータである倍率ＭとＰＳＦは、重要なパラメータであり、これらパラメータの設定を誤まると、超解像処理の出力が正しく求まらなくて良い超解像処理結果が得られないことが知られている。

また、超解像処理において、低解像度画像の枚数を増加させることによって、大きな倍率Ｍを設定することが可能であると予想されるが、低解像度画像を増加させても倍率の大きさには限界があることが経験的に知られている。この倍率の限界は、再構成型超解像処理に関しての基本的な性質である。また、低解像度画像の枚数を無限に増加させた場合に、この性質は顕著に現われることが予想される。

しかし、超解像処理において、高解像度画像の低解像度画像に対する倍率の限界が再構成型超解像処理の基本的な性質であるにもかかわらず、この分野の研究は非常に少ない。数少ない研究の例を挙げると、例えば、非特許文献１０に記載されるように、エラッド(Elad)らは、未知数の数と方程式の数の関係から、倍率と観測画像の枚数（低解像度画像の枚数）が満足しなければ成らない関係式を示している。しかし、その関係式は、単純に観測画像の総画素数（方程式の数）は、高解像度画像の画素数（未知数の数）より多くなくてはならないことを示しているに過ぎない。

また、非特許文献１１では、Box型ＰＳＦ(Point Spread Function)を仮定して、倍率Ｍが２以上の整数のとき、超解像の方程式の条件数が無限大となり、高解像度画像を安定に推定できないことが開示されている。更に、非特許文献１１では、Box型ＰＳＦの場合は、倍率Ｍを大きくするにつれて、条件数が大きくなると定性的に述べられている。

そして、非特許文献１２に記載されるように、リン（Lin）らは、Box型ＰＳＦを仮定し、摂動定理を応用することにより、条件数の逆数に対応する指標を定量的に評価している。このように、従来、再構成型超解像処理において、高解像度画像の低解像度画像に対する倍率の限界に関する研究は、Box型ＰＳＦが用いられた場合のみである。

一方、非特許文献１３に開示されるように、ケープル（Capel）は、カメラモデルのＰＳＦはGaussian型ＰＳＦで近似できることを実験的に示している。ところで、Box型ＰＳＦも、Gaussian型ＰＳＦも、どちらもカメラモデルのＰＳＦの近似でしかない。そのため、どの近似を利用すべきかといった問題が重要である。

しかしながら、上記のように、従来では、単純なBox型ＰＳＦを仮定した場合の倍率Ｍの限界に関する研究のみがなされ、より一般的なＰＳＦを仮定した場合の倍率の限界に関する研究はなされていないため、超解像処理に必要なパラメータであるＰＳＦを適切に設計する方法は、従来になかった。また、従来では、超解像処理のもう１つのパラメータである倍率の設計に関しても、単純なBox型ＰＳＦを仮定した場合のみが検討された。

以上のように、従来では、再構成型超解像処理に必要なパラメータ（つまり、カメラモデルから得られるＰＳＦ関数、高解像度画像の低解像度画像に対する倍率Ｍ）を最適に設計する超解像処理のパラメータ設計方法は、研究されていない。

本発明は、上述のような事情よりなされたものであり、本発明の目的は、再構成型超解像処理において、高解像度画像を安定に推定できるパラメータを設計できる超解像処理のパラメータ設計方法を提供することにある。

本発明は、位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法に関し、本発明の上記目的は、前記パラメータは、前記高解像度画像の前記低解像度画像に対する倍率、及びカメラモデルから得られるＰＳＦ関数であり、前記超解像処理の基本方程式の解を求める際の安定性を表す、前記基本方程式の条件数に基いて、前記高解像度画像を安定に推定できる最適化パラメータを決定し、前記基本方程式は、次の式で表される線形連立方程式であり、

ただし、

は前記複数の低解像度画像の全ての画素値を表すベクトルで、

は推定すべき前記高解像度画像の画素値を表すベクトルで、

は前記ＰＳＦ関数及び前記倍率から導かれるシステム行列を表し、（Ａ１）前記位置ずれは平行移動であり、（Ａ２）前記ＰＳＦ関数は空間不変であり、（Ａ３）前記高解像度画像の画素値は同一画素内で一定であり、（Ａ４）前記低解像度画像は位置の重なりがなく、無限に取得できるといった４つの条件を前提に前記条件数を算出することによって効果的に達成される。

また、本発明の上記目的は、前記条件数cond(Ｂ)は、次の式により算出され、

ただし、

が成立し、また、

が成立し、

は

の離散フーリエ変換であり、ｂ(ｘ,ｙ)は高解像度画像空間における連続のＰＳＦ関数を表すことにより、或いは、所定のＰＳＦ関数に対して、異なる倍率における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応する倍率を、前記超解像処理に前記所定のＰＳＦ関数が用いられた場合の最適な倍率にすることにより、或いは、所定の倍率に対して、異なるＰＳＦ関数における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応するＰＳＦ関数を、前記超解像処理に前記所定の倍率が指定された場合の最適化ＰＳＦ関数にすることによって一層効果的に達成される。

本発明では、線形連立方程式の解を求める際の安定性を表す定量的指標として利用されている条件数に注目し、任意の空間不変なＰＳＦ、任意の倍率に対して、超解像の基本方程式の条件数を算出する方法である条件数算出方法を開示し、そして、その条件数算出方法を証明した。

本発明に係る超解像処理のパラメータ設計方法によれば、任意の空間不変なＰＳＦ、任意の倍率に対して、超解像の基本方程式の条件数を算出することが可能であり、そして、それぞれの条件数を比較することにより、より良い超解像結果を得るためには、所望の倍率ではどのＰＳＦ近似が良いか、また、指定されたＰＳＦ近似では倍率をどの値に設定すれば良いかを判断することができ、よって、高解像度画像を安定に推定できる最適なパラメータを設計できるという優れた効果を奏する。

また、合成画像に対して、本発明に係る超解像処理のパラメータ方法を適用することによって、高解像度画像を安定に推定できる最適化ＰＳＦ関数、最適な倍率を設定できることを確認した。

以下、本発明を実施するための最良の形態を図面を参照して説明する。なお、本発明では、『復元したい高解像度画像の空間』を『高解像度画像空間』と称する。

まず、超解像の基本方程式について説明する。

再構成型超解像処理の基本となる連続系の画像変換方程式を下記数１に示す。

ここで、添え字ｊはｊ番目の観測を、ｙ_ｊは観測画像の座標を、ｘは高解像度画像の座標を、ｇ_ｊ(ｙ_ｊ)は連続観測画像を、ｂ_ｊ(ｙ_ｊ,ｘ)は連続のＰＳＦを、ｈ(ｘ)は連続高解像度画像を、それぞれ表す。

実際に、超解像処理が扱う方程式は、上記数１を離散化して得られた下記数２である。以下、この数２を超解像の方程式と呼ぶことにする。

ここで、

は複数の低解像度画像の全ての画素値を表すベクトルを、

は推定すべき高解像度画像の画素値を表すベクトルを、

は高解像度画像から低解像度画像への変換を表す行列（つまり、ＰＳＦ、倍率Ｍ及び低解像度画像画素のレジストレーション後の画素位置から導出されるシステム行列）を、それぞれ表す。

数２で表す超解像の方程式は、低解像度画像が無限にある場合にも拡張することができる。ここで、低解像度画像が無限の場合（低解像度画像の枚数が無限の場合）の超解像の方程式を、超解像の基本方程式と呼ぶことにする。超解像の基本方程式は、上記数２と同様であるが、

の次元と

の行の次元が無限大である。

次に、本発明の理論解析について説明する。

本発明の理論解析は以下の４つの仮定に基づく。

つまり、まず、「仮定１」として、位置ずれは平行移動である。また、「仮定２」として、カメラモデルから得られるＰＳＦ関数（以下、単に、カメラＰＳＦ、もしくは、ＰＳＦとも称する）は、空間不変である。そして、「仮定３」として、連続高解像度画像の画素値は同一画素内で一定で、つまり、高解像度画像の画素内の濃度値は一定である。最後に、「仮定４」として、低解像度画像は、位置の重なりがなく、且つ、無限に取得することができる。

本発明の大きな特徴の１つは、本発明で用いられる「仮定２」（つまり、カメラＰＳＦを単に空間不変と仮定している点）が、これまでの従来研究（非特許文献１１及び非特許文献１２参照）と大きく異なる点である。

ここで、一般性を失わず、高解像度画像の画素間隔を１として解析を行うことができる。ただし、高解像度画像の画素間隔を１と考えると、ＰＳＦは高解像度画像の低解像度画像に対する倍率Ｍをパラメータとして含むようになる。この場合の具体的なＰＳＦの例として、Box型ＰＳＦを下記数３に、Gaussian型ＰＳＦを下記数４に、それぞれ示す。

ここで、

は下記数５で定義される矩形関数である。

以下では、簡単のため、再構成型超解像処理を１次元の問題として（つまり、超解像処理が扱う低解像度画像は１次元画像であるを前提に）、本発明の説明を行うが、本発明はそれに限られることがなく、２次元の問題（つまり、超解像処理が扱う低解像度画像は２次元画像であるの場合）へ簡単に拡張できることは言うまでもない。

本発明では、超解像の基本方程式を以下のように導出する。

上述した４つの仮定（「仮定１」〜「仮定４」）を利用して、低解像度画像を無限に取得できると仮定した場合の超解像の基本方程式の導出を行う。「仮定１」と「仮定２」を利用すれば、数１で表す画像変換方程式から下記数６を得ることができる。

ここで，ｙはレジストレーション後の低解像度画像の座標である。「仮定３」を満足する連続高解像度画像は、数５で表す矩形関数を利用して、下記数７のように定式化することができる。

ここで、Ｎは１次元画像の大きさを、ｉは１次元画像の離散座標を、ｈ[ｉ]は離散高解像度画像の位置ｉにおける画像値を、それぞれ表す。数６及び数７より、数８に示される離散高解像度画像から連続低解像度画像への画像変換方程式を得ることができる。

さて、ある座標ｙを下記数９のように、ｙに最も近い整数

と残りの小数部分δに分けて考える。

このとき、小数部分δの変域は−１／２から１／２となる。数９の関係を利用して、数８を下記数１０の離散畳込み演算の形式に簡単化することができる。

ここで、

はオフセットδの離散ＰＳＦであり、下記数１１で定義される。

次に、小数部分δを固定して考えると、

は１つの離散画像を表しているとみなせる。この離散画像は、連続観測画像ｇ(ｙ)を高解像度画像の画素間隔でサンプリングした画像と同一である。ただし、サンプリング点は、高解像度画像の画素中心とδだけずれている。高解像度画像、低解像度画像及びＰＳＦに関して、それぞれ連続と離散の関係を図１に示す。

ここで、離散画像

を「オフセットδの離散低解像度画像」と呼び、

で表すことにする。数１０は、オフセットδの離散低解像度画像に対する超解像の方程式となる。数１０をベクトル表現で表すと、下記数１２となる。

ここで、

は離散畳込み演算のカーネルを表す行列である。

連続低解像度画像は、オフセットδの離散低解像度画像の線形結合として表される。従って、連続低解像度画像に対する超解像の方程式、即ち、超解像の基本方程式も同様に、オフセットδの離散低解像度画像に対する方程式の線形結合として、下記数１３のように表現することができる。

ここで、下記数１４が成立し、また、ｉ≠ｊのときに、δ_ｉ≠δ_ｊである。

以上のように、前述した４つの仮定から、数１３で表す超解像の基本方程式を導出することができた。

次に、上述したように導出された超解像の基本方程式に対する条件数の算出方法（以下、単に「超解像の条件数算出方法」とも称する）について説明する。

ここで、まず、条件数に関する数学的な性質を概説した後に、本発明のポイントである超解像の条件数算出方法を示し、証明する。

条件数は、線形連立方程式の係数行列に対して定義され、連立方程式の解を求める際の安定性の指標として知られている（非特許文献９参照）。条件数の小さい連立方程式は、解が安定に推定される（推定誤差が相対的に小さい）。また、条件数の大きい連立方程式は、解が不安定である（推定誤差が相対的に大きい）。特に、条件数が無限大のときは、その連立方程式は特異で、つまり、解が一意に求められないことを意味する。

ところで、数２で示されるような連立方程式（つまり、超解像の方程式）に対して、その条件数

、つまり、超解像の方程式の条件数（以下、単に条件数とも称する）は、下記数１５のように定義することができる。

ここで、‖・‖_２はＬ２ノルムを表す。条件数に関して、下記数１６及び数１７で表す不等式が成立し、条件数は解の相対誤差の上限を定める。

ここで、Δｈは推定された高解像度画像の誤差を、Δｇは観測された低解像度画像の誤差を、ΔＢはＰＳＦによって定まる係数行列の誤差を、それぞれ表す。

ちなみに、数１７の不等式は、ＰＳＦの近似精度が多少悪くても条件数が小さければ、結果として誤差の影響は小さく抑えられることを示している。

ところで、行列のノルムを直接計算することは一般に困難であるので、条件数を上記数１５の定義に従って計算することは現実的ではない。

そこで、本発明では、下記数１８に示される条件数の性質を利用する（非特許文献９参照）。

ここで、^＊は行列の共役転置を、

は行列

の固有値をベクトル化したものを、それぞれ表す。

次に、条件数算出方法について詳細に説明する。以下、補題を示した後に、条件数算出方法を証明する。
＜補題＞数１３の係数行列

に対して、下記数１９が成り立つ。

ここで、

は

の離散フーリエ変換を表し、｜・｜^２は要素ごとに適用されるとする。
＜補題の証明＞
先ず、オフセットδの離散低解像度画像に対する超解像の方程式の係数行列

は、下記数２０のように対角化される。

ここで、

は離散フーリエ変換の行列表現に対応するユニタリ回転行列を、

は

の要素を対角成分にもつ対角行列を、それぞれ表す。

次に、行列

の固有値を求めるため、対角化を行う。行列

の定義式と数２０の関係式より、行列

の対角化を下記数２１のように行うことができる。

ここで、行列

はユニタリ行列であるので、数２１の右辺の対角成分は、行列

の固有値となる。従って、数１９が得られる。
＜超解像の条件数算出方法＞
数１３で表す超解像の基本方程式に対する条件数

は、下記数２２により算出される。

ここで、下記数２３及び数２４が成立する。また、

をＰＳＦの平均パワースペクトラムと呼ぶことにする。

以上では、超解像処理が扱う低解像度画像が１次元画像の場合（つまり、超解像処理を１次元の問題として考えた場合）に、超解像の条件数算出方法を詳細に説明した。

なお、超解像処理を２次元の問題として考えた場合に、要するに、超解像処理が扱う低解像度画像が２次元画像である場合に、超解像の基本方程式に対する条件数

は、下記数２５により算出される。

ただし、下記数２６及び数２７が成立する。また、

は

の離散フーリエ変換であり、そして、ｂ(ｘ,ｙ)は高解像度画像空間における連続のＰＳＦ関数を表す。

＜超解像の条件数算出方法の証明＞
ここで、数２２に表される超解像の条件数算出方法を以下のように証明する。

規格化されたカメラＰＳＦを考えた場合に、下記数２８の関係を仮定することができる。

数２８を利用することにより、行列

に対する数１８の分子は、下記数２９となる。

上記数２９の右辺は無限大に発散してしまうが、条件数は数１８に示されるように最小固有値と最大固有値の比であるので、条件数は収束する。

数１８、数１９及び数２９より、下記数３０が得られる。

上記数３０を簡単にすることにより、数２２が得られる。

なお、δ_ｋの変域は−１／２から１／２であり、δ_ｋが重ならず無限あると考えた場合に、数２３を数２４の積分の形式に変形することができる。

超解像の条件数算出方法を利用して、超解像の基本方程式に対する条件数を算出するために、必要なパラメータは、次の３つのパラメータである。

つまり、「パラメータ１」として連続ＰＳＦｂ(ｘ)を、「パラメータ２」として高解像度画像の低解像度画像に対する倍率Ｍを、「パラメータ３」として低解像度画像の大きさをそれぞれ用いる。

「パラメータ１」及び「パラメータ２」が必要な理由は、条件数算出方法より明らかである。「パラメータ３」に関しては、ＰＳＦの平均パワースペクトラムを計算するために、必要である。

つまり、超解像の条件数算出方法は、任意のＰＳＦに対して、任意の倍率での条件数を算出することを可能にする。得られた条件数は、ＰＳＦ近似に対する重要な指標となる。従って、条件数を比較し、条件数のより小さなＰＳＦ近似を利用することにより、高解像度画像を安定に推定することができる。また、同一のＰＳＦ近似で倍率を変化させた場合も、条件数を比較することにより、高解像度画像を安定に推定することができる最適な倍率を設定することができ、更に、倍率の限界を示すことも可能にする。

次に、数１３で表す超解像の基本方程式をＭＬ法で解く場合を考える。この場合に、ＰＳＦの平均パワースペクトラム

がＭＬ法の評価関数の勾配とどのような関係があるかを明らかにする。つまり、ＭＬ法における勾配とＰＳＦのパワースペクトルとの関係を勾配制限として示す。

数１３で表す超解像の基本方程式に関するＭＬ法の評価関数を下記数３１のように定義する。

ここで、評価関数が発散するのを防ぐために、Ｋによる規格化を行っている。評価関数の高解像度画像

に関する微分は、下記数３２となる。

ここで、下記数３３が成立する。

行列

の逆行列が存在しないことも考えられる。しかしながら、行列

はδ_ｋまたは倍率Ｍの関数と考えられ、極限値の意味で逆行列が存在すると考えることができる。

数３２の両辺にフーリエ変換に相当する回転行列

をかけることにより、評価関数の微分の空間周波数成分を下記数３４及び数３５のように表すことができる。

ここで、

のｉ番目の要素は

の最大値である。数３５の不等式は、ＰＳＦの平均パワースペクトラム

によって、評価関数の微分の空間周波数成分が抑え込まれるということを表している。

つまり、ＰＳＦの平均パワースペクトラムのある空間周波数成分が非常に小さければ、対応する評価関数の微分の空間周波数成分も小さくなる。評価関数の微分の空間周波数成分が小さいということは、高解像度画像のその空間周波数成分がほとんど更新されないということを示している。このような制限を「勾配制限」と呼ぶことにする。

特に、ＰＳＦの平均パワースペクトラムが０であれば、その空間周波数成分において高解像度画像が全く更新されないことを意味する。ここで、ＰＳＦの平均パワースペクトラムが０である空間周波数を「特異な空間周波数」と呼ぶことにする。

以上では、一般的なＰＳＦに関して理論的な解析を行った。以下では、具体的にBox型ＰＳＦ及びGaussian型ＰＳＦを取り上げ、それぞれの条件数及び平均パワースペクトラムを算出し、解析を行う。なお、上述した一般的なＰＳＦに関する理論解析は、１次元（超解像処理が扱う低解像度画像は１次元画像である）での解析であったが、以下では２次元（超解像処理が扱う低解像度画像は２次元画像）の場合の結果を示す。

まず、Box型ＰＳＦ及びGaussian型ＰＳＦ（σ＝０.３）に関して、倍率Ｍを変化させながら、条件数を算出した。図２にその条件数の逆数を対数軸で表したときの結果を示す。条件数の逆数とは、その値が大きいほど超解像処理の基本方程式の解が安定に推定されていることを意味し、特に条件数の逆数の値がゼロの場合はその基本方程式が特異な問題であることを表す。

なお、図２において、低解像度画像のサイズを２３×２３とし、数２４の積分計算には数値積分を利用した。また、Gaussian型ＰＳＦのσは、倍率Ｍ＝１のときに、条件数がBox型ＰＳＦの条件数と等しくなるようにσ＝０.３とした。

図２から分かるように、Box型ＰＳＦの場合、２以上の整数倍率において条件数の逆数がゼロとなっていることが特徴的である。これは、Box型ＰＳＦを仮定した場合に、２以上の整数倍率は特異な問題となっていることを示しており、非特許文献１１及び非特許文献１２に記載されるように、ベイカ（Baker）ら及びリン（Lin）らの従来研究と一致している。さらに、Box型ＰＳＦの条件数の逆数の大まかな形状もLinらの定量的な解析と一致する（非特許文献１２参照）。

本発明では、上述した超解像の条件数算出方法により、Box型ＰＳＦのみならず、任意の空間一定なＰＳＦの条件数を算出できるため、例えば、Gaussian型ＰＳＦに関しても同様に条件数を算出することができる。

図２から分かるように、例えば、倍率が１から２.２４の範囲では、Gaussian型ＰＳＦの条件数の逆数がBox型ＰＳＦの条件数の逆数を上回っているため、この範囲の倍率に対してはGaussian型ＰＳＦを利用するほうが適切であることが一目瞭然である。それ以外の倍率に関しても、条件数の逆数を比較することにより、どちらのＰＳＦ近似を利用するのが適切であるかを判断することが可能である。

つまり、超解像処理において、ＰＳＦの設計に関して、同じ倍率の場合では、条件数の逆数が大きいほうのＰＳＦを用いることにより、より安定に高解像度画像を推定できる。

次に、超解像処理において、倍率Ｍの設計について説明する。まず、Box型ＰＳＦの場合について考える。図２から分かるように、例えば、倍率Ｍが２.８０から３.４４の範囲では、Ｍ＝３.４４の条件数の逆数が最も大きい。従って、Box型ＰＳＦを仮定した下で、２.８０から３.４４までの範囲の倍率Ｍが必要な場合は、３.４４倍の倍率でまず、超解像処理をして、そして必要な倍率に縮小するほうが効果的であることが分かる。

次に、Gaussian型ＰＳＦの場合について考える。図２から分かるように、条件数の逆数は、倍率に関して単調に減少している。従って、Gaussian型ＰＳＦを用いた場合に、できるだけ小さい倍率を利用するほうが効果的であることが分かる。

勾配制限は、ＰＳＦの平均パワースペクトラムが小さければ、ＭＬ法の評価関数の空間周波数成分も小さいということを示している。つまり、ＭＬ法の評価関数の勾配は、ＰＳＦの平均パワースペクトラムにより制限される。具体的に、Box型ＰＳＦ及びGaussian型ＰＳＦに関して、倍率Ｍ＝２.０,２.５,３.０の場合の平均パワースペクトラムを計算した。その結果を図３に示す。

図３より、倍率によらずGaussian型ＰＳＦの平均パワースペクトラムは、空間周波数に対して単調に減少しているのが確認される。その一方で、Box型ＰＳＦの平均パワースペクトラムは特定の空間周波数において、非常に小さな値となっている。特に、倍率がＭ＝２.０,３.０では、黒で示されている空間周波数の平均パワースペクトラムは０であり、特異な空間周波数となっている。特異な空間周波数はＭ＝２.０の場合、０.５[ＨＲＩpixel^−１]であり、Ｍ＝３.０の場合、０.３３３[ＨＲＩpixel^−１]である。勾配制限により、高解像度画像を更新する際に、この特異な空間周波数成分は、全く更新されないことが予想される。

次に、本発明における超解像の条件数算出方法及び勾配制限の有効性を確認するために、合成画像を用いて、実際に超解像処理を行った。

まず、図４に示す画像を真の高解像度画像として用意した。この画像に対して、Box型PSFを仮定し、２５６枚の低解像度画像を作成した。このとき、低解像度画像の大きさを２３×２３、平行移動の間隔は低解像度画像の画素間隔の１／１６刻みとしている。位置あわせは正確であると仮定している。

次に、勾配制限の効果が顕著に現れるように、図５に示す３種類の初期画像を用意した。図５に示される３種類の初期画像は、図３(Ａ)のBox型ＰＳＦの平均パワースペクトラムが最小値となる空間周波数成分のみを、それぞれ有している。

そして、Box型ＰＳＦとGaussian型ＰＳＦ(σ＝０.３)を仮定して、ＭＬ法による超解像処理を行った。このとき、繰返し数は２００回として、倍率はＭ＝２.０,２.５および３.０の３種類に対して超解像処理を行った。超解像処理結果を図６に示す。

まず、Box型ＰＳＦ(Ｍ＝２.５)とGaussian型ＰＳＦ(Ｍ＝２.０)の条件数の逆数に注目すると、それぞれ１.１５×１０^−２、２.４０×１０^−２である。これらの条件数の逆数は、比較的大きな値を示している。対応する超解像処理結果も、他の条件と比較した場合，初期画像によらずノイズが小さく良好な結果を示している。超解像の条件数算出方法から導かれる条件数の逆数からも、実験結果からも、Box型ＰＳＦ(Ｍ＝２.５)とGaussian型ＰＳＦ(Ｍ＝２.０)の２つのＰＳＦ条件が、他のＰＳＦ条件よりも安定であるという一致した結果が得られた。

本発明の超解像の条件数算出方法により算出される条件数も、非特許文献１１及び非特許文献１２に記載された従来の研究においても、Box型ＰＳＦの場合、倍率が２以上の整数では特異な問題となるため、解が安定に得られないということが示されている。

しかしながら、Ｍ＝２.０とＭ＝３.０のBox型ＰＳＦに対する超解像結果は、図６に示されるように、幾つかの初期画像に対しては、ノイズがなく良好な結果となっている。このような初期画像に対する依存性は、従来の研究では説明することができない。

一方、本発明の勾配制限の関係を利用することにより、この初期画像に対する依存性を簡単に説明することが可能である。例として、Box型ＰＳＦ(Ｍ＝３.０)の場合を考える。ノイズの多い結果を与えている図５(Ｃ)の初期画像は、(０.３３３,０.３３３)の空間周波数成分のみをもつ画像である。

また、図３(Ａ)に示すように、Box型ＰＳＦ(Ｍ＝３.０)の平均パワースペクトラムは、(０.３３３,０.３３３)において０であり、この空間周波数は特異な空間周波数であることが確認できる。従って、初期画像がもつ特異な空間周波数成分は、ＭＬ法では更新されることなく、超解像処理結果に現われてしまう。これが超解像処理結果のノイズである。その一方で、図５(Ａ)及び図５(Ｂ)の初期画像が有している空間周波数成分に関するBox型ＰＳＦ(Ｍ＝３.０)の平均パワースペクトラムは、０ではない。従って、初期画像は更新され、また、超解像処理の結果にも初期画像の成分は現われず、良好な結果を示す。

次に、Gaussian型ＰＳＦに注目する。Gaussian型ＰＳＦの平均パワースペクトラムは、空間周波数に対して単調に減少する。本発明では、高解像度画像の画素間隔を１と仮定しているので、空間周波数の最大値は０.５である。従って、Gaussian型ＰＳＦの平均パワースペクトラムは、空間周波数が(０.５,０.５)において最小となる。

以上のことから、空間周波数(０.５,０.５)の成分のみを有する図５(Ａ)の初期画像に対する超解像処理結果が、最も悪いことが予想される。図６の結果も、この予想と一致している。特に、図５(Ａ)の初期画像に対する超解像処理結果を、倍率を変化させながら比較すると、倍率を大きくするに従って、ノイズの大きくなっていることが確認される。この結果も、Gaussian型ＰＳＦの条件数の逆数が倍率に対して単調減少するという理論解析の結果と一致している。

以上のように、本発明では、超解像の条件数算出方法を導出し、また証明を行った。超解像の条件数算出方法によれば、任意のＰＳＦに対して、任意の倍率における条件数を算出することができる。つまり、本発明の超解像の条件数算出方法によれば、倍率ＭとＰＳＦが与えられた場合に、超解像の条件数を算出できる。従って、条件数を比較することにより、超解像処理のパラメータであるＰＳＦ及び倍率Ｍを最適に設計することができる。

また、本発明では、具体的な超解像処理方法としてＭＬ法を用いた場合、その評価関数の勾配に関する勾配制限についても検討した。さらに、本発明では、超解像の条件数算出方法及び勾配制限の有効性を、合成画像を利用して行った超解像処理により確認した。

超解像処理において、高解像度画像、低解像度画像及びＰＳＦに関してそれぞれ連続と離散との関係を説明するための模式図である。本発明におけるBox型ＰＳＦ、Gaussian型ＰＳＦをそれぞれ仮定した場合の倍率Ｍと条件数の逆数との関係を示す図である。本発明におけるBox型ＰＳＦ、Gaussian型ＰＳＦをそれぞれ仮定した下で異なる倍率のＰＳＦの平均パワースペクトラムを示す図である。本発明を適用した超解像処理において、用いられた真の高解像度画像を示す図である。本発明を適用した超解像処理において、用いられた３種類の初期画像を示す図である。本発明を適用した超解像処理の処理結果を示す図である。

Claims

位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法であって、
前記パラメータは、前記高解像度画像の前記低解像度画像に対する倍率、及びカメラモデルから得られるＰＳＦ関数であり、
前記超解像処理の基本方程式の解を求める際の安定性を表す、前記基本方程式の条件数に基いて、前記高解像度画像を安定に推定できる最適化パラメータを決定し、
前記基本方程式は、次の式で表される線形連立方程式であり、

ただし、

は前記複数の低解像度画像の全ての画素値を表すベクトルで、

は推定すべき前記高解像度画像の画素値を表すベクトルで、

は前記ＰＳＦ関数及び前記倍率から導かれるシステム行列を表し、
（Ａ１）前記位置ずれは平行移動であり、
（Ａ２）前記ＰＳＦ関数は空間不変であり、
（Ａ３）前記高解像度画像の画素値は同一画素内で一定であり、
（Ａ４）前記低解像度画像は位置の重なりがなく、無限に取得できるといった４つの条件を前提に、前記条件数を算出することを特徴とする超解像処理のパラメータ設計方法。
前記条件数cond(Ｂ)は、次の式により算出され、

ただし、

が成立し、
また、

が成立し、

の離散フーリエ変換であり、ｂ(ｘ,ｙ)は高解像度画像空間における連続のＰＳＦ関数を表す請求項１に記載の超解像処理のパラメータ設計方法。
所定のＰＳＦ関数に対して、異なる倍率における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応する倍率を、前記超解像処理に前記所定のＰＳＦ関数が用いられた場合の最適な倍率にする請求項２に記載の超解像処理のパラメータ設計方法。
所定の倍率に対して、異なるＰＳＦ関数における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応するＰＳＦ関数を、前記超解像処理に前記所定の倍率が指定された場合の最適化ＰＳＦ関数にする請求項２に記載の超解像処理のパラメータ設計方法。