JP6039236B2 - 画像推定方法、プログラム、記録媒体、画像推定装置、および画像データの取得方法 - Google Patents

Description

本発明は、任意のフォーカス位置での画像を推定する画像推定方法、プログラム、記録媒体、画像推定装置、ネットワーク機器および画像データの取得方法に関する。

バーチャルスライドシステムではバーチャルスライドとよばれるデジタル撮像装置を用いて被写体のデジタル画像を取得する。一般に医療分野では試料を光学素子（カバーグラス）で覆い固定したもの（プレパラート）が被写体として用いられる。バーチャルスライドは顕微光学系、撮像素子、情報処理装置によって構成され、プレパラートをデジタル画像化し、データを保存する。この種の装置ではプレパラートのデジタル画像のみが保存されるため、撮像後に参照できるのは撮像を行ったフォーカス位置での画像だけである。医師はフォーカス位置を変えた一連の画像から試料の立体的構造を判断することが多いため、フォーカス位置を変えた複数枚の画像を撮像しておくことが求められる。

しかし、多数の画像を取得すると、撮像時間が大幅に増大し、データ量も膨大になる。そのため、撮像枚数は最小限に抑えることが望ましい。しかし一方で、撮像枚数を減らすと、診断時に医師の求めるフォーカス位置での画像を提供できない場合が生じる。撮像枚数の削減と任意のフォーカス位置での画像の提供という２つの要望を実現する方法として、必要なフォーカス位置での画像を画像処理で推定する方法が提案されている（特許文献１、非特許文献１）。

特許文献１は、複数のフォーカス位置で取得した画像に対し、光学系の焦点ボケフィルタを利用することで画像を推定する方法を提案している。非特許文献１は、画像をフォーカス位置ｚの関数とし、ｚの多項式展開することで近似的に画像を推定する方法を提案している。

特開２００１−２２３８７４号公報

Ｋｅｎｊｉ Ｙａｍａｚｏｅ ａｎｄ Ａｎｄｒｅｗ Ｒ． Ｎｅｕｒｅｕｔｈｅｒ， "Ｍｏｄｅｌｉｎｇ ｏｆ ｔｈｒｏｕｇｈ−ｆｏｃｕｓ ａｅｒｉａｌ ｉｍａｇｅ ｗｉｔｈ ａｂｅｒｒａｔｉｏｎ ａｎｄ ｉｍａｇｉｎａｒｙ ｍａｓｋ ｅｄｇｅ ｅｆｆｅｃｔｓ ｉｎ ｏｐｔｉｃａｌ ｌｉｔｈｏｇｒａｐｈｙ ｓｉｍｕｌａｔｉｏｎ" Ａｐｐｌｉｅｄ Ｏｐｔｉｃｓ， Ｖｏｌ．５０， Ｎｏ．２０， ｐｐ．３５７０−３５７８，１０ Ｊｕｌｙ ２０１１， ＵＳＡ

特許文献１の方法は、光学系の焦点ボケフィルタを予め知る必要があり、煩雑な予備計測等が必要となる。加えて、顕微鏡のような部分コヒーレント結像系には適用できない。また、非特許文献１の方法は、展開に用いる関数がｚの多項式であるため、得られる画像は近似解である。近似の精度を高めるには、高次まで展開する必要があるため、計算時間が膨大になる。

本発明では、任意のフォーカス位置での画像を簡単かつ精度良く推定する画像推定方法、プログラム、記録媒体、画像推定装置、ネットワーク機器および画像データの取得方法を提供することを例示的な目的とする。

本発明の画像推定手法は、Ｎを２以上の整数として、撮像光学系の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、前記撮像光学系を有する撮像装置を用いて取得した被写体の画像データから、ｚ_ｊの最大値をｚ_ｍａｘ、ｚ_ｊの最小値をｚ_ｍｉｎとして、前記光軸方向の位置ｚ（ｚ_ｍｉｎ≦ｚ≦ｚ_ｍａｘ）での画像データを推定する画像推定方法であって、Ｎ枚の画像データに対して、前記光軸方向に周波数変換を行い、Ｎ枚の変換画像データを算出する画像変換ステップと、ｚ_ｍａｘ、ｚ_ｍｉｎ、ｚ及びΔｚから前記Ｎ枚の変換画像データに対して定まる複素数を、前記Ｎ枚の変換画像データに乗算し、和をとる結合ステップと、を有することを特徴とする。

本発明によれば、任意のフォーカス位置での画像を簡単かつ精度良く推定する画像推定方法、プログラム、記録媒体、画像推定装置、ネットワーク機器および画像データの取得方法を提供することができる。

本発明のバーチャルスライドのブロック図である。（実施例１、２、３、４、５、６、７、８） 本発明による被写体の画像取得から表示までのフローチャートである。（実施例１、２、３、４、５、６、７、８） 本発明による画像取得のフローチャートである。（実施例１、２、３、４、５、６、７、８） 本発明の原理を説明するための光学系である。（実施例１、２、３、４、５、６、７、８） ＣＴＦのｆ_ｙ＝０での断面である。（実施例１、２、３、４、５、６、７、８） 被写体の振幅透過率分布である。（実施例１、２、３、４、５、６、７、８） フォーカス位置を変えて取得した被写体の画像である。（実施例１、２、３、４、５、６、７、８） 本発明の画像推定方法によって得られた推定画像である。（実施例１） 本発明の画像推定方法によって画像推定を行ったフォーカス位置とＰＳＮＲである。（実施例１） 本発明の画像推定方法によって得られた推定画像である。（実施例２） 本発明の画像推定方法によって画像推定を行ったフォーカス位置とＰＳＮＲである。（実施例２） 本発明の画像推定方法のフローチャートである。（実施例３、４） 本発明の画像推定方法によって得られた推定画像である。（実施例３） 本発明の画像推定方法によって画像推定を行ったフォーカス位置とＰＳＮＲである。（実施例３） 本発明の画像推定方法によって得られた推定画像である。（実施例４） 本発明の画像推定方法によって画像推定を行ったフォーカス位置とＰＳＮＲである。（実施例４） 本発明の画像推定方法によって得られた推定画像である。（実施例５） 本発明の画像推定方法によって画像推定を行ったフォーカス位置とＰＳＮＲである。（実施例５） 本発明の画像推定方法によって得られた推定画像である。（実施例６） 本発明の画像推定方法によって画像推定を行ったフォーカス位置とＰＳＮＲである。（実施例６） ＰＳＮＲの最悪値を画像取得ステップΔｚと撮像光学系の瞳の中心遮蔽の大きさεを変えて求めたものである。（実施例７） ＰＳＮＲの最悪値が３５の場合の画像取得ステップΔｚを撮像光学系の瞳の中心遮蔽の大きさεを変えて求めたものである。（実施例７） ＰＳＮＲの最悪値を画像取得ステップΔｚと照明光学系の瞳の中心遮蔽の大きさεを変えて求めたものである。（実施例８） ＰＳＮＲの最悪値が３５の場合の画像取得ステップΔｚを照明光学系の瞳の中心遮蔽の大きさεを変えて求めたものである。（実施例８） 本発明のネットワーク機器から表示位置を指定する場合に行う動作のフローチャートである。（実施例９）

本発明は、撮像光学系を用いて撮像装置によって撮像された試料の画像を、任意のフォーカス位置における画像を推定する画像推定方法に関する。画像推定方法は、コンピュータによって実行可能なプログラムとして具現化され、記録媒体などに記録されてもよい。

画像の推定は、撮像装置において行ってもよいし、撮像装置の画像を格納する記憶装置（メモリ）に接続された画像推定装置（コンピュータ）が行ってもよい。撮像装置は画像推定手段として機能してもよい。また、取得した画像の表示は、撮像装置や画像推定装置に接続された表示手段（ディスプレイ）やこれらと（ＬＡＮ、ＷＡＮ、インターネットなどの）ネットワークを介して接続されたネットワーク機器において行ってもよい。この場合、画像推定方法はクラウドコンピューティングを利用することもでき、ネットワーク機器は位置の指定（入力）と画像データの受信を行うデスクトップ型のパーソナルコンピュータ（ＰＣ）であってディスプレイなどの表示手段に接続されてもよい。あるいは、ネットワーク機器はタッチスクリーンなどの表示手段を備えた携帯端末（ｉＰａｄ（登録商標））やラップトップＰＣ、専用端末（ＰＤＡなど）であってもよい。これによって遠隔診断が可能になる。

