JP6051393B2 - ポアソンノイズ除去装置、ポアソンノイズ除去方法、及びｘ線撮影装置 - Google Patents

Description

本発明は、画像のポアソンノイズを除去するポアソンノイズ除去装置及びポアソンノイズ除去方法並びにポアソンノイズ除去装置を備えるＸ線撮影装置に関する。

Ｘ線撮影で得られる画像に含まれるノイズはポアソンノイズが支配的であることが知られている。ポアソン分布は平均と分散が常に一致するため、Ｘ線撮影の際に被写体に照射されるＸ線量子の個数が多いほど、Ｘ線撮影で得られる画像に含まれるポアソンノイズの大きさも大きくなる。

ノイズには、信号の大きさに依存して大きさが変動するノイズと、信号の大きさに関わらず大きさが一定であるノイズとの二種類がある。ポアソンノイズは前者である。

ノイズによる絶対的な変動は分散として表現される。これに対し、ノイズの変動に対する相対的な信号強度として信号対雑音比と呼ばれる量が定義される。信号対雑音比が大きいほど、その信号に含まれるノイズの割合は小さいことを意味する。

Ｘ線撮影の際に被写体に照射されるＸ線量子の個数が多くなると、Ｘ線撮影で得られる画像に含まれるポアソンノイズの絶対的な大きさ（分散）は上述した通り大きくなるが、信号強度はそれ以上に大きくなるため、信号強度に対するポアソンノイズの相対的な大きさ（信号対雑音比の逆数）は小さくなることが知られている。したがって、Ｘ線撮影で得られる画像に含まれるポアソンノイズの分散（絶対的な大きさ）は、Ｘ線撮影の際に被写体に照射されるＸ線量子の個数が多いほど大きくなるのに対して、Ｘ線撮影で得られる画像の信号対雑音比は、Ｘ線撮影の際に被写体に照射されるＸ線量子の個数が少ないほど悪く（小さく）なる。

ここで、水中において横方向に並べた２本のアルミ棒をＸ線撮影した場合の測定画像及びその横方向へのプロファイルをそれぞれ図１９（ａ）、図１９（ｂ）に示し、水中において横方向に並べた２本のアルミ棒をＸ線撮影した場合の投影画像（測定画像を対数変換した画像）及びその横方向へのプロファイルをそれぞれ図２０（ａ）、図２０（ｂ）に示す。

アルミ棒が存在する領域について、図１９（ａ）を見ると測定画像ではノイズが小さいように見えるが、実際は図１９（ｂ）から分かるように信号対雑音比が悪い。そのため、投影画像ではアルミ棒が存在する領域周りよりもノイズが大きくなって現れる。

特開２００７−２００７６９号公報（段落００４４〜００５４） 特開２００８−００６２８８号公報 特表２０１０−５２８７６４号公報 特表２０１０−５２７７４１号公報

特許文献１では、電子顕微鏡画像に最大エントロピー法を用いた画像処理を施し、ポアソンノイズの低減を図っている。しかしながら、特許文献１で提案されている画像処理は、ポアソンノイズがビームハードニング補正処理によって非線形に増幅された場合に、充分なノイズ除去が期待できるものではない。

また、特許文献２〜４で提案されているコスト関数（エネルギー関数）を用いた画像処理は、画像処理の問題をエネルギー最小化問題に帰着させている点で、後述する本発明と共通しているものの、取り扱う問題（具体的な画像処理の内容）が本発明とは異なっている。

本発明は、上記の状況に鑑み、ポアソンノイズがビームハードニング補正処理によって非線形に増幅された場合であっても、画像のポアソンノイズを適切に除去することができるポアソンノイズ除去装置、ポアソンノイズ除去方法、及びＸ線撮影装置を提供することを目的とするものである。

上記目的を達成するために本発明に係るポアソンノイズ除去装置においては、入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去装置であって、前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を最小にするまたは近似的に最小にする前記出力画像を求める処理であり、前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子を含んでいる構成（第１の構成）としている。

このような構成によると、ポアソンノイズ除去処理で用いるエネルギー関数が、注目画素における平滑化の強さを示し注目画素毎に動的に変化する因子を含んでいるので、平滑化の強さを画素毎に動的に変化させることができる。これにより、ポアソンノイズがビームハードニング補正処理によって非線形に増幅された場合であっても、画像のポアソンノイズを適切に除去することができる。

さらに、前記エネルギー関数が、
又は
（ここで、
は前記出力画像であり、
は前記入力画像であり、iは前記注目画素のラベルであり、ijは前記注目画素の近傍に位置する画素のラベルであり、h(f_i)が前記注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子である。）
である構成（第２の構成）が好ましい。

このような構成によると、互いに近傍に位置する画素同士は同じ値を持ちやすいという傾向の強さをパラメータαによって容易に調整でき、ノイズの全体的な大きさをパラメータβによって容易に調整できる。また、上記第２の構成において、前記出力画像が、確率分布
を用いた期待値計算によって前記エネルギー関数を近似的に最小にする画像である構成（第３の構成）にしてもよい。さらに、上記第３の構成において、前記期待値計算が、前記確率分布に平均場近似を適用した計算である構成（第４の構成）にしてもよい。上記第４の構成によると、注目画素における平滑化の強さを示し注目画素毎に動的に変化する因子を含んでいるエネルギー関数の最小化問題を、確率分布を用い当該確率分布に平均場近似を適用した期待値計算によって解くことになるので、ガウス関数の公式等を利用することができ、解析的な扱いが可能となる。

