JP5889265B2 - 画像処理方法および装置並びにプログラム - Google Patents

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

Info

Publication number
JP5889265B2
JP5889265B2 JP2013230466A JP2013230466A JP5889265B2 JP 5889265 B2 JP5889265 B2 JP 5889265B2 JP 2013230466 A JP2013230466 A JP 2013230466A JP 2013230466 A JP2013230466 A JP 2013230466A JP 5889265 B2 JP5889265 B2 JP 5889265B2
Authority
JP
Japan
Prior art keywords
value
image
gradient
image processing
processing apparatus
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.)
Active
Application number
JP2013230466A
Other languages
English (en)
Other versions
JP2014225218A (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.)
GE Medical Systems Global Technology Co LLC
Original Assignee
GE Medical Systems Global Technology Co LLC
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 GE Medical Systems Global Technology Co LLC filed Critical GE Medical Systems Global Technology Co LLC
Priority to JP2013230466A priority Critical patent/JP5889265B2/ja
Priority to US14/258,081 priority patent/US9256916B2/en
Priority to CN201410162145.1A priority patent/CN104156940B/zh
Publication of JP2014225218A publication Critical patent/JP2014225218A/ja
Application granted granted Critical
Publication of JP5889265B2 publication Critical patent/JP5889265B2/ja
Active 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
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/751—Comparing pixel values or logical combinations thereof, or feature values having positional relevance, e.g. template matching
    • G06V10/7515—Shifting the patterns to accommodate for positional errors
    • 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/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • 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/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • 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/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Artificial Intelligence (AREA)
  • Health & Medical Sciences (AREA)
  • Multimedia (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Description

本発明は、同一の対象を含む複数の画像を位置合せする技術に関する。
従来、撮像装置として、例えば、磁気共鳴装置（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）したものである。つまり、標準勾配場は、画素値あるいは輝度値の大小や勾配の大きさに依存せず、勾配の方向だけを表す特徴量である。仮に、ある２つの画像において互いに対応する位置に同じ方向の標準勾配場が発生しているならば、これら２つの画像の位置は合っていると見なすことができる。したがって、この手法では、標準勾配場が示す方向の揃い具合を最適化することで、位置合せを行うことができる。
標準勾配場ｎは、数学的には以下のように記述される。
(1)式において、∇Ｉ(ｐ)は、画像Ｉの座標ｐにおける勾配ベクトルを表している。また、||∇Ｉ（ｐ）||は、その勾配ベクトルの長さ(Norm)を表している。座標ｐは、ｘ，ｙ，ｚの各成分からなる。
(1)式から理解されるように、標準勾配場ｎは、それ単独では、構造物に起因する勾配とノイズに起因する勾配とを分別できない。そこで、ノイズ（noise）に起因する勾配の影響を抑制するために、非特許文献２において、次式のような手法が提案されている。
ここで、εは画像のノイズ量に応じて任意に選択する定数、である。従来法では、(4)式に従ってノイズ項εを設定し、ノイズ項を考慮して勾配ベクトルを規格化している。
位置合せの指標となる評価関数Ｄは、次式（いずれか一方でよい）で定義される。
(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)にて開示されているフィルタ処理などが考えられる。
次に、目標画像Ｔおよび対象画像Ｒ′のそれぞれに対して、座標ｐごとに、先に求めた勾配ベクトル∇Ｔ(ｐ)、∇Ｒ′(ｐ)に基づいて、閾値判定付き標準勾配場（特徴ベクトル）ｎ(Ｔ,ｐ)、ｎ(Ｒ′,ｐ)を、次式にしたがって算出する。
ここでTHNormは、勾配ベクトルの長さ(Norm)に対する閾値である。全ての勾配ベクトルに対してベクトルのNormを算出し、Normが閾値未満ならばノイズに起因する勾配ベクトルであるとみなし、閾値判定付き標準勾配場を０ベクトルとする。これにより閾値判定付き標準勾配場には、構造物に起因するベクトル、言い換えると位置合せにおいて有意な情報を与えるベクトルのみが抽出される。なお、閾値THNormは、例えば、実験や経験を基に決定する。また、閾値THNormは、目標画像Ｔと対象画像Ｒ′とにおいて、共通の値であってもよいし、別々に独立した値であってもよい。
図３は、従来法による標準勾配場の例を示している。一方、図４は、本提案法（閾値THNormを経験的に決定）による閾値判定付き標準勾配場の例を示している。本提案法では、ノイズに起因する勾配場が除去されており、位置合せの目標となり得る臓器内の主要な構造物に起因する勾配場のみが支配的になっているのが分かる。
特徴ベクトルを算出したら、目標画像Ｔおよび対象画像Ｒ′における互いに対応する座標ｐごとに、勾配場の揃い具合を表す評価関数Ｄ（Ｔ，Ｒ′）の値（評価値）を算出する。評価関数Ｄは、次式のように定義する。
すなわち、評価関数Ｄ（Ｔ，Ｒ′）は、目標画像Ｔおよび対象画像Ｒ′における各座標ｐでの「特徴ベクトル同士の内積の二乗」（相関値）の積算値、すなわち加算値あるいは総和を、画像の大きさ、すなわち画像の総画素数で規格化したものである。
仮に、ある２つのベクトルが平行ならば内積値は最大値を、垂直ならば内積値は最小値をとる。したがって、(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を最適化することができ、より高精度な位置合せが期待できる。
なお、本実施形態では、評価関数は、特徴ベクトル同士の内積の二乗を含むものであるが、特徴ベクトル同士の内積の絶対値のＮ乗（Ｎは自然数）を含むものであってもよい。
また、評価関数は、特徴ベクトル同士の内積の二乗または内積の絶対値のＮ乗を、画像の大きさで規格化したものが好ましいが、規格化しないものであってもよい。位置合せする画像の大きさ、すなわち画像の総画素数がほとんど変わらない場合には、このようにしても、位置合せの精度をある程度保持できる。
また、本実施形態は、発明を、撮像モダリティが互いに異なる画像同士の位置合せに適用した例であるが、撮像モダリティが同一であって撮像の時相が互いに異なる画像同士の位置合せに適用することもできる。このような画像としては、例えば、手術前後の画像や造影撮影における早期相と後期相の画像などが考えられる。また、発明は、医用画像だけでなく、工業用の画像にも適用可能である。
また、本実施形態は、画像処理装置であるが、コンピュータをこのような画像処理装置として機能させるためのプログラムもまた発明の実施形態の一例である。
１ 画像処理装置
２ 画像取得部
３ 位置合せ制御条件設定部
４ 位置合せ演算部
４１ 評価部
４２ 最適化部
４３ 変換部
４４ 補間部
５ 画像出力部

Claims (18)

  1. 同一の被写体を含む第１および第２の画像における座標ごとに、該座標における画素値の勾配に係る特徴ベクトルを決定する第１の工程であって、該座標における画素値の勾配ベクトルの大きさが所定の閾値以上であるときには、前記勾配ベクトルを該勾配ベクトルの大きさで規格化したものを該座標における特徴ベクトルとして決定し、前記勾配ベクトルの大きさが前記所定の閾値未満であるときには、０（ゼロ）ベクトルを該座標における特徴ベクトルとして決定する第１の工程と、
    前記第１および第２の画像における互いに対応する座標ごとに、該座標における特徴ベクトル同士の内積の絶対値をＮ（Ｎは自然数）乗して成る値に相当する相関値を算出し、前記座標ごとに算出された相関値の積算値を含む評価値を求める第２の工程と、
    前記評価値がより大きくなるように前記第２の画像を変更する第３の工程と、を複数回繰り返し行って、前記第１および第２の画像の位置合せを行う画像処理方法であって、
    前記第１および第２の画像の少なくとも一方の対象領域における各画素での勾配強度を求め、該勾配強度のヒストグラムに基づいて前記閾値を決定する工程を有している画像処理方法。
  2. 同一の被写体を含む第１および第２の画像における座標ごとに、該座標における画素値の勾配ベクトルの大きさが所定の閾値以上であるときには、前記勾配ベクトルを該勾配ベクトルの大きさで規格化したものを該座標における特徴ベクトルとして決定し、前記勾配ベクトルの大きさが前記所定の閾値未満であるときには、０（ゼロ）ベクトルを該座標における特徴ベクトルとして決定する決定手段と、
    前記第１および第２の画像における互いに対応する座標ごとに、該座標における特徴ベクトル同士の内積の絶対値をＮ（Ｎは自然数）乗して成る値に相当する相関値を算出し、前記座標ごとに算出された相関値の積算値を含む評価値を求める算出手段と、
    前記評価値がより大きくなるように前記第２の画像を変更する変更手段と、を有し、
    前記決定手段、算出手段および変更手段による処理を複数回繰り返し行って、前記第１および第２の画像の位置合せを行う画像処理装置であって、
    前記第１および第２の画像の少なくとも一方の対象領域における各画素での勾配強度を求め、該勾配強度のヒストグラムに基づいて前記閾値を決定する手段を有している画像処理装置。
  3. 前記勾配ベクトルは、画像における各座標軸方向の１次偏微分である、請求項２に記載の画像処理装置。
  4. 前記相関値は、前記特徴ベクトル同士の内積を二乗して成る値である、請求項３に記載の画像処理装置。
  5. 前記評価値は、前記相関値の積算値を前記第１または第２の画像の大きさで規格化して成る値である、請求項２から請求項４のいずれか一項に記載の画像処理装置。
  6. 前記変更は、平行移動、回転、および変形のうち少なくとも一つを含む、請求項２から請求項５のいずれか一項に記載の画像処理装置。
  7. 前記決定手段は、前記第１および第２の画像のそれぞれに構造物強調フィルタ処理および／またはノイズ低減処理を行って得られる画像に基づいて、前記勾配ベクトルを求める、請求項２から請求項６のいずれか一項に記載の画像処理装置。
  8. 前記第１および第２の画像は、医用画像である、請求項２から請求項７のいずれか一項に記載の画像処理装置。
  9. 前記第１および第２の画像は、撮像モダリティが互いに異なる、請求項８に記載の画像処理装置。
  10. 前記第１および第２の画像の一方は、超音波画像（US画像）である、請求項９に記載の画像処理装置。
  11. 前記決定手段、算出手段および変更手段による処理を、該処理を行った回数が所定数に達するまで、または、前記評価値が実質的に収束するまで、繰り返し行う、請求項２から請求項１０のいずれか一項に記載の画像処理装置。
  12. 前記対象領域は、前記第１の画像と前記第２の画像との間に共通する解剖学的特徴点の周辺領域である、請求項２から請求項１１のいずれか一項に記載の画像処理装置。
  13. 前記対象領域は、前記少なくとも一方の画像の中央領域である、請求項２から請求項１１のいずれか一項に記載の画像処理装置。
  14. 前記対象領域は、病変部の周辺領域である、請求項２から請求項１１のいずれか一項に記載の画像処理装置。
  15. 前記対象領域は、全画像領域である、請求項２から請求項１１のいずれか一項に記載の画像処理装置。
  16. 前記閾値を決定する手段は、前記対象領域におけるノイズ量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までの範囲の重心に相当する勾配強度である。
  17. 前記σは、２〜３の値である、請求項１６に記載の画像処理装置。
  18. コンピュータを、請求項２から請求項１７のいずれか一項に記載の画像処理装置として機能させるためのプログラム。
JP2013230466A 2013-04-22 2013-11-06 画像処理方法および装置並びにプログラム Active JP5889265B2 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2013230466A JP5889265B2 (ja) 2013-04-22 2013-11-06 画像処理方法および装置並びにプログラム
US14/258,081 US9256916B2 (en) 2013-04-22 2014-04-22 Image processing method and apparatus and program
CN201410162145.1A CN104156940B (zh) 2013-04-22 2014-04-22 图像处理方法及装置

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2013088963 2013-04-22
JP2013088963 2013-04-22
JP2013230466A JP5889265B2 (ja) 2013-04-22 2013-11-06 画像処理方法および装置並びにプログラム

Publications (2)

Publication Number Publication Date
JP2014225218A JP2014225218A (ja) 2014-12-04
JP5889265B2 true JP5889265B2 (ja) 2016-03-22

Family

ID=51729035

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013230466A Active JP5889265B2 (ja) 2013-04-22 2013-11-06 画像処理方法および装置並びにプログラム

Country Status (3)

Country Link
US (1) US9256916B2 (ja)
JP (1) JP5889265B2 (ja)
CN (1) CN104156940B (ja)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107004259A (zh) * 2014-12-10 2017-08-01 皇家飞利浦有限公司 多对比度成像中的统计加权正则化
JP6374535B2 (ja) * 2015-02-02 2018-08-15 富士フイルム株式会社 操作装置、追尾システム、操作方法、及びプログラム
JP6775294B2 (ja) * 2015-11-12 2020-10-28 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像処理装置および方法並びにプログラム
US11538176B2 (en) 2015-12-15 2022-12-27 Koninklijke Philips N.V. Image processing systems and methods
EP3375399B1 (en) 2016-10-05 2022-05-25 NuVasive, Inc. Surgical navigation system
CN108280850A (zh) * 2018-03-05 2018-07-13 北京中科嘉宁科技有限公司 一种基于B样条和Levenberg-Marquardt优化的快速图像配准方法
US11430140B2 (en) * 2018-09-18 2022-08-30 Caide Systems, Inc. Medical image generation, localizaton, registration system
JP7290274B2 (ja) 2019-07-04 2023-06-13 東芝エネルギーシステムズ株式会社 荷電粒子の出射制御装置、方法及びプログラム
CN113627446B (zh) * 2021-08-18 2023-10-31 成都工业学院 基于梯度向量的特征点描述算子的图像匹配方法及***
JPWO2023054089A1 (ja) * 2021-10-01 2023-04-06

Family Cites Families (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004227519A (ja) * 2003-01-27 2004-08-12 Matsushita Electric Ind Co Ltd 画像処理方法
JP4617883B2 (ja) * 2004-01-06 2011-01-26 ソニー株式会社 画像処理装置および方法、プログラム並びに記録媒体
US8090429B2 (en) 2004-06-30 2012-01-03 Siemens Medical Solutions Usa, Inc. Systems and methods for localized image registration and fusion
JP4559844B2 (ja) * 2004-12-27 2010-10-13 株式会社東芝 画像処理装置及び画像処理方法
US20070167784A1 (en) 2005-12-13 2007-07-19 Raj Shekhar Real-time Elastic Registration to Determine Temporal Evolution of Internal Tissues for Image-Guided Interventions
KR20070110965A (ko) 2006-05-16 2007-11-21 주식회사 메디슨 초음파 영상과 외부 의료영상의 합성 영상을디스플레이하기 위한 초음파 시스템
US8358818B2 (en) 2006-11-16 2013-01-22 Vanderbilt University Apparatus and methods of compensating for organ deformation, registration of internal structures to images, and applications of same
US8098911B2 (en) 2006-12-05 2012-01-17 Siemens Aktiengesellschaft Method and system for registration of contrast-enhanced images with volume-preserving constraint
IL191615A (en) * 2007-10-23 2015-05-31 Israel Aerospace Ind Ltd A method and system for producing tie points for use in stereo adjustment of stereoscopic images and a method for identifying differences in the landscape taken between two time points
EP2131212A3 (en) 2008-06-05 2011-10-05 Medison Co., Ltd. Non-Rigid Registration Between CT Images and Ultrasound Images
US8144947B2 (en) * 2008-06-27 2012-03-27 Palo Alto Research Center Incorporated System and method for finding a picture image in an image collection using localized two-dimensional visual fingerprints
GB0913930D0 (en) * 2009-08-07 2009-09-16 Ucl Business Plc Apparatus and method for registering two medical images
JP2013020294A (ja) * 2011-07-07 2013-01-31 Sony Corp 画像処理装置と画像処理方法およびプログラムと記録媒体
CN102231191B (zh) * 2011-07-17 2012-12-26 西安电子科技大学 基于asift的多模态图像特征提取与匹配方法
JP5591309B2 (ja) * 2012-11-29 2014-09-17 キヤノン株式会社 画像処理装置、画像処理方法、プログラム、及び記憶媒体

Also Published As

Publication number Publication date
US9256916B2 (en) 2016-02-09
CN104156940A (zh) 2014-11-19
CN104156940B (zh) 2017-05-17
JP2014225218A (ja) 2014-12-04
US20140314295A1 (en) 2014-10-23

Similar Documents

Publication Publication Date Title
JP5889265B2 (ja) 画像処理方法および装置並びにプログラム
JP7203895B2 (ja) 位相相関を用いた非剛体変形の存在下での医用画像のシングル及びマルチのモダリティ位置合わせ
JP5335280B2 (ja) 位置合わせ処理装置、位置合わせ方法、プログラム、及び記憶媒体
US9230331B2 (en) Systems and methods for registration of ultrasound and CT images
JP6386752B2 (ja) 医用画像処理装置、医用画像診断装置及び画像処理方法
US8897519B2 (en) System and method for background phase correction for phase contrast flow images
Alves et al. Computer image registration techniques applied to nuclear medicine images
JP2011041656A (ja) ボリュームデータ間の対応付け方法
JP2008511395A (ja) 一連の画像における動き修正のための方法およびシステム
JP6905323B2 (ja) 画像処理装置、画像処理方法、およびプログラム
JP2005169077A (ja) 第１データセットのデータポイントを第２データセットのデータポイントと関係付ける変換を見つける方法、およびその方法を実行するコンピュータプログラム
JP2016189946A (ja) 医用画像位置合わせ装置および方法並びにプログラム
Alam et al. Evaluation of medical image registration techniques based on nature and domain of the transformation
Hameeteman et al. Carotid wall volume quantification from magnetic resonance images using deformable model fitting and learning-based correction of systematic errors
JP2016041245A (ja) 医用画像処理装置及び医用画像処理方法
Gupta et al. A hybrid framework for registration of carotid ultrasound images combining iconic and geometric features
Noorda et al. Registration of CT to pre-treatment MRI for planning of MR-HIFU ablation treatment of painful bone metastases
Carvalho et al. Joint intensity‐and‐point based registration of free‐hand B‐mode ultrasound and MRI of the carotid artery
Shah et al. Quantification and visualization of mri cartilage of the knee: A simplified approach
McLeod et al. A near-incompressible poly-affine motion model for cardiac function analysis
Carminati et al. Reconstruction of the descending thoracic aorta by multiview compounding of 3-d transesophageal echocardiographic aortic data sets for improved examination and quantification of atheroma burden
Morais et al. Dense motion field estimation from myocardial boundary displacements
US9495610B2 (en) Technique for presuming positions of organs
JP2022052210A (ja) 情報処理装置、情報処理方法及びプログラム
JP7270331B2 (ja) 医用画像診断装置及び画像処理装置

Legal Events

Date Code Title Description
A625 Written request for application examination (by other person)

Free format text: JAPANESE INTERMEDIATE CODE: A625

Effective date: 20141224

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20150430

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20150518

RD02 Notification of acceptance of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7422

Effective date: 20150602

RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7424

Effective date: 20150609

A601 Written request for extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A601

Effective date: 20150817

A601 Written request for extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A601

Effective date: 20150917

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20151007

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: 20160119

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160216

R150 Certificate of patent or registration of utility model

Ref document number: 5889265

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250