JP2005044050A - 画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体 - Google Patents

Abstract

【課題】欠損部分の推定を行わず、ＣＴから得られる原画像データを用いて、鮮明な画像を高速に生成することが可能な画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体を提供することを目的とする。
【解決手段】３次元対象物から得られた３次元ＣＴデータを収集する手順（Ｓ１０）と、前記３次元対象物から得られた３次元ＣＴデータから、補正値を生成するために用いられる閾値を設定する閾値設定手順（Ｓ１１）と、補正対象の３次元ＣＴデータと補正対象の３次元ＣＴデータの周辺の複数の３次元ＣＴデータとからなる３次元ＣＴデータブロックの平均値を求める平均値算出手順（Ｓ１２）と、前記閾値設定手順で設定された閾値と前記平均値算出手順で求められた平均値とを用いて、前記３次元ＣＴデータの値を補正する補正手順（Ｓ１３）と、補正された次元ＣＴデータを用いて、画像を出力する手順（Ｓ１４）とを有する画像処理方法である。
【選択図】 図１３

Description

【０００１】
【発明の属する技術分野】
本発明は、画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体に関する。
【０００２】
【従来の技術】
医科領域においてはＣＴ（Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ）やＭＲＩ（Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ Ｉｍａｇｉｎｇ）から得られた断層画像を用いた診断が一般的となっている。さらに近年のコンピュータハードウェア技術の発展により、これら断層画像から３次元画像を作成することが容易となり、患部の位置や形状のより正確で客観的な把握が可能となった。また、最近では３次元画像を用いた骨片の切断、合成シミュレーションや皮切開部位の距離計測など多くの研究が行われている。
【０００３】
歯科領域においては、これまで医科用ＣＴを流用してきたが装置の規模、画像の解像度などの点から歯科に最適化されているとは言えなかった。
【０００４】
これに対し、日本大学歯学部放射線学教室において新井らが、３ＤＸ（３ＤＸ Ｍｕｌｔｉ Ｉｍａｇｅ Ｍｉｃｒｏ ＣＴ）を開発した（非特許文献１、２参照。）。
【０００５】
この３ＤＸ装置は、撮影領域が局所的であることや装置自体がＸ線照射線量を低く抑えていることから、３ＤＸ装置全体でＸ線被曝線量が低減されている。撮影領域の濃度範囲は２５６階調で表現されるよう量子化されている。この３ＤＸ装置により得られた３次元投影データは従来のＣＴデータと比較して高い解像度を有し、任意方向からの断層像を３次元的に表示することが可能であり、これまで困難であった顎骨や内耳の硬組織の微細な観察が可能である。これにより耳鼻科および歯科領域の硬組織診断での有効性が確認されている。
【０００６】
この３ＤＸ装置を用いて、発明者らは、これまでに、３ＤＸから得られた３次元投影データをＺ軸（対象物の中心を通る鉛直線）を中心に多方向から２次元画像上に展開することにより断層像を作製し、そこでノイズ低減、対象物の輪郭抽出を行った後、３次元上に再配列を行う方法を考案した。これは３ＤＸから得られるデータが３次元上のボクセルデータであることに着目した手法である。
【０００７】
以下、図１を用いて、その処理概要を説明する。
（画像の切り出し）
まず、図１（Ａ）に示されている３ＤＸ装置を用いて、３次元投影データ（図１（Ｂ））に対して、Ｚ軸を中心に多方向から２次元画像に切り出す（図１（Ｃ））。画像の切り出しは断層像を積層する方法に対し逆のアルゴリズムである。ここで、単純に空間分解能である０．１２５（１ボクセル）間隔で切り出すのではなく、ノイズ低減処理を含めるために移動平均処理を行い、切り出す深度方向に８枚の画像の平均値で１枚の画像となるよう画像を切り出す。
（２次元画像処理）
次に、得られた原画像からノイズを除去する。単なる２値化処理では対象物を正確に抽出することが困難であるため、以下で述べる処理を行い、２値化輪郭画像を得る（図１（Ｄ））。
▲１▼濃度変換によるコントラスト調整
▲２▼ガウシアンフィルターを用いた輪郭検出
▲３▼２値化
▲４▼細線化処理による輪郭抽出
（３次元画像の再構築）
得られた輪郭画像を２次元上に画像を切り出す際とは逆に、切り出した順番に各方向毎に積層させ３次元上に再配列を行い、方向毎に構築した３次元画像を３次元上において合成する（図１（Ｅ））。この際、各方向によって抽出される対象物領域が異なるため、欠落している対象物の輪郭情報も他方向からの処理画像によって補間することができる。
（多軸切り出し）
なお、図１（Ｃ）のように、３次元投影データをＺ軸を中心に２次元画像に切り出す以外に、図２に示す正四面体（Ａ）、正六面体（Ｂ）、正八面体（Ｃ）、正一二面体（Ｄ）、正二十面体（Ｅ）の正多面体の各面の中心と正多面体の中心を結ぶ線を軸として２次元画像に切り出すようにしてもよい。
【０００８】
この３ＤＸ装置を用いて、これまでに５０００症例の画像診断を行い、その有効性が確認されている。一方歯科医療も高度化により３次元画像を用いた診断が要求されている。３ＤＸから得られる画像は高解像度であるためにノイズが多く、一般に骨領域の抽出に用いられている２値化処理では対象物表面が欠落し、良好な３次元画像を得ることができなかった。欠落、切断された輪郭線を修復する手法としては欠損部分または離散点の曲率を用いて滑らかに補間する方法や人間の視覚構造、主観的評価をもとに輪郭を検出する方法、ハフ変換を用いた補間方法などが提案されている。また歯科領域においては最適化に基づいたトレース図の生成方法について報告されている。
【０００９】
また、少ないデータ量にて高画質な３次元画像を生成する超音波診断装置において、収集座標のボリュームデータを所定の視線方向に基づいて定義される面に沿って分割し、複数のスライスデータを生成するスライスデータ生成し、複数のスライス面それぞれに対応する中間画像を求め、その中間画像を累積加算することにより表示画像を生成する技術知られている（特許文献１参照。）。
【００１０】
ところで、上記したこれらの手法は滑らかな連続線を検出することが可能であり有効性が確認されているが、いずれも欠損部分を何らかの方法により推定し、人工的な線で補間するものである。
【００１１】
そこで、発明者らはこれまでにノイズを多く含んだＣＴデータから対象物を抽出する新しい手法を考案し、これまでにその有効性について基本的な検討を発表した（非特許文献３、４、５、６参照。）。
【００１２】
また、医用画像処理装置において、画像再構成されたデータに対して後補正処理としてマトリックスフィルタをかける際にその画像に適するマトリックスフィルタを自動的に選択して後補正処理する技術知られている（特許文献２参照。）。
【００１３】
【特許文献１】
特開２００３−６１９５６号公報
【００１４】
【特許文献２】
特開平７−３８５号公報
【００１５】
【非特許文献１】
Ａｒａｉ Ｙ、 Ｈｏｎｄａ Ｋ、 Ｉｗａｉ Ｋ ｅｔ ａｌ ： Ｐｒａｃｔｉｃａｌ Ｍｏｄｅｌ ”３ＤＸ” ｏｆ Ｌｉｍｉｔｅｄ Ｃｏｎｅ−Ｂｅａｍ Ｘ−ｒａｙ ＣＴ ｆｏｒ Ｄｅｎｔａｌ Ｕｓｅ． ＣＡＲＳ２００１ ： ２００１、 ｐｐ．６７１−６７５
【００１６】
【非特許文献２】
小照射野Ｘ線ＣＴの実験機”３ＤＸ Ｍｕｌｔｉ Ｉｍａｇｅ Ｍｉｃｒｏ ＣＴ”の基本性能、 歯科放射．２０００；４０（２）、 ２０００、 ｐｐ．１４５−１５４
【００１７】
【非特許文献３】
別府嗣信、綱島均、新井嘉則：Ｏｒｔｈｏ−ＣＴを用いた３次元モデル構築に関する研究、 電子情報通信学会技術研究報告 ＭＩ２０００−７５、 ２００１、 ｐｐ．１４５−１４９
【００１８】
【非特許文献４】
綱島均、別府嗣信、新井嘉則：立体画像構成法（特願２０００−３５８４２０）、２０００
【００１９】
【非特許文献５】
Ｂｅｆｕ Ｓ、Ｔｓｕｎａｓｈｉｍａ Ｈ、 Ａｒａｉ Ｙ ： Ａ ｓｔｕｄｙ ｉｎ Ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ Ｉｍａｇｅ Ｐｒｏｃｅｓｓｉｎｇ Ｍｅｔｈｏｄ ｆｏｒ ３ＤＸ Ｍｕｌｔｉ Ｉｍａｇｅ Ｍｉｃｒｏ ＣＴ． ＣＡＲＳ２００１： ２００１、 ｐｐ．６６５−６７０
【００２０】
【非特許文献６】
綱島均、別府嗣信、山田鮎太、新井嘉則 ： 歯科用小型Ｘ線ＣＴにおける３次元画像構築法、Ｍｅｄ． Ｉｍａｇ． Ｔｅｃｈ． ２１：１５７−１６５、 ２００３
【００２１】
【発明が解決しようとする課題】
しかしながら、特許文献２に記載された発明は、予め、画像に対応したマトリックスフィルタを用意する必要がある。したがって、特殊な画像の場合は、その画像のために別途、マトリックスフィルタを作成し用意する必要があるという問題がある。
【００２２】
また、非特許文献３、４、５、６に記載された発明において、構築した３次元画像を検討した結果、切り出し回転中心から外周方向に向かって画像が膨張していることが確認され、対象物の寸法の再現に問題があることがわかった。
【００２３】
本発明は、上記問題に鑑みなされたものであり、欠損部分の推定を行わず、ＣＴから得られる原画像データを用いて、鮮明な画像を高速に生成することが可能な画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体を提供することを目的とするものである。
【００２４】
【課題を解決するための手段】
上記課題を解決するために、本件発明は、以下の特徴を有する課題を解決するための手段を採用している。
【００２５】
請求項１に記載された発明は、３次元対象物から得られた３次元ＣＴデータの値を補正する画像処理方法において、前記３次元対象物から得られた３次元ＣＴデータから、補正値を生成するために用いられる閾値を設定する閾値設定手順と、補正対象の３次元ＣＴデータと補正対象の３次元ＣＴデータの周辺の複数の３次元ＣＴデータとからなる３次元ＣＴデータブロックの平均値を求める平均値算出手順と、前記閾値設定手順で設定された閾値と前記平均値算出手順で求められた平均値とを用いて、前記３次元ＣＴデータの値を補正する補正手順とを有することを特徴とする画像処理方法。
【００２６】
請求項２に記載された発明は、請求項１記載の画像処理方法において、補正される前記３次元ＣＴデータの値をＶｏｘｅｌｉｎ、前記閾値をＴｈｒ_ｖｏｌ、前記平均値をＡ_Ｖｎとするとき、補正された前記３次元ＣＴデータの値、Ｖｏｘｅｌｏｕｔは、Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌ（ｘ，ｙ，ｚ）−Ｔｈｒ_ｖｏｌ（Ｔｈｒ_ｖｏｌ＞Ａ_Ｖｎ）、Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌ（ｘ，ｙ，ｚ）−Ａ_Ｖｎ（Ｔｈｒ_ｖｏｌ＜Ａ_Ｖｎ） であることを特徴とする請求項１記載の画像処理方法。
【００２７】
請求項３に記載された発明は、請求項１記載の画像処理方法において、補正される前記３次元ＣＴデータの値をＶｏｘｅｌｉｎ、前記閾値をＴｈｒ_ｖｏｌ、前記平均値をＡ_Ｖｎとするとき、補正された前記３次元ＣＴデータの値、Ｖｏｘｅｌｏｕｔは、Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌｉｎ−（Ａ_Ｖｎ−Ｔｈｒ_ｖｏｌ）であることを特徴とする。
【００２８】
請求項４に記載された発明は、請求項１記載の画像処理方法において、補正される前記３次元ＣＴデータの値をＶｏｘｅｌｉｎ、前記閾値をＴｈｒ_ｖｏｌ、前記平均値をＡ_Ｖｎとするとき、補正された前記３次元ＣＴデータの値、Ｖｏｘｅｌｏｕｔは、Ｔｈｒ_ｖｏｌ＞＝Ａ_Ｖｎの場合は、Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌｉｎ−（Ａ_Ｖｎ−Ｔｈｒ_ｖｏｌ）であり、Ｔｈｒ_ｖｏｌ＜Ａ_Ｖｎの場合は、Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌｉｎであることを特徴とする。
【００２９】
請求項５に記載された発明は、３次元対象物から得られた３次元ＣＴデータを処理する画像処理方法において、注目している３次元ＣＴデータを起点とする複数の方向毎に、注目している３次元ＣＴデータを起点として連続した所定数の３次元ＣＴデータの積算値を算出する積算値算出手順と、前記算出手順で算出された方向毎に求めた積算値の内、上位所定数の積算値の和を求める和算出手順と、所定の閾値と前記和算出手順で積算された積算値とを比較して、前記積算値が前記閾値より大きい場合に、注目している３次元ＣＴデータを処理対象のデータとすることを判定する判定手順と、判定手順とを有することを特徴とする。
【００３０】
請求項６に記載された発明は、請求項５記載の画像処理方法において、前記閾は、前記積算値算出手順で算出された全方向の積算値の平均値と、前記積算値算出手順で算出された全方向の積算値の最大値に基づいて、設定されることを特徴とする。
【００３１】
請求項７に記載された発明は、３次元対象物から得られた３次元ＣＴデータを処理する画像処理方法において、注目している３次元ＣＴデータを起点とする複数の方向毎に、注目している３次元ＣＴデータを起点として連続した所定数の３次元ＣＴデータの積算値を算出する積算値算出手順と、前記算出手順で算出された方向毎に求めた積算値の内、上位所定数の積算値の和及び下位所定数の積算値の和を求める和算出手順と、注目している３次元ＣＴデータを、上位所定数の積算値の和及び下位所定数の積算値の和に基づいて、補正する補正手順と、前記補正手順で補正された３次元ＣＴデータ値と所定の閾値とを比較して、前記補正された３次元ＣＴデータ値が前記閾値より大きい場合に、注目している３次元ＣＴデータを処理対象のデータとすることを判定する判定手順とを有することを特徴とする。
【００３２】
請求項８に記載された発明は、請求１ないし７記載の画像処理方法をコンピュータに実行させる画像処理プログラムである。
請求項９に記載された発明は、請求項８記載の画像処理プログラムを記録したコンピュータ読取可能な記録媒体である。
【発明の実施の形態】
次に、本発明の実施の形態について図面と共に説明する。
【００３３】
図３に、本発明が適用される３ＤＸ装置の例を示す。
【００３４】
撮像装置１は、３ＤＸ装置であり、放射線（Ｘ線）源１１、検知器１２、アナログ−ディジタル変換器１３、汎用コンピュータ１４から構成されている。放射線源１１は、放射線（Ｘ線）を出射しており、放射線源１１から出射された放射線は、撮像対象２１に照射される。放射線源１１からの放射線は、撮像対象２１を透過して検知器（イメージ・インテンシファイア）１２に入射される。検知器１２は、入射された放射線の強度に応じた検知信号を出力する。
【００３５】
なお、汎用コンピュータ１４は、ＨＤＤ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＦＤＤ等の記録媒体に記録された画像処理プログラムをインストールして、汎用コンピュータ１４が、画像処理を行うようにしてもよい。
【００３６】
また、アナログ−ディジタル変換器１３の処理を、汎用コンピュータ１４が、ソフトウエアにより行うようにしてもよい。その場合は、アナログ−ディジタル変換器１３が無くてもよい。
【００３７】
放射線源１１と検知器１２とは撮像対象２１を挟んで互いに対向して配置されており、Ｚ軸（対象物の中心を通る鉛直線）を中心として、１回転して３６０度方向から計５１２枚の投影データが収集される。検知器１２で検知された検知信号（アナログ信号）は、アナログ−ディジタル変換器１３に供給され、ディジタルデータに変換される。アナログ−ディジタル変換器１３でディジタルデータに変換されたデータは、汎用コンピュータ１４に供給されて画像処理が施されて、撮像対象２１の３次元データを得ることができる。
【００３８】
３ＤＸ装置により得られた３次元投影データデータは、図４（Ｂ）に示されているように、高さ４０ｍｍで直径３０ｍｍ×直径３０ｍｍの円柱型の画像形成領域を有する。また、１ボクセルは、図４（Ａ）に示されているように、０．１２５×０．１２５×０．１２５であるので、図４（Ｂ）に示された円柱型の画像は、２４０×π×（３２０／２）^２のボクセルを有する。
（濃度補正）
従来の濃度補正は、画像のヒストグラムに注目し、対象物の濃度分布領域を広げる処理を行っていた（ヒストグラム強調法）。
【００３９】
図５（Ａ）のような画像の場合、ヒストグラム強調法は、画像に対して下記の式（４）を適用して、図５（Ｂ）に示されているように、濃度値ａ−ｂ区間内のヒストグラムのみを強調するものである。
【００４０】
Ｐｉｘｅｌ_ｏｕｔ＝（ｄ−ｃ）×（Ｐｉｘｅｌ_ｉｎ−ａ）／（ｂ−ａ） ・・・・・（６）
なお、Ｐｉｘｅｌｉｎは補正前の入力画素値、Ｐｉｘｅｌｏｕｔは補正後の出力画素値を示す。
【００４１】
この方法では、顎関節の輪郭を多く含む濃度値区間（ａ、ｂ）を見つけ、次に、濃度値ａ−ｂを濃度値ｃ−ｄにマッピングし直し、部分強調を行なう。例えば、濃度値ａ（＝１００）〜ｂ（＝２００）の領域を濃度値ｃ（＝０）〜ｄ（＝２５５）の領域に広げる。
【００４２】
しかし、この手法は全ての画素に対して同じ処理パラメータ（ａ、ｂ、ｃ、ｄ）を用いて補正しているため、全てのノイズを効果的に除去することが困難であった。例えば図８の画像をこの方法により濃度補正をした図９の場合、画像の下半分は全体的に濃度値が低い（画像が暗い）が、上半分は濃度値が高い（画像が明るい）分布を示している。このことは、画像の上半分のノイズ成分の濃度値と、画像上半分のノイズ成分の濃度値が異なることを意味しており、濃度値の異なる双方の領域に分布するノイズ成分を同じ処理パラメータ（ａ、ｂ、ｃ、ｄ）で除去するのは困難である。
【００４３】
この線形濃度変換はコントラストを向上させノイズを低減することが可能であるが、画像全体のヒストグラムのみに注目しており、画像内の対象物の位置情報は含まれていない。また、３ＤＸより得られる画像の位置によって対象物の濃度値が変化するためこのため、単純な線形濃度変換処理では対象物のノイズ成分のトレード・オフが大きくノイズ除去に限界がある。そこで、位置情報を加味しながら３次元投影データ上において角画素毎に応じてパラメータを適応的に変化させる適応型の濃度変換を行った。
【００４４】
本実施の形態の濃度補正方法は、画像の全ての領域に対して同じ処理パラメータを用いるのでなく、画素毎に適応的に処理パラメータを変化させることにより、効果的にノイズ成分を除去する。
▲１▼先ず、３次元投影データ上において大津が提案した判別分析法（大津展之，判別および最小２乗規準に基づく自動しきい値選定法，電子通信学会論文誌 ’８０／４ Ｖｏｌ．Ｊ６３−Ｄ Ｎｏ．４， １９８０， ｐｐ．３４９−３５６）を用いて、領域分割のための最適閾値Ｔｈｒ_ｖｏｌを算出する。
【００４５】
なお、最適閾値Ｔｈｒ_ｖｏｌは、これに限らず、任意に設定することができる。
▲２▼次に、次式を用いて画素毎に近傍画素の平均値Ａ_Ｖｎを算出する（図６参照）。
【００４６】
【数１】