また、上記目的を達成するために本発明に係るポアソンノイズ除去方法においては、入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去方法であって、前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を最小にするまたは近似的に最小にする前記出力画像を求める処理であり、前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子を含んでいるようにする。

また、上記目的を達成するために本発明に係るＸ線撮影装置においては、被写体に対してＸ線を照射するＸ線照射部と、前記被写体を透過したＸ線を検出するＸ線検出部と、前記Ｘ線検出部の検出結果に基づく画像を入力する上記いずれかの構成のポアソンノイズ除去装置とを備える構成としている。

さらに、Ｘ線撮影で得られる投影画像に対して、第１の関数y=ｆ(x)を用いてビームハードニング補正を行うビームハードニング補正装置を備え、第１の関数が、第２の関数の逆関数であり、第２の関数が、ロジスティック曲線を表す第３の関数y=γ/[1+exp[-δ(x-ε)]]を平行移動して得られる関数であることが好ましい。

ビームハードニング補正を行われた場合、ノイズが増強されることになるので、当該構成によると、ポアソンノイズ除去装置によるポアソンノイズ除去がより一層効果的になる。

さらに、第１の関数が、
であり、
前記ポアソンノイズ除去装置が上記第２の構成のポアソンノイズ除去装置であって、
（ここで、I₀(i)は、被写体が存在しないときの測定画像のラベルがiである画素の画素値である。）
を満たすことが好ましい。

当該構成によると、ビームハードニング補正によるノイズの増強に対応した因子h(f_i)が設定されるので、ポアソンノイズ除去装置によるポアソンノイズ除去がより一層適切になる。

本発明によると、ポアソンノイズ除去処理で用いるエネルギー関数が、注目画素における平滑化の強さを示し注目画素毎に動的に変化する因子を含んでいるので、平滑化の強さを画素毎に動的に変化させることができる。これにより、ポアソンノイズがビームハードニング補正処理によって非線形に増幅された場合であっても、画像のポアソンノイズを適切に除去することができる。

本発明の一実施形態に係るＸ線撮影装置の本体部の外観を示す図である。 パノラマ撮影モードの標準軌道を示す図である。 局所ＣＴ撮影モードの軌道を示す図である。 局所ＣＴ撮影モードにおいて、撮影対象部位の中心を前歯の位置に設定した場合の撮影対象部位の中心及び画像再構成範囲を示す図である。 局所ＣＴ撮影モードにおいて、撮影対象部位の中心を左顎の位置に設定した場合の撮影対象部位の中心及び画像再構成範囲を示す図である。 局所ＣＴ撮影モードにおいて、撮影対象部位の中心を左第２小臼歯の位置に設定した場合の撮影対象部位の中心及び画像再構成範囲を示す図である。 全歯ＣＴ撮影モードの軌道を示す図である。 全顎ＣＴ撮影モードの軌道を示す図である。 画像処理装置の構成を示す図である。 測定画像及び投影画像の各例を示す図である。 第３の関数の一例を示すグラフである。 第２の関数の一例を示すグラフである。 第１の関数の一例を示すグラフである。 頭部ファントムのビームハードニング補正を適用して得られた再構成画像及びその２値化画像を示す図である。 頭部ファントムの再構成画像及びその２値化画像を示す図である。 水中において横方向に並べた２本のアルミ棒をＸ線撮影した場合のビームハードニング補正を適用して得られた投影画像及びその横方向へのプロファイルを示す図である。 比較例のノイズ除去処理を行い、ビームハードニング補正を適用して得られる再構成画像及びそれと図１４（ａ）の再構成画像との差分画像を示す図である。 本発明に係るポアソンノイズ除去処理を行い、ビームハードニング補正を適用して得られる再構成画像及びそれと図１４（ａ）の再構成画像との差分画像を示す図である。 水中において横方向に並べた２本のアルミ棒をＸ線撮影した場合の測定画像及びその横方向へのプロファイルを示す図である。 水中において横方向に並べた２本のアルミ棒をＸ線撮影した場合の投影画像及びその横方向へのプロファイルを示す図である。

本発明の実施形態について図面を参照して以下に説明する。

まず始めに、本発明の一実施形態に係るＸ線撮影装置の本体部１（以下、「Ｘ線撮影装置の本体部１」と称す）の構成について図１を参照して説明する。図１はＸ線撮影装置１の本体部の外観を示す図であり、図１（ａ）は上面図、図１（ｂ）は正面図、図１（ｃ）は側面図である。

Ｘ線撮影装置の本体部１は、歯科用あるいは耳鼻科用等のＸ線撮影装置の本体部であって、床面に載置されるベース２と、ベース２から鉛直方向に立設された下部ポール３と、鉛直方向にスライド可能に下部ポール３に接続される上部ポール４と、上部ポール４の上端部に固定されている固定アーム５と、回転可能に固定アーム５に接続される旋回アーム６と、上部ポール４の中央部に固定されており被写体（例えば歯など）を含む人体の頭部を保持する頭部保持部７とを備えている。実施形態では、固定アーム５が上部ポール４に固定されているが、例えば、Ｘ線撮影装置の本体部１を設置する部屋の壁や天井に固定アーム５が直接あるいは部屋の壁や天井との距離を調整することができる調整機構を介して取り付けられる態様であってもよい。

