(57)【要約】【課題】撮像素子によって取得された画像から低周波領域においてランダムに発生するノイズを効果的に低減すること【解決手段】光軸方向の複数の異なる位置において撮像装置を用いて撮像および統合された被写体の３次元の画像データに含まれるノイズを低減する画像処理方法は、前記３次元の画像データに対して、前記光軸方向に対して周波数変換を行い、３次元の変換画像データを算出し、前記３次元の変換画像データの絶対値を特定の周波数領域で小さくし、３次元の変調画像データを算出し、前記３次元の変調画像データに対して、前記光軸方向に対して前記周波数変換に対応する逆周波数変換を行い、３次元の逆変換画像データを算出する。前記特定の周波数領域は、｜ｆｚ｜＞０．０１ＮＡ２／λで表わされる領域の一部である。ここで、ｆｚは前記光軸方向の空間周波数、ＮＡは前記撮像光学系の被写体側の開口数、λは前記被写体からの光の波長である。【選択図】図１
2014096640.tif 000002
【特許請求の範囲】
【請求項１】
撮像光学系の光軸方向の複数の異なる位置において前記撮像光学系と撮像素子を有する撮像装置を用いて撮像および統合された被写体の３次元の画像データに含まれるノイズを、演算手段を利用して低減する画像処理方法であって、
前記３次元の画像データに対して、前記光軸方向に対して周波数変換を行い、３次元の変換画像データを算出する画像変換ステップと、
前記３次元の変換画像データの絶対値を特定の周波数領域で小さくし、３次元の変調画像データを算出する画像変調ステップと、
前記３次元の変調画像データに対して、前記光軸方向に対して前記周波数変換に対応する逆周波数変換を行い、３次元の逆変換画像データを算出する逆画像変換ステップと、
を有し、
前記特定の周波数領域は、次の条件式で表わされる領域の一部であることを特徴とする画像処理方法。
｜ｆｚ｜＞０．０１ＮＡ２／λ
ここで、ｆｚは前記光軸方向の空間周波数、ＮＡは前記撮像光学系の被写体側の開口数、λは前記被写体からの光の波長である。
【請求項２】
前記周波数変換はフーリエ変換であることを特徴とする請求項１に記載の画像処理方法。
【請求項３】
前記特定の周波数領域は、次の条件式で表される領域の一部であることを特徴とする請求項２に記載の画像処理方法。
2014096640.tif 000033
【請求項４】
前記特定の周波数領域は、次の条件式で表される領域の一部であることを特徴とする請求項２または３に記載の画像処理方法。
2014096640.tif 000034
【請求項５】
前記周波数変換はコサイン変換であることを特徴とする請求項１に記載の画像処理方法。
【請求項６】
前記特定の周波数領域が、以下の条件式で表される領域であることを特徴とする請求項５に記載の画像処理方法。
2014096640.tif 000035
【請求項７】
前記特定の周波数領域が、以下の条件式で表される領域であることを特徴とする請求項５または６に記載の画像処理方法。
2014096640.tif 000036
【請求項８】
撮像光学系の光軸方向の複数の異なる位置において前記撮像光学系と撮像素子を有する撮像装置を用いて撮像および統合された被写体の３次元の画像データに含まれるノイズを、演算手段を利用して低減する画像処理方法であって、
前記３次元の画像データに対して、前記光軸方向および光軸に垂直な２つの直交方向に対して周波数変換を行い、３次元の変換画像データを算出する画像変換ステップと、
前記３次元の変換画像データの絶対値を特定の周波数領域で小さくし、３次元の変調画像データを算出する画像変調ステップと、
前記３次元の変調画像データに対して、前記光軸方向および前記２つの直交方向に対して前記周波数変換に対応する逆周波数変換を行い、３次元の逆変換画像データを算出する逆画像変換ステップと、
を有し、
前記特定の周波数領域は、以下の３つの条件式のうち少なくとも１つで表わされる領域の一部であることを特徴とする画像処理方法。
2014096640.tif 000037
｜ｆｒ｜＞０．０１ＮＡ／λ         
｜ｆｚ｜＞０．０１ＮＡ２／λ              
ここで、ｆｚは前記光軸方向の空間周波数、ｆｒは前記撮像光学系の光軸に垂直な方向の空間周波数、ＮＡは前記撮像光学系の被写体側の開口数、λは前記被写体からの光の波長である。
【請求項９】
前記特定の周波数領域は、更に以下の条件式で表わされる領域の一部であることを特徴とする請求項８に記載の画像処理方法。
｜ｆｚ｜＞０．５ＮＡ２／λ
【請求項１０】
前記特定の周波数領域は、更に以下の条件式で表わされる領域の一部であることを特徴とする請求項８に記載の画像処理方法。
｜ｆｒ｜＞０．９ＮＡ／λ
【請求項１１】
前記周波数変換はフーリエ変換であることを特徴とする請求項８乃至１０のうちいずれか１項に記載の画像処理方法。
【請求項１２】
前記特定の周波数領域は、更に次の２つの条件式のうち少なくとも１つで表される領域の一部であることを特徴とする請求項１１に記載の画像処理方法。
2014096640.tif 000038
【請求項１３】
前記周波数変換はコサイン変換であることを特徴とする請求項８乃至１０のうちいずれか１項に記載の画像処理方法。
【請求項１４】
前記特定の周波数領域は、更に以下の３つの条件式のうち少なくとも１つで表わされる領域の一部であることを特徴とする請求項１３に記載の画像処理方法。
2014096640.tif 000039
｜ｆｒ｜＞０．０２ＮＡ／λ
｜ｆｚ｜＞０．０２ＮＡ２／λ
【請求項１５】
前記特定の周波数領域が、更に以下の２つの条件式のうち少なくとも１つで表される領域の一部であることを特徴とする請求項１３に記載の画像処理方法。
2014096640.tif 000040
【請求項１６】
以下の条件式が満たされることを特徴とする請求項１乃至１５のうちいずれか１項に記載の画像処理方法。
2014096640.tif 000041
ただし、Δｚは前記光軸方向の前記複数の異なる位置の間隔である。
【請求項１７】
前記画像変調ステップは前記３次元の変換画像データの前記特定の周波領域での値を０にするステップであることを特徴とする請求項１乃至１６のうちいずれか１項に記載の画像処理方法。
【請求項１８】
請求項１乃至１７のうちいずれか１項に記載の画像処理方法をコンピュータによって実行可能なプログラム。
【請求項１９】
請求項１８に記載のプログラムを記録した記録媒体。
【請求項２０】
請求項１乃至１７のうちいずれか１項に記載の画像処理方法を実行する演算装置を有する撮像装置。
【請求項２１】
顕微鏡であることを特徴とする請求項２０に記載の撮像装置。
【請求項２２】
カメラであることを特徴とする請求項２０に記載の撮像装置。
【請求項２３】
前記撮像装置は、３００ｎｍから１５００ｎｍの波長の範囲で画像データを取得し、
前記撮像装置は、前記波長の範囲の光で前記被写体を照明する照明ユニットを有することを特徴とする請求項２０に記載の撮像装置。
【請求項２４】
前記被写体を部分コヒーレント照明する照明光学系を有することを特徴とする請求項２０に記載の撮像装置。
撮像装置。
【請求項２５】
請求項１乃至１７のうちいずれか１項に記載の画像処理方法を実行する演算手段を有する画像処理装置。
【請求項２６】
請求項２５に記載の画像処理装置にネットワークを介して接続されたネットワーク機器であって、
前記画像処理装置が被写体の画像データを処理する位置を指定する指定手段と、
前記画像処理装置に前記位置の情報を送信し、前記画像処理装置によって画像処理された前記画像データの情報を受信する通信手段と、
を有することを特徴とするネットワーク機器。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像処理方法、プログラム、画像処理装置、撮像装置およびネットワーク機器に関する。
【背景技術】
【０００２】
バーチャルスライドシステムでは、バーチャルスライドとよばれるデジタル撮像装置を用いて被写体のデジタル画像を取得する。一般に、医療分野では試料を光学素子（カバーグラス）で覆い、固定したもの（プレパラート）が被写体として用いられる。バーチャルスライドは顕微光学系、撮像素子、情報処理装置によって構成され、プレパラートをデジタル画像化し、データを保存する。
【０００３】
医師は撮像素子で取得した画像を見て診断を行うが、ＣＣＤ、ＣＭＯＳ、ＰＭＴなどの撮像素子では画像取得時にランダムなノイズが生じるため、取得した画像から効率的にノイズを低減することが求められる。特許文献１は、画像に離散コサイン変換を施し、変換画像の高周波成分を抑制するノイズ低減方法を提案している。特許文献２は、光軸方向の複数の異なる位置における画像を取得し、画像間の類似点からノイズ成分を処理し、合焦位置での画像のノイズを低減する方法を開示している。非特許文献１は、光軸方向の複数の異なる位置における複数枚の画像に対して周波数フィルタを施すことで、高解像度の画像を得る方法を提案している。その他の従来技術としては非特許文献２がある。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開平１０−３０４３６１号公報
【特許文献２】特開２００７−２７４０６４号公報
【非特許文献】
【０００５】
【非特許文献１】ジェイ・ビー・シバリタ（J.-B. Sibarita）著、ジェン・リートドルフ（Jens Rietdorf）編「マイクロスコーピー・テクニークス（Microscopy Techniques）」（米国）スプリンガー（Springer）、２００５年、ｐ．２０１−２４４
【非特許文献２】ミン・グ（Min Gu）著「プリンシプルズ・オブ・スリーディメンジョナル・イメージング・イン・コンフォーカル・マイクロスコープス（Principles of Three-Dimensional Imaging in Confocal Microscopes）」（米国）ワールド・サイエンティフィック（World Scientific）、１９９６年
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、特許文献１の方法は、低周波領域に存在するノイズを低減することはできない。特許文献２の方法は、光軸方向の位置によらず、ノイズが一定、もしくは光軸方向の位置に対して規則的にノイズが出現することが前提としているため、ランダムに発生するノイズを低減することは困難である。また、非特許文献１に記載は、光学系によってボケた画像を復元することが目的であり、ノイズの低減を図ったものではない。
【０００７】
本発明では、撮像素子によって取得された画像から低周波領域においてランダムに発生するノイズを効果的に低減することが可能な画像処理方法、プログラム、画像処理装置、撮像装置およびネットワーク機器を提供することを例示的な目的とする。
【課題を解決するための手段】
【０００８】
本発明の画像処理方法は、撮像光学系の光軸方向の複数の異なる位置において前記撮像光学系と撮像素子を有する撮像装置を用いて撮像および統合された被写体の３次元の画像データに含まれるノイズを、演算手段を利用して低減する画像処理方法であって、前記３次元の画像データに対して、前記光軸方向に対して周波数変換を行い、３次元の変換画像データを算出する画像変換ステップと、前記３次元の変換画像データの絶対値を特定の周波数領域で小さくし、３次元の変調画像データを算出する画像変調ステップと、前記３次元の変調画像データに対して、前記光軸方向に対して前記周波数変換に対応する逆周波数変換を行い、３次元の逆変換画像データを算出する逆画像変換ステップと、を有する。前記特定の周波数領域は０．０１ＮＡ２／λ＜｜ｆｚ｜の条件式で表わされる領域の一部であることを特徴とする。ここで、ｆｚは前記光軸方向の空間周波数、ＮＡは前記撮像光学系の被写体側の開口数、λは前記被写体からの光の波長である。
【発明の効果】
【０００９】
本発明によれば、撮像素子によって取得された画像から低周波領域においてランダムに発生するノイズを効果的に低減することが可能な画像処理方法、プログラム、画像処理装置、撮像装置およびネットワーク機器を提供することができる。
【図面の簡単な説明】
【００１０】
【図１】本発明のバーチャルスライドのブロック図である。（実施例１、２、３）
【図２】本発明による被写体の画像取得のフローチャートである。（実施例１、２、３）
【図３】本発明の原理を説明するための光学系である。（実施例１、２、３）
【図４】ＣＴＦのｆｙ＝０での断面である。（実施例１、２、３）
【図５】図５（ａ）はｆｙ＝０での断面で画像のスペクトルが非０となる領域を示す図、図５（ｂ）はｆｙ＝０での断面で共焦点光学系の画像のスペクトルが非０となる領域を示す図である。（実施例１、２、３）
【図６】本発明による画像処理方法のフローチャートである。（実施例１、２、３）
【図７】被写体の振幅透過率分布である。（実施例１）
【図８】光軸方向の異なる位置で取得した被写体の画像を示す図である。（実施例１）
【図９】図８に示す画像にノイズを付加した画像を示す図である。（実施例１）
【図１０】図１０（ａ）はスペクトルの絶対値のｆｙ＝０での断面図、図１０（ｂ）はで本発明のフィルタを施した後のスペクトルの絶対値のｆｙ＝０での断面図ある。（実施例１）
【図１１】本発明のフィルタを施した後の画像を示す図である。（実施例１）
【図１２】従来のフィルタを施した後のスペクトルの絶対値のｆｙ＝０での断面図である。
【図１３】従来のフィルタを施した後の画像を示す図である。
【図１４】光軸方向の位置ｚに対するＰＳＮＲを示すグラフである。（実施例１）
【図１５】改良されたフィルタのｆｙ＝０での断面図である。（実施例１）
【図１６】図１５に示すフィルタを用いた場合の光軸方向の位置ｚに対するＰＳＮＲを示すグラフである。（実施例１）
【図１７】光軸方向の位置を変えて取得した被写体の画像を示す図である。（実施例２）
【図１８】図１７に示す画像にノイズを付加した画像を示す図である。（実施例２）
【図１９】本発明のフィルタを施した後の画像を示す図である。（実施例２）
【図２０】従来のフィルタを施した後の画像を示す図である。
【図２１】光軸方向の位置ｚに対するＰＳＮＲを示すグラフである。（実施例２）
【図２２】光軸方向の位置ｚに対するＰＳＮＲを示すグラフである。（実施例２）
【発明を実施するための形態】
【００１１】
本発明は、撮像光学系を用いて撮像装置によって撮像された試料の画像を、光軸方向における任意の位置における画像のノイズを除去する画像処理方法に関する。画像処理方法は、コンピュータによって実行可能なプログラムとして具現化され、記録媒体などに記録されてもよい。
【００１２】
画像処理は、撮像装置において行われてもよいし、撮像装置の画像を格納する記憶装置（メモリ）に接続された画像処理装置（コンピュータ）が行ってもよい。撮像装置は画像処理手段として機能してもよい。また、取得した画像の表示は、撮像装置や画像処理装置に接続された表示手段（ディスプレイ）やこれらと（ＬＡＮ、ＷＡＮ、インターネットなどの）ネットワークを介して接続されたネットワーク機器において行ってもよい。この場合、画像処理方法はクラウドコンピューティングを利用することもでき、ネットワーク機器は画像データの受信を行うパーソナルコンピュータ（ＰＣ）であってもよい。あるいは、ネットワーク機器はタッチスクリーンなどの表示手段を備えた携帯端末（ｉＰａｄ（登録商標））、専用端末（ＰＤＡなど）であってもよい。これによって遠隔診断が可能になる。
【００１３】
図１は、本発明の実施に係るバーチャルスライドのブロック図である。バーチャルスライドは、撮像ユニット（撮像装置）１００、制御ユニット２００、情報処理ユニット（画像処理装置）４００から構成される。
【００１４】
制御ユニット２００は、搬送部２０１、制御部２０２から構成される。搬送部２０１は制御部２０２の命令に基づき被写体１０３を可動ステージ１０２上へ移動する。可動ステージ１０２は制御部２０２からの命令に従い、光軸方向に動くことが可能である。可動ステージ１０２は光軸と垂直方向に移動する機能を備えていてもよい。可動ステージ１０２を用いることで光軸方向の異なる位置で画像を取得することができる。
【００１５】
撮像ユニット１００は被写体１０３の画像を取得する装置であり、照明ユニット１０１、可動ステージ１０２、光学系（撮像光学系）１０４、撮像素子１０５からなる。
【００１６】
照明ユニット１０１は光源と照明光学系を有し、可動ステージ１０２に載置された被写体１０３を照明する。光学系１０４によって撮像素子１０５上に被写体の拡大光学像が形成される。
【００１７】
撮像素子１０５は被写体の拡大光学像を光電変換する光電変換手段であり、ＣＣＤ、ＣＭＯＳ、ＰＭＴなどから構成される。このため、撮像素子１０５による画像取得時にはランダムにノイズが発生する。撮像素子１０５から出力された電気信号は画像データとして情報処理ユニット４００に伝送される。
【００１８】
情報処理ユニット４００はコンピュータ（演算手段）４０１、画像処理部４０２、データ保存部（メモリ）４０３、表示部（表示手段）４０４によって構成される。
【００１９】
画像処理部４０２は撮像素子１０５から送られた画像データをデジタル信号に変換する。このデジタル信号を輝度信号と呼ぶ。画像処理部４０２は輝度信号に変換した画像データ対して、例えば、ノイズ除去や圧縮等の画像処理を施し、コンピュータ４０１に伝送する。送られた画像データをコンピュータ４０１はデータ保存部４０３に転送する。データ保存部４０３は送られた画像データを保存する。診断時は、コンピュータ４０１がデータ保存部４０３から画像データを読み出す。読み出された画像データに対してコンピュータ４０１は画像処理を施す。変換された画像データを表示部４０４へ転送し、画像が表示される。
【００２０】
コンピュータ４０１、画像処理部４０２、データ保存部４０３、表示部４０４及び制御部２０２は１つのコンピュータを構成してもよい。また、データはネットワーク４５０を介して接続された外部の不図示のサーバーに保存してもよく、この場合、複数の人がリモートでデータにアクセスすることができる。コンピュータ４０１は、ネットワーク４５０を介して様々なネットワーク機器に接続されている。これらのネットワーク機器は、ラップトップＰＣ４６０、デスクトップＰＣ４６２、タッチスクリーン機能を有する携帯端末４６４、ＰＤＡなどの専用端末４６６を含む。
【００２１】
ネットワーク機器は、演算手段、表示手段、指定手段、通信手段、記憶手段を有する。演算手段はコンピュータ（プロセッサ）から構成され、各部を制御や必要な演算を行う。表示手段は、ネットワーク機器４６０、４６４、４６６のようにネットワーク機器の筺体と一体であってもよいし、ネットワーク機器４６２のようにネットワーク機器に接続されたディスプレイ等であってもよい。指定手段は、ユーザーに光学系１０４の光軸方向の任意の位置ｚを指定させることを可能にするタッチスクリーン、キーボード、スタイラスペン、マウス等の入力手段を含む。通信手段は、ネットワーク４５０に接続されて画像処理装置に位置ｚの情報または光軸方向のＮ点（Ｎは２以上の整数）の異なる位置の画像データを送信して、画像処理装置から位置ｚにおける被写体１０３の画像データの情報を受信する。画像データの情報は、ｊｐｅｇなどの静止画像でもよいし、それを表す後述する輝度分布（画素値の分布）を表す関数Ｉ（ｘ，ｙ，ｚ）であってもよい。記憶手段は、位置ｚの指定を可能にするためのアプリケーションプログラムを格納したメモリである。ネットワーク機器４６０、４６４、４６６は通信手段を介して受信した画像データの情報に基づいて位置ｚの画像を表示する表示手段を更に有する。
【００２２】
図２（ａ）に画像取得のフローチャートを示す。ここで、「Ｓ」はステップ（工程）の略であり、これは他のフローチャートにおいても同様である。なお、被写体１０３の拡大像を撮像素子１０５、画像処理部４０２、コンピュータ４０１を用いて画像データとして取得することを画像取得と呼ぶ。
【００２３】
Ｓ１で被写体１０３を搬送部２０１を用いて可動ステージ１０２に設置し、Ｓ２で被写体１０３の画像を光軸上の複数の異なる位置で取得する。複数の異なる位置は合焦位置でなくてもよい。取得した複数の画像に対してＳ３で画像処理部４０２は画像処理を行う。画像処理では、光学系１０４の光軸方向の複数の異なる位置において撮像素子を用いて撮像（取得）および統合した被写体の３次元の画像データに含まれるノイズを、コンピュータ４０１を利用して低減する。ノイズが低減された画像をＳ４でコンピュータ４０１がデータ保存部４０３に保存する。データは一時保存であってもよい。Ｓ３のノイズ低減処理とＳ４の順番は逆でもよい。すなわち、保存されたデータに対してノイズ低減処理を行ってもよい。
【００２４】
次に、Ｓ２における被写体１０３の画像取得に関して詳細に説明する。図２（ｂ）に画像取得のフローチャートを示す。
【００２５】
始めに、Ｓ２０１でステージ位置を光軸方向（ｚ方向）に変える範囲と移動幅を指定する。直進光の伝搬方向へのフォーカスずれ量をｚｊとし、光学系１０４に対して撮像素子１０５と共役な位置を原点とする。ｚｊは光軸方向の位置及び光軸方向のステージ位置と同義である。Ｓ２０１で画像を取得するｚｊの範囲（最大値ｚｍａｘ〜最小値ｚｍｉｎ）と画像を取得するステージの移動幅（間隔）Δｚを指定する。移動幅Δｚは一定でなくてもよい。これらの値はＳ２０１を行う前に予め測定しておいた試料の表面形状から決めてもよいし、ユーザーが任意に設定してもよい。また予め設定しておいた既定の値を用いてもよい。
【００２６】
即ち、光学系１０４の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚｊ（１≦ｊ≦Ｎ）において、光学系１０４を用いて被写体１０３を撮像することによって被写体１０３の画像データを取得する。ｚｊの最大値をｚｍａｘ、ｚｊの最小値をｚｍｉｎとする。なお、ここでは、Ｎ点における撮影中に、光学系１０４の各レンズの位置を維持し、Ｎ点の１つを合焦位置とし、他をデフォーカス位置としているが、これに限定されるものではない。
【００２７】
次に、Ｓ２０２でステージをｚｍｉｎの位置まで移動させ、Ｓ２０３で被写体の画像を取得する。画像取得が終了した後、Ｓ２０４でステージをΔｚだけ移動させ、Ｓ２０５で画像を取得する。Ｓ２０４からＳ２０５をステージ位置がｚｍａｘ以上となるまで繰り返す。以上により、光軸方向の位置を変えた複数の画像データが取得される。ステージの移動方法はこれに限らない。例えば、Ｓ２０２でステージをＺｍａｘに移動した後、前述と逆方向に−Δｚずつ移動させて複数枚の画像を取得してもよい。この結果、光学系１０４の光軸方向に間隔Δｚで決定されたＮ点の異なる位置ｚｊ（１≦ｊ≦Ｎ）において、光学系１０４を用いて被写体１０３を撮像することによって被写体１０３のＮ枚の画像データを取得する。ｚｊの最大値がｚｍａｘ、ｚｊの最小値がｚｍｉｎである。
【００２８】
本実施形態では、取得した複数の位置での画像に対して画像処理部４０２がノイズ低減を行う。その原理を説明する。撮像光学系によって撮像素子面上に形成される光学像の強度分布を与える式を説明する。
【００２９】
原理を説明するための光学系を図３に示す。図３（ａ）は、被写体５０１の像が撮像光学系５０２によって撮像素子５０３上に結像している場合の光学系を示す。図３（ｂ）は、図３（ａ）の場合に対して被写体５０１が光軸方向にｚだけずれた場合の光学系を示す。
【００３０】
始めに照明光が垂直入射の完全コヒーレント光である場合を考える。被写体５０１の像が撮像素子５０３上に結像している場合の結像式は数式１で表される。
【００３１】
【数１】
2014096640.tif 000003
【００３２】
Ｉ（ｘ，ｙ）は撮像素子５０３上での像強度分布、Ｏ（ｆｘ，ｆｙ）は物体からの回折光分布である。また、Ｐ（ｆｘ，ｆｙ）は光学系の瞳関数、ｘ、ｙは像面での光軸に垂直な直交座標、ｆｘ、ｆｙはそれぞれｘ方向の空間周波数、ｙ方向の空間周波数、はフーリエ変換を表す。の下付き添え字の２は光軸方向に直交する２方向であるｘとｙ方向に２次元フーリエ変換（周波数変換）を行うことを示し、上付き添え字の−１は逆フーリエ変換（逆周波数変換）を行うことを示す。図３（ｂ）に示すように、被写体５０１が光軸方向にｚだけ移動した場合の結像式は数式２で表される。
【００３３】
【数２】
2014096640.tif 000004
【００３４】
ここでＰ（ｆｘ，ｆｙ；ｚ）は数式３で表され、被写体のフォーカスずれと等価なデフォーカス収差を加えた瞳関数である。
【００３５】
【数３】
2014096640.tif 000005
【００３６】
λは照明光の波長、ＮＡは撮像光学系の被写体側の開口数、ｋｚ０は照明光の波数のｚ方向成分、ｋｚは回折光の波数のｚ方向成分である。ｆｒ＝（ｆｘ２＋ｆｙ２）１／２である。簡単のため倍率は１倍とする。簡単のため倍率は１倍とする。撮像光学系の瞳面の中心に円形の遮蔽部があるとし、その半径が瞳の半径のε倍であるとする。本説明では瞳面の中心付近が遮光されている場合を想定しているが、これに限ることはない。瞳面の中心付近の光強度が周囲の光強度よりも弱くなっていれば同様の効果が期待できる。また、遮蔽部は撮像光学系の瞳面に配置される必要はない。円形である必要もない。瞳面での光強度分布が中心付近で弱くなるように設計されたフィルタが撮像光学系の光路中に配置されていればよい。また、中心遮蔽がない場合は数式３の波括弧を省略した式を考えればよい。
【００３７】
波数のｚ方向成分は数式４及び数式５で表される。
【００３８】
【数４】
2014096640.tif 000006
【００３９】
【数５】
2014096640.tif 000007
【００４０】
ｆｘ０、ｆｙ０はそれぞれ入射光のｘ、ｙ方向の波数成分を２πで割った値である。ｆｒ０＝（ｆｘ０２＋ｆｙ０２）１／２である。今は垂直入射を考えているためｆｘ０＝ｆｙ０＝０である。またｃｉｒｃ関数は数式６で表される。
【００４１】
【数６】
2014096640.tif 000008
【００４２】
光軸方向の位置を変えて取得した画像データは数式２のＩ（ｘ，ｙ，ｚ）のｚを変えて取得していることに相当する。数式２では被写体が光軸方向に移動した場合を想定しているが、撮像素子を光軸方向に動かした場合も同様に議論することができる。
【００４３】
ここで、Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域を考える。Ｉ（ｘ，ｙ，ｚ）のスペクトルとはＩ（ｘ，ｙ，ｚ）を３次元フーリエ変換して得られる３次元データである。Ｉ（ｘ，ｙ，ｚ）に対して３次元フーリエ変換を施す。すると数式２は数式７に変形される。
【００４４】
【数７】
2014096640.tif 000009
【００４５】
ここで、ＣＴＦ（ｆｘ，ｆｙ，ｆｚ）は数式８で表され、瞳関数をｚに対してフーリエ変換したものであり、撮像光学系の特性を表す関数である。は畳み込み積分を示し、添え字の３はｆｘ、ｆｙ、ｆｚの３つの座標に対して畳み込み積分を行うことを表す。上付きの＊は複素共役を示す。の下付き添え字の３はｘ、ｙ、ｚの３つの座標に対する３次元フーリエ変換を行うことを示す。
【００４６】
【数８】
2014096640.tif 000010
【００４７】
ＣＴＦ（ｆｘ，ｆｙ，ｆｚ）のｆｙ＝０での断面を図４に示す。数式８よりＣＴＦはｃｉｒｃ関数とδ関数の積である。今は垂直入射を考えているため、ｆｒ０＝０の場合を考える。δ関数の特性からＣＴＦは、ｆｚ＋１／λ−（１／λ２−ｆｒ２）１／２＝０を満たす曲面でのみ非０となる。これは空間周波数空間で中心がｆｚ＝−１／λ、半径が１／λの球面を表している。更に、ｃｉｒｃ関数の特性からＣＴＦ（ｆｘ，ｆｙ，ｆｚ）はＮＡε／λ≦ｆｒ≦ＮＡ／λでのみ非０となる。よって、ＣＴＦ（ｆｘ，ｆｙ，ｆｚ）が非０となるのは図４の太線部分のみとなる。簡単のため中心遮蔽がない場合を図示している。
【００４８】
Ｏ（ｆｘ，ｆｙ）は一般に空間周波数空間の全域で非０であるため、Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域はＣＴＦ（ｆｘ，ｆｙ，ｆｚ）が非０となる領域によって決定される。Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域を明確にするため、Ｏ（ｆｘ，ｆｙ）＝１という場合を考える。すると、Ｉ（ｘ，ｙ，ｚ）のスペクトルはＣＴＦ（ｆｘ，ｆｙ，ｆｚ）の自己相関で表される。畳み込み積分の結果、Ｉ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域は図５（ａ）のハッチング領域となる。簡単のため、近軸近似を導入し中心遮蔽がないとする。すなわちＣＴＦ（ｆｘ，ｆｙ，ｆｚ）を数式９で近似する。
【００４９】
【数９】
2014096640.tif 000011
【００５０】
数式９を数式７に代入し、Ｏ（ｆｘ，ｆｙ）＝１とおいて計算を実行すると、数式７は数式１０に変形される。
【００５１】
【数１０】
2014096640.tif 000012
【００５２】
ただし、数式１１以外の領域では０である。
【００５３】
【数１１】
2014096640.tif 000013
【００５４】
すなわち、Ｉ（ｘ，ｙ，ｚ）のスペクトルは数式１１が成立する領域でのみ非０となる。
【００５５】
また、数式１１より、Ｉ（ｘ，ｙ，ｚ）のスペクトルのｆｚ方向の広がりは数式１２で示される。
【００５６】
【数１２】
2014096640.tif 000014
【００５７】
部分コヒーレント結像系を考える場合では数式２の代わりに数式１３を用いればよい。
【００５８】
【数１３】
2014096640.tif 000015
【００５９】
Ｓ（ｆｘ０，ｆｙ０）は有効光源の強度分布である。有効光源とは照明光学系に挿入される光学素子によって、被写体がない場合に結像光学系の瞳面に形成される光強度分布である。数式１１における瞳関数Ｐは数式３の瞳関数をｆｘ０及びｆｙ０だけ平行移動させたものであるため、数式１１における瞳関数Ｐから求められるＣＴＦ（ｆｘ，ｆｙ、ｆｚ）の形状は図４に示した太線と変わらない。数式１３における被積分関数の絶対値２乗の項に対して、先の議論と同様にスペクトルが非０となる領域が議論できる。従って、部分コヒーレント結像系においてもＩ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域は数式１１で表される。
【００６０】
インコヒーレント結像系はＳ（ｆｘ０，ｆｙ０）＝１という場合に等価なため、インコヒーレント結像系においてもＩ（ｘ，ｙ，ｚ）のスペクトルが非０となる領域は数式１１で表される。
【００６１】
以上より、通常の光学系では照明光のコヒーレンスにかかわらずＩ（ｘ、ｙ、ｚ）のスペクトルが非０となる領域は数式１１で表される。
【００６２】
共焦点光学系で被写体の画像を取得した場合に、画像のスペクトルが非０となる領域を考察する。
【００６３】
共焦点光学系では、点照明、点検出を行う。被写体を透過した光、または被写体から反射した光を観察する際の結像式は数式１４となる。複号同順で−／＋が透過光を、＋／−が反射光を観察する共焦点光学系の場合である。
【００６４】
【数１４】
2014096640.tif 000016
【００６５】
染色された被写体からの蛍光光を観察する共焦点光学系では、照明光の吸収と蛍光の発生という過程によりインコヒーレント結像となるため結像式が数式１５となる。
【００６６】
【数１５】
2014096640.tif 000017
【００６７】
Ｏｆは被写体の発光効率のフーリエ変換である。共焦点光学系の場合も同様にＯ（ｆｘ，ｆｙ）＝１として考えることで、Ｉ（ｘ、ｙ、ｚ）のスペクトルが非０となる領域を導くことができる。Ｏ（ｆｘ，ｆｙ）＝１とすると数式１４及び数式１５は数式１６に変形される。これがＩ（ｘ、ｙ、ｚ）のスペクトルが非０となる領域である。数式１６は数式７の右辺でＯ（ｆｘ，ｆｙ）＝１とした式の自己相関を計算していることに相当する。自己相関によりＩ（ｘ、ｙ、ｚ）のスペクトルが非０となる領域は数式１１よりも広くなる。数式１６を計算した結果得られる共焦点光学系で取得した画像のスペクトルが非０となる領域をｆｙ＝０の断面で図５（ｂ）に示す。この計算は非特許文献２を参照されたい。図５（ｂ）よりＩ（ｘ、ｙ、ｚ）のスペクトルが非０となるのは数式１７と数式１８を満たす領域内である。
【００６８】
【数１６】
2014096640.tif 000018
【００６９】
【数１７】
2014096640.tif 000019
【００７０】
【数１８】
2014096640.tif 000020
【００７１】
以上のように、Ｉ（ｘ，ｙ，ｚ）のスペクトルは、光学系によって定まる領域でのみ非０となる。もし、所定の領域以外で値が存在すれば、ノイズなどに起因した真の画像には存在しないスペクトル成分が画像に混在していることを意味する。従って、取得した画像のスペクトルに対して所定の領域をそのまま残し、それ以外の領域の値を０にする、すなわち３次元のローパスフィルタを施せばノイズを低減することができる。
【００７２】
ここまで述べたノイズ低減方法（画像処理方法）を図６に示すフローチャートに従って説明する。
【００７３】
まず、Ｓ３０１で取得した複数の画像に対して３次元フーリエ変換を行い、画像の３次元のスペクトルを取得する。Ｓ３０１は、３次元の画像データに対して、少なくとも光軸方向（ここでは光軸方向と光軸に垂直な２つの直交方向）に対して周波数変換を行い、３次元の変換画像データを算出する画像変換ステップとして機能する。但し、周波数変換の例はフーリエ変換に限定されず、コサイン変換などであってもよい。
【００７４】
取得したスペクトルに対してＳ３０２で３次元のローパスフィルタを施す。Ｓ３０２は、３次元の変換画像データの絶対値を特定の周波数領域で小さくし、３次元の変調画像データを算出する画像変調ステップとして機能する。真の画像の情報は数式１１や数式１２、また数式１７及び数式１８で定まる特定の領域に存在する。そのため、ローパスフィルタとして数式１１や数式１２、また数式１７及び数式１８の条件式によって表される領域を非０とした周波数フィルタを用いれば、真の画像の情報を失うことなくノイズの低減が可能である。
【００７５】
例えば、数式１１であれば、３次元の変換画像データの絶対値を小さくすべき特定の周波数領域は、例えば、以下の条件式で表される領域の一部となる。
【００７６】
【数１９】
2014096640.tif 000021
【００７７】
あるいは、数式１２、数式１７、１８であれば、特有の周波数領域は次の条件式で表される領域の一部となる。
【００７８】
【数２０】
2014096640.tif 000022
【００７９】
【数２１】
2014096640.tif 000023
【００８０】
【数２２】
2014096640.tif 000024
【００８１】
フィルタをかけたスペクトルに対して、Ｓ３０３で３次元逆フーリエ変換を行う。Ｓ３０３は、３次元の変調画像データに対して、少なくとも光軸方向（ここでは光軸方向と光軸に垂直な２つの直交方向）に対して周波数変換に対応する逆周波数変換を行い、３次元の逆変換画像データを算出する逆画像変換ステップとして機能する。
【００８２】
以上により、ノイズが低減された画像を得ることができる。
【００８３】
周波数フィルタによって、画像のスペクトルの内、数式１１や数式１２、また数式１７及び数式１８で定まる特定の領域以外の成分はノイズとして除去される。これは、数式１１や数式１２、また数式１７及び数式１８で定まる特定の領域以外の領域が広ければ広いほど、ノイズ低減効果が高いことを意味する。画像のスペクトルは｜ｆｚ｜≦１／２Δｚの範囲で得られる。一方で、画像のスペクトルは、通常の光学系では数式１２、共焦点顕微鏡では数式１８で示される範囲に限定される。ノイズ除去効果を十分に得るためには、標準顕微鏡の場合１／２Δｚ≧ＮＡ２／２λ、共焦点顕微鏡の場合は１／２Δｚ≧ＮＡ２／λであることが望ましい。従って、画像を取得するステップ幅Δｚは通常の光学系の場合、数式１９を、共焦点顕微鏡の場合は数式２０を満たしていることが望ましい。
【００８４】
【数２３】
2014096640.tif 000025
【００８５】
【数２４】
2014096640.tif 000026
【００８６】
また、生体試料を透過した光や、人物、風景などから反射した光など一般的な被写体からの回折光分布は周波数が大きくなると値が小さくなる。そのため、光学系で定まる領域よりも更に限定したローパスフィルタを施してもよい。
【００８７】
図６に示すフローチャートではフーリエ変換を用いて３次元のスペクトルを取得したが、変換の方法はこれに限ることはない。例えば、コサイン変換を用いてもよい。ただし、コサイン変換で用いる周波数はフーリエ変換の場合に比べて２倍異なるので、数式１１、数式１２、数式１７、数式１８で示された条件式を２倍だけ変更する必要がある。すなわち、上記４つの数式に現れるｆｒ及びｆｚをｆｒ／２及びｆｚ／２に変更した数式がコサイン変換の場合に用いる条件式となる。また、コサイン変換は、基底の周期がフーリエ変換の場合に比べて２倍大きくなるため、より高い精度でノイズを低減することができる。
【実施例１】
【００８８】
以下、本発明の実施例１におけるノイズ低減手法について説明する。
【００８９】
シミュレーションによって得られた被写体の画像を示す。撮像光学系の被写体側の開口数ＮＡを０．７、結像倍率を４０倍、照明光を波長が５５０ｎｍで単色、撮像素子の画素サイズを４μｍ、同画素の面積開口率を５０％とする。照明は部分コヒーレント照明とし、有効光源の形状を円形とする。照明光学系の開口数をＮＡｉとした時、ＮＡｉ／ＮＡ＝０．７を仮定する。
【００９０】
被写体の振幅透過率を図７に示す。被写体として幅０．５μｍのラインチャートを用いた。計算上では被写体を２次元の物体と仮定している。もちろん、被写体に厚みがあっても本手法は適用可能である。光軸方向の位置を−４．０μｍから４．０μｍまで０．２５μｍごとに変化させて画像の輝度分布を計算した。図８に光軸方向の位置を変えて得られる画像の一例を示す。画像の軸は撮像素子面上での位置を示す。
【００９１】
画像にランダムノイズを付加する。ランダムノイズとして１０％ガウシアンノイズを仮定する。ガウシアンノイズとは数式２１に従うノイズである。
ノイズ＝ａ（１＋Ｇａｕｓｉａｎ（σ，ｍ））  （２１）
ここで、Ｇａｕｓｉａｎ（σ，ｍ）とは数ｉが出現する頻度Ｐ（ｉ）が標準偏差σ、平均ｍのガウス分布に従うガウス乱数である。ａは定数である。１０％ガウシアンノイズとは式のａが画像の最大輝度の１０％の値、σが１、ｍが０である場合を指す。光軸方向の位置を変えて取得した一連の画像にノイズを付加した後の画像の一例を図９に示す。図８と比較すると、画像全体にノイズが付加されていることが分かる。
【００９２】
次に、本実施例の画像処理方法について説明する。
【００９３】
光軸方向の位置ｚでの画像の輝度分布をＩ（ｘ，ｙ，ｚ）とする。Ｉ（ｘ，ｙ，ｚ）に対してｘ、ｙ、ｚに関する３次元の離散フーリエ変換を施しＩ’（ｆｘ，ｆｙ，ｆｚ）を取得する。図１０（ａ）に得られたＩ’（ｆｘ，ｆｙ，ｆｚ）の絶対値をｆｙ＝０での断面で示す。Ｉ’（ｆｘ，ｆｙ，ｆｚ）は図５（ａ）に示す領域で値が大きく、それ以外の領域ではランダムな信号が存在している。このランダムな信号がノイズ成分と考えられる。よって、図５（ａ）に示す領域を非０とし、それ以外の領域を０としたフィルタとＩ’（ｆｘ，ｆｙ，ｆｚ）との要素ごとの積を取ることでノイズを低減する。フィルタとＩ’（ｆｘ，ｆｙ，ｆｚ）との要素ごとの積を取ることをフィルタリングと呼ぶ。
【００９４】
用いるフィルタとして数式１１を満たす領域を１、それ以外を０とするフィルタを採用する。もちろんフィルタはこれに限ることはない。
【００９５】
フィルタリング後のスペクトルをＩ’’（ｆｘ，ｆｙ，ｆｚ）とする。図１０（ｂ）にＩ’’（ｆｘ，ｆｙ，ｆｚ）の絶対値をｆｙ＝０での断面で示す。図５（ａ）に示した領域のみで値が有限で、それ以外の領域では値が０となるスペクトルが得られる。フィルタリングによってＩ’（ｆｘ，ｆｙ，ｆｚ）にみられていた図５（ａ）で示した領域以外にあったランダムな信号がなくなっている。
【００９６】
Ｉ’’（ｆｘ，ｆｙ，ｆｚ）をｆｘ、ｆｙ、ｆｚに対して３次元逆フーリエ変換する。その結果得られた一連の画像の一例を図１１に示す。図９と比較すると、ノイズが少なくなっていることが分かる。
【００９７】
本実施例を従来技術と比較する。従来の画像処理方法では１枚１枚の画像に対して２次元のローパスフィルタを施している。すなわち、スペクトルＩ’（ｆｘ，ｆｙ，ｆｚ）の（ｆｘ２＋ｆｙ２）１／２＞２ＮＡ／λの領域の値を０とする。従来技術でフィルタリングした後のスペクトルをＩ’’２（ｆｘ，ｆｙ，ｆｚ）とする。図１２にＩ’’２（ｆｘ，ｆｙ，ｆｚ）の絶対値をｆｙ＝０での断面で示す。図１０（ａ）と比較すると図１２ではｆｘ＞２ＮＡ／λの領域にあったランダムな信号が除去されているが、図１０（ｂ）と比較するとｆｚの高周波領域でランダムな信号が残存する。
【００９８】
Ｉ’’２（ｆｘ，ｆｙ，ｆｚ）をｆｘ、ｆｙ、ｆｚに対して３次元逆フーリエ変換を行い、得られた画像を図１３に示す。図９と比較すると取得画像からノイズが低減されていることが分かるが、図１１と比較すると、本実施例の方がノイズ低減効果が高いことが分かる。
【００９９】
画像の違いを定量化するためにノイズを付加する前の画像を基準としたＰＳＮＲ（Ｐｅａｋ  Ｓｉｇｎａｌ  ｔｏ  Ｎｏｉｓｅ  Ｒａｔｉｏ）を計算する。ＰＳＮＲとは数式２２と数式２３を用いて計算される値で、画像の類似性を定量化したものである。０以上の値をとり、値が大きいほど画像の類似度が高いことを示す。今回の場合、ＰＳＮＲが大きいほどノイズ低減の効果が大きいことになる。
【０１００】
【数２５】
2014096640.tif 000027
【０１０１】
【数２６】
2014096640.tif 000028
【０１０２】
ここで、Ｉｒｓｔ（ｘ，ｙ，ｚ）はノイズ低減処理を施した画像である。数式２２は８−ｂｉｔ画像に対応した式である。
【０１０３】
図１４に位置ｚにおけるＰＳＮＲを示す。点線がノイズを付加した画像に対するＰＳＮＲ、一点破線が従来の画像処理方法で得られた画像に対するＰＳＮＲ、実線が本発明の画像処理方法で得られた画像に対するＰＳＮＲである。すべてのｚにおいて実線が最も値が大きく、一点破線よりも値が大きい。これは本実施例が従来技術よりもノイズ低減効果が高いことを示している。
【０１０４】
更にノイズ低減効果の高いフィルタについて説明する。
【０１０５】
一般的な被写体からの回折光分布は空間周波数が大きくなると値が小さくなる。そのため、空間周波数の高いスペクトル成分はノイズに起因する信号よりも値が小さくなる場合がある。その結果、数式１１を満たす領域であっても画像を劣化させる要因が含まれている場合がある。よって、数式１１で定めたフィルタよりも更に領域の狭いローパスフィルタを用いることでノイズ低減効果を高めることができる。ただし、この場合、物体に起因する回折光成分まで低減してしまうため、画像の劣化が同時に起こる。フィルタをどこまで小さくするかは、処理後の画像からユーザーが画像の劣化を許容できる範囲で定める必要がある。
【０１０６】
解析的に求めた領域よりも狭いフィルタを作成する。新たなフィルタとして数式１１、数式２４、数式２５の３つすべてを満たす領域を非０とするフィルタを採用する。ただし、ｆｒは光軸に垂直な方向の空間周波数であり、ｆｚは光軸方向の空間周波数である。
Ｃｒ≧｜ｆｒ｜  （２４）
Ｃｚ≧｜ｆｚ｜  （２５）
Ｃｒ及びＣｚはユーザーが任意に設定可能な定数である。例えば、Ｃｒ＝０．９ＮＡ／λ、Ｃｚ＝０．５ＮＡ２／λ（数式１２’）として作成したフィルタを図１５にｆｙ＝０の断面で示す。図５（ａ）の場合と比べてｆｘ方向に非０の領域が狭くなっていることが分かる。
【０１０７】
このフィルタを用いてノイズ低減を行った画像に対するＰＳＮＲを図１６に示す。太線がＣｒ＝０．９ＮＡ／λとしたフィルタを用いた場合であり、分かり易くするため矢印で示す。細い実線が数式１１のみから定めたフィルタを用いた場合、点線がフィルタをかけない場合である。太線は細い実線よりも値が大きくなっている。よって、数式２４と数式２５の条件を加えたフィルタを用いると、よりノイズ低減効果を高めることができる。
【０１０８】
適切なＣｒとＣｚは被写体によって異なる。例えば、Ｃｒ＝０．０１ＮＡ／λ、好ましくはＣｒ＝０．９ＮＡ／λである。同様に、Ｃｚ＝０．０１ＮＡ２／λ、好ましくはＣｚ＝０．５ＮＡ２／λである。従って、３次元の変換画像データの絶対値を小さくすべき特定の周波数領域は、例えば、｜ｆｚ｜＞０．０１ＮＡ２／λまたは｜ｆｚ｜＞０．５ＮＡ２／λの条件式で表される領域の一部となる。あるいは、特定の周波数領域は、例えば、｜ｆｒ｜＞０．０１ＮＡ／λまたは｜ｆｚ｜＞０．９ＮＡ／λの条件式で表される領域の一部となる。
【０１０９】
本実施例では被写体として幅０．５μｍのラインチャートを用いたが、被写体はこれに限ることはない。例えば、生体の組織片や半導体の基盤、または風景や人物でもよい。
【０１１０】
本実施例では近軸近似を仮定した数式１０を利用してフィルタを設定したが、近軸近似を用いずに求めたフィルタを用いてもよい。例えば、数式７をコンピュータで計算し、非０の領域を求めてフィルタを設定しても良い。
【０１１１】
本実施例ではフーリエ変換を使用したが、これに限ることはない。例えば、コサイン変換を用いてもよい。ただし、コサイン変換では、用いる周波数がフーリエ変換の場合に比べて２倍異なるので、数式１１、数式２４、数式２５を用いる際に、これらの式に現れるｆｒ及びｆｚをｆｒ／２及びｆｚ／２に変更する必要がある。
【０１１２】
例えば、周波数変換がコサイン変換の場合には、数式１１’、数式１２’、数式１７’、数式１８’は次式のようになる。
【０１１３】
【数２７】
2014096640.tif 000029
【０１１４】
【数２８】
2014096640.tif 000030
【０１１５】
【数２９】
2014096640.tif 000031
【０１１６】
【数３０】
2014096640.tif 000032
【０１１７】
本実施例では波長を５５０ｎｍの単色光と仮定したが、これに限ることはない。例えば、発光ダイオードを複数利用することで多色光による照明及び撮像を行ってもよいし、ハロゲンランプ等を利用することで、波長に広がりを持った照明光で照明及び撮像を行ってもよい。この時利用する波長は生体イメージングや一般的な写真撮影で利用される波長領域、すなわち３００ｎｍ〜１５００ｎｍに含まれていることが望ましい。この場合、撮像素子も前記波長領域に対して応答することが望ましい。また、撮像装置は前記波長領域の分光データを取得してもよい。用いる光が単色光でない場合、フィルタの設定に用いる波長は用いた光の波長の内、最も短い波長を用いれば、画像を劣化させることなくノイズ低減を実現することができる。
【０１１８】
本実施例で撮像装置の結像倍率を４０倍としたが、これに限ることはない、例えば、縮小光学系でもよい。撮像光学系が縮小光学系で照明光学系がインコヒーレント系の場合は一般的なカメラによる撮像を行っている場合に相当する。すなわち、一般的なカメラを用いた写真撮影においても、本発明は実施可能である。
【０１１９】
本実施例では収差がない場合を仮定しているが、これに限ることはない。
【０１２０】
本実施例では取得した画像に対して像ずれの補正等の画像処理を行っていないが、これらの画像処理を施してもよい。
【実施例２】
【０１２１】
以下、本発明の実施例２におけるノイズ低減手法について説明する。実施例２では共焦点光学系で被写体の蛍光を観察した場合を想定する。
【０１２２】
図２（ｂ）に示すフローチャートに従って取得した被写体の画像をシミュレーションによって示す。撮像光学系の被写体側開口数を０．７、照明光を波長が５５０ｎｍで単色、被写体からの蛍光を波長が５５０ｎｍで単色、ｘ、ｙ方向のサンプリングピッチを物体上で１００ｎｍとする。被写体の発光効率の分布を図７に示す分布で定めた。計算上では被写体を２次元の物体と仮定している。もちろん、被写体に厚みがあっても本手法は適用可能である。光軸方向の位置は−４．０μｍから４．０μｍまで０．２５μｍごとに変化させて画像の輝度分布を計算した。図１７に光軸方向の位置を変えて得られる画像の一例を示す。画像の軸は物体面上での値を示す。
【０１２３】
画像に対して１０％ガウシアンノイズを仮定する。ノイズを付加した後の画像の一例を図１８に示す。図１７と比較すると、画像全体にノイズが付加されていることが分かる。
【０１２４】
ノイズ低減手法について説明する。
【０１２５】
位置ｚでの画像の輝度分布をＩ（ｘ，ｙ，ｚ）とする。Ｉ（ｘ，ｙ，ｚ）に対してｘ、ｙ、ｚに関する３次元の離散フーリエ変換を施しＩ’（ｆｘ，ｆｙ，ｆｚ）を取得する。数式１７及び数式１８を満たす領域を非０としたローパスフィルタを施す。もちろん、フィルタはこれに限ることはない。
【０１２６】
フィルタリングによって得られる画像の一例を図１９に示す。図１８と比較すると、ノイズが少なくなっていることが分かる。また、図２０に従来技術によって得られる画像の一例を示す。従来技術では１枚１枚の画像に対してローパスフィルタを施すため、数式１７を満たす領域を非０とするフィルタを用いた。図１９と図２０を比較すると図１９の方がノイズが少なく、本実施例の方がノイズ低減効果が高いことが分かる。
【０１２７】
画像の違いを定量化するためにノイズを付加する前の画像とのＰＳＮＲを計算する。図２１に位置ｚとＰＳＮＲを示す。点線がノイズを付加した画像に対するＰＳＮＲ、一点破線が従来のノイズ低減手法で得られた画像に対するＰＳＮＲ、実線が本発明のノイズ低減手法で得られた画像に対するＰＳＮＲである。実線が最も値が大きく、一点破線よりも値が大きい。これは本実施例が従来技術よりもノイズ低減効果が高いことを示している。
【０１２８】
自然光の回折光分布は一般に空間周波数が大きくなると値が小さくなる。そのため、空間周波数の高いスペクトル成分はノイズに起因する信号よりも値が小さくなる場合がある。よって、数式１７及び数式１８で定めたフィルタよりも更に領域の狭いローパスフィルタを用いることでノイズ低減効果を高めることができる。一例として、Ｃｒ＝１．１ＮＡ／λ、Ｃｚ＝ＮＡ２／λとした時に数式２４と数式２５を満たす領域を非０としたフィルタを用いた結果を示す。図２２にノイズ低減後の画像に対するＰＳＮＲを示す。太線がＣｒ＝１．１ＮＡ／λとしたフィルタを用いた場合であり、分かり易くするため矢印で示す。先の結果よりも、ノイズ低減効果が高くなっていることが分かる。
【０１２９】
この場合、物体に起因する回折光成分まで低減してしまうため、画像の劣化が同時に起こる。フィルタをどこまで小さくするかは、処理後の画像からユーザーが画像の劣化を許容できる範囲で定める必要がある。
【０１３０】
本実施例では被写体として幅０．５μｍのラインチャートを用いたが、被写体はこれに限ることはない。例えば、生体の組織片や半導体の基盤、または風景や人物でもよい。
【０１３１】
本実施例では蛍光を観察する共焦点光学系を想定しているが、反射光を観察する共焦点光学系をおよび透過光を観察する共焦点光学系においても同様の方法でノイズ低減を行うことができる。
【実施例３】
【０１３２】
以下、画像処理装置にネットワーク４５０を介して接続されたネットワーク機器４６０、４６２、４６４、４６６（以下、「４６０」で代表させる）を利用してノイズが除去された画像を取得する場合の動作について説明する。この場合、ユーザーは、既にノイズが除去された状態で不図示のサーバーに保存されている画像データの送信を依頼することができる。あるいは、ユーザーは、不図示のサーバーに保存されている、所定位置（例えば、合焦位置またはデフォーカス位置）の画像データを含むＮ枚の画像データを指定してノイズが除去された所定位置の画像データの送信を依頼する。あるいは、ユーザーは、所定位置の画像データを含むＮ枚の画像データを送信してノイズが除去された所定位置の画像データの送信を依頼する。
【０１３３】
まず、ユーザーはネットワーク機器４６０を操作してネットワーク４５０を介して画像処理装置または不図示のサーバーにアクセスする。不図示のサーバーには不図示の記憶装置（ストレージ）が接続されている。この場合、画像処理装置は画像データを保存しておらず、処理すべき画像データは記憶装置に保存され、画像処理装置は、処理すべき画像データを記憶装置から取得する。ネットワーク機器４６０が画像処理装置にアクセスする場合、処理すべき画像データを保存する記憶装置は画像処理装置に直接に、または、ＬＡＮ等のネットワークを介して接続されている。以下、ネットワーク機器４６０は、画像処理装置である情報処理ユニット４００にアクセスするものとする。
【０１３４】
ネットワーク機器４６０は情報処理ユニット４００にアクセスすると、ユーザー名とパスワードの入力を求められ、これによって認証が行われる。
【０１３５】
認証されたユーザー（例えば、医師）は、例えば、患者のＩＤを入力することによってデータ保存部４０３に保存されている患者の情報にアクセスする。ここでは、データ保存部４０３が、Ｎを２以上の整数として、撮像光学系の光軸方向のＮ点の異なる位置ｚｊ（１≦ｊ≦Ｎ）において、その撮像光学系を有する撮像装置を用いて撮像（取得）した被写体（患者）のＮ枚の画像データを保存しているものとする。そこで、ユーザーは、Ｎ枚の画像データのサムネイルやファイル名から少なくとも１枚をノイズ除去対象としての画像データとして選択する。次に、ユーザーは、ノイズ除去のボタンをクリックする。これに応答して、上述したノイズ除去（画像処理）がコンピュータ４０１と画像処理部４０２によって行われる。
【０１３６】
その後、ユーザーは、この画像をネットワーク上で閲覧するか、ダウンロードして保存するか、別の画像を選択するか、をメニューに従って選択する。これによって、ユーザーはノイズ除去された画像データを取得することができる。
【産業上の利用可能性】
【０１３７】
本発明の画像処理方法は、デジタル顕微鏡やデジタルカメラの用途に適用することができる。
【符号の説明】
【０１３８】
１０３…被写体、１０４…光学系（撮像光学系）、１０５…撮像素子、４００…情報処理ユニット（画像処理装置）、４０１…コンピュータ（演算手段）、４０２…画像処理部
【図２】
2014096640.tif 000042
【図３】
2014096640.tif 000043
【図４】
2014096640.tif 000044
【図５】
2014096640.tif 000045
【図６】
2014096640.tif 000046
【図１４】
2014096640.tif 000047
【図１６】
2014096640.tif 000048
【図２１】
2014096640.tif 000049
【図２２】
2014096640.tif 000050
【図１】
2014096640.tif 000051
【図７】
2014096640.tif 000052
【図８】
2014096640.tif 000053
【図９】
2014096640.tif 000054
【図１０】
2014096640.tif 000055
【図１１】
2014096640.tif 000056
【図１２】
2014096640.tif 000057
【図１３】
2014096640.tif 000058
【図１５】
2014096640.tif 000059
【図１７】
2014096640.tif 000060
【図１８】
2014096640.tif 000061
【図１９】
2014096640.tif 000062
【図２０】
2014096640.tif 000063