ここではＮの値を２１とし注目画素（Ｖｏｘｅｌ（ｘ，ｙ，ｚ））の２１×２１×２１近傍の平均画素値を算出した。Ｎとして、２７等の他の値を用いてもよい。
【００４７】
次いで、この近傍平均値を画素毎に独立に算出し、以下の補正処理を行う。Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌ（ｘ，ｙ，ｚ）−Ｔｈｒ_ｖｏｌ（Ｔｈｒ_ｖｏｌ＞Ａ_Ｖｎ） （８）
Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌ（ｘ，ｙ，ｚ）−Ａ_Ｖｎ（Ｔｈｒ_ｖｏｌ＜Ａ_Ｖｎ） （９）
なお、Ｖｏｘｅｌｏｕｔは補正後の出力画素値を示す。
（判別処理の変形）
▲１▼３ＤＸから得られた３次元ＣＴデータの平均画素値Ａを算出する。
▲２▼平均画素値Ａから補正値を生成するための閾値Ｔｈｒ_ｖｏを設定する。
【００４８】
Ｔｈｒ_ｖｏｌ＝ｋ１×Ａ ・・・・（１０）
但し、０＜ｋ１＝＜１
なお、ｋ１は、過去のデータに基づいて、最適な値を予め設定しておく。また、その都度、ｋ１を設定するようにしてもよい。
▲２▼次に、上記式（７）を用いて画素毎に近傍画素の平均値Ａ_ＶＮを算出する．
▲３▼判別処理
（その１）
近傍画素の平均値Ａ_ＶＮとＣＴデータ領域全体の平均値Ａの差Ｃを算出する。
Ａ_ＶＮ−Ａ＝Ｃ ・・・（１１）
Ｃを濃度補正量として次式を用いて補正する。
【００４９】
Ｖｏｘｅｌｏｕｔ ＝ Ｖｏｘｅｌｉｎ −Ｃ ・・・（１２）
Ｖｏｘｅｌｉｎは補正前の入力画素値、Ｖｏｘｅｌｏｕｔは補正後の出力画素値を示す。
（その２）
ｉ）Ａ＞Ａ_ＶＮ 近傍の平均値Ａ_ＶＮとＣＴデータ領域全体の平均値Ａの差Ｃを算出する。
Ａ_ＶＮ−Ａ＝Ｃ ・・・（１３）
Ｃを濃度補正量として次式を用いて補正する。
【００５０】
Ｖｏｘｅｌｏｕｔ ＝ Ｖｏｘｅｌｉｎ − Ｃ ・・・（１４）
Ｖｏｘｅｌｉｎは補正前の入力画素値、Ｖｏｘｅｌｏｕｔは補正後の出力画素値を示す。この場合、濃度値が増加の方向へ補正される（画像が明るくなる）。
【００５１】
ｉｉ）Ａ＜Ａ_ＶＮ 濃度値を補正せず、次式の様にそのまま出力する。
【００５２】
Ｖｏｘｅｌｏｕｔ ＝ Ｖｏｘｅｌｉｎ ・・・（１５）
この場合、明るい画像が、補正されることなく出力する。
【００５３】
（その２）の処理は、図７（Ｂ）に示されているように、Ａ＜Ａ_ＶＮの場合は、画像が暗くなる処理を行わないが、（その１）の処理は、図７（Ａ）に示されているように、Ａ＜Ａ_ＶＮの場合は、画像を暗くする処理を行う。
【００５４】
この処理は画像全体の情報と、局所領域の位置に関する情報（今回の場合平均値）とを判別処理することにより、ノイズ成分の濃度値分布の位置におけるバラツキを低減することが可能である。
【００５５】
図１０は、本方法によって、濃度補正されたものである。図１０と図９とを比較すると、下顎頭、下顎窩の近接部分のノイズ成分（左上部の「もやもや」）を除去することが可能となることが分かる。これは近接部分の近傍平均値が判別分析法による最適閾値より値が高いために、注目点の近傍領域の濃度変化に濃度変化に敏感となったためであると考えられる。
【００５６】
一方、画像の背景部分や対象物の領域内においては、近接部分の近傍平均値が判別分析法による最適閾値より低い値となっている。このため３次元投影データ全体の濃度分布から得られた情報を用いた処理となっており、３次元投影データに定常的に存在する一定濃度値領域を効率よく除去することが可能であった。
【００５７】
本手法は３次元上の画像全体の情報と、画像の注目点を中心とする局所領域とを比較しながら適応的に処理する手法である。これによって画像全体に対し濃度値の低い部分は一定値で変換し、濃度値の高い部分に対しては、その注目点の周辺情報に強く依存した変数で変換され画像全体では主観的に良好に変換することが可能となった。
（境界面の検出）
３ＤＸから得られた３次元ＣＴデータ上で、各画素毎にその点が境界点もしくは面の一部かを判定していく。
▲１▼注目３次元ＣＴデータ点をから、図１１及び図１２に示されている２６方向の処理方向を設定する。
【００５８】
なお、図１１及び図１２に示されている２６方向は、注目３次元ＣＴデータ点Ｐを中心にした立法体を考えた場合、
（１）注目３次元ＣＴデータ点Ｐから、８つの「立法体の角」の方向（１，２，７、８、１１、１２、１５、１６）
（２）注目３次元ＣＴデータ点Ｐから、１２の「辺の中点」の方向（９、１０、１３、１４、１７、１８、１９、２０、２１、２２，２３、２４）
（３）注目３次元ＣＴデータ点Ｐから、６つの「面の中心」の方向（３、４、５、６、２５、２６）である。
なお、注目３次元ＣＴデータ点からの方向はこれに限らず実施できる。例えば、正多角形の各面の中心の方向であってもよい。また、図１１の方向において、方向の付加又は削除を行ってもよい。
▲２▼各方向毎に４ｖｏｘｅｌを追跡し画素値を加算（線積分）する。２６方向の加算値ｆ［０］〜ｆ［２５］が算出される。
【００５９】
なお、線積分の数は、４に限らず実施することができる。
▲３▼２６個の加算値ｆ［０］〜ｆ［２５］を大きい順に並び替える。並び替えの後、上８方向の値の平均値を算出する。
【００６０】
なお、上８方向の値の平均値を算出したが、これに限らず実施することができる。例えば、上８方向の値の合計を用いてもよい。
▲４▼８方向の平均値が、閾値よりも高い場合、この画素は対象物の情報とみなし画像出力する。
【００６１】
この閾値Ｔｈｒ_Ｆとして、前記加算値ｆ［０］〜ｆ［２５］の平均値Ｆ_ＡＶと、前記加算値ｆ［０］〜ｆ［２５］の最大値Ｆ_ＭＡＸに基づいて、閾値が設定される。
【００６２】
つまり、例えば、次式で閾値Ｔｈｒ_Ｆを求める。
Ｔｈｒ_Ｆ＝ｋ２×（Ｆ_ＭＡＸ−Ｆ_ＡＶ） ・・・（１６）
本方法により、エッジが鮮明になった。
【００６３】
なお、上記▲３▼における「上８方向の値の平均」法に代えて、「上８方向の値の平均」と「下８方向の値の平均」の差を求めて、その差に基づいて補正して、補正した値が所定の値より大きい場合に、境界面として認定するようにしてもよい。
【００６４】
つまり、上８方向の値の平均値をＦ８ｍａｘとし、下８方向の値の平均値をＦ８ｍｉｎとし、その差に応じて、注目３次元ＣＴデータ点の値Ｖｏｘｅｌを補正し、その補正値の大きさに基づいて、対象物の情報とみなすか否かを判断する。
【００６５】
具体的には、強調係数Ｋ３（０＜ｋ３）を用いて、注目３次元ＣＴデータ点の値Ｖｏｘｅｌを、次の式（１６）に基づいて、強調された補正Ｖｏｘｅｌ ｅを得る。
【００６６】
Ｖｏｘｅｌ ｅ＝Ｖｏｘｅｌ＋ｋ３（Ｆ８ｍａｘ−Ｆ８ｍｉｎ） ・・・・・・（１６）
この強調されたＶｏｘｅｌ ｅが、所定の大きさを超えるか否かで、対象物の情報とみなすか否かを判断する。
【００６７】
なお、注目３次元ＣＴデータ点が、丁度、面の境界上にあるとき、Ｆ８ｍａｘは大きな値を持ち、Ｆ８ｍｉｎは小さな値を持つので、「Ｆ８ｍａｘ−Ｆ８ｍｉｎ」は、大きな値を持つ。したがって、式（１６）により得られた強調されたＶｏｘｅｌ ｅに基づいて、判断することにより、的確な、境界面の判断を行うことができる。
【００６８】
また、本手法は各方向４ｖｏｘｅｌ、注目点を中心に９×９×９の領域から特徴点を抽出する処理である。しかし、本来９×９×９領域から特徴点を抽出するには膨大な計算量であり処理時間に課題がある。しかしながら、本手法では注目点から２６方向のベクトルの線積分を求めソート（並び替え）を行い特徴点を抽出する処理であり非常に効率よく計算することが可能である。
（画像処理方法）
図１３に、３次元対象物から得られた３次元ＣＴデータの値を補正する画像処理方法の処理フローを示す。
【００６９】
先ず、３ＤＸ装置において、所定の軸を中心として、３次元対象物を１回転して３次元ＣＴデータが収集される。（Ｓ１０）
次いで、３次元対象物から得られた３次元ＣＴデータから、補正値を生成するために用いられる閾値を設定する（Ｓ１１）。
【００７０】
大津が提案した最小２乗規準に基づく自動しきい値選定法又は式（１０）に基づいて、閾値を設定する。
【００７１】
次いで、補正対象の３次元ＣＴデータと補正対象の３次元ＣＴデータの周辺の複数の３次元ＣＴデータとからなる３次元ＣＴデータブロックの平均値を求める（Ｓ１２）。
【００７２】
式（７）を用いて、画素毎に近傍画素の平均値Ａ_Ｖｎを算出する。
【００７３】
前記閾値設定手順で設定された閾値と前記平均値算出手順で求められた平均値とを用いて、３次元ＣＴデータの値を補正する（Ｓ１３）。
【００７４】
３次元ＣＴデータの値を補正は、式（１）及び式（２）、式（３）、式（４）及び式（５）に基づいて、行う。
その後、補正された３次元ＣＴデータを用いて、画像を出力する（Ｓ１４）。
【００７５】
図１４に、３次元対象物から得られた３次元ＣＴデータにおいて、注目している３次元ＣＴデータを処理対象のデータとするか否かを決定する処理フローを示す。
【００７６】
先ず、３ＤＸ装置において、所定の軸を中心として、３次元対象物を１回転して３６０度方向から計５１２枚の投影データが収集される。（Ｓ２０）
次いで、注目している３次元ＣＴデータを起点とする複数の方向を設定する。例えば、注目３次元ＣＴデータ点をから、図１１及び図１２に示されている２６方向の処理方向を設定する（Ｓ２１）。
次いで、複数の方向毎に、注目している３次元ＣＴデータを起点として連続した所定数の３次元ＣＴデータの積算値を算出する。例えば、各方向毎に４ｖｏｘｅｌを追跡し画素値を加算（線積分）する（Ｓ２２）。２６方向の加算値ｆ［０］〜ｆ［２５］が算出される。
【００７７】
次いで、２６個の加算値ｆ［０］〜ｆ［２５］を大きい順に並び替える。並び替えの後、上８方向の値の平均値を算出する（Ｓ２３）。
【００７８】
次いで、注目している３次元ＣＴデータを処理対象のデータとするか否かを決定する（Ｓ２４）。例えば、ステップＳ２３で算出された平均値と所定の閾値とを比較して、所定の閾値より大きい場合に、注目している３次元ＣＴデータを処理対象のデータとする。
【００７９】
補正された３次元ＣＴデータを用いて、画像を出力する（Ｓ２５）。
【００８０】
図１５に、「上８方向の値の平均」と「下８方向の値の平均」の差を求めて、その差に基づいて補正して、補正した値が所定の値より大きい場合に、境界面として認定する方法の処理フローを説明する。
【００８１】
先ず、３ＤＸ装置において、所定の軸を中心として、３次元対象物を１回転して３６０度方向から計５１２枚の投影データが収集される。（Ｓ３０）
次いで、注目している３次元ＣＴデータを起点とする複数の方向を設定する。例えば、注目３次元ＣＴデータ点をから、図１１及び図１２に示されている２６方向の処理方向を設定する（Ｓ３１）。
【００８２】
次いで、複数の方向毎に、注目している３次元ＣＴデータを起点として連続した所定数の３次元ＣＴデータの積算値を算出する。例えば、各方向毎に４ｖｏｘｅｌを追跡し画素値を加算（線積分）する（Ｓ３２）。２６方向の加算値ｆ［０］〜ｆ［２５］が算出される。
【００８３】
次いで、２６個の加算値ｆ［０］〜ｆ［２５］を大きい順に並び替える。並び替えの後、上８方向の値の平均値と、下８方向の値の平均を算出し（Ｓ３３）、注目３次元ＣＴデータ点の値Ｖｏｘｅｌを、上記式（１６）に基づいて補正し、強調された補正値Ｖｏｘｅｌ ｅを得る（Ｓ３４）。
【００８４】
次いで、強調された補正値Ｖｏｘｅｌ ｅに基づいて、この注目している３次元ＣＴデータを処理対象のデータとするか否かを決定する（Ｓ３５）。例えば、所定の閾値より大きい場合に、注目している３次元ＣＴデータを処理対象のデータとする。
【００８５】
補正された３次元ＣＴデータを用いて、画像を出力する（Ｓ３５）。
【００８６】
本実施の形態によれば、濃度変換の際，画像全体の最適しきい値，局所領域での画素平均値を用い判別処理を行うことにより効果的にノイズを低減することが可能であった。特に関節腔などの近接部分におけるノイズ除去が顕著である。
【００８７】
また、本実施の形態によれば、ＣＴから得られる原画像データを用いて、２次元データに変換することなく処理を行うので、高速な処理が可能である。
【発明の効果】
上述の如く本発明によれば、欠損部分の推定を行わず、ＣＴから得られる原画像データを用いて、鮮明な画像を高速に生成することが可能な画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体を提供するこができる。
【００８８】
【図面の簡単な説明】
【図１】従来の画像処理方法の処理概要を説明するための図である。
【図２】正多面体である。
【図３】本発明が適用される３ＤＸ装置の例を説明するための図である。
【図４】３次元投影データデータを説明するための図である。
【図５】ヒストグラム強調法を説明するための図である。
【図６】近傍画素を説明するための図である。
【図７】判別処理の比較のための図である。
【図８】原画像である。
【図９】従来の方法によって、濃度補正をした場合の画像である。
【図１０】本実施の形態の方法によって、濃度補正をした場合の画像である。
【図１１】２６の方向を説明するための図（その１）である。
【図１２】２６の方向を説明するための図（その２）である。
【図１３】３次元対象物から得られた３次元ＣＴデータの値を補正する画像処理方法の処理フローである。
【図１４】３次元対象物から得られた３次元ＣＴデータにおいて、注目している３次元ＣＴデータを処理対象のデータとするか否かを決定する処理フローである。
【図１５】「上８方向の値の平均」と「下８方向の値の平均」の差を求めて、その差に基づいて補正して、補正した値が所定の値より大きい場合に、境界面として認定する方法の処理フローである。
【符号の説明】
１ ３ＤＸ装置
１１ 放射線源
１２ 検知器
１３ アナログ−ディジタル変換器
１４ 汎用コンピュータ
２１ 撮像対象