旋回アーム６は、被写体に対してＸ線を照射するＸ線照射部８と、被写体を透過したＸ線を検出するＸ線検出部９とを対向して配置している。本実施形態では、Ｘ線検出部９として、照射されたＸ線に応じて電気信号を生成する変換素子が二次元状に配置されている二次元Ｘ線検出器を用いる。

Ｘ線撮影装置の本体部１の撮影モードは特に限定されないが、例えば、パノラマ撮影モードやＣＴ撮影モードを挙げることができる。パノラマ撮影モードでは、Ｘ線照射部８及びＸ線検出部９が歯列弓の形状に沿った所定の軌跡を描くように、旋回アーム６の旋回軸を旋回軸に垂直な方向（Ｘ方向、Ｙ方向）に移動させ、旋回アーム６を旋回軸回りに旋回させながら断層撮影を行う。ＣＴ撮影モードでは、頭部の対象撮影領域（画像再構成範囲）を中心にして旋回アーム６を回転させながら、対象撮影領域（画像再構成範囲）の断層撮影を行う。

ここで、パノラマ撮影モードについて図２を参照してより詳細に説明する。図２はパノラマ撮影モードの標準（成人用）軌道を示している。パノラマ撮影モードの標準（成人用）軌道では、Ｘ線照射部８及びＸ線検出部９が仮想歯列弓２０１の形状に沿った所定の軌跡を描いてＸ線ビームの軌跡が包絡線状の軌跡２０２になるように、Ｘ線照射部８及びＸ線検出部９が配置されている旋回アーム６を、撮影開始位置Ｐ１から図２に示す軌道に沿って撮影終了位置Ｐ２まで移動させる。撮影開始位置Ｐ１と撮影終了位置Ｐ２との間における旋回アーム６の旋回角度は約２２０度である。なお、撮影終了位置Ｐ２を除く図２に示された旋回アーム６の位置は被写体の撮影領域の左半分における各撮影位置である。Ｘ線照射部８のＸ線焦点８Ａから射出されるＸ線は、Ｘ線照射部８に設けられているＸ線絞り８Ｂによって絞られ、Ｘ線検出部９上でのＸ線ビーム幅Ｗが調整される。

パノラマ撮影モードは、上述した標準（成人用）軌道の他に、小児用軌道、直行軌道、顎関節撮影軌道、上顎洞撮影軌道などを有していることが好ましい。小児用軌道は、仮想歯列弓２０１の形状が小さくなる点が標準軌道と異なっている。直行軌道は、各撮影位置でのＸ線ビームが患者歯列弓２０３の歯と歯の間を通過するようにしている点が標準軌道と異なっている。顎関節撮影軌道（側面）は、Ｘ線照射部８及びＸ線検出部９が仮想歯列弓２０１の両端部分（顎関節撮影可能部分）の形状に沿った所定の軌跡を描くように旋回アーム６を移動させる点が標準軌道と異なっている。顎関節撮影軌道（正面）は、Ｘ線照射部８及びＸ線検出部９が仮想線２０４の形状に沿った所定の軌跡を描くように旋回アーム６を移動させる点が標準軌道と異なっている。上顎洞撮影軌道は、Ｘ線照射部８及びＸ線検出部９が仮想線２０５の形状に沿った所定の軌跡を描くように旋回アーム６を移動させる点が標準軌道と異なっている。

続いて、ＣＴ撮影モードについて図３〜図８を参照してより詳細に説明する。なお、図３〜図８において図２と同一の部分には同一の符号を付す。

局所ＣＴ撮影モードは、歯顎領域内の上下歯牙領域全体よりも狭い特定の領域を撮影対象とするＣＴ撮影モードである。局所ＣＴ撮影モードの画像再構成範囲は例えば直径５１ｍｍ高さ５５ｍｍの円柱形状の空間領域である。図３は局所ＣＴ撮影モードの軌道を示している。局所ＣＴ撮影モードでは、図３に示すように、Ｘ線検出部９の中心がＸ線照射部８と旋回アーム６の旋回軸中心２０６とを結ぶラインの延長線上にくるように旋回アーム６を旋回させながら複数の撮影位置で撮影が行われる。また、局所ＣＴ撮影モードでは、通常、図３に示すように、旋回アーム６の旋回軸中心２０６は定位置になっている。なお、図３には撮影位置として４箇所が図示されているが、これはあくまで例示であり撮影位置は図示された箇所に限定されるものではない。

局所ＣＴ撮影モードは、後述する全歯ＣＴ撮影モードや全顎ＣＴ撮影モードに比べてＸ線検出部９上でのＸ線ビーム幅Ｗが狭いため、Ｘ線検出部９のサイズが小さくても実施可能である。

なお、局所ＣＴ撮影モードでは、撮影対象部位（関心領域）の中心を何処に設定するかに応じて旋回アーム６の旋回軸中心２０６の位置を変えるようにしており、通常、図３に示すように、撮影対象部位（関心領域）の中心と旋回アーム６の旋回軸中心２０６の位置とが一致するように位置調整がなされる。局所ＣＴ撮影モードにおける撮影対象部位（関心領域）の中心は任意に設定することができる。図３に示した位置設定の他にも、例えば、図４に示すように撮影対象部位（関心領域）の中心２０８を仮想歯列弓２０１上の前歯の位置に設定することもでき、図５に示すように撮影対象部位（関心領域）の中心２０８を仮想歯列弓２０１上の左顎の位置に設定することもでき、図６に示すように撮影対象部位（関心領域）の中心２０８を仮想歯列弓２０１上の右第２小臼歯の位置に設定することもでき、その他種々の位置設定が可能である。

