JP2020038129A - 平均光子数の推定方法及び平均光子数の推定装置 - Google Patents

Abstract

【課題】平均光子数を精度よく推定することが可能な推定方法、及び平均光子数の推定装置を提供する。【解決手段】電子増倍型ＣＣＤ１１を用いて撮像された画像データから、画素毎に平均光子数を推定する方法であって、電子増倍型ＣＣＤ１１の特性に基づいた画素値の確率分布を算出する準備工程Ｓ１と、準備工程Ｓ１で算出された確率分布を用いて、画像データから画素毎に平均光子数を推定する推定工程Ｓ４と、を含み、推定工程Ｓ１は、ピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑの画素値の確率分布に基づいて、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑが、それぞれ、ｌ個、ｋ個、ｍ個の光子を検出した場合におけるピーク画素Ｐの画素値の確率分布を算出する。【選択図】図２

Description

本発明は、平均光子数の推定方法に関する。

従来、測定対象からの光を検出装置によって検出し、検出装置の出力値に基づいて測定対象から放出された光子の数を推定する定量方法が知られている。例えば、特許文献１には、光電子増倍管から出力される電流信号に基づいて、統計的な推定手法である最尤推定法を用いて光子数を推定することが記載されている。

極微弱光を検出するための検出装置としては、例えば非特許文献１に記載の電子増倍型ＣＣＤ（ＥＭＣＣＤ：Electron Multiplying Charge Coupled Devices）を搭載した撮像装置が知られている。

特開平９−１９６７５２号公報

A. G. Basden et al.,"Photon counting strategies with low-light-level CCDs" Mon. Not. R. Astron. Soc., 345, 985-991 (2003)

ところで、本発明者らは、極微弱光下でのイメージングのために、上記の電子増倍型ＣＣＤを用いて撮像された画像の画素値から平均光子数を推定することを検討している。しかしながら、電子増倍型ＣＣＤを用いて撮像された画像の画素値に対して、特許文献１に記載の最尤推定法を適用しても、平均光子数を精度よく推定することは困難であった。

本発明の一側面は、平均光子数を精度よく推定することが可能な推定方法、及び平均光子数の推定装置を提供することを目的とする。

本発明の一側面は、電子増倍型ＣＣＤを用いて撮像された画像データから、画素毎に平均光子数を推定する方法であって、電子増倍型ＣＣＤの特性に基づいた画素値の確率分布を算出する準備工程と、準備工程で算出された確率分布を用いて、画像データから画素毎に平均光子数を推定する推定工程と、を含み、準備工程は、電子増倍型ＣＣＤの１画素で１光子が検出された準備用画像データを取得する第１工程と、第１工程で取得した準備用画像データから、１光子が検出されたピーク画素と、電子増倍型ＣＣＤにおける電荷の転送方向においてピーク画素の上流側に隣接する第１隣接画素と、転送方向においてピーク画素の下流側に隣接する第２隣接画素とを含む少なくとも３つの連続する画素の画素値をそれぞれ取得する第２工程と、第２工程で取得したピーク画素、第１隣接画素、及び第２隣接画素の画素値に基づいてピーク画素で検出される光子の数が２以上且つ所定数以下であるそれぞれの場合について、当該ピーク画素、第１隣接画素、及び第２隣接画素に発生する画素値の確率分布をそれぞれ算出する第３工程と、第３工程で算出したピーク画素、第１隣接画素、及び第２隣接画素の画素値の確率分布に基づいて、第１隣接画素、ピーク画素、及び第２隣接画素が、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合におけるピーク画素の画素値の確率分布を算出する第４工程と、を有し、推定工程は、画像データを取得する第５工程と、第４工程で算出したピーク画素の画素値の確率分布を用いて、画像データの各画素に対する平均光子数の尤度を算出する第６工程と、平均光子数の尤度に基づいた目的関数を最大化することにより、平均光子数を推定する第７工程と、を有する。

この平均光子数の推定方法は、電子増倍型ＣＣＤの特性に基づいた画素値の確率分布を算出する準備工程と、準備工程で算出された確率分布を用いて、画像データから画素毎に平均光子数を推定する推定工程とを含んでいる。準備工程では、１光子が検出されたピーク画素と、電荷の転送方向においてピーク画素の上流側に隣接する第１隣接画素と、転送方向においてピーク画素の下流側に隣接する第２隣接画素とを含む少なくとも３つの連続する画素の画素値の確率分布に基づいて、第１隣接画素、ピーク画素、及び第２隣接画素が、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合におけるピーク画素の画素値の確率分布を算出する。これにより、準備工程で算出されたピーク画素の画素値の確率分布に、電子増倍型ＣＣＤにおける電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映される。推定工程では、電子増倍型ＣＣＤの特性による影響が反映された確率分布を用いて平均光子数の尤度を算出し、この尤度に基づいた目的関数を最大化することにより、平均光子数を推定する。このように、この平均光子数の推定方法では、電子増倍型ＣＣＤの特性の１つである電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映された統計モデルが用いられるので、平均光子数を精度よく推定することが可能である。

推定工程は、画像データの隣接する画素間の平均光子数の差分の二乗和を算出する工程を更に有し、第７工程において、平均光子数の差分の二乗和を目的関数に反映させてもよい。この場合、平均光子数を更に精度よく推定することができる。

第６工程において、平均光子数がポアソン分布に従うと仮定して尤度を算出してもよい。この場合においても、平均光子数を精度よく推定することができる。

第３工程において、準備用画像データのうち、光子が検出されていない画素の画素値の分布にガウス分布をフィッティングすることによって電子増倍型ＣＣＤの読み出しノイズの確率分布を更に算出し、第４工程において、第３工程で算出したピーク画素、第１隣接画素、及び第２隣接画素の画素値の確率分布と、読み出しノイズの確率分布とに基づいて、ピーク画素の画素値の確率分布を算出してもよい。この場合、準備工程で算出されたピーク画素の画素値の確率分布に、電子増倍型ＣＣＤにおける電荷の転送に伴う損失による画素値への影響に加え、電子増倍型ＣＣＤの読み出しノイズによる画素値への影響も反映される。したがって、平均光子数を更に精度よく推定することができる。

準備工程において、ピーク画素と、第１隣接画素と、第２隣接画素とを含む４つ以上の連続する画素の画素値に基づいてピーク画素の画素値の確率分布を算出してもよい。この場合、準備工程で算出されるピーク画素の画素値の確率分布の精度向上を図ることができる。したがって、平均光子数を更に精度よく推定することができる。

