(57)【要約】【課題】  再構成型超解像処理において、高解像度画像を安定に推定できるパラメータを設計できる超解像処理のパラメータ設計方法を提供する。【解決手段】  位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法であって、前記パラメータは、前記高解像度画像の前記低解像度画像に対する倍率、及びカメラモデルから得られるＰＳＦ関数であり、前記超解像処理の基本方程式の解を求める際の安定性を表す定量的指標に基いて、前記高解像度画像を安定に推定できる最適化パラメータを決定する。また、前記定量的指標として、前記基本方程式の条件数を用いる。【選択図】  図１
2006195856.tif 000002
【特許請求の範囲】
【請求項１】
位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法であって、
前記パラメータは、前記高解像度画像の前記低解像度画像に対する倍率、及びカメラモデルから得られるＰＳＦ関数であり、
前記超解像処理の基本方程式の解を求める際の安定性を表す定量的指標に基いて、前記高解像度画像を安定に推定できる最適化パラメータを決定することを特徴とする超解像処理のパラメータ設計方法。
【請求項２】
前記基本方程式は、次の式で表される線形連立方程式であり、
2006195856.tif 000088
ただし、
2006195856.tif 000089
は前記複数の低解像度画像の全ての画素値を表すベクトルで、
2006195856.tif 000090
は推定すべき前記高解像度画像の画素値を表すベクトルで、
2006195856.tif 000091
は前記ＰＳＦ関数及び前記倍率から導かれるシステム行列を表し、
（Ａ１）前記位置ずれは平行移動であり、
（Ａ２）前記ＰＳＦ関数は空間不変であり、
（Ａ３）前記高解像度画像の画素値は同一画素内で一定であり、
（Ａ４）前記低解像度画像は位置の重なりがなく、無限に取得できるといった４つの条件を前提に、前記定量的指標を算出する請求項１に記載の超解像処理のパラメータ設計方法。
【請求項３】
前記定量的指標として、前記基本方程式の条件数を用いる請求項２に記載の超解像処理のパラメータ設計方法。
【請求項４】
前記条件数cond(Ｂ)は、次の式により算出され、
2006195856.tif 000092
ただし、
2006195856.tif 000093
が成立し、また、
2006195856.tif 000094
が成立し、
2006195856.tif 000095
は
2006195856.tif 000096
の離散フーリエ変換であり、ｂ(ｘ,ｙ)は高解像度画像空間における連続のＰＳＦ関数を表す請求項３に記載の超解像処理のパラメータ設計方法。
【請求項５】
所定のＰＳＦ関数に対して、異なる倍率における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応する倍率を、前記超解像処理に前記所定のＰＳＦ関数が用いられた場合の最適な倍率にする請求項４に記載の超解像処理のパラメータ設計方法。
【請求項６】
所定の倍率に対して、異なるＰＳＦ関数における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応するＰＳＦ関数を、前記超解像処理に前記所定の倍率が指定された場合の最適化ＰＳＦ関数にする請求項４に記載の超解像処理のパラメータ設計方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法に関し、特に、再構成型超解像処理において、高解像度画像を安定に推定できるパラメータを設計できる超解像処理のパラメータ設計方法に関する。
【背景技術】
【０００２】
超解像処理とは、位置ずれのある複数の低解像度画像（観測画像）より一つの高解像度画像を推定する処理である。高解像度画像から低解像度画像への画像変換は、点広がり関数(ＰＳＦ)と呼ばれる積分核との積分変換として、定式化される。超解像処理は、この積分変換の逆変換に相当する処理である。
【０００３】
これまでにも、例えば、非特許文献１〜非特許文献６に開示されるように、数多くの超解像処理方法が提案されている。その中でも、再構成型超解像処理と分類される方法が広く利用されている。例えば、非特許文献２に開示されているＭＬ（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型ＰＳＦを仮定した場合のみが検討された。
【００１２】
以上のように、従来では、再構成型超解像処理に必要なパラメータ（つまり、カメラモデルから得られるＰＳＦ関数、高解像度画像の低解像度画像に対する倍率Ｍ）を最適に設計する超解像処理のパラメータ設計方法は、研究されていない。
【００１３】
本発明は、上述のような事情よりなされたものであり、本発明の目的は、再構成型超解像処理において、高解像度画像を安定に推定できるパラメータを設計できる超解像処理のパラメータ設計方法を提供することにある。
【課題を解決するための手段】
【００１４】
本発明は、位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理に必要なパラメータを設計するための超解像処理のパラメータ設計方法に関し、本発明の上記目的は、前記パラメータは、前記高解像度画像の前記低解像度画像に対する倍率、及びカメラモデルから得られるＰＳＦ関数であり、前記超解像処理の基本方程式の解を求める際の安定性を表す定量的指標に基いて、前記高解像度画像を安定に推定できる最適化パラメータを決定することにより、或いは、前記基本方程式は、次の式で表される線形連立方程式であり、
2006195856.tif 000003
ただし、
2006195856.tif 000004
は前記複数の低解像度画像の全ての画素値を表すベクトルで、
2006195856.tif 000005
は推定すべき前記高解像度画像の画素値を表すベクトルで、
2006195856.tif 000006
は前記ＰＳＦ関数及び前記倍率から導かれるシステム行列を表し、（Ａ１）前記位置ずれは平行移動であり、（Ａ２）前記ＰＳＦ関数は空間不変であり、（Ａ３）前記高解像度画像の画素値は同一画素内で一定であり、（Ａ４）前記低解像度画像は位置の重なりがなく、無限に取得できるといった４つの条件を前提に前記定量的指標を算出することにより、或いは、前記定量的指標として、前記基本方程式の条件数を用いることによって効果的に達成される。
【００１５】
また、本発明の上記目的は、前記条件数cond(Ｂ)は、次の式により算出され、
2006195856.tif 000007
ただし、
2006195856.tif 000008
が成立し、また、
2006195856.tif 000009
が成立し、
2006195856.tif 000010
は
2006195856.tif 000011
の離散フーリエ変換であり、ｂ(ｘ,ｙ)は高解像度画像空間における連続のＰＳＦ関数を表すことにより、或いは、所定のＰＳＦ関数に対して、異なる倍率における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応する倍率を、前記超解像処理に前記所定のＰＳＦ関数が用いられた場合の最適な倍率にすることにより、或いは、所定の倍率に対して、異なるＰＳＦ関数における条件数をそれぞれ算出し、算出された条件数の中で、最小値となる条件数に対応するＰＳＦ関数を、前記超解像処理に前記所定の倍率が指定された場合の最適化ＰＳＦ関数にすることによって一層効果的に達成される。
【発明の効果】
【００１６】
本発明では、線形連立方程式の解を求める際の安定性を表す定量的指標として利用されている条件数に注目し、任意の空間不変なＰＳＦ、任意の倍率に対して、超解像の基本方程式の条件数を算出する方法である条件数算出方法を開示し、そして、その条件数算出方法を証明した。
【００１７】
本発明に係る超解像処理のパラメータ設計方法によれば、任意の空間不変なＰＳＦ、任意の倍率に対して、超解像の基本方程式の条件数を算出することが可能であり、そして、それぞれの条件数を比較することにより、より良い超解像結果を得るためには、所望の倍率ではどのＰＳＦ近似が良いか、また、指定されたＰＳＦ近似では倍率をどの値に設定すれば良いかを判断することができ、よって、高解像度画像を安定に推定できる最適なパラメータを設計できるという優れた効果を奏する。
【００１８】
また、合成画像に対して、本発明に係る超解像処理のパラメータ方法を適用することによって、高解像度画像を安定に推定できる最適化ＰＳＦ関数、最適な倍率を設定できることを確認した。
【発明を実施するための最良の形態】
【００１９】
以下、本発明を実施するための最良の形態を図面を参照して説明する。なお、本発明では、『復元したい高解像度画像の空間』を『高解像度画像空間』と称する。
【００２０】
まず、超解像の基本方程式について説明する。
【００２１】
再構成型超解像処理の基本となる連続系の画像変換方程式を下記数１に示す。
【００２２】
【数１】
2006195856.tif 000012
ここで、添え字ｊはｊ番目の観測を、ｙｊは観測画像の座標を、ｘは高解像度画像の座標を、ｇｊ(ｙｊ)は連続観測画像を、ｂｊ(ｙｊ,ｘ)は連続のＰＳＦを、ｈ(ｘ)は連続高解像度画像を、それぞれ表す。
【００２３】
実際に、超解像処理が扱う方程式は、上記数１を離散化して得られた下記数２である。以下、この数２を超解像の方程式と呼ぶことにする。
【００２４】
【数２】
2006195856.tif 000013
ここで、
2006195856.tif 000014
は複数の低解像度画像の全ての画素値を表すベクトルを、
2006195856.tif 000015
は推定すべき高解像度画像の画素値を表すベクトルを、
2006195856.tif 000016
は高解像度画像から低解像度画像への変換を表す行列（つまり、ＰＳＦ、倍率Ｍ及び低解像度画像画素のレジストレーション後の画素位置から導出されるシステム行列）を、それぞれ表す。
【００２５】
数２で表す超解像の方程式は、低解像度画像が無限にある場合にも拡張することができる。ここで、低解像度画像が無限の場合（低解像度画像の枚数が無限の場合）の超解像の方程式を、超解像の基本方程式と呼ぶことにする。超解像の基本方程式は、上記数２と同様であるが、
2006195856.tif 000017
の次元と
2006195856.tif 000018
の行の次元が無限大である。
【００２６】
次に、本発明の理論解析について説明する。
【００２７】
本発明の理論解析は以下の４つの仮定に基づく。
【００２８】
つまり、まず、「仮定１」として、位置ずれは平行移動である。また、「仮定２」として、カメラモデルから得られるＰＳＦ関数（以下、単に、カメラＰＳＦ、もしくは、ＰＳＦとも称する）は、空間不変である。そして、「仮定３」として、連続高解像度画像の画素値は同一画素内で一定で、つまり、高解像度画像の画素内の濃度値は一定である。最後に、「仮定４」として、低解像度画像は、位置の重なりがなく、且つ、無限に取得することができる。
【００２９】
本発明の大きな特徴の１つは、本発明で用いられる「仮定２」（つまり、カメラＰＳＦを単に空間不変と仮定している点）が、これまでの従来研究（非特許文献１１及び非特許文献１２参照）と大きく異なる点である。
【００３０】
ここで、一般性を失わず、高解像度画像の画素間隔を１として解析を行うことができる。ただし、高解像度画像の画素間隔を１と考えると、ＰＳＦは高解像度画像の低解像度画像に対する倍率Ｍをパラメータとして含むようになる。この場合の具体的なＰＳＦの例として、Box型ＰＳＦを下記数３に、Gaussian型ＰＳＦを下記数４に、それぞれ示す。
【００３１】
【数３】
2006195856.tif 000019
【００３２】
【数４】
2006195856.tif 000020
ここで、
2006195856.tif 000021
は下記数５で定義される矩形関数である。
【００３３】
【数５】
2006195856.tif 000022
以下では、簡単のため、再構成型超解像処理を１次元の問題として（つまり、超解像処理が扱う低解像度画像は１次元画像であるを前提に）、本発明の説明を行うが、本発明はそれに限られることがなく、２次元の問題（つまり、超解像処理が扱う低解像度画像は２次元画像であるの場合）へ簡単に拡張できることは言うまでもない。
【００３４】
本発明では、超解像の基本方程式を以下のように導出する。
【００３５】
上述した４つの仮定（「仮定１」〜「仮定４」）を利用して、低解像度画像を無限に取得できると仮定した場合の超解像の基本方程式の導出を行う。「仮定１」と「仮定２」を利用すれば、数１で表す画像変換方程式から下記数６を得ることができる。
【００３６】
【数６】
2006195856.tif 000023
ここで，ｙはレジストレーション後の低解像度画像の座標である。「仮定３」を満足する連続高解像度画像は、数５で表す矩形関数を利用して、下記数７のように定式化することができる。
【００３７】
【数７】
2006195856.tif 000024
ここで、Ｎは１次元画像の大きさを、ｉは１次元画像の離散座標を、ｈ[ｉ]は離散高解像度画像の位置ｉにおける画像値を、それぞれ表す。数６及び数７より、数８に示される離散高解像度画像から連続低解像度画像への画像変換方程式を得ることができる。
【００３８】
【数８】
2006195856.tif 000025
さて、ある座標ｙを下記数９のように、ｙに最も近い整数
2006195856.tif 000026
と残りの小数部分δに分けて考える。
【００３９】
【数９】
2006195856.tif 000027
このとき、小数部分δの変域は−１／２から１／２となる。数９の関係を利用して、数８を下記数１０の離散畳込み演算の形式に簡単化することができる。
【００４０】
【数１０】
2006195856.tif 000028
ここで、
2006195856.tif 000029
はオフセットδの離散ＰＳＦであり、下記数１１で定義される。
【００４１】
【数１１】
2006195856.tif 000030
次に、小数部分δを固定して考えると、
2006195856.tif 000031
は１つの離散画像を表しているとみなせる。この離散画像は、連続観測画像ｇ(ｙ)を高解像度画像の画素間隔でサンプリングした画像と同一である。ただし、サンプリング点は、高解像度画像の画素中心とδだけずれている。高解像度画像、低解像度画像及びＰＳＦに関して、それぞれ連続と離散の関係を図１に示す。
【００４２】
ここで、離散画像
2006195856.tif 000032
を「オフセットδの離散低解像度画像」と呼び、
2006195856.tif 000033
で表すことにする。数１０は、オフセットδの離散低解像度画像に対する超解像の方程式となる。数１０をベクトル表現で表すと、下記数１２となる。
【００４３】
【数１２】
2006195856.tif 000034
ここで、
2006195856.tif 000035
は離散畳込み演算のカーネルを表す行列である。
【００４４】
連続低解像度画像は、オフセットδの離散低解像度画像の線形結合として表される。従って、連続低解像度画像に対する超解像の方程式、即ち、超解像の基本方程式も同様に、オフセットδの離散低解像度画像に対する方程式の線形結合として、下記数１３のように表現することができる。
【００４５】
【数１３】
2006195856.tif 000036
ここで、下記数１４が成立し、また、ｉ≠ｊのときに、δｉ≠δｊである。
【００４６】
【数１４】
2006195856.tif 000037
以上のように、前述した４つの仮定から、数１３で表す超解像の基本方程式を導出することができた。
【００４７】
次に、上述したように導出された超解像の基本方程式に対する条件数の算出方法（以下、単に「超解像の条件数算出方法」とも称する）について説明する。
【００４８】
ここで、まず、条件数に関する数学的な性質を概説した後に、本発明のポイントである超解像の条件数算出方法を示し、証明する。
【００４９】
条件数は、線形連立方程式の係数行列に対して定義され、連立方程式の解を求める際の安定性の指標として知られている（非特許文献９参照）。条件数の小さい連立方程式は、解が安定に推定される（推定誤差が相対的に小さい）。また、条件数の大きい連立方程式は、解が不安定である（推定誤差が相対的に大きい）。特に、条件数が無限大のときは、その連立方程式は特異で、つまり、解が一意に求められないことを意味する。
【００５０】
ところで、数２で示されるような連立方程式（つまり、超解像の方程式）に対して、その条件数
2006195856.tif 000038
、つまり、超解像の方程式の条件数（以下、単に条件数とも称する）は、下記数１５のように定義することができる。
【００５１】
【数１５】
2006195856.tif 000039
ここで、‖・‖２はＬ２ノルムを表す。条件数に関して、下記数１６及び数１７で表す不等式が成立し、条件数は解の相対誤差の上限を定める。
【００５２】
【数１６】
2006195856.tif 000040
【００５３】
【数１７】
2006195856.tif 000041
ここで、Δｈは推定された高解像度画像の誤差を、Δｇは観測された低解像度画像の誤差を、ΔＢはＰＳＦによって定まる係数行列の誤差を、それぞれ表す。
【００５４】
ちなみに、数１７の不等式は、ＰＳＦの近似精度が多少悪くても条件数が小さければ、結果として誤差の影響は小さく抑えられることを示している。
【００５５】
ところで、行列のノルムを直接計算することは一般に困難であるので、条件数を上記数１５の定義に従って計算することは現実的ではない。
【００５６】
そこで、本発明では、下記数１８に示される条件数の性質を利用する（非特許文献９参照）。
【００５７】
【数１８】
2006195856.tif 000042
ここで、＊は行列の共役転置を、
2006195856.tif 000043
は行列
2006195856.tif 000044
の固有値をベクトル化したものを、それぞれ表す。
【００５８】
次に、条件数算出方法について詳細に説明する。以下、補題を示した後に、条件数算出方法を証明する。
＜補題＞数１３の係数行列
2006195856.tif 000045
に対して、下記数１９が成り立つ。
【００５９】
【数１９】
2006195856.tif 000046
ここで、
2006195856.tif 000047
は
2006195856.tif 000048
の離散フーリエ変換を表し、｜・｜２は要素ごとに適用されるとする。
＜補題の証明＞
先ず、オフセットδの離散低解像度画像に対する超解像の方程式の係数行列
2006195856.tif 000049
は、下記数２０のように対角化される。
【００６０】
【数２０】
2006195856.tif 000050
ここで、
2006195856.tif 000051
は離散フーリエ変換の行列表現に対応するユニタリ回転行列を、
2006195856.tif 000052
は
2006195856.tif 000053
の要素を対角成分にもつ対角行列を、それぞれ表す。
【００６１】
次に、行列
2006195856.tif 000054
の固有値を求めるため、対角化を行う。行列
2006195856.tif 000055
の定義式と数２０の関係式より、行列
2006195856.tif 000056
の対角化を下記数２１のように行うことができる。
【００６２】
【数２１】
2006195856.tif 000057
ここで、行列
2006195856.tif 000058
はユニタリ行列であるので、数２１の右辺の対角成分は、行列
2006195856.tif 000059
の固有値となる。従って、数１９が得られる。
＜超解像の条件数算出方法＞
数１３で表す超解像の基本方程式に対する条件数
2006195856.tif 000060
は、下記数２２により算出される。
【００６３】
【数２２】
2006195856.tif 000061
ここで、下記数２３及び数２４が成立する。また、
2006195856.tif 000062
をＰＳＦの平均パワースペクトラムと呼ぶことにする。
【００６４】
【数２３】
2006195856.tif 000063
【００６５】
【数２４】
2006195856.tif 000064
以上では、超解像処理が扱う低解像度画像が１次元画像の場合（つまり、超解像処理を１次元の問題として考えた場合）に、超解像の条件数算出方法を詳細に説明した。
【００６６】
なお、超解像処理を２次元の問題として考えた場合に、要するに、超解像処理が扱う低解像度画像が２次元画像である場合に、超解像の基本方程式に対する条件数
2006195856.tif 000065
は、下記数２５により算出される。
【００６７】
【数２５】
2006195856.tif 000066
ただし、下記数２６及び数２７が成立する。また、
2006195856.tif 000067
は
2006195856.tif 000068
の離散フーリエ変換であり、そして、ｂ(ｘ,ｙ)は高解像度画像空間における連続のＰＳＦ関数を表す。
【００６８】
【数２６】
2006195856.tif 000069
【００６９】
【数２７】
2006195856.tif 000070
＜超解像の条件数算出方法の証明＞
ここで、数２２に表される超解像の条件数算出方法を以下のように証明する。
【００７０】
規格化されたカメラＰＳＦを考えた場合に、下記数２８の関係を仮定することができる。
【００７１】
【数２８】
2006195856.tif 000071
数２８を利用することにより、行列
2006195856.tif 000072
に対する数１８の分子は、下記数２９となる。
【００７２】
【数２９】
2006195856.tif 000073
上記数２９の右辺は無限大に発散してしまうが、条件数は数１８に示されるように最小固有値と最大固有値の比であるので、条件数は収束する。
【００７３】
数１８、数１９及び数２９より、下記数３０が得られる。
【００７４】
【数３０】
2006195856.tif 000074
上記数３０を簡単にすることにより、数２２が得られる。
【００７５】
なお、δｋの変域は−１／２から１／２であり、δｋが重ならず無限あると考えた場合に、数２３を数２４の積分の形式に変形することができる。
【００７６】
超解像の条件数算出方法を利用して、超解像の基本方程式に対する条件数を算出するために、必要なパラメータは、次の３つのパラメータである。
【００７７】
つまり、「パラメータ１」として連続ＰＳＦｂ(ｘ)を、「パラメータ２」として高解像度画像の低解像度画像に対する倍率Ｍを、「パラメータ３」として低解像度画像の大きさをそれぞれ用いる。
【００７８】
「パラメータ１」及び「パラメータ２」が必要な理由は、条件数算出方法より明らかである。「パラメータ３」に関しては、ＰＳＦの平均パワースペクトラムを計算するために、必要である。
【００７９】
つまり、超解像の条件数算出方法は、任意のＰＳＦに対して、任意の倍率での条件数を算出することを可能にする。得られた条件数は、ＰＳＦ近似に対する重要な指標となる。従って、条件数を比較し、条件数のより小さなＰＳＦ近似を利用することにより、高解像度画像を安定に推定することができる。また、同一のＰＳＦ近似で倍率を変化させた場合も、条件数を比較することにより、高解像度画像を安定に推定することができる最適な倍率を設定することができ、更に、倍率の限界を示すことも可能にする。
【００８０】
次に、数１３で表す超解像の基本方程式をＭＬ法で解く場合を考える。この場合に、ＰＳＦの平均パワースペクトラム
2006195856.tif 000075
がＭＬ法の評価関数の勾配とどのような関係があるかを明らかにする。つまり、ＭＬ法における勾配とＰＳＦのパワースペクトルとの関係を勾配制限として示す。
【００８１】
数１３で表す超解像の基本方程式に関するＭＬ法の評価関数を下記数３１のように定義する。
【００８２】
【数３１】
2006195856.tif 000076
ここで、評価関数が発散するのを防ぐために、Ｋによる規格化を行っている。評価関数の高解像度画像
2006195856.tif 000077
に関する微分は、下記数３２となる。
【００８３】
【数３２】
2006195856.tif 000078
ここで、下記数３３が成立する。
【００８４】
【数３３】
2006195856.tif 000079
行列
2006195856.tif 000080
の逆行列が存在しないことも考えられる。しかしながら、行列
2006195856.tif 000081
はδｋまたは倍率Ｍの関数と考えられ、極限値の意味で逆行列が存在すると考えることができる。
【００８５】
数３２の両辺にフーリエ変換に相当する回転行列
2006195856.tif 000082
をかけることにより、評価関数の微分の空間周波数成分を下記数３４及び数３５のように表すことができる。
【００８６】
【数３４】
2006195856.tif 000083
【００８７】
【数３５】
2006195856.tif 000084
ここで、
2006195856.tif 000085
のｉ番目の要素は
2006195856.tif 000086
の最大値である。数３５の不等式は、ＰＳＦの平均パワースペクトラム
2006195856.tif 000087
によって、評価関数の微分の空間周波数成分が抑え込まれるということを表している。
【００８８】
つまり、ＰＳＦの平均パワースペクトラムのある空間周波数成分が非常に小さければ、対応する評価関数の微分の空間周波数成分も小さくなる。評価関数の微分の空間周波数成分が小さいということは、高解像度画像のその空間周波数成分がほとんど更新されないということを示している。このような制限を「勾配制限」と呼ぶことにする。
【００８９】
特に、ＰＳＦの平均パワースペクトラムが０であれば、その空間周波数成分において高解像度画像が全く更新されないことを意味する。ここで、ＰＳＦの平均パワースペクトラムが０である空間周波数を「特異な空間周波数」と呼ぶことにする。
【００９０】
以上では、一般的なＰＳＦに関して理論的な解析を行った。以下では、具体的に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型ＰＳＦをそれぞれ仮定した下で異なる倍率のＰＳＦの平均パワースペクトラムを示す図である。
【図４】本発明を適用した超解像処理において、用いられた真の高解像度画像を示す図である。
【図５】本発明を適用した超解像処理において、用いられた３種類の初期画像を示す図である。
【図６】本発明を適用した超解像処理の処理結果を示す図である。
【図１】
2006195856.tif 000097
【図２】
2006195856.tif 000098
【図３】
2006195856.tif 000099
【図４】
2006195856.tif 000100
【図５】
2006195856.tif 000101
【図６】
2006195856.tif 000102