全歯ＣＴ撮影モードは、上下歯牙領域全体を撮影対象とするＣＴ撮影モードである。全歯ＣＴ撮影モードの画像再構成範囲は例えば直径９７ｍｍ高さ１００ｍｍの円柱形状の空間領域である。図７は全歯ＣＴ撮影モードの軌道を示している。全歯ＣＴ撮影モードでは、図７に示すように、Ｘ線検出部９の中心がＸ線照射部８と旋回アーム６の旋回軸中心２０６とを結ぶラインの延長線上にくるように旋回アーム６を旋回させながら複数の撮影位置で撮影が行われる。また、全歯ＣＴ撮影モードでは、通常、図７に示すように、旋回アーム６の旋回軸中心２０６は定位置になっている。なお、図７には撮影位置として４箇所が図示されているが、これはあくまで例示であり撮影位置は図示された箇所に限定されるものではない。

全歯ＣＴ撮影モードは、上述した局所ＣＴ撮影モードに比べて撮影対象が広範囲になりＸ線検出部９上でのＸ線ビーム幅Ｗが広くなるため、その広いＸ線ビーム幅Ｗに見合ったＸ線検出部９のサイズを必要とする。

全顎ＣＴ撮影モードは、歯顎領域の全ての範囲を撮影対象とするＣＴ撮影モードである。全顎ＣＴ撮影モードの画像再構成範囲は例えば直径１６１ｍｍ高さ１００ｍｍの円柱形状の空間領域である。図８は全顎ＣＴ撮影モードの軌道を示している。全顎ＣＴ撮影モードでは、図８に示すように、Ｘ線検出部９の中心がＸ線照射部８と旋回アーム６の旋回軸中心２０６とを結ぶラインの延長線上からずれるように旋回アーム６を旋回させながら複数の撮影位置で撮影が行われる。また、全顎ＣＴ撮影モードでは、通常、図８に示すように、旋回アーム６の旋回軸中心２０６は定位置になっている。なお、図８には撮影位置として４箇所が図示されているが、これはあくまで例示であり撮影位置は図示された箇所に限定されるものではない。

全顎ＣＴ撮影モードは、Ｘ線検出部９の中心をＸ線照射部８と旋回アーム６の旋回軸中心２０６とを結ぶラインの延長線上からずらして撮影を行っているので、上述した全歯ＣＴ撮影モードよりも画像再構成範囲２０７を拡大することができる。したがって、Ｘ線検出部９のサイズアップを抑えながら歯顎領域の全ての範囲を撮影対象とすることができる。

なお、全顎ＣＴ撮影モードにおいて、Ｘ線検出部９をサイズアップして、Ｘ線検出部９上でのＸ線ビーム幅Ｗを図８に示す場合よりも拡大し、画像再構成範囲を例えば直径２３０ｍｍ高さ１６４ｍｍの円柱形状の空間領域にすることで、歯顎領域の全ての範囲のみならず、頭頸部領域の全ての範囲を撮影対象とすることも可能である。

上述したパノラマ撮影モード及びＣＴ撮影モードでは、撮影時に患者歯列弓２０３が想定した位置（図２、図３、図７、図８に図示した位置）に存在することで、撮影者が意図していた通りの撮影を行うことができる。患者歯列弓２０３の想定した位置への位置合わせを容易に実現する方法としては、例えば、光ビームを利用する方法を挙げることができる。当該光ビームとしては、例えば、頭の正中線の位置を示す正中線光ビーム、眼窩下縁と外耳道を結ぶ線の位置を示す水平線光ビーム、犬歯の位置（断層撮影の基準位置）を示す断層基準線光ビームなどがあり、これらの光ビームの出力部をＸ線撮影装置に設け、これらの光ビームを参考にして患者が頭の位置を微調整するとよい。

また、旋回アーム６から離れた位置に設置するセファロ用ユニット（不図示）を用い、セファロ撮影モードでの撮影が行えるようにしてもよい。セファロ用ユニットは、被写体を透過 したＸ線を検出して、被写体をセファロ撮影するためのセファロ用Ｘ線検出部と、頭部を固定するための頭部固定部とを備える。セファロ撮影は、歯科矯正の診断等に用いられ、頭部規格Ｘ線撮影法（セファロ撮影法）を用いて撮影する。セファロ撮影では、例えば、頭部固定部のイヤーロッドを頭部の左右の外耳孔部に挿入して固定し、旋回アーム６に設けられたＸ線照射部８からＸ線を照射して、被写体を透過したＸ線をセファロ用Ｘ線検出部で検出する。

本発明の一実施形態に係るＸ線撮影装置は、Ｘ線撮影装置の本体部１の他に、図９に示す画像処理装置１０も備えている。