本発明の一側面に係る平均光子数の推定装置は、電子増倍型ＣＣＤと、電子増倍型ＣＣＤに接続されると共に、当該電子増倍型ＣＣＤを用いて撮像された画像データから画素毎に平均光子数を推定する推定演算部と、を備え、推定演算部は、電子増倍型ＣＣＤの１画素で１光子が検出された準備用画像データ、及び画像データを取得する画像データ取得手段と、準備用画像データから、１光子が検出されたピーク画素と、電子増倍型ＣＣＤにおける電荷の転送方向においてピーク画素の上流側に隣接する第１隣接画素と、転送方向においてピーク画素の下流側に隣接する第２隣接画素とを含む少なくとも３つの連続する画素の画素値をそれぞれ取得する画素値取得手段と、画素値取得手段によって取得されたピーク画素、第１隣接画素、及び第２隣接画素の画素値に基づいて、ピーク画素で検出される光子の数が２以上且つ所定数以下であるそれぞれの場合について、当該ピーク画素、第１隣接画素、及び第２隣接画素に発生する画素値の確率分布をそれぞれ算出する第１確率分布算出手段と、第１確率分布算出手段によって算出されたピーク画素、第１隣接画素、及び第２隣接画素の画素値の確率分布に基づいて、第１隣接画素、ピーク画素、及び第２隣接画素が、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合におけるピーク画素の画素値の確率分布を算出する第２確率分布算出手段と、第２確率分布算出手段によって算出されたピーク画素の画素値の確率分布を用いて、画像データの各画素に対する平均光子数の尤度を算出する尤度算出手段と、尤度算出手段によって算出された尤度に基づいた目的関数を最大化することにより、平均光子数を推定する推定手段と、を有する。

この平均光子数の推定装置の推定演算部は、１光子が検出されたピーク画素と、電荷の転送方向においてピーク画素の上流側に隣接する第１隣接画素と、転送方向においてピーク画素の下流側に隣接する第２隣接画素とを含む少なくとも３つの連続する画素の画素値の確率分布をそれぞれ算出する第１確率分布算出手段と、第１確率分布算出手段によって算出されたピーク画素、第１隣接画素、及び第２隣接画素の画素値の確率分布に基づいて、第１隣接画素、ピーク画素、及び第２隣接画素が、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合におけるピーク画素の画素値の確率分布を算出する第２確率分布算出手段とを有する。これにより、第２確率分布算出手段によって算出されたピーク画素の画素値の確率分布に、電子増倍型ＣＣＤにおける電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映される。推定演算部の尤度算出手段は、電子増倍型ＣＣＤの特性による影響が反映された確率分布を用いて平均光子数の尤度を算出し、推定演算部の推定手段は、尤度算出手段によって算出された尤度に基づいた目的関数を最大化することにより、平均光子数を推定する。このように、この平均光子数の推定装置では、電子増倍型ＣＣＤの特性の１つである電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映された統計モデルが用いられているので、平均光子数を精度よく推定することが可能である。

推定演算部は、画像データの隣接する画素間の平均光子数の差分の二乗和を算出する差分算出手段を更に備え、推定手段は、平均光子数の差分の二乗和を目的関数に反映させてもよい。この場合、平均光子数を更に精度よく推定することができる。

尤度算出手段は、平均光子数がポアソン分布に従うと仮定して尤度を算出してもよい。この場合においても、平均光子数を精度よく推定することができる。

推定演算部は、準備用画像データのうち、光子が検出されていない画素の画素値の分布にガウス分布をフィッティングすることによって電子増倍型ＣＣＤの読み出しノイズの確率分布を算出する読み出しノイズ算出手段を更に備え、第２確率分布算出手段は、第１確率分布算出手段によって算出されたピーク画素、第１隣接画素、及び第２隣接画素の画素値の確率分布と、読み出しノイズ算出手段によって算出された読み出しノイズの確率分布とに基づいて、ピーク画素の画素値の確率分布を算出してもよい。この場合、第２確率分布算出手段で算出されたピーク画素の画素値の確率分布に、電子増倍型ＣＣＤにおける電荷の転送に伴う損失による画素値への影響に加え、電子増倍型ＣＣＤの読み出しノイズによる画素値への影響も反映される。したがって、平均光子数を更に精度よく推定することができる。

画素値取得手段は、ピーク画素と、第１隣接画素と、第２隣接画素とを含む４つ以上の連続する画素の画素値を取得し、第１確率分布算出手段及び第２確率分布算出手段は、４つ以上の連続する画素の画素値に基づいて確率分布を算出してもよい。この場合、第２確率分布算出手段によって算出されるピーク画素の画素値の確率分布の精度向上を図ることができる。したがって、平均光子数を更に精度よく推定することができる。

本発明によれば、平均光子数を精度よく推定することが可能な推定方法、及び平均光子数の推定装置が提供される。

一実施形態に係る平均光子数の推定装置の概略構成図である。 一実施形態に係る平均光子数の推定方法を示すフローチャートである。 図２の準備工程における処理を示すフローチャートである。 図２の推定工程における処理を示すフローチャートである。 準備用画像データにおけるピーク画素、第１隣接画素、及び第２隣接画素を示す図である。 ピーク画素の画素値と、第１隣接画素の画素値との関係をプロットしたグラフである。 ピーク画素の画素値と、第２隣接画素の画素値との関係をプロットしたグラフである。 ピーク画素で検出された光子の数がｋ個である場合の第１隣接画素、ピーク画素、第２隣接画素の画素値の確率分布を示すグラフである。 本実施形態に係る平均光子数の推定方法による推定結果を説明するための図である。 本実施形態に係る平均光子数の推定方法による推定結果を説明するためのグラフである。 本実施形態に係る平均光子数の推定方法による推定結果を説明するためのグラフである。 比較例に係る平均光子数の推定方法による推定結果を説明するためのグラフである。

以下、図面を参照して種々の実施形態について詳細に説明する。なお、各図面において同一又は相当の部分に対しては同一の符号を付し、重複する説明を省略する。

図１は、本実施形態に係る平均光子数の推定装置の概略構成図である。図１に示される平均光子数の推定装置１は、例えば、細胞等の測定対象から発せられる極微弱光のイメージングにおいて、電子増倍型ＣＣＤを用いて撮像された画像データから平均光子数を推定するための装置である。図１に示されるように、平均光子数の推定装置１は、電子増倍型ＣＣＤ（ＥＭＣＣＤ：Electron Multiplying Charge Coupled Devices）１１が搭載された撮像部１０と、撮像部１０に接続された推定演算部２０とを備えている。

