JP6039371B2 - 画像処理方法、プログラム、画像処理装置、および撮像装置 - Google Patents

画像処理方法、プログラム、画像処理装置、および撮像装置 Download PDF

Info

Publication number
JP6039371B2
JP6039371B2 JP2012245701A JP2012245701A JP6039371B2 JP 6039371 B2 JP6039371 B2 JP 6039371B2 JP 2012245701 A JP2012245701 A JP 2012245701A JP 2012245701 A JP2012245701 A JP 2012245701A JP 6039371 B2 JP6039371 B2 JP 6039371B2
Authority
JP
Japan
Prior art keywords
image
image processing
processing method
frequency
image data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2012245701A
Other languages
English (en)
Other versions
JP2014096640A (ja
JP2014096640A5 (ja
Inventor
輝 江口
輝 江口
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Canon Inc
Original Assignee
Canon Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Canon Inc filed Critical Canon Inc
Priority to JP2012245701A priority Critical patent/JP6039371B2/ja
Priority to US14/070,684 priority patent/US9014506B2/en
Publication of JP2014096640A publication Critical patent/JP2014096640A/ja
Publication of JP2014096640A5 publication Critical patent/JP2014096640A5/ja
Application granted granted Critical
Publication of JP6039371B2 publication Critical patent/JP6039371B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G02—OPTICS
    • G02B—OPTICAL ELEMENTS, SYSTEMS OR APPARATUS
    • G02B21/00—Microscopes
    • G02B21/36—Microscopes arranged for photographic purposes or projection purposes or digital imaging or video purposes including associated control and data processing arrangements
    • G02B21/365—Control or image processing arrangements for digital or video microscopes
    • G02B21/367—Control or image processing arrangements for digital or video microscopes providing an output produced by processing a plurality of individual source images, e.g. image tiling, montage, composite images, depth sectioning, image comparison
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration using non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20052—Discrete cosine transform [DCT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04N—PICTORIAL COMMUNICATION, e.g. TELEVISION
    • H04N23/00—Cameras or camera modules comprising electronic image sensors; Control thereof
    • H04N23/80—Camera processing pipelines; Components thereof
    • H04N23/81—Camera processing pipelines; Components thereof for suppressing or minimising disturbance in the image signal generation

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Optics & Photonics (AREA)
  • Image Processing (AREA)
  • Microscoopes, Condenser (AREA)
  • Studio Devices (AREA)

Description

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

Claims (27)

  1. 撮像光学系を介して被写体を撮像することで得られる３次元の画像データに含まれるノイズを低減する画像処理方法であって、
    前記画像データに対して、前記撮像光学系の光軸方向における周波数変換を行うことで、３次元の変換画像データを算出する画像変換ステップと、
    前記変換画像データの絶対値を第１の周波数領域において小さくすることで、３次元の変調画像データを算出する画像変調ステップと、
    前記変調画像データに対して、前記光軸方向における前記周波数変換に対応する逆周波数変換を行うことで、３次元の逆変換画像データを算出する逆画像変換ステップと、
    を有し、
    前記光軸方向における空間周波数をｆｚ、前記撮像光学系の前記被写体の側の開口数をＮＡ、前記被写体からの光の波長をλ、とするとき、前記第１の周波数領域は、
    ｜ｆｚ｜＞０．０１ＮＡ２／λ
    なる条件式で表わされる領域に含まれることを特徴とする画像処理方法。
  2. 前記画像データは、前記被写体の前記光軸方向における複数の異なる位置を撮像することで得られる、複数の画像を統合することにより取得されることを特徴とする請求項１に記載の画像処理方法。
  3. 前記周波数変換はフーリエ変換であることを特徴とする請求項１または２に記載の画像処理方法。
  4. 前記第１の周波数領域は、次の条件式で表される領域に含まれることを特徴とする請求項３に記載の画像処理方法。
  5. 前記第１の周波数領域は、次の条件式で表される領域に含まれることを特徴とする請求項３または４に記載の画像処理方法。
  6. 前記周波数変換はコサイン変換であることを特徴とする請求項１または２に記載の画像処理方法。
  7. 前記第１の周波数領域が、以下の条件式で表される領域であることを特徴とする請求項６に記載の画像処理方法。
  8. 前記第１の周波数領域が、以下の条件式で表される領域であることを特徴とする請求項６または７に記載の画像処理方法。
  9. 撮像光学系を介して被写体を撮像することで得られる３次元の画像データに含まれるノイズを、低減する画像処理方法であって、
    前記画像データに対して、前記撮像光学系の光軸方向および光軸に垂直な２つの直交方向における周波数変換を行うことで、３次元の変換画像データを算出する画像変換ステップと、
    前記変換画像データの絶対値を第１の周波数領域において小さくすることで、３次元の変調画像データを算出する画像変調ステップと、
    前記変調画像データに対して、前記光軸方向および前記２つの直交方向における前記周波数変換に対応する逆周波数変換を行うことで、３次元の逆変換画像データを算出する逆画像変換ステップと、
    を有し、
    前記光軸方向における空間周波数をｆｚ、前記撮像光学系の光軸に垂直な方向における空間周波数をｆｒ、前記撮像光学系の被写体側の開口数をＮＡ、前記被写体からの光の波長をλ、とするとき、前記第１の周波数領域は、以下の３つの条件式のうち少なくとも１つで表わされる領域に含まれることを特徴とする画像処理方法。

    ｜ｆｒ｜＞０．０１ＮＡ／λ
    ｜ｆｚ｜＞０．０１ＮＡ２／λ
  10. 前記画像データは、前記被写体の前記光軸方向における複数の異なる位置を撮像することで得られる、複数の画像を統合することにより取得されることを特徴とする請求項９に記載の画像処理方法。
  11. 前記第１の周波数領域は、更に以下の条件式で表わされる領域に含まれることを特徴とする請求項９または１０に記載の画像処理方法。
    ｜ｆｚ｜＞０．５ＮＡ２／λ
  12. 前記第１の周波数領域は、更に以下の条件式で表わされる領域に含まれることを特徴とする請求項９に記載の画像処理方法。
    ｜ｆｒ｜＞０．９ＮＡ／λ
  13. 前記周波数変換はフーリエ変換であることを特徴とする請求項９乃至１２のうちいずれか１項に記載の画像処理方法。
  14. 前記第１の周波数領域は、更に次の２つの条件式のうち少なくとも１つで表される領域に含まれることを特徴とする請求項１３に記載の画像処理方法。
  15. 前記周波数変換はコサイン変換であることを特徴とする請求項９乃至１２のうちいずれか１項に記載の画像処理方法。
  16. 前記第１の周波数領域は、更に以下の３つの条件式のうち少なくとも１つで表わされる領域に含まれることを特徴とする請求項１５に記載の画像処理方法。

    ｜ｆｒ｜＞０．０２ＮＡ／λ
    ｜ｆｚ｜＞０．０２ＮＡ２／λ
  17. 前記第１の周波数領域が、更に以下の２つの条件式のうち少なくとも１つで表される領域に含まれることを特徴とする請求項１５に記載の画像処理方法。
  18. 前記光軸方向の前記複数の異なる位置の間隔をΔｚとするとき、以下の条件式が満たされることを特徴とする請求項１乃至１７のうちいずれか１項に記載の画像処理方法。
  19. 前記画像変調ステップは前記変換画像データの前記第１の周波領域での値を０にするステップであることを特徴とする請求項１乃至１８のうちいずれか１項に記載の画像処理方法。
  20. 請求項１乃至１９のうちいずれか１項に記載の画像処理方法をコンピュータによって実行可能なプログラム。
  21. 請求項２０に記載のプログラムを記録した記録媒体。
  22. 請求項１乃至１９のうちいずれか１項に記載の画像処理方法を実行する演算装置を有する撮像装置。
  23. 顕微鏡であることを特徴とする請求項２２に記載の撮像装置。
  24. カメラであることを特徴とする請求項２２に記載の撮像装置。
  25. 前記撮像装置は、３００ｎｍから１５００ｎｍの波長の範囲で画像データを取得し、
    前記撮像装置は、前記波長の範囲の光で前記被写体を照明する照明ユニットを有することを特徴とする請求項２２に記載の撮像装置。
  26. 前記被写体を部分コヒーレント照明する照明光学系を有することを特徴とする請求項２２に記載の撮像装置。
    撮像装置。
  27. 請求項１乃至１９のうちいずれか１項に記載の画像処理方法を実行する演算手段を有する画像処理装置。
JP2012245701A 2012-11-07 2012-11-07 画像処理方法、プログラム、画像処理装置、および撮像装置 Expired - Fee Related JP6039371B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2012245701A JP6039371B2 (ja) 2012-11-07 2012-11-07 画像処理方法、プログラム、画像処理装置、および撮像装置
US14/070,684 US9014506B2 (en) 2012-11-07 2013-11-04 Image processing method, program, image processing apparatus, image-pickup optical apparatus, and network device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2012245701A JP6039371B2 (ja) 2012-11-07 2012-11-07 画像処理方法、プログラム、画像処理装置、および撮像装置

Publications (3)

Publication Number Publication Date
JP2014096640A JP2014096640A (ja) 2014-05-22
JP2014096640A5 JP2014096640A5 (ja) 2015-12-17
JP6039371B2 true JP6039371B2 (ja) 2016-12-07

Family

ID=50622444

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2012245701A Expired - Fee Related JP6039371B2 (ja) 2012-11-07 2012-11-07 画像処理方法、プログラム、画像処理装置、および撮像装置

Country Status (2)

Country Link
US (1) US9014506B2 (ja)
JP (1) JP6039371B2 (ja)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6362363B2 (ja) * 2014-03-10 2018-07-25 キヤノン株式会社 画像推定方法、プログラム、記録媒体および画像推定装置
JP6408893B2 (ja) * 2014-12-16 2018-10-17 オリンパス株式会社 ３次元位置情報取得方法及び３次元位置情報取得装置
CN104992414B (zh) * 2015-06-19 2017-11-24 中国科学技术大学 三维离子速度聚焦影像的处理方法
CN106204497B (zh) * 2016-07-20 2018-12-25 长安大学 一种基于smooth平滑曲线和拟合曲线的路面裂缝提取算法
US10181181B2 (en) * 2017-05-10 2019-01-15 Southwest Research Institute Denoising with three dimensional fourier transform for three dimensional images, including image sequences

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH10304361A (ja) * 1997-04-25 1998-11-13 Matsushita Electric Ind Co Ltd 画像符号化装置およびノイズ成分検出方法並びにノイズ除去方法
DE69914473D1 (de) * 1998-04-15 2004-03-04 Vincent Lauer Mikroskop für die erzeugung einer dreidimensionalen darstellung von einem objekt und bilder, die mittels dieses mikroskops erzeugt sind
JP3961729B2 (ja) * 1999-03-03 2007-08-22 株式会社デンソー 全焦点撮像装置
JP2001257932A (ja) * 2000-03-09 2001-09-21 Denso Corp 撮像装置
US7312875B2 (en) * 2003-04-23 2007-12-25 Ut-Battelle Llc Two-wavelength spatial-heterodyne holography
EP1801754B8 (en) * 2004-10-14 2011-09-14 Lightron Co. Ltd. Degradation information restoring method and device
US7619563B2 (en) * 2005-08-26 2009-11-17 Step Communications Corporation Beam former using phase difference enhancement
JP2007274064A (ja) * 2006-03-30 2007-10-18 Matsushita Electric Ind Co Ltd ノイズ低減装置
WO2008011510A2 (en) * 2006-07-21 2008-01-24 Zygo Corporation Compensation of systematic effects in low coherence interferometry
JP2009038439A (ja) * 2007-07-31 2009-02-19 Tokyo Institute Of Technology 空間フィルタ処理を行う撮像方法及びその撮像装置
JP2011053901A (ja) * 2009-09-01 2011-03-17 Ricoh Co Ltd 文書画像データ提供装置、文書画像データ提供システム、文書画像データ提供方法、文書画像データ提供プログラム、背景処理プログラム

Also Published As

Publication number Publication date
US20140126806A1 (en) 2014-05-08
JP2014096640A (ja) 2014-05-22
US9014506B2 (en) 2015-04-21

Similar Documents

Publication Publication Date Title
Luo et al. Pixel super-resolution using wavelength scanning
Han et al. Pixelation effect removal from fiber bundle probe based optical coherence tomography imaging
JP6039371B2 (ja) 画像処理方法、プログラム、画像処理装置、および撮像装置
JP4865930B2 (ja) 構造化された照射および均一な照射の両方を用いて光学的に切片化された画像を生成するためのシステムおよび方法
JP6112872B2 (ja) 撮像システム、画像処理方法、および撮像装置
US10269103B2 (en) Image processing apparatus, image processing method, and image processing system
Goud et al. Low cost digital holographic microscope for 3-D cell imaging by integrating smartphone and DVD optical head
Thapa et al. Comparison of super-resolution algorithms applied to retinal images
Huang et al. Real-time reference A-line subtraction and saturation artifact removal using graphics processing unit for high-frame-rate Fourier-domain optical coherence tomography video imaging
JP2015052663A (ja) 画像処理方法、画像処理装置、撮像装置およびプログラム
Watanabe et al. Graphics processing unit accelerated intensity-based optical coherence tomography angiography using differential frames with real-time motion correction
Liu et al. Honeycomb pattern removal for fiber bundle endomicroscopy based on a two-step iterative shrinkage thresholding algorithm
Pérez et al. Lightfield recovery from its focal stack
US20190261843A1 (en) Extended depth of field intraoral imaging apparatus
JP6661491B2 (ja) 画像処理装置および画像処理方法
JP6039236B2 (ja) 画像推定方法、プログラム、記録媒体、画像推定装置、および画像データの取得方法
Marks et al. Cone-beam tomography with a digital camera
JP6362363B2 (ja) 画像推定方法、プログラム、記録媒体および画像推定装置
Myakinin et al. A complex noise reduction method for improving visualization of SD-OCT skin biomedical images
Zhu et al. Phase retrieval for interferometry imaging from microlens array
Boroomand et al. Lateral resolution enhancement via imbricated spectral domain optical coherence tomography in a maximum-a-posterior reconstruction framework
JP2014002636A (ja) 画像推定方法、画像推定プログラムおよび画像推定装置
Reunov et al. Comparison of Methods for Removing Noise in an Image Obtained in a Specular Microscope at a Wavelength of 13.84 nm
LeBlanc et al. Joint camera blur and pose estimation from aliased data
Chen et al. An OCT image denoising method based on fractional integral

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20151102

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20151102

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20160707

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20160712

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20160909

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20161004

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20161104

R151 Written notification of patent or utility model registration

Ref document number: 6039371

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151

LAPS Cancellation because of no payment of annual fees