(57)【要約】【課題】高速かつ高画質で被写体の画像を再構成可能な画像処理方法、画像処理装置、撮像装置および画像処理プログラムを提供すること。【解決手段】符号化マスクを用いた被写体の撮像に基づく符号化データを取得するステップと、符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出するステップと、推定画像に対する複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出するステップと、誤差画像を用いて推定画像を更新するステップと、を有する。【選択図】図３
2018206274.tif 000002
【特許請求の範囲】
【請求項１】
符号化マスクを用いた被写体の撮像に基づく符号化データを取得するステップと、
前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出するステップと、
前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出するステップと、
前記誤差画像を用いて前記推定画像を更新するステップと、を有することを特徴とする画像処理方法。
【請求項２】
前記符号化マスクの透過率分布は、前記複素関数の実部および虚部の重み付き平均により表されることを特徴とする請求項１に記載の画像処理方法。
【請求項３】
前記推定画像を更新するステップは、
前記符号化データから前記誤差画像を減算することで減算画像を算出するステップと、
前記減算画像に対する前記複素関数マスクによる逆畳み込み演算を行うことで前記推定画像を算出するステップと、を有することを特徴とする請求項１または２に記載の画像処理方法。
【請求項４】
前記複素関数は、複素関数フィルタの逆フーリエ変換で表され、
前記逆畳み込み演算は、
フーリエ変換によって、前記符号化データまたは前記符号化データから前記誤差画像を減算することで算出される減算画像から変換データを算出するステップと、
前記変換データに対して、前記複素関数フィルタで除算した後、逆フーリエ変換を行うステップと、を有することを特徴とする請求項１から３のいずれか１項に記載の画像処理方法。
【請求項５】
前記複素関数は、複素関数フィルタの逆フーリエ変換で表され、
前記複素関数フィルタは、２次元の各方向の空間周波数をそれぞれξ，ηとするとき、正の実関数ａ（ξ,η）および実関数Φ（ξ,η）を用いて、ａ（ξ,η）×ｅｘｐ（ｊΦ（ξ,η））と表され、
前記正の実関数ａ（ξ,η）は、前記符号化データのサンプリング間隔をΔ、前記正の実関数ａ（ξ,η）の最大値をａｍａｘとするとき、
（ξ２＋η２）１／２≦１／２Δ
なる条件式を満足する領域内で、
ａ（ξ,η）≧０．１５ａｍａｘ
なる条件式を満足することを特徴とする請求項１から４のいずれか１項に記載の画像処理方法。
【請求項６】
前記複素関数マスクは、同心円状のパターンを有することを特徴とする請求項１から５のいずれか１項に画像処理方法。
【請求項７】
前記同心円状のパターンは、複素定数Ａと実定数ξ０を用いてＡｅｘｐ（２πｊξ０（ｘ２＋ｙ２））で表されることを特徴とする請求項６に記載の画像処理方法。
【請求項８】
前記符号化マスクの透過率分布は、ａ≦ｃを満たす実定数ａ，ｂ，ｃを用いて、ｃ＋ａ×ｃｏｓ（２πξ０（ｘ２＋ｙ２）＋ｂ）で表されることを特徴とする請求項７に記載の画像処理方法。
【請求項９】
前記逆畳み込み演算は、前記符号化データに対して、同心円状の複素パターンを乗じた後、フーリエ変換を行うステップを有することを特徴とする請求項６から８のいずれか1項に記載の画像処理方法。
【請求項１０】
前記同心円状の複素パターンは、Ａｅｘｐ（−２πｊξ０（ｘ２＋ｙ２））であることを特徴とする請求項９に記載の画像処理方法。
【請求項１１】
前記実関数Φ（ξ,η）は、空間周波数ξ，ηの少なくとも一方の２次関数を含むことを特徴とする請求項５に記載の画像処理方法。
【請求項１２】
前記実関数Φ（ξ,η）は、空間周波数ξ，ηの少なくとも一方の３次関数を含むことを特徴とする請求項５に記載の画像処理方法。
【請求項１３】
符号化マスクを用いた被写体の撮像に基づく符号化データを取得する取得手段と、
前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出する第１の算出手段と、
前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出する第２の算出手段と、
前記誤差画像を用いて前記推定画像を更新する更新手段と、を有することを特徴とする画像処理装置。
【請求項１４】
撮像素子と、
前記撮像素子に対して被写体側に配置された符号化マスクと、
前記撮像素子による前記符号化マスクを用いた被写体の撮像に基づく符号化データを取得する取得手段と、
前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出する第１の算出手段と、
前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出する第２の算出手段と、
前記誤差画像を用いて前記推定画像を更新する更新手段と、を有することを特徴とする撮像装置。
【請求項１５】
コンピュータに、
符号化マスクを用いた被写体の撮像に基づく符号化データを取得するステップと、
前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出するステップと、
前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出するステップと、
前記誤差画像を用いて前記推定画像を更新するステップと、を有することを特徴とする画像処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、符号化マスクを用いた被写体を撮像することで生成された符号化データに基づく演算処理によって被写体の画像を再構成する画像処理方法に関する。
【背景技術】
【０００２】
カメラや望遠鏡などの撮像装置は、レンズなどの光学素子を用いて被写体の像を撮像素子の撮像面に形成し、像の輝度分布を画像として取得する。撮像面上に形成された像は被写体の輝度分布をほぼそのまま再現しているため、撮像装置は特段の後処理なしに被写体の画像を取得することができる。しかしながら、理想的な画像を取得するためには高度に設計・加工されたレンズを使用する必要があり、そのようなレンズを備えた撮像装置は必然的に高価になってしまう。
【０００３】
特許文献１，２では、撮像装置の価格を抑えるために、レンズを用いることなく被写体の画像を再構成する方法が提案されている。上記方法では、撮像素子の前側（被写体側）に、特殊な形状の符号化マスクを配置し、取得された透過像から逆問題を計算機上で解くことで被写体の画像を再構成している。逆問題を解く方法として、特許文献１では、符号化マスクを、Ｔｏｅｐｌｉｔｚ行列を用いて分解し、ウィーナーフィルタおよびＬａｎｄｗｅｂｅｒ法を用いて再構成演算を行う方法が開示されている。また、特許文献２では、Ｔｉｋｈｏｎｏｖ正則化を用いた方法や圧縮検出を用いた方法などが開示されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】米国特許第９４４５１１５号明細書
【特許文献２】特表２０１６−５１０９１０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
通常、符号化マスクの周波数スペクトルには、振幅がゼロ近傍まで落ち込むゼロ点が存在する。そのため、周波数領域の情報を復元する場合、復元によってノイズが増幅されてしまい、良好な画像を算出することが困難である。一般的に、このようなノイズの増幅を抑えるために、ウィーナーフィルタを用いて回復処理が行われるが、回復処理はゼロ点を避けて行われるため、ゼロ点が存在する領域では実質的に回復効果が得られない。そのため、特許文献１の方法では、回復効果が得られない周波数成分が存在する。また、特許文献２の方法では、最適化問題を解くための計算負荷が大きく、画像を再構成するまでの時間が長くなってしまう。
【０００６】
本発明は、高速かつ高画質で被写体の画像を再構成可能な画像処理方法、画像処理装置、撮像装置および画像処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明の一側面としての画像処理方法は、符号化マスクを用いた被写体の撮像に基づく符号化データを取得するステップと、前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出するステップと、前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出するステップと、前記誤差画像を用いて前記推定画像を更新するステップと、を有することを特徴とする。
【０００８】
また、本発明の他の側面としての画像処理装置は、符号化マスクを用いた被写体の撮像に基づく符号化データを取得する取得手段と、前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出する第１の算出手段と、前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出する第２の算出手段と、前記誤差画像を用いて前記推定画像を更新する更新手段と、を有することを特徴とする。
【０００９】
また、本発明の他の側面としての撮像装置は、撮像素子と、前記撮像素子に対して被写体側に配置された符号化マスクと、前記撮像素子による前記符号化マスクを用いた被写体の撮像に基づく符号化データを取得する取得手段と、前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出する第１の算出手段と、前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出する第２の算出手段と、前記誤差画像を用いて前記推定画像を更新する更新手段と、を有することを特徴とする。
【００１０】
また、本発明の他の側面としての画像処理プログラムは、コンピュータに、符号化マスクを用いた被写体の撮像に基づく符号化データを取得するステップと、前記符号化データに対する、複素関数で表される複素関数マスクによる逆畳み込み演算を行うことで推定画像を算出するステップと、前記推定画像に対する前記複素関数マスクの共役マスクとの畳み込み演算を行うことで誤差画像を算出するステップと、前記誤差画像を用いて前記推定画像を更新するステップと、を有することを特徴とする。
【発明の効果】
【００１１】
本発明によれば、高速かつ高画質で被写体の画像を再構成可能な画像処理方法、画像処理装置、撮像装置および画像処理プログラムを提供することができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の実施形態に係る画像処理装置を有する撮像装置の概略図である。
【図２】被写体と符号化データとの関係図である。
【図３】本実施形態の画像処理方法を示すフローチャートである。
【図４】被写体を示す図である。
【図５】実施例１の符号化マスクを示す図である。
【図６】実施例１の符号化データを示す図である。
【図７】実施例１の再構成画像である。
【図８】実施例１の繰り返し回数に対応するＰＳＮＲを示す図である。
【図９】実施例２の符号化マスクを示す図である。
【図１０】実施例２の再構成画像を示す図である。
【図１１】実施例２の繰り返し回数に対応するＰＳＮＲを示す図である。
【図１２】実施例３の符号化マスクを示す図である。
【図１３】実施例３の再構成画像を示す図である。
【図１４】実施例３の繰り返し回数に対応するＰＳＮＲを示す図である。
【発明を実施するための形態】
【００１３】
以下、本発明の好ましい実施の形態を、添付の図面に基づいて詳細に説明する。各図において、同一の部材については同一の参照番号を付し、重複する説明は省略する。
【００１４】
図１は、本発明の実施形態に係る画像処理装置１０４を有する撮像装置１００の概略図である。撮像装置１００は、符号化マスク１０１、撮像素子１０２、データ保持装置１０３および画像処理装置１０４を有する。
【００１５】
被写体１から発せられた光は、符号化マスク１０１を透過して撮像素子１０２に到達する。撮像素子１０２は、符号化マスク１０１によって符号化された光強度分布（符号化データ）を取得する。データ保持装置１０３は、撮像素子１０２により取得された符号化データを保存する。画像処理装置１０４は、データ保持装置１０３から読み込んだ符号化データに対して、マスク形状に対応する画像処理を行う。
【００１６】
以下、画像処理装置１０４が実行する画像処理方法について説明する。まず、図２を参照して、本発明の原理を説明する。図２は、被写体１と符号化データとの関係図である。
【００１７】
被写体１上の点から発せられた光は符号化マスク１０１をほぼ一様に照明し、撮像素子１０２の撮像面上には符号化マスク１０１の影が投影される。被写体１上の光軸Ｏ上の点でない点１ａから発せられた光は、符号化マスク１０１を斜めに照明する。そのため、撮像面上には、符号化マスク１０１に対して、光軸Ｏと物点１ａから発せられた光とがなす角θ、および符号化マスク１０１と撮像素子１０２との間隔Ｄによって定まるシフト量δだけずれた位置に符号化マスク１０１の影２０１が投影される。撮像素子１０２は、各物点から発せられた光によって形成された符号化マスク１０１の投影像が足し合わされたものを符号化データとして取得する。被写体１の輝度分布をずれ量δの関数ｏ（δ）、符号化マスク１０１の透過率を光軸Ｏに垂直な方向の位置座標ｘの関数ｔ（ｘ）とすると、撮像素子１０２によって取得される符号化データｉ（ｘ）は以下の式（１）で表すことができる。
【００１８】
【数１】
2018206274.tif 000003
【００１９】
ここで、＊は畳み込み積分を意味する。画像処理装置１０４は、式（１）を満足する被写体１の輝度分布ｏ（δ）を符号化マスク１０１の透過率ｔ（ｘ）および符号化データｉ（ｘ）から算出する。しかしながら、一般にこの演算は困難とされている。式（１）の両辺をフーリエ変換すると、畳み込みの定理により式（１）は以下の式（２）で表される。
【００２０】
Ｉ（ξ）＝Ｏ（ξ）×Ｔ（ξ）                               （２）
Ｉ（ξ），Ｏ（ξ），Ｔ（ξ）はそれぞれ符号化データｉ（ｘ）、被写体１の輝度分布ｏ（δ）、符号化マスク１０１の透過率ｔ（ｘ）のフーリエ変換であり、ξは空間周波数を表す。式（２）より、Ｉ（ξ）をＴ（ξ）で除算することでＯ（ξ）を取得することができる。しかしながら、一般に符号化マスクの周波数スペクトルＴは特定の周波数で値が非常に小さくなるゼロ点を有するため、ゼロ割によってその周波数領域での情報が過剰に増幅されたり、ノイズが増幅されてしまったりする。そのため、従来の方法では、ウィーナーフィルタ等を用いることで周波数領域でのゼロ割を避けていた。しかしながら、ウィーナーフィルタを用いた場合、回復効果が得られない周波数領域が生じるため、十分な画質の再生像を取得することができない。
【００２１】
本実施形態では、このような問題を回避する方法について説明する。以下、任意の複素関数ｇ（ｘ）が任意の実関数Φ（ξ）を用いて以下の式（３），（４）で定義される場合について説明する。
【００２２】
Ｇ（ξ）＝ｅｘｐ（ｊΦ（ξ））                            （３）
ｇ（ｘ）＝ＦＴ−１［Ｇ（ξ）］                             （４）
ここで、ＦＴはフーリエ変換を意味し、ＦＴ−１は逆フーリエ変換を意味する。以下の説明では、Ｇ（ξ）を複素関数フィルタ、ｇ（ｘ）を複素関数マスクという。本実施形態では、複素関数マスクｇ（ｘ）の実部を用いて、以下の式（５）で表される透過率分布を有する符号化マスクが作成される。
【００２３】
ｔ（ｘ）＝（１＋Ｒｅ［ｇ（ｘ）］）／２                            （５）
本実施形態では、符号化データを取得するために透過率分布ｔ（ｘ）を使用するが、再構成演算を行う場合にはゼロ割を避けるため、式（５）を変形した以下の式（６）を使用する。
【００２４】
ｔ（ｘ）＝（１＋（ｇ（ｘ）＋ｇ（ｘ）＊）／２）／２         （６）
式（６）において、＊は複素共役を意味する。以下の説明では、ｇ（ｘ）＊を共役マスクという。式（１）に式（６）を代入し、変形することで、式（７）が得られる。
【００２５】
ｉ（ｘ）−Ｃ−ｏ（ｘ）＊ｇ（ｘ）＊／４＝ｏ（ｘ）＊ｇ（ｘ）／４    （７）
式（７）の左辺の第３項を書き換えると、式（８）が得られる。
【００２６】
ｉ（ｘ）−Ｃ−Ｎ（ｘ）＝ｏ（ｘ）＊ｇ（ｘ）／４             （８）
ここで、Ｃは座標に依らない定数であり、Ｎ（ｘ）は誤差画像であり、以下の式（９）で表される。
【００２７】
Ｎ（ｘ）＝ｏ（ｘ）＊ｇ（ｘ）＊／４                         （９）
式（８）を解くためには誤差画像Ｎ（ｘ）を知る必要があるが、式（９）より誤差画像Ｎ（ｘ）は被写体１の輝度分布ｏ（ｘ）によって決まるため、真の誤差画像Ｎ（ｘ）をあらかじめ知ることはできない。そこで、本実施形態では、まず、誤差画像Ｎ（ｘ）を０として暫定的な推定画像ｏ０（ｘ）を算出する。式（８）において誤差画像Ｎ（ｘ）を０とすると、式（８）は以下の式（１０）で表される。
【００２８】
ｉ（ｘ）−Ｃ＝ｏ０（ｘ）＊ｇ（ｘ）／４                    （１０）
式（３）より複素関数マスクｇ（ｘ）は空間スペクトルにゼロ点を持たないため、式（１０）で符号化データｉ（ｘ）に対して逆畳み込み演算を行うことで、推定画像ｏ０（ｘ）を容易に取得することができる。式（９）より、推定画像ｏ０（ｘ）と共役マスクｇ（ｘ）＊との畳み込み演算によって新たな誤差画像Ｎ０（ｘ）を算出することができる。式（８）の誤差画像Ｎ（ｘ）を算出された誤差画像Ｎ０（ｘ）で置き換えれば、逆畳み込み演算によって新たな推定画像ｏ１（ｘ）を算出することができる。そして、式（９）を用いて推定画像ｏ１（ｘ）に基づいて新たな誤差画像Ｎ１（ｘ）を算出することができる。
【００２９】
以上のように、誤差画像Ｎ（ｘ）および推定画像ｏ（ｘ）を交互に推定することで、より高画質な推定画像を算出することができる。さらに、畳み込み演算と逆畳み込み演算は畳み込みの定理から、ＦＦＴ（高速フーリエ変換）を用いて高速に処理が可能である。特に、複素関数マスクｇ（ｘ）との逆畳み込み演算は、ＦＦＴを用いて画像をスペクトルに変換した後、複素関数フィルタＧ（ξ）で除算し、逆ＦＦＴを行うことで高速に実行可能である。
【００３０】
なお、本実施形態では、簡単のため１次元の場合について説明したが、２次元の場合にも容易に拡張することができる。
【００３１】
以下、図３を参照して、画像処理装置１０４が実行する画像処理方法について説明する。図３は、本実施形態の画像処理方法を示すフローチャートである。本実施形態の画像処理方法は、ソフトウエアおよびハードウエア上で動作する画像処理プログラムにしたがって実行される。画像処理プログラムは、例えば、画像処理装置１０４内に格納されていてもよいし、コンピュータが読み取り可能な記録媒体に記録されていてもよい。また、本実施形態では画像処理装置１０４が画像処理方法を実行するが、パーソナルコンピュータ（ＰＣ）や専用の装置が画像処理装置として本実施形態の画像処理方法を実行してもよい。また、本実施形態の画像処理プログラムに対応する回路を設け、回路を動作させることで本実施形態の画像処理方法を実行してもよい。
【００３２】
ステップＳ１では、画像処理装置１０４は、符号化マスクを用いた被写体１の撮像に基づく符号化データを取得する。すなわち、画像処理装置１０４は、取得手段として機能する。
【００３３】
ステップＳ２では、画像処理装置１０４は、符号化データからオフセットを減算する。
【００３４】
ステップＳ３では、画像処理装置１０４は、符号化データに対して、逆畳み込み演算を行うことで推定画像を算出する。すなわち、画像処理装置１０４は、第１の算出手段として機能する。具体的には、まず、画像処理装置１０４は、符号化データに対してＦＦＴを実行し、変換データ（スペクトルデータ）を算出する。次に、画像処理装置１０４は、変換データに対して、複素関数フィルタＧ（ξ）で除算した後、逆ＦＦＴを実行することで、推定画像を算出する。
【００３５】
ステップＳ４では、画像処理装置１０４は、繰り返し回数が所定の回数に到達したかどうかを判定する。繰り返し回数が所定の回数より多い場合、処理を終了し、推定画像を被写体の画像として取得する。また、繰り返し回数が所定の回数より少ない場合、ステップＳ５に進む。なお、繰り返し回数が所定の回数と等しい場合にどちらのステップに進むかは任意である。
【００３６】
ステップＳ５では、画像処理装置１０４は、式（９）より、誤差画像を算出する。すなわち、画像処理装置１０４は、第２の算出手段として機能する。
【００３７】
ステップＳ６では、画像処理装置１０４は、符号化データから誤差画像を減算することで減算画像を算出する。その後、ステップＳ３に戻り、画像処理装置１０４は、減算画像に対して、逆畳み込み演算を行うことで新たに推定画像を算出する。すなわち、画像処理装置１０４は、更新手段として機能する。
【００３８】
以上の方法によって、符号化データから被写体の画像を高速かつ高画質で再構成することができる。
【００３９】
以下の各実施例では、被写体として図４に示される画像を使用して再構成画像を取得する場合について説明する。
【実施例１】
【００４０】
本実施例では、符号化マスク１０１として、図５に示される符号化マスクを使用する。本実施例の符号化マスクは、同心円状のパターンを有する。本実施例では、実関数Φ（ξ，η）を−π（ξ２＋η２）／２ξ０とし、式（３）〜（５）より符号化マスク１０１の透過率分布ｔ（ｘ）は１＋ｃｏｓ（２πξ０（ｘ２＋ｙ２））／２と算出される。
【００４１】
まず、画像処理装置１０４は、符号化マスク１０１を用いて生成された図６に示される符号化データをデータ保持装置１０３から読み込む。図６に示される符号化データは、ξ０＝０．０１７［１／ｐｉｘｅｌ］、サンプリング数として３００×３００を用いてシミュレーションによって生成される。
【００４２】
次に、画像処理装置１０４は、図３のフローチャートに従って被写体の画像を算出する。本実施例では、繰り返し回数の上限を５０回としている。図７は、繰り返し回数ｋが１回、１１回、２１回の場合に取得された推定画像を示している。図８は、繰り返し回数に対応する、画像評価指標の１つであるＰＳＮＲを示している。図７に示されるように、繰り返し回数が増加することで、画質が向上する。また、図８に示されるように、繰り返し回数が増加することで、ＰＳＮＲが増加する。
【実施例２】
【００４３】
本実施例では、符号化マスク１０１として、図９に示される符号化マスクを使用する。本実施例では、実関数Φ（ξ、η）を−π（ξ２＋０．５η２）／２ξ０とし、符号化マスク１０１の透過率分布ｔ（ｘ）は式（３）〜（５）を用いて算出される。実関数Φは楕円形状を表すため、符号化マスク１０１は楕円形状を有する。
【００４４】
画像処理装置１０４は、符号化マスク１０１を用いて実施例１と同様のシミュレーションにより生成された符号化データを取得した後、図３のフローチャートに従って被写体の画像を算出する。本実施例では、繰り返し回数の上限を５０回としている。図１０は、繰り返し回数ｋが１回、１１回、２１回の場合に取得された推定画像を示している。図１１は、繰り返し回数に対応するＰＳＮＲを示している。図１０に示されるように、繰り返し回数が増加することで、画質が向上する。また、図１１に示されるように、繰り返し回数が増加することで、ＰＳＮＲが増加する。
【実施例３】
【００４５】
本実施例では、符号化マスク１０１として、図１２に示される符号化マスクを使用する。本実施例では、実関数Φ（ξ、η）を−π（｜ξ｜３＋｜η｜３）／ξ０で表される３次関数とし、符号化マスク１０１の透過率分布ｔ（ｘ）は式（３）〜（５）を用いて算出される。なお、本実施例では、実関数Φ（ξ，η）は、３次関数であるが、２次関数などのｎ次関数であってもよい。
【００４６】
画像処理装置１０４は、符号化マスク１０１を用いて実施例１と同様のシミュレーションにより生成された符号化データを取得した後、図３のフローチャートに従って被写体の画像を算出する。本実施例では、繰り返し回数の上限を５０回としている。図１３は、繰り返し回数ｋが１回、１１回、２１回の場合に取得された推定画像を示している。図１４は、繰り返し回数に対応するＰＳＮＲを示している。図１３に示されるように、繰り返し回数が増加することで、画質が向上する。また、図１４に示されるように、繰り返し回数が増加することで、ＰＳＮＲが増加する。
【実施例４】
【００４７】
被写体１の輝度分布ｏ（ｘ）を算出する際に用いる逆畳み込み処理は、演算の一部をあらかじめ解析的に実行しておくことで更に高速な処理を行うことも可能である。
【００４８】
本実施例では、実関数Φ（ξ）を−πξ２／２ξ０とする。簡単のため１次元での結果を示す。式（３）、式（４）および式（１０）より逆畳み込み演算によって算出される推定画像ｏ０（ｘ）は、畳み込みの定理により以下の式（１１）で表される。
【００４９】
【数２】
2018206274.tif 000004
【００５０】
ここで、Ｉ‘（ξ）は符号化データのフーリエ変換であり、Ｃ’は複素定数である。
【００５１】
ｅｘｐ（４πｊξ０ｘｘ‘）を乗じて積分を行うことはフーリエ変換を行うことと等価であるため、この畳み込み演算は１回のＦＦＴで算出することができる。数値的に逆畳み込み演算を行うためには、少なくとも画像のスペクトルを算出する場合と、フィルタが乗じられたスペクトルを逆フーリエ変換する場合の最低２回のＦＦＴを要する。しかしながら、本実施例で示したように、演算の一部をあらかじめ解析的に実行しておくことで、複素パターンを乗じた後、１回のＦＦＴで画像を再構成することが可能となり、演算の高速化を実現することができる。
［他の実施形態］
本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサーがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。
【００５２】
以上、本発明の好ましい実施形態について説明したが、本発明はこれらの実施形態に限定されず、その要旨の範囲内で種々の変形及び変更が可能である。
【００５３】
本実施形態では、誤差画像の初期値を０としているが、本発明はこれに限定されない。初期値を任意の値に設定し、図３のフローチャートにおけるステップＳ３の前に減算しておいてもよい。
【００５４】
本実施形態では、式（３）において、複素関数フィルタＧ（ξ）として全周波数領域で絶対値が１となる複素関数を用いているが、本発明はこれに限定されない。例えば、ａ（ξ、η）×ｅｘｐ（ｊΦ（ξ、η））で表される、絶対値に減衰や変調を有する関数であってもよい。特に、絶対値を表す関数ａ（ξ、η）がゼロ点を持たなければ、高い回復効果を得ることができる。
【００５５】
本実施形態では、空間周波数に対して特段の限定を行っていないが、実際に作成できるマスクには加工できる最小サイズが決まっており、これにより符号化マスクが有することが可能な最大周波数が決まってしまう。したがって、厳密に言えば、全周波数領域で絶対値が１となる複素関数フィルタは存在しえない。しかしながら、本発明を実施する場合、対象とする周波数領域内で関数ａ（ξ、η）にゼロ点が存在しなければ、本発明の効果を得ることができる。例えば、以下の式（１２）で表される、符号化データのサンプリングピッチ（サンプリング間隔）Δの逆数から定まる周波数よりも小さな周波数領域内でゼロ点がなければよい。
【００５６】
ξｍａｘ＝１／２Δ                                  （１２）
本実施形態では、周波数スペクトルの値がゼロとなる箇所をゼロ点と称してきたが、ゼロ点として取り扱われるべき箇所は周波数スペクトルの値が厳密にゼロとなる箇所のみではない。特に、数値計算でデータを処理する場合、厳密なゼロ点というのは現れないことが多い。しかしながら、このような場合においても、他の領域に比べて周波数スペクトルの値が小さい領域は存在しており、その領域がゼロ点として称される箇所となる。例えば、ＦＦＴなどの数値計算で得られる周波数スペクトルの絶対値の最大値に対して、周波数スペクトルの値が十分に小さい点および領域はゼロ点として扱ってもよい。具体的には、周波数スペクトルの絶対値の最大値に対して、周波数スペクトルの値が０．１５倍以下の点および領域はゼロ点として扱ってもよい。また、式（５）で定義される符号化マスクの透過率には、負となる箇所を避けるために定数によるオフセットが足されている。このような符号化マスクの場合、周波数スペクトルの原点を除く領域での絶対値の最大値に対して、周波数スペクトルの値が０．１５倍以下の点および領域をゼロ点として扱えばよい。
【００５７】
本実施形態では、式（５）を用いて、複素関数マスクｇ（ｘ）から符号化マスクを生成する方法を説明したが、本発明はこれに限定されない。例えば、複素関数マスクｇ（ｘ）の振幅が大きい場合、規格化などを施した後に式（５）を用いて符号化マスクを生成してもよいし、式（５）において加える定数を１とは異なる値で調整してもよい。また、複素関数マスクｇ（ｘ）の実部ではなく虚部から符号化マスクを生成してもよい。更に言えば、符号化マスクの形状は、複素関数ｇ（ｘ）の実部または虚部を陽に含む形状である必要はない。符号化マスクは、複素関数ｇ（ｘ）の実部および虚部から新たに算出される別の関数を透過率分布として有するように生成されてもよい。例えば、実部および虚部を重み付き平均することで得られる新たな関数を透過率分布として有するように生成されてもよい。この場合、式（７）は、以下の式（１３）に変更される。
【００５８】
ｉ（ｘ）−Ｃ−Ｃ１×ｏ（ｘ）＊ｇ（ｘ）＊＝Ｃ２×ｏ（ｘ）＊ｇ（ｘ）   （１３）
Ｃ１、Ｃ２は先の重みによって定まる複素係数である。式（１３）を用いることで、本実施形態で説明した効果を得ることができる。
【００５９】
本発明では、周波数スペクトルを算出する際、または周波数スペクトルから実空間でのデータを算出する際、フーリエ変換と逆フーリエ変換のどちらを用いてもよい。すべての演算において、一貫した組み合わせになっていればよい。変換方法の違いによって本発明の趣旨は変わるものではない。
【符号の説明】
【００６０】
１０１  符号化マスク
【図１】
2018206274.tif 000005
【図２】
2018206274.tif 000006
【図３】
2018206274.tif 000007
【図４】
2018206274.tif 000008
【図５】
2018206274.tif 000009
【図６】
2018206274.tif 000010
【図７】
2018206274.tif 000011
【図８】
2018206274.tif 000012
【図９】
2018206274.tif 000013
【図１０】
2018206274.tif 000014
【図１１】
2018206274.tif 000015
【図１２】
2018206274.tif 000016
【図１３】
2018206274.tif 000017
【図１４】
2018206274.tif 000018