画像処理装置１０は、ＲＯＭ１０２やＨＤＤ１０７に格納されているプログラムに従って画像処理装置１０全体を制御するＣＰＵ１０１と、固定的なプログラムやデータを記録するＲＯＭ１０２と、作業メモリを提供するＲＡＭ１０３と、Ｘ線撮影装置の本体部１内に格納されＸ線撮影装置の本体部１の各部を制御する制御部（不図示）との間で通信を行うための通信インターフェース部１０４と、画像データを一時的に記憶するＶＲＡＭ１０５と、ＶＲＡＭ１０５に記憶された画像データに基づいて画像を表示する表示部１０６と、前記制御部及びＣＰＵ１０１が協働してＸ線撮影動作を制御するための撮影制御プログラム、再構成画像を生成するための画像再構成処理プログラム、ポアソンノイズ除去処理を行うためのポアソンノイズ除去処理プログラム等の各種プログラム、各種プログラムを実行する際に用いられる各種パラメータの設定値、並びに、再構成画像データ等の各種データを記憶するＨＤＤ１０７と、キーボード、ポインティングデバイス等の入力部１０８とを備えている。

画像処理装置１０と前記制御部との通信方法は、有線通信でもよく、無線通信でもよく、有線と無線を組み合わせた通信であってもよい。画像処理装置１０としては、例えば、パーソナルコンピュータを挙げることができる。なお、画像処理装置１０は、画像処理以外に、Ｘ線撮影装置の本体部１の遠隔操作、画像表示も行う。ＨＤＤ１０７に記憶されている各プログラムは、画像処理装置１０にプリインストールされていてもよく、光ディスク等の記憶媒体に格納された形態で流通されて画像処理装置１０にインストールされてもよく、ネットワークを介して流通されて画像処理装置１０にインストールされてもよい。

ポアソンノイズ除去処理プログラムを実行すると、画像処理装置１０はポアソンノイズ除去装置として機能する。ポアソンノイズ除去処理は画像再構成処理中に割り込んで実施される。

Ｘ線検出部９から出力され、画像処理装置１０が受信する画像は、測定画像と呼ばれる。測定画像の各画素値は、Ｘ線検出部９の各検出素子に到達したＸ線量子の個数に比例した値を表しているとみなすことができる。

画像再構成処理プログラムで使用される再構成アルゴリズムは、Ｘ線量子の個数ではなく、Ｘ線の線減弱係数の線積分を用いて断面を再構成する。このため、画像再構成処理において、Ｘ線量子の個数の分布を示す測定画像を、線減弱係数の線積分の分布を示す画像（投影画像）に変換する必要がある。以下、測定画像から投影画像への変換（対数変換）について詳述する。

画素のラベルをi、被写体が存在しないときの測定画像の画素値をI₀(i)、今回再構成を所望する被写体を置いたときの測定画像の画素値をI(i)、Ｘ線の経路をs、位置sにおける線減弱係数をμ(s)とすると、次の式（１）が成り立つ。この式（１）は、Ｘ線の減弱の様子を表す方程式である。

式（１）の左辺が測定画像の画素値であるのに対し、右辺の指数関数の中身は投影画像の画素値を正負反転したものである。よって、式（１）を指数関数の中身に関して解くことで、測定画像から投影画像への変換（対数変換）の変換式が求まる。よって、対数変換の変換式は次の式（２）のようになる。

ここで、測定画像の一例を図１０（ａ）に示し、図１０（ａ）の測定画像を対数変換して得られる投影画像を図１０（ｂ）に示す。

画像処理装置１０は、測定画像または投影画像をポアソンノイズ除去装置の入力画像として扱う。当該ポアソンノイズ除去装置が入力画像に施すポアソンノイズ除去処理は、当該ポアソンノイズ除去装置の出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を近似的に最小にする前記出力画像を求める処理である。そして、前記エネルギー関数は、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子を含んでいることを特徴としている。

本実施形態では、前記エネルギー関数を次の式（３）としている。

ここで、
は前記出力画像であり、
は前記入力画像であり、iは前記注目画素のラベルであり、ijは前記注目画素の近傍に位置する画素のラベルであり、h(f_i)が前記注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子である。

式（３）で表されるエネルギー関数の第１項は、互いに近傍に位置する画素同士は同じ値を持ちやすいという仮定に基づいて設計されている項であり、パラメータαは互いに近傍に位置する画素同士は同じ値を持ちやすいという傾向の強さを表すパラメータである。また、式（３）で表されるエネルギー関数の第２項は，因子h(f_i)に従って平滑化の強さが変化するという仮定に基づいて設計されている項であり、パラメータβはノイズの全体的な大きさを表すパラメータである。

式（３）で表されるエネルギー関数は、因子h(f_i)の値が小さいほど平滑化が強くなり、因子h(f_i)の値が大きいほど平滑化が弱くなるように設計されている。因子h(f_i)は注目画素毎に動的に変化する因子であるため、式（３）で表されるエネルギー関数を最小にする出力画像を求めようとすると、必然的に画像内の場所（画素）に応じて平滑化の強さが調節されることになり、画像のポアソンノイズを適切に除去することができるノイズ除去アルゴリズムが構築できる。

因子h(f_i)の具体的な形としては様々なものが考えられるが、図１９（ｂ）及び図２０（ｂ）から分かるように、少なくとも次の方針[１]、[２]には従うことになる。
[１] 入力画像が測定画像である場合（測定画像を対数変換する前にポアソンノイズ除去処理を実施する場合）は、注目画素の近傍に位置する画素の画素値が小さいほど、注目画素の因子h(f_i)が小さくなるようにする。
[２] 入力画像が投影画像である場合（測定画像を対数変換して投影画像を得た後にポアソンノイズ除去処理を実施する場合）は、注目画素の近傍に位置する画素の画素値が小さいほど、注目画素の因子h(f_i)が大きくなるようにする。