Claims (9)

1. ３次元対象物から得られた３次元ＣＴデータの値を補正する画像処理方法において、
   前記３次元対象物から得られた３次元ＣＴデータから、補正値を生成するために用いられる閾値を設定する閾値設定手順と、
   補正対象の３次元ＣＴデータと補正対象の３次元ＣＴデータの周辺の複数の３次元ＣＴデータとからなる３次元ＣＴデータブロックの平均値を求める平均値算出手順と、
   前記閾値設定手順で設定された閾値と前記平均値算出手順で求められた平均値とを用いて、前記３次元ＣＴデータの値を補正する補正手順と、
   を有することを特徴とする画像処理方法。
2. 補正される前記３次元ＣＴデータの値をＶｏｘｅｌｉｎ、前記閾値をＴｈｒ_ｖｏｌ、前記平均値をＡ_Ｖｎとするとき、補正された前記３次元ＣＴデータの値、Ｖｏｘｅｌｏｕｔは、
   Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌ（ｘ，ｙ，ｚ）−Ｔｈｒ_ｖｏｌ（Ｔｈｒ_ｖｏｌ＞Ａ_Ｖｎ） （１）
   Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌ（ｘ，ｙ，ｚ）−Ａ_Ｖｎ（Ｔｈｒ_ｖｏｌ＜Ａ_Ｖｎ） （２）
   であることを特徴とする請求項１記載の画像処理方法。