図１は本発明の実施に係るバーチャルスライドのブロック図である。バーチャルスライドは、撮像ユニット（撮像装置）１００、制御ユニット２００、情報処理ユニット（画像推定装置）４００から構成される。

制御ユニット２００は搬送部２０１、制御部２０２から構成される。搬送部２０１は制御部２０２の命令に基づき被写体１０３を可動ステージ１０２上へ移動する。可動ステージ１０２は制御部２０２からの命令に従い光軸方向に動くことが可能である。可動ステージ１０２は光軸と垂直方向に移動する機能を備えていてもよい。可動ステージ１０２を用いることでフォーカス位置を変えて画像を取得することができる。

撮像ユニット１００は被写体１０３の画像を取得する装置であり、照明系１０１、可動ステージ１０２、光学系（撮像光学系）１０４、撮像素子１０５からなる。

照明系１０１により可動ステージ１０２に載置された被写体１０３が照明され、光学系１０４によって撮像素子１０５上に被写体の拡大光学像が形成される。撮像素子１０５は被写体の拡大光学像を光電変換する光電変換手段（光電変換素子）である。撮像素子１０５から出力された電気信号は画像データとして情報処理ユニット４００に伝送される。

情報処理ユニット４００はコンピュータ（演算手段）４０１、画像処理部４０２、データ保存部（メモリ）４０３、表示部（表示手段）４０４によって構成される。

画像処理部４０２は撮像素子１０５から送られた画像データをデジタル信号に変換する。このデジタル信号を輝度信号と呼ぶ。画像処理部４０２は輝度信号に変換した画像データ対して、例えば、ノイズ除去や圧縮等の画像処理を施し、コンピュータ４０１に伝送する。送られた画像データをコンピュータ４０１はデータ保存部４０３に転送する。データ保存部４０３は送られた画像データを保存する。診断時は、コンピュータ４０１が保存部４０３から画像データを読み出す。読み出された画像データに対してコンピュータ４０１は画像処理を施し、ユーザーが指定したフォーカス位置の画像データへ変換する。変換された画像データを表示部４０４へ転送し、画像が表示される。

コンピュータ４０１、画像処理部４０２、データ保存部４０３、表示部４０４及び制御部２０２は１つのコンピュータに含まれていてもよい。また、データはネットワーク４５０を介して接続された外部の不図示のサーバーに保存してもよく、この場合、複数の人がリモートでデータにアクセスすることができる。コンピュータ４０１は、ネットワーク４５０を介して様々なネットワーク機器に接続されている。これらのネットワーク機器は、ラップトップＰＣ４６０、デスクトップＰＣ４６２、タッチスクリーン機能を有する携帯端末４６４、ＰＤＡなどの専用端末４６６を含む。

ネットワーク機器は、演算手段、表示手段、指定手段、通信手段、記憶手段を有する。演算手段はコンピュータ（プロセッサ）から構成され、各部を制御や必要な演算を行う。表示手段は、ネットワーク機器４６０、４６４、４６６のようにネットワーク機器の筺体と一体であってもよいし、ネットワーク機器４６２のようにネットワーク機器に接続されたディスプレイ等であってもよい。指定手段は、ユーザーに光学系１０４の光軸方向の任意の位置ｚを指定させることを可能にするタッチスクリーン、キーボード、スタイラスペン、マウス等の入力手段を含む。通信手段は、ネットワーク４５０に接続されて画像推定装置に位置ｚの情報を送信し、画像推定装置から位置ｚにおける被写体１０３の画像データの情報を受信する。画像データの情報は、ｊｐｅｇなどの静止画像でもよいし、それを表す後述する輝度分布（画素値の分布）を表す関数Ｉ（ｘ，ｙ，ｚ）であってもよい。記憶手段は、位置ｚの指定を可能にするためのアプリケーションプログラムを格納したメモリである。ネットワーク機器４６０、４６４、４６６は通信手段を介して受信した画像データの情報に基づいて位置ｚの画像を表示する表示手段を更に有する。

図２に被写体の画像取得から表示までのフローチャートを示す。ここで、「Ｓ」はステップ（工程）の略であり、これは他のフローチャートにおいても同様である。なお、被写体１０３の拡大像を撮像素子１０５、画像処理部４０２、コンピュータ４０１を用いて画像データとして取得することを画像取得と呼ぶ。

Ｓ１で被写体１０３を搬送部２０１を用いて可動ステージ１０２に設置し、Ｓ２で被写体１０３の画像を複数のフォーカス位置で取得する。取得した一連の画像をＳ３でデータ保存部４０３に保存する。診断時には、保存されたデータをＳ４で読み込む。被写体１０３の画像取得と診断が同時に行われる場合、データは一時保存でもよい。Ｓ５で所望のフォーカス位置をユーザーが設定する。フォーカス位置は予め設定されていた値を読み込んでもよいし、予備計測で得た試料の表面形状から算出してもよい。設定したフォーカス位置に対してＳ６で画像推定を行う。推定された画像をＳ７で表示部４０４に送り画像を表示する。以上の流れで被写体の画像取得と表示が行われる。

次に、Ｓ２における被写体１０３の画像取得に関して詳細に説明する。図３に画像取得のフローチャートを示す。

始めに、Ｓ２０１でステージ位置を光軸方向（ｚ方向）に変える範囲と移動幅を指定する。直進光の伝搬方向へのフォーカスずれ量をｚ_ｊとし、光学系１０４に対して撮像素子１０５と共役な位置を原点とする。ｚ_ｊはフォーカス位置及び光軸方向のステージ位置と同義である。Ｓ２０１で画像を取得するｚ_ｊの範囲ｚ_ｍａｘ〜ｚ_ｍｉｎと画像を取得するステージの移動幅（間隔）Δｚを指定する。

即ち、光学系１０４の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、光学系１０４を用いて被写体１０３を撮像することによって被写体１０３の画像データを取得する。ｚ_ｊの最大値をｚ_ｍａｘ、ｚ_ｊの最小値をｚ_ｍｉｎとする。

これらの値はＳ２０１を行う前に予め測定しておいた試料の表面形状から決めてもよいし、ユーザーが任意に設定してもよい。また予め設定しておいた既定の値を用いてもよい。次に、Ｓ２０２で可動ステージ１０２をｚ_ｍｉｎの位置まで移動させ、Ｓ２０３で被写体１０３の画像を取得する。画像取得が終了した後、Ｓ２０４で可動ステージ１０２をΔｚ移動させ、Ｓ２０５で画像を取得する。Ｓ２０４からＳ２０５を可動ステージ１０２の位置がｚ_ｍａｘ以上となるまで繰り返す。以上によりフォーカス位置を変えた複数の（即ち、Ｎ枚の）被写体１０３の画像データが取得される。

なお、可動ステージ１０２の移動方法はこれに限ることはない。例えば、Ｓ２０２で可動ステージ１０２をｚ_ｍａｘに移動した後、−Δｚ移動させて画像を取得してもよい。

本発明では、任意のフォーカス位置での画像を、取得した複数のフォーカス位置での画像データから推定する。その原理を説明する。

撮像光学系によって撮像素子面上に形成される光学像の強度分布を与える式を説明する。

原理を説明するための光学系を図４に示す。図４（ａ）は被写体５０１の像が撮像光学系５０２によって撮像素子５０３上に結像している場合の光学系を示す。図４（ｂ）は図４（ａ）の場合に対して被写体５０１が光軸方向にｚだけずれた場合の光学系を示す。

始めに照明光が垂直入射の完全コヒーレント光である場合を考える。被写体５０１の像が撮像素子５０３上に結像している場合の結像式は次式で表される。

Ｉ（ｘ，ｙ）は撮像素子５０３上での像強度分布、Ｏ（ｆ_ｘ，ｆ_ｙ）は物体からの回折光分布である。また、Ｐ（ｆ_ｘ，ｆ_ｙ）は光学系の瞳関数、ｘ、ｙは像面での光軸に垂直な直交座標、ｆ_ｘ、ｆ_ｙはそれぞれｘ方向の空間周波数、ｙ方向の空間周波数、はフーリエ変換を表す。の下付き添え字の２は光軸方向に直交する２方向であるｘとｙ方向に２次元フーリエ変換（周波数変換）を行うことを示し、上付き添え字の−１は逆フーリエ変換（逆周波数変換）を行うことを示す。図４（ｂ）に示したように、被写体５０１が光軸方向にｚだけ移動した場合の結像式は数式２で表される。

ここでＰ（ｆ_ｘ，ｆ_ｙ；ｚ）は数式３で表され、被写体のフォーカスずれと等価なデフォーカス収差を加えた瞳関数である。