因子h(f_i)の具体例としては、下記のようなものを挙げることができる。因子h(f_i)の決定方法としては、例えば、複数の因子h(f_i)を試行錯誤して、その中から最も適切にポアソンノイズを除去するものを選定するようにすればよい。

入力画像が測定画像である場合、最も適切にポアソンノイズを除去すると考えられる因子h(f_i)を次に示す。画像処理装置１０がビームハードニング補正を行わない場合には定数aを大きな値にしてy_iとx_iとがほぼ同一になるようにすればよく、画像処理装置１０がビームハードニング補正処理を行う場合にはビームハードニング補正の強さに応じて定数aの値を調整するとよい。
ここで、I₀(i)は、被写体が存在しないときの測定画像のラベルがiである画素の画素値である。

次に、式（３）で表されるエネルギー関数を最小にする前記出力画像を近似的に求める手法について説明する。

まずは、式（３）で表されるエネルギー関数の符号を逆転させ、それを指数関数の中身に放り込み、さらにこれを規格化することで次の式（４）で表される確率分布が得られる。なお、式（３）で表されるエネルギー関数が低い値を取るほど、式（４）で表される確率分布は高い確率を示すことになる。

式（４）で表される確率分布は多変数の確率分布であり、因子h(f_i)の中に確率変数が含まれるため、解析的な扱いができない。そこで、解析的な扱いを可能にするために、式（４）に対して平均場近似を適用する。平均場近似とは、注目画素以外の全ての画素の値を期待値m_ijに固定し、次の式（５）で与えられる近似確率分布に近似することをいう。

平均場近似を適用する前に、注目画素の画素値が注目画素の近傍に位置する画素の画素値の平均に近い値になるという仮定の下、因子h(f_i)を次の式（６）のように近似する。
ここで、Jは注目画素の近傍に位置する画素の個数である。
上記の式（６）に対して平均場近似を適用すると、注目画素の近傍に位置するJ個の画素の画素値f_ijが全て期待値m_ijに固定されるため、因子h(f_i)が次の式（７）に示すように単なる定数となる。
したがって、式（４）に対して平均場近似を適用すると、次の式（８）で表される１変数の確率分布となる。

平均場近似を適用した式（８）は１変数の確率分布であり、因子
は単なる定数であるため、ガウス積分の公式などが適用でき、解析的な扱いができるようになる。よって、あとは各画素の期待値を次の式（９）で表される連立方程式に従って計算すればよい。この連立方程式は反復法を用いて解くことができる。全ての画素値の期待値が収束したとき、ノイズが除去された画像（エネルギー関数を最小にする出力画像）が得られる。すなわち、画像処理装置１０は、式（３）で表されるエネルギー関数を近似的に最小にする出力画像を求める処理として、式（９）で表される連立方程式を反復法を用いて解く処理を実行している。なお、上述した式（３）で表されるエネルギー関数が２つのパラメータα、βを導入しているのに対し、式（９）で表される連立方程式ではその２つのパラメータα、βを１つのパラメータAにまとめることができる。

なお、画像処理装置１０は、ビームハードニング補正処理を行うためのビームハードニング補正処理プログラムをＨＤＤ１０７に記憶し、ビームハードニング補正装置としても機能することが好ましい。これにより、ビームハードニングによる誤差に起因するアーチファクトの出現を抑えることができる。

１回の処理で再構成画像の軟組織撮影領域も硬組織撮影領域もビームハードニング補正が可能であるようにするために、ここでは、画像処理装置１０が、投影画像に対して、第１の関数y=ｆ(x)を用いてビームハードニング補正を行うようにする。ここで、第１の関数は、第２の関数の逆関数である。また、第２の関数は、ロジスティック曲線を表す第３の関数y=γ/[1+exp[-δ(x-ε)]]を平行移動して得られる関数である。

また、本実施形態では、第３の関数の各定数γ、δ、εをγ=2a、δ=2/a、ε=0に設定している。さらに、本実施形態では、第２の関数が、原点を通る関数であって、尚且つ、第２の関数の原点での接線の傾きが1であるという条件、第２の関数の変曲点が原点であるという条件を満たす関数となるように、第３の関数をy方向に-a平行移動して得られる関数を第２の関数にしている。よって、本実施形態では、第３の関数は次の式（１０）で表され、第２の関数は次の式（１１）で表され、第１の関数は次の式（１２）で表される。そして、式（１０）で定数aの設定値を5にしたときのグラフは図１１に示すようになり、式（１１）で定数aの設定値を5にしたときのグラフは図１２に示すようになり、式（１２）で定数aの設定値を5にしたときのグラフは図１３に示すようになる。図１３に示すグラフのx軸はビームハードニング補正前の投影画像の画素値を示す軸であり、図１３に示すグラフのy軸はビームハードニング補正後の投影画像の画素値を示す軸である。

投影画像の画素値は零または正の値であるので、本実施形態では第１の関数の原点または第１象限を用いてビームハードニング補正を行っている。上述した第１の関数を用いたビームハードニング補正によって、投影画像の画素値は大きいほど、より大きくなるように変換される。また、上述した第１の関数を用いたビームハードニング補正が軟組織撮影領域にも硬組織撮影領域にも適しているので、上述した第１の関数を用いたビームハードニング補正のみ、つまり１回の処理で再構成画像の軟組織撮影領域も硬組織撮影領域もビームハードニング補正が可能である。

