(57)【要約】【課題】撮像モダリティが異なる画像間での位置合せにおいて、ノイズによる勾配の影響を効果的に抑制する。【解決手段】第１及び第２の画像T,R上の座標pごとに、画素値の勾配に係る特徴ベクトルnを決定する第１の工程であって、座標pにおける画素値の勾配ベクトル∇T(p),∇R’(p)の大きさが閾値THNorm以上であれば、当該勾配ベクトルをその大きさで規格化したものを座標pの特徴ベクトルnとし、勾配ベクトルの大きさが閾値THNorm未満であれば、０ベクトルを座標pの特徴ベクトルnとする第１の工程と、第１及び第２の画像T,R’における互いに対応する座標pごとに、座標pの特徴ベクトル同士の内積の絶対値をN（Nは自然数）乗して成る相関値を算出し、これら相関値の積算値を含む評価値Dを求める第２の工程と、評価値Dが増大するよう第２の画像R’を変更する第３の工程と、を複数回繰返し行う。【選択図】図１
2014225218.tif 000002
【特許請求の範囲】
【請求項１】
同一の被写体を含む第１および第２の画像における座標ごとに、該座標における画素値の勾配に係る特徴ベクトルを決定する第１の工程であって、該座標における画素値の勾配ベクトルの大きさが所定の閾値以上であるときには、前記勾配ベクトルを該勾配ベクトルの大きさで規格化したものを該座標における特徴ベクトルとして決定し、前記勾配ベクトルの大きさが前記所定の閾値未満であるときには、０（ゼロ）ベクトルを該座標における特徴ベクトルとして決定する第１の工程と、
前記第１および第２の画像における互いに対応する座標ごとに、該座標における特徴ベクトル同士の内積の絶対値をＮ（Ｎは自然数）乗して成る値に相当する相関値を算出し、前記座標ごとに算出された相関値の積算値を含む評価値を求める第２の工程と、
前記評価値がより大きくなるように前記第２の画像を変更する第３の工程と、を複数回繰り返し行って、前記第１および第２の画像の位置合せを行う画像処理方法。
【請求項２】
同一の被写体を含む第１および第２の画像における座標ごとに、該座標における画素値の勾配ベクトルの大きさが所定の閾値以上であるときには、前記勾配ベクトルを該勾配ベクトルの大きさで規格化したものを該座標における特徴ベクトルとして決定し、前記勾配ベクトルの大きさが前記所定の閾値未満であるときには、０（ゼロ）ベクトルを該座標における特徴ベクトルとして決定する決定手段と、
前記第１および第２の画像における互いに対応する座標ごとに、該座標における特徴ベクトル同士の内積の絶対値をＮ（Ｎは自然数）乗して成る値に相当する相関値を算出し、前記座標ごとに算出された相関値の積算値を含む評価値を求める算出手段と、
前記評価値がより大きくなるように前記第２の画像を変更する変更手段と、を有し、
前記決定手段、算出手段および変更手段による処理を複数回繰り返し行って、前記第１および第２の画像の位置合せを行う画像処理装置。
【請求項３】
前記勾配ベクトルは、画像における各座標軸方向の１次偏微分である、請求項２に記載の画像処理装置。
【請求項４】
前記相関値は、前記特徴ベクトル同士の内積を二乗して成る値である、請求項３に記載の画像処理装置。
【請求項５】
前記評価値は、前記相関値の積算値を前記第１または第２の画像の大きさで規格化して成る値である、請求項２から請求項４のいずれか一項に記載の画像処理装置。
【請求項６】
前記変更は、平行移動、回転、および変形のうち少なくとも一つを含む、請求項２から請求項５のいずれか一項に記載の画像処理装置。
【請求項７】
前記決定手段は、前記第１および第２の画像のそれぞれに対して、構造物強調フィルタ処理および／またはノイズ低減処理を行った画像に基づいて、前記勾配ベクトルを求める、請求項２から請求項６のいずれか一項に記載の画像処理装置。
【請求項８】
前記第１および第２の画像は、医用画像である、請求項２から請求項７のいずれか一項に記載の画像処理装置。
【請求項９】
前記第１および第２の画像は、撮像モダリティが互いに異なる、請求項８に記載の画像処理装置。
【請求項１０】
前記第１および第２の画像の一方は、超音波画像（US画像）である、請求項９に記載の画像処理装置。
【請求項１１】
前記決定手段、算出手段および変更手段による処理を、該処理を行った回数が所定数に達するまで、または、前記評価値が実質的に収束するまで、繰り返し行う、請求項２から請求項１０のいずれか一項に記載の画像処理装置。
【請求項１２】
前記第１および第２の画像の少なくとも一方の対象領域におけるノイズ量を推定し、該ノイズ量に基づいて前記閾値を決定する手段をさらに有する、請求項２から請求項１１のいずれか一項に記載の画像処理装置。
【請求項１３】
前記対象領域は、前記第１の画像と前記第２の画像との間に共通する解剖学的特徴点の周辺領域である、請求項１２に記載の画像処理装置。
【請求項１４】
前記対象領域は、前記少なくとも一方の画像の中央領域である、請求項１２に記載の画像処理装置。
【請求項１５】
前記対象領域は、病変部の周辺領域である、請求項１２に記載の画像処理装置。
【請求項１６】
前記対象領域は、全画像領域である、請求項１２に記載の画像処理装置。
【請求項１７】
前記閾値を決定する手段は、前記対象領域における各画素での勾配強度を求め、該勾配強度のヒストグラムに基づいて前記閾値を決定する、請求項１２に記載の画像処理装置。
【請求項１８】
前記閾値を決定する手段は、ノイズ量Nqtを、以下の（ａ）〜（ｆ）式のいずれか１つ、または、複数の組合せの平均値から算出する、請求項１２から請求項１７のいずれか一項に記載の画像処理装置。
Nqt = Mfmax + σ×HWHML                     （ａ）
Nqt = Mfmax + σ×HWHMR                     （ｂ）
Nqt = Mfmax + σ×(HWHML+HWHMR)/2   （ｃ）
Nqt = σ×Mfmax                      （ｄ）
Nqt = σ×Mmom1                      （ｅ）
Nqt = σ×Mmom2                      （ｆ）
ただし、
Mfmaxは、前記ヒストグラム分布上に現れるピークのうち、最も勾配強度が低いピークにおいて最頻値を与える勾配強度であり、
HWHMLは、前記ヒストグラム分布上におけるMfmaxから見て低値側の半値半幅であり、
HWHMRは、前記ヒストグラム分布上におけるMfmaxから見て高値側の半値半幅であり、
HWHMmom1は、前記ヒストグラム分布上における勾配強度が0からHWHMRまでの範囲の重心に相当する勾配強度であり、
HWHMmom2は、前記ヒストグラム分布上における勾配強度がHWHMLからHWHMRまでの範囲の重心に相当する勾配強度である。
【請求項１９】
コンピュータを、請求項２から請求項１８のいずれか一項に記載の画像処理装置として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、撮像モダリティ（modality）が互いに異なる複数の画像間での位置合せ技術に関する。
【背景技術】
【０００２】
従来、撮像装置として、例えば、磁気共鳴装置（MR）、放射線断層撮影装置（CT）、超音波装置（US）などが知られている。これらの撮像装置には、その撮像モダリティごとにそれぞれ長所/短所が存在し、ある特定の撮像モダリティによる画像だけでは診断における精度が不足する場合がある。このため、近年、特定の撮像モダリティによる画像だけでなく、複数の異なる撮像モダリティによる画像を用いて診断することで、診断精度の向上を図る試みが盛んに行われている。
【０００３】
複数の異なる撮像モダリティによる画像を用いた診断では、撮像モダリティごとに画像の座標系が異なる。そのため、これらの座標系の相違や臓器の変動・変形に起因する位置ズレを補正する技術、すなわち、画像間での位置合せ（Registration）技術が重要である。
【０００４】
撮像モダリティが互いに異なる複数の画像間での位置合せ手法としては、相互情報量（Mutual Information）を用いる手法が最も一般的である（例えば、非特許文献１等参照）。この手法は、広義においては、画像の輝度値に基づく手法（intensity based method）である。つまり、相互情報量を用いて位置合せを行うには、対象画像間で輝度値に関連性があることが前提となる。
【０００５】
しかしながら、US画像においては、音響陰影（Acoustic shadow）が発生し、高反射体の後方の輝度値は本来の値よりも低下する。また血管の輝度値も血管の走行方向に依って変化する。このため、例えば、MR画像とUS画像との位置合せ、あるいは、CT画像とUS画像との位置合せにおいては、輝度値の関連性に乏しい状況がしばしば発生する。このような画像に対しては、相互情報量のような輝度値に基づく手法ではなく、画像の特徴に基づく手法（feature based method）の方が適切である。
【０００６】
画像の特徴に基づく手法の代表例として、標準勾配場（Normalized Gradient Field：NGF）を用いる手法が提案されている（例えば、非特許文献２等参照）。標準勾配場とは、画像上の座標において各方向ｘ，ｙ，ｚの１次偏微分、すなわち勾配ベクトル（Gradient Vector）を算出した後、その勾配ベクトルをその勾配ベクトルの長さ（Vector
Norm）で規格化（normalized）したものである。つまり、標準勾配場は、画素値あるいは輝度値の大小や勾配の大きさに依存せず、勾配の方向だけを表す特徴量である。仮に、ある２つの画像において互いに対応する位置に同じ方向の標準勾配場が発生しているならば、これら２つの画像の位置は合っていると見なすことができる。したがって、この手法では、標準勾配場が示す方向の揃い具合を最適化することで、位置合せを行うことができる。
【０００７】
標準勾配場ｎは、数学的には以下のように記述される。
【０００８】
2014225218.tif 000003
(1)式において、∇Ｉ(ｐ)は、画像Ｉの座標ｐにおける勾配ベクトルを表している。また、||∇Ｉ（ｐ）||は、その勾配ベクトルの長さ(Norm)を表している。座標ｐは、ｘ，ｙ，ｚの各成分からなる。
【０００９】
(1)式から理解されるように、標準勾配場ｎは、それ単独では、構造物に起因する勾配とノイズに起因する勾配とを分別できない。そこで、ノイズ（noise）に起因する勾配の影響を抑制するために、非特許文献２において、次式のような手法が提案されている。
【００１０】
2014225218.tif 000004
ここで、εは画像のノイズ量に応じて任意に選択する定数、である。従来法では、(4)式に従ってノイズ項εを設定し、ノイズ項を考慮して勾配ベクトルを規格化している。
【００１１】
位置合せの指標となる評価関数Ｄは、次式（いずれか一方でよい）で定義される。
【００１２】
2014225218.tif 000005
(5),(6)式において、Ｔ,Ｒは、それぞれ位置合せの対象となる２つの画像を表している。(5)式は、２つの標準勾配場ｎ(Ｒ,ｐ)、ｎ(Ｔ,ｐ)の内積の二乗の積算値を算出している。(6)式は、２つの標準勾配場ｎ(Ｒ,ｐ)、ｎ(Ｔ,ｐ)の外積の二乗の積算値を算出している。
【先行技術文献】
【特許文献】
【００１３】
【非特許文献１】IEEE Trans. on Med. Imaging, 16:187-198, 1997
【非特許文献２】Proceeding of SPIE Vol.7261, 72610G-1, 2009
【発明の概要】
【発明が解決しようとする課題】
【００１４】
前述の標準勾配場を用いる画像の位置合せ手法では、ノイズに起因するベクトルは標準勾配場上で微小な値をとる。しかし、(5),(6)式のように評価関数はすべての座標における内積または外積の積算値をとるため、一つ一つは微小な値であってもその絶対数が多ければ評価関数に与える影響は無視できない。そのため、ノイズに起因する勾配(gradient)の影響を十分に排除できない場合がある。
【００１５】
このような事情により、同一の被写体を含む２画像間の位置合せにおいて、ノイズに起因する勾配(gradient)の影響をより効果的に抑制することが可能な位置合せ手法が望まれている。
【課題を解決するための手段】
【００１６】
第１の観点の発明は、
同一の被写体を含む第１および第２の画像における座標ごとに、該座標における画素値の勾配に係る特徴ベクトルを決定する第１の工程であって、該座標における画素値の勾配ベクトルの大きさが所定の閾値以上であるときには、前記勾配ベクトルを該勾配ベクトルの大きさで規格化したものを該座標における特徴ベクトルとして決定し、前記勾配ベクトルの大きさが前記所定の閾値未満であるときには、０（ゼロ）ベクトルを該座標における特徴ベクトルとして決定する第１の工程と、
前記第１および第２の画像における互いに対応する座標ごとに、該座標における特徴ベクトル同士の内積の絶対値をＮ（Ｎは自然数）乗して成る値に相当する相関値を算出し、前記座標ごとに算出された相関値の積算値を含む評価値を求める第２の工程と、
前記評価値がより大きくなるように前記第２の画像を変更する第３の工程と、を複数回繰り返し行って、前記第１および第２の画像の位置合せを行う画像処理方法を提供する。
【００１７】
第２の観点の発明は、
同一の被写体を含む第１および第２の画像における座標ごとに、該座標における画素値の勾配ベクトルの大きさが所定の閾値以上であるときには、前記勾配ベクトルを該勾配ベクトルの大きさで規格化したものを該座標における特徴ベクトルとして決定し、前記勾配ベクトルの大きさが前記所定の閾値未満であるときには、０（ゼロ）ベクトルを該座標における特徴ベクトルとして決定する決定手段と、
前記第１および第２の画像における互いに対応する座標ごとに、該座標における特徴ベクトル同士の内積の絶対値をＮ（Ｎは自然数）乗して成る値に相当する相関値を算出し、前記座標ごとに算出された相関値の積算値を含む評価値を求める算出手段と、
前記評価値がより大きくなるように前記第２の画像を変更する変更手段と、を有し、
前記決定手段、算出手段および変更手段による処理を複数回繰り返し行って、前記第１および第２の画像の位置合せを行う画像処理装置を提供する。
【００１８】
第３の観点の発明は、
前記勾配ベクトルが、画像における各座標軸方向の１次偏微分である、上記第２の観点の画像処理装置を提供する。
【００１９】
第４の観点の発明は、
前記相関値が、前記特徴ベクトル同士の内積を二乗して成る値である、上記第３の観点の画像処理装置を提供する。
【００２０】
第５の観点の発明は、
前記評価値が、前記相関値の積算値を前記第１または第２の画像の大きさで規格化して成る値である、上記第２の観点から第４の観点のいずれか一つの観点の画像処理装置を提供する。
【００２１】
第６の観点の発明は、
前記変更が、平行移動、回転、および変形のうち少なくとも一つを含む、上記第２の観点から第５の観点のいずれか一つの観点の画像処理装置を提供する。
【００２２】
第７の観点の発明は、
前記決定手段が、前記第１および第２の画像のそれぞれに対して、構造物強調フィルタ処理および／またはノイズ低減処理を行った画像に基づいて、前記勾配ベクトルを求める、上記第２の観点から第６の観点のいずれか一つの観点の画像処理装置を提供する。
【００２３】
第８の観点の発明は、
前記第１および第２の画像が、医用画像である、上記第２の観点から第７の観点のいずれか一つの観点の画像処理装置を提供する。
【００２４】
第９の観点の発明は、
前記第１および第２の画像が、撮像モダリティが互いに異なる、上記第８の観点の画像処理装置を提供する。
【００２５】
第１０の観点の発明は、
前記第１および第２の画像の一方が、超音波画像（US画像）である、上記第９の観点の画像処理装置を提供する。
【００２６】
第１１の観点の発明は、
前記決定手段、算出手段および変更手段による処理を、該処理を行った回数が所定数に達するまで、または、前記評価値が実質的に収束するまで、繰り返し行う、上記第２の観点から第１０の観点のいずれか一つの観点の画像処理装置を提供する。
【００２７】
第１２の観点の発明は、
前記第１および第２の画像の少なくとも一方の対象領域におけるノイズ量を推定し、該ノイズ量に基づいて前記閾値を決定する手段をさらに有する、上記第２の観点から第１１の観点のいずれか一つの観点の画像処理装置を提供する。
【００２８】
第１３の観点の発明は、
前記対象領域が、前記第１の画像と前記第２の画像との間に共通する解剖学的特徴点の周辺領域である、上記第１２の観点の画像処理装置を提供する。
【００２９】
第１４の観点の発明は、
前記対象領域が、前記少なくとも一方の画像の中央領域である、上記第１２の観点の画像処理装置を提供する。
【００３０】
第１５の観点の発明は、
前記対象領域が、病変部の周辺領域である、上記第１２の観点の画像処理装置を提供する。
【００３１】
第１６の観点の発明は、
前記対象領域が、全画像領域である、上記第１２の観点のの画像処理装置を提供する。
【００３２】
第１７の観点の発明は、
前記閾値を決定する手段が、前記対象領域における各画素での勾配強度を求め、該勾配強度のヒストグラムに基づいて前記閾値を決定する、上記第１２から第１６の観点のいずれか一つの観点の画像処理装置を提供する。
【００３３】
第１８の観点の発明は、
前記閾値を決定する手段が、ノイズ量Nqtを、以下の（ａ）〜（ｆ）式のいずれか１つ、または、複数の組合せの平均値から算出する、上記第１２の観点から第１６の観点のいずれか一つの観点の画像処理装置を提供する。
Nqt = Mfmax + σ×HWHML                     （ａ）
Nqt = Mfmax + σ×HWHMR                     （ｂ）
Nqt = Mfmax + σ×(HWHML+HWHMR)/2   （ｃ）
Nqt = σ×Mfmax                      （ｄ）
Nqt = σ×Mmom1                      （ｅ）
Nqt = σ×Mmom2                      （ｆ）
ただし、
Mfmaxは、前記ヒストグラム分布上に現れるピークのうち、最も勾配強度が低いピークにおいて最頻値を与える勾配強度であり、
HWHMLは、前記ヒストグラム分布上におけるMfmaxから見て低値側の半値半幅であり、
HWHMRは、前記ヒストグラム分布上におけるMfmaxから見て高値側の半値半幅であり、
HWHMmom1は、前記ヒストグラム分布上における勾配強度が0からHWHMRまでの範囲の重心に相当する勾配強度であり、
HWHMmom2は、前記ヒストグラム分布上における勾配強度がHWHMLからHWHMRまでの範囲の重心に相当する勾配強度である。
【００３４】
第１９の観点の発明は、
コンピュータ(computer)を、上記第２の観点から第１８の観点のいずれか一つの観点の画像処理装置として機能させるためのプログラム(program)を提供する。
【発明の効果】
【００３５】
上記観点の発明によれば、勾配ベクトルの大きさが所定の閾値未満のときは、特徴ベクトルは、０ベクトルになる。そのため、特徴ベクトルは、主要な構造物による成分のみで構成されることになる。画像間での位置合せの評価値は、座標ごとの特徴ベクトル同士の内積の絶対値をＮ乗して成る相関値の積算値を含むので、ノイズに起因する勾配ベクトルの影響が効果的に抑制される。
【図面の簡単な説明】
【００３６】
【図１】本実施形態に係る画像処理装置の構成を概略的に示す図である。
【図２】位置合せ演算の概念図である。
【図３】従来法による標準勾配場の例を示す図である。
【図４】本提案法による閾値判定付き標準勾配場の例を示す図である。
【図５】本実施形態に係る画像処理装置による画像位置合せ処理のフロー図である。
【図６】従来法と本提案法とによる画像位置合せの例を示す図である。
【図７】目標画像Ｔ及び対象画像Ｒ′の少なくとも一方に基づいて画像のノイズ量を推定する処理のフロー図である。
【図８】勾配強度画像における勾配強度のヒストグラムの例を示す図である。
【図９】勾配強度のヒストグラムにおいて解析する特徴量を示す第１の図である。
【図１０】勾配強度のヒストグラムにおいて解析する特徴量を示す第２の図である。
【図１１】従来法による標準勾配場の例と本提案法（閾値THNormをノイズ量の自動推定結果から決定）による閾値判定付き標準勾配場の例とを示す。
【図１２】従来法と本提案法（閾値THNormをノイズ量の自動推定結果から決定）とによる画像位置合せの他の例を示す。
【発明を実施するための形態】
【００３７】
以下、発明の実施形態について説明する。なお、これにより、発明が限定されるものではない。
【００３８】
図１に、本実施形態に係る画像処理装置１の構成を概略的に示す。同図に示すように、画像処理装置１は、画像取得部２と、位置合せ条件設定部３と、位置合せ演算部４と、画像出力部５とを有している。
【００３９】
画像取得部２は、撮像モダリティが互いに異なる２つの画像を取得する。通常は、ユーザの操作に応じてこれら２つの画像が入力され取得される。画像取得部２は、これら２つの画像のうち一方を、位置合せ処理において固定される目標画像Ｔに設定し、他方を、位置合せ処理において変更される対象画像Ｒに設定する。なお、ここでは、取得した元の対象画像をＲ、位置合せ処理中の対象画像をＲ′、位置合せ処理後の対象画像をＲ″で表すことにする。
【００４０】
位置合せ条件設定部３は、ユーザ(user)の操作に応じて、位置合せ演算の条件である最適化パラメータ(parameter)を設定する。最適化パラメータの設定については、後ほど詳述する。
【００４１】
位置合せ演算部４は、設定された最適化パラメータに従って位置合せ演算を行う。位置合せ演算は、対象画像Ｒ′に移動や変形などの微小な変更を加える処理と、目標画像Ｔと対象画像Ｒ′との位置合せの評価値を算出する処理とを、その評価値が改善される方向に繰り返し行うことにより成される。
【００４２】
本実施形態では、図１に示すように、位置合せ演算部４は、評価部（決定手段、算出手段）４１と、最適化部（変更手段）４２と、変換部（変更手段）４３と、補間部（変更手段）４４とを有している。なお、これは、位置合せ演算部を図式化する際によく用いられる構成図であり、評価部４１は、Cost Function、最適化部４２は、Optimizer、変換部４３は、Transformer、補間部４４は、Image Interpolatorとも呼ばれる。
【００４３】
位置合せ演算部４は、これら評価部４１〜補間部４４による処理を複数回繰り返し行うことにより、位置合せ演算を行う。なお、評価部４１は、発明における決定手段および算出手段の一例である。また、最適化部４２、変換部４３および補間部４４は、発明における変更手段の一例である。以下、これら評価部４１〜補間部４４の各部の処理について、画像を参照しながら詳しく説明する。
【００４４】
図２に、位置合せ演算の概念図を示す。図２では、位置合せの対象となる目標画像Ｔおよび対象画像Ｒの例として、腹部（肝臓）におけるUS画像およびMR画像を提示している。
【００４５】
評価部４１は、目標画像Ｔおよび対象画像Ｒ′のそれぞれに対して、座標ｐごとに、ｘ，ｙ，ｚの各方向の１次偏微分、すなわち勾配ベクトルを算出する。このとき、勾配ベクトルの算出前に、目標画像Ｔおよび対象画像Ｒ′のそれぞれに対して、血管などの構造物を強調するフィルタ処理と、ノイズを低減するフィルタ処理との少なくとも１つを行うことが好ましい。なお、これらのフィルタ処理としては、一般的な既存の構造物強調フィルタ処理やノイズ低減フィルタ処理のほか、特にUS画像に対してはUS画像特有のフィルタ処理を用いることができる。US画像特有のフィルタ処理としては、例えば、非特許文献３、Real-Time
Speckle Reduction and Coherence Enhancement in Ultrasound Imaging via Nonlinear
Anisotropic Diffusion, Khaled Z. Abd-Elmoniem,
Student Member, IEEE, Abou-Bakr M. Youssef, and
Yasser M. Kadah*, Member, IEEE(IEEE TRANSACTIONS ON
BIOMEDICAL ENGINEERING, VOL. 49, NO. 9, SEPTEMBER 2002)にて開示されているフィルタ処理などが考えられる。
【００４６】
次に、目標画像Ｔおよび対象画像Ｒ′のそれぞれに対して、座標ｐごとに、先に求めた勾配ベクトル∇Ｔ(ｐ)、∇Ｒ′(ｐ)に基づいて、閾値判定付き標準勾配場（特徴ベクトル）ｎ(Ｔ,ｐ)、ｎ(Ｒ′,ｐ)を、次式にしたがって算出する。
【００４７】
2014225218.tif 000006
【００４８】
ここでTHNormは、勾配ベクトルの長さ(Norm)に対する閾値である。全ての勾配ベクトルに対してベクトルのNormを算出し、Normが閾値未満ならばノイズに起因する勾配ベクトルであるとみなし、閾値判定付き標準勾配場を０ベクトルとする。これにより閾値判定付き標準勾配場には、構造物に起因するベクトル、言い換えると位置合せにおいて有意な情報を与えるベクトルのみが抽出される。なお、閾値THNormは、例えば、実験や経験を基に決定する。また、閾値THNormは、目標画像Ｔと対象画像Ｒ′とで共通の値を用いてもよいし、別々に独立した値を用いてもよい。
【００４９】
図３は、従来法による標準勾配場の例を示している。一方、図４は、本提案法（閾値THNormを経験的に決定）による閾値判定付き標準勾配場の例を示している。本提案法では、ノイズに起因する勾配場が除去されており、位置合せの目標となり得る臓器内の主要な構造物に起因する勾配場のみが支配的になっているのが分かる。
【００５０】
特徴ベクトルを算出したら、目標画像Ｔおよび対象画像Ｒ′における互いに対応する座標ｐごとに、勾配場の揃い具合を表す評価関数Ｄ（Ｔ，Ｒ′）の値（評価値）を算出する。評価関数Ｄは、次式のように定義する。
【００５１】
2014225218.tif 000007
すなわち、評価関数Ｄ（Ｔ，Ｒ′）は、目標画像Ｔおよび対象画像Ｒ′における各座標ｐでの「特徴ベクトル同士の内積の二乗」（相関値）の積算値、すなわち加算値あるいは総和を、画像の大きさ、すなわち画像の総画素数で規格化したものである。
【００５２】
仮に、ある２つのベクトルが平行ならば内積値は最大値を、垂直ならば内積値は最小値をとる。したがって、(9)式の評価関数Ｄが最大値となるように２つの画像の位置関係を最適化していけば、その２つの画像の勾配場が平行に近づいていき、位置合せが達成される。
【００５３】
なお、評価部４１による評価関数の値の算出は、対象画像Ｒ′が変更され、更新される度に行われる。
【００５４】
最適化部４２は、算出された評価関数の値を基に、対象画像Ｒ′の変換パラメータの決定と変更とを繰返し行い、変換パラメータを最適化する。この変換パラメータは、対象画像Ｒ′の平行移動、回転、変形などの変更における向きと変更量を規定するパラメータである。変換パラメータの最適化のルールは、例えば、以下の通りである。
【００５５】
初回は、変換パラメータを事前に設定された方向にシフト(shift)させる。２回目以降は、評価関数Ｄの値が前回に比べて改善している場合には、前回と同方向に変換パラメータをシフトさせ、悪化している場合には、逆方向に変換パラメータをシフトさせる。変換パラメータをシフトさせていくときのステップ(step)幅は、事前に設定された「最大ステップ幅」から始める。そして、逆方向に変換パラメータをシフトさせるときに、その時点でのステップ幅に、事前に設定された「緩和係数」を乗じて成るステップ幅を新たなステップ幅とする。
【００５６】
なお、最適化部４２による変換パラメータの最適化の方法としては、例えば、Gradient Descent法、Conjugate gradient法、LBFGS法、Gauss-Newton法、Levenberg-Marquardt法、Amoeba法などを用いることができる。
【００５７】
変換部４３は、最適化部４２で決定した変換パラメータにしたがって、対象画像Ｒ′に変換処理を行う。
【００５８】
なお、変換部４３による対象画像Ｒ′に対する変換処理の方法としては、例えば、Euler Transform（オイラー変換）、Affine Transform（アフィン変換）、B-Spline Transform（Bスプライン変換）、Versor Transform（ベルソル変換）、Quaternion Transform（クォータニオン変換,または4元数変換）などを用いることができる。
【００５９】
補間部４４は、変換部４３で変換処理が成された画像に対して補間処理を行い、各座標における画素の画素値を算出する。対象画像Ｒ′に対して変換処理が行われた後は、１画素幅単位の平行移動など特殊な場合を除き、座標の位置に画素が重ならない場合が多い。そこで、このような補間処理により、各座標の画素値を求める必要がある。この補間処理が完了すると、対象画像Ｒ′は更新されたことになる。
【００６０】
なお、補間部４４による補間処理の方法としては、例えば、Nearest Neighbor Interpolation（最近傍補間）、Linear
Interpolation（線形補間）、B-Spline Interpolation（Bスプライン補間）などを用いることができる。
【００６１】
位置合せ演算部４は、評価部４１〜補間部４４による処理を複数回繰返し行う。そして、評価関数の前回の値と今回の値との差分が、事前に設定された「許容誤差(tolerance)」以下になった場合、ステップ幅が、事前に設定された「最小ステップ幅」未満になった場合、または、変換パラメータをシフトする繰返し演算の回数が、事前に設定された「繰返し演算上限回数」に達した場合には、この繰返し演算を終了させる。
【００６２】
位置合せ条件設定部３は、位置合せ演算条件として、最適化部４２による変換パラメータをシフトする繰返し演算における最適化パラメータを設定する。最適化パラメータとしては、上記の「最大ステップ幅」、「最小ステップ幅」、「繰返し演算上限回数」、「許容誤差」、および「緩和係数」などが含まれる。
【００６３】
画像出力部５は、位置合せ演算終了後に、目標画像Ｔと、位置合せされた対象画像Ｒ″とを出力する。
【００６４】
以下、本実施形態に係る画像処理装置による画像位置合せ処理の流れについて説明する。図５は、本実施形態に係る画像処理装置による画像位置合せ処理のフロー図である。
【００６５】
ステップＳ１では、画像取得部２が、位置合せする２つの画像を取得する。そして、取得したこれら２つの画像の一方を目標画像Ｔとして設定し、他方を対象画像Ｒとして設定する。
【００６６】
ステップＳ２では、位置合せ条件設定部３が、最適化パラメータを設定する。最適化パラメータとしては、最大ステップ幅、最小ステップ幅、繰返し演算の上限回数、位置合せの許容誤差、緩和係数などが含まれる。
【００６７】
ステップＳ３では、最適化部４２が、対象画像Ｒ′の変換パラメータ（平行移動、回転移動、変形など）を決定する。評価関数Ｄの値が前回に比べて改善している場合は、前回と同方向に変換パラメータをシフトさせる。逆に悪化している場合には、逆方向に変換パラメータをシフトさせる。初回は、適当な変換パラメータを決定する。
【００６８】
ステップＳ４では、変換部４３が、決定された変換パラメータにしたがって、対象画像Ｒ′に変換処理を行う。
【００６９】
ステップＳ５では、補間部４４が、対象画像Ｒ′に対して補間処理を行う。
【００７０】
ステップＳ６では、評価部４１が、評価関数Ｄ（Ｔ，Ｒ′）の値、すなわち評価値を算出する。
【００７１】
ステップＳ７では、位置合せ演算部４が、繰返し演算の終了判定を行う。繰返し回数が所定数に達した場合、または、評価関数の値が実質的に収束した場合（許容誤差に到達する、または変換パラメータをシフト量が最小ステップ幅未満に到達した場合）には、繰返し演算を終了し、ステップＳ８に進む。それ以外の場合には、ステップＳ３に戻る。
【００７２】
ステップＳ８では、画像出力部５が、最終的な変換パラメータを対象画像に適用し、位置合せ後画像を出力する。
【００７３】
図６に、従来法と本提案法（閾値THNormを経験的に決定）とによる画像位置合せの例を示す。図６において、左から順にAxial, Sagittal, Coronalの各断面を、上から順に目標画像（この例ではUS画像）、対象画像（この例ではMR画像）、従来法による位置合せ画像、本提案法による位置合せ画像である。図中の矢印は、従来法と本提案法で位置合せ精度の向上が顕著に視認できる箇所を示している。図６に示すように、本提案法は、従来法に比べて位置合せの精度が向上する。また、位置合せの計算結果が所望の精度に到達しているかどうかを容易に判別できる。
【００７４】
（変形例）
上記閾値THNormは、目標画像Ｔ及び対象画像Ｒ′の少なくとも一方に基づくノイズ量の推定結果を基に決定してもよい。
【００７５】
図７は、目標画像Ｔ及び対象画像Ｒ′の少なくとも一方に基づいて画像のノイズ量を推定する処理のフロー図である。
【００７６】
ステップＴ１では、目標画像Ｔ及び対象画像Ｒ′の少なくとも一方において、ノイズ量を推定する対象領域を設定する。画像中には診断対象の臓器だけではなく、その他の臓器が含まれることもあるため、これらを除外するように対象領域を設定することが望ましい。また、臓器あるいは医用画像においては、ノイズ量が一定ではなく画像位置によって異なる場合もあるので、ノイズ量を推定する対象領域を、重点的に診断したい範囲（例えば病変部の周辺など）に限定することで、より効果的にノイズ量を推定することができる。上記対象領域として切り出す範囲は、例えば、画像中央に位置する第１の部分範囲、病変部周辺に位置する第２の部分範囲、あるいは、目標画像Ｔと対象画像Ｒ′との間で共通する解剖学的特徴点の周辺に位置する第３の部分範囲などとすることができる。第１の部分範囲は、例えば、ＦＯＶに対応する範囲のうち面積がＦＯＶ全体の１/２〜１/３に相当する画像中央の範囲などとすることができる。第２の部分範囲は、例えば、腫瘍などの病変部の周辺であって、体積が５０〜１００ｍｍ3となる範囲などとすることができる。また、第３の部分範囲は、例えば、骨形状や血管形状もしくは血管分岐など共通する解剖学的特徴点の周辺であって、体積が５０〜１００ｍｍ3となる範囲などとすることができる。対象領域の設定は、手動であっても自動であってもよい。自動設定の場合には、必要に応じて、病変部や解剖学的特徴点を検出する処理を実行する。なお、本ステップを実行しなくても、以降の処理は実行可能であるが、正確なノイズ量を推定には、本ステップを実行することが望ましい。本ステップを実行しない場合、上記対象領域は、画像全体の領域とする。
【００７７】
ステップＴ２では、ステップＴ１で設定した対象領域に対してスムージング処理を行う。スムージング処理は、ガウシアンフィルタなどの既知のフィルタリングでよい。
【００７８】
ステップＴ３では、ステップＴ２で生成したスムージング処理済み画像に対して、各画像位置における勾配強度（Gradient Magnitude）を算出し、各画像位置に対応した画素の画素値をその画像位置での勾配強度とする勾配強度画像を生成する。勾配強度は、よく知られているように、ｘ，ｙ，ｚの各方向の１次偏微分、すなわち勾配ベクトル（Gradient Vector）の大きさを算出することで得られる。
【００７９】
ステップＴ４では、ステップ３で生成した勾配強度画像に対して、各画像位置での勾配強度のヒストグラムを算出し、ヒストグラムの特徴量を解析する。図８に、勾配強度のヒストグラムの例を示す。図８のヒストグラムにおいて、横軸は勾配強度を、縦軸は頻度の相対度数をそれぞれ表している。ノイズに起因する勾配の場合、その勾配強度は相対的に低値であり、その度数はヒストグラム上で低値に集中する。一方、構造物に起因する勾配の場合、その勾配強度は相対的に中値/高値をとり、その度数は広範囲に分布する。よって、図８のヒストグラムにおいて、勾配強度が低い領域に現れる山は、ノイズに起因する勾配に対応していると考えることができる。
【００８０】
ここで解析する特徴量は、以下の通りである（図９，図１０参照）。
・ヒストグラム分布上に現れるいくつかのピークのうち、最も勾配強度が低いピークにおいて最頻値を与える勾配強度（図９のMfmax）
・ヒストグラム分布上におけるMfmaxから見て低値側の半値半幅（図９のHWHML）
・ヒストグラム分布上におけるMfmaxから見て高値側の半値半幅（図９のHWHMR）
・ヒストグラム分布上における勾配強度が0からHWHMRまでの範囲の重心に相当する勾配強度（図１０のHWHMmom1）
・ヒストグラム分布上における勾配強度がHWHMLからHWHMRまでの範囲の重心に相当する勾配強度（図１０のHWHMmom2）
【００８１】
ステップＴ５では、ノイズ量を算出する。ノイズ量Nqtは、以下の（ａ）〜（ｆ）式のいずれか１つ、または、複数の組合せの平均値から算出される。
Nqt = Mfmax + σ×HWHML                     （ａ）
Nqt = Mfmax + σ×HWHMR                     （ｂ）
Nqt = Mfmax + σ×(HWHML+HWHMR)/2   （ｃ）
Nqt = σ×Mfmax                      （ｄ）
Nqt = σ×Mmom1                      （ｅ）
Nqt = σ×Mmom2                      （ｆ）
【００８２】
ここで、σは任意の定数であり、２．０〜３．０程度の値が望ましい。後述の画像の位置合せ（Image Registration）のプロセスにおいて、ユーザが比較的大きな構造物のみを重点的に位置合せすることを意図する場合には、σを大き目の値に設定することが望ましく、ユーザが比較的小さな構造物も含めて位置合せすることを意図する場合には、σを小さ目の値に設定することが望ましい。
【００８３】
閾値THNormは、この推定したノイズ量Nqtに基づいて決定するが、本例では、ノイズ量Nqtをそのまま閾値THNormとして用いる。
【００８４】
なお、閾値THNormを目標画像Ｔ及び対象画像Ｒ′で共通に用いる場合には、（１）目標画像Ｔにおけるノイズ量の推定結果から閾値THNormを決定する、（２）目標画像Ｒ′におけるノイズ量の推定結果から閾値THNormを決定する、（３）目標画像Ｔにおけるノイズ量の推定結果と目標画像Ｒ′におけるノイズ量の推定結果を総合的に用いて閾値THNormを決定する等の方法により閾値を決定する。また、閾値THNormを目標画像Ｔと対象画像Ｒ′とで別々に決定する場合には、目標画像Ｔにおけるノイズ量の推定結果から目標画像Ｔ用の閾値THNormを決定し、対象画像Ｒ′におけるノイズ量の推定結果から対象画像Ｒ′用の閾値THNormを決定する。
【００８５】
図１１は、従来法による標準勾配場の例と本提案法（閾値THNormをノイズ量の自動推定結果から決定）による閾値判定付き標準勾配場の例とを示す。この図では、左端の画像（腹部US画像）における標準勾配場を算出した例を示している。図１１から理解されるように、従来法では構造物に由来する勾配場だけでなくノイズに由来する勾配場も含まれているが、本提案法では構造物に由来する勾配場のみが的確に描出されている（図１１の白線参照）。
【００８６】
図１２に、従来法と本提案法（閾値THNormをノイズ量の自動推定結果から決定）とによる画像位置合せの他の例を示す。図１２中の画像は、それぞれ、目標画像（この例では腹部US画像）、対象画像（この例では腹部MR画像）、従来法による位置合わせ画像、目標画像と従来法による位置合わせ画像とのFusion表示、本提案法による位置合わせ画像、目標画像と本提案法による位置合わせ画像とのFusion表示をそれぞれ表示している。US画像、MR画像ともに3Dスキャンされた画像であるが、図１２では、代表的な１断面のみを表示している。従来法の場合には白線で示した血管の視認がやや困難であるが、本提案法では明瞭に視認できる。また、図１２中の白点線で示されるように、従来法の場合には、目標画像と位置合せ画像との間に補正しきれていない位置ずれ（misregistration）が観察されるのに対し、本提案法では、ほとんど位置ずれ無しに位置合せが達成されている。
【００８７】
このように、本実施形態に係る画像処理装置によれば、以下に説明する通り、精度の高い位置合せが可能になる。
【００８８】
従来法では、ノイズに起因する勾配ベクトルは、標準勾配場上で微小な値をとる。しかし、(5),(6)式のように、評価関数はすべての座標における内積または外積の積算値をとる。そのため、内積または外積一つ一つは微小な値であっても、その絶対数が多ければ評価関数に与える影響は無視できない。一方、本提案法では、勾配ベクトルの大きさが閾値THNorm未満のときは、閾値判定付き標準勾配場は、０ベクトルになる。そのため、閾値判定付き標準勾配場は、主要な構造物による成分のみで構成されており、ノイズに起因する勾配ベクトルの影響が効果的に抑制される。よって、本提案法では、位置合せ精度の向上が期待できる。
【００８９】
また、従来法では、(5),(6)式のように、２つの標準勾配場の内積または外積の二乗の積算値を算出しており、評価関数の値は、画像の大きさ（総画素数）に依存する。このため、評価関数の絶対値から位置合せの精度を予測することが困難である。本提案法では、評価関数が画像の総画素数で規格化されているため、評価関数の絶対値と位置合せの精度が相関性を持つ。仮に、位置合せする２つの画像が完全に同一であるならば、(8)式の評価関数Ｄの値は1.0となる。また仮に、位置合せする２つの画像において、50％の画素が同一であり、残り50％の画素にまったく類似性がなければ、(8)式の評価関数Ｄの値は0.5となる。したがって、あらかじめユーザが望む位置合せ精度を指定しておけば、位置合せの計算結果が所望の精度を満たしているかどうかを、簡単に判別することができる。
【００９０】
マルチモダリティの画像による診断では必ずしも２つの画像は同時期に撮影されているとは限らず、有意な時差が存在する可能性も考えられる。例えば、一方の画像が術前に撮影され、他方の画像が術後に撮影されている状況も想定される。この場合、一方の画像上に存在する病変が、他方の画像に必ずしも存在するとは限らない。よって一方の画像上にのみ病変に起因する勾配ベクトルが存在する場合には、このような勾配ベクトルは評価関数の計算からは除外されることが望ましい。しかし、従来法では、一方の画像上にのみ病変に起因する勾配ベクトルが存在する場合、他方の画像における当該病変に対応する座標には、ノイズ起因の勾配ベクトルが存在している。したがって、これらの標準勾配場の内積値は何らかの値を持ってしまい、評価関数に影響を与えてしまう。これに対して本提案法では、他方の画像における当該病変に対応する座標には、０ベクトルが存在しているので、内積値は０となり、評価関数に影響を与えない。よって、本提案法では、位置合せ精度の向上が期待できる。
【００９１】
また、位置合せを行う２つの画像の組合せとしては、MR画像とUS画像の組合せだけでなく、CT画像とUS画像や、MR画像とCT画像の組合せなど、あらゆる撮像モダリティの画像に適用できる。ただし、標準勾配場を用いる従来法による位置合せ手法は、特に、位置合せ対象の画像にUS画像が含まれる場合に有効な手法であることから、本提案法による位置合せ手法も、US画像を含む位置合せに対して特に有効である。
【００９２】
また、閾値THNormを上述の変形例のようにノイズ量の自動推定結果から決定するようにすれば、画像ごとに閾値THNormを最適化することができ、より高精度な位置合せが期待できる。
【００９３】
なお、本実施形態では、評価関数は、特徴ベクトル同士の内積の二乗を含むものであるが、特徴ベクトル同士の内積の絶対値のＮ乗（Ｎは自然数）を含むものであってもよい。
【００９４】
また、評価関数は、特徴ベクトル同士の内積の二乗または内積の絶対値のＮ乗を、画像の大きさで規格化したものが好ましいが、規格化しないものであってもよい。位置合せする画像の大きさ、すなわち画像の総画素数がほとんど変わらない場合には、このようにしても、位置合せの精度をある程度保持できる。
【００９５】
また、本実施形態は、発明を、撮像モダリティが互いに異なる画像同士の位置合せに適用した例であるが、撮像モダリティが同一であって撮像の時相が互いに異なる画像同士の位置合せに適用することもできる。このような画像としては、例えば、手術前後の画像や造影撮影における早期相と後期相の画像などが考えられる。また、発明は、医用画像だけでなく、工業用の画像にも適用可能である。
【００９６】
また、本実施形態は、画像処理装置であるが、コンピュータをこのような画像処理装置として機能させるためのプログラムもまた発明の実施形態の一例である。
【符号の説明】
【００９７】
１  画像処理装置
２  画像取得部
３  位置合せ制御条件設定部
４  位置合せ演算部
４１  評価部
４２  最適化部
４３  変換部
４４  補間部
５  画像出力部
【図１】
2014225218.tif 000008
【図５】
2014225218.tif 000009
【図８】
2014225218.tif 000010
【図９】
2014225218.tif 000011
【図１０】
2014225218.tif 000012
【図２】
2014225218.tif 000013
【図３】
2014225218.tif 000014
【図４】
2014225218.tif 000015
【図６】
2014225218.tif 000016
【図７】
2014225218.tif 000017
【図１１】
2014225218.tif 000018
【図１２】
2014225218.tif 000019