λは照明光の波長、ＮＡは撮像光学系の被写体側の開口数、ｋ_ｚ０は照明光の波数のｚ方向成分、ｋ_--ｚは回折光の波数のｚ方向成分である。ｆ_ｒ＝（ｆ_ｘ ^２＋ｆ_ｙ ^２）^１／２である。簡単のため倍率は１倍とする。撮像光学系の瞳面の中心に円形の遮蔽部があるとし、その半径が瞳の半径のε倍であるとする。本実施例では瞳面の中心付近が遮光されている場合を想定しているが、これに限ることはない。瞳面の中心付近の光強度が周囲の光強度よりも弱くなっていれば同様の効果が期待できる。また、遮蔽部は撮像光学系の瞳面に配置される必要はない。円形である必要もない。瞳面での光強度分布が中心付近で弱くなるように設計されたフィルタが撮像光学系の光路中に配置されていればよい。また、中心遮蔽がない場合は数式３の波括弧を省略した式を考えればよい。

波数のｚ方向成分は数式４及び数式５で表される。

ｆ_ｘ０、ｆ_ｙ０はそれぞれ入射光のｘ、ｙ方向の波数成分を２πで割った値である。ｆ_ｒ０＝（ｆ_ｘ０ ^２＋ｆ_ｙ０ ^２）^１／２である。今は垂直入射を考えているためｆ_ｘ０＝ｆ_ｙ０＝０である。またｃｉｒｃ関数は数式６で表される。

フォーカス位置を変えて取得した画像データは数式２のＩ（ｘ，ｙ，ｚ）をｚを変えて取得していることに相当する。数式２では被写体が光軸方向に移動した場合を想定しているが、撮像素子を光軸方向に動かした場合も同様に議論することができる。

ここで、Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域を考える。Ｉ（ｘ，ｙ，ｚ）を各座標に対してフーリエ変換する。すると数式２は数式７に変形される。

ここで、ＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）は数式８で表され、瞳関数をｚに対してフーリエ変換したものであり、撮像光学系の特性を表す関数である。は畳み込み積分を示し、添え字の３はｆ_ｘ、ｆ_ｙ、ｆ_ｚの３つの座標に対して畳み込み積分を行うことを表す。上付きの＊は複素共役を示す。の下付き添え字の３はｘ、ｙ、ｚの３つの座標に対する３次元フーリエ変換を行うことを示す。

ＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）のｆ_ｙ＝０での断面を図５に示す。数式８よりＣＴＦはｃｉｒｃ関数とδ関数の積である。δ関数の特性からＣＴＦは、ｆ_ｚ＋１／λ−（１／λ^２−ｆ_ｒ ^２）^１／２＝０を満たす曲面でのみ非０となる。これは空間周波数空間で中心がｆ_ｚ＝−１／λ、半径が１／λの球面を表している。更に、ｃｉｒｃ関数の特性からＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）はＮＡε／λ≦ｆ_ｒ≦ＮＡ／λでのみ非０となる。よって、ＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）が非０となるのは図５の太線部分のみとなる。図５よりＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）が非０となるｆ_ｚの領域は数式９に限られる。

Ｏ（ｆ_ｘ，ｆ_ｙ）は一般に空間周波数空間の全域で非０となるため、Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域はＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）が非０となる領域によって決定される。Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域を明確にするため、Ｏ（ｆ_ｘ，ｆ_ｙ）＝１という場合を考える。すると、Ｉ（ｘ，ｙ，ｚ）のスペクトルはＣＴＦ（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）の自己相関で表される。自己相関の結果、Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となるｆ_ｚの領域は数式１０となる。

この結果は任意のＯ（ｆ_ｘ，ｆ_ｙ）に対して成立する。

部分コヒーレント結像系を考える場合では数式２の代わりに数式１１を用いればよい。

Ｓ（ｆ_ｘ０，ｆ_ｙ０）は有効光源の強度分布である。有効光源とは試料がない場合に結像光学系の瞳面に形成される光強度分布である。数式１１における瞳関数Ｐは数式３の瞳関数をｆ_ｘ０及びｆ_ｙ０だけ平行移動させたものであるため、数式１１における瞳関数Ｐから求められるＣＴＦ（ｆ_ｘ，ｆ_ｙ、ｆ_ｚ）の形状は図５に示した太線と変わらない。よって、数式１１における被積分関数の絶対値２乗の項に対して、完全コヒーレントの場合と同様にスペクトルが非０となるｆ_ｚの領域が定まる。従って、部分コヒーレント結像系においてもＩ（ｘ，ｙ，ｚ）のスペクトルが非０となるｆ_ｚの領域は数式１０で表される。

インコヒーレント結像系はＳ（ｆ_ｘ０，ｆ_ｙ０）＝１という場合に等価なため、インコヒーレント結像系においてもＩ（ｘ，ｙ，ｚ）のスペクトルが非０となるｆ_ｚの領域は数式１０で表される。

以上より、撮像方式にかかわらずＩ（ｘ，ｙ，ｚ）のスペクトルが非０となるｆ_ｚの領域は数式１０に限られる。

任意のフォーカス位置での画像を推定する手法を説明する。以下ではＩ（ｘ，ｙ，ｚ）をフォーカス位置ｚで取得した画像データの輝度分布として取り扱う。また、簡単のためｚ_ｍｉｎ＝０として考える。ｚ_ｍｉｎ≠０の場合は以下に示す式のｚをｚ−ｚ_ｍｉｎに置き換えればよい。

Ｉ（ｘ，ｙ，ｚ）を数式１２に従ってｚに対し、ｉを虚数単位としてフーリエ級数展開（複素形フーリエ級数展開）する。

Ｉ’（ｘ，ｙ；ｆ_ｚｎ）はフォーカス位置を変えて取得した画像データＩ（ｘ，ｙ，ｚ_ｊ）から数式１３に示した離散フーリエ変換によって求まるフーリエ係数（複素フーリエ係数）である。ここで、ｚ方向の空間周波数ｆ_ｚは画像を取得するｚの範囲が有限であるため、空間周波数ｆ_ｚは離散的なｆ_ｚｎ（ｆ_ｚｎ＝ｎ／（ｚ_ｍａｘ−ｚ_ｍｉｎ＋Δｚ））となり、ｚ_ｊの範囲ｚ_ｍａｘ−ｚ_ｍｉｎより数式１４によって定まる。ｎはフーリエ係数Ｉ’（ｘ，ｙ；ｆ_ｚｎ）を指定する整数であり、−Ｎ／２≦ｎ＜Ｎ／２を満たす整数である。

Ｎは撮像枚数であり、２以上の整数である。Ｉ’（ｘ，ｙ；ｆ_ｚｎ）は数式１０を満たす領域でのみ非０となるため数式１２は数式１５に示した−ｎ_０〜ｎ_０の範囲の級数と同一である。

ここで、ｎ_０は数式１６を満たす最大のｎである。

数式１６は両辺を（ｚ_ｍａｘ−ｚ_ｍｉｎ＋Δｚ）で割ると数式１４から次式のようになる。

数式１５はフォーカス位置を変えて取得した画像データＩ（ｘ，ｙ，ｚ_ｊ）から数式１３に従ってフーリエ係数を求め、それを係数としてｅｘｐ（−２πｉｆ_ｚｎｚ）を線型結合すれば、任意のフォーカス位置ｚでの画像が推定できることを意味している。即ち、２次元データであるフーリエ係数Ｉ’（ｘ，ｙ；ｆ_ｚｎ）を線型結合することで画像推定が可能となる。

また、サンプリング定理より取得した画像の間隔Δｚから、ｆ_ｚｎの最大値ｆ_ｚｍａｘが数式１８と定まる。

ｆ_ｚｍａｘが数式１０の右辺よりも小さいと折り返し雑音が発生し画像の推定精度が低下する。画像を精度よく推定するためにはｆ_ｚｍａｘが数式１０の右辺よりも大きくなければならない。よって、画像を誤差なく推定するため要求されるΔｚの条件は数式１８と数式１０より数式１９のように定まる。

画像を取得する幅Δｚを数式１９の範囲で設定しておけば、任意のフォーカス位置での画像が精度良く推定可能となる。数式１９より瞳に中心遮蔽がある場合はない場合に比べてΔｚの大きさをＭ倍に広げてよいことを示している。ここでＭは数式２０で表わされる。