ここで、Ｘ線撮影装置の被写体が頭部ファントムであり、上述したポアソンノイズ除去処理は行わず、ビームハードニング補正を適用して得られる再構成画像を図１４（ａ）に示し、図１４（ａ）を２値化した画像を図１４（ｂ）に示す。これに対して、上述したポアソンノイズ除去処理もビームハードニング補正も適用しないで得られる再構成画像を図１５（ａ）に示し、図１５（ａ）を２値化した画像を図１５（ｂ）に示す。２値化した画像である図１５（ｂ）を観察すると、前歯部と臼歯部とで画素値が異なることが分かる。本来は前歯部と臼歯部とでは画素値は同じはずであるが、ビームハードニングの影響で画素値が変動している。このことは、図１５（ａ）にビームハードニングによる誤差に起因するアーチファクトが出現していることを示している。一方、２値化した画像である図１４（ｂ）を観察すると、前歯部と臼歯部との大部分で画素値が同じであることが分かる。このことは、図１４（ａ）にビームハードニングによる誤差に起因するアーチファクトが殆ど出現していないことを示している。上述した図１４と図１５との比較結果から明らかな通り、ビームハードニング補正を適用することで、ビームハードニングによる誤差に起因するアーチファクトを低減できていることが分かる。

なお、定数aの設定値を変更することでビームハードニング補正の強さを調整することができるので、定数aの設定値が変更可能であることが好ましい。

ビームハードニング補正を適用することで、ビームハードニングによる誤差に起因するアーチファクトを低減できるものの、その反面ノイズが増強されてしまう。例えば、水中において横方向に並べた２本のアルミ棒をＸ線撮影した場合の上述したポアソンノイズ除去処理は行わず、ビームハードニング補正を適用して得られる投影画像及びその横方向へのプロファイルをそれぞれ図１６（ａ）、図１６（ｂ）に示す。図１６（ｂ）と図２０（ｂ）とを比較すれば明らかな通り、ビームハードニング補正を適用することで、ノイズが増強されてしまうことが分かる。

このようにビームハードニング補正を行われた場合、ノイズが増強されることになるので、本発明に係るポアソンノイズ除去がより一層効果的になる。ポアソンノイズ除去はノイズを早く除去しておく観点からビームハードニング補正の前に行うようにすれば良いが、ビームハードニング補正の後に行うことも可能である。そして、式（１２）で表される関数でビームハードニング補正が行われた場合、因子h(f_i)を次のように設定し、因子h(f_i)中の定数aの値と式（１２）で表される関数中の定数aの値とを一致させることが望ましい。
ここで、I₀(i)は、被写体が存在しないときの測定画像のラベルがiである画素の画素値である。x _iはf_iを対数変換したものであり、y _iはx_iを式（１２）で表される関数に代入して得られるものであり、h(f_i)はy _iを逆対数変換して得られるものであるため、式（１２）で表される関数を用いたビームハードニング補正で増強されるノイズを適切に除去することができる。

ここで、平滑化の強さが画像の場所に依らず一定であるノイズ除去処理（比較例のノイズ除去処理）を行い、ビームハードニング補正を適用して得られる再構成画像を図１７（ａ）に示し、図１７（ａ）の再構成画像と、ノイズ除去処理は行わず、ビームハードニング補正を適用して得られる図１４（ａ）の再構成画像との差分画像を図１７（ｂ）に示す。また、本発明に係るポアソンノイズ除去処理を行い、ビームハードニング補正を適用して得られる再構成画像を図１８（ａ）に示し、図１８（ａ）の再構成画像と、ノイズ除去処理は行わず、ビームハードニング補正を適用して得られる図１４（ａ）の再構成画像との差分画像を図１８（ｂ）に示す。

差分画像において、歯列の輪郭が写っていれば、それだけ歯列がぼやけたことを意味し、歯列の輪郭以外で像が写っていれば、それだけノイズが除去できたことを意味している。図１７（ｂ）と図１８（ｂ）とを比較すれば明らかな通り、平滑化の強さが画像の場所に依らず一定であるノイズ除去処理（比較例のノイズ除去処理）よりも、本発明に係るポアソンノイズ除去処理の方が、歯列のエッジを保存し、かつノイズも除去できている。つまり、比較例のノイズ除去処理よりも、本発明に係るポアソンノイズ除去処理の方が、ポアソンノイズを適切に除去できている。

なお、上述した実施形態では、式（３）で表されるエネルギー関数を用いたが、他のエネルギー関数、例えば次の式で表されるエネルギー関数を用いてもよい。

１ 本発明の一実施形態に係るアーム型Ｘ線撮影装置の本体部
２ ベース
３ 下部ポール
４ 上部ポール
５ 固定アーム
６ 旋回アーム
７ 頭部保持部
８ Ｘ線照射部
８Ａ Ｘ線焦点
８Ｂ Ｘ線絞り
９ Ｘ線検出部
１０ 画像処理装置
１０１ ＣＰＵ
１０２ ＲＯＭ
１０３ ＲＡＭ
１０４ 通信インターフェース部
１０５ ＶＲＡＭ
１０６ 表示部
１０７ ＨＤＤ
１０８ 入力部
２０１ 仮想歯列弓
２０２ 包絡仮想歯列弓
２０３ 患者歯列弓
２０４、２０５ 仮想線
２０６ 旋回アームの旋回軸中心
２０７ 画像再構成範囲
Ｐ１ 撮影開始位置
Ｐ２ 撮影終了位置
Ｗ Ｘ線検出部上でのＸ線ビーム幅