撮像部１０は、平均光子数の推定に用いられる画像データを取得する部分である。撮像部１０は推定演算部２０とデータ通信可能に接続されており、電子増倍型ＣＣＤ１１を用いて撮像された画像のデータを推定演算部２０に送信する。電子増倍型ＣＣＤ１１は、極微弱光を高いＳ／Ｎ比で検出するために用いられる検出装置である。電子増倍型ＣＣＤ１１は、一般的なＣＣＤが有する水平シフトレジスタ（不図示）に加え、入射光から変換された電子を増幅させる機能を有する電子増倍レジスタ（不図示）を有している。これにより、電子増倍型ＣＣＤ１１は、電子増倍レジスタに通常よりも高い電圧を印加して電子の転送を行うことにより、電子を所定の増倍率で増倍することが可能である。

推定演算部２０は、撮像部１０からのデータを受け取り、電子増倍型ＣＣＤを用いて撮像された画像データから画素毎に平均光子数を推定する部分である。推定演算部２０は、物理的には、プロセッサであるＣＰＵ（Central Processing Unit）と、記録媒体であるＲＡＭ（Random Access Memory）及びＲＯＭ（Read Only Memory）と、通信モジュールと、ディスプレイ、マウス、キーボード等の入出力デバイスとを含んだコンピュータ等の装置である。また、推定演算部２０は、機能的な構成として、画像データ取得手段２１と、画素値取得手段２２と、読み出しノイズ算出手段２３と、第１確率分布算出手段２４と、第２確率分布算出手段２５と、尤度算出手段２６と、差分算出手段２７と、推定手段２８と、を有している。

以下、図１〜４を参照して、推定演算部２０による処理を含めた平均光子数の推定装置１の動作を説明すると共に、本実施形態に係る平均光子数の推定方法について説明する。図２は、本実施形態に係る平均光子数の推定方法を示すフローチャートである。図３は、図２の平均光子数の推定方法の準備工程を示すフローチャートである。図４は、図２の平均光子数の推定方法の推定工程を示すフローチャートである。また、図５は、準備用画像データにおけるピーク画素、第１隣接画素、及び第２隣接画素を示す図である。図６は、ピーク画素の画素値と、第１隣接画素の画素値との関係をプロットしたグラフである。図７は、ピーク画素の画素値と、第２隣接画素の画素値との関係をプロットしたグラフである。図８は、ピーク画素で検出された光子の数がｋ個である場合の第１隣接画素、ピーク画素、第２隣接画素の画素値の確率分布を示すグラフである。

図２に示されるように、本実施形態に係る平均光子数の推定方法では、まず、電子増倍型ＣＣＤ１１の特性に基づいた画素値の確率分布を予め算出し（準備工程Ｓ１）、その結果をハードディスク等の記憶装置に保存する（工程Ｓ２）。その後、実際に平均光子数の推定を行う際に、保存された準備工程Ｓ１の結果をＲＡＭ等のメモリに読み出し（工程Ｓ３）、準備工程Ｓ１で算出された確率分布を用いて、画像データから画素毎に平均光子数を推定する（推定工程Ｓ４）。その後、次の画像データが存在するか否かを判断する（Ｓ５）。工程Ｓ５の判断の結果、次の画像データが存在した場合（Ｓ５−ＹＥＳ）には、再度推定工程Ｓ４を行い、画素毎に平均光子数を推定する。工程Ｓ５の判断の結果、次の画像データが存在しなかった場合（Ｓ５−ＮＯ）には、平均光子数の推定を終了する。

次に、図３を参照して、電子増倍型ＣＣＤ１１の特性に基づいた画素値の確率分布を算出する準備工程Ｓ１について詳細に説明する。準備工程Ｓ１では、まず、推定演算部２０の画像データ取得手段２１が、準備用画像データを取得する（第１工程Ｓ１１）。ここで、準備用画像データとは、電子増倍型ＣＣＤ１１の１画素で１光子が検出された箇所を含む画像のデータである。このような準備用画像データは、例えば、光源の輝度を低下させて、１画素で検出される光子数が高々１となるように入射光強度を調整した後、この入射光強度で電子増倍型ＣＣＤ１１を一様に照射して撮像することによって取得することができる。入射光強度は、例えば、全画素数（ここでは５１２×５１２画素）に対して、光子を検出した箇所の割合が１％程度となるように調整され得る。第１工程Ｓ１１では、準備用画像データを複数取得する。本実施形態では、１画素で１光子が検出された箇所が２９００個程度含まれている画像２０４８枚を１セットとし、このようなセットを１０セット取得する。この場合、１画素で１光子が検出された箇所の総数は、約６０００万となる。

次に、図５に示されるように、推定演算部２０の画素値取得手段２２（図１参照）が、第１工程Ｓ１１で取得した準備用画像データから、１光子が検出されたピーク画素Ｐと、第１隣接画素Ｏと、第２隣接画素Ｑと、を含む３つの連続する画素の画素値ｖ，ｕ，ｗをそれぞれ取得する（第２工程Ｓ１２）。第１隣接画素Ｏは、電子増倍型ＣＣＤ１１における電荷の転送方向Ａにおいてピーク画素Ｐの上流側に隣接する画素である。第２隣接画素Ｑは、転送方向Ａにおいてピーク画素Ｐの下流側に隣接する画素である。図５に示される例では、準備用画像データの一部において、３組のピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑの画素値が取得されている。第１隣接画素Ｏの画素値ｕ、ピーク画素Ｐの画素値ｖ、第２隣接画素Ｑの画素値ｗの組合わせ（ｕ，ｖ，ｗ）は、例えば、（６３，２２１，１４）、（３１，３４４，５）、（４７，３４８，４６）である。

なお、１光子が検出されたピーク画素Ｐの判別は、例えば画素値に基づいてなされ得る。１光子が検出された画素における画素値は、例えば、１５０より大きくなるので、このような画素値を有する画素をピーク画素Ｐと判断することができる。本実施形態では、画素値が１５０より大きい点に加え、転送方向Ａにおいて隣接する画素の画素値より大きい点、及び、転送方向Ａにおいて、他のピーク画素Ｐから１０画素以上離れている点を判別条件としている。上記の３つの条件を全て満たす画素をピーク画素Ｐと判別する。

次に、推定演算部２０の読み出しノイズ算出手段２３が、電子増倍型ＣＣＤ１１の読み出しノイズの確率分布を算出する（工程Ｓ１３）。具体的には、準備用画像データからオフセット値を差し引いた後、準備用画像データの画素のうち、光子を検出していない画素の画素値の分布に式（１）のガウス分布をフィッティングすることにより、読み出しノイズｚの確率分布の算出に必要となるσの値が求められる。式（１）のおけるσは標準偏差である。σの値としては、例えば、σ＝１６．５が得られる。