ただし、生体試料の画像や、人物、風景などを撮像した画像のスペクトルは一般にｆ_ｘ、ｆ_ｙ、ｆ_ｚの絶対値が大きくなると値が小さくなる。そのため、数式１９よりもさらに大きなΔｚで画像を取得しても、一定の誤差を許容すれば画像の推定は可能である。この場合、Δｚはユーザーが求める画像推定精度から事前にまたは撮像時に設定すればよい。

本実施形態では撮像光学系の瞳に中心遮蔽がある場合を説明したが、照明光学系の瞳に中心遮蔽を配置した場合にも同様の効果が期待できる。

数式１２はｚ方向のみのフーリエ級数展開だが、同時にｘ、ｙ方向に対してもフーリエ級数展開をしてもよいため、次のような展開も可能である。

係数Ｉ’（ｆ_ｘｓ，ｆ_ｙｔ，ｆ_ｚｎ）はＩ（ｘ_ｐ，ｙ_ｑ，ｚ_ｊ）から数式２２より求められるフーリエ係数である。

数式２１は数式１２より煩雑で計算時間も多くなるが、画像処理においてローパスフィルタ等を施す場合は数式２１の形で計算をする方が便利である。

画像推定に用いる関数としてｓｉｎｃ関数を採用してもよい。すなわち画像を数式２３で計算してもよい。

ここで、ｓｉｎｃ関数とはｓｉｎ（πｘ）／πｘで表される関数である。Ｂは画像を取得した間隔Δｚから定まる定数で数式２４である。

数式２３は無限級数であるが、取得した画像データはＮ枚であるため和は有限枚の画像に対してしかとれない。そこで、画像を取得したｚの範囲外では画像が周期的に繰り返されると仮定して、数式２３を数式２５に変更する。この場合、次式で示すように、Ｎ枚の画像データに、ｓｉｎｃ（（ｚ−ｚ_ｊ）／Δｚ）を乗算し、和をとることが行われる。

ここでＴは画像がｚ方向に繰り返される周期でＴ＝ｚ_ｍａｘ−ｚ_ｍｉｎ＋Δｚである。Ｎ_ｖは０以上の整数で、ユーザーが自由に設定できる値である。

画像推定方法として、数式１２に示したフーリエ級数展開を用いたが、級数展開の方法はこれに限ることはない。他にも、コサイン変換を用いることで推定が可能である。その方法を簡単に述べる。

コサイン変換を用いる方法では、Ｉ（ｘ，ｙ，ｚ）に対する展開係数Ｉ’（ｘ，ｙ；ｆ_ｚｎ）を数式２６に示す離散コサイン変換を用いて求め、空間周波数ｆ_ｚｎを数式１４に従って求める。

先のフーリエ級数展開を用いる方法と異なり、ここではｎを０≦ｎ≦Ｎ−１を満たす整数とする。求められたＩ’（ｘ，ｙ；ｆ_ｚｎ）を数式２５に従って線形結合することで、位置ｚでの画像が求められる。

また、コサイン変換を用いる場合、式２７に示されるｃｏｓの括弧内の係数が式１２で示されるｅｘｐの括弧内の係数と２だけ異なるため、数式１７に相当する条件は次式で表される。

この式を満たす最大のｎをｎ_０とすれば、式１５を導出した過程と同様の議論を行うことで、式２７は式２９と等価となる。

以上の手順で任意のｚでの画像を推定することができる。

また、画像推定に用いる式は数式２６、数式２７に限ることはない。ｃｏｓを用いた類似の変換方法はすべて本手法に含まれる。また、ｓｉｎを用いても同様の議論は可能である。

以下、実施例１における任意のフォーカス位置での画像推定方法について説明する。図３で示したフローチャートに従って取得した被写体の画像をシミュレーションによって示す。撮像光学系の開口数を０．７、結像倍率を１０倍、波長を５５０ｎｍで単色、中心に円形の遮蔽部がある光学系を仮定し、円形遮蔽部の半径は瞳の絞りの３割とする。照明は部分コヒーレント照明とし、有効光源は輪帯状とする。照明光学系の開口数をＮＡ_ｉ、有効光源の中心遮蔽部に対する開口数をＮＡ_ｉ＿ｉｎとした時、ＮＡ_ｉ／ＮＡ＝０．７、ＮＡ_ｉ＿ｉｎ／ＮＡ＝０．３を仮定する。

被写体の振幅透過率を図６に示す。被写体としてｙ方向の幅が０．５μｍのラインチャートを用いた。計算上では被写体を２次元の物体と仮定している。もちろん、被写体に厚みがあっても本手法は適用可能である。フォーカス位置は−４μｍ（Ｚ_ｍｉｎ）から４μｍ（Ｚ_ｍａｘ）まで１．０μｍごとに変化させた。図７にフォーカス位置を、Ｚ＝２μｍ、−１μｍ、０μｍと変えて取得した被写体の画像の一例を示す。図７の軸は物体側の座標を示している。

次に、画像の推定方法について説明する。フォーカス位置ｚでの画像の輝度分布をＩ（ｘ，ｙ，ｚ）とする。Ｉ（ｘ，ｙ，ｚ）に対して数式１３に従い、ｚ方向に対して周波数変換である離散フーリエ変換を施し、Ｎ枚の変換画像データであるＩ’（ｘ，ｙ；ｆ_ｚ）を取得する（画像変換ステップ）。また、Ｚ_ｍａｘとＺ_ｍｉｎから数式１４に従って複数の変換画像データに対する周波数ｆ_ｚｎを算出する（周波数算出ステップ）。

続いて、観察したいフォーカス位置ｚを指定する。例えば、ｚ＝−１．５μｍとｚ＝−０．５μｍを指定する。これらの値を数式１２に代入し計算を行う。即ち、Ｎ枚の変換画像データに対してｆ_ｚｎ及びｚから定まる複素数をＮ枚の変換画像データに乗算し、和をとる（結合ステップ）。

得られた画像を図８（ａ）に示す。比較対象として、推定を行わずに取得した場合の画像を図８（ｂ）に示す。両者を見比べると、視認できる差はない。画像の違いを定量化するためにＰＳＮＲ（Ｐｅａｋ Ｓｉｇｎａｌ ｔｏＮｏｉｓｅ Ｒａｔｉｏ）を計算する。ＰＳＮＲとは数式３０と数式３１を用いて計算される値で、画像の類似性を定量化したものである。０以上の値をとり、値が大きいほど画像の類似度が高いことを示す。一般的にはＰＳＮＲが３５ｄＢを超えると高品質な画像とされる。

ここで、Ｉ_ｅｓｔ（ｘ，ｙ，ｚ）は推定した画像データである。数式３０は８−ｂｉｔ画像に対応した式である。

図９に推定を行ったフォーカス位置（ｚ［μｍ］）（横軸）とＰＳＮＲ［ｄＢ］（縦軸）を示す。フォーカス位置を−４μｍから４μｍまで０．１２５μｍ刻みで変化させて推定を行った。ただし、実際に画像を取得した箇所は除いている。推定した全領域でＰＳＮＲが３５ｄＢを超えており推定画像が実際に取得した画像を良く再現していることを示している。

生体試料の画像や、人物、風景などを撮像した画像のスペクトルは一般にｆ_ｘ、ｆ_ｙ、ｆ_ｚの絶対値が大きくなると値が小さくなる。このことからΔｚをλ／２（１−（１−ＮＡ^２）^１／２）より大きくしても、高い精度で画像が推定できると考えられる。この結果に関しては実施例７で説明する。

本実施例では、被写体からの光の波長を５５０ｎｍの単色光と仮定したが、これに限ることはない。例えば、発光ダイオードを複数利用することで多色光による照明及び撮像を行ってもよいし、ハロゲンランプ等を利用することで波長に広がりを持った照明光で照明及び撮像を行ってもよい。この時利用する波長は、生体イメージングや一般的な写真撮影で利用される波長領域、即ち、３００ｎｍ〜１５００ｎｍに含まれていることが望ましい。この場合、撮像素子も前記波長領域に対して応答することが望ましい。また、撮像装置は波長領域の分光データを取得してもよい。

本実施例で撮像装置の結像倍率を１０倍としたが、これに限ることはない、例えば、縮小光学系でもよい。撮像光学系が縮小光学系で照明光学系がインコヒーレント系の場合は一般的なカメラによる撮像を行っている場合に相当する。すなわち、一般的なカメラを用いた写真撮影においても、本発明は実施可能である。