Claims (13)

1. 入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去装置であって、
   前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を最小にするまたは近似的に最小にする前記出力画像を求める処理であり、
   前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する第１因子を含んでおり、
   前記第１因子が、増幅前のノイズの大きさに対応する平滑化の強さを入力とし、増幅後のノイズの大きさに対応する平滑化の強さを出力とする第２因子を含んでいることを特徴とするポアソンノイズ除去装置。
2. 入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去装置であって、
   前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を近似的に最小にする前記出力画像を求める処理であり、
   前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子を含んでおり、
   前記出力画像が、前記エネルギー関数を含んだ確率分布を用いた期待値計算によって前記エネルギー関数を近似的に最小にする画像であることを特徴とするポアソンノイズ除去装置。
3. 前記期待値計算が、前記確率分布に平均場近似を適用した計算である請求項２に記載のポアソンノイズ除去装置。
4. 入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去装置であって、
   前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を近似的に最小にする前記出力画像を求める処理であり、
   前記エネルギー関数が、
   
   又は
   
   （ここで、
   は前記出力画像であり、
   は前記入力画像であり、iは注目画素のラベルであり、ijは前記注目画素の近傍に位置する画素のラベルであり、h(f _i )が前記注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子である。）
   であることを特徴とするポアソンノイズ除去装置。
5. 前記出力画像が、確率分布
   
   を用いた期待値計算によって前記エネルギー関数を近似的に最小にする画像である請求項４に記載のポアソンノイズ除去装置。
6. 前記期待値計算が、前記確率分布に平均場近似を適用した計算である請求項５に記載のポアソンノイズ除去装置。
7. 入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去方法であって、
   前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を最小にするまたは近似的に最小にする前記出力画像を求める処理であり、
   前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する第１因子を含んでおり、
   前記第１因子が、増幅前のノイズの大きさに対応する平滑化の強さを入力とし、増幅後のノイズの大きさに対応する平滑化の強さを出力とする第２因子を含んでいることを特徴とするポアソンノイズ除去方法。
8. 入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去方法であって、
   前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を近似的に最小にする前記出力画像を求める処理であり、
   前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子を含んでおり、
   前記出力画像が、前記エネルギー関数を含んだ確率分布を用いた期待値計算によって前記エネルギー関数を近似的に最小にする画像であることを特徴とするポアソンノイズ除去方法。
9. 入力画像にポアソンノイズ除去処理を施して出力画像を生成するポアソンノイズ除去方法であって、
   前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を最小にするまたは近似的に最小にする前記出力画像を求める処理であり、
   前記エネルギー関数が、
   
   又は
   
   （ここで、
   は前記出力画像であり、
   は前記入力画像であり、iは注目画素のラベルであり、ijは前記注目画素の近傍に位置する画素のラベルであり、h(f _i )が前記注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子である。）
   であることを特徴とするポアソンノイズ除去方法。
10. 被写体に対してＸ線を照射するＸ線照射部と、
    前記被写体を透過したＸ線を検出するＸ線検出部と、
    前記Ｘ線検出部の検出結果に基づく画像を入力する請求項１〜６のいずれか一項に記載のポアソンノイズ除去装置とを備えることを特徴とするＸ線撮影装置。
11. Ｘ線撮影で得られる投影画像に対して、第１の関数y=ｆ(x)を用いてビームハードニング補正を行うビームハードニング補正装置を備え、
    第１の関数が、第２の関数の逆関数であり、
    第２の関数が、ロジスティック曲線を表す第３の関数y=γ/[1+exp{-δ(x-ε)}]を平行移動して得られる関数である請求項１０に記載のＸ線撮影装置。
12. 被写体に対してＸ線を照射するＸ線照射部と、
    前記被写体を透過したＸ線を検出するＸ線検出部と、
    前記Ｘ線検出部の検出結果に基づく画像を入力するポアソンノイズ除去装置と、
    Ｘ線撮影で得られる投影画像に対して、第１の関数y=ｆ(x)を用いてビームハードニング補正を行うビームハードニング補正装置とを備えるＸ線撮影装置であって、
    前記ポアソンノイズ除去装置は、前記Ｘ線検出部の検出結果に基づく画像にポアソンノイズ除去処理を施して出力画像を生成する装置であり、
    前記ポアソンノイズ除去処理が、前記出力画像がノイズが存在しない理想画像に近づくほど、その値が小さくなるように定式化されたエネルギー関数を用い、前記エネルギー関数を最小にするまたは近似的に最小にする前記出力画像を求める処理であり、
    前記エネルギー関数が、注目画素における平滑化の強さを示し前記注目画素毎に動的に変化する因子を含んでいる関数であり、
    第１の関数が、第２の関数の逆関数であり、
    第２の関数が、ロジスティック曲線を表す第３の関数y=γ/[1+exp{-δ(x-ε)}]を平行移動して得られる関数であることを特徴とするＸ線撮影装置。
13. 第１の関数が、
    
    であり、
    前記ポアソンノイズ除去装置が請求項４に記載のポアソンノイズ除去装置であって、
    
    （ここで、I ₀ (i)は、被写体が存在しないときの測定画像のラベルがiである画素の画素値である。）
    を満たす請求項１１に記載のＸ線撮影装置。