ラドンデータに基づくイメージ関数を処理する方法並びに撮像方法
【課題】ラドンデータに基づいてイメージ関数を処理する。
【課題手段】複数の画素で構成され所定の画素サイズを有するデジタル画像をイメージ関数ｆの近似値として提供するための画像処理方法であって、複数の所定の射影方向νに応じて測定された複数の射影関数を含むラドンデータから前記イメージ関数ｆを提供する工程と、所望のデジタル画像を表す画素値を前記イメージ関数から決定する工程とを含み、前記画素値は各画素の画素サイズ内で少なくとも２つのイメージ関数値に基づいて決定される。さらに、前記画像処理方法を用いた撮像方法及び撮像装置について記述する。
【課題手段】複数の画素で構成され所定の画素サイズを有するデジタル画像をイメージ関数ｆの近似値として提供するための画像処理方法であって、複数の所定の射影方向νに応じて測定された複数の射影関数を含むラドンデータから前記イメージ関数ｆを提供する工程と、所望のデジタル画像を表す画素値を前記イメージ関数から決定する工程とを含み、前記画素値は各画素の画素サイズ内で少なくとも２つのイメージ関数値に基づいて決定される。さらに、前記画像処理方法を用いた撮像方法及び撮像装置について記述する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の所定の射影方向に応じて測定した複数の射影関数を含むラドンデータに基づくイメージ関数を処理する方法に関する。さらに本発明は、上記のイメージ関数処理方法を適用して検査領域を撮像するための方法及び装置に関する。
【背景技術】
【０００２】
試料の非破壊検査は、材料科学、健康診断、考古学、建築技術、セキュリティ問題に関わる技術等、多様な技術分野において重要である。例えばコンピュータ断層撮影法（ＣＴ）による試料の画像を得る一方法は、試料平面に異なる射影方向からＸ線を照射し、次に異なる方向で測定した減衰データに依拠して試料平面を再構成する。測定した減衰データ全体は、ラドン空間におけるいわゆるラドンデータで記述することができる。
【０００３】
今日ではラドンデータに関して多様な再構成方法が知られており、特に反復再構成法やフィルター（補正）逆射影法が挙げられる。ラドンデータからイメージ関数を再構成するため、さらに改良された方法が欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記述されている。この方法はディスク上で直交多項式展開を使うものであり、この方法（以下、「ＯＰＥＤアルゴリズム」という。）を用いることにより、検査領域を表すイメージ関数は、複数の所定の射影方向に応じて検査領域を透過して測定した射影関数の値を乗じた多項式の和として、ラドンデータから決定される。
【０００４】
ＯＰＥＤ式再構成法で得られたイメージ関数は連続関数である。画像を視覚化するための典型的な装置は、デジタル化して出力するので、この連続関数は視覚化された画像を表示するため離散化される。一例として、コンピュータのディスプレイにイメージ関数を表示したり、デジタルデータ用プリンタでイメージ関数をプリントアウトしたりするには、ディスプレイやプリンタの（ディスプレイ）解像度に従ってイメージ関数を離散化することが必要となる。一例として、連続したイメージ関数は画素数５１２×５１２又は１０２４×１０２４のパターンで離散化される。濃淡値に変換されたイメージ関数の値が、ぞれぞれ各画素に割り当てられる。
【０００５】
しかしながら、従来の撮像技術には連続したイメージ関数の離散化において共通の欠点があった。通常、視覚化された画像の局部的分解能（解像度）は、イメージ関数の局部的な特長である空間周波数に比べて本質的に小さい。一方で視覚化された画像の画素サイズはイメージ関数の局部的な特長であるサイズに比べて本質的に大きい。従来の離散化は、各画素においてディラック関数（デルタ関数、単位インパルス関数）を用いたイメージ関数の畳み込み（コンボリューション）を含んでいる。このため、連続したイメージ関数の局部的な特長が視覚化された画像にひずみを生じる。その結果、イメージ関数の局部的な特長としてイメージ関数の高周波成分が、離散化のパターンに匹敵する特有の周期を有するのであれば、アーティファクトが生じることがある。
【０００６】
上述した欠点は、従来のＣＴ撮像のみならず、例えば中性子透過撮像、超音波断層法等、ラドンデータの収集に基づくあらゆる既存のＯＰＥＤ式再構成法や画像処理法においても関連することである。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】欧州特許出願公開第０４０３１０４３号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明の目的は、従来の画像処理法の欠点を排除して、ラドンデータに基づいてイメージ関数を処理するための改良方法を提供することにある。本発明の目的は、ラドンデータの収集に基づく従来の撮像方法の欠点を克服するように改良された撮像方法を提供することにある。また本発明の他の目的は、測定したラドンデータを再構成し処理して検査領域を撮像するように改良した装置を提供することにある。
【課題を解決するための手段及び発明の効果】
【０００９】
上述の目的は、独立請求項の特徴を備える方法又は装置で達成される。本発明の好適な実施形態及び適用例は、従属請求項に記述される。
【００１０】
第１発明によれば、イメージ関数におけるデジタルの近似値を提供する画像処理方法は、各画素の画素サイズ内で少なくとも２つのイメージ関数で決定される画素値を演算する工程を含む。 本質的な利点として、従来より用いられるディラック関数でサンプリングする位置においてイメージ関数の局部的な特長が確率的に発生することによる影響を抑制又は完全に排除可能な、改良された画素値を取得できる。画素値はそれぞれの画素サイズ内でイメージ関数に関する少なくとも２つの局部的な特長によって決定される。本発明者らは、ディラック関数を用いた単一のサンプリングに起因して発生する疑似情報が、画素サイズ内で２つのイメージ関数値を考慮するだけで低減できるという知見を得た。少なくとも２つのイメージ関数値を考慮することで、画像サイズ内でさらなる統合が得られ、これによりアーティファクトを排除することができる。
【００１１】
本発明によれば、連続的に定義したイメージ関数がデジタル化される。「連続的イメージ関数」とは、イメージ関数又は準連続的イメージ関数いずれの解析的表現も含む。準連続的イメージ関数はイメージ関数値に関して離散化されたデータセットで表わされ、それぞれのデータセットは得るべきイメージ関数のデジタル表現の画素サイズに比べて、本質的に小さな局部領域（すなわちエリア）を表す。一例として、局部領域はデジタル画像の画素サイズと比べて、少なくとも１／２、好ましくは少なくとも１／１０（例えば１／１００）である。対応する多数のイメージ関数値は、デジタル化された各画素内におけるイメージ関数の経路を表す。
【００１２】
第２発明は、対象物の検査領域を撮像するための撮像方法である。この撮像方法は、撮像装置でラドンデータを測定する工程と、ラドンデータから導いたイメージ関数を再構成する工程と、このイメージ関数に前記第１発明に係る画像処理方法を適用する工程とを含む。本発明の撮像方法における本質的な利点は、イメージ関数の改良された離散化を主たる画像データの収集と組み合わせることによって得られる。
【００１３】
第３発明によれば、対象物の検査領域を撮像するために配置した撮像装置は、射影関数を測定する測定装置と、この測定された射影関数から導いたイメージ関数を再構成する再構成回路と、前記第１発明に係る方法でイメージ関数を処理する画像処理回路とを備える。
【００１４】
本発明の画像処理方法における好ましい実施形態によれば、画素値は画素サイズ内で２つ以上のイメージ関数値、例えば少なくとも４つ、特に９つ以上のイメージ関数値で決定される。したがって、各画素内におけるイメージ関数の局部的な特長の統合を一層向上させることができる。特に好ましい実施形態によれば、各画素の画素サイズ内で所定のフィルタ関数でイメージ関数を畳み込むことから得られるイメージ関数の近似値を求めることで画素値が決定できるよう、画素サイズ内で考えられるイメージ関数値の数はさらに増加される。好適には、フィルタ関数からイメージ関数における全ての局部的な特長にわたる統合がもたらされ、それにより、アーティファクトの削減を促進できる。一例として、このフィルタ関数は各画素の画素サイズ全体を受け持つステップ関数を備えている。したがって、画素サイズ内のイメージ関数値全てが、各画素に割り当てられた画素値に寄与する。
【００１５】
本発明により処理したイメージ関数は、連続したイメージ関数をもたらす従来の任意の再構成方法に基づくラドンデータから再構成できる。イメージ関数が前述のＯＰＥＤアルゴリズムで決定されれば、さらにメリットを享受できる。ＯＰＥＤアルゴリズムから詳細な局部的特長を有する連続したイメージ関数がもたらされるとともに、イメージ関数に関する本発明の離散化によって画質が特に強力に向上し、かつ画像エラーが抑制される。したがって、本発明によりデジタル化したイメージ関数は、ラドンデータの射影関数の値を乗じた多項式の和として提供されることが好ましい。これに応じて、本発明の撮像装置における再構成回路は、測定された射影関数の値を乗じた多項式の和としてイメージ関数を再構成するように適合させることが好ましい。イメージ関数を離散化させるには余分な処理時間を要するので、検査中の領域に関するデジタル画像を表示するために必要な全体時間が増加することもあり得る。好適には、高速フーリエ変換法（ＦＦＴ法）を以下のようにしてイメージ関数の再構成に実装することで、このような処理時間の増加分を補償できる。
【００１６】
本発明の撮像装置で測定した射影関数は離散した射影プロファイルを有しており、この離散した各射影プロファイルは同一の射影方向（ｖ）を有する複数の射影線（ｊ）に対応した射影値γ（ν，ｊ）を有する。欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に概説されるように、イメージ関数ｆは以下の射影の和で表すことができる。
【００１７】
【数１】
【００１８】
（但し、Ｓｋ，ν、φνは次式で表現され、かつＵkはチェビシェフの多項式である。）
【００１９】
【数２】
【００２０】
【数３】
【００２１】
本発明の好ましい実施形態（いわゆる高速ＯＰＥＤアルゴリズム）によれば、イメージ関数ｆは正弦変換法で中間和
【００２２】
【数４】
【００２３】
【数５】
【００２４】
を決定する工程と、
【００２５】
【数６】
【００２６】
として射影の和を決定する工程とで求められる。
【００２７】
正弦変換法を適用すると、イメージ関数の再構成が本質的に高速化される。最良の例では、本発明による再構成は（２ｍ）２の演算ではなく、２ｍｌｏｇ２２ｍの演算だけで済む。典型的には、この実施形態では処理時間（例えばｍ＝５１２）を百分の一に低減することができる。さらに、投射の和を算出するためのハードウェア回路の設計も改善される。
【００２８】
好ましくは、中間和は以下により求められる。
【００２９】
【数７】
【００３０】
特に中間和は、和の後続値である（ｌ＋１）π／（２ｍ＋２），ｌ＝０，１，…，２ｍ−１間で補間処理により得られるθｌ（ｘ，ｙ）で決定できる。
【００３１】
好ましくは、上述した正弦変換法の実行、すなわちイメージ関数の決定は、高速フーリエ変換法を用いて次式を決定する工程を含む。
【００３２】
【数８】
【００３３】
高速ＯＰＥＤアルゴリズムは、ディスク内で任意のコンパクトな部分集合に均一に収束することが証明されており、高速ＯＰＥＤアルゴリズムが正確にしかも従来のＯＰＥＤアルゴリズムに匹敵する画像を再構成できることが数値テストで示されている。補間工程を有する高速ＯＰＥＤアルゴリズムは遙かに高速であるのみならず、従来のＯＰＥＤアルゴリズムで得られる主な利点も享受できる。
【００３４】
高速ＯＰＥＤアルゴリズムの実行は、前記第１発明に係る画像処理方法と組み合わせる例のみに限定されない。事実、一般に高速ＯＰＥＤアルゴリズムは、例えば欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載されるように、得られたラドンベースのデータを再構成するために用いることができる。したがって、高速ＯＰＥＤアルゴリズムは、本発明においてさらに独立した態様をなす。
【００３５】
本発明の撮像方法は、対象物を撮像装置の測定装置内に配置する工程と、対象物を前記複数の所定の射影方向に沿って指向されたエネルギー入力を受けさせる工程と、射影プロファイルを測定する工程を含むことが好ましい。少なくとも一の対象物と測定装置は、ヘリカル状の射影データを得るため対象物にエネルギー入力に受けさせる工程中に所定の方向に並進移動させることができる。
【００３６】
好ましくは、イメージ関数のデジタル画像は所望の視覚化された画像として表示される。なお一般に「表示する」とは、例えば画像ディスプレイ上に又は印刷形態で表示する等、任意の方法で視覚化された表現形態を意味する。
【００３７】
一般に、本発明はラドンデータ又はラドン様データから画像を生成するために用いることができる。このような画像生成が、医療関係の撮像（例えばＣＴ、ＰＥＴ、ＳＰＥＣＴ、ガンマカメラによる撮像等）における数多くの使用例のように多様な用途に利用できることは、本発明の本質的な利点となる。これに限られず、超音波断層撮像法、光断層撮影法、産業用検査又は生物学的研究のための任意の多次元撮像等のように、さらに多くの用途も考えられる。好ましくは、イメージ関数は、Ｘ線ＣＴ装置、超音波断層撮影装置、ＰＥＴ撮像装置、光断層撮影装置、ガンマ線撮像装置、ＳＰＥＣＴ撮像装置、中性子透過検出システム、又は電気インピーダンス断層撮影装置で測定されたラドンデータから決定される。したがって、検査中の対象物は、生体有機物又はその一部、流体組成物、固形物質、工作物及び／又はセキュリティ上の理由で検査すべき対象物を含むことが好ましい。
【００３８】
本発明の他の側面として、多項式の和からなるデータであって、本発明に係る方法を実行するために設置される撮像装置と相互作用可能な、電子的に読み取り可能なデータを有するデジタル保存用媒体又はコンピュータプログラム製品、コンピュータで読み取り可能な媒体に存在し、本発明の方法を実行するためのプログラムコードを有するコンピュータプログラム、本発明の方法を実行するためのプログラム命令を含むコンピュータ読み取り可能な記憶媒体が挙げられる。
【００３９】
本発明における更なる詳細及び利点を、添付図面を参照しながら以下説明する。
【図面の簡単な説明】
【００４０】
【図１】本発明に係るイメージ関数処理方法の実施形態を示すフローチャートである。
【図２】本発明に係る撮像方法の実施形態を示すフローチャートである。
【図３−４】本発明に係る画素値を決定する概略図である。
【図５】本発明に従って使用される画像再構成アルゴリズムに含まれる補間工程を示す概略図である。
【図６】本発明に係る撮像装置の実施形態を示す概略図である。
【発明を実施するための形態】
【００４１】
本発明を、Ｘ線ベースのコンピュータ断層撮影の分野に応用した好ましい適用例について、以下説明する。本発明は、（例えばＶＩＳ又はＩＲ範囲における例えば中性子又は光のような）他種のエネルギー入力ビーム及び／又は（例えば投影Ｘ線撮影のような）他種の投影型撮像法を適用して類似の態様で実行可能である。さらに、好ましい実施形態に関する以下の説明では、主としてイメージ関数の離散化に言及する。本発明の実行に用いるＣＴ装置の詳細については、従来のＣＴ装置から既知であるので詳細説明を省略する。ＯＰＥＤアルゴリズムに関する詳細及び項については欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載してあり、その記載内容を参照して本願に包含するよう援用する。
１．離散化した画素値の決定
【００４２】
図１によれば、本発明に係るイメージ関数の処理方法における基本的な工程は、イメージ関数を用意する工程１００と、離散化した画素値を決定する工程２００と、（必要に応じて）離散化した画素値に基づいてデジタル画像を表示する工程３００とを含む。図２は、本発明に係る撮像方法に関する実施形態の工程を示す。工程１００は、イメージ関数の処理方法を具体的に適用することによって実行できる。一例として、工程１００は特に図２に示す二次工程１１０、１２０、１３０を含む。これらの二次工程は、欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載されるように、ラドンデータを測定しＯＰＥＤベースで再構成することを概略的に示している。あるいは、工程１００は、分析的表現として又は準連続的データセットとしてイメージ関数を用意することを含むこともできる。一方で工程３００は本発明に不可欠な特長というわけではない。デジタル画像を表示する代わりに、工程２００で決定された画素値を例えば撮像装置又は他の画像処理装置におけるデータ記憶装置に格納することもできる。
【００４３】
特に、二次工程１１０、１２０、１３０は、所定の投射方向で複数のＸ線ビームを検査領域に透過させる工程と、これら複数のＸ線ビームで測定した減衰値を含む射影プロファイルを決定する工程と、この射影プロファイルから導いたイメージ関数ｆを再構成する工程とを含む。
【００４４】
工程２００に関する更なる詳細を、図３及び図４の概略図を参照して説明する。これらの図において、図３はイメージ関数ｆを表すように配置した複数の画素１１を備える画素配列１０を示す。図４Ａは例示的なイメージ関数ｆを有する３つの特定画素１１、１２、１３の拡大図である。
【００４５】
典型的には、画素配列１０で適切な画素数を選択する特定の方法は、イメージ関数を提供するのに使用する撮像装置における測定条件によって決定される。一例として、３０ｃｍの検査領域に関するラドンデータを収集するようにしたＣＴ装置を撮像装置が備えており、また検出装置が０．５ｍｍの分解能を有している場合は、検査領域のデジタル画像を示すためには約６００×６００画素１１を用いればよい。画素配列１０は、例えば、約２００μｍの典型的な画素サイズのＬＥＤディスプレイのようなディスプレイ装置の画素、又は画素配列を備え、この画素配列は例えば画素サイズが８０μｍのデジタルプリンタでプリントアウトできる。
【００４６】
従来の方法によれば、例えば各画素の左側境界に位置するディラック関数でイメージ関数ｆを畳み込むことにより、画素値１’、２’、３’（図４Ｂ）が得られる。その結果、画素値１’、２’、３’は画素境界におけるイメージ関数値と同一となる。このような結果は、画素値がイメージ関数の局部的な特長のみで決定されるので、ある確率特性を有する。図４Ｃに示すように各画素のサイズ内の少なくとも２つのイメージ関数値に基づいて画素値を決定することで、従来の方法を置き換えることができる。
【００４７】
図４Ｃによれば、画素値１、２、３は、例えば、各画素の画素サイズ内のステップ関数を備えるフィルタ関数でイメージ関数ｆを畳み込むことによって決定される。その結果、画素値１、２、３は、各画素１１、１２、１３内のイメージ関数値に関して改良された平均近似値を示す。デジタル化した画素値の算出は同様の方法で２つ又はそれ以上の画像寸法において得られる。
２．イメージ関数の再構成の高速化（高速ＯＰＥＤ）
２．１数学的背景
【００４８】
各画素においてイメージ関数をさらに統合するには余分な処理時間が必要となる。この余分な処理時間は、好ましくは本発明に従い、以下に記載する改良ＯＰＥＤアルゴリズム（いわゆる高速ＯＰＥＤ）を用いてイメージ関数を再構成することで補償される。
【００４９】
欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載されるＯＰＥＤアルゴリズムによれば、イメージ関数の近似値は、測定した投射値から以下の二重和Ａ２ｍとして直接算出できる。
【００５０】
【数９】
【００５１】
上式においては、各々以下の通りである（Ｕｋは第２種のチェビシェフ多項式である）。
【００５２】
【数１０】
【００５３】
【数１１】
【００５４】
【数１２】
【００５５】
Ａ２ｍの構造は直接的な方法で一度だけ高速フーリエ変換（ＦＦＴ）を実行するために使われる。以下の定義式を用いて、Ｓｋ，νの値を離散正弦変換のためのＦＦＴで求めることができる。
【００５６】
【数１３】
【００５７】
したがって、項Ａ２ｍは次式のように書き直せる。
【００５８】
【数１４】
【００５９】
したがって、ＯＰＥＤアルゴリズムの主要な工程は上記の二重和を求めることにあり、逆射影の工程として考えることができる。Ｎ＝２ｍ＋１とすれば、マトリクスＳｋ，νを求めるにはＯ（Ｎ２ｌｏｇＮ）の演算（フロップス「flops」)を必要とする。二重和を求めるにはＯ（Ｎ２）の演算を必要とする。ここで、Ｍ×Ｍの格子を求めるにはＯ（Ｎ２（Ｍ２＋ｌｏｇＮ）が必要となる。特にＭ〜Ｎであれば、Ｏ（Ｎ４）が必要となる。主たる演算に必要なのは格子点における二重和を求めることにある。換言すれば、主として必要なことは逆射影の工程にある。
【００６０】
フィルター逆射影法のアルゴリズムとは異なり、式（４）の和にはＦＦＴを介して求められる離散した畳み込みが含まれていない。ただ、式（２）におけるＵｋは、次式で表現できる。
【００６１】
【数１５】
【００６２】
【数１６】
【００６３】
内側の和（inner sum）は特定点において、離散正弦変換のためのＦＦＴで求めることができる。このことは、補間工程を導入してＦＦＴにより高速で数値を求めることが可能であることを示している。さらに詳述すれば、式（５）の内側の和の後で次式を定義できる。
【００６４】
【数１７】
【００６５】
ただし０＜θ＜π………式(6)
図５において、黒丸はθνを表わし、白丸はθ（ｘ，ｙ）を意味する。
【００６６】
離散正弦変換用のＦＦＴは、以下の計算を効率的に行うために利用できる。
【００６７】
【数１８】
すなわち、式（５）における内側の和は、θν（ｘ，ｙ）＝ξιのとき、ＦＦＴで求めることができる。
【００６８】
補間工程は、線形補間を用いるのが好ましい。所与の（ｘ，ｙ）については、θν（ｘ，ｙ）がξ１とξ１＋１の間にあって、α1,yとαl+1v 間の線形補間の値を式（５）の内側の和への近似値として用いるように、整数ｌが選択される。線形補間は、次式によって与えられる。
【００６９】
【数１９】
【００７０】
但しξ１≦θ≦ξ１＋１である。次にＯＰＥＤアルゴリズムの高速実行（図２の工程１３０）は以下のステップを含む。
【００７１】
ステップ１：各ν＝０、・・・２ｍについて、各ｋ＝０、・・・２ｍの演算に以下のＦＦＴが用いられる（ｍ：正の整数）。
【００７２】
【数２０】
【００７３】
ステップ２：各Ｌ＝０、１、・・・２ｍ−１について、次式を演算するためにＦＦＴが用いられる。
【００７４】
【数２１】
【００７５】
ステップ３：検査領域（半径ｃｏｓ（π／（２ｍ＋１））のディスク）内の各再構成点（ｘ，ｙ）について、整数ｌは次式のように決定される。
【００７６】
【数２２】
【００７７】
【数２３】
【００７８】
最後に、イメージ関数ｆは次式に従って算出される。
【００７９】
【数２４】
【００８０】
但し、ｕν＝（２ｍ＋１）θν／π−（ｌ＋１）………式(7)
【００８１】
以上から、ｓｉｎθvがアルゴリズムの最終ステップで分母に現れることが判る。但し、ｃｏｓθv（ｘ，ｙ）＝ｘｃｏｓφν＋ｙｓｉｎφν＝１のときのみθv＝０となる。これは（ｘ，ｙ）＝（ｃｏｓφν，ｓｉｎφν）のときのみ生じる。点（ｃｏｓφν，ｓｉｎφν）は領域Ｂの境界上にあるので、これらの点を再構成点とみなす必要はない。実際、検査領域は普通、単位ディスク内にあり。したがって、演算はＢ内の小さなディスク内の点に限定することができる。これは確実に、演算において有効桁の消失を起こすのに最終段階のｓｉｎθvの値が小さすぎることがないようにする。さらにまた、（ｘ，ｙ）が半径ｃｏｓξ０＝ｃｏｓπ／（２ｍ＋１）であるディスクに限定される場合、ステップ３におけるｌの選択が該ディスク内のすべての（ｘ，ｙ）に対して固有であることも保証する。
【００８２】
アルゴリズムはＦＦＴを２度用いて、ステップ３における最終和は単一の和となり、それを求めるにはＯ（Ｎ）の演算を要する。したがって、Ｍ〜Ｎを用いてＭ×Ｍの格子上で数値を求める上で必要なものはＯ（Ｎ３）であり、これはＦＢＰ（filtered back projection）アルゴリズムと同じ桁数である。
２．２ 演算と結果
【００８３】
数値演算には、パッケージソフト「ＦＦＴＷ」（http://www.fftw.org/）における離散正弦変換用のＦＦＴを用いた。数値例はいわゆるシェップ−ローガン（Ｓｈｅｐｐ−Ｌｏｇａｎ）の頭部ファントム上で行われる。この頭部ファントムは非常に優れた解析用ファントムであり、画像内のあらゆる楕円の境界近傍で画像がジャンプ（jump）を含み、境界上のジャンプも含む。画像を表す関数はステップ関数であり、これは楕円の境界において連続していない。
【００８４】
画像の再構成は、補間工程のないＯＰＥＤアルゴリズムと高速ＯＰＥＤアルゴリズムとをそれぞれ使って行い、後者は上述したように補間工程を含んでいる。いずれの場合も、Ｓｋ、νはＦＦＴで演算した。
【００８５】
再構成の実行は、インテル社製ＣＰＵ「Ｘｅｎｏｎ（商標）」を２個（各々３０６５ＭＨｚ）と４ＧＢのＲＡＭ（コードはＣ言語で書き込まれる）を備えたＣＥＬＳＩＵＳコンピュタＲ６１０型で行った。ＯＰＥＤアルゴリズムを使った場合、再構成に３４４秒要し、このうち９５％以上は逆射影工程に費やされた。これに対して高速ＯＰＥＤアルゴリズムを用いた場合は、所要時間は僅か１３秒となり、２６倍以上の速度向上が得られた。また、これら２つの画像の見た目は殆ど変わらなかった。
【００８６】
再構成の誤差を求めたところ、相対的な最小二乗平均誤差では高速ＯＰＥＤが僅かに劣り、平均誤差では僅かに勝っていた。誤差の桁数は同じであり、誤差の違いは殆ど無視できるレベルであった。
３．撮像装置
【００８７】
図６に、一実施の形態に係る撮像装置１００を概略的に示す。撮像装置１００は、エネルギー生成器２００及び検出装置３００を備えた測定装置と、この測定装置に接続した再構成回路４００とを含む。さらに保持装置５００を備えており、これは例えばＣＴシステムで知られているようなキャリアテーブルや、検査中の対象物を測定装置に配置したり、対象物の幾何学的形状をエネルギー生成器２００及び検出装置３００に対して調整したりする他の任意のキャリア又は物体保持器である。制御装置、ディスプレイ装置（図示せず）等のような他の構成部品も、これら自体は従来技術から既知の態様で設けられる。
【００８８】
エネルギー生成器２００は、従来のＣＴ装置から判るように、線源キャリア２２０（例えばガイドレール）上に配置される、例えば移動可能なＸ線管のようなエネルギー入力源２１０を含む。検出装置は、エネルギー入力源２１０とは対向する関係で線源キャリア２２０上に可動配置されるセンサ配列３１０を含む。このような構造で、（図面の平面と平行な）ＲＯＩ（Region Of Investigation）を通過する射影方向は、保持装置５００の回りで構成部品２１０、３１０の組合せを回転させることにより設定できる。線源キャリア２２０は対象物の回りでエネルギー生成器２００及び検出装置３００を回転可能にする円環として示される。変形例によれば（変更することにより）、線源キャリアは楕円形状や他の形状であってもよい。これは検査対象物の幾何学的形状に適合させるという点で特長を示すことができる。
【００８９】
測定装置２００、３００に接続した再構成回路４００は、測定したイメージ関数の値を乗じた多項式の和としてイメージ関数ｆを再構成するように適合させてある。再構成回路４００は、上記第１発明に係る画像処理方法でイメージ関数ｆを離散化させるよう構成した画像処理回路４１０を含む。画像処理回路４１０には、例えば本発明を実施するためのプログラム命令を含んだコンピュータで読み取り可能な記憶媒体を備えるマイクロプロセッサが実装される。
【００９０】
上記の詳細な説明、図面、及び特許請求の範囲で開示した本発明の特長は、単独もしくは組み合わせて本発明を様々な実施形態で実現することができる重要なものである。
【符号の説明】
【００９１】
１・・・画素値
２・・・画素値
３・・・画素値
１’・・・画素値
２’・・・画素値
３’・・・画素値
１０・・・画素配列
１１・・・画素
１２・・・画素
１３・・・画素
１００・・・撮像装置
２００・・・エネルギー生成器
２１０・・・エネルギー入力源
２２０・・・線源キャリア
３００・・・検出装置
３１０・・・センサ配列
４００・・・再構成回路
４１０・・・画像処理回路
５００・・・保持装置
ｆ・・・イメージ関数
ν・・・射影方向
【技術分野】
【０００１】
本発明は、複数の所定の射影方向に応じて測定した複数の射影関数を含むラドンデータに基づくイメージ関数を処理する方法に関する。さらに本発明は、上記のイメージ関数処理方法を適用して検査領域を撮像するための方法及び装置に関する。
【背景技術】
【０００２】
試料の非破壊検査は、材料科学、健康診断、考古学、建築技術、セキュリティ問題に関わる技術等、多様な技術分野において重要である。例えばコンピュータ断層撮影法（ＣＴ）による試料の画像を得る一方法は、試料平面に異なる射影方向からＸ線を照射し、次に異なる方向で測定した減衰データに依拠して試料平面を再構成する。測定した減衰データ全体は、ラドン空間におけるいわゆるラドンデータで記述することができる。
【０００３】
今日ではラドンデータに関して多様な再構成方法が知られており、特に反復再構成法やフィルター（補正）逆射影法が挙げられる。ラドンデータからイメージ関数を再構成するため、さらに改良された方法が欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記述されている。この方法はディスク上で直交多項式展開を使うものであり、この方法（以下、「ＯＰＥＤアルゴリズム」という。）を用いることにより、検査領域を表すイメージ関数は、複数の所定の射影方向に応じて検査領域を透過して測定した射影関数の値を乗じた多項式の和として、ラドンデータから決定される。
【０００４】
ＯＰＥＤ式再構成法で得られたイメージ関数は連続関数である。画像を視覚化するための典型的な装置は、デジタル化して出力するので、この連続関数は視覚化された画像を表示するため離散化される。一例として、コンピュータのディスプレイにイメージ関数を表示したり、デジタルデータ用プリンタでイメージ関数をプリントアウトしたりするには、ディスプレイやプリンタの（ディスプレイ）解像度に従ってイメージ関数を離散化することが必要となる。一例として、連続したイメージ関数は画素数５１２×５１２又は１０２４×１０２４のパターンで離散化される。濃淡値に変換されたイメージ関数の値が、ぞれぞれ各画素に割り当てられる。
【０００５】
しかしながら、従来の撮像技術には連続したイメージ関数の離散化において共通の欠点があった。通常、視覚化された画像の局部的分解能（解像度）は、イメージ関数の局部的な特長である空間周波数に比べて本質的に小さい。一方で視覚化された画像の画素サイズはイメージ関数の局部的な特長であるサイズに比べて本質的に大きい。従来の離散化は、各画素においてディラック関数（デルタ関数、単位インパルス関数）を用いたイメージ関数の畳み込み（コンボリューション）を含んでいる。このため、連続したイメージ関数の局部的な特長が視覚化された画像にひずみを生じる。その結果、イメージ関数の局部的な特長としてイメージ関数の高周波成分が、離散化のパターンに匹敵する特有の周期を有するのであれば、アーティファクトが生じることがある。
【０００６】
上述した欠点は、従来のＣＴ撮像のみならず、例えば中性子透過撮像、超音波断層法等、ラドンデータの収集に基づくあらゆる既存のＯＰＥＤ式再構成法や画像処理法においても関連することである。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】欧州特許出願公開第０４０３１０４３号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明の目的は、従来の画像処理法の欠点を排除して、ラドンデータに基づいてイメージ関数を処理するための改良方法を提供することにある。本発明の目的は、ラドンデータの収集に基づく従来の撮像方法の欠点を克服するように改良された撮像方法を提供することにある。また本発明の他の目的は、測定したラドンデータを再構成し処理して検査領域を撮像するように改良した装置を提供することにある。
【課題を解決するための手段及び発明の効果】
【０００９】
上述の目的は、独立請求項の特徴を備える方法又は装置で達成される。本発明の好適な実施形態及び適用例は、従属請求項に記述される。
【００１０】
第１発明によれば、イメージ関数におけるデジタルの近似値を提供する画像処理方法は、各画素の画素サイズ内で少なくとも２つのイメージ関数で決定される画素値を演算する工程を含む。 本質的な利点として、従来より用いられるディラック関数でサンプリングする位置においてイメージ関数の局部的な特長が確率的に発生することによる影響を抑制又は完全に排除可能な、改良された画素値を取得できる。画素値はそれぞれの画素サイズ内でイメージ関数に関する少なくとも２つの局部的な特長によって決定される。本発明者らは、ディラック関数を用いた単一のサンプリングに起因して発生する疑似情報が、画素サイズ内で２つのイメージ関数値を考慮するだけで低減できるという知見を得た。少なくとも２つのイメージ関数値を考慮することで、画像サイズ内でさらなる統合が得られ、これによりアーティファクトを排除することができる。
【００１１】
本発明によれば、連続的に定義したイメージ関数がデジタル化される。「連続的イメージ関数」とは、イメージ関数又は準連続的イメージ関数いずれの解析的表現も含む。準連続的イメージ関数はイメージ関数値に関して離散化されたデータセットで表わされ、それぞれのデータセットは得るべきイメージ関数のデジタル表現の画素サイズに比べて、本質的に小さな局部領域（すなわちエリア）を表す。一例として、局部領域はデジタル画像の画素サイズと比べて、少なくとも１／２、好ましくは少なくとも１／１０（例えば１／１００）である。対応する多数のイメージ関数値は、デジタル化された各画素内におけるイメージ関数の経路を表す。
【００１２】
第２発明は、対象物の検査領域を撮像するための撮像方法である。この撮像方法は、撮像装置でラドンデータを測定する工程と、ラドンデータから導いたイメージ関数を再構成する工程と、このイメージ関数に前記第１発明に係る画像処理方法を適用する工程とを含む。本発明の撮像方法における本質的な利点は、イメージ関数の改良された離散化を主たる画像データの収集と組み合わせることによって得られる。
【００１３】
第３発明によれば、対象物の検査領域を撮像するために配置した撮像装置は、射影関数を測定する測定装置と、この測定された射影関数から導いたイメージ関数を再構成する再構成回路と、前記第１発明に係る方法でイメージ関数を処理する画像処理回路とを備える。
【００１４】
本発明の画像処理方法における好ましい実施形態によれば、画素値は画素サイズ内で２つ以上のイメージ関数値、例えば少なくとも４つ、特に９つ以上のイメージ関数値で決定される。したがって、各画素内におけるイメージ関数の局部的な特長の統合を一層向上させることができる。特に好ましい実施形態によれば、各画素の画素サイズ内で所定のフィルタ関数でイメージ関数を畳み込むことから得られるイメージ関数の近似値を求めることで画素値が決定できるよう、画素サイズ内で考えられるイメージ関数値の数はさらに増加される。好適には、フィルタ関数からイメージ関数における全ての局部的な特長にわたる統合がもたらされ、それにより、アーティファクトの削減を促進できる。一例として、このフィルタ関数は各画素の画素サイズ全体を受け持つステップ関数を備えている。したがって、画素サイズ内のイメージ関数値全てが、各画素に割り当てられた画素値に寄与する。
【００１５】
本発明により処理したイメージ関数は、連続したイメージ関数をもたらす従来の任意の再構成方法に基づくラドンデータから再構成できる。イメージ関数が前述のＯＰＥＤアルゴリズムで決定されれば、さらにメリットを享受できる。ＯＰＥＤアルゴリズムから詳細な局部的特長を有する連続したイメージ関数がもたらされるとともに、イメージ関数に関する本発明の離散化によって画質が特に強力に向上し、かつ画像エラーが抑制される。したがって、本発明によりデジタル化したイメージ関数は、ラドンデータの射影関数の値を乗じた多項式の和として提供されることが好ましい。これに応じて、本発明の撮像装置における再構成回路は、測定された射影関数の値を乗じた多項式の和としてイメージ関数を再構成するように適合させることが好ましい。イメージ関数を離散化させるには余分な処理時間を要するので、検査中の領域に関するデジタル画像を表示するために必要な全体時間が増加することもあり得る。好適には、高速フーリエ変換法（ＦＦＴ法）を以下のようにしてイメージ関数の再構成に実装することで、このような処理時間の増加分を補償できる。
【００１６】
本発明の撮像装置で測定した射影関数は離散した射影プロファイルを有しており、この離散した各射影プロファイルは同一の射影方向（ｖ）を有する複数の射影線（ｊ）に対応した射影値γ（ν，ｊ）を有する。欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に概説されるように、イメージ関数ｆは以下の射影の和で表すことができる。
【００１７】
【数１】
【００１８】
（但し、Ｓｋ，ν、φνは次式で表現され、かつＵkはチェビシェフの多項式である。）
【００１９】
【数２】
【００２０】
【数３】
【００２１】
本発明の好ましい実施形態（いわゆる高速ＯＰＥＤアルゴリズム）によれば、イメージ関数ｆは正弦変換法で中間和
【００２２】
【数４】
【００２３】
【数５】
【００２４】
を決定する工程と、
【００２５】
【数６】
【００２６】
として射影の和を決定する工程とで求められる。
【００２７】
正弦変換法を適用すると、イメージ関数の再構成が本質的に高速化される。最良の例では、本発明による再構成は（２ｍ）２の演算ではなく、２ｍｌｏｇ２２ｍの演算だけで済む。典型的には、この実施形態では処理時間（例えばｍ＝５１２）を百分の一に低減することができる。さらに、投射の和を算出するためのハードウェア回路の設計も改善される。
【００２８】
好ましくは、中間和は以下により求められる。
【００２９】
【数７】
【００３０】
特に中間和は、和の後続値である（ｌ＋１）π／（２ｍ＋２），ｌ＝０，１，…，２ｍ−１間で補間処理により得られるθｌ（ｘ，ｙ）で決定できる。
【００３１】
好ましくは、上述した正弦変換法の実行、すなわちイメージ関数の決定は、高速フーリエ変換法を用いて次式を決定する工程を含む。
【００３２】
【数８】
【００３３】
高速ＯＰＥＤアルゴリズムは、ディスク内で任意のコンパクトな部分集合に均一に収束することが証明されており、高速ＯＰＥＤアルゴリズムが正確にしかも従来のＯＰＥＤアルゴリズムに匹敵する画像を再構成できることが数値テストで示されている。補間工程を有する高速ＯＰＥＤアルゴリズムは遙かに高速であるのみならず、従来のＯＰＥＤアルゴリズムで得られる主な利点も享受できる。
【００３４】
高速ＯＰＥＤアルゴリズムの実行は、前記第１発明に係る画像処理方法と組み合わせる例のみに限定されない。事実、一般に高速ＯＰＥＤアルゴリズムは、例えば欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載されるように、得られたラドンベースのデータを再構成するために用いることができる。したがって、高速ＯＰＥＤアルゴリズムは、本発明においてさらに独立した態様をなす。
【００３５】
本発明の撮像方法は、対象物を撮像装置の測定装置内に配置する工程と、対象物を前記複数の所定の射影方向に沿って指向されたエネルギー入力を受けさせる工程と、射影プロファイルを測定する工程を含むことが好ましい。少なくとも一の対象物と測定装置は、ヘリカル状の射影データを得るため対象物にエネルギー入力に受けさせる工程中に所定の方向に並進移動させることができる。
【００３６】
好ましくは、イメージ関数のデジタル画像は所望の視覚化された画像として表示される。なお一般に「表示する」とは、例えば画像ディスプレイ上に又は印刷形態で表示する等、任意の方法で視覚化された表現形態を意味する。
【００３７】
一般に、本発明はラドンデータ又はラドン様データから画像を生成するために用いることができる。このような画像生成が、医療関係の撮像（例えばＣＴ、ＰＥＴ、ＳＰＥＣＴ、ガンマカメラによる撮像等）における数多くの使用例のように多様な用途に利用できることは、本発明の本質的な利点となる。これに限られず、超音波断層撮像法、光断層撮影法、産業用検査又は生物学的研究のための任意の多次元撮像等のように、さらに多くの用途も考えられる。好ましくは、イメージ関数は、Ｘ線ＣＴ装置、超音波断層撮影装置、ＰＥＴ撮像装置、光断層撮影装置、ガンマ線撮像装置、ＳＰＥＣＴ撮像装置、中性子透過検出システム、又は電気インピーダンス断層撮影装置で測定されたラドンデータから決定される。したがって、検査中の対象物は、生体有機物又はその一部、流体組成物、固形物質、工作物及び／又はセキュリティ上の理由で検査すべき対象物を含むことが好ましい。
【００３８】
本発明の他の側面として、多項式の和からなるデータであって、本発明に係る方法を実行するために設置される撮像装置と相互作用可能な、電子的に読み取り可能なデータを有するデジタル保存用媒体又はコンピュータプログラム製品、コンピュータで読み取り可能な媒体に存在し、本発明の方法を実行するためのプログラムコードを有するコンピュータプログラム、本発明の方法を実行するためのプログラム命令を含むコンピュータ読み取り可能な記憶媒体が挙げられる。
【００３９】
本発明における更なる詳細及び利点を、添付図面を参照しながら以下説明する。
【図面の簡単な説明】
【００４０】
【図１】本発明に係るイメージ関数処理方法の実施形態を示すフローチャートである。
【図２】本発明に係る撮像方法の実施形態を示すフローチャートである。
【図３−４】本発明に係る画素値を決定する概略図である。
【図５】本発明に従って使用される画像再構成アルゴリズムに含まれる補間工程を示す概略図である。
【図６】本発明に係る撮像装置の実施形態を示す概略図である。
【発明を実施するための形態】
【００４１】
本発明を、Ｘ線ベースのコンピュータ断層撮影の分野に応用した好ましい適用例について、以下説明する。本発明は、（例えばＶＩＳ又はＩＲ範囲における例えば中性子又は光のような）他種のエネルギー入力ビーム及び／又は（例えば投影Ｘ線撮影のような）他種の投影型撮像法を適用して類似の態様で実行可能である。さらに、好ましい実施形態に関する以下の説明では、主としてイメージ関数の離散化に言及する。本発明の実行に用いるＣＴ装置の詳細については、従来のＣＴ装置から既知であるので詳細説明を省略する。ＯＰＥＤアルゴリズムに関する詳細及び項については欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載してあり、その記載内容を参照して本願に包含するよう援用する。
１．離散化した画素値の決定
【００４２】
図１によれば、本発明に係るイメージ関数の処理方法における基本的な工程は、イメージ関数を用意する工程１００と、離散化した画素値を決定する工程２００と、（必要に応じて）離散化した画素値に基づいてデジタル画像を表示する工程３００とを含む。図２は、本発明に係る撮像方法に関する実施形態の工程を示す。工程１００は、イメージ関数の処理方法を具体的に適用することによって実行できる。一例として、工程１００は特に図２に示す二次工程１１０、１２０、１３０を含む。これらの二次工程は、欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載されるように、ラドンデータを測定しＯＰＥＤベースで再構成することを概略的に示している。あるいは、工程１００は、分析的表現として又は準連続的データセットとしてイメージ関数を用意することを含むこともできる。一方で工程３００は本発明に不可欠な特長というわけではない。デジタル画像を表示する代わりに、工程２００で決定された画素値を例えば撮像装置又は他の画像処理装置におけるデータ記憶装置に格納することもできる。
【００４３】
特に、二次工程１１０、１２０、１３０は、所定の投射方向で複数のＸ線ビームを検査領域に透過させる工程と、これら複数のＸ線ビームで測定した減衰値を含む射影プロファイルを決定する工程と、この射影プロファイルから導いたイメージ関数ｆを再構成する工程とを含む。
【００４４】
工程２００に関する更なる詳細を、図３及び図４の概略図を参照して説明する。これらの図において、図３はイメージ関数ｆを表すように配置した複数の画素１１を備える画素配列１０を示す。図４Ａは例示的なイメージ関数ｆを有する３つの特定画素１１、１２、１３の拡大図である。
【００４５】
典型的には、画素配列１０で適切な画素数を選択する特定の方法は、イメージ関数を提供するのに使用する撮像装置における測定条件によって決定される。一例として、３０ｃｍの検査領域に関するラドンデータを収集するようにしたＣＴ装置を撮像装置が備えており、また検出装置が０．５ｍｍの分解能を有している場合は、検査領域のデジタル画像を示すためには約６００×６００画素１１を用いればよい。画素配列１０は、例えば、約２００μｍの典型的な画素サイズのＬＥＤディスプレイのようなディスプレイ装置の画素、又は画素配列を備え、この画素配列は例えば画素サイズが８０μｍのデジタルプリンタでプリントアウトできる。
【００４６】
従来の方法によれば、例えば各画素の左側境界に位置するディラック関数でイメージ関数ｆを畳み込むことにより、画素値１’、２’、３’（図４Ｂ）が得られる。その結果、画素値１’、２’、３’は画素境界におけるイメージ関数値と同一となる。このような結果は、画素値がイメージ関数の局部的な特長のみで決定されるので、ある確率特性を有する。図４Ｃに示すように各画素のサイズ内の少なくとも２つのイメージ関数値に基づいて画素値を決定することで、従来の方法を置き換えることができる。
【００４７】
図４Ｃによれば、画素値１、２、３は、例えば、各画素の画素サイズ内のステップ関数を備えるフィルタ関数でイメージ関数ｆを畳み込むことによって決定される。その結果、画素値１、２、３は、各画素１１、１２、１３内のイメージ関数値に関して改良された平均近似値を示す。デジタル化した画素値の算出は同様の方法で２つ又はそれ以上の画像寸法において得られる。
２．イメージ関数の再構成の高速化（高速ＯＰＥＤ）
２．１数学的背景
【００４８】
各画素においてイメージ関数をさらに統合するには余分な処理時間が必要となる。この余分な処理時間は、好ましくは本発明に従い、以下に記載する改良ＯＰＥＤアルゴリズム（いわゆる高速ＯＰＥＤ）を用いてイメージ関数を再構成することで補償される。
【００４９】
欧州特許出願公開第０４０３１０４３．５号明細書（特許文献１）に記載されるＯＰＥＤアルゴリズムによれば、イメージ関数の近似値は、測定した投射値から以下の二重和Ａ２ｍとして直接算出できる。
【００５０】
【数９】
【００５１】
上式においては、各々以下の通りである（Ｕｋは第２種のチェビシェフ多項式である）。
【００５２】
【数１０】
【００５３】
【数１１】
【００５４】
【数１２】
【００５５】
Ａ２ｍの構造は直接的な方法で一度だけ高速フーリエ変換（ＦＦＴ）を実行するために使われる。以下の定義式を用いて、Ｓｋ，νの値を離散正弦変換のためのＦＦＴで求めることができる。
【００５６】
【数１３】
【００５７】
したがって、項Ａ２ｍは次式のように書き直せる。
【００５８】
【数１４】
【００５９】
したがって、ＯＰＥＤアルゴリズムの主要な工程は上記の二重和を求めることにあり、逆射影の工程として考えることができる。Ｎ＝２ｍ＋１とすれば、マトリクスＳｋ，νを求めるにはＯ（Ｎ２ｌｏｇＮ）の演算（フロップス「flops」)を必要とする。二重和を求めるにはＯ（Ｎ２）の演算を必要とする。ここで、Ｍ×Ｍの格子を求めるにはＯ（Ｎ２（Ｍ２＋ｌｏｇＮ）が必要となる。特にＭ〜Ｎであれば、Ｏ（Ｎ４）が必要となる。主たる演算に必要なのは格子点における二重和を求めることにある。換言すれば、主として必要なことは逆射影の工程にある。
【００６０】
フィルター逆射影法のアルゴリズムとは異なり、式（４）の和にはＦＦＴを介して求められる離散した畳み込みが含まれていない。ただ、式（２）におけるＵｋは、次式で表現できる。
【００６１】
【数１５】
【００６２】
【数１６】
【００６３】
内側の和（inner sum）は特定点において、離散正弦変換のためのＦＦＴで求めることができる。このことは、補間工程を導入してＦＦＴにより高速で数値を求めることが可能であることを示している。さらに詳述すれば、式（５）の内側の和の後で次式を定義できる。
【００６４】
【数１７】
【００６５】
ただし０＜θ＜π………式(6)
図５において、黒丸はθνを表わし、白丸はθ（ｘ，ｙ）を意味する。
【００６６】
離散正弦変換用のＦＦＴは、以下の計算を効率的に行うために利用できる。
【００６７】
【数１８】
すなわち、式（５）における内側の和は、θν（ｘ，ｙ）＝ξιのとき、ＦＦＴで求めることができる。
【００６８】
補間工程は、線形補間を用いるのが好ましい。所与の（ｘ，ｙ）については、θν（ｘ，ｙ）がξ１とξ１＋１の間にあって、α1,yとαl+1v 間の線形補間の値を式（５）の内側の和への近似値として用いるように、整数ｌが選択される。線形補間は、次式によって与えられる。
【００６９】
【数１９】
【００７０】
但しξ１≦θ≦ξ１＋１である。次にＯＰＥＤアルゴリズムの高速実行（図２の工程１３０）は以下のステップを含む。
【００７１】
ステップ１：各ν＝０、・・・２ｍについて、各ｋ＝０、・・・２ｍの演算に以下のＦＦＴが用いられる（ｍ：正の整数）。
【００７２】
【数２０】
【００７３】
ステップ２：各Ｌ＝０、１、・・・２ｍ−１について、次式を演算するためにＦＦＴが用いられる。
【００７４】
【数２１】
【００７５】
ステップ３：検査領域（半径ｃｏｓ（π／（２ｍ＋１））のディスク）内の各再構成点（ｘ，ｙ）について、整数ｌは次式のように決定される。
【００７６】
【数２２】
【００７７】
【数２３】
【００７８】
最後に、イメージ関数ｆは次式に従って算出される。
【００７９】
【数２４】
【００８０】
但し、ｕν＝（２ｍ＋１）θν／π−（ｌ＋１）………式(7)
【００８１】
以上から、ｓｉｎθvがアルゴリズムの最終ステップで分母に現れることが判る。但し、ｃｏｓθv（ｘ，ｙ）＝ｘｃｏｓφν＋ｙｓｉｎφν＝１のときのみθv＝０となる。これは（ｘ，ｙ）＝（ｃｏｓφν，ｓｉｎφν）のときのみ生じる。点（ｃｏｓφν，ｓｉｎφν）は領域Ｂの境界上にあるので、これらの点を再構成点とみなす必要はない。実際、検査領域は普通、単位ディスク内にあり。したがって、演算はＢ内の小さなディスク内の点に限定することができる。これは確実に、演算において有効桁の消失を起こすのに最終段階のｓｉｎθvの値が小さすぎることがないようにする。さらにまた、（ｘ，ｙ）が半径ｃｏｓξ０＝ｃｏｓπ／（２ｍ＋１）であるディスクに限定される場合、ステップ３におけるｌの選択が該ディスク内のすべての（ｘ，ｙ）に対して固有であることも保証する。
【００８２】
アルゴリズムはＦＦＴを２度用いて、ステップ３における最終和は単一の和となり、それを求めるにはＯ（Ｎ）の演算を要する。したがって、Ｍ〜Ｎを用いてＭ×Ｍの格子上で数値を求める上で必要なものはＯ（Ｎ３）であり、これはＦＢＰ（filtered back projection）アルゴリズムと同じ桁数である。
２．２ 演算と結果
【００８３】
数値演算には、パッケージソフト「ＦＦＴＷ」（http://www.fftw.org/）における離散正弦変換用のＦＦＴを用いた。数値例はいわゆるシェップ−ローガン（Ｓｈｅｐｐ−Ｌｏｇａｎ）の頭部ファントム上で行われる。この頭部ファントムは非常に優れた解析用ファントムであり、画像内のあらゆる楕円の境界近傍で画像がジャンプ（jump）を含み、境界上のジャンプも含む。画像を表す関数はステップ関数であり、これは楕円の境界において連続していない。
【００８４】
画像の再構成は、補間工程のないＯＰＥＤアルゴリズムと高速ＯＰＥＤアルゴリズムとをそれぞれ使って行い、後者は上述したように補間工程を含んでいる。いずれの場合も、Ｓｋ、νはＦＦＴで演算した。
【００８５】
再構成の実行は、インテル社製ＣＰＵ「Ｘｅｎｏｎ（商標）」を２個（各々３０６５ＭＨｚ）と４ＧＢのＲＡＭ（コードはＣ言語で書き込まれる）を備えたＣＥＬＳＩＵＳコンピュタＲ６１０型で行った。ＯＰＥＤアルゴリズムを使った場合、再構成に３４４秒要し、このうち９５％以上は逆射影工程に費やされた。これに対して高速ＯＰＥＤアルゴリズムを用いた場合は、所要時間は僅か１３秒となり、２６倍以上の速度向上が得られた。また、これら２つの画像の見た目は殆ど変わらなかった。
【００８６】
再構成の誤差を求めたところ、相対的な最小二乗平均誤差では高速ＯＰＥＤが僅かに劣り、平均誤差では僅かに勝っていた。誤差の桁数は同じであり、誤差の違いは殆ど無視できるレベルであった。
３．撮像装置
【００８７】
図６に、一実施の形態に係る撮像装置１００を概略的に示す。撮像装置１００は、エネルギー生成器２００及び検出装置３００を備えた測定装置と、この測定装置に接続した再構成回路４００とを含む。さらに保持装置５００を備えており、これは例えばＣＴシステムで知られているようなキャリアテーブルや、検査中の対象物を測定装置に配置したり、対象物の幾何学的形状をエネルギー生成器２００及び検出装置３００に対して調整したりする他の任意のキャリア又は物体保持器である。制御装置、ディスプレイ装置（図示せず）等のような他の構成部品も、これら自体は従来技術から既知の態様で設けられる。
【００８８】
エネルギー生成器２００は、従来のＣＴ装置から判るように、線源キャリア２２０（例えばガイドレール）上に配置される、例えば移動可能なＸ線管のようなエネルギー入力源２１０を含む。検出装置は、エネルギー入力源２１０とは対向する関係で線源キャリア２２０上に可動配置されるセンサ配列３１０を含む。このような構造で、（図面の平面と平行な）ＲＯＩ（Region Of Investigation）を通過する射影方向は、保持装置５００の回りで構成部品２１０、３１０の組合せを回転させることにより設定できる。線源キャリア２２０は対象物の回りでエネルギー生成器２００及び検出装置３００を回転可能にする円環として示される。変形例によれば（変更することにより）、線源キャリアは楕円形状や他の形状であってもよい。これは検査対象物の幾何学的形状に適合させるという点で特長を示すことができる。
【００８９】
測定装置２００、３００に接続した再構成回路４００は、測定したイメージ関数の値を乗じた多項式の和としてイメージ関数ｆを再構成するように適合させてある。再構成回路４００は、上記第１発明に係る画像処理方法でイメージ関数ｆを離散化させるよう構成した画像処理回路４１０を含む。画像処理回路４１０には、例えば本発明を実施するためのプログラム命令を含んだコンピュータで読み取り可能な記憶媒体を備えるマイクロプロセッサが実装される。
【００９０】
上記の詳細な説明、図面、及び特許請求の範囲で開示した本発明の特長は、単独もしくは組み合わせて本発明を様々な実施形態で実現することができる重要なものである。
【符号の説明】
【００９１】
１・・・画素値
２・・・画素値
３・・・画素値
１’・・・画素値
２’・・・画素値
３’・・・画素値
１０・・・画素配列
１１・・・画素
１２・・・画素
１３・・・画素
１００・・・撮像装置
２００・・・エネルギー生成器
２１０・・・エネルギー入力源
２２０・・・線源キャリア
３００・・・検出装置
３１０・・・センサ配列
４００・・・再構成回路
４１０・・・画像処理回路
５００・・・保持装置
ｆ・・・イメージ関数
ν・・・射影方向
【特許請求の範囲】
【請求項１】
複数の画素で構成され所定の画素サイズを有するデジタル画像を、イメージ関数(f)の近似値として生成するための画像処理方法であって、
複数の所定の射影方向(v)に応じて測定された複数の射影関数を含むラドンデータからイメージ関数(f)を取得する工程と、
所望のデジタル画像を表す画素値をイメージ関数から決定する工程と
を含み、
画素値は各画素の画素サイズ内で少なくとも２つのイメージ関数値に基づいて決定されることを特徴とする画像処理方法。
【請求項２】
請求項１に記載の画像処理方法であって、各画素の画素サイズ内でフィルタ関数を用いてイメージ関数を畳み込むことにより画素値が決定されることを特徴とする画像処理方法。
【請求項３】
請求項２に記載の画像処理方法であって、各画素の画素サイズ内でステップ関数を用いてイメージ関数を畳み込むことにより画素値を決定することを特徴とする画像処理方法。
【請求項４】
請求項１〜３のいずれか一に記載の画像処理方法であって、各画素の画素サイズ内でイメージ関数の平均値として画素値が決定されることを特徴とする画像処理方法。
【請求項５】
請求項１〜４のいずれか一に記載の画像処理方法であって、射影関数の値を乗じた多項式の和としてイメージ関数(f)が提供されることを特徴とする画像処理方法。
【請求項６】
請求項１〜５のいずれか一に記載の画像処理方法であって、
射影関数は離散した射影プロファイルを備えており、離散した各射影プロファイルは同一の射影方向(v)を有する複数の射影線(j)に対応した射影値(γ(v,j)）を備え、
イメージ関数(f)は射影和
【数２５】
（但し、
【数２６】
【数２７】
であって、Ｕｋは直交リッジ多項式）として表わされ、
イメージ関数(f)は正弦波変換法を用いて中間和
【数２８】
（但し、θv=arccos(xcosφv+ysinφv)）を決定する工程と、
射影和を
【数２９】
として決定する工程とで決定されることを特徴とする画像処理方法。
【請求項７】
請求項６に記載の画像処理方法であって、中間和が
【数３０】
ｌ＝０，１，…，２ｍ−１で決定されることを特徴とする画像処理方法。
【請求項８】
請求項６に記載の画像処理方法であって、中間和が後続値
【数３１】
ｌ＝０，１，…，２ｍ−１間の補間法によって得られる
【数３２】
で決定されることを特徴とする画像処理方法。
【請求項９】
請求項８に記載の画像処理方法であって、線形補間法で得られる
【数３３】
で中間和が決定されることを特徴とする画像処理方法。
【請求項１０】
請求項６〜９のいずれか一に記載の画像処理方法であって、高速フーリエ変換法を用いて
【数３４】
を決定する工程を含むことを特徴とする画像処理方法。
【請求項１１】
請求項１〜１０のいずれか一に記載の画像処理方法であって、
Ｘ線コンピュータ断層撮影装置、
超音波断層撮影装置、
ＰＥＴ撮像装置
ＰＥＴ撮像、
光断層撮影、
ガンマ線撮像装置
ＳＰＥＣＴ撮像装置、
中性子透過検出システム、又は
電気インピーダンス断層撮影装置
において測定されたラドンデータでイメージ関数(f)が決定されることを特徴とする画像処理方法。
【請求項１２】
対象物(1)内の検査領域(2)を撮像するための撮像方法であって、
所定の射影方向(v)において検査領域(2)を通過して複数のエネルギー入力を指向させる工程と、
複数のエネルギー入力ビームで測定された減衰値を含む射影プロファイルを決定する工程と、
射影プロファイルから導いたイメージ関数(f)を再構成する工程と、
請求項１〜１１のいずれか一に係る画像処理方法をイメージ関数(f)に適用する工程と、
を含むことを特徴とする撮像方法。
【請求項１３】
請求項１２に記載の撮像方法であって、さらに、
所望の視覚化された画像として、イメージ関数(f)のデジタル画像を表示する工程を含むことを特徴とする撮像方法。
【請求項１４】
請求項１２又は１３に記載の撮像方法であって、射影プロファイルを決定する過程が
測定装置(200, 300)に対象物(1)を配置する工程と、
対象物(1)が複数の所定の射影方向(v)に沿って指向させたエネルギー入力を受ける工程と、
射影プロファイルを測定する工程と、
を含むことを特徴とする撮像方法。
【請求項１５】
請求項１４に記載の撮像方法であって、ヘリカル状の射影データを得るために対象物(1)がエネルギー入力を受ける工程中に、少なくとも一の対象物(1)と測定装置(200, 300)を所定の方向に並進移動させることを特徴とする撮像方法。
【請求項１６】
請求項１２〜１５のいずれか一に記載の撮像方法であって、対象物(1)が
生体有機物又はその一部、
流体組成物、
固形物質、
工作物及び／又は、
セキュリティ上の理由で検査すべき対象物
を含むことを特徴とする撮像方法。
【請求項１７】
対象物(1)における検査領域(2)を撮像するための撮像装置(100)であって、該撮像装置が
複数の所定の射影方向(v)に応じて射影関数を測定するようにした測定装置(200, 300)と、
測定した射影関数に基づいてイメージ関数(f)を再構成するようにした再構成回路(400)であって、測定装置(200, 300)に接続した再構成回路(400)と、
を備え、さらに、
請求項１〜１６のいずれか一に係る方法でイメージ関数(f)を処理するようにした画像処理回路(410)も備えることを特徴とする撮像装置。
【請求項１８】
請求項１７に記載の撮像装置であって、再構成回路(400)は、測定された射影関数の値を乗じた多項式の和としてイメージ関数(f)を再構成するようにしたことを特徴とする撮像装置。
【請求項１９】
多項式の和を含む電子的に読み取り可能なデータを有するデジタル保存用媒体又はコンピュータプログラム製品であって、請求項１〜１３のいずれか一に係る方法を実行するための撮像装置(100)内で前記データが演算単位と相互作用可能であることを特徴とする媒体又はコンピュータプログラム製品。
【請求項２０】
請求項１〜１３のいずれか一に係る方法を実行するためのプログラムコードを有し、コンピュータ読み取り可能な媒体に存在するコンピュータプログラム。
【請求項２１】
請求項１〜１３のいずれか一に係る方法を実行するためのプログラム命令を含むコンピュータ読み取り可能な記憶媒体を備える装置。
【請求項１】
複数の画素で構成され所定の画素サイズを有するデジタル画像を、イメージ関数(f)の近似値として生成するための画像処理方法であって、
複数の所定の射影方向(v)に応じて測定された複数の射影関数を含むラドンデータからイメージ関数(f)を取得する工程と、
所望のデジタル画像を表す画素値をイメージ関数から決定する工程と
を含み、
画素値は各画素の画素サイズ内で少なくとも２つのイメージ関数値に基づいて決定されることを特徴とする画像処理方法。
【請求項２】
請求項１に記載の画像処理方法であって、各画素の画素サイズ内でフィルタ関数を用いてイメージ関数を畳み込むことにより画素値が決定されることを特徴とする画像処理方法。
【請求項３】
請求項２に記載の画像処理方法であって、各画素の画素サイズ内でステップ関数を用いてイメージ関数を畳み込むことにより画素値を決定することを特徴とする画像処理方法。
【請求項４】
請求項１〜３のいずれか一に記載の画像処理方法であって、各画素の画素サイズ内でイメージ関数の平均値として画素値が決定されることを特徴とする画像処理方法。
【請求項５】
請求項１〜４のいずれか一に記載の画像処理方法であって、射影関数の値を乗じた多項式の和としてイメージ関数(f)が提供されることを特徴とする画像処理方法。
【請求項６】
請求項１〜５のいずれか一に記載の画像処理方法であって、
射影関数は離散した射影プロファイルを備えており、離散した各射影プロファイルは同一の射影方向(v)を有する複数の射影線(j)に対応した射影値(γ(v,j)）を備え、
イメージ関数(f)は射影和
【数２５】
（但し、
【数２６】
【数２７】
であって、Ｕｋは直交リッジ多項式）として表わされ、
イメージ関数(f)は正弦波変換法を用いて中間和
【数２８】
（但し、θv=arccos(xcosφv+ysinφv)）を決定する工程と、
射影和を
【数２９】
として決定する工程とで決定されることを特徴とする画像処理方法。
【請求項７】
請求項６に記載の画像処理方法であって、中間和が
【数３０】
ｌ＝０，１，…，２ｍ−１で決定されることを特徴とする画像処理方法。
【請求項８】
請求項６に記載の画像処理方法であって、中間和が後続値
【数３１】
ｌ＝０，１，…，２ｍ−１間の補間法によって得られる
【数３２】
で決定されることを特徴とする画像処理方法。
【請求項９】
請求項８に記載の画像処理方法であって、線形補間法で得られる
【数３３】
で中間和が決定されることを特徴とする画像処理方法。
【請求項１０】
請求項６〜９のいずれか一に記載の画像処理方法であって、高速フーリエ変換法を用いて
【数３４】
を決定する工程を含むことを特徴とする画像処理方法。
【請求項１１】
請求項１〜１０のいずれか一に記載の画像処理方法であって、
Ｘ線コンピュータ断層撮影装置、
超音波断層撮影装置、
ＰＥＴ撮像装置
ＰＥＴ撮像、
光断層撮影、
ガンマ線撮像装置
ＳＰＥＣＴ撮像装置、
中性子透過検出システム、又は
電気インピーダンス断層撮影装置
において測定されたラドンデータでイメージ関数(f)が決定されることを特徴とする画像処理方法。
【請求項１２】
対象物(1)内の検査領域(2)を撮像するための撮像方法であって、
所定の射影方向(v)において検査領域(2)を通過して複数のエネルギー入力を指向させる工程と、
複数のエネルギー入力ビームで測定された減衰値を含む射影プロファイルを決定する工程と、
射影プロファイルから導いたイメージ関数(f)を再構成する工程と、
請求項１〜１１のいずれか一に係る画像処理方法をイメージ関数(f)に適用する工程と、
を含むことを特徴とする撮像方法。
【請求項１３】
請求項１２に記載の撮像方法であって、さらに、
所望の視覚化された画像として、イメージ関数(f)のデジタル画像を表示する工程を含むことを特徴とする撮像方法。
【請求項１４】
請求項１２又は１３に記載の撮像方法であって、射影プロファイルを決定する過程が
測定装置(200, 300)に対象物(1)を配置する工程と、
対象物(1)が複数の所定の射影方向(v)に沿って指向させたエネルギー入力を受ける工程と、
射影プロファイルを測定する工程と、
を含むことを特徴とする撮像方法。
【請求項１５】
請求項１４に記載の撮像方法であって、ヘリカル状の射影データを得るために対象物(1)がエネルギー入力を受ける工程中に、少なくとも一の対象物(1)と測定装置(200, 300)を所定の方向に並進移動させることを特徴とする撮像方法。
【請求項１６】
請求項１２〜１５のいずれか一に記載の撮像方法であって、対象物(1)が
生体有機物又はその一部、
流体組成物、
固形物質、
工作物及び／又は、
セキュリティ上の理由で検査すべき対象物
を含むことを特徴とする撮像方法。
【請求項１７】
対象物(1)における検査領域(2)を撮像するための撮像装置(100)であって、該撮像装置が
複数の所定の射影方向(v)に応じて射影関数を測定するようにした測定装置(200, 300)と、
測定した射影関数に基づいてイメージ関数(f)を再構成するようにした再構成回路(400)であって、測定装置(200, 300)に接続した再構成回路(400)と、
を備え、さらに、
請求項１〜１６のいずれか一に係る方法でイメージ関数(f)を処理するようにした画像処理回路(410)も備えることを特徴とする撮像装置。
【請求項１８】
請求項１７に記載の撮像装置であって、再構成回路(400)は、測定された射影関数の値を乗じた多項式の和としてイメージ関数(f)を再構成するようにしたことを特徴とする撮像装置。
【請求項１９】
多項式の和を含む電子的に読み取り可能なデータを有するデジタル保存用媒体又はコンピュータプログラム製品であって、請求項１〜１３のいずれか一に係る方法を実行するための撮像装置(100)内で前記データが演算単位と相互作用可能であることを特徴とする媒体又はコンピュータプログラム製品。
【請求項２０】
請求項１〜１３のいずれか一に係る方法を実行するためのプログラムコードを有し、コンピュータ読み取り可能な媒体に存在するコンピュータプログラム。
【請求項２１】
請求項１〜１３のいずれか一に係る方法を実行するためのプログラム命令を含むコンピュータ読み取り可能な記憶媒体を備える装置。
【図１】
【図２】
【図３−４】
【図５】
【図６】
【図２】
【図３−４】
【図５】
【図６】
【公表番号】特表２０１０−５１８５２４（Ｐ２０１０−５１８５２４Ａ）
【公表日】平成２２年５月２７日（２０１０．５．２７）
【国際特許分類】
【出願番号】特願２００９−５４９３７２（Ｐ２００９−５４９３７２）
【出願日】平成２０年２月８日（２００８．２．８）
【国際出願番号】ＰＣＴ／ＥＰ２００８／０００９８３
【国際公開番号】ＷＯ２００８／０９８７１１
【国際公開日】平成２０年８月２１日（２００８．８．２１）
【出願人】（３９８０６１２４５）ヘルムホルツ・ツェントルム・ミュンヒェン・ドイチェス・フォルシュンクスツェントルム・フューア・ゲズントハイト・ウント・ウムベルト（ゲーエムベーハー） (19)
【氏名又は名称原語表記】Ｈｅｌｍｈｏｌｔｚ Ｚｅｎｔｒｕｍ Ｍｕｅｎｃｈｅｎ Ｄｅｕｔｓｃｈｅｓ Ｆｏｒｓｃｈｕｎｇｓｚｅｎｔｒｕｍ ｆｕｅｒ Ｇｅｓｕｎｄｈｅｉｔ ｕｎｄ Ｕｍｗｅｌｔ （ＧｍｂＨ）
【出願人】（５０７３８６２８７）ザ・ステイト・オブ・オレゴン・アクティング・バイ・アンド・スルー・ザ・ステイト・ボード・オブ・ハイアー・エドゥケイション・オン・ビハーフ・オブ・ザ・ユニバーシティー・オブ・オレゴン (3)
【Ｆターム（参考）】
【公表日】平成２２年５月２７日（２０１０．５．２７）
【国際特許分類】
【出願日】平成２０年２月８日（２００８．２．８）
【国際出願番号】ＰＣＴ／ＥＰ２００８／０００９８３
【国際公開番号】ＷＯ２００８／０９８７１１
【国際公開日】平成２０年８月２１日（２００８．８．２１）
【出願人】（３９８０６１２４５）ヘルムホルツ・ツェントルム・ミュンヒェン・ドイチェス・フォルシュンクスツェントルム・フューア・ゲズントハイト・ウント・ウムベルト（ゲーエムベーハー） (19)
【氏名又は名称原語表記】Ｈｅｌｍｈｏｌｔｚ Ｚｅｎｔｒｕｍ Ｍｕｅｎｃｈｅｎ Ｄｅｕｔｓｃｈｅｓ Ｆｏｒｓｃｈｕｎｇｓｚｅｎｔｒｕｍ ｆｕｅｒ Ｇｅｓｕｎｄｈｅｉｔ ｕｎｄ Ｕｍｗｅｌｔ （ＧｍｂＨ）
【出願人】（５０７３８６２８７）ザ・ステイト・オブ・オレゴン・アクティング・バイ・アンド・スルー・ザ・ステイト・ボード・オブ・ハイアー・エドゥケイション・オン・ビハーフ・オブ・ザ・ユニバーシティー・オブ・オレゴン (3)
【Ｆターム（参考）】
[ Back to top ]