本実施例では収差がない場合を仮定しているが、これに限ることはない。収差が存在しても本実施例は適用可能であることを発明者がシミュレーションによって確認している。本実施例では取得した画像に対してノイズ除去や像ずれの補正等の画像処理を行っていないが、これらの画像処理を施してもよい。本実施例では画像を等間隔に取得した場合を示しているが、これに限ることはない。

また、周波数算出ステップは、ｆ_ｚｎを数式１４に従って算出するものに限定されない。

本実施例では画像を取得する位置ｚ_ｊの間隔Δｚを等間隔としているが、これに限ることはない。実際の計測ではステージの駆動誤差等により等間隔とならない場合が生じるが、この場合でも本手法が適用可能であることを本発明者は確認している。

以下、実施例２における任意のフォーカス位置での画像推定方法について説明する。図３で示したフローチャートに従って取得した被写体の画像をシミュレーションによって求める。計算条件は実施例１と同じであり、得られる画像も図７と同一である。取得した画像の輝度分布Ｉ（ｘ，ｙ，ｚ）を数式２２に従いｘ、ｙ、ｚに対して３次元離散フーリエ変換を施しＩ’（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）を取得する。Ｉ’（ｆ_ｘ，ｆ_ｙ，ｆ_ｚ）を用いて推定画像を数式２１によって計算する。推定するフォーカス位置は実施例１と同じである。

推定した画像の一例を図１０（ａ）に示す。比較対象として、推定を行わずに取得した場合の画像を図１０（ｂ）に示す。両者を見比べると、視認できる差はない。画像の違いを定量化するためにＰＳＮＲを計算する。

図１１に推定を行ったフォーカス位置（ｚ［μｍ］）（横軸）とＰＳＮＲ［ｄＢ］（縦軸）を示す。推定したフォーカス位置は実施例１と同じであり、画像を取得した箇所は除いている。推定した全領域でＰＳＮＲが３５ｄＢを超えており、推定画像が実際に取得した画像を良く再現していることが分かる。また、Δｚをλ／２（１−（１−ＮＡ^２）^１／２）の１．５倍の値まで大きくしても、十分な画質が得られることを確認している。

実施例２ではｘとｙに関するフーリエ変換を行うため計算時間が多くなるが、ノイズ除去などの画像処理を同時に行う場合は好適に使用することができる。ノイズ除去としては、画像に対してローパスフィルタを施す方法が知られている。これは画像のスペクトルに対して高周波成分を減衰させるフィルタを乗算することで実施できる。実施例２ではｘとｙにフーリエ変換を行うことで、画像のスペクトルを求めているため、画像推定の過程でノイズ除去を施すことが可能となる。

以下、実施例３における任意のフォーカス位置での画像推定方法について説明する。図３で示したフローチャートに従って取得した被写体の画像をシミュレーションによって求める。計算条件は実施例１と同じであり、得られる画像も図７と同一である。

実施例３では図１２に示したフローチャートに従って画像を推定する。取得した画像データに対してＳ６０１で数式１３に従いｚに対して離散フーリエ変換（周波数変換）を施し、変換画像データを取得する（画像変換ステップ）。次に、Ｓ６０２で取得した変換画像データに対して、｜ｆ_ｚ｜＞１／２Δの領域（周波数ｆ_ｚｎの最大値よりも高周波領域と最小値よりも低周波領域）に０行列（０値）を付加し、変換画像データを拡張する。この拡張操作をゼロパディングと呼ぶ（ゼロパディングステップ）。ゼロパディングされたデータに対してＳ６０３でｆ_ｚに対して逆フーリエ変換（逆周波数変換）を行う（逆周波数変換ステップ）。

なお、画像変換ステップの周波数変換はフーリエ変換または逆フーリエ変換であれば足りる。周波数算出ステップは、ｆ_ｚｎを数式１４に従って算出する。以上により、取得したΔｚより小さい間隔で画像が等間隔に再現される。Ｓ６０１及びＳ６０３では高速フーリエ変換のアルゴリズムを用いることで計算が速くできる。

推定した画像の一例を図１３（ａ）に示す。比較対象として、推定を行わずに取得した場合の画像を図１３（ｂ）に示す。両者を見比べると、視認できる差はない。画像の違いを定量化するためにＰＳＮＲを計算する。

図１４に推定を行うフォーカス位置（ｚ［μｍ］）（横軸）とＰＳＮＲ［ｄＢ］（縦軸）を示す。推定した全領域で推定したフォーカス位置は実施例１と同じであり、画像を取得した箇所は除いている。推定した全領域でＰＳＮＲが３５ｄＢを超えており、推定画像が実際に取得した画像を良く再現していることが分かる。

また、Δｚをλ／２（１−（１−ＮＡ^２）^１／２）の１．５倍の値まで大きくしても、十分な画質が得られることを確認している。

以上のようにして、撮像を行っていないフォーカス位置での画像を推定することができる。実施例３の場合は複数の推定画像を高速逆フーリエ変換により一度に生成するため、ｘｚ断面の観察を行いたい場合などに好適に使用することができる。

以下、実施例４における任意のフォーカス位置での画像推定方法について説明する。図３で示したフローチャートに従って取得した被写体の像をシミュレーションによって求める。計算条件は実施例１と同じであり、得られる画像も図７と同一である。

実施例４においても図１２に示したフローチャートに従って画像を推定する。ただし、実施例４ではＳ６０１の離散フーリエ変換をｘ、ｙ、ｚの３つの座標に対して実行する。またステップＳ６０３の逆フーリエ変換もｆ_ｘ、ｆ_ｙ、ｆ_ｚの３つの座標に対して実行する。

推定した画像の一例を図１５（ａ）に示す。比較対象として、推定を行わずに取得した場合の画像を図１５（ｂ）に示す。両者を見比べると、視認できる差はない。画像の違いを定量化するためにＰＳＮＲを計算する。

図１６に推定を行うフォーカス位置（ｚ［μｍ］）（横軸）とＰＳＮＲ［ｄＢ］（縦軸）を示す。推定したフォーカス位置は実施例１と同じであり、画像を取得した箇所は除いている。推定した全領域でＰＳＮＲが３５ｄＢを超えており、推定画像が実際に取得した画像を良く再現していることが分かる。

また、Δｚをλ／２（１−（１−ＮＡ^２）^１／２）の１．５倍の値まで大きくしても、十分な画質が得られることを確認している。

実施例４の場合は複数の推定画像を逆フーリエ変換により一度に生成するため、ｘｚ断面の観察を行いたい場合などに好適に使用することができる。実施例４ではｘとｙに関するフーリエ変換を行うため計算時間が多くなるが、ノイズ除去などの画像処理を同時に行う場合は好適に使用することができる。ノイズ除去としては、画像に対してローパスフィルタを施す方法が知られている。これは画像のスペクトルに対して高周波成分を減衰させるフィルタを乗算することで実施できる。実施例４ではｘとｙにフーリエ変換を行うことで、画像のスペクトルを求めているため、画像推定の過程でノイズ除去を施すことが可能となる。

以下、実施例５における任意のフォーカス位置での画像推定方法について説明する。図３で示したフローチャートに従って取得した被写体の像をシミュレーションによって求める。計算条件は実施例１と同じであり、得られる画像も図７と同一である。

取得した画像の輝度分布Ｉ（ｘ，ｙ，ｚ）を用いて、推定画像を数式２５によって計算する。推定するフォーカス位置は実施例１と同じであり、Ｎ_ｖ＝２とした。推定した画像の一例を図１７（ａ）に示す。比較対象として、推定を行わずに取得した場合の画像を図１７（ｂ）に示す。両者を見比べると、視認できる差はない。画像の違いを定量化するためにＰＳＮＲを計算する。

図１８に推定を行うフォーカス位置（ｚ［μｍ］）（横軸）とＰＳＮＲ［ｄＢ］（縦軸）を示す。推定したフォーカス位置は実施例１と同じであり、画像を取得した箇所は除いている。推定した全領域でＰＳＮＲが３５ｄＢを超えており、推定画像が実際に取得した画像を良く再現していることが分かる。

また、Δｚをλ／２（１−（１−ＮＡ^２）^１／２）の１．５倍の値まで大きくしても、十分な画質が得られることを確認している。実施例５ではフーリエ変換を伴わないため、より高速に画像を推定することが可能である。

以下、実施例６における任意のフォーカス位置での画像推定方法について説明する。図３で示したフローチャートに従って取得した被写体の画像をシミュレーションによって求める。計算条件は実施例１と同じであり、得られる画像も図７と同一である。取得した画像の輝度分布Ｉ（ｘ，ｙ，ｚ）を数式２６に従ってｚに対して離散コサイン変換を施しＩ’（ｘ，ｙ；ｆ_ｚ）を取得する。Ｉ’（ｘ，ｙ；ｆ_ｚ）を用いて推定画像を数式２７によって計算する。推定するフォーカス位置は実施例１と同じである。