次に、第２工程Ｓ１２で取得したピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑの画素値ｖ，ｕ，ｗに基づいて、ピーク画素Ｐで検出される光子の数が２以上且つ所定数以下であるそれぞれの場合について、ピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑに発生する画素値の確率分布をそれぞれ算出する（第３工程Ｓ１４）。第３工程Ｓ１４は、推定演算部２０の第１確率分布算出手段２４によって行われる。

第３工程Ｓ１４における確率分布の算出方法について具体的に説明する。図６は、ピーク画素Ｐの画素値ｖと、第１隣接画素Ｏの画素値ｕとの関係を示す二次元ヒストグラムである。図７は、ピーク画素Ｐの画素値ｖと、第２隣接画素Ｑの画素値ｗとの関係を示す２次元ヒストグラムである。図６及び図７の二次元ヒストグラムは、第２工程Ｓ１２で取得した第１隣接画素Ｏの画素値ｕ、ピーク画素Ｐの画素値ｖ、第２隣接画素Ｑの画素値ｗの複数の組合わせ（ｕ，ｖ，ｗ）に基づいて作成され得る。図６及び図７においては、グレースケールの色が濃い座標ほど発生頻度が高く、色が薄い座標ほど発生頻度が低いことを示している。画素値ｕ，ｖ，ｗのそれぞれに読み出しノイズｚが付加される前の画素値を、それぞれ、画素値ｒ，ｓ，ｔとすると、画素値ｕ，ｖ，ｗは以下の式（２）のように表される。

ここで、ｘ＝ｒ＋ｓ＋ｔとすると、画素値ｒ，ｓ，ｔについて、それぞれ以下の式（３）〜（５）に示される関係式が導かれる。



式（３）〜（５）におけるα、β、γ、εは、図６及び図７のヒストグラムのデータに式（３）〜（５）を当てはめて求まる定数である。本実施形態における図６及び図７のヒストグラムのデータを用いた場合、α＝１９．１、β＝２１６．０、γ＝０．０４６７、ε＝０．０３７８が得られる。

ここで、ｘの確率分布は、以下の式（６）に示される指数分布で表され得る。式（６）におけるμは、電子増倍型ＣＣＤ１１における電子の増倍率である。

電子増倍型ＣＣＤ１１における電子の増倍率μは、第２工程Ｓ１２で取得した連続する３つの画素（第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑ）の画素値ｕ，ｖ，ｗの和のヒストグラムに対し、式（６）をフィッティングすることによって求めることができる。本実施形態の場合、増倍率μの値として、μ＝２１２．０が得られる。

次に、以上の式（１）〜（６）に示される第２工程Ｓ１２で取得したデータの解析結果に基づいて、読み出しノイズｚが付加される前における第１隣接画素Ｏの画素値ｒの確率分布ｆ_１（ｒ）、ピーク画素Ｐの画素値ｓの確率分布ｇ_１（ｓ）、及び第２隣接画素Ｑの画素値ｔの確率分布ｈ_１（ｔ）を算出する。まず、ｘの累積確率をΩ（ｘ）とすると、Ω（ｘ）は下記の式（７）で示される。

準備用画像データのある画素において１光子が検出された場合（すなわちピーク画素Ｐとなる画素）の、この画素に対して転送方向Ａの上流側の画素（すなわち第１隣接画素Ｏとなる画素）の、読み出しノイズｚが付加される前の画素値ｒの確率分布ｆ_１（ｒ）は、任意の非負整数値ｒに対して、下記の式（８）によって表される。なお、Ｒ^-１（ｒ）は、Ｒ（ｘ）の逆関数である。同様に、読み出しノイズｚが付加される前のピーク画素Ｐの画素値ｓの確率分布ｇ_１（ｓ）は式（９）によって表され、読み出しノイズｚが付加される前の第２隣接画素Ｑの画素値ｔの確率分布ｈ_１（ｔ）は式（１０）によって表される。

次に、式（８）〜（１０）に示される第１隣接画素Ｏの画素値ｒの確率分布ｆ_１（ｒ）、ピーク画素Ｐの画素値ｓの確率分布ｇ_１（ｓ）、及び第２隣接画素Ｑの画素値ｔの確率分布ｈ_１（ｔ）に基づいて、ピーク画素Ｐで検出された光子の数が２以上且つ所定数以下であるそれぞれの場合について、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑに発生する画素値の確率分布をそれぞれ算出する。ピーク画素Ｐで検出された光子の数がｋ個（ｋは２以上且つ所定数以下の整数）である場合の第１隣接画素Ｏ、ピーク画素Ｐ、第２隣接画素Ｑの画素値の確率分布は、それぞれ、畳み込み積分を使った以下の漸化式（１１）〜（１３）によって算出することができる。本実施形態では、ｋ＝２０までの各確率分布ｆ_ｋ（ｒ）、ｇ_ｋ（ｓ）、ｈ_ｋ（ｔ）を算出しており、その結果は図８（ａ）〜（ｃ）にそれぞれ示されている。なお、ｋの値の上限は２０に限定されず、測定環境に応じて適宜設定され得る。

なお、ピーク画素Ｐにおいて光子が検出されていない場合、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑの画素値は何れも０となる。この場合の第１隣接画素Ｏの画素値ｒ、ピーク画素Ｐの画素値ｓ、第２隣接画素Ｑの画素値ｔの確率分布は、それぞれ、ｆ_０（ｒ）＝δ（ｒ）、ｇ_０（ｓ）＝δ（ｓ）、ｈ_０（ｔ）＝δ（ｔ）となる。ここで、δ（・）はディラックのデルタ関数である。

次に、第３工程で算出した第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑの画素値の確率分布ｆ_ｋ（ｒ）、ｇ_ｋ（ｓ）、ｈ_ｋ（ｔ）に基づいて、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑが、それぞれ、ｌ個、ｋ個、ｍ個の光子を検出した場合におけるピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する（第４工程Ｓ１５）。ここで、ｌ、ｋ、ｍはそれぞれ１以上の整数である。第４工程Ｓ１５は、推定演算部２０の第２確率分布算出手段２５によって行われる。ここで、ｌ個の光子を検出した画素（すなわち第１隣接画素Ｏ）を基準とした場合、ｋ個の光子を検出した画素（ピーク画素Ｐ）は、ｌ個の光子を検出した画素の第２隣接画素に相当する。また、ｍ個の光子を検出した画素（すなわち第２隣接画素Ｑ）を基準とした場合、ｋ個の光子を検出した画素（ピーク画素Ｐ）は、ｍ個の光子を検出した画素の第１隣接画素に相当する。したがって、ピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）は、式（１）に示される読み出しノイズｚを与えるガウス分布ｎ（ｚ）、及び、第３工程Ｓ１４で算出した式（１１）〜（１３）に示される各確率分布ｆ_ｋ（ｒ）、ｇ_ｋ（ｓ）、ｈ_ｋ（ｔ）を用いて、以下の式（１４）に示す畳み込み積分によって算出される。本実施形態では、例えば、０から２０までのｋ、ｌ、ｍの各組合わせについて、確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する。以上の工程により、電子増倍型ＣＣＤ１１の特性に基づいた画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）が算出され、準備工程が完了する。

