(57)【要約】【課題】波面伝達媒質を通して観測される少数の観測波面を用いて、標本から波面伝達媒質に向かう未知の波面を再構成する。【解決手段】照明光源１０１からの波面を整形して標本に向けて出射する空間光変調器１０２と、内部で光を複数回、散乱または反射させることにより入射する平面波の角度に応じた波面を出射する特性を有し、標本からの波面を伝達する波面伝達媒質１０３と、波面伝達媒質から出射した波面を観測する検出器１０４と、該検出器により観測された波面である観測波面を用いて標本から波面伝達媒質に向かう波面である標本波面を再構成する演算を行う波面演算装置１０５とを含む撮像システムにおいて、波面演算装置は、１つの観測波面と、該観測波面の推定の誤差に関する誤差項と、標本波面の推定の振幅に関する正則化項とを含む目的関数を用意し、該目的関数を最小化することで標本波面を再構成する。【選択図】図１
2016057172.tif 000002
【特許請求の範囲】
【請求項１】
照明光源からの波面を整形して標本に向けて出射する空間光変調器と、
内部で光を複数回、散乱または反射させることにより、入射する平面波の角度に応じた波面を出射する特性を有し、前記標本からの波面を伝達する波面伝達媒質と、
前記波面伝達媒質から出射した波面を観測する検出器と、
前記検出器により観測された波面である観測波面を用いて、前記標本から前記波面伝達媒質に向かう波面である標本波面を再構成する演算を行う波面演算装置とを含む撮像システムに用いられる前記波面演算装置であって、
１つの前記観測波面と、該観測波面の推定の誤差に関する誤差項と、前記標本波面の推定の振幅に関する正則化項とを含む目的関数を用意し、
前記目的関数を最小化することで前記標本波面を再構成することを特徴とする波面演算装置。
【請求項２】
前記正則化項は、前記標本波面の推定の振幅のTotal Variationであり、
前記誤差項および前記正則化項を勾配法によって交互に最小化することにより、前記標本波面を再構成することを特徴とする請求項１に記載の波面演算装置。
【請求項３】
前記標本波面の推定の振幅を２次元配列から辞書式順序により１次元ベクトルに変換した場合において、該推定の振幅のｉ番目の画素値をｘｉとし、該推定の振幅の画像サイズをＭ×Ｎ画素とするとき、前記Total Variationは、
2016057172.tif 000015
なる式で表されることを特徴とする請求項２に記載の波面演算装置。
【請求項４】
前記観測波面および前記誤差項を勾配法により最小化した後に前記標本波面の振幅を強閾値処理することで前記正則化項を最小化する処理を繰り返すことにより、前記標本波面を再構成することを特徴とする請求項１に記載の波面演算装置。
【請求項５】
前記強閾値処理は、
前記標本波面の振幅をスパース変換してスパース係数を得るステップと、
前記スパース係数に対する閾値処理を行うステップと、
前記閾値処理が行われたスパース係数を逆スパース変換するステップとを含むことを特徴とする請求項４に記載の波面演算装置。
【請求項６】
照明光源からの波面を整形して標本に向けて出射する空間光変調器と、
内部で光を複数回、散乱または反射させることにより、入射する平面波の角度に応じた波面を出射する特性を有し、前記標本からの波面を伝達する波面伝達媒質と、
前記波面伝達媒質から出射した波面を観測する検出器と、
前記検出器により観測された波面である観測波面を用いて、前記標本から前記波面伝達媒質に向かう波面である標本波面を再構成する演算を行う請求項１から５のいずれか一項に記載の波面演算装置とを有することを特徴とする撮像システム。
【請求項７】
照明光源からの波面を整形して標本に向けて出射する空間光変調器と、
内部で光を複数回、散乱または反射させることにより、入射する平面波の角度に応じた波面を出射する特性を有し、前記標本からの波面を伝達する波面伝達媒質と、
前記波面伝達媒質から出射した波面を観測する検出器と、
前記検出器により観測された波面である観測波面を用いて、前記標本から前記波面伝達媒質に向かう波面である標本波面を再構成する演算を行う波面演算装置とを含む撮像システムに用いられる前記波面演算装置としてのコンピュータに前記演算を行わせるコンピュータプログラムであって、
前記コンピュータに、
１つの前記観測波面と、該観測波面の推定の誤差に関する誤差項と、前記標本波面の推定の振幅に関する正則化項とを含む目的関数を用意させ、
前記目的関数を最小化することで前記標本波面を再構成させることを特徴とする波面演算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、標本からの波面を散乱媒質等の波面伝達が可能な媒質（以下、波面伝達媒質という）を通して観測することで得られた観測波面を用いて、標本からの波面を再構成する技術に関する。
【背景技術】
【０００２】
波面伝達媒質を通った観測波面を用いて標本からの波面を再構成することが可能な方法として、非特許文献１および非特許文献２に開示されたものがある。
【０００３】
非特許文献１には、波面伝達媒質の透過行列を用いて、波面伝達媒質を介して観測される１つの出射波面から未知の入射波面を再構成する方法が開示されている。透過行列とは、波面伝達媒質の伝達特性を表す行列であり、波面伝達媒質への入射波面と波面伝達媒質からの出射波面とを関係付ける。なお、再構成は観測された出射波面に対して、透過行列の（擬似）逆行列をかけることで行われる。この方法は、濁レンズイメージング（Turbid Lens Imaging：ＴＬＩ）と呼ばれる。
【０００４】
非特許文献２には、非特許文献１と同様にＴＬＩを用いた方法として、観測された複数の出射波面のそれぞれに対してＴＬＩによる再構成を行い、得られた結果を足し合わせる方法が開示されている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Youngwoon Choi, Taeseok Daniel Yang, Christopher Fang-Yen, Pilsung Kang, Kyoung Jin Lee, Ramachandra R. Dasari, Michael S. Feld, Wonshik Choi, "Overcoming the Diffraction Limit Using Multiple Light Scattering in a Highly Disordered Medium", Physical Review Letter, American Physical Society, 2011, Vol.107, Issue2
【非特許文献２】Youngwoon Choi, Changhyeong Yoon, Moonseok Kim, Taeseok Daniel Yang, Christopher Fang-Yen, Ramachandra R. Dasari, Kyoung Jin Lee, Wonshik Choi, "Scanner-Free and Wide-Field Endoscopic Imaging by Using a Single Multimode Optical Fiber", Physical Review Letter, American Physical Society, 2012, Vol.109, p.203901-1-203901-5
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、非特許文献１にて開示されたように単にＴＬＩを用いる方法では、イメージセンサによる観測過程で生じるノイズにロバストな再構成を行うことができないという問題がある。これに対して、非特許文献２にて開示された方法では、複数の観測波面（出射波面）のそれぞれに対するＴＬＩによる再構成結果を足し合わせることで、観測過程で生じるノイズにロバストな再構成を行うことができる。しかしながら、この方法では極めて多くの回数（例えば５００回）もの出射波面の観測が必要であり、全観測に必要な時間も長くなる。しかも、複数の出射波面の取得には、ガルバノミラーが必要であり、装置の構成が複雑になる。
【０００７】
本発明は、波面伝達媒質を通して観測される少数の観測波面を用いて、標本から波面伝達媒質に向かう未知の波面を、観測過程で生じるノイズにロバストに再構成することができるようにした波面演算装置を提供する。
【課題を解決するための手段】
【０００８】
本発明の一側面としての波面演算装置は、照明光源からの波面を整形して標本に向けて出射する空間光変調器と、内部で光を複数回、散乱または反射させることにより、入射する平面波の角度に応じた波面を出射する特性を有し、標本からの波面を伝達する波面伝達媒質と、波面伝達媒質から出射した波面を観測する検出器と、該検出器により観測された波面である観測波面を用いて、標本から波面伝達媒質に向かう波面である標本波面を再構成する演算を行う波面演算装置とを含む撮像システムに用いられる。該波面演算装置は、１つの観測波面と、該観測波面の推定の誤差に関する誤差項と、標本波面の推定の振幅に関する正則化項とを含む目的関数を用意し、該目的関数を最小化することで標本波面を再構成することを特徴とする。
【０００９】
また、本発明の他の一側面としての波面演算プログラムは、照明光源からの波面を整形して標本に向けて出射する空間光変調器と、内部で光を複数回、散乱または反射させることにより、入射する平面波の角度に応じた波面を出射する特性を有し、標本からの波面を伝達する波面伝達媒質と、波面伝達媒質から出射した波面を観測する検出器と、該検出器により観測された波面である観測波面を用いて、標本から波面伝達媒質に向かう波面である標本波面を再構成する演算を行う波面演算装置とを含む撮像システムに用いられる波面演算装置としてのコンピュータに上記演算を行わせるコンピュータプログラムである。該プログラムは、コンピュータに、１つの観測波面と、該観測波面の推定の誤差に関する誤差項と、標本波面の推定の振幅に関する正則化項とを含む目的関数を用意させ、該目的関数を最小化することで標本波面を再構成させることを特徴とする。
【発明の効果】
【００１０】
本発明によれば、波面伝達媒質を通して検出器により観測される１つの観測波面を用いて、標本から波面伝達媒質に向かう未知の標本波面を、検出器による観測過程で生じるノイズにロバストに再構成することができる。
【図面の簡単な説明】
【００１１】
【図１】本発明の基本実施例である波面演算装置の構成を示すブロック図。
【図２】基本実施例の波面演算方法の手順を示すフローチャート。
【図３】実施例１の波面演算方法の手順を示すフローチャート。
【図４】実施例１による波面演算結果を示す図。
【図５】実施例２による波面演算結果を示す図。
【図６】従来の波面演算結果を示す図。
【発明を実施するための形態】
【００１２】
以下、本発明の実施例について図面を参照しながら説明する。
【００１３】
まず、具体的な実施例の説明に先立って、本発明の基本的な実施例である撮像システムについて説明する。図１には、該撮像システムの構成を示している。撮像システム１００は、照明光源１０１と、空間光変調器１０２と、波面伝達媒質（図には擬ランダム媒質と記す）１０３と、検出器１０４と、波面演算装置１０５とを有する。照明光源１０１、空間光変調器１０２および検出器１０４と波面演算装置１０５とは、バス配線１０６によって接続されている。図中の破線は、上記各構成要素間を伝搬する波面を模式的に示している。
【００１４】
また、空間光変調器１０２と検出器１０４との間には、標本が配置されている。標本は、平面波で照明された際に、該標本の情報を含んだ波面を透過または反射させるものあればどのようなものでもよい。
【００１５】
照明光源１０１は、例えばＨｅ−Ｎｅレーザーにより構成され、標本を照明するための照明光を発する。
【００１６】
空間光変調器１０２は、例えばＬＣＯＳ（Liquid Crystal On Silicon）により構成され、照明光源１０１からの波面を整形する。ＬＣＯＳは、シリコン基板上の液晶を電圧で駆動し、入射した波面（位相）を任意の波面に整形するデバイスである。なお、空間光変調器１０２は、ＬＣＯＳに限らず、照明光源１０１からの入射波面を整形して標本を照明する平面波を出射させるものであればどのようなものでもよい。また、前述した照明光源１０１もＨｅ−Ｎｅレーザーに限らず、空間光変調器１０２と組み合わせて標本を平面波で照明できる光源であればどのようなものでもよい。
【００１７】
なお、図１に示したように照明光源１０１と空間光変調器１０２を波面伝達媒質１０３に対して標本側に配置するのではなく、検出器１０４側に設置し、波面伝達媒質１０３を通して標本を照明する構成も採ることができる。この構成でも、標本からの波面は、波面伝達媒質１０３を通して検出器１０４により観測される。この場合は、波面伝達媒質１０３の伝達特性まで考慮して標本を平面波で照明できるように、空間光変調器１０２で照明光源１０１からの入射波面を適切に整形すればよい。
【００１８】
波面伝達媒質１０３は、例えばマルチモードファイバで構成され、標本からの波面を標本側から検出器１０４側（以下、観測側ともいう）に伝達する。本実施例で用いる波面伝達媒質１０３は、内部で光を複数回、散乱または反射させることにより、入射する平面波の角度に応じた波面を出射する特性を有する媒質であり、マルチモードファイバ以外には、拡散板等も含まれる。この媒質に平面波を入射させると、一見ランダムなスペックルパターンが得られるため、この媒質を図１中に記載したように「擬ランダム媒質」ということもできる。ここにいう擬ランダム媒質は、論文等では一般的に散乱媒質（scatter mediumまたはturbid medium）と呼ばれているが、マルチモードファイバと散乱とは容易に結び付かないため、擬ランダム媒質との呼称がより適切である。
【００１９】
マルチモードファイバは、屈折率が異なるコアとクラッドと呼ばれる光学部材を用い、コアの外周を覆うようにクラッドを配置した光伝送路である。コアの内部で全反射された光は、複数のモードを有してコアの内部を伝送される。なお、マルチモードファイバも、入射する平面波の入射角度に応じた波面を出射することが非特許文献２で開示されており、擬ランダム媒質の１つである。また、波面伝達媒質（擬ランダム媒質）１０３は、マルチモードファイバに限らず、上述した特性を有する媒質であればどのようなものでもよい。
【００２０】
検出器１０４は、例えば干渉計で構成され、標本からの波面を波面伝達媒質１０３を通して観測する。干渉計は、未知の入射波面と既知の参照波面とを干渉させて得られる強度分布から未知の入射波面を推定するための装置である。なお、検出器１０４は干渉計に限らず、標本からの波面を波面伝達媒質１０３を通して観測できるものであればどのようなものでもよい。
【００２１】
波面演算装置１０５は、パーソナルコンピュータやマイクロコンピュータ等のコンピュータにより構成され。バス配線１０６を介して照明光源１０１、空間光変調器１０２および検出器１０４と相互に通信可能である。波面演算装置１０５は、波面伝達媒質１０３を通して検出器１０４で観測された標本からの波面（既知の波面：以下、観測波面という）を入力とし、数値演算によって標本から波面伝達媒質１０３に向かう波面（未知の波面：以下、標本波面という）を再構成する。波面演算装置１０５は、再構成の結果を出力または保存する。
【００２２】
なお、波面演算装置１０５に入力される観測波面は、検出器１０４から直接入力される観測波面だけではなく、検出器１０４により観測されて半導体メモリや光ディスク等の記憶媒体に保存され、該記憶媒体から読み出された観測波面であってもよい。
【００２３】
波面演算装置１０５は、図示はしないが、制御部であるＣＰＵや、記憶部であるＲＡＭ、入力デバイスであるキーボードやマウス、出力デバイスであるディスプレイやＵＳＢインターフェース等を有する。
【００２４】
このように構成される撮像システムの具体例として、内視鏡が挙げられる。内視鏡では、波面伝達媒質１０３としてマルチモードファイバを用い、該マルチモードファイバよりも観測側に照明光源１０１と空間光変調器１０２を設置する。そして、照明光源１０１からの波面をマルチモードファイバの伝達特性を考慮して空間光変調器１０２で整形した後、マルチモードファイバを通して体内の標本（例えば、血管や乳管の内壁）を照明する。次に、マルチモードファイバを通して観測波面を検出器１０４により観測する。最後に、観測波面を波面演算装置１０５に入力して標本波面を再構成する。これにより、従来では観察することが困難であった、狭い空間内部の標本を観察することができる。
【００２５】
次に、波面演算装置１０５が、観測波面から標本波面を再構成する処理（波面演算方法）を、図２に示すフローチャートを用いて説明する。コンピュータである波面演算装置１０５は、コンピュータプログラムである波面演算プログラムに従って以下の処理を行う。このことは、後述する具体的な実施例１，２においても同じである。
【００２６】
なお、以下の説明では、照明光源１０１の波長を単波長であるとしているが、多波長の場合は、以下の処理を波長毎に行えばよい。多波長の例としては、複数の単波長の光源を含む照明光源１０１からの照明光によって標本を照明するという構成が考えられる。また多波長の照明光源の構成はこれに限らず、多波長で標本が照明できればどのようなものでもよい。
【００２７】
また、以下では、１つの観測波面を用いる場合について説明するが、これは標本波面の再構成には１つの観測波面があれば十分であることを意味する一方、複数の観測波面を用いることを制限する意味ではない。例えば、観測波面が複数ある場合には、以下の処理を観測波面毎に行えばよい。
【００２８】
ステップＳ２０１では、波面演算装置１０５は、透過行列を準備する。透過行列とは、波面伝達媒質１０３にその標本側から様々な入射角度で平面波を入射させ、波面伝達媒質１０３を通して観測側で観測される波面を並べた行列である。波面を観測する検出器１０４の総画素数をＭとし、入射する平面波の入射角度数をＮとすると、透過行列のサイズはＭ×Ｎである。波面は振幅と位相とを有する複素数であるので、透過行列の要素は複素数となる。また、透過行列は観測の前に予め取得しておく。また、入射する平面波の波長は、照明光源１０１の波長と同じであるとする。また、マルチモードファイバのように波面伝達媒質１０３が変形する場合は、変形に伴って透過行列が変化するので、波面伝達媒質１０３の変形と関連付けて複数の透過行列を取得しておく。そして、観察の際には、複数の透過行列のうち適切な透過行列を用いればよい。また、ユーザが適切な透過行列を波面演算装置１０５に入力することで、このステップに代えてもよい。
【００２９】
ステップＳ２０２では、波面演算装置１０５は、ステップＳ２０１で準備した透過行列と、波面伝達媒質１０３を通して観測される観測波面と、標本波面の推定の角度スペクトルとを含む目的関数を用意する。そして、この目的関数を標本波面の推定の角度スペクトルに関して最小化することで標本波面を再構成する。目的関数は、例えば以下の式（１）で与えられる。
【００３０】
【数１】
2016057172.tif 000003
【００３１】
ただし、ｙは波面伝達媒質１０３を通して観測される観測波面であり、ｘは標本波面の推定の角度スペクトルであり、Φは透過行列である。
【００３２】
式（１）における第一項Ψ１は、波面伝達媒質１０３を通して観測される観測波面とその観測波面の推定との誤差を表す誤差項である。第二項Ψ２は、標本波面の推定の尤もらしさを表す正則化項である。μは、これら２つの項の寄与を調整する正則化パラメータである。
【００３３】
また、透過行列ΦのサイズをＭ×Ｎとすると、波面伝達媒質１０３を通して観測される観測波面ｙのサイズはＭ×１のベクトルであり、標本波面の推定の角度スペクトルｘのサイズはＮ×１のベクトルである。ただし、観測波面は検出器１０４を構成する２次元アレイセンサを用いて取得されるので、そのサイズは通常２次元配列となるが、ここでは２次元配列を辞書式順序で並び替えた１次元ベクトルとして表現している。２次元配列を１次元ベクトルへ変換するルールとして、辞書式順序での並び替え以外のルールを用いてもよい。また、標本波面の角度スペクトルに関しても同様である。
【００３４】
標本波面とその角度スペクトルの関係は、以下の式（２）で与えられる。標本波面が推定であったとしても、以下の関係は変わらない。
【００３５】
【数２】
2016057172.tif 000004
【００３６】
ただし、ｘは前述した標本波面の推定の角度スペクトルであり、Ｘは標本波面であり、ＦＴはフーリエ変換を表す演算子である。
【００３７】
また、式（１）中の観測波面とその推定との誤差は、例えばＬ２ノルムで計算される。なお、Ｌ２ノルムは以下の式（３）で与えられる。
【００３８】
【数３】
2016057172.tif 000005
【００３９】
ただし、ａ，ｂはサイズＭ×１の任意のベクトルである。なお、本実施例では誤差関数としてＬ２ノルムを用いるが、誤差を計算できる関数であればどのような誤差関数を用いてもよい。
【００４０】
以下、より具体的な実施例について説明する。
【実施例１】
【００４１】
実施例１では、ＴＶ正則化（TV Regularization）と呼ばれる画像回復手法を用いて観測波面から標本波面を再構成する処理を、図３に示すフローチャートを用いて説明する。なお、ＴＶはTotal Variation（全変動）の略である。
【００４２】
ステップＳ３０１では、波面演算装置１０５は、ステップＳ２０１と同様に、透過行列を準備する。
【００４３】
ステップＳ３０２では、波面演算装置１０５は初期化を行う。具体的には、標本波面の推定の角度スペクトルを設定する。本実施例では、一様乱数を用いて初期化を行うが、他の初期化方法を用いてもよい。
【００４４】
ステップＳ３０３では、波面演算装置１０５は、以下の式（４）に従って目的関数の誤差項を最小化し、標本波面の推定の角度スペクトルを更新する。なお、最小化は一度に行わず、少しずつ推定解を更新していく手順を反復する逐次近似（勾配法）によって行う。
【００４５】
【数４】
2016057172.tif 000006
【００４６】
ただし、∇は勾配を表す演算子であり、αは１回の反復における更新の大きさを表すステップ幅である。←は右辺の計算結果を左辺の変数に代入する記号である。他の記号は式（１）と同じである。
【００４７】
ステップＳ３０４では、波面演算装置１０５は、標本波面の推定の角度スペクトルから、標本波面の推定の振幅と位相を、以下の式（５）を用いて求める。
【００４８】
【数５】
2016057172.tif 000007
【００４９】
ただし、Ｘａｍｐは標本波面の推定の振幅であり、Ｘｐｈａｓｅは標本波面の推定の位相である。ＩＦＴは逆フーリエ変換であり、ａｂｓは絶対値をとる演算子である。ｉｍａｇは波面Ｘの虚部をとる演算子であり、ｒｅａｌは波面Ｘの実部をとる演算子である。．／は要素毎の除算である。その他の記号は式（１），（２）と同である。
【００５０】
ステップＳ３０５では、波面演算装置１０５は、以下の式（６）に従って目的関数の正則化項を最小化し、標本波面の推定の振幅を更新する。本ステップでも、ステップＳ３０２同様に最小化を逐次近似（勾配法）により行う。
【００５１】
【数６】
2016057172.tif 000008
【００５２】
ただし、ＴＶはＴＶ正則化項であり、∇は勾配を計算する演算子である。βは１回の反復における更新の大きさを表すステップ幅である。他の記号は式（４），（５）と同じである。ＴＶ正則化項は、以下の式（７）で与えられる。
【００５３】
【数７】
2016057172.tif 000009
【００５４】
ただし、ｃは画素サイズがＭ×Ｎ画素の２次元配列を辞書式順序によって１次元ベクトルに変換したときの要素であり、ｃｉはｃのｉ番目の要素である。
【００５５】
ステップＳ３０６では、波面演算装置１０５は、以下の式（８）に従って標本波面の推定の角度スペクトルを更新する。
【００５６】
【数８】
2016057172.tif 000010
【００５７】
ただし、ＸａｍｐはステップＳ３０５で更新された標本波面の推定の振幅であり、ＸｐｈａｓｅはステップＳ３０４で求められた標本波面の推定の位相である。ｉは虚数単位である。他の記号は式（４）と同じである。
【００５８】
ステップＳ３０７では、波面演算装置１０５は、以下の式（９）を用いて、標本波面の推定の角度スペクトルから標本波面の推定を求める。
【００５９】
【数９】
2016057172.tif 000011
【００６０】
ただし、ＩＦＴは逆フーリエ変換を表す演算子であり、他の記号は式（１），（２）と同じである。
【００６１】
ステップＳ３０８では、波面演算装置１０５は、ステップＳ３０３〜ステップＳ３０７の処理の反復を終了するか否か、つまりは反復回数の上限値に達したか否かを判断する。反復回数の上限値に達するまで、ステップＳ３０３での誤差項の最小化とステップＳ３０５での正則化項の最小化とが交互に行われる。
【００６２】
本実施例では、反復回数の上限値に達したことを終了条件とするが、他のどのような終了条件を用いてもよい。例えば、前回の反復と今回の反復における標本波面の推定の変化量が所定値以下であることを用いてもよい。終了条件を満たさなかった場合は、ステップＳ３０３に戻る。
【００６３】
以上の処理を行うことで観測波面から標本波面を再構成した結果の例を図４に示す。なお、正解の標本波面の振幅と位相の画像サイズはともに５０×５０画素であり、透過行列のサイズは２５００×２５００である。また、ステップ幅α，βはそれぞれ、１Ｅ−５および１Ｅ−６である。反復回数は５０００回である。また、観測波面は理想的なものではなく、検出器１０４のノイズを模して、最大画素値に対して０．５％の標準偏差を持つ白色ガウスノイズを加えた波面である。また、１画素のサイズは０．５μｍであるが、この数値は便宜的に定めたものであり、特別な意味はない。
【００６４】
図４（ａ）は正解の標本波面の振幅を示す。画素値は０〜１に規格化している。
【００６５】
図４（ｂ）は正解の標本波面の位相を示す。−π〜πの範囲の一様乱数で定めたものである。
【００６６】
図４（ｃ），（ｄ）はそれぞれ、標本波面の振幅と位相の初期解である。初期解の設定方法はステップＳ３０２に示した通りである。
【００６７】
図４（ｅ），（ｆ）はそれぞれ、標本波面の振幅と位相の推定解である。
【００６８】
再構成精度の評価指標は、二乗平均平方根誤差（Root Mean Square Error：RMSE）を用いた以下の式（１０）で与えられる。
【００６９】
【数１０】
2016057172.tif 000012
【００７０】
ただし、ｄは正解の標本波面の振幅または位相であり、ｅは標本波面の推定解の振幅または位相である。ｄｉ，ｅｉはそれぞれ、ｄとｅのｉ番目の要素である。Ｍ×Ｎはｄとｅのサイズである。簡単に言うと、ＲＭＳＥが０に近いほど正解と推定解との誤差が小さい、すなわち再構成精度が良いと言える。
【００７１】
本実施例において、標本波面の正解の振幅とその推定解のＲＭＳＥは０．１６１３であり、標本波面の正解の位相とその推定解のＲＭＳＥは０．６３７３であった。
【００７２】
次に、従来のＴＬＩによって観測波面から標本波面を再構成した結果の例を図５に示す。ここでの再構成は、簡単に言うと、波面伝達媒質１０３を介して観測される観測波面に対して透過行列の（擬似）逆行列を掛けることで行われる。また、ここでは、透過行列が正方行列であるので、逆行列を計算した。
【００７３】
実施例１と同様に、正解の振幅と位相の画像サイズは５０×５０画素であり、透過行列のサイズは２５００×２５００である。また、標本波面の正解と初期解の定め方も実施例１と同様であるので、ここでは示さず、標本波面の推定解のみ示す。
【００７４】
図５（ａ），（ｂ）はそれぞれ、標本波面の推定解の振幅と位相である。標本波面の正解の振幅とその推定解のＲＭＳＥは０．６５４０であり、標本波面の正解の位相とその推定解のＲＭＳＥは１．３３４８であった。
【００７５】
このことから、実施例１により、従来よりも小さい標本波面の推定解のＲＭＳＥが得られることが分かる。すなわち、実施例１の処理は、従来に比べて、標本波面の高い再構成精度を得ることができる。
【実施例２】
【００７６】
実施例２では、強閾値処理（または強閾値反復処理）と呼ばれる画像回復手法を用いて観測波面から標本波面を再構成する処理について説明する。なお、強閾値処理は一般にＩＨＴ（Iterative Hard Thresholding）として知られている。
【００７７】
強閾値処理は、標本波面の振幅をスパース変換してスパース係数を得るステップと、スパース係数に対する閾値処理を行うステップと、閾値処理が行われたスパース係数を逆スパース変換するステップとを含む。
【００７８】
本実施例の処理は、実施例１のステップＳ３０５に相当するステップＳ３０５′において実施例１で用いたＴＶ正則化に代えてＩＨＴを用いる点で実施例１と異なるが、他のステップは同じである。ここでは、ステップＳ３０５′についてのみ説明する。
【００７９】
ステップＳ３０５′では、波面演算装置１０５は、以下の式（１１）に従って、目的関数の正則化項を最小化し、標本波面の推定の振幅を更新する。実施例１のステップＳ３０２と同様に、最小化は逐次近似により行う。
【００８０】
【数１１】
2016057172.tif 000013
【００８１】
ただし、ψは標本波面の推定の振幅を非ゼロ要素が少ないスパースな係数に変換するスパース変換演算子であり、Ｔθは閾値演算子である。他の記号は式（４），（５）と同じである。スパース変換演算子ψとしては、離散コサイン変換や離散ウェーブレット変換等がある。本実施例では、離散ウェーブレット変換を用いるが、スパースな係数に変換できれば、他のどのようなスパース変換演算子を用いてもよい。また、閾値演算子Ｔθは、以下の式（１２）で与えられる。
【００８２】
【数１２】
2016057172.tif 000014
【００８３】
ただし、ｆは任意のスパース係数であり、ｆｉはｆのｉ番目の要素である。θは任意の閾値である。
【００８４】
本実施例において観測波面から標本波面を再構成した結果の例を図６に示す。実施例１と同様に、標本波面の正解の振幅と位相の画像サイズはともに５０×５０画素であり、透過行列のサイズは２５００×２５００である。反復回数は５０００回である。また、ステップＳ３０３での逐次近似におけるステップ幅は１Ｅ−５である。また、ステップＳ３０５′での閾値は１Ｅ−４．５である。標本波面の正解と初期解の定め方は実施例１と同様であるので、ここでは示さず、標本波面の推定解のみを示す。
【００８５】
図６（ａ），（ｂ）はそれぞれ、標本波面の推定解の振幅と位相である。標本波面の正解の振幅とその推定解のＲＭＳＥは０．１４９３であり、標本波面の正解の位相とその推定解のＲＭＳＥは０．６１０１であった。
【００８６】
このことから、実施例２によっても、従来よりも小さい標本波面の推定解のＲＭＳＥが得られることが分かる。すなわち、実施例２の処理も、従来に比べて、標本波面の高い再構成精度を得ることができる。
【００８７】
以上説明した各実施例は代表的な例にすぎず、本発明の実施に際しては、各実施例に対して種々の変形や変更が可能である。
【符号の説明】
【００８８】
１０１  照明光源
１０２  空間光変調素子
１０３  波面伝達媒質
１０４  検出器
１０５  波面演算装置
【図１】
2016057172.tif 000016
【図２】
2016057172.tif 000017
【図３】
2016057172.tif 000018
【図４】
2016057172.tif 000019
【図５】
2016057172.tif 000020
【図６】
2016057172.tif 000021