推定した画像の一例を図１９（ａ）に示す。比較対象として、推定を行わずに取得した場合の画像を図１９（ｂ）に示す。両者を見比べると、視認できる差はない。画像の違いを定量化するためにＰＳＮＲを計算する。

図２０に推定を行ったフォーカス位置（ｚ［μｍ］）（横軸）とＰＳＮＲ［ｄＢ］（縦軸）を示す。推定したフォーカス位置は実施例１と同じであり、画像を取得した箇所は除いている。推定した全領域でＰＳＮＲが３５ｄＢを超えており、推定画像が実際に取得した画像を良く再現していることが分かる。また、Δｚをλ／２（１−（１−ＮＡ^２）^１／２）の１．５倍の値まで大きくしても、十分な画質が得られることを確認している。

実施例６ではｚに関する周波数変換の方法としてコサイン変換を採用している。コサイン変換は実数で計算できるため、フーリエ変換を用いる方法に比べ使用するメモリが少なくて済む。また、コサイン変換の特性からＩ（ｘ，ｙ，ｚ_ｍｉｎ）とＩ（ｘ，ｙ，ｚ_ｍａｘ）の差分が大きい場合、例えば光学系に大きな収差がある場合でも高い精度で推定が可能となる。

以下、実施例７における任意のフォーカス位置での画像推定方法について説明する。実施例７では画像取得間隔と推定精度の関係、及び撮像光学系の瞳に配置する中心遮蔽の効果について説明する。この場合、照明光学系は瞳面の中心部近傍において光強度が周辺部に比べて低い領域を有することになる。

図３で示したフローチャートに従って取得した被写体の画像をシミュレーションによって求める。計算条件として。撮像光学系の開口数を０．７、結像倍率を１０倍、波長を５５０ｎｍで単色とし、撮像光学系の中心に配置した円形の遮蔽部の半径を瞳の大きさのε倍とする。εは０から１までの値をとり０の時は中心遮蔽が存在しない場合に相当する。照明は部分コヒーレント照明とし、有効光源は円形とする。照明光学系の開口数をＮＡ_ｉとし、ＮＡ_ｉ／ＮＡ＝０．７を仮定する。以上の条件を用いて、図６に示した被写体に対して撮像される画像をシミュレーションによって求める。画像を取得するフォーカス位置の範囲は−８μｍ（Ｚ_ｍｉｎ）から８μｍ（Ｚ_ｍａｘ）までとし、画像を取得する間隔をΔｚとする。取得した画像から、位置ｚでの画像を実施例１に示した方法を用いて推定する。推定を行うｚは−８μｍ（Ｚ_ｍｉｎ）から８μｍ（Ｚ_ｍａｘ）の範囲で０．１μｍごとに指定した。各ｚに対して推定された画像を推定を行わずに取得した場合の画像と比較するためＰＳＮＲを計算した。得られたＰＳＮＲから最悪値（縦軸）を抽出し、その値を取得間隔Δｚ（横軸）に対して示したグラフが図２１である。図２１ではΔｚをλ／２（１−（１−ＮＡ^２）^１／２）で規格化している。εを変えた結果をマークの違いで分けて示している。

まずε＝０の場合を考える。Δｚが増加するとＰＳＮＲが減少する。これは取得間隔を広げると推定誤差が大きくなることを示す。画像取得の際には、図１９に示したΔｚとＰＳＮＲの関係から、ユーザーが求める精度目標が達成できるようにΔｚを設定すればよい。例えば、良好な画像をＰＳＮＲが３５ｄＢ以上とすれば、Δｚを１．４以下にすればよい。規格化しない大きさに変換すると、Δｚは１．４×λ／２（１−（１−ＮＡ^２）^１／２）以下にすればよいことになる。

また、図２１よりεが大きい方がΔｚに対してＰＳＮＲ最悪値の減衰が緩やかになる。すなわち撮像光学系の瞳に配置される中心遮蔽が大きい方が、Δｚを大きくしても推定精度を高く保つことができる。

更にΔｚの許容範囲に関して考察する。図２２にＰＳＮＲの最悪値が３５となるΔｚ（縦軸）をε（横軸）に対して示す。図２２の理論予測値はε＝０での結果に対して、数式２０で示されるＭを乗じた結果を示している。図２２よりεが大きくなると、ＰＳＮＲの最悪値が３５となるΔｚが増加する。すなわち、撮像光学系の瞳の中心遮蔽が大きくなると、推定誤差が許容値に収まるΔｚが大きくなるということを示している。さらに、理論予測値はシミュレーション結果をほぼ再現している。すなわち、ε＝０でのΔｚの許容範囲が１．４×λ／２（１−（１−ＮＡ^２）^１／２）以下であったことを考えると、ε≠０の時のΔｚの許容範囲は数式３２となる。

この式はε＝０を包含している。

以下、実施例８における任意のフォーカス位置での画像推定方法について説明する。実施例８では照明光学系の瞳に中心遮蔽を配置した場合、すなわち有効光源を輪帯にした場合の効果について説明する。

図３で示したフローチャートに従って取得した被写体の画像をシミュレーションによって求める。計算条件として。撮像光学系の開口数を０．７、結像倍率を１０倍、波長を５５０ｎｍで単色とし、瞳は円形とする。照明は部分コヒーレント照明とし、有効光源は照明光学系の瞳に配置される中心遮蔽によって輪帯状となっているとする。照明光学系の開口数をＮＡ_ｉ、中心遮蔽部に対する開口数をεとした時、ＮＡ_ｉ／ＮＡ＝０．７を仮定する。

以上の条件を用いて、図６に示した被写体に対して撮像される画像をシミュレーションによって求める。画像を取得するフォーカス位置の範囲、取得間隔Δｚ、推定する位置ｚは実施例７と同じである。更に実施例７と同様に、各ｚに対して推定された画像と推定を行わずに取得した場合の画像からＰＳＮＲを計算し、得られたＰＳＮＲから最悪値（縦軸）を抽出した。その値を取得間隔Δｚ（横軸）に対して示したグラフが図２３である。図２３ではΔｚをλ／２（１−（１−ＮＡ^２）^１／２）で規格化している。εを変えた結果をマークの違いで分けて示している。

図２３よりεが大きい方がΔｚに対してＰＳＮＲ最悪値の減衰が緩やかになる。すなわち照明光学系の瞳に配置される中心遮蔽が大きい方が、Δｚを大きくしても推定精度を高く保つことができる。

Δｚの許容範囲に関して考察する。図２４にＰＳＮＲの最悪値が３５となるΔｚ（縦軸）をε（横軸）に対して示す。図２４の理論予測値はε＝０での結果に対して、数式２０で示されるＭを乗じた結果を示している。図２２よりεが大きくなると、ＰＳＮＲの最悪値が３５となるΔｚが増加する。すなわち、照明光学系の瞳に配置される中心遮蔽が大きくなると、推定誤差が許容値に収まるΔｚが大きくなるということを示している。さらに、理論予測値はシミュレーション結果をほぼ再現している。すなわち、ε＝０でのΔｚの許容範囲が１．４×λ／２（１−（１−ＮＡ^２）^１／２）以下であったことを考えると、ε≠０の時のΔｚの許容範囲は式２７と同じになる。

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

なお、本発明を適用可能な撮像装置は、被写体の拡大画像データを取得する顕微鏡であってもよいし、被写体の縮小画像データを取得するカメラであってもよい。また、撮像装置が画像データを取得する波長の範囲が３００ｎｍから１５００ｎｍの範囲に含まれており、被写体を照明する照明光の波長が３００ｎｍから１５００ｎｍの範囲に含まれていてもよい。

以下、画像推定装置にネットワーク４５０を介して接続されたネットワーク機器４６０、４６２、４６４、４６６（以下、「４６０」で代表させる）から表示位置ｚを指定する場合の動作について説明する。ユーザーは、ネットワーク機器４６０から直ちに表示装置ｚを入力することができないので、まず、表示すべき位置ｚを決定するために画像データを確認することが必要である。このため、ネットワーク機器４６０は、位置ｚの指定（入力）と画像データの受信を行う機能を有すると共に、ディスプレイなどの画像表示手段を有するか、これに接続されている。

