知財求人 - 知財ポータルサイト「IP Force」
▶ 独立行政法人情報通信研究機構の特許一覧 ▶ 大学共同利用機関法人自然科学研究機構の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024030022
(43)【公開日】2024-03-07
(54)【発明の名称】顕微鏡画像補正方法及び顕微鏡画像補正プログラム
(51)【国際特許分類】
G02B 21/00 20060101AFI20240229BHJP
G01N 21/64 20060101ALI20240229BHJP
【ＦＩ】
G02B21/00
G01N21/64 E
【審査請求】未請求
【請求項の数】4
【出願形態】ＯＬ
(21)【出願番号】P 2022132546
(22)【出願日】2022-08-23
【国等の委託研究の成果に係る記載事項】（出願人による申告）令和３年度、国立研究開発法人科学技術振興機構、戦略的創造研究推進事業「時空間を一括取得する超高速超解像光センサー 超高速超解像光センサーの顕微イメージングへの応用研究」委託研究、産業技術力強化法第１７条の適用を受ける特許出願
(71)【出願人】
【識別番号】301022471
【氏名又は名称】国立研究開発法人情報通信研究機構
(71)【出願人】
【識別番号】504261077
【氏名又は名称】大学共同利用機関法人自然科学研究機構
(74)【代理人】
【識別番号】110001195
【氏名又は名称】弁理士法人深見特許事務所
(72)【発明者】
【氏名】松田 厚志
(72)【発明者】
【氏名】根本 知己
(72)【発明者】
【氏名】大友 康平
(72)【発明者】
【氏名】高橋 泰伽
【テーマコード（参考）】
2G043
2H052
【Ｆターム（参考）】
2G043AA03
2G043EA01
2G043FA01
2G043HA01
2G043HA02
2G043HA09
2G043JA02
2G043LA03
2G043NA01
2H052AA07
2H052AA09
2H052AC04
2H052AC15
2H052AC34
2H052AF14
2H052AF21
2H052AF25
(57)【要約】
【課題】顕微鏡で得られた原画像を補正して、より高いＳＮ比を有する補正画像を算出することができる顕微鏡画像補正方法を提供する。
【解決手段】顕微鏡画像補正方法は、試料の原画像から試料の累乗根画像を算出するステップと、試料の累乗根画像から試料の累乗根画像の空間周波数スペクトルを算出するステップと、試料の累乗根画像の空間周波数スペクトルに高周波除去フィルタを掛けて、試料の累乗根画像の補正空間周波数スペクトルを算出するステップと、試料の累乗根画像の補正空間周波数スペクトルから試料の補正累乗根画像を算出するステップと、試料の補正累乗根画像から試料の補正画像を算出するステップとを備える。
【選択図】図４
【解決手段】顕微鏡画像補正方法は、試料の原画像から試料の累乗根画像を算出するステップと、試料の累乗根画像から試料の累乗根画像の空間周波数スペクトルを算出するステップと、試料の累乗根画像の空間周波数スペクトルに高周波除去フィルタを掛けて、試料の累乗根画像の補正空間周波数スペクトルを算出するステップと、試料の累乗根画像の補正空間周波数スペクトルから試料の補正累乗根画像を算出するステップと、試料の補正累乗根画像から試料の補正画像を算出するステップとを備える。
【選択図】図４
【特許請求の範囲】
【請求項1】
試料の原画像から前記試料の累乗根画像を算出するステップと、
前記試料の前記累乗根画像から前記試料の前記累乗根画像の空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記空間周波数スペクトルに高周波除去フィルタを掛けて、前記試料の前記累乗根画像の補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記補正空間周波数スペクトルから前記試料の補正累乗根画像を算出するステップと、
前記試料の前記補正累乗根画像から前記試料の補正画像を算出するステップとを備える、顕微鏡画像補正方法。
前記試料の前記累乗根画像から前記試料の前記累乗根画像の空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記空間周波数スペクトルに高周波除去フィルタを掛けて、前記試料の前記累乗根画像の補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記補正空間周波数スペクトルから前記試料の補正累乗根画像を算出するステップと、
前記試料の前記補正累乗根画像から前記試料の補正画像を算出するステップとを備える、顕微鏡画像補正方法。
【請求項2】
試料の原画像から前記試料の累乗根画像を算出するステップと、
前記試料の前記累乗根画像から前記試料の前記累乗根画像の第１空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１空間周波数スペクトルに高周波除去フィルタを掛けて、前記試料の前記累乗根画像の第１補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１補正空間周波数スペクトルを、前記試料の前記原画像を得る光学系の点像分布関数の累乗根をフーリエ変換することによって算出される前記光学系の光学伝達関数でデコンボリューションすることによって、前記試料の前記累乗根画像の第２補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第２補正空間周波数スペクトルから前記試料の補正累乗根画像を算出するステップと、
前記試料の前記補正累乗根画像を累乗することによって前記試料の暫定補正画像を算出するステップと、
前記試料の前記暫定補正画像から前記試料の前記暫定補正画像の第２空間周波数スペクトルを算出するステップと、
前記試料の前記暫定補正画像の前記第２空間周波数スペクトルを前記光学系の補正光学伝達関数でデコンボリューションすることによって、前記試料の前記原画像の第３補正空間周波数スペクトルを算出するステップとを備え、前記補正光学伝達関数は、第１補正点像分布関数をフーリエ変換することによって算出され、前記第１補正点像分布関数は、第２補正点像分布関数を累乗することによって算出され、前記第２補正点像分布関数は、前記試料の前記原画像を得る前記光学系の前記点像分布関数の累乗根のフーリエ変換と前記高周波除去フィルタとの積を、前記光学伝達関数でデコンボリューションし、さらに逆フーリエ変換することによって算出され、
前記試料の前記原画像の前記第３補正空間周波数スペクトルから前記試料の補正画像を算出するステップとを備える、顕微鏡画像補正方法。
前記試料の前記累乗根画像から前記試料の前記累乗根画像の第１空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１空間周波数スペクトルに高周波除去フィルタを掛けて、前記試料の前記累乗根画像の第１補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１補正空間周波数スペクトルを、前記試料の前記原画像を得る光学系の点像分布関数の累乗根をフーリエ変換することによって算出される前記光学系の光学伝達関数でデコンボリューションすることによって、前記試料の前記累乗根画像の第２補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第２補正空間周波数スペクトルから前記試料の補正累乗根画像を算出するステップと、
前記試料の前記補正累乗根画像を累乗することによって前記試料の暫定補正画像を算出するステップと、
前記試料の前記暫定補正画像から前記試料の前記暫定補正画像の第２空間周波数スペクトルを算出するステップと、
前記試料の前記暫定補正画像の前記第２空間周波数スペクトルを前記光学系の補正光学伝達関数でデコンボリューションすることによって、前記試料の前記原画像の第３補正空間周波数スペクトルを算出するステップとを備え、前記補正光学伝達関数は、第１補正点像分布関数をフーリエ変換することによって算出され、前記第１補正点像分布関数は、第２補正点像分布関数を累乗することによって算出され、前記第２補正点像分布関数は、前記試料の前記原画像を得る前記光学系の前記点像分布関数の累乗根のフーリエ変換と前記高周波除去フィルタとの積を、前記光学伝達関数でデコンボリューションし、さらに逆フーリエ変換することによって算出され、
前記試料の前記原画像の前記第３補正空間周波数スペクトルから前記試料の補正画像を算出するステップとを備える、顕微鏡画像補正方法。
【請求項3】
試料の原画像から前記試料の累乗根画像を算出するステップと、
前記試料の前記累乗根画像から前記試料の前記累乗根画像の第１空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１空間周波数スペクトルを、前記試料の前記原画像を得る光学系の点像分布関数の累乗根をフーリエ変換することによって算出される前記光学系の光学伝達関数でデコンボリューションすることによって、前記試料の前記累乗根画像の第１補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１補正空間周波数スペクトルから前記試料の補正累乗根画像を算出するステップと、
前記試料の前記補正累乗根画像を累乗することによって前記試料の暫定補正画像を算出するステップと、
前記試料の前記暫定補正画像から前記試料の前記暫定補正画像の第２空間周波数スペクトルを算出するステップと、
前記試料の前記暫定補正画像の前記第２空間周波数スペクトルを前記光学系の補正光学伝達関数でデコンボリューションすることによって、前記試料の前記原画像の第２補正空間周波数スペクトルを算出するステップとを備え、前記補正光学伝達関数は、第１補正点像分布関数をフーリエ変換することによって算出され、前記第１補正点像分布関数は、第２補正点像分布関数を累乗することによって算出され、前記第２補正点像分布関数は、前記試料の前記原画像を得る前記光学系の前記点像分布関数の累乗根のフーリエ変換と高周波除去フィルタとの積を、前記光学伝達関数でデコンボリューションし、さらに逆フーリエ変換することによって算出され、
前記試料の前記原画像の前記第２補正空間周波数スペクトルから前記試料の補正画像を算出するステップとを備える、顕微鏡画像補正方法。
前記試料の前記累乗根画像から前記試料の前記累乗根画像の第１空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１空間周波数スペクトルを、前記試料の前記原画像を得る光学系の点像分布関数の累乗根をフーリエ変換することによって算出される前記光学系の光学伝達関数でデコンボリューションすることによって、前記試料の前記累乗根画像の第１補正空間周波数スペクトルを算出するステップと、
前記試料の前記累乗根画像の前記第１補正空間周波数スペクトルから前記試料の補正累乗根画像を算出するステップと、
前記試料の前記補正累乗根画像を累乗することによって前記試料の暫定補正画像を算出するステップと、
前記試料の前記暫定補正画像から前記試料の前記暫定補正画像の第２空間周波数スペクトルを算出するステップと、
前記試料の前記暫定補正画像の前記第２空間周波数スペクトルを前記光学系の補正光学伝達関数でデコンボリューションすることによって、前記試料の前記原画像の第２補正空間周波数スペクトルを算出するステップとを備え、前記補正光学伝達関数は、第１補正点像分布関数をフーリエ変換することによって算出され、前記第１補正点像分布関数は、第２補正点像分布関数を累乗することによって算出され、前記第２補正点像分布関数は、前記試料の前記原画像を得る前記光学系の前記点像分布関数の累乗根のフーリエ変換と高周波除去フィルタとの積を、前記光学伝達関数でデコンボリューションし、さらに逆フーリエ変換することによって算出され、
前記試料の前記原画像の前記第２補正空間周波数スペクトルから前記試料の補正画像を算出するステップとを備える、顕微鏡画像補正方法。
【請求項4】
請求項１から請求項３のいずれか一項に記載の前記顕微鏡画像補正方法の各ステップをプロセッサに実行させる、顕微鏡画像補正プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、顕微鏡画像補正方法及び顕微鏡画像補正プログラムに関する。
【背景技術】
【0002】
特表２０１８－５３３７６９号公報（特許文献１）は、二光子励起顕微鏡法（ＴＰＭ）を開示している。
【先行技術文献】
【特許文献】
【0003】
【特許文献1】特表２０１８－５３３７６９号公報
【発明の概要】
【発明が解決しようとする課題】
【0004】
本開示の目的は、顕微鏡で得られた原画像を補正して、より高いＳＮ比を有する補正画像を算出することができる顕微鏡画像補正方法及び顕微鏡画像補正プログラムを提供することである。
【課題を解決するための手段】
【0005】
本開示の第一の局面の顕微鏡画像補正方法は、試料の原画像から試料の累乗根画像を算出するステップと、試料の累乗根画像から試料の累乗根画像の空間周波数スペクトルを算出するステップと、試料の累乗根画像の空間周波数スペクトルに高周波除去フィルタを掛けて、試料の累乗根画像の補正空間周波数スペクトルを算出するステップと、試料の累乗根画像の補正空間周波数スペクトルから試料の補正累乗根画像を算出するステップと、試料の補正累乗根画像から試料の補正画像を算出するステップとを備える。
【0006】
本開示の第二の局面の顕微鏡画像補正方法は、試料の原画像から試料の累乗根画像を算出するステップと、試料の累乗根画像から試料の累乗根画像の第１空間周波数スペクトルを算出するステップと、試料の累乗根画像の第１空間周波数スペクトルに高周波除去フィルタを掛けて、試料の累乗根画像の第１補正空間周波数スペクトルを算出するステップと、試料の累乗根画像の第１補正空間周波数スペクトルを、試料の原画像を得る光学系の点像分布関数の累乗根をフーリエ変換することによって算出される光学系の光学伝達関数でデコンボリューションすることによって、試料の累乗根画像の第２補正空間周波数スペクトルを算出するステップと、試料の累乗根画像の第２補正空間周波数スペクトルから試料の補正累乗根画像を算出するステップと、試料の補正累乗根画像を累乗することによって、試料の暫定補正画像を算出するステップと、試料の暫定補正画像から試料の暫定補正画像の第２空間周波数スペクトルを算出するステップと、試料の暫定補正画像の第２空間周波数スペクトルを光学系の補正光学伝達関数でデコンボリューションすることによって、試料の原画像の第３補正空間周波数スペクトルを算出するステップを備える。光学系の補正光学伝達関数は、第１補正点像分布関数をフーリエ変換することによって算出される。第１補正点像分布関数は、第２補正点像分布関数を累乗することによって算出される。第２補正点像分布関数は、試料の原画像を得る光学系の点像分布関数の累乗根のフーリエ変換と高周波除去フィルタとの積を、光学系の光学伝達関数でデコンボリューションし、さらに逆フーリエ変換することによって算出される。本開示の第二の局面の顕微鏡画像補正方法は、試料の原画像の第３補正空間周波数スペクトルから試料の補正画像を算出するステップを備える。
【0007】
本開示の第三の局面の顕微鏡画像補正方法は、試料の原画像から試料の累乗根画像を算出するステップと、試料の累乗根画像から試料の累乗根画像の第１空間周波数スペクトルを算出するステップと、試料の累乗根画像の第１空間周波数スペクトルを、試料の原画像を得る光学系の点像分布関数の累乗根をフーリエ変換することによって算出される光学系の光学伝達関数でデコンボリューションすることによって、試料の累乗根画像の第１補正空間周波数スペクトルを算出するステップと、試料の累乗根画像の第１補正空間周波数スペクトルから試料の補正累乗根画像を算出するステップと、試料の補正累乗根画像を累乗することによって、試料の暫定補正画像を算出するステップと、試料の暫定補正画像から試料の暫定補正画像の第２空間周波数スペクトルを算出するステップと、試料の暫定補正画像の第２空間周波数スペクトルを光学系の補正光学伝達関数でデコンボリューションすることによって、試料の原画像の第２補正空間周波数スペクトルを算出するステップとを備える。補正光学伝達関数は、第１補正点像分布関数をフーリエ変換することによって算出される。第１補正点像分布関数は、第２補正点像分布関数を累乗することによって算出される。第２補正点像分布関数は、試料の原画像を得る光学系の点像分布関数の累乗根のフーリエ変換と高周波除去フィルタとの積を、光学系の光学伝達関数でデコンボリューションし、さらに逆フーリエ変換することによって算出される。本開示の第三の局面の顕微鏡画像補正方法は、試料の原画像の第２補正空間周波数スペクトルから試料の補正画像を算出するステップを備える。
【0008】
本開示の顕微鏡画像補正プログラムは、本開示の第一局面から第三局面のいずれかの顕微鏡画像補正方法の各ステップをプロセッサに実行させる。
【発明の効果】
【0009】
本開示の顕微鏡画像補正方法及び本開示の顕微鏡画像補正プログラムによれば、顕微鏡で得られた原画像を補正して、より高いＳＮ比を有する補正画像を算出することができる。
【図面の簡単な説明】
【0010】
【図1】実施の形態１から実施の形態３の顕微鏡システムを示す概略図である。
【図2】実施の形態１から実施の形態３の顕微鏡システムに含まれるコンピュータの機能的構成を説明するブロック図である。
【図3】実施の形態１の画像処理部の機能的構成を説明するブロック図である。
【図4】実施の形態１の顕微鏡で得られた原画像の補正方法のフローチャートを示す図である。
【図5】実施の形態１及び実施の形態２のデノイズ処理のフローチャートを示す図である。
【図6】実施の形態１の補正画像を算出するステップのフローチャートを示す図である。
【図7】二光子励起顕微鏡もしくは三光子励起顕微鏡のような多光子励起顕微鏡、共焦点顕微鏡または超解像顕微鏡の原画像及びその累乗根画像、光検出器が単一素子でありピンホールを持たない一光子励起走査型顕微鏡の一光子励起画像、並びに、それらのフーリエ変換（空間周波数スペクトル）の関係を示す図である。
【図8】実施の形態２の画像処理部の機能的構成を説明するブロック図である。
【図9】実施の形態２の顕微鏡で得られた原画像の補正方法のフローチャートを示す図である。
【図10】実施の形態２の第２デコンボリューション処理ステップのフローチャートを示す図である。
【図11】実施の形態３の画像処理部の機能的構成を説明するブロック図である。
【図12】実施の形態３の顕微鏡で得られた原画像の補正方法のフローチャートを示す図である。
【図13】実施の形態３の第１デコンボリューション処理ステップのフローチャートを示す図である。
【図14】実施の形態３の第２デコンボリューション処理ステップのフローチャートを示す図である。
【図15】第１原画像を示す図である。
【図16】実施例１－１の補正画像を示す図である。
【図17】実施例２－１の補正画像を示す図である。
【図18】実施例３－１の補正画像を示す図である。
【図19】実施例１－１から実施例３－２の画像のＳＮ比の改善率を表すグラフを示す図である。
【発明を実施するための形態】
【0011】
以下、実施の形態を説明する。なお、同一の構成には同一の参照番号を付し、その説明は繰り返さない。
【0012】
（実施の形態１）
【0013】
図１から図３を参照して、実施の形態１の顕微鏡システム１を説明する。図１に示されるように、顕微鏡システム１は、顕微鏡２と、コンピュータ３と、表示装置４とを備える。顕微鏡２は、試料１２の原画像を得ることができる顕微鏡である。本実施の形態では、顕微鏡２は、二光子励起顕微鏡である。
【0014】
＜顕微鏡の構成＞
【0015】
図１に示されるように、顕微鏡２は、照明光学系５と、試料ステージ１３と、観察光学系１５とを含む。
【0016】
照明光学系５は、例えば、光源６と、走査ミラー８と、レンズ９ａ，９ｂ，９ｃ，９ｄと、ダイクロイックミラー１０と、対物レンズ１１とを含む。
【0017】
光源６は、光７を出力する。光源６は、例えば、チタンサファイアレーザのようなフェムト秒超短パルスレーザである。光７は、例えば、フェムト秒パルスレーザビームである。光７は、レンズ９ａを通って、レンズ９ｂに入射する。レンズ９ａ，９ｂによって、光７は拡げられる。光７は、レンズ９ｂでコリメートされて、走査ミラー８に入射する。走査ミラー８は、光７を反射するとともに、光７を走査する。走査ミラー８は、例えば、ガルバノミラーである。光７は、レンズ９ｃを通って、レンズ９ｄに入射する。光７は、レンズ９ｄでコリメートされて、ダイクロイックミラー１０に入射する。ダイクロイックミラー１０は、光７を試料１２に向けて反射する。光７は、対物レンズ１１で集光されて、試料１２の一部分を照射する。光７の集光点は、走査ミラー８によって、試料１２において、試料１２における光７の入射方向に垂直な面内に走査される。
【0018】
試料１２は、試料ステージ１３によって支持されている。試料１２は、特に限定されないが、例えば、生細胞、分子、タンパク質またはＤＮＡである。試料１２は、例えば、蛍光色素分子のような蛍光物質によって標識されている。光７が試料１２に照射されると、光７は試料１２を標識する蛍光物質に吸収される。蛍光物質は、光７を二光子吸収して励起され、光１４を放射する。光１４は、例えば、蛍光光のようなインコヒーレント光である。こうして、試料１２から光１４が放射される。
【0019】
観察光学系１５は、例えば、対物レンズ１１と、ダイクロイックミラー１０と、光検出器１６とを含む。試料１２から放射された光１４は、対物レンズ１１によってコリメートされて、ダイクロイックミラー１０に入射する。光１４は、ダイクロイックミラー１０を通って、光検出器１６に入射する。光検出器１６は、試料１２の各部分から放射される光１４の強度を検出する。光検出器１６は、例えば、光電子増倍管（ＰＭＴ）である。
【0020】
図１を参照して、コンピュータ３は、顕微鏡２に通信可能に接続されている。コンピュータ３は、顕微鏡２の制御、試料１２の原画像を構築するための画像構築処理、及び、試料１２の原画像を補正するための画像処理などを行う。コンピュータ３は、例えば、パーソナルコンピュータ、マイクロコンピュータ、クラウドサーバまたはスマートデバイス（スマートフォンまたはタブレット端末など）である。表示装置４は、コンピュータ３に通信可能に接続されている。表示装置４は、例えば、液晶表示装置または有機ＥＬ表示装置である。表示装置４は、例えば、顕微鏡２の動作状態、試料１２の原画像、及び、試料１２の補正画像などを表示する。
【0021】
＜コンピュータ３の構成＞
【0022】
図１を参照して、コンピュータ３は、プロセッサ２０と、記憶装置２２とを含む。記憶装置２２は、例えば、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、または、ハードディスクである。記憶装置２２は、プロセッサ２０により実行されるコンピュータプログラム、試料１２の原画像、及び、試料１２の補正画像などを格納する。コンピュータプログラムは、例えば、試料１２の原画像を取得するための画像取得プログラムと、試料１２の原画像を補正して試料１２の補正画像を算出するための画像補正プログラムとを含む。画像取得プログラム及び画像補正プログラムは、非一時的なコンピュータ可読記憶媒体に格納されてもよい。
【0023】
プロセッサ２０は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、または、ＦＰＧＡ（field-programmable gate array）などによって構成されている。プロセッサ２０は、記憶装置２２に格納されているコンピュータプログラムを実行する。プロセッサ２０がコンピュータプログラムを実行することにより、コンピュータ３は、図２及び図３に示される機能を実現する。図２及び図３を参照して、コンピュータ３の機能構成を説明する。
【0024】
図２に示されるように、コンピュータ３は、原画像取得部２５と、画像処理部３０とを含む。原画像取得部２５は、プロセッサ２０が画像取得プログラムを実行することによって、実現される。画像処理部３０は、顕微鏡２で得られた試料１２の原画像を補正して、試料１２の補正画像を算出する。画像処理部３０は、プロセッサ２０が画像補正プログラムを実行することによって、実現される。
【0025】
図３に示されるように、画像処理部３０は、累乗根画像算出部３１と、デノイズ処理部３２と、補正画像算出部３５とを含む。累乗根画像算出部３１は、試料１２の原画像から累乗根画像を算出するステップ（ステップＳ１、図４を参照）を行う。デノイズ処理部３２は、デノイズ処理ステップ（ステップＳ２、図４及び図５を参照）を行う。補正画像算出部３５は、試料１２の補正画像を算出するステップ（ステップＳ３、図４及び図６を参照）を行う。
【0026】
＜試料１２の原画像を取得する方法＞
【0027】
試料１２の原画像を取得する方法は、原画像取得部２５によって実行される。具体的には、原画像取得部２５は、光源６を制御して、光源６に光７を出射させる。光７は、レンズ９ａを通って、レンズ９ｂに入射する。レンズ９ａ，９ｂによって、光７は拡げられる。光７は、レンズ９ｂでコリメートされて、走査ミラー８に入射する。光７は、走査ミラー８で反射される。光７は、レンズ９ｃ，９ｄを通って、ダイクロイックミラー１０で反射される。光７は、対物レンズ１１で集光されて、試料１２の一部分を照射する。原画像取得部２５は、走査ミラー８を制御して、試料１２において、光７の集光点を、試料１２における光７の入射方向に垂直な面内に走査する。光検出器１６は、試料１２の各部分から放射される光１４の強度を検出する。原画像取得部２５は、試料１２の各部分から放射されて光検出器１６で検出された光１４の強度から、試料１２の原画像として、試料１２の二次元原画像を構築する。例えば、顕微鏡２が二光子励起顕微鏡である場合、試料１２の原画像は、二光子励起画像Ｉ２である。また、顕微鏡２が三光子励起顕微鏡である場合、試料１２の原画像は、三光子励起画像Ｉ３である。
【0028】
原画像取得部２５は、さらに対物レンズ１１を制御して、試料１２において、光７の集光点を、試料１２における光７の入射方向に移動させてもよい。原画像取得部２５は、試料１２の二次元原画像を積層して、試料１２の原画像として、試料１２の三次元原画像を構築してもよい。
【0029】
＜顕微鏡２で得られた試料１２の原画像の補正方法＞
【0030】
図４から図７を参照して、顕微鏡２で得られた試料１２の原画像から、試料１２の補正画像を算出する方法を説明する。試料１２の補正画像を算出する方法は、画像処理部３０によって実行される。以下では、主に、顕微鏡２が二光子励起顕微鏡であり、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合について説明する。
【0031】
試料１２の原画像（例えば、二光子励起画像Ｉ２）を得るために試料１２に照射される光７の波長と同一波長を有する光を用いて試料１２が一光子励起されることによって得られる試料１２の画像を、一光子励起画像Ｉ１（図７を参照）とする。原画像の累乗根によって与えられかつ一光子励起画像Ｉ１と解像度が等価である画像を、累乗根画像Ｉ１’（図７を参照）と定義する。試料の原画像Ｉｍがｍ光子励起画像である時、累乗根はｍ乗根を意味し、累乗はｍ乗を意味し、かつ、累乗根画像Ｉ１’は原画像Ｉｍを（１／ｍ）乗した画像を意味する（なお、ｍは自然数）。
【0032】
図４に示されるように、試料１２の原画像から累乗根画像Ｉ１’を算出する（ステップＳ１）。試料１２の原画像の累乗根をとることによって、累乗根画像Ｉ１’が算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、累乗根画像Ｉ１’は、式（１）によって算出される。すなわち、二光子励起画像Ｉ２の各画素の光強度の平方根を算出することによって、累乗根画像Ｉ１’の各画素の光強度が算出される。
【数1】
【0033】
例えば、試料１２の原画像が試料１２の三光子励起画像Ｉ３である場合、累乗根画像Ｉ１’は、式（２）によって算出される。すなわち、三光子励起画像Ｉ３の各画素の光強度の三乗根を算出することによって、累乗根画像Ｉ１’の各画素の光強度が算出される。
【数2】
【0034】
【0035】
具体的には、式（３）に示されるように、累乗根画像Ｉ１’をフーリエ変換して、累乗根画像Ｉ１’の空間周波数スペクトルＦ１’を算出する（ステップＳ２－１）。本明細書において、ＦＴは、フーリエ変換を表す。
【数3】
【0036】
式（４）に示されるように、累乗根画像Ｉ１’の空間周波数スペクトルＦ１’に高周波除去フィルタＧ（ｋx、ｋy、ｋz）（以下、「高周波除去フィルタＧ」ということがある。）を掛けることによって、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ａ’を算出する（ステップＳ２－２）。そのため、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ａ’では、高周波ノイズが除去される。
【数4】
【0037】
高周波除去フィルタＧは、例えば、水平方向の波数ｋrと垂直方向（ｚ方向）の波数ｋｚの二軸において、式（５）によって与えられる。垂直方向は、試料１２における光７の入射方向を意味する。水平方向は、試料１２における光７の入射方向に垂直な方向を意味する。
【数5】
【0038】
ｋｒは、式（６）によって与えられる。ｘ方向は、試料１２における光７の入射方向に垂直な方向のうちの一方向を意味する。ｙ方向は、試料１２における光７の入射方向に垂直な方向のうち、ｘ方向に垂直な方向を意味する。ｋｘはｘ方向の波数を意味し、ｋｙはｙ方向の波数を意味する。
【数6】
【0039】
Ｋは、式（７）によって与えられる。
【数7】
【0040】
式（５）における角度αは、対物レンズ１１を通過して試料１２に入射する光の最大入射角である（図１を参照）。角度αは、対物レンズ１１の開口数ＮＡと、対物レンズ１１から試料１２までの媒質の屈折率ｎとによって規定される。具体的には、角度αは、式（８）によって与えられる。
【数8】
【0041】
【0042】
具体的には、式（９）に示されるように、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ａ’を逆フーリエ変換して、累乗根補正画像Ｉ１ａ’を算出する（ステップＳ３－１）。本明細書において、ＦＴ－１は、逆フーリエ変換を表す。
【数9】
【0043】
累乗根補正画像Ｉ１ａ’から試料１２の補正画像Ｉ２ｆを算出する（ステップＳ３－２）。具体的には、累乗根補正画像Ｉ１ａ’を累乗することによって、試料１２の補正画像Ｉ２ｆが算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、式（１０）に示されるように、累乗根補正画像Ｉ１ａ’を二乗することによって、試料１２の補正画像Ｉ２ｆを算出する。累乗根補正画像Ｉ１ａ’の各画素の光強度の二乗を算出することによって、試料１２の補正画像Ｉ２ｆの各画素の光強度が算出される。高周波除去フィルタＧによって、試料１２の補正画像Ｉ２ｆでは、高周波ノイズが除去される。そのため、補正画像Ｉ２ｆのＳＮ比（コントラスト比）は向上する。
【数10】
【0044】
試料１２の原画像が試料１２の三光子励起画像Ｉ３である場合、式（１１）に示されるように、累乗根補正画像Ｉ１ａ’を三乗することによって、試料１２の補正画像Ｉ３ｆを算出する。累乗根補正画像Ｉ１ａ’の各画素の光強度の三乗を算出することによって、試料１２の補正画像Ｉ３ｆの各画素の光強度が算出される。高周波除去フィルタＧによって、試料１２の補正画像Ｉ３ｆでは、高周波ノイズが除去される。そのため、補正画像Ｉ３ｆのＳＮ比（コントラスト比）は向上する。
【数11】
【0045】
本実施の形態の効果を説明する。
【0046】
本実施の形態の顕微鏡画像補正方法は、試料１２の原画像（例えば、二光子励起画像Ｉ２）から試料１２の累乗根画像Ｉ１’を算出するステップ（ステップＳ１）と、試料１２の累乗根画像Ｉ１’から試料１２の累乗根画像Ｉ１’の空間周波数スペクトルＦ１’を算出するステップ（ステップＳ２－１）と、試料１２の累乗根画像Ｉ１’の空間周波数スペクトルＦ１’に高周波除去フィルタＧを掛けて、試料１２の累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ａ’を算出するステップ（ステップＳ２－２）と、試料１２の累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ａ’から試料１２の累乗根補正画像Ｉ１ａ’を算出するステップ（ステップＳ３－１）と、試料１２の累乗根補正画像Ｉ１ａ’から試料１２の補正画像Ｉ２ｆを算出するステップ（ステップＳ３－２）とを備える。
【0047】
本実施の形態の顕微鏡画像補正プログラムは、本実施の形態の顕微鏡画像補正方法の各ステップをプロセッサ２０に実行させる。
【0048】
図７に示されるように、照明光学系５を用いて試料１２を一光子励起すると、試料１２のうち光７が通る部分の全体において一光子励起が生じる（図７の画像Ｉ１の実空間の点像を参照）。そのため、試料１２の全体から光が放射される。光検出器が単一素子でありかつピンホールを持たない一光子励起走査型顕微鏡では、試料１２の全体から放射される光を合算して検出する。これに対し、二光子励起顕微鏡もしくは三光子励起顕微鏡のような多光子励起顕微鏡、共焦点顕微鏡または超解像顕微鏡（Image Scanning Microscopy（ＩＳＭ））などによって検出される光１４は、試料１２のうち対物レンズ１１の焦点付近の部分のみに限られている（図７の原画像Ｉｍの実空間の点像を参照）。その理由は、これらの顕微鏡の点像分布関数が、ピンホールを持たない単純な一光子励起走査型顕微鏡の点像分布関数を累乗（ｍ乗）した関係にあるため、照明光の強度が強い部分では光検出器１６によって検出される光の量がさらに増加するのに対し、照明光の強度が弱い部分では光検出器１６によって検出される光の量がさらに減少するためである。
【0049】
その結果、原画像Ｉｍの解像度は、画像Ｉ１の解像度よりも大きくなる。言い換えると、空間周波数領域では、原画像Ｉｍの遮断周波数は、画像Ｉ１の遮断周波数よりも大きくなる（図７の原画像Ｉｍの空間周波数及び画像Ｉ１の空間周波数を参照）。これらの画像の遮断周波数は、照明光学系５によって与えられる。原画像Ｉｍの累乗根により得られる累乗根画像Ｉ１’の解像度は、画像Ｉ１の解像度と同等である（図７の実空間の点像の欄を参照）。したがって、画像Ｉ１と同様に、累乗根画像Ｉ１’の遮断周波数は、原画像Ｉｍの遮断周波数よりも小さくなり、画像Ｉ１の遮断周波数と同等になる（図７の空間周波数の欄を参照）。
【0050】
累乗根画像Ｉ１’における試料１２からの信号成分は、試料１２の原画像Ｉｍにおける試料１２からの信号成分よりも、より小さな空間周波数領域に分布している。これに対し、累乗根画像Ｉ１’におけるノイズ成分は、試料１２の原画像Ｉｍにおけるノイズ成分と同様に、全空間周波数領域にわたってランダムに分布している。そのため、累乗根画像Ｉ１’の空間周波数スペクトルＦ１’に高周波除去フィルタＧを掛けるというデノイズ処理ステップ（Ｓ２）を行うことによって、より広い空間周波数領域に分布するノイズを除去することができる。こうして、より高いＳＮ比（コントラスト比）を有する補正画像を算出することができる。
【0051】
（実施の形態２）
【0052】
図１、図２及び図８を参照して、実施の形態２の顕微鏡システム１を説明する。実施の形態２の顕微鏡システム１は、実施の形態１の顕微鏡システム１と同様に構成されているが、コンピュータ３の機能構成が異なっている。
【0053】
図８に示されるように、本実施の形態の画像処理部３０は、実施の形態１の画像処理部３０と同様に構成されているが、第１デコンボリューション処理部３３と、第２デコンボリューション処理部３４とをさらに含む点で、実施の形態１の画像処理部３０と異なっている。
【0054】
累乗根画像算出部３１は、試料１２の原画像（例えば、二光子励起画像Ｉ２）から累乗根画像Ｉ１’を算出するステップ（ステップＳ１、図９を参照）を行う。デノイズ処理部３２は、デノイズ処理ステップ（ステップＳ２、図９を参照）を行う。第１デコンボリューション処理部３３は、第１デコンボリューション処理ステップ（ステップＳ３、図９を参照）を行う。第２デコンボリューション処理部３４は、第２デコンボリューション処理ステップ（ステップＳ４、図９及び図１０を参照）を行う。補正画像算出部３５は、試料１２の補正画像を算出するステップ（ステップＳ５、図９を参照）を行う。
【0055】
図５、図９及び図１０を参照して、実施の形態２の試料１２の原画像の補正方法を説明する。試料１２の補正画像を算出する方法は、画像処理部３０によって実行される。以下では、主に、顕微鏡２が二光子励起顕微鏡であり、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合について説明する。実施の形態２の試料１２の原画像の補正方法は、実施の形態１の試料１２の原画像の補正方法と同様であるが、主に、第１デコンボリューション処理ステップ（ステップＳ３）と、第２デコンボリューション処理ステップ（ステップＳ４）とをさらに備えている点において、実施の形態１の試料１２の原画像の補正方法と異なっている。
【0056】
試料１２に照射される光７は、振幅Ｅ及び位相θを有する電界Ｅｅｘｐ（ｉθ）として表現される。光７の強度｜Ｅ２｜として光７が試料１２に吸収されて、試料１２が光７によって励起される。試料１２による光７の強度の吸収において、光７の電界の位相θの情報が失われ、光７の電界の振幅Ｅが二乗される。そのため、試料１２の原画像の空間周波数スペクトルでは、照明光学系５の遮断周波数に近い空間周波数成分の振幅成分が減少して、試料１２の原画像のＳＮ比（コントラスト比）が低下する。試料１２の原画像が試料１２の一光子励起画像Ｉ１である場合、照明光学系５の遮断周波数に近い空間周波数成分の振幅成分の減少は一回発生するが、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、照明光学系５の遮断周波数に近い空間周波数成分の振幅成分の減少は二回発生する。そこで、試料１２の二光子励起画像Ｉ２において、照明光学系５の遮断周波数に近い空間周波数成分の振幅成分を回復させるために、第１デコンボリューション処理ステップ（ステップＳ３）と、第２デコンボリューション処理ステップ（ステップＳ４）とを行う。
【0057】
【0058】
図９に示されるように、デノイズ処理（ステップＳ２）の後に、第１デコンボリューション処理ステップ（ステップＳ３）を行う。
【0059】
具体的には、式（１２）に示されるように、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ａ’を、光学系（照明光学系５）の光学伝達関数ＯＴＦ１でデコンボリューションして、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ｂ’を算出する。そのため、照明光学系５の遮断周波数未満の周波数を有するノイズが低減される。第１デコンボリューション処理は、例えば、ＯＴＦ１をウィナーフィルタとして用いるデコンボリューション処理を含む。式（１２）におけるｗ１は、ウィナー定数を表す。
【数12】
【0060】
光学伝達関数ＯＴＦ１は、試料１２の原画像（例えば、二光子励起画像Ｉ２）を得る光学系（照明光学系５）の点像分布関数ＰＳＦ２の累乗根をフーリエ変換することによって算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、ＯＴＦ１は、式（１３）によって与えられる。
【数13】
【0061】
本実施の形態の目的は、顕微鏡２で得られた試料１２の原画像のＳＮ比（コントラスト比）を向上させることである。そのため、第１デコンボリューション処理は、光学伝達関数ＯＴＦ１で行われてもよいし、光学伝達関数ＯＴＦ１の絶対値である変調伝達関数ＭＴＦ１で行われてもよい。
【0062】
図１０に示されるように、第１デコンボリューション処理ステップ（ステップＳ３）の後に、第２デコンボリューション処理ステップ（ステップＳ４）を行う。
【0063】
具体的には、式（１４）に示されるように、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ｂ’を逆フーリエ変換して、累乗根補正画像Ｉ１ｂ’を算出する（ステップＳ４－１）。
【数14】
【0064】
累乗根補正画像Ｉ１b’を累乗することによって、試料１２の暫定補正画像Ｉ２ｂを算出する（ステップＳ４－２）。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、式（１５）に示されるように、累乗根補正画像Ｉ１b’を二乗することによって、試料１２の暫定補正画像Ｉ２ｂを算出する。具体的には、累乗根補正画像Ｉ１ｂ’の各画素の光強度の二乗を算出することによって、試料１２の暫定補正画像Ｉ２ｂの各画素の光強度が算出される。試料１２の暫定補正画像Ｉ２ｂでは、デノイズ処理（ステップＳ２）によって照明光学系５の遮断周波数以上の周波数を有する高周波ノイズが除去されているとともに、第１デコンボリューション処理ステップ（ステップＳ３）によって照明光学系５の遮断周波数未満の周波数を有するノイズが低減されている。
【数15】
【0065】
式（１６）に示されるように、試料１２の暫定補正画像Ｉ２ｂをフーリエ変換して、試料１２の暫定補正画像Ｉ２ｂの空間周波数スペクトルＦ２ｂを算出する（ステップＳ４－３）。
【数16】
【0066】
式（１７）に示されるように、試料１２の暫定補正画像Ｉ２ｂの空間周波数スペクトルＦ２ｂを光学系（照明光学系５）の補正光学伝達関数ＯＴＦ２ｂでデコンボリューションすることによって、試料１２の原画像の補正空間周波数スペクトルＦ２ｃを算出する（ステップＳ４－４）。第２デコンボリューション処理は、例えば、ＯＴＦ２ｂをウィナーフィルタとして用いるデコンボリューション処理を含む。式（１７）におけるｗ２は、ウィナー定数を表す。
【数17】
【0067】
本実施の形態の目的は、顕微鏡２で得られた試料１２の原画像のＳＮ比（コントラスト比）を向上させることである。そのため、第２デコンボリューション処理は、補正光学伝達関数ＯＴＦ２ｂで行われてもよいし、補正光学伝達関数ＯＴＦ２ｂの絶対値である補正変調伝達関数ＭＴＦ２ｂで行われてもよい。
【0068】
補正光学伝達関数ＯＴＦ２ｂは、補正点像分布関数ＰＳＦ２ｂをフーリエ変換することによって算出される。具体的には、補正光学伝達関数ＯＴＦ２ｂは、式（１８）によって与えられる。
【数18】
【0069】
補正点像分布関数ＰＳＦ２ｂは、補正点像分布関数ＰＳＦ１ｂを累乗することによって算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、ＰＳＦ２ｂは、式（１９）によって与えられる。
【数19】
【0070】
補正点像分布関数ＰＳＦ１ｂは、試料１２の原画像を得る光学系（照明光学系５）の点像分布関数ＰＳＦ２の累乗根のフーリエ変換と高周波除去フィルタＧとの積を、光学伝達関数ＯＴＦ１でデコンボリューションし、さらに逆フーリエ変換することによって算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、ＰＳＦ１ｂは、式（２０）によって与えられる。
【数20】
【0071】
本実施の形態の目的は、顕微鏡２で得られた試料１２の原画像のＳＮ比（コントラスト比）を向上させることである。そのため、補正点像分布関数ＰＳＦ１ｂを算出するために、点像分布関数ＰＳＦ２の累乗根のフーリエ変換と高周波除去フィルタＧとの積は、光学伝達関数ＯＴＦ１でデコンボリューションされてもよいし、光学伝達関数ＯＴＦ１の絶対値である変調伝達関数ＭＴＦ１でデコンボリューションされてもよい。
【0072】
【0073】
本実施の形態の効果は、実施の形態１の効果に加えて、以下の効果を奏する。
【0074】
本実施の形態の顕微鏡画像補正方法は、試料１２の原画像（例えば、二光子励起画像Ｉ２）から試料１２の累乗根画像Ｉ１’を算出するステップ（ステップＳ１）と、試料１２の累乗根画像Ｉ１’から試料１２の累乗根画像Ｉ１’の第１空間周波数スペクトル（空間周波数スペクトルＦ１’）を算出するステップ（ステップＳ２－１）と、試料１２の累乗根画像Ｉ１’の第１空間周波数スペクトルに高周波除去フィルタＧを掛けて、試料１２の累乗根画像Ｉ１’の第１補正空間周波数スペクトル（補正空間周波数スペクトルＦ１ａ’）を算出するステップ（ステップＳ２－２）と、試料１２の累乗根画像Ｉ１’の第１補正空間周波数スペクトルを、試料１２の原画像を得る光学系（照明光学系５）の点像分布関数の累乗根をフーリエ変換することによって算出される光学系の光学伝達関数（光学伝達関数ＯＴＦ１）でデコンボリューションすることによって、試料１２の累乗根画像Ｉ１’の第２補正空間周波数スペクトル（補正空間周波数スペクトルＦ１ｂ’）を算出するステップ（ステップＳ３）と、試料１２の累乗根画像Ｉ１’の第２補正空間周波数スペクトルから試料１２の累乗根補正画像Ｉ１ｂ’を算出するステップ（ステップＳ４－１）と、試料１２の累乗根補正画像Ｉ１ｂ’を累乗することによって試料１２の暫定補正画像Ｉ２ｂを算出するステップ（ステップＳ４－２）と、試料１２の暫定補正画像Ｉ２ｂから試料１２の暫定補正画像Ｉ２ｂの第２空間周波数スペクトル（空間周波数スペクトルＦ２ｂ）を算出するステップ（ステップＳ４－３）と、試料１２の暫定補正画像Ｉ２ｂの第２空間周波数スペクトルを光学系の補正光学伝達関数（補正光学伝達関数ＯＴＦ２ｂ）でデコンボリューションすることによって、試料１２の原画像の第３補正空間周波数スペクトル（補正空間周波数スペクトルＦ２ｃ）を算出するステップ（ステップＳ４－４）とを備える。補正光学伝達関数は、第１補正点像分布関数（補正点像分布関数ＰＳＦ２ｂ）をフーリエ変換することによって算出される。第１補正点像分布関数は、第２補正点像分布関数（補正点像分布関数ＰＳＦ１ｂ）を累乗することによって算出される。第２補正点像分布関数は、試料１２の原画像を得る光学系の点像分布関数の累乗根のフーリエ変換と高周波除去フィルタＧとの積を、光学系の光学伝達関数（光学伝達関数ＯＴＦ１）でデコンボリューションし、さらに逆フーリエ変換することによって算出される。本実施の形態の顕微鏡画像補正方法は、試料１２の原画像の第３補正空間周波数スペクトルから試料１２の補正画像Ｉ２ｆを算出するステップ（ステップＳ５）とを備える。
【0075】
本実施の形態の顕微鏡画像補正プログラムは、本実施の形態の顕微鏡画像補正方法の各ステップをプロセッサ２０に実行させる。
【0076】
照明光学系５の光７が試料１２により吸収される際に光７の電界の振幅Ｅが二乗されることに起因して、試料１２の原画像の空間周波数スペクトルのうち照明光学系５の遮断周波数に近い空間周波数成分の振幅成分が減少する。そのため、試料１２の原画像のＳＮ比（コントラスト比）が低下する。本実施の形態では、第１デコンボリューション処理ステップ（ステップＳ３）及び第２デコンボリューション処理ステップ（ステップＳ４）を行っているため、照明光学系５の遮断周波数に近い空間周波数成分の振幅成分を回復させることができる。試料１２の補正画像Ｉ２ｆのＳＮ比（コントラスト比）が向上する。
【0077】
（実施の形態３）
【0078】
図１、図２及び図１１を参照して、実施の形態３の顕微鏡システム１を説明する。実施の形態３の顕微鏡システム１は、実施の形態２の顕微鏡システム１と同様に構成されているが、コンピュータ３の機能構成が異なっている。
【0079】
図１１に示されるように、本実施の形態の画像処理部３０は、実施の形態２の画像処理部３０と同様に構成されているが、デノイズ処理部３２が省略されている点で、実施の形態２の画像処理部３０と異なっている。
【0080】
累乗根画像算出部３１は、試料１２の原画像（例えば、二光子励起画像Ｉ２）から累乗根画像Ｉ１’を算出するステップ（ステップＳ１、図１２を参照）を行う。第１デコンボリューション処理部３３は、第１デコンボリューション処理ステップ（ステップＳ２、図１２及び図１３を参照）を行う。第２デコンボリューション処理部３４は、第２デコンボリューション処理ステップ（ステップＳ３、図１２及び図１４を参照）を行う。補正画像算出部３５は、試料１２の補正画像を算出するステップ（ステップＳ４、図１２を参照）を行う。
【0081】
図１２から図１４を参照して、実施の形態３の試料１２の原画像の補正方法を説明する。試料１２の補正画像を算出する方法は、画像処理部３０によって実行される。以下では、主に、顕微鏡２が二光子励起顕微鏡であり、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合について説明する。実施の形態３の試料１２の原画像の補正方法は、実施の形態２の試料１２の原画像の補正方法と同様であるが、主に、実施の形態２のデノイズ処理ステップ（Ｓ３）を備えていない点において、実施の形態２の試料１２の原画像の補正方法と異なっている。
【0082】
具体的には、図１２を参照して、本実施の形態のステップＳ１は、実施の形態２のステップＳ１と同じである。
【0083】
図１２に示されるように、累乗根画像Ｉ１’を算出するステップ（ステップＳ１）の後に、第１デコンボリューション処理ステップ（ステップＳ２）を行う。
【0084】
具体的には、図１３を参照して、実施の形態１の式（３）に示されるように、累乗根画像Ｉ１’をフーリエ変換して、累乗根画像Ｉ１’の空間周波数スペクトルＦ１を算出する（ステップＳ２－１）。本実施の形態のステップＳ２－１は、実施の形態１のステップＳ２－１と同じである。
【0085】
式（２２）に示されるように、累乗根画像Ｉ１’の空間周波数スペクトルＦ１を、光学系（照明光学系５）の光学伝達関数ＯＴＦ１でデコンボリューションして、累乗根画像Ｉ１’の補正空間周波数スペクトルＦ１ｄ’を算出する（ステップＳ２－２）。実施の形態２において既に記載したとおり、光学伝達関数ＯＴＦ１は、試料１２の原画像（例えば、二光子励起画像Ｉ２）を得る光学系（照明光学系５）の点像分布関数ＰＳＦ２の累乗根をフーリエ変換することによって算出される。第１デコンボリューション処理は、例えば、ＯＴＦ１をウィナーフィルタとして用いるデコンボリューション処理を含む。式（２２）におけるｗ１は、ウィナー定数を表す。
【数22】
【0086】
本実施の形態の目的は、顕微鏡２で得られた試料１２の原画像のＳＮ比（コントラスト比）を向上させることである。そのため、第１デコンボリューション処理は、光学伝達関数ＯＴＦ１で行われてもよいし、光学伝達関数ＯＴＦ１の絶対値である変調伝達関数ＭＴＦ１で行われてもよい。
【0087】
図１２に示されるように、第１デコンボリューション処理ステップ（ステップＳ２）の後に、第２デコンボリューション処理ステップ（ステップＳ３）を行う。
【0088】
【0089】
式（２４）に示されるように、累乗根補正画像Ｉ１ｄ’を累乗することによって、試料１２の暫定補正画像Ｉ２ｄを算出する（ステップＳ３－２）。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、累乗根補正画像Ｉ１ｄ’を二乗して、試料１２の暫定補正画像Ｉ２ｄを算出する。具体的には、累乗根補正画像Ｉ１ｄ’の各画素の光強度の二乗を算出することによって、試料１２の暫定補正画像Ｉ２ｄの各画素の光強度が算出される。試料１２の暫定補正画像Ｉ２ｄでは、第１デコンボリューション処理ステップ（ステップＳ２）によって照明光学系５の遮断周波数未満の周波数を有するノイズが低減されている。
【数24】
【0090】
式（２５）に示されるように、試料１２の暫定補正画像Ｉ２ｄをフーリエ変換して、試料１２の暫定補正画像Ｉ２ｄの空間周波数スペクトルＦ２ｄを算出する（ステップＳ３－３）。
【数25】
【0091】
式（２６）に示されるように、試料１２の暫定補正画像Ｉ２ｄの空間周波数スペクトルＦ２ｄを光学系の補正光学伝達関数ＯＴＦ２ｂでデコンボリューションすることによって、試料１２の原画像の補正空間周波数スペクトルＦ２ｅを算出する（ステップＳ３－４）。第２デコンボリューション処理は、例えば、ＯＴＦ２ｂをウィナーフィルタとして用いるデコンボリューション処理を含む。式（２６）におけるｗ２は、ウィナー定数を表す。
【数26】
【0092】
本実施の形態の目的は、顕微鏡２で得られた試料１２の原画像のＳＮ比（コントラスト比）を向上させることである。そのため、第２デコンボリューション処理は、補正光学伝達関数ＯＴＦ２ｂで行われてもよいし、補正光学伝達関数ＯＴＦ２ｂの絶対値である補正変調伝達関数ＭＴＦ２ｂで行われてもよい。
【0093】
補正光学伝達関数ＯＴＦ２ｂは、補正点像分布関数ＰＳＦ２ｂをフーリエ変換することによって算出される。具体的には、補正光学伝達関数ＯＴＦ２ｂは、式（２７）によって与えられる。式（２７）におけるＰＳＦ２ｂは、試料１２を二光子励起する照明光学系５の点像分布関数に第１デコンボリューション処理を行った補正点像分布関数を表す。
【数27】
【0094】
補正点像分布関数ＰＳＦ２ｂは、補正点像分布関数ＰＳＦ１ｂを累乗することによって算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、ＰＳＦ２ｂは、式（２８）によって与えられる。
【数28】
【0095】
補正点像分布関数ＰＳＦ１ｂは、試料１２の原画像を得る光学系の点像分布関数ＰＳＦ２の累乗根のフーリエ変換と高周波除去フィルタＧとの積を、光学伝達関数（光学伝達関数ＯＴＦ１）でデコンボリューションし、さらに逆フーリエ変換することによって算出される。例えば、試料１２の原画像が試料１２の二光子励起画像Ｉ２である場合、ＰＳＦ１ｂは、式（２９）によって与えられる。
【数29】
【0096】
本実施の形態の目的は、顕微鏡２で得られた試料１２の原画像のＳＮ比（コントラスト比）を向上させることである。そのため、補正点像分布関数ＰＳＦ１ｂを算出するために、点像分布関数ＰＳＦ２の累乗根のフーリエ変換と高周波除去フィルタＧとの積は、光学伝達関数ＯＴＦ１でデコンボリューションされてもよいし、光学伝達関数ＯＴＦ１の絶対値である変調伝達関数ＭＴＦ１でデコンボリューションされてもよい。
【0097】
【0098】
本実施の形態の効果は、実施の形態２の効果と同様の以下の効果を奏する。
【0099】
本実施の形態の顕微鏡画像補正方法は、試料１２の原画像（例えば、二光子励起画像Ｉ２）から試料１２の累乗根画像Ｉ１’を算出するステップ（ステップＳ１）と、試料１２の累乗根画像Ｉ１’から試料１２の累乗根画像Ｉ１’の第１空間周波数スペクトル（空間周波数スペクトルＦ１’）を算出するステップ（ステップＳ２－１）と、試料１２の累乗根画像Ｉ１’の第１空間周波数スペクトルを、試料１２の原画像を得る光学系（照明光学系５）の点像分布関数の累乗根をフーリエ変換することによって算出される光学系の光学伝達関数（光学伝達関数ＯＴＦ１）でデコンボリューションすることによって、試料１２の累乗根画像Ｉ１’の第１補正空間周波数スペクトル（補正空間周波数スペクトルＦ１ｄ’）を算出するステップ（ステップＳ２－２）と、試料１２の累乗根画像Ｉ１’の第１補正空間周波数スペクトルから試料１２の累乗根補正画像Ｉ１ｄ’を算出するステップ（ステップＳ３－１）と、試料１２の累乗根補正画像Ｉ１ｄ’を累乗することによって試料１２の暫定補正画像Ｉ２ｄを算出するステップ（ステップＳ３－２）と、試料１２の暫定補正画像Ｉ２ｄから試料１２の暫定補正画像Ｉ２ｄの第２空間周波数スペクトル（空間周波数スペクトルＦ２ｄ）を算出するステップ（ステップＳ３－３）と、試料１２の暫定補正画像Ｉ２ｄの第２空間周波数スペクトルを光学系の補正光学伝達関数（補正光学伝達関数ＯＴＦ２ｂ）でデコンボリューションすることによって、試料１２の原画像の第２補正空間周波数スペクトル（補正空間周波数スペクトルＦ２ｅ）を算出するステップ（ステップＳ３－４）とを備える。補正光学伝達関数は、第１補正点像分布関数（補正点像分布関数ＰＳＦ２ｂ）をフーリエ変換することによって算出される。第１補正点像分布関数は、第２補正点像分布関数（補正点像分布関数ＰＳＦ１ｂ）を累乗することによって算出される。第２補正点像分布関数は、試料１２の原画像を得る光学系の点像分布関数の累乗根のフーリエ変換と高周波除去フィルタＧとの積を、光学系の光学伝達関数（光学伝達関数ＯＴＦ１）でデコンボリューションし、さらに逆フーリエ変換することによって算出される。本実施の形態の顕微鏡画像補正方法は、試料１２の原画像の第２補正空間周波数スペクトルから試料１２の補正画像Ｉ２ｆを算出するステップ（ステップＳ４）とを備える。
【0100】
本実施の形態の顕微鏡画像補正プログラムは、本実施の形態の顕微鏡画像補正方法の各ステップをプロセッサ２０に実行させる。
【0101】
照明光学系５の光７が試料１２により吸収される際に光７の電界の振幅Ｅが二乗されることに起因して、試料１２の原画像の空間周波数スペクトルのうち照明光学系５の遮断周波数に近い空間周波数成分の振幅成分が減少する。そのため、試料１２の原画像のＳＮ比（コントラスト比）が低下する。本実施の形態では、第１デコンボリューション処理ステップ（ステップＳ２）及び第２デコンボリューション処理ステップ（ステップＳ３）を行っているため、照明光学系５の遮断周波数に近い空間周波数成分の振幅成分を回復させることができる。試料１２の補正画像Ｉ２ｆのＳＮ比（コントラスト比）が向上する。
【0102】
（実施例）
【0103】
【0104】
試料１２の原画像として、ＳＮ比（コントラスト比）が１０の第１原画像（図１５を参照）と、ＳＮ比が２００の第２原画像とを準備した。第１原画像及び第２原画像は、コンピュータ上で作成された点像に、コンピュータ上で作成されたガウシアン分布のノイズを追加することによって、作成されている。
【0105】
実施例１－１（図１６を参照）は、実施の形態１の一つの実施例であり、第１原画像を実施の形態１のようにデノイズ処理（図４を参照）して算出された補正画像である。実施例１－２は、実施の形態１の別の実施例であり、第２原画像を実施の形態１のようにデノイズ処理（図４を参照）して算出された補正画像である。
【0106】
実施例２－１（図１７を参照）は、実施の形態２の一つの実施例であり、第１原画像を実施の形態２のようにデノイズ処理、第１デコンボリューション処理及び第２デコンボリューション処理（図９を参照）して算出された補正画像である。実施例２－２は、実施の形態２の別の実施例であり、第２原画像を実施の形態２のようにデノイズ処理、第１デコンボリューション処理及び第２デコンボリューション処理（図９を参照）して算出された補正画像である。
【0107】
実施例３－１（図１８を参照）は、実施の形態３の一つの実施例であり、第１原画像を実施の形態３のように第１デコンボリューション処理及び第２デコンボリューション処理（図１２を参照）して算出された補正画像である。実施例３－２は、実施の形態３の別の実施例であり、第２原画像を実施の形態３のように第１デコンボリューション処理及び第２デコンボリューション処理（図１２を参照）して算出された補正画像である。
【0108】
図１９に、実施例１－１から実施例３－２の画像のＳＮ比の改善率を示す。実施例１－１の画像のＳＮ比の改善率は、実施例１－１の補正画像のＳＮ比を第１原画像のＳＮ比（すなわち１０）で除算することによって算出される。実施例１－２の画像のＳＮ比の改善率は、実施例１－２の補正画像のＳＮ比を第２原画像のＳＮ比（すなわち２００）で除算することによって算出される。実施例２－１の画像のＳＮ比の改善率は、実施例２－１の補正画像のＳＮ比を第１原画像のＳＮ比で除算することによって算出される。実施例２－２の画像のＳＮ比の改善率は、実施例２－２の補正画像のＳＮ比を第２原画像のＳＮ比で除算することによって算出される。実施例３－１の画像のＳＮ比の改善率は、実施例３－１の補正画像のＳＮ比を第１原画像のＳＮ比で除算することによって算出される。実施例３－２の画像のＳＮ比の改善率は、実施例３－２の補正画像のＳＮ比を第２原画像のＳＮ比で除算することによって算出される。
【0109】
実施例１－１及び実施例１－２のＳＮ比の改善率に示されるように、原画像に対して実施の形態１のデノイズ処理（図４を参照）を行うことによって、補正画像のＳＮ比が向上する。実施例３－１及び実施例３－２のＳＮ比の改善率に示されるように、原画像に対して実施の形態３の第１デコンボリューション処理及び第２デコンボリューション処理（図１２を参照）を行うことによって、補正画像のＳＮ比が大幅に向上する。実施例２－１及び実施例２－２のＳＮ比の改善率に示されるように、原画像に対して実施の形態２のデノイズ処理、第１デコンボリューション処理及び第２デコンボリューション処理（図９を参照）を行うことによって、補正画像のＳＮ比がさらに向上する。
【0110】
今回開示された実施の形態１から実施の形態３及びそれらの変形例はすべての点で例示であって制限的なものではないと考えられるべきである。例えば、顕微鏡２は、二光子励起顕微鏡以外の多光子励起顕微鏡、共焦点顕微鏡または超解像顕微鏡などであってもよい。試料１２の原画像は、二光子励起顕微鏡以外の多光子励起顕微鏡、共焦点顕微鏡または超解像顕微鏡などによって得られてもよい。本開示の範囲は、上記した説明ではなく特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることを意図される。
【符号の説明】
【0111】
１ 顕微鏡システム、２ 顕微鏡、３ コンピュータ、４ 表示装置、５ 照明光学系、６ 光源、７ 光、８ 走査ミラー、９ａ，９ｂ，９ｃ，９ｄ レンズ、１０ ダイクロイックミラー、１１ 対物レンズ、１２ 試料、１３ 試料ステージ、１４ 光、１５ 観察光学系、１６ 光検出器、２０ プロセッサ、２２ 記憶装置、２５ 原画像取得部、３０ 画像処理部、３１ 累乗根画像算出部、３２ デノイズ処理部、３３ 第１デコンボリューション処理部、３４ 第２デコンボリューション処理部、３５ 補正画像算出部。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