続いて、推定工程Ｓ４について詳細に説明する。推定工程Ｓ４では、まず、電子増倍型ＣＣＤ１１を用いて、平均光子数の推定に用いられる画像データを取得する（第５工程Ｓ４１）。平均光子数の推定に用いられる画像データの取得は、画像データ取得手段２１によって行われ、画像データの画素値の取得は、画素値取得手段２２によって行われる。以下の説明では、第５工程Ｓ４１で取得した画像データを｛ｙ_ｉ，ｊ｝と標記し、画像データ｛ｙ_ｉ，ｊ｝に含まれる１の画素を（ｉ，ｊ）、画素（ｉ，ｊ）における画素値をｙ_ｉ，ｊとして説明する場合がある。なお、ｉは電子増倍型ＣＣＤ１１における電子の転送方向Ａにおける画素の位置を示し、ｊは転送方向Ａに交差する方向における画素の位置を示している。

次に、準備工程Ｓ１で算出した画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を用いて、画像データ｛ｙ_ｉ，ｊ｝の各画素（ｉ，ｊ）に対する平均光子数の尤度を算出する（第６工程Ｓ４２）。第６工程Ｓ４２は、推定演算部２０の尤度算出手段２６によって行われる。第６工程Ｓ４２では、まず、画素（ｉ，ｊ）において検出される光子数の平均がλ_ｉ，ｊであるときに、画素（ｉ，ｊ）で光子が検出されたことによる画素値ｙ_ｉ，ｊへの寄与分をｓとし、寄与分ｓの確率分布Ｇ（ｓ）を計算する。画素（ｉ，ｊ）で検出される光子数がポアソン分布に従うと仮定すると、確率分布Ｇ（ｓ）は下記の式（１５）によって算出される。ここで、式（１５）におけるＮは、準備工程Ｓ１において確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する際に用いたｌ、ｋ、ｍの上限値である。本実施形態においては、Ｎの値は例えば２０である。

次に、画素（ｉ−１，ｊ）において検出される光子数の平均がλ_{ｉ−１，ｊ}であるときに、画素（ｉ−１，ｊ）で光子が検出されたことによる画素値ｙ_ｉ，ｊへの寄与分をｔとし、寄与分ｔの確率分布Ｈ（ｔ）を計算する。画素（ｉ−１，ｊ）で検出される光子数がポアソン分布に従うと仮定すると、確率分布Ｈ（ｔ）は下記の式（１６）によって算出される。

更に、画素（ｉ＋１，ｊ）において検出される光子数の平均がλ_{ｉ＋１，ｊ}であるときに、画素（ｉ＋１，ｊ）で光子が検出されたことによる画素値ｙ_ｉ，ｊへの寄与分をｒとし、寄与分ｒの確率分布Ｆ（ｒ）を計算する。画素（ｉ＋１，ｊ）で検出される光子数がポアソン分布に従うと仮定すると、確率分布Ｆ（ｒ）は下記の式（１７）によって算出される。

以上の式（１５）〜（１７）により、連続する３つの画素（ｉ−１，ｊ）、（ｉ，ｊ）、（ｉ＋１，ｊ）において検出される光子数の平均が、それぞれ、λ_{ｉ−１，ｊ}、λ_ｉ，ｊ、λ_{ｉ＋１，ｊ}である場合、画素（ｉ，ｊ）の画素値ｙ_ｉ，ｊの確率分布はｙ_ｉ，ｊ＝ｓ＋ｔ＋ｒ＋ｚであることから、下記の式（１８）で表される。ここで、ｚはガウス分布ｎ（ｚ）に従う読み出しノイズである。

そして、画像データ｛ｙ_ｉ，ｊ｝が得られる確率は下記の式（１９）で表され、画像データ｛ｙ_ｉ，ｊ｝に対する対数尤度は、下記の式（２０）で表される。

次に、画像データ｛ｙ_ｉ，ｊ｝の隣接する画素間の平均光子数λ_ｉ，ｊの差分の二乗和を算出する（工程Ｓ４３）。本実施形態では、転送方向Ａにおいて隣り合う画素間、及び転送方向Ａに直交する方向において隣り合う画素間の平均光子数λ_ｉ，ｊの差分の二乗和をそれぞれ算出する。この工程Ｓ４３は、差分算出手段２７によって行われる。隣り合う画素間の平均光子数λ_ｉ，ｊの差分の二乗和は、下記の式（２１）によって算出することができる。

次に、第６工程Ｓ４２で算出した式（２０）によって表される平均光子数の尤度、及び、工程Ｓ４３で算出した式（２１）によって表される平均光子数λ_ｉ，ｊの差分の二乗和に基づいて、目的関数を求める。本実施形態では、隣接する画素間では、検出される光子数は大きくは異ならないと仮定して、差分の二乗和ψをペナルティとして対数尤度に反映させたものを目的関数とし、この目的関数を最大化することによって平均光子数画像｛λ_ｉ，ｊ｝を推定する（第７工程Ｓ４４）。目的関数は、下記の式（２２）によって表される。なお、本実施形態における実際の計算としては、式（２２）に−１を乗算して目的関数とし、共役勾配法により最小化して計算を行う。この計算式は、下記の式（２３）によって示される。以上の工程により、画像データ｛ｙ_ｉ，ｊ｝から平均光子数画像｛λ_ｉ，ｊ｝が推定され、推定工程Ｓ４が完了する。

続いて、図９〜１１を参照して、本実施形態に係る平均光子数の推定装置１及び推定方法による平均光子数画像の推定結果について説明する。図９〜１１は、本実施形態に係る平均光子数の推定方法による推定結果を説明するためのグラフである。図９（ａ）は、測定対象となるグレースケールパターンを示している。図９（ｂ）は、電子増倍型ＣＣＤ１１により得られた画像データを示している。図９（ｃ）は、図９（ｂ）の画像データから推定された平均光子数画像を示している。図９（ｂ）に示される画像データは、パルス点灯するＬＥＤからの極微弱光によって照明された図９（ａ）に示されるパターンを、電子増倍型ＣＣＤ１１を用いて撮像したものである。なお、本実施形態で用いられた電子増倍型ＣＣＤ１１の画素数は５１２×５１２画素であるが、図９（ｂ）では、中央部の２５６×２５６画素の画像を示している。