図２５に、本実施例のネットワーク機器４６０が行う動作のフローチャートを示す。このフローチャートはアプリケーションプログラムとしてネットワーク機器４６０にインストールされていてもよいし、画像推定装置またはサーバーにインストールされて、そこにアクセスしたネットワーク機器４６０が実行してもよい。

まず、ユーザーはネットワーク機器４６０を操作してネットワーク４５０を介して画像推定装置または不図示のサーバーにアクセスする。不図示のサーバーには不図示の記憶装置（ストレージ）が接続されている。この場合には、画像推定装置は画像データを保存しておらず、推定すべき画像データは記憶装置に保存され、画像推定装置は、推定すべき画像データを記憶装置から取得する。ネットワーク機器４６０が画像推定装置にアクセスする場合には推定すべき画像データを保存する記憶装置は画像推定装置に直接に、または、ＬＡＮ等のネットワークを介して接続されている。以下、ネットワーク機器４６０は、画像推定装置である情報処理ユニット４００にアクセスするものとする。

ネットワーク機器４６０は情報処理ユニット４００にアクセスすると、ユーザー名とパスワードの入力を求められ、これによって認証が行われる。

認証されたユーザー（例えば、医師）は、例えば、患者のＩＤを入力することによってデータ保存部４０３に保存されている患者の情報にアクセスする。データ保存部４０３には、Ｎを２以上の整数として、撮像光学系の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、その撮像光学系を有する撮像装置を用いて取得した被写体（患者）のＮ枚の画像データが保存されている。そこで、ユーザーは、Ｎ枚の画像データのサムネイルやファイル名から少なくとも１枚を選択する。その後、ユーザーは、この画像をネットワーク上で閲覧するか、ダウンロードして保存するか、別の画像を選択するか、をメニューに従って選択する。なお、取得すべき画像は、リストから順次送られてくる画像であってもよい。また、ユーザーは、既に推定された画像データがあればそれを選択してもよい。

ここでは、ユーザーはネットワーク機器４６０において画像データのダウンロードを選択し、画像を取得する（Ｓ５１０）。必要があれば、この時点で、情報処理ユニット４００とのネットワーク接続を一旦切断してもよい。

続いて、ネットワーク機器４６０は、受信した画像データを表示手段に表示させる（Ｓ５０２）。これにより、ユーザーは表示手段に表示された画像を確認することができる。この時点で診断に値する画像が表示されていれば（Ｓ５０３のＹｅｓ）、以降の動作は不要となる。画像データが表示される際には、Ｎ枚の画像データのフォルダのＩＤと、閲覧している画像データのファイルのＩＤも表示される。

但し、一般には、表示される画像は所望の部位にフォーカスが合っているとは限らない（Ｓ５０３のＮｏ）。そのため、ネットワーク機器４６０は、再び、情報処理ユニット４００にアクセスし、認証を受け、保存されている他の画像データがあれば、それを閲覧またはダウンロードする。しかし、ｚ位置の画像データがなければメニューから画像データの推定ボタンをクリックする。

するとネットワーク機器４６０の表示手段には、Ｎ枚の画像データのフォルダのＩＤを入力するか、患者のＩＤに基づいて検索してＮ枚の画像データのフォルダを選択するかのメニューが表示される。ユーザーは、ネットワーク機器４６０の指定手段を介して、Ｎ枚の画像データのフォルダを指定すると、次に、ネットワーク機器４６０の表示手段にはＮ枚の画像データの中で最も近い１枚または２枚の画像データの選択を促すメニューが表示される。例えば、ユーザーが２枚の画像データの間にある画像データを取得したいとして、２枚の画像データのファイルを選択すると、それらの画像データのサムネイルまたはファイル名の下にそれらのｚ_ｊ位置が表示される。その際、ネットワーク機器４６０の表示手段には、入力可能なｚ位置の範囲（即ち、最小値ｚ_ｍｉｎ≦ｚ≦最大値ｚ_ｍａｘ）の情報が併せて表示されてもよい。

次に、ユーザーが所望のｚ位置を入力する。このｚ位置の指定ステップにおいては、ユーザーは、ネットワーク機器４６０の指定手段を介して、その２つのｚ_ｊ位置の間の位置を所望のｚ位置の候補として入力し（Ｓ５０４）、エンターキーを押して画像推定装置へ送信する（Ｓ５０５）。その際、ネットワーク機器４６０は入力されたｚ位置がその最大値ｚ_ｍａｘまたは最小値ｚ_ｍｉｎを超えていればエラーメッセージを表示してもよい。あるいは、入力されたｚ位置が選択された２枚の画像データの中間でない場合にはエラーメッセージを表示してもよい。あるいは、ユーザーが２枚の画像データを選択すると、その中間のｚ位置を２つのｚ_ｊ位置の平均値などによって計算して自動的に決定してもよい。ユーザーが１枚の画像データを選択してｚ位置を入力した場合に、それが隣りの画像データのｚ_ｊ位置に近かったり、隣りの画像データのｚ_ｊ位置を超えていたりした場合にはその旨が表示されてもよい。また、これに加えて隣の画像データを（そのｚ_ｊ位置と共に）表示するかどうか、ｚ位置を変更するかなどを表示してもよい。

ネットワーク機器４６０は画像推定装置が画像の推定を完了するまで待機する（Ｓ５０６）。画像推定が終了したら、ネットワーク機器４６０は推定された画像を、ネットワーク４５０を通じて受信し（Ｓ５０７）、表示手段に表示させる（Ｓ５０８）。ユーザーは表示された画像を観察し、所望の画像が得られていれば処理を終了する（Ｓ５０３のＹｅｓ）。もし、他のｚ位置での画像を望む場合は（Ｓ５０３のＮｏ）、Ｓ５０４からＳ５０８までの動作を繰り返す。この場合は、２枚の画像は今回推定された画像データと、前回推定された画像データのどちらかであってもよい。以上のフローでネットワーク機器４６０はユーザーが求めるｚ位置での画像を提供する。

本発明の画像推定方法は、試料の画像を照明光学、撮像光学系及びデジタルセンサを用いて取得する撮像装置に好適に使用することが可能であり、特にデジタル顕微鏡やデジタルカメラにおいて有用である。

１０３…被写体、１０４…光学系（撮像光学系）、１０５…撮像素子（光電変換手段）、４００…情報処理ユニット（画像推定装置）、４０１…コンピュータ（演算手段）、４５０…ネットワーク、４６０〜４６４…ネットワーク機器

Claims (22)

1. Ｎを２以上の整数として、撮像光学系の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、前記撮像光学系を有する撮像装置を用いて取得した被写体の画像データから、ｚ_ｊの最大値をｚ_ｍａｘ、ｚ_ｊの最小値をｚ_ｍｉｎとして、前記光軸方向の位置ｚ（ｚ_ｍｉｎ≦ｚ≦ｚ_ｍａｘ）での画像データを推定する画像推定方法であって、
   Ｎ枚の画像データに対して、前記光軸方向に周波数変換を行い、Ｎ枚の変換画像データを算出する画像変換ステップと、
   ｚ_ｍａｘ、ｚ_ｍｉｎ、ｚ及びΔｚから前記Ｎ枚の変換画像データに対して定まる複素数を、前記Ｎ枚の変換画像データに乗算し、和をとる結合ステップと、
   を有することを特徴とする画像推定方法。
2. 前記周波数変換はフーリエ変換または逆フーリエ変換であり、
   前記画像推定方法は、ｚ_ｍａｘ、ｚ_ｍｉｎ及びΔｚから、前記Ｎ枚の変換画像データに対応する周波数ｆ_ｚｎを算出する周波数算出ステップを含み、
   前記周波数算出ステップは、前記周波数ｆ_ｚｎを次式に従って算出することを特徴とし、
   前記複素数は虚数単位をｉとしてｅｘｐ（−２πｉｆ_ｚｎｚ）またはｅｘｐ（２πｉｆ_ｚｎｚ）であること、
   を特徴とする請求項１に記載の画像推定方法。
   
3. 前記Ｎ枚の変換画像データに対して、ｚ_ｍａｘ、ｚ_ｍｉｎ及びΔｚから、前記Ｎ枚の変換画像データに対応する周波数ｆ_ｚｎを算出する周波数算出ステップと、
   前記周波数ｆ_ｚｎの最大値よりも高周波領域に０値を付加し、前記周波数ｆ_ｚｎの最小値よりも低周波領域に０値を与えることで前記Ｎ枚の変換画像データを拡張するゼロパディングステップと、
   を更に有し、
   前記結合ステップが前記ゼロパディングステップで拡張された画像データに対して光軸方向に前記周波数変換に対応した逆周波数変換を行う逆周波数変換ステップである、
   ことを特徴とする請求項１に記載の画像推定方法。
