(57)【要約】【課題】少ないデータ量で任意のフォーカス位置での画像を簡単かつ精度良く推定する【解決手段】光学系１０４の光軸方向に所定の間隔ΔＺを置いた複数の位置において、光学系を介して撮像されることにより被写体１０３の画像データを取得し、複数の位置とは異なる位置での画像データを、コンピュータ４０１を利用して推定する画像推定方法は、画像データを取得するＳ２０１と、取得した画像データに基づいて、横カットオフ周波数ｆｒｃを算出するＳ２０８と、算出された横カットオフ周波数に基づいて所定の間隔を算出するＳ２１０と、を有する。【選択図】図１
2015169852.tif 000002
【特許請求の範囲】
【請求項１】
撮像光学系の光軸方向に所定の間隔を置いた複数の位置において、前記撮像光学系を介して被写体を撮像することにより画像データを取得し、前記複数の位置とは異なる位置での画像データを、コンピュータを利用して推定する画像推定方法であって、
画像データを取得する画像取得ステップと、
前記画像取得ステップで取得した前記画像データに基づいて、前記光軸方向に垂直な方向の横カットオフ周波数を算出する周波数解析ステップと、
前記周波数解析ステップで算出された前記横カットオフ周波数に基づいて、前記所定の間隔を算出する間隔算出ステップと、
を有することを特徴とする画像推定方法。
【請求項２】
前記周波数解析ステップで算出された前記横カットオフ周波数に基づいて、前記光軸方向の縦カットオフ周波数を算出する周波数変換ステップを更に有し、
前記間隔算出ステップは、前記周波数変換ステップによって算出された縦カットオフ周波数に基づいて前記所定の間隔を算出することを特徴とする請求項１に記載の画像推定方法。
【請求項３】
前記画像取得ステップは、前記画像データを撮像または読み込みによって取得することを特徴とする請求項１または２に記載の画像推定方法。
【請求項４】
前記周波数解析ステップは、
前記画像取得ステップで取得した前記画像データにフィルタをかけるフィルタリングステップと、
前記フィルタリングステップによって得られた前記画像データに基づいて評価値を算出する評価値算出ステップと、
を含み、
前記フィルタの特性を変更しながら前記フィルタリングステップと前記評価値算出ステップを繰り返すことによって、前記評価値が、推定される画像データの画質の許容値を表す基準値以上の前記横カットオフ周波数の最小値を算出することを特徴とする請求項１乃至３のうちいずれか１項に記載の画像推定方法。
【請求項５】
前記フィルタは、前記画像データの空間周波数の絶対値が所定値よりも大きい領域で前記画像データのスペクトルを減衰させ、前記フィルタの特性を変更することは前記所定値を変更することであり、前記横カットオフ周波数は、前記所定値から算出されることを特徴とする請求項４に記載の画像推定方法。
【請求項６】
前記評価値が前記基準値と等しくなる前記横カットオフ周波数を内挿法によって算出することを特徴とする請求項５に記載の画像推定方法。
【請求項７】
前記周波数変換ステップは、前記撮像光学系の開口数をＮＡ、前記被写体からの光の波長をλ、前記横カットオフ周波数をｆｒｃ、前記縦カットオフ周波数をｆｚｃとすると、
ｆｒｃがＮＡ／λよりも小さい場合は次式を解くことによってｆｚｃを算出し、
2015169852.tif 000012
ｆｒｃがＮＡ／λ以上の場合は次式を解くことによってｆｚｃを算出することを特徴とする請求項２に記載の画像推定方法。
2015169852.tif 000013
【請求項８】
前記周波数変換ステップは、前記撮像光学系の開口数をＮＡ、前記被写体からの光の波長をλ、前記横カットオフ周波数をｆｒｃ、前記縦カットオフ周波数をｆｚｃとすると、
ｆｒｃがＮＡ／λよりも小さい場合は次式を解くことによってｆｚｃを算出し、
2015169852.tif 000014
ｆｒｃがＮＡ／λ以上の場合は次式を解くことによってｆｚｃを算出することを特徴とする請求項２に記載の画像推定方法。
2015169852.tif 000015
【請求項９】
前記横カットオフ周波数に基づいて算出される上限値と、前記上限値に３０％を掛けた下限値の範囲内で前記所定の間隔を選択することを特徴とする請求項１乃至８のうちいずれか１項に記載の画像推定方法。
【請求項１０】
算出された前記所定の間隔の情報をデータベースに保存するステップを更に有することを特徴とする請求項１乃至９のうちいずれか１項に記載の画像推定方法。
【請求項１１】
請求項１乃至１０のうちいずれか１項に記載の画像推定方法をコンピュータによって実行させるためのプログラム。
【請求項１２】
請求項１１に記載のプログラムを記録した記録媒体。
【請求項１３】
請求項１乃至１０のうちいずれか１項に記載の画像推定方法を実行する演算手段を有する画像推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、任意のフォーカス位置での画像を推定する画像推定方法、プログラム、記録媒体および画像推定装置に関する。
【背景技術】
【０００２】
バーチャルスライドシステムではバーチャルスライドとよばれるデジタル撮像装置を用いて被写体のデジタル画像を取得する。一般に医療分野では試料を光学素子（カバーグラス）で覆い固定したもの（プレパラート）が被写体として用いられる。バーチャルスライドは顕微光学系、撮像素子、情報処理装置によって構成され、プレパラートをデジタル画像化し、データを保存する。この種の装置ではプレパラートのデジタル画像のみが保存されるため、撮像後に参照できるのは撮像を行ったフォーカス位置での画像だけである。医師はフォーカス位置を変えた一連の画像から試料の立体的構造を判断することが多いため、フォーカス位置を変えた複数枚の画像（スタック画像）を撮像しておくことが求められる。
【０００３】
しかし、多数の画像を取得すると、撮像時間が大幅に増大し、データ量も膨大になる。そのため、撮像枚数は最小限に抑えることが望ましい。一方で、撮像枚数を減らすと、診断時に医師の求めるフォーカス位置での画像を提供できない場合が生じる。撮像枚数の削減と任意のフォーカス位置での画像の提供という２つの要望を実現する方法として、必要なフォーカス位置での画像を画像処理で推定する方法が提案されている。
【０００４】
特許文献１は、複数のフォーカス位置で取得した画像に対し、光学系の焦点ボケフィルタを利用することで画像を推定する方法を提案している。非特許文献１は、画像をフォーカス位置ｚの関数とし、ｚの多項式展開することで近似的に画像を推定する方法を提案している。
【０００５】
撮像枚数を最小限に抑えるためには、画像を取得するフォーカス位置の間隔をより大きくすることが望ましい。しかし、画像の取得間隔を広くすると画像の推定精度が低下してしまう。より少ない枚数から医師の望む画像を提示するためには、推定精度が保証される範囲で画像の取得間隔を最大することが求められる。
【０００６】
画像の取得間隔の決定法として、非特許文献１及び特許文献２では画像の取得間隔を光学系から定まる値にすることを提案している。非特許文献２は、カットオフ周波数の算出方法について説明している。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００１−２２３８７４号公報
【特許文献２】特開２０１２−４２５２５号公報
【非特許文献】
【０００８】
【非特許文献１】Ｋｅｎｊｉ  Ｙａｍａｚｏｅ  ａｎｄ  Ａｎｄｒｅｗ  Ｒ．  Ｎｅｕｒｅｕｔｈｅｒ，  ”Ｍｏｄｅｌｉｎｇ  ｏｆ  ｔｈｒｏｕｇｈ−ｆｏｃｕｓ  ａｅｒｉａｌ  ｉｍａｇｅ  ｗｉｔｈ  ａｂｅｒｒａｔｉｏｎ  ａｎｄ  ｉｍａｇｉｎａｒｙ  ｍａｓｋ  ｅｄｇｅ  ｅｆｆｅｃｔｓ  ｉｎ  ｏｐｔｉｃａｌ  ｌｉｔｈｏｇｒａｐｈｙ  ｓｉｍｕｌａｔｉｏｎ”  Ａｐｐｌｉｅｄ  Ｏｐｔｉｃｓ，  Ｖｏｌ．５０，  Ｎｏ．２０，  ｐｐ．３５７０−３５７８，１０  Ｊｕｌｙ  ２０１１，  ＵＳＡ
【非特許文献２】Ｃ．  Ｊ．  Ｒ．  Ｓｈｅｐｐａｒｄ，  Ｙ．  Ｋａｗａｔａ  ａｎｄ  Ｓ．  Ｋａｗａｔａ，  ”Ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ  ｔｒａｎｓｆｅｒ  ｆｕｎｃｔｉｏｎｓ  ｆｏｒ  ｈｉｇｈ−ａｐｅｒｔｕｒｅ  ｓｙｓｔｅｍｓ”，  Ｊｏｕｒｎａｌ  ｏｆ  ｔｈｅ  Ｏｐｔｉｃａｌ  Ｓｏｃｉｅｔｙ  ｏｆ  Ａｍｅｒｉｃａ  Ａ，  Ｖｏｌ．１１，  Ｎｏ．２，  ｐｐ．５９３−５９８，１  Ｆｅｂ  １９９４，  ＵＳＡ
【発明の概要】
【発明が解決しようとする課題】
【０００９】
特許文献１の方法は、光学系の焦点ボケフィルタを予め知る必要があり、煩雑な予備計測等が必要となる。加えて、顕微鏡のような部分コヒーレント結像系には適用できない。また、非特許文献１の方法は、展開に用いる関数がｚの多項式であるため、得られる画像は近似解である。近似の精度を高めるには、高次まで展開する必要があるため、計算時間が膨大になる。更に、画像の取得間隔を光学系の情報だけから決めると、光学系によって決まる最も微細な解像力に合わせることになる。その結果、必要以上に狭い間隔で画像を取得することになり、データ量の削減効果が小さくなってしまう。
【００１０】
本発明は、少ないデータ量で任意のフォーカス位置での画像を簡単かつ精度良く推定することが可能な画像推定方法、プログラム、記録媒体および画像推定装置を提供することを例示的な目的とする。
【課題を解決するための手段】
【００１１】
本発明の画像推定方法は、撮像光学系の光軸方向に所定の間隔を置いた複数の位置において、前記撮像光学系を介して被写体を撮像することにより画像データを取得し、前記複数の位置とは異なる位置での画像データを、コンピュータを利用して推定する画像推定方法であって、画像データを取得する画像取得ステップと、前記画像取得ステップで取得した前記画像データに基づいて、前記光軸方向に垂直な方向の横カットオフ周波数を算出する周波数解析ステップと、前記周波数解析ステップで算出された前記横カットオフ周波数に基づいて、前記所定の間隔を算出する間隔算出ステップと、を有することを特徴とする。
【発明の効果】
【００１２】
本発明によれば、少ないデータ量で任意のフォーカス位置での画像を簡単かつ精度良く推定することができる。
【図面の簡単な説明】
【００１３】
【図１】本発明のバーチャルスライドのブロック図である。
【図２】本発明による被写体の画像取得から表示までのフローチャートである。
【図３】本発明による画像の取得間隔Δｚを決定するフローチャートである。
【図４】スタック画像が持つ３次元スペクトルを説明する図である。
【図５】被写体の振幅透過率分布である。
【図６】ｆｒ０とＰＳＮＲの関係を例図する図である。
【図７】画像の推定結果に対する評価量である。
【発明を実施するための形態】
【００１４】
本実施形態は、撮像光学系の光軸方向に所定の間隔ΔＺを置いた複数の置において、撮像光学系を介して撮像されることにより被写体の画像データを取得し、複数の位置とは異なる位置での画像データを推定する画像推定方法に関する。また、本実施形態は、画像取得間隔である所定の間隔ΔＺを適切に設定する方法を提案する。画像推定方法は、コンピュータによって実行可能なプログラムとして具現化され、記録媒体などに記録されてもよい。
【００１５】
画像の推定は、撮像装置において実行されてもよいし、撮像装置の画像を格納する記憶装置（メモリ）に接続された画像推定装置（コンピュータ）において実行されてもよい。撮像装置は画像推定手段として機能してもよい。画像の推定はクラウドコンピューティングを利用して実行することもできる。画像の推定に必要な情報、例えばフォーカス位置などの入力や推定結果の受信、及び既存データの受信はクラウド（ネットワーク）上のコンピュータに接続されたネットワーク機器を用いて行ってもよい。ネットワーク機器はフォーカス位置などを指定（入力）できる入力装置や画像データの受信を行う通信手段を備えたコンピュータである。例えば、デスクトップ型のパーソナルコンピュータ（ＰＣ）や携帯端末、あるいはラップトップＰＣや専用端末（ＰＤＡなど）であってもよい。これらはディスプレイなどの表示装置を同時に備えていてもよいし、外部の表示装置に接続されていてもよい。これにより遠隔診断が可能になる。
【００１６】
図１は、本実施形態に係るバーチャルスライドのブロック図である。バーチャルスライドは、撮像ユニット（撮像装置）１００、制御ユニット２００、情報処理ユニット（画像推定装置）４００から構成される。
【００１７】
制御ユニット２００は、搬送部２０１と制御部２０２から構成される。搬送部２０１は、制御部２０２の命令に基づいて、被写体１０３を可動ステージ１０２上へ移動する。可動ステージ１０２は、制御部２０２からの命令に従って光軸方向に動くことが可能である。可動ステージ１０２は光軸と垂直な方向に移動する機能を備えていてもよい。可動ステージ１０２を用いることでフォーカス位置を変えて画像を取得することができる。
【００１８】
撮像ユニット１００は、被写体１０３の画像を取得する装置であり、照明系１０１、可動ステージ１０２、光学系（撮像光学系）１０４、撮像素子１０５からなる。
【００１９】
照明系１０１により可動ステージ１０２に載置された被写体１０３が照明され、光学系１０４によって撮像素子１０５上に被写体の拡大光学像が形成される。撮像素子１０５は被写体の拡大光学像を光電変換する光電変換手段（光電変換素子）である。撮像素子１０５から出力された電気信号は画像データとして情報処理ユニット４００に伝送される。
【００２０】
情報処理ユニット４００は、コンピュータ（演算手段）４０１、画像処理部４０２、データ保存部（メモリ）４０３、表示部（表示手段）４０４によって構成される。
【００２１】
画像処理部４０２は撮像素子１０５から送られた画像データをデジタル信号に変換する。このデジタル信号を輝度信号と呼ぶ。画像処理部４０２は、輝度信号に変換した画像データ対して、例えば、ノイズ除去や圧縮等の画像処理を施し、コンピュータ４０１に伝送する。コンピュータ４０１は、受信した画像データをデータ保存部４０３に転送する。データ保存部４０３は、受信した画像データを保存する。
【００２２】
診断時は、コンピュータ４０１がデータ保存部４０３から画像データを読み出す。読み出された画像データに対してコンピュータ４０１は画像処理を施し、ユーザーが指定したフォーカス位置の画像データへ変換する。変換された画像データを表示部４０４へ転送し、画像が表示される。
【００２３】
コンピュータ４０１、画像処理部４０２、データ保存部４０３、表示部４０４及び制御部２０２は１つのコンピュータに含まれていてもよい。また、データはネットワーク４５０を介して接続された外部の不図示のサーバーに保存してもよく、この場合、複数の人がリモートでデータにアクセスすることができる。コンピュータ４０１は、ネットワーク４５０を介して様々なネットワーク機器に接続されている。これらのネットワーク機器は、ラップトップＰＣ４６０、デスクトップＰＣ４６２、タッチスクリーン機能を有する携帯端末４６４、ＰＤＡなどの専用端末４６６を含む。
【００２４】
ネットワーク機器は、演算手段、表示手段、指定手段、通信手段、記憶手段を有する。演算手段はコンピュータ（プロセッサ）から構成され、各部の制御や必要な演算を行う。表示手段は、ネットワーク機器４６０、４６４、４６６のようにネットワーク機器の筺体と一体であってもよいし、ネットワーク機器４６２のようにネットワーク機器に接続されたディスプレイ等であってもよい。指定手段は、ユーザーに可動ステージ１０２の光学系１０４の光軸方向に対する任意の位置ｚ（フォーカス位置と同義）を指定させることを可能にするタッチスクリーン、キーボード、スタイラスペン、マウス等の入力手段を含む。通信手段は、ネットワーク４５０に接続されて画像推定装置に位置ｚの情報を送信し、画像推定装置から位置ｚにおける被写体１０３の画像データの情報を受信する。画像データの情報は、ｊｐｅｇなどの静止画像でもよい。記憶手段は、位置ｚの指定を可能にするためのアプリケーションプログラムを格納したメモリである。ネットワーク機器４６０、４６４、４６６は通信手段を介して受信した画像データの情報に基づいて位置ｚの画像を表示する表示手段を更に有する。
【００２５】
図２に被写体１０３の設置から表示までのフローチャートを示す。ここで、「Ｓ」はステップ（工程）の略であり、これは他のフローチャートにおいても同様である。なお、被写体１０３の拡大像を撮像素子１０５、画像処理部４０２、コンピュータ４０１を用いて画像データとして取得することを画像取得と呼ぶ。フォーカス位置を変える所定の間隔をΔｚとする。
【００２６】
Ｓ１では、搬送部２０１が被写体１０３を可動ステージ１０２に設置する。Ｓ２では、被写体１０３に対して最適なΔｚを画像処理装置が決定する。Ｓ２の詳細は後述する。Ｓ３では、被写体１０３のスタック画像を取得する。取得したスタック画像をＳ４でデータ保存部４０３に保存する。診断時には、保存されたデータをＳ５で読み込む。被写体１０３の画像取得と診断が同時に行われる場合、データは一時保存でもよい。Ｓ６では、所望のフォーカス位置をユーザーが設定する。フォーカス位置は予め設定されていた値を読み込んでもよいし、予備計測で得た被写体の表面形状から算出してもよい。設定したフォーカス位置に対してＳ７で画像推定を行う。推定された画像をＳ８で表示部４０４に送り画像を表示する。以上の流れで被写体１０３に対する画像取得と表示が行われる。
【００２７】
画像の推定精度が保証されるΔｚの大きさは被写体に依存する。この理由については後述する。Ｓ２では、被写体の画像のスペクトルを解析することでΔｚを決定する。
【００２８】
Ｓ２におけるΔｚの決定方法を図３に示すフローチャートを用いて示す。
【００２９】
まず、Ｓ２０１（画像取得ステップ）で、被写体１０３の画像データを、撮像ユニット１００を用いて取得する。Ｓ２０２では、取得した画像に離散コサイン変換を行い、スペクトルを算出する。Ｓ２０３では、空間周波数の絶対値がｆｒ０より大きな領域のスペクトルを０にする（減衰させる）。フィルタの特性である所定値ｆｒ０の初期値は光学系の開口数をＮＡ、被写体からの光の波長をλとしたときに、２ＮＡ／λよりも大きな値であれば特に限定されない。Ｓ２０４では、逆離散コサイン変換を行い、劣化画像を算出する。Ｓ２０２からＳ２０４によってローパスフィルタがかかった劣化画像が算出されたことになるため、Ｓ２０２からＳ２０４は、Ｓ２０１で取得した画像データにフィルタをかけるフィルタリングステップとして機能する。
【００３０】
この劣化画像に対してステップＳ２０５（評価値算出ステップ）で評価量を算出する。評価量としては画像の輝度分散値、ＳＳＩＭ（Ｓｔｒｕｃｔｕｒａｌ  Ｓｉｍｉｌａｒｉｔｙ）、ＰＳＮＲ（Ｐｅａｋ  Ｓｉｇｎａｌ  ｔｏＮｏｉｓｅ  Ｒａｔｉｏ）等の評価値や、スペクトルの積分値などの評価値が挙げられる。もちろん画像の評価量としてはこれらに限ることはない。
【００３１】
Ｓ２０６では、評価量を基準値と比較する。基準値は推定画像の画質を決めるための値である。基準値はユーザーが自由に設定してもよいし、予め保存されていた値を読み込んでもよい。好ましくは、画像が良好とされる値、例えばＰＳＮＲならば３５［ｄＢ］程度の値を設定するのが良い。評価量が基準値よりも大きい場合には、ｆｒ０を小さくすることによってフィルタ特性を変更し（Ｓ２０７）、Ｓ２０３からステップＳ２０６までを繰り返す。ｆｒ０を小さくする量としては画像のサイズから決まる周波数空間のサンプリング幅程度の大きさが好ましい。
【００３２】
評価量が、推定される画像の画質の許容値を表す基準値を下回ったら繰り返しを止め、Ｓ２０８（周波数解析ステップ）へと進む。１つ前の繰り返しで用いたｆｒ０が画像の劣化が許容できる空間周波数（カットオフ周波数）と考えることができる。ここまでのステップはコサイン変換を基にした理論を用いていたが、この後のステップはフーリエ変換を基にした理論を用いる。コサイン変換で用いる周波数はフーリエ変換の場合に比べて２倍異なるので、繰り返し計算が終わる１つ前のｆｒ０を２で割った値を光軸と垂直な方向のカットオフ周波数（横カットオフ周波数）ｆｒｃとする（Ｓ２０８）。Ｓ２０８では、フィルタの特性を変更しながらＳ２０３〜Ｓ２０５を繰り返すことによって、評価値が基準値以上の横カットオフ周波数の最小値を算出する。最小値は、離散的に更新されたｆｒ０から算出されてもよいし、評価値が基準値に等しくなるときのｆｒｃを内挿法によって求めてもよい。
【００３３】
Ｓ２０９では、Ｓ２０８までで得られたｆｒｃに基づいて光軸方向のカットオフ周波数（縦カットオフ周波数）ｆｚｃを算出する周波数変換ステップである。ｆｚｃは数式１または数式２を解くことで求めることができる。ＮＡは光学系１０４の開口数であり、λは光学系１０４に使用される光の波長である。ｆｒｃがＮＡ／λよりも小さい場合は数式１を解き、ｆｒｃがＮＡ／λ以上の場合は数式２を解くことによってｆｚｃを算出する。
【００３４】
【数１】
2015169852.tif 000003
【００３５】
【数２】
2015169852.tif 000004
【００３６】
本式は非特許文献２に記載の数式１３を基に導出している。Ｓ２１０（間隔算出ステップ）では、Ｓ２０９で算出されたｆｚｃから画像取得間隔Δｚを算出する。サンプリング定理よりカットオフ周波数がｆｚｃであるデータを完全に取得するために必要な最大のΔｚは次式で定まる。
【００３７】
【数３】
2015169852.tif 000005
【００３８】
以上の手順で定められたΔｚを用いてスタック画像を取得すれば、データ量を最大限抑えることが可能となる。なお、本実施形態では、中間値である縦カットオフ周波数ｆｚｃを求めているが、数式３のｆｚｃに数式１または２を代入することによって中間値の算出を省略することができる。即ち、本実施形態は、算出されたｆｒｃに基づいて、所定の間隔ΔＺを算出すれば足りる。
【００３９】
ここで被写体によって必要なΔｚが変化する理由について説明する。非特許文献２によればスタック画像が持つ３次元のスペクトルは図４に示された白い領域にのみ存在する。ここでｆｒは光軸と垂直な方向の空間周波数、ｆｚは光軸方向の空間周波数である。ｆｒはＮＡ／λ、ｆｚは２（１−（１−ＮＡ２）１／２）／λで規格化されている。
【００４０】
また、データ量の増加が許容されれば、算出されたΔＺよりも多少小さい値が所定の間隔ΔＺとして採用されてもよい。例えば、ｆｒｃに基づいて算出される上限値と、上限値に３０％を掛けた下限値の範囲内でΔＺを選択してもよい。
【００４１】
ｆｒの絶対値がＮＡ／λより小さい時、白い領域はｆｒの絶対値が増加するにつれてｆｚ方向に広がる。言いかえれば、ｆｒの絶対値が小さい時、３次元のスペクトルのｆｚ方向への広がりは小さくなる。すなわち、被写体の持つｆｒ方向のスペクトルが存在する領域が低周波になればなるほど、ｆｚ方向のスペクトルの広がりは小さくなる。これを実空間で考えると、被写体の光軸と垂直な方向の構造が粗ければ粗いほど、光軸方向の画像の広がりは大きくなる、すなわちΔｚを大きくすることができる。以上が被写体に依存してΔｚが異なる理由である。
【００４２】
今、被写体の横カットオフ周波数ｆｒｃが図４のＡの線であるとすると、ｆｚｃは白い領域の境界線との交点（図４黒点）から縦カットオフ周波数ｆｚｃを求めることができる。数式１は白領域の第一象限における境界線を表すため、Ｓ２０８までで求めたｆｒｃを数式１に代入することでｆｚｃを計算できる。ただし、ｆｚはｆｒがＮＡ／λの時に最大となるため、ｆｒｃがＮＡ／λよりも大きい時は数式２を用いる必要がある。ｆｚｃが分かればサンプリング定理から必要なΔｚは数式３よりただちに導出されるため、被写体の画像を解析することで適切なΔｚを定めることができる。
【００４３】
本方法では横カットオフ周波数ｆｒｃをＳ２０１からＳ２０８によって定めることを例として挙げているが、先のステップの本質的な点は被写体の持つ横カットオフ周波数ｆｒｃを取得することである。それゆえ、ｆｒｃを求める方法は例示された方法に限られるものではない。
【００４４】
例えば、Ｓ２０１で被写体の画像を撮像する代わりに被写体と同じ部位で症例が同じ試料の画像を読み込んでもよい。同じ部位で同じ症例であれば被写体が持つスペクトルは類似している可能性が高いため、撮像をせずともｆｒｃを推定することが可能となる。既存の画像を使うことができるため、医者が既に所有している画像データを活用することが可能である。また、観察対象の典型的な大きさから横カットオフ周波数を算出してもよい。この場合は観察対象を円や楕円などモデル化し、モデル画像をＳ２０１で読み込んでもよい。
【００４５】
Ｓ７で行う画像推定方法を説明する。まず原理を説明する。今スタック画像の輝度値分布をＩ（ｘ，ｙ，ｚ）とする。Ｉ（ｘ，ｙ，ｚ）を数式４に従ってｚに対してフーリエ級数展開（複素形フーリエ級数展開）する。
【００４６】
【数４】
2015169852.tif 000006
【００４７】
ｉは虚数単位である。空間周波数ｆｚはｆｚｎ＝ｎ／（ｚｍａｘ−ｚｍｉｎ＋Δｚ）で表わされる離散的な値である。ｚｍａｘは画像を取得する最大のｚ、ｚｍｉｎは画像を取得する最小のｚ、ｚｊは画像を取得した位置でのｚ、ｎはフーリエ係数Ｉ’（ｘ，ｙ；ｆｚｎ）を指定する整数であり、−Ｎ／２≦ｎ＜Ｎ／２を満たす整数である。Ｎは撮像枚数であり、２以上の整数である。Ｉ’（ｘ，ｙ；ｆｚｎ）はフォーカス位置を変えて取得した画像データＩ（ｘ，ｙ，ｚｊ）から求められるフーリエ係数（複素フーリエ係数）である。これらのフーリエ係数はｅｘｐ（−２πｉｆｚｎｚ）の直交性から、数式５に示した離散フーリエ変換によって求めることができる。
【００４８】
【数５】
2015169852.tif 000007
【００４９】
数式４に所望のｚを代入すれば、所望の位置での画像が算出される。よって、画像推定の手順は次のようになる。まず、取得したスタック画像（複数の位置の画像データ）の輝度値分布から数式５に示した離散フーリエ変換（光軸方向への周波数変換）によってフーリエ係数Ｉ’（ｘ，ｙ；ｆｚｎ）（変換画像データ）を算出する。続いて、得られたフーリエ係数Ｉ’（ｘ，ｙ；ｆｚｎ）に対して画像を算出する位置ｚと空間周波数ｆｚｎから定まる複素数である係数ｅｘｐ（−２πｉｆｚｎｚ）を乗算し、−Ｎ／２≦ｎ＜Ｎ／２で和を取る。以上の演算を行うことで所望の位置ｚでの画像を算出することができる。
【００５０】
画像を推定する方法は先に述べた方法に限ることはない。例えば、数式４で級数展開に用いた関数としてｃｏｓ関数やｓｉｎ関数を用いてもよい。例えば、複数の位置の画像データに、ｓｉｎｃ（（ｚ−ｚｊ）／Δｚ）を乗算し、和をとるなどである。また、非特許文献１記載の方法を用いることも可能である。
【００５１】
Ｓ２０８はＳ２０７までに行った周波数変換と数式１を導出すために用いた周波数変換が異なることに起因するものである。そのため、ｆｒｃを算出する方法に応じて演算内容は変化する。例えば、Ｓ２０２で離散フーリエ変換を行えばＳ２０８は不要となる。
【００５２】
本実施形態で算出されたΔｚは被写体の情報とともに記録装置に保存してもよい。こうすれば同一の部位の画像を取得する際に、新たにΔｚの算出を行う必要がなくなる。
【００５３】
本実施例ではｆｒｃをｆｚｃに変換する式として数式１または数式２を用いたがこれに限ることはない。例えば、近軸近似が十分に成立する（ＮＡ＜＜１）光学系では数式１及び数式２を近似した数式６及び数式７を利用することも可能である。
【００５４】
【数６】
2015169852.tif 000008
【００５５】
【数７】
2015169852.tif 000009
【００５６】
本実施例ではＳ２０８でｆｒｃをｆｒ０の１／２としたがこれに限ることはない。１／２が必要となるのは各ステップにおいて基にしている周波数変換が異なることが原因である。他の周波数変換によって本発明を実施する際は利用する周波数変換に応じて各ステップでの周波数を対応付ければよい。Ｓ２１０でΔｚを算出するために数式３を用いているが、このステップで用いる式も使用する周波数変換に応じて変化させればよい。
【００５７】
以下の実施例では、本発明の撮像装置、及び、撮像方法、本発明の効果、更なる応用方法について、図を用いて詳細に説明する。
【実施例１】
【００５８】
以下、実施例１における任意のフォーカス位置での画像推定方法について説明する。本実施例では光学系の条件として撮像光学系の開口数ＮＡを０．２５、結像倍率を１０倍、波長λを６５０ｎｍで単色、照明は部分コヒーレント照明とし、照明光学系の開口数をＮＡｉとして、ＮＡｉ／ＮＡ＝０．７を仮定する。図２及び図３に示したフローチャートに従ってΔｚを決定し、任意の位置での画像を表示する方法を具体的に説明する。本実施例では数値シミュレーションによってその効果を示す。
【００５９】
まず、計算に用いる被写体として図５に示す振幅透過率を設定する（Ｓ１）。ｘ及びｙは光軸と垂直な方向に対する直交した２つの座標である。続いてΔｚを設定する（Ｓ２）。
【００６０】
Δｚを設定するために図５の被写体を光学系１０４を介して撮像し画像データを取得する（Ｓ２０１）。得られた画像に対して空間周波数の絶対値がｆｒ０より大きい領域の値を０にするローパスフィルタを施し劣化画像を算出する（Ｓ２０２乃至Ｓ２０４）。ｆｒ０の初期値として画像のサンプリングピッチから定まる最大の空間周波数を用いる。評価値として劣化画像とフィルタリング前の画像とのＰＳＮＲを計算する（Ｓ２０５）。
【００６１】
ここで、ＰＳＮＲとは数式８と数式９を用いて計算される値で、画像の類似性を定量化したものである。０以上の値をとり、値が大きいほど画像の類似度が高いことを示す。一般的にはＰＳＮＲが３５［ｄＢ］を超えると高品質な画像とされる。
【００６２】
【数８】
2015169852.tif 000010
【００６３】
【数９】
2015169852.tif 000011
【００６４】
ここで、Ｉ（ｘ，ｙ）は真の画像データの輝度値分布、Ｉｅｓｔ（ｘ，ｙ）は推定した画像データの輝度値分布、Ｎｐ、Ｎｑはｘ方向のサンプリング数、ｙ方向のサンプリング数である。数式８は８−ｂｉｔ画像に対応した式である。
【００６５】
算出されたＰＳＮＲは基準値３５［ｄＢ］より大きいため（Ｓ２０６）、ｆｒ０を小さくして（Ｓ２０７）、再計算を行う（Ｓ２０３乃至Ｓ２０５）。
【００６６】
この計算を繰り返して得られたｆｒ０（横軸）に対するＰＳＮＲ（縦軸）を図６に示す。図６の右端の点が２回目の繰り返し計算で得られた結果である。ｆｒ０を小さくしていくとＰＳＮＲが減少し、ｆｒ０が０．４７４［１／μｍ］となった時点で基準値３５［ｄＢ］を下回る。
【００６７】
１つめの繰り返しで用いたｆｒ０を２で割った値、すなわち０．２５７［１／μｍ］を横カットオフ周波数ｆｒｃとする（Ｓ２０８）。得られたｆｒｃはＮＡ／λ＝０．３８５より小さいため数式１を用いて縦カットオフ周波数ｆｚｃを算出する（Ｓ２０９）。計算の結果、ｆｚ０は０．０４３５［１／μｍ］となる。得られたｆｚｃから画像取得間隔Δｚを数式３によって算出すると１１．５［μｍ］である（Ｓ２１０）。
【００６８】
算出されたΔｚを用いてスタック画像を取得し、保存する（Ｓ３、Ｓ４）。得られたスタック画像から、スタック画像を取得していないｚ位置での画像を推定する（Ｓ５乃至Ｓ７）。
【００６９】
画像推定の効果を検証するため、スタック画像を取得したフォーカス範囲−４６〜４６［μｍ］の内で２［μｍ］刻みに計算した。推定結果を評価するため、推定された画像と真の画像とのＰＳＮＲを計算する。計算結果を図７に示す。横軸はｚ座標であり、縦軸はＰＳＮＲである。画像を取得した位置でのＰＳＮＲは推定結果に対する値よりも十分に大きいため表示していない。図７よりすべてのｚにおいてＰＳＮＲが３５［ｄＢ］を上回っている。これは、すべてのフォーカス位置で良好な画像が推定可能であることを示している。
【００７０】
更に、Ｓ２１０までで定まったΔｚ＝１１．５［μｍ］は、光学理論だけから定まる値λ／２（１−（１−ＮＡ）１／２）＝１０．２［μｍ］よりも大きい。これは、画像の特性を利用することでデータ削減効果を高めることが可能であることを示している。
【００７１】
本実施例では、被写体からの光の波長を６５０ｎｍの単色光と仮定したが、これに限ることはない。例えば、発光ダイオードを複数利用することで多色光による照明及び撮像を行ってもよいし、ハロゲンランプ等を利用することで波長に広がりを持った照明光で照明及び撮像を行ってもよい。また、撮像装置は波長領域の分光データを取得してもよい。
【００７２】
本実施例で撮像装置の結像倍率を１０倍としたが、これに限ることはない、例えば、縮小光学系でもよい。撮像光学系が縮小光学系で照明光学系がインコヒーレント系の場合は一般的なカメラによる撮像を行っている場合に相当する。すなわち、一般的なカメラを用いた写真撮影においても、本発明は実施可能である。
【００７３】
本実施例では取得した画像に対してノイズ除去や像ずれの補正等の画像処理を行っていないが、これらの画像処理を施してもよい。
【００７４】
本実施例では画像を等間隔に取得した場合を示しているが、これに限ることはない。例えばスタック画像を取得する最中にフォーカス位置に応じてΔｚの設定を変えることも可能である。
【００７５】
本実施例では１枚の画像からΔｚを算出したが、これに限ることはない。複数枚の画像を取得し、各々画像に対して算出されるΔｚの最小値を取る方法なども考えうる。
【００７６】
本実施例では被写体の画像からΔｚを算出したが、これに限ることはない。被写体と同一の部位の画像を医師が既に所有している画像を読み込み、その画像から算出してもよい。
【００７７】
本実施例で算出したΔｚは、観察した試料の情報、例えば試料の部位の名称や症例など共にデータ保存部に保存されてもよい。同じ部位または同じ症例の標本を観察する場合、記憶しておいたΔｚを読み込み、これを用いてもよい。こうすることで、本特許に示した計算を省略することができる。また、本発明で求めたΔｚと部位の名称、症例などのデータをデータベース化し、通信すれば、他の顕微鏡システムで活用することができる。なお、カテゴリー分けはどのようにしてもよい。
【００７８】
以上、本発明の実施例について説明したが、本発明はこれらの実施例に限定されず、その要旨の範囲内で種々の変形及び変更が可能である。上述したレンズ鏡筒を有する撮像装置も本発明の一側面を構成する。
【産業上の利用可能性】
【００７９】
本発明の画像推定方法は、試料の画像を照明光学、撮像光学系及びデジタルセンサを用いて取得する撮像装置に好適に使用することが可能であり、特にデジタル顕微鏡やデジタルカメラにおいて有用である。
【符号の説明】
【００８０】
１０３…被写体、１０４…光学系（撮像光学系）、４００…情報処理ユニット（画像推定装置）、４０１…コンピュータ（演算手段）
【図７】
2015169852.tif 000016
【図１】
2015169852.tif 000017
【図２】
2015169852.tif 000018
【図３】
2015169852.tif 000019
【図４】
2015169852.tif 000020
【図５】
2015169852.tif 000021
【図６】
2015169852.tif 000022