3. 補正される前記３次元ＣＴデータの値をＶｏｘｅｌｉｎ、前記閾値をＴｈｒ_ｖｏｌ、前記平均値をＡ_Ｖｎとするとき、補正された前記３次元ＣＴデータの値、Ｖｏｘｅｌｏｕｔは、
   Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌｉｎ−（Ａ_Ｖｎ−Ｔｈｒ_ｖｏｌ） ・・・（３）
   であることを特徴とする請求項１記載の画像処理方法。
4. 補正される前記３次元ＣＴデータの値をＶｏｘｅｌｉｎ、前記閾値をＴｈｒ_ｖｏｌ、前記平均値をＡ_Ｖｎとするとき、補正された前記３次元ＣＴデータの値、Ｖｏｘｅｌｏｕｔは、
   Ｔｈｒ_ｖｏｌ＞＝Ａ_Ｖｎの場合は、
   Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌｉｎ−（Ａ_Ｖｎ−Ｔｈｒ_ｖｏｌ） ・・・（４）
   であり、
   Ｔｈｒ_ｖｏｌ＜Ａ_Ｖｎの場合は、
   Ｖｏｘｅｌｏｕｔ＝Ｖｏｘｅｌｉｎ ・・・（５）
   であることを特徴とする請求項１記載の画像処理方法。