図１０は、図９（ｃ）のラインＬ１上の推定値のプロファイルを示している。ラインＬ１は、ｊ＝７８の一列の画素上に位置している。ここでは、ＬＥＤの点灯パルス幅を２００μ秒、４００μ秒、８００μ秒として、それぞれ１００枚の平均光子数画像の推定を行い、ｊ＝７８の一列の画素の推定値がプロットされている。分散したそれぞれの点が各１００枚の平均光子数画像からの推定値を示しており、一点鎖線、点線、実線は、それぞれ、パルス幅が２００μ秒、４００μ秒、８００μ秒の場合の平均値を示している。ＬＥＤの点灯パルスとしては矩形パルスが用いられているので、電子増倍型ＣＣＤ１１によって検出される光強度は、ＬＥＤのパルス幅に比例する。図１０に示される結果から、照明用ＬＥＤの点灯時間（すなわち光強度）に比例した平均光子数画像が推定されていることが確認できる。

図１１は、図１０のラインＬ２上の推定値の平均値とＬＥＤのパルス幅との関係を示している。ラインＬ２は、ｉ＝７８の一列の画素上に位置している。すなわち、図１１は、画素（７８，７８）における推定値の平均値とＬＥＤのパルス幅との関係を示している。図１１に示されるように、画素（７８，７８）における推定値の平均値とパルス幅とは、良好な直線関係にある。これにより、本実施形態に係る平均光子数の推定装置１及び推定方法によって、平均光子数を定量的に精度よく推定できていることが確認できる。

次に、図１２を参照して、比較例に係る平均光子数の推定方法による推定結果について説明する。図１２は、比較例に係る平均光子数の推定方法による推定結果を説明するための図である。比較例に係る平均光子数の推定方法では、電子増倍型ＣＣＤ１１の１画素の画素値の分布のみに基づいて、その画素における平均光子数を推定している点で、本発明に係る平均光子数の推定方法と相違している。

図１２は、図１１と同様に、ＬＥＤの点灯パルス幅を２００μ秒、４００μ秒、８００μ秒として、それぞれ８１７５枚の画像を取得し、画素（７８，７８）における画素値の分布から比較例に係る平均光子数の推定方法により推定された平均光子数とＬＥＤのパルス幅との関係を示したものである。図１２に示されるように、比較例に係る平均光子数の推定方法では、画素（７８，７８）における平均光子数の推定値とパルス幅とは直線関係になく、平均光子数を正確に推定できていないことが確認できる。図１２に示される結果と図１１に示される結果とを比較すると、本実施形態に係る平均光子数の推定方法では、比較例に係る平均光子数の推定方法より正確に平均光子数を推定できていることが確認できる。

以上説明したように、本実施形態に係る平均光子数の推定方法は、電子増倍型ＣＣＤ１１の特性に基づいた画素値の確率分布を算出する準備工程Ｓ１と、準備工程Ｓ１で算出された確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を用いて、画像データから画素毎に平均光子数を推定する推定工程Ｓ４とを含んでいる。準備工程Ｓ１では、１光子が検出されたピーク画素Ｐと、電荷の転送方向Ａにおいてピーク画素Ｐの上流側に隣接する第１隣接画素Ｏと、転送方向Ａにおいてピーク画素Ｐの下流側に隣接する第２隣接画素Ｑとを含む少なくとも３つの連続する画素の画素値の確率分布に基づいて、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑが、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合におけるピーク画素Ｐの画素値の確率分布を算出する。これにより、準備工程Ｓ１で算出されたピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）に、電子増倍型ＣＣＤ１１における電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映される。推定工程Ｓ４では、電子増倍型ＣＣＤ１１の特性による影響が反映された確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を用いて平均光子数の尤度を算出し、この尤度に基づいた目的関数を最大化することにより、平均光子数を推定する。このように、この平均光子数の推定方法では、電子増倍型ＣＣＤ１１の特性の１つである電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映された統計モデルが用いられるので、平均光子数を精度よく推定することが可能である。

また、推定工程Ｓ４は、画像データの隣接する画素間の平均光子数の差分の二乗和ψを算出する工程Ｓ４３を更に有し、第７工程Ｓ４４において、平均光子数の差分の二乗和ψを目的関数に反映させる。これにより、平均光子数を更に精度よく推定することができる。

また、本実施形態に係る平均光子数の推定方法では、第６工程Ｓ４２において、平均光子数がポアソン分布に従うと仮定して尤度を算出する。この場合においても、平均光子数を精度よく推定することができる。

また、本実施形態に係る平均光子数の推定方法では、準備用画像データのうち、光子が検出されていない画素の画素値の分布にガウス分布をフィッティングすることによって電子増倍型ＣＣＤ１１の読み出しノイズの確率分布を更に算出し（工程Ｓ１３）、第４工程Ｓ１５において、第３工程Ｓ１４で算出したピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑの画素値の確率分布と、読み出しノイズの確率分布とに基づいて、ピーク画素の画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する。これにより、準備工程Ｓ１で算出されたピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）に、電子増倍型ＣＣＤ１１における電荷の転送に伴う損失による画素値への影響に加え、電子増倍型ＣＣＤ１１の読み出しノイズによる画素値への影響も反映される。したがって、平均光子数を更に精度よく推定することができる。

本実施形態に係る平均光子数の推定装置１の推定演算部２０は、１光子が検出されたピーク画素Ｐと、電荷の転送方向Ａにおいてピーク画素Ｐの上流側に隣接する第１隣接画素Ｏと、転送方向Ａにおいてピーク画素Ｐの下流側に隣接する第２隣接画素Ｑとを含む少なくとも３つの連続する画素の画素値の確率分布をそれぞれ算出する第１確率分布算出手段２４と、第１確率分布算出手段２４によって算出されたピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑの画素値の確率分布に基づいて、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑが、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合におけるピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する第２確率分布算出手段２５とを有する。これにより、第２確率分布算出手段２５によって算出されたピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）に、電子増倍型ＣＣＤ１１における電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映される。推定演算部２０の尤度算出手段２６は、電子増倍型ＣＣＤ１１の特性による影響が反映された確率分布を用いて平均光子数の尤度を算出し、推定演算部２０の推定手段２８は、尤度算出手段２６によって算出された尤度に基づいた目的関数を最大化することにより、平均光子数を推定する。このように、この平均光子数の推定装置１では、電子増倍型ＣＣＤ１１の特性の１つである電荷の転送に伴う損失（Charge Transfer Inefficiency）による画素値への影響が反映された統計モデルが用いられているので、平均光子数を精度よく推定することが可能である。

また、推定演算部２０は、画像データの隣接する画素間の平均光子数の差分の二乗和ψを算出する差分算出手段２７を更に備え、推定手段２８は、平均光子数の差分の二乗和ψを目的関数に反映させる。これにより、平均光子数を更に精度よく推定することができる。