4. 前記周波数変換はフーリエ変換または逆フーリエ変換であり、
   前記周波数算出ステップは、前記周波数ｆ_ｚｎを次式に従って算出すること、
   を特徴とする請求項３に記載の画像推定方法。
   
5. 前記周波数変換が光軸と垂直方向に対する周波数変換を含み、
   前記逆周波数変換が前記光軸と垂直方向に対して前記周波数変換に対応する逆周波数変換を含むこと、
   を特徴とする請求項３または４に記載の画像推定方法。
6. Ｎを２以上の整数として、撮像光学系の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、前記撮像光学系を有する撮像装置を用いて取得した被写体の画像データから、ｚ_ｊの最大値をｚ_ｍａｘ、ｚ_ｊの最小値をｚ_ｍｉｎとして、前記光軸方向の位置ｚ（ｚ_ｍｉｎ≦ｚ≦ｚ_ｍａｘ）での画像データを推定する画像推定方法であって、
   Ｎ枚の画像データに、ｓｉｎｃ（（ｚ−ｚ_ｊ）／Δｚ）を乗算し、和をとること、
   を特徴とする画像推定方法。
7. 前記被写体からの光の波長をλ、前記撮像光学系の被写体側の開口数をＮＡとすると、Δｚが次式で表されること、
   を特徴とする請求項１乃至６のうちいずれか１項に画像推定方法。
   
8. 前記被写体からの光の波長をλ、前記撮像光学系の被写体側の開口数をＮＡとすると、Δｚが次式で表されること、
   を特徴とする請求項１乃至６のうちいずれか１項に記載の画像推定方法。
   
9. 前記周波数変換はコサイン変換であり、
   前記画像推定方法はｚ_ｍａｘ、ｚ_ｍｉｎ及びΔｚから、前記Ｎ枚の変換画像データに対応する周波数ｆ_ｚｎを算出する周波数算出ステップを含み、
   前記周波数算出ステップは、前記周波数ｆ_ｚｎを次式に従って算出し、
   前記複素数は次式によって定まる実数であること、
   を特徴とする請求項１に記載の画像推定方法。
   
10. 前記撮像光学系は瞳面の中心部近傍において光強度が周辺部に比べて低い領域を有すること、
    を特徴とする請求項１乃至９のうちいずれか１項に記載の画像推定方法。
11. 前記領域は円形であり、
    前記撮像光学系の被写体側の開口数をＮＡとし、前記領域の半径を前記撮像光学系の瞳の半径のε倍とし、前記被写体からの光の波長をλとすると、Δｚが次式で表されること、
    を特徴とした請求項１０に記載の画像推定方法。
    
12. 前記領域は円形であることを特徴とし、
    前記撮像光学系の被写体側の開口数をＮＡとし、前記領域の半径を前記撮像光学系の瞳の半径のε倍とし、前記被写体からの光の波長をλとすると、Δｚが次式で表されること、
    を特徴とした請求項１０に記載の画像推定方法。
    
13. 前記撮像装置は前記被写体を照明する照明光学系を有し、
    前記照明光学系は瞳面の中心部近傍において光強度が周辺部に比べて低い領域を有することを特徴とする、
    請求項１乃至８のうちいずれか１項に記載の画像推定方法。
14. 前記領域は円形であり、
    前記撮像光学系の被写体側の開口数をＮＡとし、前記領域の半径を前記撮像光学系の瞳の半径のε倍とし、前記被写体からの光の波長をλとすると、Δｚが次式で表されること、
    を特徴とした請求項１３に記載の画像推定方法。
    
15. 前記領域は円形であり、
    前記撮像光学系の被写体側の開口数をＮＡとし、前記領域の半径を前記撮像光学系の瞳の半径のε倍とし、前記被写体からの光の波長をλとすると、Δｚが次式で表されること、
    を特徴とした請求項１３に記載の画像推定方法。
    
16. Ｎを２以上の整数として、撮像光学系の光軸方向に沿って間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、前記撮像光学系を有する撮像装置を用いて取得した被写体のＮ枚の画像データから、ｚ_ｊの最大値をｚ_ｍａｘ、ｚ_ｊの最小値をｚ_ｍｉｎとして、前記光軸方向の位置ｚ（ｚ_ｍｉｎ≦ｚ≦ｚ_ｍａｘ）の画像データを推定する画像推定方法であって、
    前記位置ｚの前記画像データを表す輝度分布を前記光軸方向とそれに直交する２方向からなる３次元の関数で表したものをＩ（ｘ，ｙ，ｚ）、前記被写体の前記Ｎ枚の画像データを表す輝度分布をＩ（ｘ，ｙ，ｚ_ｊ）、前記光軸方向の離散的な空間周波数をｆ_ｚｎ（ｆ_ｚｎ＝ｎ／（ｚ_ｍａｘ−ｚ_ｍｉｎ＋Δｚ））、ｎを−Ｎ／２≦ｎ＜Ｎ／２を満たす整数、前記撮像光学系の被写体側の開口数をＮＡ、前記被写体からの光の波長をλ、次式を満たす最大のｎをｎ_０、前記撮像光学系の瞳面の中心に円形の遮蔽部があるとした場合の前記遮蔽部の半径が瞳の半径のε倍であるとし、前記遮蔽部がない場合はεが０であるとすると、
    
    Ｉ（ｘ，ｙ，ｚ）を−ｎ_０〜ｎ_０の範囲で複素形フーリエ級数展開することによって取得し、前記複素形フーリエ級数展開の複素フーリエ係数はＩ（ｘ，ｙ，ｚ_ｊ）の離散フーリエ変換によって求まること、
    を特徴とする画像推定方法。
17. Ｎを２以上の整数として、撮像光学系の光軸方向に沿って間隔Δｚで決定されたＮ点の異なる位置ｚ_ｊ（１≦ｊ≦Ｎ）において、前記撮像光学系を有する撮像装置を用いて取得した被写体のＮ枚の画像データから、ｚ_ｊの最大値をｚ_ｍａｘ、ｚ_ｊの最小値をｚ_ｍｉｎとして、前記光軸方向の位置ｚ（ｚ_ｍｉｎ≦ｚ≦ｚ_ｍａｘ）の画像データを推定する画像推定方法であって、
    前記位置ｚの前記画像データを表す輝度分布を前記光軸方向とそれに直交する２方向からなる３次元の関数で表したものをＩ（ｘ，ｙ，ｚ）、前記被写体の前記Ｎ枚の画像データを表す輝度分布をＩ（ｘ，ｙ，ｚ_ｊ）、前記光軸方向の離散的な空間周波数をｆ_ｚｎ（ｆ_ｚｎ＝ｎ／（ｚ_ｍａｘ−ｚ_ｍｉｎ＋Δｚ））、ｎを０≦ｎ＜Ｎ−１を満たす整数、前記撮像光学系の被写体側の開口数をＮＡ、前記被写体からの光の波長をλ、次式を満たす最大のｎをｎ_０、前記撮像光学系の瞳面の中心に円形の遮蔽部があるとした場合の前記遮蔽部の半径が瞳の半径のε倍であるとし、前記遮蔽部がない場合はεが０であるとすると、
    
    Ｉ（ｘ，ｙ，ｚ）を１からｎ_０の範囲で級数展開することによって取得し、前記級数展開の係数はＩ（ｘ，ｙ，ｚ_ｊ）の離散コサイン変換によって求まること、
    を特徴とする画像推定方法。
18. 請求項１乃至１７のうちいずれか１項に記載の画像推定方法をコンピュータによって実行可能なプログラム。
19. 請求項１８に記載のプログラムを記録した記録媒体。
20. 請求項１乃至１７のうちいずれか１項に記載の画像推定方法を実行する演算手段を有する画像推定装置。
21. 請求項２０に記載の画像推定装置によって推定された被写体の画像データを取得する方法であって、
    被写体のＮ枚の画像データおよび前記画像推定装置によって既に推定された画像データのうちの少なくとも１枚の画像データの光軸方向の位置を表示するステップと、
    前記画像推定装置が推定すべき、前記被写体の前記画像データの前記光軸方向の位置を指定するステップと、
    を有することを特徴とする方法。
22. 請求項２０に記載の画像推定装置にネットワークを介して接続されたネットワーク機器が前記方法を実行することを特徴とする請求項２１に記載の方法。