JPWO2008059982A1 - 画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置 - Google Patents

画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置 Download PDF

Info

Publication number
JPWO2008059982A1
JPWO2008059982A1 JP2008544217A JP2008544217A JPWO2008059982A1 JP WO2008059982 A1 JPWO2008059982 A1 JP WO2008059982A1 JP 2008544217 A JP2008544217 A JP 2008544217A JP 2008544217 A JP2008544217 A JP 2008544217A JP WO2008059982 A1 JPWO2008059982 A1 JP WO2008059982A1
Authority
JP
Japan
Prior art keywords
image
sectional image
cross
projection
projection image
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.)
Granted
Application number
JP2008544217A
Other languages
English (en)
Other versions
JP5190825B2 (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.)
Kyoto Institute of Technology NUC
Original Assignee
Kyoto Institute of Technology NUC
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 Kyoto Institute of Technology NUC filed Critical Kyoto Institute of Technology NUC
Priority to JP2008544217A priority Critical patent/JP5190825B2/ja
Publication of JPWO2008059982A1 publication Critical patent/JPWO2008059982A1/ja
Application granted granted Critical
Publication of JP5190825B2 publication Critical patent/JP5190825B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/037—Emission tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/448—Computed tomography involving metal artefacts, streaking artefacts, beam hardening or photon starvation

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Surgery (AREA)
  • Public Health (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Veterinary Medicine (AREA)
  • Biophysics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Pulmonology (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Abstract

何らかの形で求めた暫定の断面像ｆ（ｘ，ｙ）に対し、ｆ（ｘ，ｙ）から計算した投影像と測定された投影像との差を含む、評価関数Ｅを定義し、Ｅがおおむね減少するようにｆ（ｘ，ｙ）を変更することで、投影像に対応する断面像を求める計算機トモグラフィー装置及びプログラム。通常の計算機トモグラフィーで必要とされる逆投影操作が本質的に必要ではないことを特徴とし、メタルアーティファクトやエリアシングアーティファクトなどの除去・低減に特に効果がある。

Description

本発明は、対象物の放射線投影像から断面像を再構成する技術に関する。
計算機トモグラフィー（ＣＴ）は、Ｘ線などの光線を使って、様々な角度からの物体の投影像を撮影し、その後、計算によって断面像を得る技術である。図１に典型的なＸ線ＣＴ装置の模式図を示す。
典型的なＣＴ装置では、Ｘ線光源が回転移動することにより、測定対象の様々な角度からの投影像を取得する仕組みになっている。こうして得られた投影像から計算操作によって断面像を得るのがＣＴである。通常は、Ｆｉｌｔｅｒｅｄ Ｂａｃｋ−ｐｒｏｊｅｃｔｉｏｎ法（ＦＢＰ）によって投影像から断面像を変換操作によって求める。ＦＢＰでは、投影像にフィルタ（基本的には微分フィルタ）を施した後、もともとの投影方向に「逆投影」という操作を施すことで断面像を得る。このとき、微分フィルタはノイズや誤差を増幅し、アーティファクト（本来存在しない誤差や虚像）を生みやすいという特徴がある。さらに、逆投影操作はそのアーティファクトを断面像全体に伝播させる効果がある。そのため、ＣＴにおけるアーティファクトは問題のある部分だけに留まらず、断面像全体を損なうという点で致命的になることが多い。
アーティファクトの多くは、ＦＢＰに含まれるフィルタ操作や逆投影操作が原因なのだから、ＦＢＰを用いなければ、アーティファクトのない断面像が得られるはずである。ＦＢＰ以外で断面像を計算する方法として、Ａｌｇｅｂｒａｉｃ Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ Ｔｅｃｈｎｉｑｕｅ（ＡＲＴ）が歴史的に重要である。ＡＲＴはＥＢＰが考案される前から利用されるほど歴史のある方法である。ＡＲＴでは断面像の計算をフィッティング問題と捉え、断面像をパラメータ、投影像をフィッティング対象として、断面像から計算した投影像（ｐ）が、実験で求めた投影像（ｐ０）に一致するように逐次的に断面像を修正する。ＡＲＴでは（ｐ−ｐ０）が０になるように断面像を漸近的に改良してゆく点が特徴となっている。ＡＲＴはＦＢＰを全く用いないかわりに、断面像を求める計算に時間がかかるため、今では特殊な用途（地震波の解析など）にしか用いられない。ＡＲＴはＦＢＰほど極端なアーティファクトは現れないが、得られる断面像はＦＢＰの方が自然であることが多い。
さらに、フィルタ操作や逆操作に由来しないアーティファクトの要因として、投影像におけるデータの欠落・過少も挙げられる。データが少なければ結果として得られる断面像の画質も低下するのは当たり前である。ＦＢＰではデータの欠落・過少も致命的なアーティファクトを生み出すことが知られており、大きな問題と認識されている。フィッティングに基づくＣＴ（ＡＲＴを含む）はデータの欠落や過少に対してＦＢＰよりは強いとされている。しかしながら、ＣＴにおけるデータの欠落は、極端に「条件が悪い」問題とされており、ＡＲＴなどを用いても改善が得られにくい。ＡＲＴでは（ｐ−ｐ０）を誤差として考えるのであるが、条件が悪いときはフィッティングに失敗しやすい。いわゆる二乗誤差（ｐ−ｐ０）２の最小化を行った方がフィッティングとしては安定する。二乗誤差を最小化する方法としては、最小二乗法が最も標準的である。最小二乗法では、一辺がパラメータの個数だけある正方行列の逆行列を求める。ＣＴにおけるパラメータは断面像を構成するそれぞれの画素値であるため、パラメータの個数は莫大になる。１０００×１０００画素の断面像だと、パラメータの数は１００万個になり、行列の要素数は１兆個と莫大である。従って、正攻法の最小二乗法では行列が大きくなりすぎで破綻してしまう。そこで、正攻法の最小二乗法を避け、Ｓｉｍｕｌｔａｎｅｏｕｓ ｉｔｅｒａｔｉｖｅ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｔｅｃｈｎｉｑｕｅ（ＳＩＲＴ）やＩｔｅｒａｔｉｖｅ ｌｅａｓｔ ｓｑｕａｒｅｔｅｃｈｎｉｑｕｅ（ＩＬＳＴ）が考案された。ＡＲＴと同様に断面像の計算をフィッティング問題として捉えるのであるが、先に述べたように、正攻法の最小二乗法を避けなければならないため、ＳＩＲＴもＩＬＳＴも計算の途中に投影操作の逆演算としてＦＢＰが利用される。そのため、ＦＢＰに起因するフィルタ操作と逆投影にまつわる問題を根本的には解決しない。おそらくそれが理由となり、ＳＩＲＴやＩＬＳＴでは各種アーティファクトの「低減」効果が報告されているに留まる。
Ｙａｚｄｉ Ｍ，Ｇｉｎｇｒａｓ Ｌ，ＢｅａｕｌｉｅｕＬ．Ａｎ ａｄａｐｔｉｖｅ ａｐｐｒｏａｃｈ ｔｏ ｍｅｔａｌ ａｒｔｉｆａｃｔ ｒｅｄｕｃｔｉｏｎ ｉｎ ｈｅｌｉｃａｌ ｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｙ ｆｏｒ ｒａｄｉａｔｉｏｎ ｔｈｅｒａｐｙ ｔｒｅａｔｍｅｎｔ ｐｌａｎｎｉｎｇ：ｅｘｐｅｒｉｍｅｎｔａｌ ａｎｄ ｃｌｉｎｉｃａｌｓｔｕｄｉｅｓ．Ｉｎｔ Ｊ Ｒａｄｉａｔ Ｏｎｃｏｌ Ｂｉｏｌ Ｐｈｙｓ，６２：１２２４−１２３１，２００５．
以上の議論から、アーティファクトの無い断面像を得るには、ＦＢＰを用いない、データの欠落に対して強い、という性質が重要であることがわかる。そこで、ＡＲＴと同様にＦＢＰを用いず、フィッティングの評価関数として二乗誤差を利用しようというのは、誰でも思いつくことである。にもかかわらず、ＣＴにおけるアーティファクトの問題が根本的に解決していないのは、この単純な考えの実現は簡単ではないからである。第一の問題は、断面像の計算はフィッティングの問題としては大規模すぎて、計算に時間がかかることである。適切な高速化が必須である。第二の問題は、最小二乗法はフィッティングアルゴリズムとしては「条件が悪い」問題に対して弱いことである。そのため、最小二乗法に基づく既存のアルゴリズム（ＩＬＳＴなど）の延長では、アーティファクトの問題を解決できないと思われる。第三の問題は、ＳＩＲＴやＩＬＳＴがそうであるように、既存のアルゴリズムではＦＢＰを完全に排除することができていないことである。
本発明の第一の新規性は、フィッティング方法としてＳｉｍｕｌａｔｅｄ Ａｎｎｅａｌｉｎｇ法（ＳＡ）を採用したことである。ＳＡは公知の技術であるが、フィッティングに時間がかかることが知られており、ＣＴと本質的に相性が悪い。しかしながら、敢えてＳＡをＣＴに適用することで、ＦＢＰを全く用いずとも、二乗誤差の最小化による断面像の計算ができるようになった。また、ＳＡはデータの欠落など条件の悪いフィッティングに対してかなり安定であるという性質があり、この点においてもアーティファクトの根本的な解決に有利である。以上の事柄を鑑み、ＳＡをＣＴに適用することで、アーティファクトの根本的解決に効果を見出した点が、本発明の重要な新規性である。単純にＳＡのＣＴへの適用を考えたとき、断面像から投影像を求める計算を無数に繰り返すことになる。まともに計算すると、断面像から投影像を求める計算は、ＦＢＰとほぼ同じであり、ＳＡをＣＴに適用すると、ＦＢＰに比べて数百万倍の時間がかかってしまうことになる。現在の高速な計算機を以ってしても、「年」単位の計算時間を要してしまう。そこで、本発明では、式変形によって計算量を大幅に低減することで、この問題を解決した。この工夫は、ＳＡをＣＴに適用する際に必要な技術であり、本発明の重要な要素である。
第二の新規性は、フィッティングの評価関数として、二乗誤差だけでなく、アーティファクトを積極的に破壊するスムージング項とエントロピー項を導入した点である。従来のＣＴでは、ｐとｐ０との差あるいは二乗誤差だけを対象としていた。この場合、アーティファクトを御座なりにした方が良いフィッティング結果（小さな誤差）が達成できるという可能性が常にあり、現実の多くのケースではそのような条件になっていた。すなわち、アーティファクトを別のアーティファクトで相殺する状態になっていたのである。そのため、アーティファクトはなかなか消えないのである。本発明においても二乗誤差だけを評価関数とした場合には、アーティファクトは完全には消えなかった。この事実は、アーティファクトのない断面像を得るには、二乗誤差以外の要請が望ましいことを意味している。本発明では、統計力学の基本にのっとり、エントロピー項やスムージング項を考案した。これらの項は、断面像は滑らかで自然な濃淡画像であるべきだ、という当たり前の要請を数式の形で表現したものである。エントロピー項は、アーティファクトを破壊し、断面像全体の画質を均質にするように、フィッティングを誘導する。スムージング項は、エントロピー項で誘導される断面像のザラツキを押さえる効果がある。これらの項を導入することで、アーティファクトが消え、かつ自然な断面像が得られるようになった。エントロピー項とスムージング項は、それぞれ単独でもアーティファクトを低減する効果があるが、２つを組み合わせるとさらに効果的である。
なお、本発明においては、一般の定義よりも広い概念として、Ｘ線・可視光・電波を含む電磁波、電子や荷電粒子からなる粒子線、媒体の振動である音波等を全て含めて「放射線」と呼ぶことにする。
本発明の第一の効果として、データの欠落が要因になっているアーティファクトの低減効果が挙げられる。データの欠落が問題になる場合として、観察対象に不透明部位がある場合、投影角度に制限がある場合、投影角度の刻みが不均一である場合、コーンビームやヘリカルスキャンなどの三次元ＣＴの場合が挙げられる。
特に効果があることがわかっているのは、観察対象に不透明部位がある場合に現れるメタルアーティファクトの低減効果である。メタルアーティファクトとは、観察対象中にＸ線に対して不透明な部位（多くの場合、金属部分）が存在する場合、ＣＴで得られる断面像全体（不透明部位だけでない）が破壊的に乱されることを指す。メタルアーティファクトの原因は、不透明部分で投影像の輝度が不連続に変化すること、不透明部分の情報が欠落することである。ＦＢＰにおいて微分フィルタが適用されると、輝度の不連続な変化は異常値となり、逆投影操作によって、不透明部位を中心とした筋状のアーティファクトを生み出す。さらに、情報の欠落によって不透明部位とは直接関係のない部分に、予期しない明暗が生じる。ＦＢＰを用いず、情報の欠落に対して安定なＳＡを利用した本発明がメタルアーティファクトの除去に有効であろうことは想像に難くない。
実用上、重要なのは、コーンビームやヘリカルスキャンへの応用であろう。両者は３次元ＣＴと呼ばれ、現在急速に普及しつつある技術である。しかしながら、３次元ＣＴには特有のアーティファクトが現れることが知られており、その原因はわかっているのだが解決には至っていない。その原因は、コーンビームの場合には、データの欠落である。コーンビームでは完全な断面像を得るための条件を満たすことができず、データの欠落に由来するアーティファクトが現れる。ヘリカルスキャンでのアーティファクトの根本的な原因は、逆投影操作にある。ヘリカルスキャンでは系の幾何学的異方性（らせん）がフィルタ操作・逆投影操作に影響し、風車アーティファクトが現れる。本発明は、フィルタ操作・逆投影操作を必要とせず、データの欠落にも強いので、３次元ＣＴにおける諸問題を解決できる。
投影角度が不均一である場合として挙げられるのは、地震波による地球内部のＣＴや、人工衛星からの電波を用いた大気状態のＣＴによる解析である。これらはＦＢＰが利用できない典型として知られている。本発明を利用することで、解析精度の向上を期待できる。
本発明の第二の効果は、投影像撮影の迅速化、Ｘ線被爆量の低減である。ＳＡはデータの欠落だけでなく、過少にも安定であるという特徴があり、本発明も同様である。ＣＴの場合、データの過少に相当するのは、投影角度の数が少ない場合である。すなわち、本発明を利用すると、従来のＣＴより投影角度の数を少なくすることができる。投影角度の数とは放射線による投影像の撮影枚数である。撮影枚数は、撮影時間と被爆量に比例するので、撮影枚数が少なければ、撮影時間の短縮、被爆量の低減につながる。
また、データの過少の別のパターンとして、投影像の画質が不良（Ｓ／Ｎ比が低い）な場合が挙げられる。本発明はそのような場合にも比較的安定であることがわかっている。投影像の画質の低下を許容できるとなると、これも撮影時間の短縮・被爆量の低減につながる。あるいは、Ｓ／Ｎ比が極端に低いＳＰＥＣＴ，ＰＥＴにおいて本発明は画質の向上に寄与できるだろう。
本発明の第三の効果は、断面像の輝度値の決定精度が高いことである。本発明では測定した投影像をかなり忠実に再現する断面像が得られる。再現の精度はＦＢＰより２桁程度高い。これは二乗誤差を最小化するというフィッティングアルゴリズムの恩恵である。輝度値の決定精度が高いことで、輝度値の定量性が保証され、輝度値を用いた密度測定が可能になる。これは骨密度測定などにおける測定精度の向上につながる。また、異常のある部位（腫瘍のある臓器など）の検出精度の向上にも寄与するだろう。
本発明の第四の効果は、本発明で得られる断面像のコントラストがＦＢＰよりも高いことである。第三の効果で述べたように、本発明は輝度値の決定精度が高いのであるが、その副次的な効果として、断面像のコントラストが高くなるのである。コントラストが高いと見掛けの空間分解能も高くなる傾向がある。その結果、本発明で得られる断面像は、従来のものより高画質になる。この効果は、特殊な条件・用途のＣＴではなく、普通のＣＴにも本発明が役立つという点で、特筆すべきものである。
図１は、Ｘ線ＣＴ装置の模式図である。
図２は、断面像ｆ（ｘ，ｙ）、投影像ｐ（ｒ，θ）の関係を示す模式図である。
図３は、典型的なｐ０（ｒ，θ）の例（横軸に検出器のチャンネル位置、縦軸に投影角度をとった画像）である。
図４は、本発明の基本手順を示すフローチャートである。
図５（ａ）は、手順（Ｉ）〜（ＶＩ）のフィッティング経過の模式図である。図５（ｂ）は、手順（１）〜（９）のフィッティング経過の模式図である。
図６は、［非特許文献１］に掲載されている手法での結果を示す図である。
図７は、図６（ａ）からシミュレーションによって求めたｓｉｎｏｇｒａｍである（白い部分は不可視領域）。
図８（ａ）は、図６（ｃ）をさらに拡大した図である。図８（ｂ）は、本発明での結果を示す図である。
図９は、仮想エネルギーＥの各項の効果を表した模式図である。
図１０は、角度制限がある場合のアーティファクトと本発明でのアーティファクト低減効果を示す図である。
（第１の実施形態）
本発明の実施形態によるＸ線ＣＴ装置の概略構成を図１に示す。このＸ線ＣＴ装置は、撮影部とコンピュータとを備えている。撮影部は、Ｘ線を測定対象に照射するための光源と、測定対象を透過したＸ線を検出する検出器とを備えており、測定対象の周囲を多数の方向からＸ線により投影して投影像を取得する。コンピュータは、Ｘ線ＣＴ装置全体を制御する制御部と、撮影部により得られたＸ線投影像に基づいて測定対象の関心領域の断面像を生成する画像再構成処理部とを備えている。なお、図１に示した構成はこれ以降の各実施形態に共通のものであり、各実施形態では画像再構成処理部において行われる処理がそれぞれ異なっている。
本実施形態の画像再構成処理部は、投影像から断面像を求めるためのフィッティング方法としてＳｉｍｕｌａｔｅｄ Ａｎｎｅａｌｉｎｇ法（ＳＡ）を採用している。そこで、まず、Ｓｉｍｕｌａｔｅｄ Ａｎｎｅａｌｉｎｇ法（ＳＡ）のフレームワークを説明する。ＳＡはＭｏｎｔｅ Ｃａｒｌｏ法から派生したフィッティングアルゴリズムで、乱数に基づいてフィッティングを進める点と、熱力学を模して仮想エネルギーや仮想温度を取り扱う点、が特徴となっている。ＳＡ自体は公知の技術であり以下のステップ（ｉ）〜（ｖｉ）により行われる。
（ｉ）パラメータを乱数に基づいて１つ選択し、さらに乱数に基づいてそのパラメータを変更する（乱数）。
（ｉｉ）変更後の評価を行う。評価関数として仮想的なエネルギーＥを考える。通常のＳＡではＥは二乗誤差の総和である。変更前後のＥの変化をΔＥとする（評価）。
（ｉｉｉ）評価が改善される（ΔＥ＜０）なら、その変更を受け入れる（変更）。
（ｉｖ）ｅｘｐ（−ΔＥ／Ｔ）の確率で改悪を受け入れる。
（ｖ）Ｔを少しずつ減じる。
（ｖｉ）（ｉ）から繰り返す。
ＳＡでは、（ｉｖ）にあるように、ボルツマン統計に従って改悪を受け入れることで、局所解から脱出できる可能性を担保している。そのため、局所解にとらわれることなく、大域解に辿りつくことができ、条件の悪いフィッティングに対しても安定に機能する。また、Ｔを徐々に減じることで、大域解にソフトランディングするようになっている。（ｉ）から（ｖ）まで一通り実行することを１モンテカルロステップと呼ぶ。ＳＡではモンテカルロステップを無数に繰り返すことで、フィッティングが進行する。計算に必要な時間は、１モンテカルロステップに必要な時間×必要なモンテカルロステップ数である。必要なモンテカルロステップ数は、パラメータの数・自由度に比例する。
次に、請求項に沿った形で本実施形態の説明を行う。ＣＴをフィッティング問題と考えた場合、フィッティングパラメータは断面像（ｆ（ｘ，ｙ））である。フィッティング対象となるデータは、測定した投影像（ｐ０（ｒ，θ））である。ｒは投影像を撮影した１次元検出器のチャンネル位置、θは投影角度である。座標の定義を図３に示す。ｐ０（ｒ，θ）は、角度θを変えながら投影像を測定することで得られるデータセットであり、横軸にｒ、縦軸にθを取ると、２次元の画像とみなすことができる。このようなデータをｓｉｎｏｇｒａｍと呼ぶ。典型的なｓｉｎｏｇｒａｍを図３に示す。
極論すれば、ＣＴはｓｉｎｏｇｒａｍから断面像を求める技術である。本実施形態では、暫定の断面像ｆ（ｘ，ｙ）と測定した投影像ｐ０（ｒ，θ）との二乗誤差を考えるが、二乗誤差を計算するために、数１０１を使って、ｆ（ｘ，ｙ）から計算で投影像（ｐ（ｒ，θ））を求める。
数１０１は図２において、ｓの方向に和をとっており、ｆ（ｘ，ｙ）のｓ方向の投影像の計算であることがわかる。こうして求めたｐ（ｒ，θ）を用いて、二乗誤差Ｈは以下のように計算できる。
通常のフィッティングでは、仮想エネルギーＥとして、数１０２をそのまま利用するが、本実施形態では、Ｈ以外にスムージング項とエントロピー項を導入している点が特徴となっている。Ｅの定義は次の通りである。
Ｔは仮想温度（温度パラメータ）、Ｓはエントロピー、σは画素値の標準偏差、ｃはスムージング項の強さを表す係数である。ＴＳがエントロピー項、ｃσがスムージング項である。Ｓやσの定義・計算法は後述する。本実施形態では、これらの数式を使い、図４に示すような手順で断面像を計算する。
・ステップ（ａ）：何らかの形で求めた暫定の断面像ｆ（ｘ，ｙ）に対し、評価関数（以下「エネルギー」と呼ぶ）（Ｅ０）を求める（ＳＴ１０およびＳＴ２０）。
・ステップ（ｂ）：乱数で（ｘ０，ｙ０）およびΔμを選択し、断面像ｆ（ｘ，ｙ）の一部を改変する（ＳＴ３０）。
・ステップ（ｃ）：改変後の断面像ｆ（ｘ，ｙ）に対して、評価関数Ｅ１を求める（ＳＴ４０およびＳＴ５０）。
・ステップ（ｄ）：前記エネルギー（Ｅ０）と前記エネルギー（Ｅ１）との差（ΔＥ）を求める（ＳＴ６０）。
・ステップ（ｅ）：前記改変を採用するか否かを前記差（ΔＥ）および温度パラメータ（Ｔ）を用いた受理関数（典型的にはボルツマン統計）により判定する（ＳＴ７０〜ＳＴ１１０）。
・ステップ（ｆ）：前記ステップ（ａ）に戻る（ＳＴ１２０）。
・ステップ（ｇ）：前記ステップ（ａ）〜（ｆ）の繰り返しが所定回数に達するごとに仮想温度（Ｔ）の値を変更する（ＳＴ１３０）。
・ステップ（ｈ）：前記ステップ（ｅ）における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、ＳＴ８０，ＳＴ１００が実行された場合を「成功」、ＳＴ１１０が実行された場合を「失敗」とし、成功の確率が１０％（この値は適宜調整可能）を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
ステップ（ａ）からステップ（ｈ）の一連の操作が請求項１２に対応し、本実施形態ではこれら一連の操作を図１の画像再構成処理部が行う。
ＳＡの基本手順（ｉ）〜（ｖｉ）との対応を考えると、ステップ（ｂ）が（ｉ）、ステップ（ａ）（ｃ）（ｄ）が（ｉｉ）、ステップ（ｅ）に（ｉｉｉ）と（ｉｖ）が含まれ、ステップ（ｇ）が（ｖ）、ステップ（ｆ）が（ｖｉ）にそれぞれ対応する。このことから、本実施形態はＣＴにＳＡを忠実に当てはめたものであることがわかる。
なお、上記ステップ（ｈ）における処理終了の判定ついては、図１の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
（第２の実施形態）
数１０３に示されたＥを計算するためには、数１０１と数１０２を経ることになり、ｓ，ｒ，θに関して級数和を実行することになる。３重の積分（級数和）を計算するので、大変時間がかかる。すなわち、１モンテカルロステップに必要な計算時間が長いことを意味する。さらに、ＣＴではパラメータの数が莫大であり、その結果、ＳＡの実施に必要な総計算時間は、現在の計算機でも「年単位」になってしまう。
そこで、本実施形態では、Ｅを計算するのではなく、変更を行ったときの変化分ΔＥだけを主に計算する。
今、暫定断面像ｆ（ｘ，ｙ）に関する変更をΔｆ（ｘ，ｙ）とする。Δｆ（ｘ，ｙ）は座標（ｘ０，ｙ０）においてのみΔμの値を持ち、それ以外は０であるような断面像である。そのΔｆ（ｘ，ｙ）の投影像Δｐ（ｒ，θ）は数１０１と同じ方法で計算できる。このΔｐ（ｒ，θ）を用いて、ΔＨは以下のように計算できる。
これを書き下して、
となる。今、Δｆ（ｘ，ｙ）は（ｘ０，ｙ０）においてのみ値を持つので、Δｐ（ｒ，θ）はｒ（θ）＝ｘ０ｃｏｓθ＋ｙ０ｓｉｎθにおいてのみ値Δμを持ち、それ以外は０である。従って、数１０６の｛｝の中身も、ｒ（θ）＝ｘ０ｃｏｓθ＋ｙ０ｓｉｎθにおいてのみ値を持つ。そのため、級数和はｒとθの両方に対して行う必要はなく、以下のように表すことができる。
数１０７ではθに関する級数和だけになっている点が重要である。これにより、計算量を大幅に削減することができる。ｐやｐ０はデジタル画像であるため、数１０７の計算にあたっては、通常は、ｒ（θ）に関する内挿を行う必要がある。そのため、数１０７の｛｝内をさらに展開することはできない。しかしながら、誤差を容認して数１０７を展開すると次式が得られる。
ただし、Ｍは投影角度の数である。数１０７の代わりに数１０８を用いるとΔＨの値に１％程度の誤差が生じる。しかしながら、数１０８は数１０７より高速であり、本発明での利用価値は高い。
次にΔσの計算を説明する。σは座標（ｘ０，ｙ０）付近の輝度値の標準偏差である。（ｘ０，ｙ０）付近として、（ｘ０，ｙ０）の周囲ｄ×ｄ画素を考える（実施例ではｄ＝５を用いている）。それらの画素の標準偏差は次式で求めることができる。
ただし、
である。数１０９〜１１１を用いると、Δσは次のように計算することができる。
ただし、ｆｉ，ｆｊは変更前および後のｆ（ｘ０，ｙ０）の値である。
次にΔＳの計算であるが、その前にエントロピーＳの定義を行う。一般に計算機が取り扱う画像はデジタル画像であり、画素の座標（ｘ，ｙ）だけでなく、画素の値もデジタルの値となっている。そこで、１つの画素を１つの量子と考え、画素値を量子状態とみなすことにする。すると画像は複数の量子からなる系であると考えることができる。統計力学に則ると、そのような系に対するエントロピーは次のように定義される。
Ｎは画像中の画素の総数であり、Ｎｉは画素値がデジタル値でｉである画素の総数である。ｋは通常の物理ではボルツマン定数であるが、本発明は物理とは直接関係がないため任意である。本発明では、σと同様にＳも（ｘ０，ｙ０）の周囲ｄ×ｄ画素で定義されるとした。
数１１３で定義されたＳに対してΔＳを考える。変更によって、画素値はデジタル値でｉからｊに変更されるとする。すると、ΔＳは次式のようになる。
すなわち、変更によって、Ｎｉが１減って、Ｎｊが１増えるのである。数１１４を展開すると、非常に簡単な式が得られる。
以上をまとめて、本実施形態の手順を書き下すと次のようになる。
（Ｉ）乱数を使って（ｘ０，ｙ０）およびΔμを選択する。
（ＩＩ）数１０４、１０８（または１０７）、１１２、１１５を使ってΔＥを計算する。
（ＩＩＩ）評価が改善される（ΔＥ＜０）なら、ｆ（ｘ０，ｙ０）にΔμを足す。
（ＩＶ）ΔＥ＞０の時にも、ｅｘｐ（−ΔＥ／Ｔ）の確率で、ｆ（ｘ０，ｙ０）にΔμを足す。
（Ｖ）Ｔを少しずつ減じる。
（ＶＩ）（Ｉ）から繰り返す。
（ＶＩＩ）所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。
ここでは、（ＩＩＩ），（ＩＶ）においてΔμが足された場合を「成功」、それ以外を「失敗」とし、成功の確率が所定値（例えば１０％、この値は適宜調整可能）を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
手順（Ｉ）〜（ＶＩ）はＳＡの基本手順（ｉ）〜（ｖｉ）を忠実に踏襲していることがわかる。なお、手順（Ｉ）〜（ＶＩＩ）の処理は図１の画像再構成処理部において行われる。
なお、上記手順（ＶＩＩ）における処理終了の判定ついては、図１の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
（第３の実施形態）
次に、上で述べた方法の重要な変形について述べる。上で述べた方法では数１０７あるいは数１０８を利用し、θに関する級数和を求めている。そのため、画素１つの変更を検討するのにＭ回の繰り返し演算を要する。本実施形態はこの部分の計算量をさらに減らすものである。
まず、ｐ（ｒ，θ）の逆投影像ｇ（ｘ，ｙ）を考える。
数１１６では、フィルタ操作をしていないため、ｇ（ｘ，ｙ）は断面像ｆ（ｘ，ｙ）に似ているが、ぼやけたような画像になる。ｇ（ｘ，ｙ）を用いると、数１０８は次のようになる。
ただし、ｇ０（ｘ，ｙ）はｐ０（ｒ，θ）の逆投影像で、数１１６と同様の方法で計算される。数１１７では繰り返し演算がないという点で数１０８より優れている。ｇ０（ｘ，ｙ）は計算過程で全く変化しないため、予め計算しておくことができる。一方、ｇ（ｘ，ｙ）はｆ（ｘ，ｙ）の１点を変更する度にわずかずつ変化するため、厳密には（Ｉ）〜（ＩＶ）の手順を１回経るごとに再計算が必要である。しかしながら、ｆ（ｘ，ｙ）の変更に伴うｇ（ｘ，ｙ）の変化は少ないという仮定をすると、別の手順を適用することができ、それが本実施形態である。以下に手順を示す。なお、以下の手順（１）〜（９）の処理は図１の画像再構成処理部において行われる。
（１）ｐ０（ｒ，θ）からｇ０（ｘ，ｙ）を求める。
（２）ｆ（ｘ，ｙ）からｐ（ｒ，θ）を求め、さらにｐ（ｒ，θ）からｇ（ｘ，ｙ）を求める。
（３）ｆ（ｘ，ｙ）に対する変更値を画素値とする画像Δμ（ｘ，ｙ）を乱数で生成する。
（４）各画素値について数１１７を適用し、ΔＨ（ｘ，ｙ）という画像を計算する。
（５）同様に各画素値について数１１２、数１１５を適用し、Δσ（ｘ，ｙ）、ΔＳ（ｘ，ｙ）という画像を計算する。
（６）数１０４に基づいてΔＥ（ｘ，ｙ）を計算する。
（７）ΔＥが正である座標（ｘ，ｙ）について、Δμ（ｘ，ｙ）を０にする。
（８）ｆ（ｘ，ｙ）にΔμ（ｘ，ｙ）を足す。
（９）Ｔをα倍（α＜１）し、（２）から繰り返す。
（１０）所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、（７）においてΔμ（ｘ，ｙ）が０にされた場合を「失敗」、されなかった場合を「成功」とし、成功の確率が所定値（例えば１０％、この値は適宜調整可能）を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
なお、上記手順（１０）における処理終了の判定ついては、図１の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
手順（Ｉ）〜（ＶＩＩ）と、手順（１）〜（１０）を図５に模式的に示した。図５ではパラメータの値（本件では画素の値）を軸とし、仮想エネルギーＥを濃淡で表している（Ｅが小さいほど濃い）。図示の都合により、パラメータの数（画素の数）は２個の場合を取り扱っている。図５（ａ）に示す手順（Ｉ）〜（ＶＩＩ）では、ジグザグに折れ曲がりながら、時にでたらめな方向に進みながら、最終的にはＥの最小値を探索してゆく。一方、図５（ｂ）に示した手順（１）〜（１０）では、Ｅの勾配を考慮しながら、Ｅの最小値を探索してゆく。図５（ｂ）の方が効率が良いということが直感的に理解できるだろう。数１１７による計算量の低減及び、図５（ｂ）に示した探索効率の向上により、本実施形態は計算速度的に有利である。手順（１）〜（１０）はＳＡの基本アルゴリズム（ｉ）〜（ｖｉ）を忠実に踏襲するものではないが、（Ｉ）〜（ＶＩＩ）の手順の自然な発展形である。
（その他の実施形態）
上述の各実施形態で説明した処理を実行する画像再構成処理部は、これらの処理をコンピュータに実行させるためのプログラム、当該プログラムがインストールされたコンピュータ、さらには、これらの処理を実行する専用ＬＳＩなどを使って実現可能である。
実施例としては、シミュレーションを用いたメタルアーティファクト除去効果を取り扱う。まず、比較のために、［非特許文献１］に掲載されている手法での結果を図６（非特許文献のｆｉｇ．５）に示す。
図６において、左のカラム（ａ）（ｄ）（ｇ）はＭｅｔａｌ Ａｒｔｉｆａｃｔがないオリジナルのファントム画像である。中央（ｂ）（ｅ）（ｈ）は、オリジナルのファントムの所定の位置に不可視領域を設定し、通常のＦＢＰ法により再構成した画像である。右（ｃ）（ｆ）（ｉ）は、上記と同様に不可視領域を設定し［非特許文献１］の手法により再構成した画像であり、従来技術で最もうまくＭｅｔａｌ Ａｒｔｉｆａｃｔが除去できている例の１つである。
そこで、［非特許文献１］から図６（ａ）の画像を取り込み、文献と同じ位置に不可視領域を設定し、シミュレーションで投影像を計算した。結果を図７に示す。図７が本実施例におけるｐ０（ｒ，θ）に対応する。図７について、本発明によって再構成した結果および、比較のための図６（ｃ）の拡大図を図８に示した。効果は一目瞭然である。本実施形態（ｂ）ではＭｅｔａｌ Ａｒｔｉｆａｃｔの特徴を見つけ出すことすら困難である。
ただし本実施形態（ｂ）ではエッジ付近が若干ボケている。これは「断面像を滑らかにする因子」（スムージング項ｃσ）がやや強いためである。Ｍｅｔａｌ Ａｒｔｉｆａｃｔを低減するには、その因子をやや強めに設定する必要があることがわかっている。このように本実施形態によれば、仮想エネルギーＥを計算する際の各因子のバランスには改良の余地があるものの、Ｍｅｔａｌ Ａｒｔｉｆａｃｔ除去の効果は非常に高い。
なお、参考までに、スムージング項ｃσ、エントロピー項ＴＳを設定せずに本実施形態のアルゴリズムを実行した場合の結果を図９に示す。図９（ａ）はスムージング項ｃσなし（エントロピー項ＴＳのみ）の場合の結果、図９（ｂ）はエントロピー項ＴＳなし（スムージング項ｃσのみ）の結果、図９（ｃ）はエントロピー項ＴＳ、スムージング項ｃσともになしの場合の結果を示す。図９と図８（ｂ）とを比較すると、本実施形態においては、エネルギーＥの因子としてスムージング項ｃσとエントロピー項ＴＳの両方を設定したほうが、いずれか一方のみを設定した場合、ともに設定しない場合よりもＭｅｔａｌ Ａｒｔｉｆａｃｔ除去の効果が高いことが分かる。スムージング項だけでは、不透明を中心に放射状に広がる縞模様がアーティファクトして残るという特徴がある。一方、エントロピー項のみでは、断面像がざらつきＳ／Ｎが低いという特徴がある。また、最終的な画質は、徐冷過程におけるスムージング項の係数ｃと仮想温度Ｔの比が重要であるようである。
次に、回転角度制限がある場合の実施例を図１０に示す。図１０（ａ）は、ＦＢＰのみ（通常のＣＴ）による再構成結果を示す。図１０（ｂ）は、ＩＬＳＴ（最小二乗法に基づく画像復元法）を適用した結果を示す。図１０（ｃ）は、本実施形態による再構成結果を示す。これらを比較すると、角度制限時に現れるアーティファクトに対して本実施形態が有効であることが分かる。角度制限によるアーティファクトには、円状の領域が両端の尖ったアーモンド状の領域になり、アーモンド状領域の尖っていない胴回りの輝度が反転するという特徴がある。本実施形態によると、それら両方の特徴が低減される。
本発明は、Ｍｅｔａｌ Ａｒｔｉｆａｃｔに対して顕著な効果があるので、Ｍｅｔａｌ Ａｒｔｉｆａｃｔが深刻な分野、たとえば歯のＣＴ、金属製インプラントを含むＣＴに特に有用である。
また、本発明は情報欠落のある投影像のセットに対して一般的に有効である。情報欠落が顕著な場合の一つに、投影角度制限がある場合を挙げることができる。投影角度制限が問題になる分野は、三次元電子顕微鏡法、マンモグラフィーのＣＴ版、平行移動ＣＴ法（特開２００６−７１４７２号公報）などである。情報欠落が問題になる別の場合として、コーンビームＣＴ、ヘリカルスキャンＣＴなど３次元ＣＴが挙げられる。本発明は、３次元ＣＴに現れるアーティファクトを除去する目的にも有効である。
さらには、情報が極端に少ない系に関しても利用可能である。例えば、蛍光Ｘ線ＣＴ、地震波による地球内部のＣＴなどに有用である。
また、本発明の副次的な効果として、断面像における輝度値（Ｘ線では線吸収係数に対応する）を従来法よりも高い精度で決定できることが挙げられる。この効果は骨密度測定等における精度改善に応用可能である。
本発明を用いると、従来の方法よりコントラストの高い再構成像を得ることができる。そのため、アーティファクト等が問題とならない通常のＸ線ＣＴにおいても、本発明の利用価値は高い。また、本発明はデータの過少に対しても安定であることから、投影像測定の時間短縮、ひいてはＸ線被爆量の低減に有効である。これらをあわせて考えると、本発明は、既存のすべてのＸ線ＣＴ技術を置き換える可能性を秘めている。
本発明は、対象物の放射線投影像から断面像を再構成する技術に関する。
計算機トモグラフィー(ＣＴ)は、Ｘ線などの光線を使って、様々な角度からの物体の投影像を撮影し、その後、計算によって断面像を得る技術である。図１に典型的なＸ線ＣＴ装置の模式図を示す。
典型的なＣＴ装置では、Ｘ線光源が回転移動することにより、測定対象の様々な角度からの投影像を取得する仕組みになっている。こうして得られた投影像から計算操作によって断面像を得るのがＣＴである。通常は、Filtered Back-projection法(FBP)によって投影像から断面像を変換操作によって求める。FBPでは、投影像にフィルタ(基本的には微分フィルタ)を施した後、もともとの投影方向に「逆投影」という操作を施すことで断面像を得る。このとき、微分フィルタはノイズや誤差を増幅し、アーティファクト(本来存在しない誤差や虚像)を生みやすいという特徴がある。さらに、逆投影操作はそのアーティファクトを断面像全体に伝播させる効果がある。そのため、ＣＴにおけるアーティファクトは問題のある部分だけに留まらず、断面像全体を損なうという点で致命的になることが多い。
アーティファクトの多くは、FBPに含まれるフィルタ操作や逆投影操作が原因なのだから、FBPを用いなければ、アーティファクトのない断面像が得られるはずである。FBP以外で断面像を計算する方法として、Algebraic Reconstruction Technique(ART)が歴史的に重要である。ARTはFBPが考案される前から利用されるほど歴史のある方法である。ARTでは断面像の計算をフィッティング問題と捉え、断面像をパラメータ、投影像をフィッティング対象として、断面像から計算した投影像(ｐ)が、実験で求めた投影像(ｐ０)に一致するように逐次的に断面像を修正する。ARTでは(ｐ−ｐ０)が０になるように断面像を漸近的に改良してゆく点が特徴となっている。ARTはFBPを全く用いないかわりに、断面像を求める計算に時間がかかるため、今では特殊な用途(地震波の解析など)にしか用いられない。ARTはFBPほど極端なアーティファクトは現れないが、得られる断面像はFBPの方が自然であることが多い。
さらに、フィルタ操作や逆操作に由来しないアーティファクトの要因として、投影像におけるデータの欠落・過少も挙げられる。データが少なければ結果として得られる断面像の画質も低下するのは当たり前である。FBPではデータの欠落・過少も致命的なアーティファクトを生み出すことが知られており、大きな問題と認識されている。フィッティングに基づくＣＴ(ARTを含む)はデータの欠落や過少に対してFBPよりは強いとされている。しかしながら、ＣＴにおけるデータの欠落は、極端に「条件が悪い」問題とされており、ARTなどを用いても改善が得られにくい。ARTでは(ｐ−ｐ０)を誤差として考えるのであるが、条件が悪いときはフィッティングに失敗しやすい。いわゆる二乗誤差(ｐ−ｐ０)２の最小化を行った方がフィッティングとしては安定する。二乗誤差を最小化する方法としては、最小二乗法が最も標準的である。最小二乗法では、一辺がパラメータの個数だけある正方行列の逆行列を求める。ＣＴにおけるパラメータは断面像を構成するそれぞれの画素値であるため、パラメータの個数は莫大になる。1000×1000画素の断面像だと、パラメータの数は100万個になり、行列の要素数は1兆個と莫大である。従って、正攻法の最小二乗法では行列が大きくなりすぎで破綻してしまう。そこで、正攻法の最小二乗法を避け、Simultaneous iterative reconstruction technique (SIRT)やIterative least squaretechnique (ILST)が考案された。ARTと同様に断面像の計算をフィッティング問題として捉えるのであるが、先に述べたように、正攻法の最小二乗法を避けなければならないため、SIRTもILSTも計算の途中に投影操作の逆演算としてFBPが利用される。そのため、FBPに起因するフィルタ操作と逆投影にまつわる問題を根本的には解決しない。おそらくそれが理由となり、SIRTやILSTでは各種アーティファクトの「低減」効果が報告されているに留まる。
Yazdi M, Gingras L, BeaulieuL. An adaptive approach to metal artifact reduction in helical computedtomography for radiation therapy treatment planning: experimental and clinicalstudies. Int J Radiat Oncol Biol Phys,62: 1224-1231, 2005.
以上の議論から、アーティファクトの無い断面像を得るには、FBPを用いない、データの欠落に対して強い、という性質が重要であることがわかる。そこで、ARTと同様にFBPを用いず、フィッティングの評価関数として二乗誤差を利用しようというのは、誰でも思いつくことである。にもかかわらず、ＣＴにおけるアーティファクトの問題が根本的に解決していないのは、この単純な考えの実現は簡単ではないからである。第一の問題は、断面像の計算はフィッティングの問題としては大規模すぎて、計算に時間がかかることである。適切な高速化が必須である。第二の問題は、最小二乗法はフィッティングアルゴリズムとしては「条件が悪い」問題に対して弱いことである。そのため、最小二乗法に基づく既存のアルゴリズム(ILSTなど)の延長では、アーティファクトの問題を解決できないと思われる。第三の問題は、SIRTやILSTがそうであるように、既存のアルゴリズムではFBPを完全に排除することができていないことである。
本発明の第一の新規性は、フィッティング方法としてSimulated Annealing法(ＳＡ)を採用したことである。ＳＡは公知の技術であるが、フィッティングに時間がかかることが知られており、ＣＴと本質的に相性が悪い。しかしながら、敢えてＳＡをＣＴに適用することで、FBPを全く用いずとも、二乗誤差の最小化による断面像の計算ができるようになった。また、ＳＡはデータの欠落など条件の悪いフィッティングに対してかなり安定であるという性質があり、この点においてもアーティファクトの根本的な解決に有利である。以上の事柄を鑑み、ＳＡをＣＴに適用することで、アーティファクトの根本的解決に効果を見出した点が、本発明の重要な新規性である。単純にＳＡのＣＴへの適用を考えたとき、断面像から投影像を求める計算を無数に繰り返すことになる。まともに計算すると、断面像から投影像を求める計算は、FBPとほぼ同じであり、ＳＡをＣＴに適用すると、FBPに比べて数百万倍の時間がかかってしまうことになる。現在の高速な計算機を以ってしても、「年」単位の計算時間を要してしまう。そこで、本発明では、式変形によって計算量を大幅に低減することで、この問題を解決した。この工夫は、ＳＡをＣＴに適用する際に必要な技術であり、本発明の重要な要素である。
第二の新規性は、フィッティングの評価関数として、二乗誤差だけでなく、アーティファクトを積極的に破壊するスムージング項とエントロピー項を導入した点である。従来のＣＴでは、ｐとｐ０との差あるいは二乗誤差だけを対象としていた。この場合、アーティファクトを御座なりにした方が良いフィッティング結果(小さな誤差)が達成できるという可能性が常にあり、現実の多くのケースではそのような条件になっていた。すなわち、アーティファクトを別のアーティファクトで相殺する状態になっていたのである。そのため、アーティファクトはなかなか消えないのである。本発明においても二乗誤差だけを評価関数とした場合には、アーティファクトは完全には消えなかった。この事実は、アーティファクトのない断面像を得るには、二乗誤差以外の要請が望ましいことを意味している。本発明では、統計力学の基本にのっとり、エントロピー項やスムージング項を考案した。これらの項は、断面像は滑らかで自然な濃淡画像であるべきだ、という当たり前の要請を数式の形で表現したものである。エントロピー項は、アーティファクトを破壊し、断面像全体の画質を均質にするように、フィッティングを誘導する。スムージング項は、エントロピー項で誘導される断面像のザラツキを押さえる効果がある。これらの項を導入することで、アーティファクトが消え、かつ自然な断面像が得られるようになった。エントロピー項とスムージング項は、それぞれ単独でもアーティファクトを低減する効果があるが、２つを組み合わせるとさらに効果的である。
なお、本発明においては、一般の定義よりも広い概念として、X線・可視光・電波を含む電磁波、電子や荷電粒子からなる粒子線、媒体の振動である音波等を全て含めて「放射線」と呼ぶことにする。
本発明の第一の効果として、データの欠落が要因になっているアーティファクトの低減効果が挙げられる。データの欠落が問題になる場合として、観察対象に不透明部位がある場合、投影角度に制限がある場合、投影角度の刻みが不均一である場合、コーンビームやヘリカルスキャンなどの三次元ＣＴの場合が挙げられる。
特に効果があることがわかっているのは、観察対象に不透明部位がある場合に現れるメタルアーティファクトの低減効果である。メタルアーティファクトとは、観察対象中にＸ線に対して不透明な部位(多くの場合、金属部分)が存在する場合、ＣＴで得られる断面像全体(不透明部位だけでない)が破壊的に乱されることを指す。メタルアーティファクトの原因は、不透明部分で投影像の輝度が不連続に変化すること、不透明部分の情報が欠落することである。FBPにおいて微分フィルタが適用されると、輝度の不連続な変化は異常値となり、逆投影操作によって、不透明部位を中心とした筋状のアーティファクトを生み出す。さらに、情報の欠落によって不透明部位とは直接関係のない部分に、予期しない明暗が生じる。FBPを用いず、情報の欠落に対して安定なＳＡを利用した本発明がメタルアーティファクトの除去に有効であろうことは想像に難くない。
実用上、重要なのは、コーンビームやヘリカルスキャンへの応用であろう。両者は３次元ＣＴと呼ばれ、現在急速に普及しつつある技術である。しかしながら、３次元ＣＴには特有のアーティファクトが現れることが知られており、その原因はわかっているのだが解決には至っていない。その原因は、コーンビームの場合には、データの欠落である。コーンビームでは完全な断面像を得るための条件を満たすことができず、データの欠落に由来するアーティファクトが現れる。ヘリカルスキャンでのアーティファクトの根本的な原因は、逆投影操作にある。ヘリカルスキャンでは系の幾何学的異方性(らせん)がフィルタ操作・逆投影操作に影響し、風車アーティファクトが現れる。本発明は、フィルタ操作・逆投影操作を必要とせず、データの欠落にも強いので、３次元ＣＴにおける諸問題を解決できる。
投影角度が不均一である場合として挙げられるのは、地震波による地球内部のＣＴや、人工衛星からの電波を用いた大気状態のＣＴによる解析である。これらはFBPが利用できない典型として知られている。本発明を利用することで、解析精度の向上を期待できる。
本発明の第二の効果は、投影像撮影の迅速化、Ｘ線被爆量の低減である。ＳＡはデータの欠落だけでなく、過少にも安定であるという特徴があり、本発明も同様である。ＣＴの場合、データの過少に相当するのは、投影角度の数が少ない場合である。すなわち、本発明を利用すると、従来のＣＴより投影角度の数を少なくすることができる。投影角度の数とは放射線による投影像の撮影枚数である。撮影枚数は、撮影時間と被爆量に比例するので、撮影枚数が少なければ、撮影時間の短縮、被爆量の低減につながる。
また、データの過少の別のパターンとして、投影像の画質が不良(Ｓ/Ｎ比が低い)な場合が挙げられる。本発明はそのような場合にも比較的安定であることがわかっている。投影像の画質の低下を許容できるとなると、これも撮影時間の短縮・被爆量の低減につながる。あるいは、Ｓ/Ｎ比が極端に低いSPECT，PETにおいて本発明は画質の向上に寄与できるだろう。
本発明の第三の効果は、断面像の輝度値の決定精度が高いことである。本発明では測定した投影像をかなり忠実に再現する断面像が得られる。再現の精度はFBPより２桁程度高い。これは二乗誤差を最小化するというフィッティングアルゴリズムの恩恵である。輝度値の決定精度が高いことで、輝度値の定量性が保証され、輝度値を用いた密度測定が可能になる。これは骨密度測定などにおける測定精度の向上につながる。また、異常のある部位（腫瘍のある臓器など）の検出精度の向上にも寄与するだろう。
本発明の第四の効果は、本発明で得られる断面像のコントラストがFBPよりも高いことである。第三の効果で述べたように、本発明は輝度値の決定精度が高いのであるが、その副次的な効果として、断面像のコントラストが高くなるのである。コントラストが高いと見掛けの空間分解能も高くなる傾向がある。その結果、本発明で得られる断面像は、従来のものより高画質になる。この効果は、特殊な条件・用途のＣＴではなく、普通のＣＴにも本発明が役立つという点で、特筆すべきものである。
（第１の実施形態）
本発明の実施形態によるＸ線ＣＴ装置の概略構成を図１に示す。このＸ線ＣＴ装置は、撮影部とコンピュータとを備えている。撮影部は、Ｘ線を測定対象に照射するための光源と、測定対象を透過したＸ線を検出する検出器とを備えており、測定対象の周囲を多数の方向からＸ線により投影して投影像を取得する。コンピュータは、Ｘ線ＣＴ装置全体を制御する制御部と、撮影部により得られたＸ線投影像に基づいて測定対象の関心領域の断面像を生成する画像再構成処理部とを備えている。なお、図１に示した構成はこれ以降の各実施形態に共通のものであり、各実施形態では画像再構成処理部において行われる処理がそれぞれ異なっている。
本実施形態の画像再構成処理部は、投影像から断面像を求めるためのフィッティング方法としてSimulated Annealing法(ＳＡ)を採用している。そこで、まず、Simulated Annealing法(ＳＡ)のフレームワークを説明する。ＳＡはMonte Carlo法から派生したフィッティングアルゴリズムで、乱数に基づいてフィッティングを進める点と、熱力学を模して仮想エネルギーや仮想温度を取り扱う点、が特徴となっている。ＳＡ自体は公知の技術であり以下のステップ(i)〜(vi)により行われる。
（i）パラメータを乱数に基づいて１つ選択し、さらに乱数に基づいてそのパラメータを変更する(乱数)。
（ii）変更後の評価を行う。評価関数として仮想的なエネルギーＥを考える。通常のＳＡではＥは二乗誤差の総和である。変更前後のＥの変化をΔＥとする(評価)。
（iii）評価が改善される(ΔＥ＜０)なら、その変更を受け入れる(変更)。
（iv）exp(−ΔＥ/Ｔ)の確率で改悪を受け入れる。
（v）Ｔを少しずつ減じる。
（vi）(i)から繰り返す。
ＳＡでは、(iv)にあるように、ボルツマン統計に従って改悪を受け入れることで、局所解から脱出できる可能性を担保している。そのため、局所解にとらわれることなく、大域解に辿りつくことができ、条件の悪いフィッティングに対しても安定に機能する。また、Ｔを徐々に減じることで、大域解にソフトランディングするようになっている。(i)から(v)まで一通り実行することを１モンテカルロステップと呼ぶ。ＳＡではモンテカルロステップを無数に繰り返すことで、フィッティングが進行する。計算に必要な時間は、１モンテカルロステップに必要な時間×必要なモンテカルロステップ数である。必要なモンテカルロステップ数は、パラメータの数・自由度に比例する。
次に、請求項に沿った形で本実施形態の説明を行う。ＣＴをフィッティング問題と考えた場合、フィッティングパラメータは断面像(ｆ(ｘ,ｙ))である。フィッティング対象となるデータは、測定した投影像(ｐ０(ｒ,θ))である。ｒは投影像を撮影した１次元検出器のチャンネル位置、θは投影角度である。座標の定義を図２に示す。ｐ０(ｒ,θ)は、角度θを変えながら投影像を測定することで得られるデータセットであり、横軸にｒ、縦軸にθを取ると、２次元の画像とみなすことができる。このようなデータをsinogramと呼ぶ。典型的なsinogramを図３に示す。
極論すれば、ＣＴはsinogramから断面像を求める技術である。本実施形態では、暫定の断面像ｆ(ｘ,ｙ)と測定した投影像ｐ０(ｒ,θ)との二乗誤差を考えるが、二乗誤差を計算するために、数１０１を使って、ｆ(ｘ,ｙ)から計算で投影像(ｐ(ｒ,θ))を求める。
数１０１は図２において、ｓの方向に和をとっており、ｆ(ｘ,ｙ)のｓ方向の投影像の計算であることがわかる。こうして求めたｐ(ｒ,θ)を用いて、二乗誤差Ｈは以下のように計算できる。
通常のフィッティングでは、仮想エネルギーＥとして、数１０２をそのまま利用するが、本実施形態では、Ｈ以外にスムージング項とエントロピー項を導入している点が特徴となっている。Ｅの定義は次の通りである。
Ｔは仮想温度(温度パラメータ)、Ｓはエントロピー、σは画素値の標準偏差、ｃはスムージング項の強さを表す係数である。ＴＳがエントロピー項、ｃσがスムージング項である。Ｓやσの定義・計算法は後述する。本実施形態では、これらの数式を使い、図４に示すような手順で断面像を計算する。
・ステップ(ａ)：何らかの形で求めた暫定の断面像ｆ(ｘ,ｙ)に対し、評価関数(以下「エネルギー」と呼ぶ)(Ｅ０)を求める(ＳＴ１０およびＳＴ２０)。
・ステップ(ｂ)：乱数で(ｘ０,ｙ０)およびΔμを選択し、断面像ｆ(ｘ,ｙ)の一部を改変する(ＳＴ３０)。
・ステップ(ｃ)：改変後の断面像ｆ(ｘ,ｙ)に対して、評価関数Ｅ１を求める(ＳＴ４０およびＳＴ５０)。
・ステップ(ｄ)：前記エネルギー(Ｅ０)と前記エネルギー(Ｅ１)との差(ΔＥ)を求める(ＳＴ６０)。
・ステップ(ｅ)：前記改変を採用するか否かを前記差(ΔＥ)および温度パラメータ(Ｔ)を用いた受理関数(典型的にはボルツマン統計)により判定する(ＳＴ７０〜ＳＴ１１０)。
・ステップ(ｆ)：前記ステップ(ａ)に戻る(ＳＴ１２０)。
・ステップ(ｇ)：前記ステップ(ａ)〜(ｆ)の繰り返しが所定回数に達するごとに仮想温度(Ｔ)の値を変更する(ＳＴ１３０)。
・ステップ(ｈ)：前記ステップ(ｅ)における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、ＳＴ８０,ＳＴ１００が実行された場合を「成功」、ＳＴ１１０が実行された場合を「失敗」とし、成功の確率が１０％(この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
ステップ(ａ)からステップ(ｈ)の一連の操作が請求項１２に対応し、本実施形態ではこれら一連の操作を図１の画像再構成処理部が行う。
ＳＡの基本手順(i)〜(vi)との対応を考えると、ステップ(ｂ)が（i）、ステップ(ａ)(ｃ)(ｄ)が(ii)、ステップ(ｅ)に(iii)と(iv)が含まれ、ステップ(ｇ)が(v)、ステップ(ｆ)が(vi)にそれぞれ対応する。このことから、本実施形態はＣＴにＳＡを忠実に当てはめたものであることがわかる。
なお、上記ステップ(ｈ)における処理終了の判定ついては、図１の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
（第２の実施形態）
数１０３に示されたＥを計算するためには、数１０１と数１０２を経ることになり、ｓ，ｒ，θに関して級数和を実行することになる。３重の積分(級数和)を計算するので、大変時間がかかる。すなわち、１モンテカルロステップに必要な計算時間が長いことを意味する。さらに、ＣＴではパラメータの数が莫大であり、その結果、ＳＡの実施に必要な総計算時間は、現在の計算機でも「年単位」になってしまう。
そこで、本実施形態では、Ｅを計算するのではなく、変更を行ったときの変化分ΔＥだけを主に計算する。
今、暫定断面像ｆ(ｘ,ｙ)に関する変更をΔｆ(ｘ,ｙ)とする。Δｆ(ｘ,ｙ)は座標(ｘ０,ｙ０)においてのみΔμの値を持ち、それ以外は０であるような断面像である。そのΔｆ(ｘ,ｙ)の投影像Δｐ(ｒ,θ)は数１０１と同じ方法で計算できる。このΔｐ(ｒ,θ)を用いて、ΔＨは以下のように計算できる。
これを書き下して、
となる。今、Δｆ(ｘ,ｙ)は(ｘ０,ｙ０)においてのみ値を持つので、Δｐ(ｒ,θ)はｒ(θ)＝ｘ０cosθ＋ｙ０sinθにおいてのみ値Δμを持ち、それ以外は０である。従って、数１０６の｛｝の中身も、ｒ(θ)＝ｘ０cosθ＋ｙ０sinθにおいてのみ値を持つ。そのため、級数和はｒとθの両方に対して行う必要はなく、以下のように表すことができる。
数１０７ではθに関する級数和だけになっている点が重要である。これにより、計算量を大幅に削減することができる。ｐやｐ０はデジタル画像であるため、数１０７の計算にあたっては、通常は、ｒ(θ)に関する内挿を行う必要がある。そのため、数１０７の｛ ｝内をさらに展開することはできない。しかしながら、誤差を容認して数１０７を展開すると次式が得られる。
ただし、Mは投影角度の数である。数１０７の代わりに数１０８を用いるとΔＨの値に１％程度の誤差が生じる。しかしながら、数１０８は数１０７より高速であり、本発明での利用価値は高い。
次にΔσの計算を説明する。σは座標(ｘ０,ｙ０)付近の輝度値の標準偏差である。(ｘ０,ｙ０)付近として、(ｘ０,ｙ０)の周囲ｄ×ｄ画素を考える(実施例ではｄ=５を用いている)。それらの画素の標準偏差は次式で求めることができる。
ただし、
である。数１０９〜１１１を用いると、Δσは次のように計算することができる。
ただし、ｆｉ，ｆｊは変更前および後のｆ(ｘ０,ｙ０)の値である。
次にΔＳの計算であるが、その前にエントロピーＳの定義を行う。一般に計算機が取り扱う画像はデジタル画像であり、画素の座標(ｘ,ｙ)だけでなく、画素の値もデジタルの値となっている。そこで、１つの画素を１つの量子と考え、画素値を量子状態とみなすことにする。すると画像は複数の量子からなる系であると考えることができる。統計力学に則ると、そのような系に対するエントロピーは次のように定義される。
Ｎは画像中の画素の総数であり、Ｎｉは画素値がデジタル値でｉである画素の総数である。ｋは通常の物理ではボルツマン定数であるが、本発明は物理とは直接関係がないため任意である。本発明では、σと同様にＳも(ｘ０,ｙ０)の周囲ｄ×ｄ画素で定義されるとした。
数１１３で定義されたＳに対してΔＳを考える。変更によって、画素値はデジタル値でｉからｊに変更されるとする。すると、ΔＳは次式のようになる。
すなわち、変更によって、Ｎｉが１減って、Ｎｊが１増えるのである。数１１４を展開すると、非常に簡単な式が得られる。
以上をまとめて、本実施形態の手順を書き下すと次のようになる。
（I）乱数を使って(ｘ０,ｙ０)およびΔμを選択する。
（II）数１０４、１０８(または１０７)、１１２、１１５を使ってΔＥを計算する。
（III）評価が改善される(ΔＥ＜０)なら、ｆ(ｘ０,ｙ０)にΔμを足す。
（IV）ΔＥ＞０の時にも、exp(−ΔＥ/Ｔ)の確率で、ｆ(ｘ０,ｙ０)にΔμを足す。
（V）Ｔを少しずつ減じる。
（VI）（I）から繰り返す。
（VII）所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、（III）,（IV）においてΔμが足された場合を「成功」、それ以外を「失敗」とし、成功の確率が所定値(例えば１０％、この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
手順(I)〜(VI)はＳＡの基本手順(i)〜(vi)を忠実に踏襲していることがわかる。なお、手順(I)〜(VII)の処理は図１の画像再構成処理部において行われる。
なお、上記手順(VII)における処理終了の判定ついては、図１の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
（第３の実施形態）
次に、上で述べた方法の重要な変形について述べる。上で述べた方法では数１０７あるいは数１０８を利用し、θに関する級数和を求めている。そのため、画素１つの変更を検討するのにM回の繰り返し演算を要する。本実施形態はこの部分の計算量をさらに減らすものである。
まず、ｐ(ｒ,θ)の逆投影像ｇ(ｘ,ｙ)を考える。
数１１６では、フィルタ操作をしていないため、ｇ(ｘ,ｙ)は断面像ｆ(ｘ,ｙ)に似ているが、ぼやけたような画像になる。ｇ(ｘ,ｙ)を用いると、数１０８は次のようになる。
ただし、ｇ０(ｘ,ｙ)はｐ０(ｒ,θ)の逆投影像で、数１１６と同様の方法で計算される。数１１７では繰り返し演算がないという点で数１０８より優れている。ｇ０(ｘ,ｙ)は計算過程で全く変化しないため、予め計算しておくことができる。一方、ｇ(ｘ,ｙ)はｆ(ｘ,ｙ)の１点を変更する度にわずかずつ変化するため、厳密には(I)〜(IV)の手順を１回経るごとに再計算が必要である。しかしながら、ｆ(ｘ,ｙ)の変更に伴うｇ(ｘ,ｙ)の変化は少ないという仮定をすると、別の手順を適用することができ、それが本実施形態である。以下に手順を示す。なお、以下の手順(１)〜(９)の処理は図１の画像再構成処理部において行われる。
（１）ｐ０(ｒ,θ)からｇ０(ｘ,ｙ)を求める。
（２）ｆ(ｘ,ｙ)からｐ(ｒ,θ)を求め、さらにｐ(ｒ,θ)からｇ(ｘ,ｙ)を求める。
（３）ｆ(ｘ,ｙ)に対する変更値を画素値とする画像Δμ(ｘ,ｙ)を乱数で生成する。
（４）各画素値について数１１７を適用し、ΔＨ(ｘ,ｙ)という画像を計算する。
（５）同様に各画素値について数１１２、数１１５を適用し、Δσ(ｘ,ｙ)、ΔＳ(ｘ,ｙ)という画像を計算する。
（６）数１０４に基づいてΔＥ(ｘ,ｙ)を計算する。
（７）ΔＥが正である座標(ｘ,ｙ)について、Δμ(ｘ,ｙ)を0にする。
（８）ｆ(ｘ,ｙ)にΔμ(ｘ,ｙ)を足す。
（９）Ｔをα倍(α＜１)し、(２)から繰り返す。
（10）所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了する。ここでは、（７）においてΔμ(ｘ,ｙ)が０にされた場合を「失敗」、されなかった場合を「成功」とし、成功の確率が所定値(例えば１０％、この値は適宜調整可能)を下回ると終了する。そして、終了時点において採用されている推定断面像を測定対象の断面像とし、これをコンピュータのディスプレイに表示したり記録媒体に記録したりする。
なお、上記手順(１０)における処理終了の判定ついては、図１の画像再構成処理部において行うのではなく、コンピュータのディスプレイに推定断面像を逐次表示し、これをユーザが目で見て処理の終了をコンピュータに指示できるようにしてもよい。
手順(I)〜(VII)と、手順(１)〜(１０)を図５に模式的に示した。図５ではパラメータの値(本件では画素の値)を軸とし、仮想エネルギーＥを濃淡で表している(Ｅが小さいほど濃い)。図示の都合により、パラメータの数(画素の数)は２個の場合を取り扱っている。図５(ａ)に示す手順(I)〜(VII)では、ジグザグに折れ曲がりながら、時にでたらめな方向に進みながら、最終的にはＥの最小値を探索してゆく。一方、図５(ｂ)に示した手順(１)〜(１０)では、Ｅの勾配を考慮しながら、Ｅの最小値を探索してゆく。図５(ｂ)の方が効率が良いということが直感的に理解できるだろう。数１１７による計算量の低減及び、図５(ｂ)に示した探索効率の向上により、本実施形態は計算速度的に有利である。手順(１)〜(１０)はＳＡの基本アルゴリズム(i)〜(vi)を忠実に踏襲するものではないが、(I)〜(VII)の手順の自然な発展形である。
（その他の実施形態）
上述の各実施形態で説明した処理を実行する画像再構成処理部は、これらの処理をコンピュータに実行させるためのプログラム、当該プログラムがインストールされたコンピュータ、さらには、これらの処理を実行する専用ＬＳＩなどを使って実現可能である。
（実施例）
実施例としては、シミュレーションを用いたメタルアーティファクト除去効果を取り扱う。まず、比較のために、[非特許文献１]に掲載されている手法での結果を図６(非特許文献のfig.5)に示す。
図６において、左のカラム(ａ)(ｄ)(ｇ)はMetal Artifactがないオリジナルのファントム画像である。中央(ｂ)(ｅ)(ｈ)は、オリジナルのファントムの所定の位置に不可視領域を設定し、通常のFBP法により再構成した画像である。右(ｃ)(ｆ)(ｉ)は、上記と同様に不可視領域を設定し[非特許文献１]の手法により再構成した画像であり、従来技術で最もうまくMetal Artifactが除去できている例の１つである。
そこで、[非特許文献１]から図６(ａ)の画像を取り込み、文献と同じ位置に不可視領域を設定し、シミュレーションで投影像を計算した。結果を図７に示す。図７が本実施例におけるｐ０(ｒ,θ)に対応する。図７について、本発明によって再構成した結果および、比較のための図６(ｃ)の拡大図を図８に示した。効果は一目瞭然である。本実施形態(ｂ)ではMetal Artifactの特徴を見つけ出すことすら困難である。
ただし本実施形態(ｂ)ではエッジ付近が若干ボケている。これは「断面像を滑らかにする因子」(スムージング項ｃσ)がやや強いためである。Metal Artifactを低減するには、その因子をやや強めに設定する必要があることがわかっている。このように本実施形態によれば、仮想エネルギーＥを計算する際の各因子のバランスには改良の余地があるものの、Metal Artifact除去の効果は非常に高い。
なお、参考までに、スムージング項ｃσ、エントロピー項ＴＳを設定せずに本実施形態のアルゴリズムを実行した場合の結果を図９に示す。図９(ａ)はスムージング項ｃσなし(エントロピー項ＴＳのみ)の場合の結果、図９(ｂ)はエントロピー項ＴＳなし(スムージング項ｃσのみ)の結果、図９(ｃ)はエントロピー項ＴＳ、スムージング項ｃσともになしの場合の結果を示す。図９と図８(ｂ)とを比較すると、本実施形態においては、エネルギーＥの因子としてスムージング項ｃσとエントロピー項ＴＳの両方を設定したほうが、いずれか一方のみを設定した場合、ともに設定しない場合よりもMetal Artifact除去の効果が高いことが分かる。スムージング項だけでは、不透明を中心に放射状に広がる縞模様がアーティファクトして残るという特徴がある。一方、エントロピー項のみでは、断面像がざらつきＳ/Ｎが低いという特徴がある。また、最終的な画質は、徐冷過程におけるスムージング項の係数ｃと仮想温度Ｔの比が重要であるようである。
次に、回転角度制限がある場合の実施例を図１０に示す。図１０(ａ)は、FBPのみ(通常のＣＴ)による再構成結果を示す。図１０(ｂ)は、ILST(最小二乗法に基づく画像復元法)を適用した結果を示す。図１０(ｃ)は、本実施形態による再構成結果を示す。これらを比較すると、角度制限時に現れるアーティファクトに対して本実施形態が有効であることが分かる。角度制限によるアーティファクトには、円状の領域が両端の尖ったアーモンド状の領域になり、アーモンド状領域の尖っていない胴回りの輝度が反転するという特徴がある。本実施形態によると、それら両方の特徴が低減される。
本発明は、Metal Artifactに対して顕著な効果があるので、Metal Artifactが深刻な分野、たとえば歯のＣＴ、金属製インプラントを含むＣＴに特に有用である。
また、本発明は情報欠落のある投影像のセットに対して一般的に有効である。情報欠落が顕著な場合の一つに、投影角度制限がある場合を挙げることができる。投影角度制限が問題になる分野は、三次元電子顕微鏡法、マンモグラフィーのＣＴ版、平行移動ＣＴ法(特開2006-71472号公報)などである。情報欠落が問題になる別の場合として、コーンビームＣＴ、ヘリカルスキャンＣＴなど３次元ＣＴが挙げられる。本発明は、３次元ＣＴに現れるアーティファクトを除去する目的にも有効である。
さらには、情報が極端に少ない系に関しても利用可能である。例えば、蛍光Ｘ線ＣＴ、地震波による地球内部のＣＴなどに有用である。
また、本発明の副次的な効果として、断面像における輝度値(Ｘ線では線吸収係数に対応する)を従来法よりも高い精度で決定できることが挙げられる。この効果は骨密度測定等における精度改善に応用可能である。
本発明を用いると、従来の方法よりコントラストの高い再構成像を得ることができる。そのため、アーティファクト等が問題とならない通常のＸ線ＣＴにおいても、本発明の利用価値は高い。また、本発明はデータの過少に対しても安定であることから、投影像測定の時間短縮、ひいてはＸ線被爆量の低減に有効である。これらをあわせて考えると、本発明は、既存のすべてのＸ線ＣＴ技術を置き換える可能性を秘めている。
Ｘ線ＣＴ装置の模式図である。 断面像ｆ(ｘ,ｙ)、投影像ｐ(ｒ,θ)の関係を示す模式図である。 典型的なｐ０(ｒ,θ)の例（横軸に検出器のチャンネル位置、縦軸に投影角度をとった画像）である。 本発明の基本手順を示すフローチャートである。 (ａ)は、手順(I)〜(VI)のフィッティング経過の模式図である。(ｂ)は、手順(１)〜(９)のフィッティング経過の模式図である。 ［非特許文献１］に掲載されている手法での結果を示す図である。 図６(ａ)からシミュレーションによって求めたsinogramである(白い部分は不可視領域)。 (ａ)は、図６(ｃ)をさらに拡大した図である。(ｂ)は、本発明での結果を示す図である。 仮想エネルギーＥの各項の効果を表した模式図である。 角度制限がある場合のアーティファクトと本発明でのアーティファクト低減効果を示す図である。

Claims (14)

  1. 対象物に放射線を照射して得られた投影像（以下「放射線投影像」と呼ぶ）から前記対象物の断面像を求める装置であって、
    前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数（以下「エネルギー」と呼ぶ）（Ｅ０）を求める手段（ａ）と、
    前記現在の推定断面像の一部を改変する手段（ｂ）と、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー（Ｅ１）を求める手段（ｃ）と、
    前記エネルギー（Ｅ０）と前記エネルギー（Ｅ１）との差（ΔＥ）を求める手段（ｄ）と、
    前記改変を受理するか否かを前記差（ΔＥ）と受理確率を制御する温度パラメータ（Ｔ）とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させる手段（ｅ）と、
    前記手段（ａ）〜（ｅ）による一連の処理の繰り返しが所定回数に達するごとに前記温度パラメータ（Ｔ）の値を変更する手段（ｆ）と、
    を備えることを特徴とする画像再構成装置。
  2. 請求項１において、
    前記手段（ａ）は、
    前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域の標準偏差との和を含むエネルギー（Ｅ０）を求め、
    前記手段（ｃ）は、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域の標準偏差との和を含むエネルギー（Ｅ１）を求める、
    ことを特徴とする画像再構成装置。
  3. 請求項１において、
    前記手段（ａ）は、
    前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域のエントロピーとの和を含むエネルギー（Ｅ０）を求め、
    前記手段（ｃ）は、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域のエントロピーとの和を含むエネルギー（Ｅ１）を求める、
    ことを特徴とする画像再構成装置。
  4. 請求項１において、
    前記手段（ａ）は、
    前記現在の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記現在の推定断面像の局所領域の標準偏差と、前記現在の推定断面像の局所領域のエントロピーとの和を含むエネルギー（Ｅ０）を求め、
    前記手段（ｃ）は、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差と、前記改変後の推定断面像の局所領域の標準偏差と、前記改変後の推定断面像の局所領域のエントロピーとの和を含むエネルギー（Ｅ１）を求める、
    ことを特徴とする画像再構成装置。
  5. 請求項１において、
    前記手段（ａ），（ｃ），（ｄ）に代えて、
    ［数１］によりΔＨを算出し、算出したΔＨを構成要素として含むΔＥ（ΔＥ＝ΔＨ＋…）を求める手段（ｈ）を備える、
    ここで、前記対象物の現在の推定断面像をｆ（ｘ，ｙ）、前記手段（ｂ）による改変部分をΔｆ（ｘ，ｙ）とすると、Δｆ（ｘ，ｙ）は座標（ｘ０，ｙ０）においてのみΔμの値を持ち、それ以外は０であるような断面像である。ｐ（ｒ，θ）は、前記対象物の現在の推定断面像から演算により得られた投影像である。ｐ０（ｒ，θ）は、前記対象物の放射線投影像である。ｒは、投影像を撮影した１次元検出器のチャンネル位置である。θは、投影角度である。ｒ（θ）＝ｘ０ｃｏｓθ＋ｙ０ｓｉｎθである。
    ことを特徴とする画像再構成装置。
  6. 請求項１において、
    前記手段（ａ），（ｃ），（ｄ）に代えて、
    ［数２］によりΔＨを算出し、算出したΔＨを構成要素として含むΔＥ（ΔＥ＝ΔＨ＋…）を求める手段（ｈ）を備える、
    ここで、前記対象物の現在の推定断面像をｆ（ｘ，ｙ）、前記手段（ｂ）による改変部分をΔｆ（ｘ，ｙ）とすると、Δｆ（ｘ，ｙ）は座標（ｘ０，ｙ０）においてのみΔμの値を持ち、それ以外は０であるような断面像である。ｐ（ｒ，θ）は、前記対象物の現在の推定断面像から演算により得られた投影像である。ｐ０（ｒ，θ）は、前記対象物の放射線投影像である。ｒは、投影像を撮影した１次元検出器のチャンネル位置である。θは、投影角度である。ｒ（θ）＝ｘ０ｃｏｓθ＋ｙ０ｓｉｎθである。Ｍは、投影角度の数である。
    ことを特徴とする画像再構成装置。
  7. 請求項５または６において、
    前記手段（ｈ）は、
    ［数３］によりΔσを算出し、算出したΔσに係数ｃを掛けたもの（ｃΔσ）と前記ΔＨとの和を構成要素として含むΔＥ（ΔＥ＝ΔＨ＋ｃΔσ＋…）を求める、
    なお、σは座標（ｘ０，ｙ０）の周囲ｄ×ｄ画素の輝度値の標準偏差であり［数４］により表される。ｆｉ，ｆｊはそれぞれ前記手段（ｂ）による改変前後のｆ（ｘ０，ｙ０）の値である。
    ただし、
    である。
    ことを特徴とする画像再構成装置。
  8. 請求項５または６において、
    前記手段（ｈ）は、
    ［数７］によりΔＳを算出し、算出したΔＳに温度パラメータ（Ｔ）を掛けたもの（−ＴΔＳ）と前記ΔＨとの和を構成要素として含むΔＥ（ΔＥ＝ΔＨ−ＴΔＳ＋…）を求める、
    なお、Ｓは座標（ｘ０，ｙ０）の周囲ｄ×ｄ画素の局所領域画像のエントロピーであり［数８］により表される。
    また、
    Ｎ：前記局所領域画像中の画素の総数、
    Ｎｉ：画素値がデジタル値でｉである画素の総数、
    Ｎｊ：画素値がデジタル値でｊである画素の総数、
    ｋ：定数、
    であり、前記手段（ｂ）による改変により画素値がデジタル値でｉからｊに変更されるものとする。
    ことを特徴とする画像再構成装置。
  9. 請求項１において、
    前記手段（ｅ），（ｆ）に代えて、
    前記改変を受理するか否かを前記差（ΔＥ）と受理確率を制御する温度パラメータ（Ｔ）とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させることを予約する手段（ｅ１）と、
    前記手段（ａ）〜（ｄ）および（ｅ１）による一連の処理の繰り返しが所定回数に達するごとに、前記手段（ｅ１）における予約を前記現在の推定断面像に反映させ、前記温度パラメータ（Ｔ）の値を変更する手段（ｆ１）と、
    を備えることを特徴とする画像再構成装置。
  10. 対象物に放射線を照射して得られた投影像から前記対象物の断面像を求める装置であって、
    前記対象物の放射線投影像ｐ０（ｒ，θ）の逆投影像ｇ０（ｘ，ｙ）をフィルタなしの逆投影操作により計算する手段（ｍ１）と、
    前記対象物の現在の推定断面像ｆ（ｘ，ｙ）から投影像ｐ（ｒ，θ）を計算し、さらに当該投影像ｐ（ｒ，θ）の逆投影像ｇ（ｘ，ｙ）をフィルタなしの逆投影操作により計算する手段（ｍ２）と、
    前記対象物の現在の推定断面像ｆ（ｘ，ｙ）に対する変更値を画素値とする画像Δμ（ｘ，ｙ）を生成する手段（ｍ３）と、
    各画素値について［数９］を適用して、画像ΔＨ（ｘ，ｙ）を生成する手段（ｍ４）と、
    ここで、
    Ｍ：投影角度の数、
    である。
    前記ΔＨ（ｘ，ｙ）を用いてΔＥ（ｘ，ｙ）を計算する手段（ｍ５）と、
    ここで、
    ΔＥ（ｘ，ｙ）は、評価関数Ｅ０（ｘ，ｙ）とＥ１（ｘ，ｙ）との差である。
    Ｅ０（ｘ，ｙ）は、前記推定断面像ｆ（ｘ，ｙ）から演算により得られる投影像ｐ（ｒ，θ）と前記放射線投影像ｐ０（ｒ，θ）との差を含む評価関数である。
    Ｅ１（ｘ，ｙ）は、前記推定断面像ｆ（ｘ，ｙ）に前記手段（ｍ３）により得られた画像Δμ（ｘ，ｙ）を足し合わせたもの｛ｆ（ｘ，ｙ）＋Δμ（ｘ，ｙ）｝から演算により得られる投影像｛ｐ（ｒ，θ）＋Δｐ（ｒ，θ）｝と前記放射線投影像ｐ０（ｒ，θ）との差を含む評価関数である。
    前記ΔＥが正である座標（ｘ，ｙ）について、前記Δμ（ｘ，ｙ）を０にする手段（ｍ６）と、
    前記推定断面像ｆ（ｘ，ｙ）に前記手段（ｍ６）により得られた画像Δμ（ｘ，ｙ）を足し合わせたものを新たな推定断面像ｆ（ｘ，ｙ）として前記手段（ｍ２）〜（ｍ６）による処理を繰り返す手段（ｍ７）と、
    を備えることを特徴とする画像再構成装置。
  11. 対象物に放射線を照射して得られた投影像から前記対象物の断面像を求める装置であって、
    前記対象物の放射線投影像ｐ０（ｒ，θ）の逆投影像ｇ０（ｘ，ｙ）を［数１０］により求める手段（ｍ１）と、
    ここで、
    ｒ：投影像を撮影した１次元検出器のチャンネル位置、
    θ：投影角度、
    である。
    前記対象物の現在の推定断面像ｆ（ｘ，ｙ）から演算により投影像ｐ（ｒ，θ）を求め、さらに当該投影像ｐ（ｒ，θ）の逆投影像ｇ（ｘ，ｙ）を［数１１］により求める手段（ｍ２）と、
    前記対象物の現在の推定断面像ｆ（ｘ，ｙ）に対する変更値を画素値とする画像Δμ（ｘ，ｙ）を生成する手段（ｍ３）と、
    各画素値について［数１２］を適用して、画像ΔＨ（ｘ，ｙ）を生成する手段（ｍ４）と、
    ここで、
    Ｍ：投影角度の数、
    である。
    各画素値について［数１３］を適用して、画像Δσ（ｘ，ｙ）を生成する手段（ｍ５）と、
    なお、σは座標（ｘ０，ｙ０）の周囲ｄ×ｄ画素の輝度値の標準偏差であり［数１４］により表される。ｆｉ，ｆｊはそれぞれ変更前後のｆ（ｘ０，ｙ０）の値である。
    ただし、
    である。
    各画素値について［数１７］を適用して、画像ΔＳ（ｘ，ｙ）を生成する手段（ｍ６）と、
    なお、Ｓは座標（ｘ０，ｙ０）の周囲ｄ×ｄ画素の局所領域画像のエントロピーであり［数１８］により表される。
    また、
    Ｎ：前記局所領域画像中の画素の総数、
    Ｎｉ：画素値がデジタル値でｉである画素の総数、
    Ｎｊ：画素値がデジタル値でｊである画素の総数、
    ｋ：定数、
    であり、画素値がデジタル値でｉからｊに変更されるものとする。
    ［数１９］に基づいてΔＥ（ｘ，ｙ）を計算する手段（ｍ７）と、
    ここで、
    ｃ：係数
    Ｔ：仮想温度（温度パラメータ）
    である。
    前記ΔＥが正である座標（ｘ，ｙ）について、前記Δμ（ｘ，ｙ）を０にする手段（ｍ８）と、
    前記推定断面像ｆ（ｘ，ｙ）に前記手段（ｍ８）により得られた画像Δμ（ｘ，ｙ）を足したものを新たな推定断面像ｆ（ｘ，ｙ）とする手段（ｍ９）と、
    前記Ｔをα倍（α＜１）し、前記手段（ｍ２）〜（ｍ９）による処理を繰り返す手段（ｍ１０）と、
    を備えることを特徴とする画像再構成装置。
  12. 対象物に放射線を照射して得られた投影像（以下「放射線投影像」と呼ぶ）から前記対象物の断面像を求める方法であって、
    前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数（以下「エネルギー」と呼ぶ）（Ｅ０）を求めるステップ（ａ）と、
    前記現在の推定断面像の一部を改変するステップ（ｂ）と、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー（Ｅ１）を求めるステップ（ｃ）と、
    前記エネルギー（Ｅ０）と前記エネルギー（Ｅ１）との差（ΔＥ）を求めるステップ（ｄ）と、
    前記改変を受理するか否かを前記差（ΔＥ）と受理確率を制御する温度パラメータ（Ｔ）とを用いた受理関数に基づいて判定するステップ（ｅ）と、
    前記判定結果を前記現在の推定断面像に反映させて前記ステップ（ａ）に戻るステップ（ｆ）と、
    前記ステップ（ａ）〜（ｆ）の繰り返しが所定回数に達するごとに前記温度パラメータ（Ｔ）の値を変更するステップ（ｇ）と、
    前記ステップ（ｅ）における判定結果が所定の停止条件を満たしているか否かを判定し、満たしていれば処理を終了するステップ（ｈ）と、
    を備えることを特徴とする画像再構成方法。
  13. 対象物に放射線を照射して得られた投影像（以下「放射線投影像」と呼ぶ）から前記対象物の断面像を求めるためのコンピュータプログラムであって、
    コンピュータに、
    前記対象物の現在の推定断面像から演算により得られた投影像と前記放射線投影像との差を含む評価関数（以下「エネルギー」と呼ぶ）（Ｅ０）を求めるステップ（ａ）、
    前記現在の推定断面像の一部を改変するステップ（ｂ）、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー（Ｅ１）を求めるステップ（ｃ）、
    前記エネルギー（Ｅ０）と前記エネルギー（Ｅ１）との差（ΔＥ）を求めるステップ（ｄ）、
    前記改変を受理するか否かを前記差（ΔＥ）と受理確率を制御する温度パラメータ（Ｔ）とを用いた受理関数に基づいて判定するステップ（ｅ）、
    前記判定結果を前記現在の推定断面像に反映させて前記ステップ（ａ）に戻るステップ（ｆ）、
    前記ステップ（ａ）〜（ｆ）の繰り返しが所定回数に達するごとに前記温度パラメータ（Ｔ）の値を変更するステップ（ｇ）、
    を実行させるための画像再構成プログラム。
  14. 対象物に放射線を照射して投影像を得る手段（Ａ）と、
    前記投影像から前記対象物の断面像を求める手段（Ｂ）とを備えており、
    前記手段（Ｂ）は、
    前記対象物の現在の推定断面像から演算により得られた投影像と前記対象物に放射線を照射して得られた投影像（以下「放射線投影像」と呼ぶ）との差を含む評価関数（以下「エネルギー」と呼ぶ）（Ｅ０）を求める手段（ｂ１）と、
    前記現在の推定断面像の一部を改変する手段（ｂ２）と、
    前記改変後の推定断面像から演算により得られた投影像と前記放射線投影像との差を含むエネルギー（Ｅ１）を求める手段（ｂ３）と、
    前記エネルギー（Ｅ０）と前記エネルギー（Ｅ１）との差（ΔＥ）を求める手段（ｂ４）と、
    前記改変を受理するか否かを前記差（ΔＥ）と受理確率を制御する温度パラメータ（Ｔ）とを用いた受理関数に基づいて判定し、判定結果を前記現在の推定断面像に反映させる手段（ｂ５）と、
    前記手段（ｂ１）〜（ｂ５）による一連の処理の繰り返しが所定回数に達するごとに前記温度パラメータ（Ｔ）の値を変更する手段（ｂ６）とを備える、
    ことを特徴とするＣＴ装置。
JP2008544217A 2006-11-13 2007-11-13 画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置 Expired - Fee Related JP5190825B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2008544217A JP5190825B2 (ja) 2006-11-13 2007-11-13 画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2006307058 2006-11-13
JP2006307058 2006-11-13
JP2008544217A JP5190825B2 (ja) 2006-11-13 2007-11-13 画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置
PCT/JP2007/072339 WO2008059982A1 (fr) 2006-11-13 2007-11-13 Dispositif de reconfiguration d'image, procédé de reconfiguration d'image, programme de reconfiguration d'image, appareil de tomodensitométrie

Publications (2)

Publication Number Publication Date
JPWO2008059982A1 true JPWO2008059982A1 (ja) 2010-03-04
JP5190825B2 JP5190825B2 (ja) 2013-04-24

Family

ID=39401784

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2008544217A Expired - Fee Related JP5190825B2 (ja) 2006-11-13 2007-11-13 画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置

Country Status (6)

Country Link
US (1) US8090182B2 (ja)
EP (1) EP2092891B1 (ja)
JP (1) JP5190825B2 (ja)
KR (1) KR20090079994A (ja)
CN (1) CN101553171B (ja)
WO (1) WO2008059982A1 (ja)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5543448B2 (ja) * 2008-06-27 2014-07-09 アール．ヤーリッシュ ウォルフラム 高効率コンピュータ断層撮影
US8660330B2 (en) * 2008-06-27 2014-02-25 Wolfram Jarisch High efficiency computed tomography with optimized recursions
DE102010034918A1 (de) * 2010-08-20 2012-02-23 Siemens Aktiengesellschaft Verfahren und Vorrichtung zum Bereitstellen von Güteeinformation für eine Röntgenbildgebung
FR2969793B1 (fr) 2010-12-23 2013-01-11 Gen Electric Reconstruction tomographique d'un objet en mouvement
DE102012206714A1 (de) * 2011-08-10 2013-02-14 Friedrich-Alexander-Universität Erlangen-Nürnberg Verfahren, Recheneinheit, CT-System und C-Bogen-System zur Reduktion von Metallartefakten in CT-Bilddatensätzen
KR101245536B1 (ko) * 2011-10-25 2013-03-21 한국전기연구원 저밀도 촬영상 ｃｔ 영상 재구성에서 줄 인공물 억제 방법
JP6122269B2 (ja) * 2011-12-16 2017-04-26 キヤノン株式会社 画像処理装置、画像処理方法、及びプログラム
US9741127B1 (en) * 2013-05-10 2017-08-22 Shimadzu Corporation Image processing device
US9495770B2 (en) * 2013-08-14 2016-11-15 University Of Utah Research Foundation Practical model based CT construction
WO2016002084A1 (ja) * 2014-07-04 2016-01-07 株式会社島津製作所 画像再構成処理方法
CN105243678B (zh) * 2015-09-23 2018-01-09 倪昕晔 放疗中基于mvcbct和kvct的金属伪影去除方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS5917334A (ja) * 1982-07-21 1984-01-28 株式会社東芝 心拍連動画像診断装置
JPS6017568A (ja) 1983-07-11 1985-01-29 Hitachi Ltd 画像処理方法および装置
US4969110A (en) * 1988-08-01 1990-11-06 General Electric Company Method of using a priori information in computerized tomography
JP3597918B2 (ja) * 1995-09-11 2004-12-08 株式会社日立メディコ Ｘ線ｃｔ装置
JPH09149902A (ja) * 1995-12-01 1997-06-10 Hitachi Medical Corp Ｘ線断層撮影方法および装置
US5907594A (en) 1997-07-01 1999-05-25 Analogic Corporation Reconstruction of volumetric images by successive approximation in cone-beam computed tomography systems
US6639965B1 (en) * 1999-09-30 2003-10-28 General Electric Company Methods and apparatus for cardiac imaging with conventional computed tomography
JP3707347B2 (ja) * 2000-04-07 2005-10-19 株式会社島津製作所 Ｘ線ｃｔ装置の画像処理方法及びｘ線ｃｔ装置並びにｘ線ｃｔ撮影用記録媒体
US6873678B2 (en) * 2000-12-28 2005-03-29 Ge Medical Systems Global Technology Company Llc Methods and apparatus for computed tomographic cardiac or organ imaging
US7477928B2 (en) * 2002-05-17 2009-01-13 Ge Medical Systems Global Technology Company, Llc Method and system for associating an EKG waveform with a CT image
JP2006017568A (ja) * 2004-07-01 2006-01-19 Ricoh Elemex Corp 超音波流量計および受信回路
US20110142316A1 (en) * 2009-10-29 2011-06-16 Ge Wang Tomography-Based and MRI-Based Imaging Systems

Also Published As

Publication number Publication date
WO2008059982A1 (fr) 2008-05-22
EP2092891A1 (en) 2009-08-26
EP2092891B1 (en) 2014-09-17
EP2092891A4 (en) 2011-08-31
KR20090079994A (ko) 2009-07-22
CN101553171A (zh) 2009-10-07
US8090182B2 (en) 2012-01-03
US20090279768A1 (en) 2009-11-12
CN101553171B (zh) 2011-11-16
JP5190825B2 (ja) 2013-04-24

Similar Documents

Publication Publication Date Title
JP5190825B2 (ja) 画像再構成装置、画像再構成方法、画像再構成プログラム、ｃｔ装置
Xu et al. A practical cone-beam CT scatter correction method with optimized Monte Carlo simulations for image-guided radiation therapy
US9934597B2 (en) Metal artifacts reduction in cone beam reconstruction
Jia et al. GPU-based iterative cone-beam CT reconstruction using tight frame regularization
US8000435B2 (en) Method and system for error compensation
US8594407B2 (en) Plane-by-plane iterative reconstruction for digital breast tomosynthesis
US8442353B2 (en) Incorporation of mathematical constraints in methods for dose reduction and image enhancement in tomography
US8724889B2 (en) Method and apparatus for CT image reconstruction
US10282872B2 (en) Noise reduction in tomograms
US20200151921A1 (en) Methods for metal artifact reduction in cone beam reconstruction
Maier et al. Fast simulation of x-ray projections of spline-based surfaces using an append buffer
Heußer et al. Prior‐based artifact correction (PBAC) in computed tomography
CA2729607A1 (en) Incorporation of mathematical constraints in methods for dose reduction and image enhancement in tomography
Xu et al. Statistical iterative reconstruction to improve image quality for digital breast tomosynthesis
Dillon et al. Evaluating reconstruction algorithms for respiratory motion guided acquisition
Jin et al. Bone-induced streak artifact suppression in sparse-view CT image reconstruction
Park et al. A fully GPU-based ray-driven backprojector via a ray-culling scheme with voxel-level parallelization for cone-beam CT reconstruction
US7269244B2 (en) Methods and apparatus for generating thick images in cone beam volumetric CT
Matej et al. Analytic TOF PET reconstruction algorithm within DIRECT data partitioning framework
JP2017221339A (ja) Ｘ線ｃｔ画像再構成方法およびコンピュータプログラム
Qiu et al. New iterative cone beam CT reconstruction software: parameter optimisation and convergence study
KR101493683B1 (ko) 콘-빔 기반 반응선 재구성을 이용한 초고해상도 ｐｅｔ 영상 재구성 장치 및 그 방법
Sun Rigid motion correction for head CT imaging
Ni Reduce Penumbra for 3D X-Ray Reconstruction
Shen et al. Image Quality Enhancement for Digital Breast Tomosynthesis: High-Density Object Artifact Reduction

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20101115

RD02 Notification of acceptance of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7422

Effective date: 20120416

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

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20130121

R150 Certificate of patent or registration of utility model

Ref document number: 5190825

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20160208

Year of fee payment: 3

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

LAPS Cancellation because of no payment of annual fees