また、尤度算出手段２６は、平均光子数がポアソン分布に従うと仮定して尤度を算出する。この場合においても、平均光子数を精度よく推定することができる。

また、推定演算部２０は、準備用画像データのうち、光子が検出されていない画素の画素値の分布にガウス分布をフィッティングすることによって電子増倍型ＣＣＤ１１の読み出しノイズの確率分布を算出する読み出しノイズ算出手段２３を更に備え、第２確率分布算出手段２５は、第１確率分布算出手段２４によって算出されたピーク画素Ｐ、第１隣接画素Ｏ、及び第２隣接画素Ｑの画素値の確率分布と、読み出しノイズ算出手段２３によって算出された読み出しノイズの確率分布とに基づいて、ピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する。これにより、第２確率分布算出手段２５で算出されたピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）に、電子増倍型ＣＣＤ１１における電荷の転送に伴う損失による画素値への影響に加え、電子増倍型ＣＣＤ１１の読み出しノイズによる画素値への影響も反映される。したがって、平均光子数を更に精度よく推定することができる。

以上、本発明の実施形態について説明してきたが、本発明は上記の実施形態に限定されず、種々の変更を行うことができる。例えば、上記の実施形態では、第１隣接画素Ｏ、ピーク画素Ｐ、及び第２隣接画素Ｑの連続する３つの画素に基づいて確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する例について説明したが、準備工程Ｓ１において、第１隣接画素Ｏと、ピーク画素Ｐと、第２隣接画素Ｑとを含む４つ以上の連続する画素の画素値に基づいてピーク画素の画素値の確率分布を算出してもよい。一例として、転送方向Ａにおいて第１隣接画素Ｏの上流側に隣接する第３隣接画素Ｒを更に含む４つの画素に基づいて確率分布を算出する場合、第３工程Ｓ１４において、第１確率分布算出手段２４によって第３隣接画素Ｒの画素値の確率分布ε_ｋ（ｑ）を求める。そして、第４工程Ｓ１５において、第２確率分布算出手段２５によって、転送方向Ａに沿った４つの４つの連続した画素が、それぞれ、ｌ個、ｋ個、ｍ個、ｏ個の光子を検出した場合におけるｋ個の光子を検出した画素の画素値の確率分布φ_{k，ｌ，ｍ，ｏ}（ｙ）を下記の式（２４）により算出する。なお、ｌ個の光子を検出した画素を基準とした場合、ｋ個の光子を検出した画素は、ｌ個の光子を検出した画素の第２隣接画素に相当する。また、ｍ個の光子を検出した画素を基準とした場合、ｋ個の光子を検出した画素は、ｍ個の光子を検出した画素の第１隣接画素に相当する。さらに、ｏ個の光子を検出した画素を基準とした場合、ｋ個の光子を検出した画素は、ｏ個の光子を検出した画素の第３隣接画素に相当する。したがって、ピーク画素Ｐの画素値の確率分布φ_{k，ｌ，ｍ，ｏ}（ｙ）は下記の式（２４）によって表される。その後、第６工程Ｓ４２において、尤度算出手段２６によって対数尤度を算出する。この場合の対数尤度は、下記の式（２５）及び式（２６）によって算出され得る。この場合、準備工程Ｓ１において第２確率分布算出手段２５によって算出されるピーク画素Ｐの画素値の確率分布の精度向上を図ることができる。したがって、平均光子数を更に精度よく推定することができる。

また、上記の実施形態では、推定工程Ｓ４が画像データの隣接する画素間の平均光子数の差分の二乗和を算出する工程Ｓ４３を更に有し、第７工程Ｓ４４において、平均光子数の差分の二乗和ψをペナルティとして目的関数に反映させる例について説明したが、目的関数に反映させる関数は差分の二乗和に限定されない。例えば、下記の式（２７）によって表される関数をペナルティとして目的関数に反映させてもよい。なお、式（２７）におけるξは、ペナルティの増加傾向の変曲点の位置を調整するパラメータである。上記の式（２１）で求められる差分の二乗和ψをペナルティとする場合、差分が大きくなるほどペナルティも大きくなる。これに対し、下記の式（２７）で表される関数をペナルティとする場合、ξに比べて差分が十分に小さいときには差分の二乗に比例したペナルティとなり、差分がξよりも大きくなると、ペナルティの増加は式（２１）に比べて緩やかになる。このため、式（２７）で表される関数は、画像データが平均光子数が急峻に変化するエッジを有する場合に、好適に適用され得る。

また、上記の実施形態では、準備工程Ｓ１によって確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出した後に推定工程Ｓ４を行う例について説明したが、同一の電子増倍型ＣＣＤ１１を用いて撮像した複数の画像データから平均光子数を推定する場合には、確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を１度だけ算出し、その後推定工程Ｓ４のみを行ってもよい。

また、上記の実施形態では、第６工程Ｓ４２において、尤度算出手段２６によって平均光子数がポアソン分布に従うと仮定して尤度を算出する例について説明したが、平均光子数が他の確率分布に従うと仮定してもよい。

また、上記の実施形態では、ガウス分布をフィッティングすることによって電子増倍型ＣＣＤ１１の読み出しノイズの確率分布を算出する例について説明したが、読み出しノイズは、他の確率分布をフィッティングすることによって算出されてもよい。

また、第４工程Ｓ１５において、第３工程Ｓ１４で算出した画素値の確率分布と、読み出しノイズの確率分布とに基づいて、ピーク画素Ｐの画素値の確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出する例について説明したが、第４工程Ｓ１５では、第３工程Ｓ１４で算出した画素値の確率分布のみに基づいて確率分布φ_{ｋ，ｌ，ｍ}（ｙ）を算出してもよい。

１…推定装置、１０…撮像部、１１…電子増倍型ＣＣＤ、２０…推定演算部、２１…画像データ取得手段、２２…画素値取得手段、２３…読み出しノイズ算出手段、２４…第１確率分布算出手段、２５…第２確率分布算出手段、２６…尤度算出手段、２７…差分算出手段、２８…推定手段、Ａ…転送方向、Ｏ…第１隣接画素、Ｐ…ピーク画素、Ｑ…第２隣接画素、Ｓ１…準備工程、Ｓ４…推定工程、Ｓ１１…第１工程、Ｓ１２…第２工程、Ｓ１４…第３工程、Ｓ１５…第４工程、Ｓ４１…第５工程、Ｓ４２…第６工程、Ｓ４４…第７工程。

Claims (10)