5. ３次元対象物から得られた３次元ＣＴデータを処理する画像処理方法において、
   注目している３次元ＣＴデータを起点とする複数の方向毎に、注目している３次元ＣＴデータを起点として連続した所定数の３次元ＣＴデータの積算値を算出する積算値算出手順と、
   前記算出手順で算出された方向毎に求めた積算値の内、上位所定数の積算値の和を求める和算出手順と、
   所定の閾値と前記和算出手順で積算された積算値とを比較して、前記積算値が前記閾値より大きい場合に、注目している３次元ＣＴデータを処理対象のデータとすることを判定する判定手順と、
   を有することを特徴とする画像処理方法。
6. 閾値設定手順において、前記閾は、前記積算値算出手順で算出された全方向の積算値の平均値と、前記積算値算出手順で算出された全方向の積算値の最大値に基づいて、設定されることを特徴とする請求項４記載の画像処理方法。
7. ３次元対象物から得られた３次元ＣＴデータを処理する画像処理方法において、
   注目している３次元ＣＴデータを起点とする複数の方向毎に、注目している３次元ＣＴデータを起点として連続した所定数の３次元ＣＴデータの積算値を算出する積算値算出手順と、
   前記算出手順で算出された方向毎に求めた積算値の内、上位所定数の積算値の和及び下位所定数の積算値の和を求める和算出手順と、
   注目している３次元ＣＴデータを、上位所定数の積算値の和及び下位所定数の積算値の和に基づいて、補正する補正手順と、
   前記補正手順で補正された３次元ＣＴデータ値と所定の閾値とを比較して、前記補正された３次元ＣＴデータ値が前記閾値より大きい場合に、注目している３次元ＣＴデータを処理対象のデータとすることを判定する判定手順と、
   を有することを特徴とする画像処理方法。
8. 請求１ないし７記載の画像処理方法をコンピュータに実行させる画像処理プログラム。
9. 請求項８記載の画像処理プログラムを記録したコンピュータ読取可能な記録媒体。

Images