1. 電子増倍型ＣＣＤを用いて撮像された画像データから、画素毎に平均光子数を推定する方法であって、
   前記電子増倍型ＣＣＤの特性に基づいた画素値の確率分布を算出する準備工程と、
   前記準備工程で算出された前記確率分布を用いて、前記画像データから画素毎に平均光子数を推定する推定工程と、を含み、
   前記準備工程は、
   前記電子増倍型ＣＣＤの１画素で１光子が検出された準備用画像データを取得する第１工程と、
   前記第１工程で取得した準備用画像データから、１光子が検出されたピーク画素と、前記電子増倍型ＣＣＤにおける電荷の転送方向において前記ピーク画素の上流側に隣接する第１隣接画素と、前記転送方向において前記ピーク画素の下流側に隣接する第２隣接画素とを含む少なくとも３つの連続する画素の画素値をそれぞれ取得する第２工程と、
   前記第２工程で取得した前記ピーク画素、前記第１隣接画素、及び前記第２隣接画素の画素値に基づいて前記ピーク画素で検出される光子の数が２以上且つ所定数以下であるそれぞれの場合について、当該ピーク画素、前記第１隣接画素、及び前記第２隣接画素に発生する画素値の確率分布をそれぞれ算出する第３工程と、
   前記第３工程で算出した前記ピーク画素、前記第１隣接画素、及び前記第２隣接画素の画素値の確率分布に基づいて、前記第１隣接画素、前記ピーク画素、及び前記第２隣接画素が、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合における前記ピーク画素の画素値の確率分布を算出する第４工程と、を有し、
   前記推定工程は、
   前記画像データを取得する第５工程と、
   前記第４工程で算出した前記ピーク画素の画素値の確率分布を用いて、前記画像データの各画素に対する前記平均光子数の尤度を算出する第６工程と、
   前記平均光子数の前記尤度に基づいた目的関数を最大化することにより、前記平均光子数を推定する第７工程と、を有する、平均光子数の推定方法。
2. 前記推定工程は、前記画像データの隣接する画素間の前記平均光子数の差分の二乗和を算出する工程を更に有し、
   前記第７工程において、前記平均光子数の差分の二乗和を前記目的関数に反映させる、請求項１に記載の平均光子数の推定方法。
3. 前記第６工程において、前記平均光子数がポアソン分布に従うと仮定して前記尤度を算出する、請求項１又は２に記載の平均光子数の推定方法。
4. 前記第３工程において、前記準備用画像データのうち、光子が検出されていない画素の画素値の分布にガウス分布をフィッティングすることによって前記電子増倍型ＣＣＤの読み出しノイズの確率分布を更に算出し、
   前記第４工程において、前記第３工程で算出した前記ピーク画素、前記第１隣接画素、及び前記第２隣接画素の画素値の確率分布と、前記読み出しノイズの確率分布とに基づいて、前記ピーク画素の画素値の確率分布を算出する、請求項１〜３の何れか一項に記載の平均光子数の推定方法。
5. 前記準備工程において、前記ピーク画素と、前記第１隣接画素と、前記第２隣接画素とを含む４つ以上の連続する画素の画素値に基づいて前記ピーク画素の画素値の確率分布を算出する、請求項１〜４の何れか一項に記載の平均光子数の推定方法。
6. 電子増倍型ＣＣＤが搭載された撮像部と、
   前記撮像部に接続されると共に、前記電子増倍型ＣＣＤを用いて撮像された画像データから画素毎に平均光子数を推定する推定演算部と、を備え、
   前記推定演算部は、
   前記電子増倍型ＣＣＤの１画素で１光子が検出された準備用画像データ、及び前記画像データを取得する画像データ取得手段と、
   前記準備用画像データから、１光子が検出されたピーク画素と、前記電子増倍型ＣＣＤにおける電荷の転送方向において前記ピーク画素の上流側に隣接する第１隣接画素と、前記転送方向において前記ピーク画素の下流側に隣接する第２隣接画素とを含む少なくとも３つの連続する画素の画素値をそれぞれ取得する画素値取得手段と、
   前記画素値取得手段によって取得された前記ピーク画素、前記第１隣接画素、及び前記第２隣接画素の画素値に基づいて、前記ピーク画素で検出される光子の数が２以上且つ所定数以下であるそれぞれの場合について、当該ピーク画素、前記第１隣接画素、及び前記第２隣接画素に発生する画素値の確率分布をそれぞれ算出する第１確率分布算出手段と、
   前記第１確率分布算出手段によって算出された前記ピーク画素、前記第１隣接画素、及び前記第２隣接画素の画素値の確率分布に基づいて、前記第１隣接画素、前記ピーク画素、及び前記第２隣接画素が、それぞれ、ｌ個、ｋ個、ｍ個（ｌ、ｋ、ｍは１以上の整数）の光子を検出した場合における前記ピーク画素の画素値の確率分布を算出する第２確率分布算出手段と、
   前記第２確率分布算出手段によって算出された前記ピーク画素の画素値の確率分布を用いて、前記画像データの各画素に対する前記平均光子数の尤度を算出する尤度算出手段と、
   前記尤度算出手段によって算出された前記尤度に基づいた目的関数を最大化することにより、前記平均光子数を推定する推定手段と、を有する、平均光子数の推定装置。
7. 前記推定演算部は、前記画像データの隣接する画素間の前記平均光子数の差分の二乗和を算出する差分算出手段を更に備え、
   前記推定手段は、前記平均光子数の差分の二乗和を前記目的関数に反映させる、請求項６に記載の平均光子数の推定装置。
8. 前記尤度算出手段は、前記平均光子数がポアソン分布に従うと仮定して前記尤度を算出する、請求項６又は７に記載の平均光子数の推定装置。
9. 前記推定演算部は、前記準備用画像データのうち、光子が検出されていない画素の画素値の分布にガウス分布をフィッティングすることによって前記電子増倍型ＣＣＤの読み出しノイズの確率分布を算出する読み出しノイズ算出手段を更に備え、
   前記第２確率分布算出手段は、前記第１確率分布算出手段によって算出された前記ピーク画素、前記第１隣接画素、及び前記第２隣接画素の画素値の確率分布と、前記読み出しノイズ算出手段によって算出された前記読み出しノイズの確率分布とに基づいて、前記ピーク画素の画素値の確率分布を算出する、請求項６〜８の何れか一項に記載の平均光子数の推定装置。
10. 前記画素値取得手段は、前記ピーク画素と、前記第１隣接画素と、前記第２隣接画素とを含む４つ以上の連続する画素の画素値を取得し、
    前記第１確率分布算出手段及び前記第２確率分布算出手段は、前記４つ以上の連続する画素の画素値に基づいて確率分布を算出する、請求項６〜９の何れか一項に記載の平均光子数の推定装置。