JP2016534482A - 二次尤度汎関数を用いてデータの統計モデリングを行う方法およびシステム - Google Patents

二次尤度汎関数を用いてデータの統計モデリングを行う方法およびシステム Download PDF

Info

Publication number
JP2016534482A
JP2016534482A JP2016543996A JP2016543996A JP2016534482A JP 2016534482 A JP2016534482 A JP 2016534482A JP 2016543996 A JP2016543996 A JP 2016543996A JP 2016543996 A JP2016543996 A JP 2016543996A JP 2016534482 A JP2016534482 A JP 2016534482A
Authority
JP
Japan
Prior art keywords
qlf
data
parameters
parameter
computer
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
JP2016543996A
Other languages
English (en)
Other versions
JP6453892B2 (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.)
Imagerecon LLC
Original Assignee
Imagerecon LLC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Imagerecon LLC filed Critical Imagerecon LLC
Publication of JP2016534482A publication Critical patent/JP2016534482A/ja
Application granted granted Critical
Publication of JP6453892B2 publication Critical patent/JP6453892B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10104Positron emission tomography [PET]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10108Single photon emission computed tomography [SPECT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10116X-ray image
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/424Iterative

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Complex Calculations (AREA)
  • Nuclear Medicine (AREA)

Abstract

対象物体を記述するデータおよび複数のパラメータを含んでいてノイズ部分を有している入力信号をソースから受信し、初期パラメータのグループを選択し、二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数(pdf)を推定し、二次尤度汎関数(QLF)を計算し、当該QLFに基づいて当該データに対する当該初期パラメータの適合度を評価し、所定の条件が満たされるまで新たなパラメータのグループを選択して当該新たなパラメータのグループの適合度を評価することにより当該QLFを反復的に最適化することによりコンピュータプロセッサ内に対象物体のモデルを構築する方法およびシステムを提供する。受容可能な適合度に達したならば、最適化されたパラメータを用いて構築された対象物体のモデルの出力を表示することができる。【選択図】図3

Description

関連出願の相互参照
本出願は、2013年9月18日出願の米国仮特許出願第61/879,629号明細書を優先権主張するものである。
本発明は、対数尤度関数を代替すべくデータの予測、特に推定された係数の個数がデータ点の個数と同等かまたはより大きい非漸近的統計的推定問題において統計的推定を行う方法、より具体的には二次尤度汎関数に関する。
検出画像、受信信号、世論調査、検索、または他の任意のデータ収集方法からの測定データはノイズまたは不確実性の影響下にあり得る。当該方法は、データを過度に解釈することなく、なるべく多くの情報をデータから抽出する能力を向上させる必要がある。
統計的推定は、ノイズが多いデータ、すなわちデータの実現値毎に変化する未知のランダムな成分を有し、且つ統計的にしか特徴付けできないデータから情報を抽出する必要がある場合に用いられる。統計的推定の目的は、データ内の信号、すなわちデータの実現値毎に変化しない再現可能な成分をモデリングすることである。信号は、データに対するモデルの適合度により値が決定されるパラメータの関数としてモデリングされる。
ノイズが多いデータをモデリングする処理には広範な用途があり、いくつか例を挙げれば対象追跡、信号処理、画像処理(CT、SPECT、PET、X線等の医療画像処理を含む)、市場調査、サプライチェーン管理、在庫管理、および金融市場が含まれるがこれらに限定されない。履歴データから抽出された情報は往々にして将来的な動向および/または付随するリスクを予測するために用いられる。基本的なモデルが正しい限り、予測および/またはリスク解析の品質は推定の精度に応じて決定される。
統計的推定が広義の最適化の分野に関係し、これには非統計的な方法が含まれ、ノイズが多いデータが関与せず、最適化は単に最良のパラメータの発見する問題に過ぎないことを指摘せねばならない。
統計的推定の起源は、最小二乗(ガウス、1809年)法まで2世紀遡り、後に最大尤度(ML)方(フィッシャー1912、1922年)に発展した。統計的に独立したデータの確率分布関数(pdf)を与えられたならば、ML法は、モデルが前提として、当該データの条件付き確率を最大化するかまたは、同等的に、対数尤度関数(LLF)を最小化する。
Figure 2016534482
ここに、xはn次元空間であっていずれの次元も連続的、離散的、更には圏論的であってよく、θはモデルパラメータ、xは観測位置、f(x,θ)はpdf、δはn次元のディラック(1958年)デルタ関数である。積分記号は、連続的な次元についての積分、および離散的および圏論的次元にわたる和を表し、但し実用的には積分を和で近似してよいものと理解されたい。
多くの用途において、pdfの積分的正規化は典型的には単位元に固定されている。この場合、式(1)の右辺第1項は定数であって省略可能であり、次式が得られる。
Figure 2016534482
追加的な項が式(1)に含まれているのは、例えば検出器により観測された事象の発生率または製品の販売率等、観測値の正規化自体が関心対象の場合を考慮に入れるためである。その場合、式(1)はポアソン(1837年)分布の非ビン化LLFである。
ML法には3個の明確な利点がある。
1.データのビン化を必要とせずに確率分布関数(pdf)を直接推定する。
2.データ点の個数がパラメータの個数を大幅に上回る漸近的極限において、ML法により推定パラメータの分散は、競合統計値の分散以下である。当該漸近的極限において、推定パラメータの分散は、その推定に用いたデータ点の個数に反比例する。所定の精度に対して、ML法は従って、他の方法よりも小さいサンプルからパラメータを推定することができるため、サンプリング効率がより高い。代替的な推定量の効率は実際、ML分散と他の推定量の分散との比として定義される。これにより、ML推定量の効率が定義から単位元に設定される一方、競合者の効率は単位元以下の分数に設定される。
3.推定パラメータの不確実性の共分散行列は、情報行列(フィッシャー、1922年)からの漸近的極限において容易に計算され、LLFの二次偏導関数のヘッセン行列は最小値をとる。
非漸近領域において、パラメータの個数がデータ点の個数と同等か上回る場合、モデルがランダムな統計的ノイズを再現可能な信号として扱うことを回避すべく解を制約する必要がある(例:Puetter、Gosnell&Yahil、2005参照)。一般的な方式は、信号を基底関数の一般的な線形「ノンパラメトリック」結合として表してその係数を推定するものである。(基底関数を特徴付ける更なる非線形パラメータがあってもよい)。推定により、有意な係数値を得ると共に、有意でない係数をゼロにするかまたは少なくとも最小化することにより制約することが目標である。このようにして、信号をノイズから分離することが望まれる。
最も信頼性が高いパラメータ化は最も保守的なものであり、最小限の複雑さまたはオッカムの剃刀としても知られる、入力データと整合する最も単純な基本的パラメータを求めるものである。単純性は文脈依存であるが、大多数の連続的なアプリケーションの場合、最も単純な解は最も滑らかなものである。PIXON(登録商標)法は、データにより許容される最大且つ空間的に適合可能な平滑化を課すことによりこの解を与える(本明細書に引用するPinya&Puetter、1993;Puetter&Yahil 1999;Puetter et al2005;米国特許第5,912,993号、第6,353,688号、第6,490,374号、第6,895,125号、第6,993,204号、第8,014,580号、第8,086,011号、第8,090,179号、第8,160,340号、第8,396,313号;米国公開特許第2012/0263393号を参照)。ALGEBRON(商標)法は、金融システムにおける測およびリスク査定(米国特許第7,328,182号)等、連続的な空間に制約されない離散的な問題に対して設計された同等な技術である。
フィッシャー(1912、1922年)の先駆的な業績以来、統計的推定の一般的傾向はMLおよびそのLLF推定量を利用するものであった。しかし、MLには多くの重大な短所があって、その有用性が制約される場合がある。
1.MLは漸近的に効率的なだけである。非漸近領域の解に追加的な制約が適用された場合、ML法の効率はもはや保証されない。
2.パラメータの共分散行列は、非漸近領域の情報行列から推定できない。無制約ML法では実際、信号として扱われるノイズを増幅することにより顕著なアーチファクトが生じることが多い。解を制約することで上述のようなアーチファクトを減らすことができるが、残余の精度は従って情報行列ではなく制約条件により決定される度合が高い。
3.LLFは一般に、パラメータθの二次関数ではなく、パラメータの決定に計算の労力を掛けることは、特に大規模な問題の場合に、漸近的サンプリング効率が多少上がるにせよ引き合わない。
4.パラメータに関するLLFの勾配は、pdf f(x,θ)に反比例する項を有している。
Figure 2016534482
低pdf領域のデータは、恐らく異常値(赤色データ)を含み、従ってパラメータの推定に大きなバイアスおよび/または変動が生じる恐れがある。
上述の制約に鑑みて、ML法は、非漸近領域のノンパラメトリック推定の他の推定量に比べて特段の利点がない。精度および計算効率の方が、もはや存在しないMLの「利点」よりも重要である。
本発明によれば、画像または信号の再構築、データ解析等において、特にパラメータの個数がデータ点の個数と同等か上回る場合に、ノイズまたは他の不確実性を有する入力データをモデリングする方法を提供する。特に、二次尤度汎関数(QLF)は、統計的推定問題における対数尤度関数を代替すべく用いられる。(a)QLFがフィッシャー一致を示し、従って統計的推定に適していること、および(b)大多数の関心対象問題について漸近的極限における効率が対数尤度関数よりも10〜15%小さいだけであることが示されている。従って、対数尤度関数が、非漸近領域において本質的な利点を有しておらず、当該分野でQLFの方が精度および計算効率並びに低確率事象に対する低感度により理想的な推定量であることが指摘される。以下の詳細な記述により、データ空間における関数の効率的な積分を含む、QLFを構築および最適化する明示的なプロシージャを提供する。この新技術の潜在的用途の例を示す。
LLFと同様に、QLFはデータのビン化を必要とせず、二次形式であることで解の計算および制約の適用が簡素化され、低pdf領域においてデータに対する感度が低い。漸近的極限においても、一般に用いられる基底関数に対するサンプリング効率においてQLFはLLFに大きく見劣りしない。線形非パラメータpdfを考慮する場合、QLFの利点が明らかになる。
本開示は統計的推定に焦点を合わせているが、本明細書に示すいくつかの技術が一般的な最適化問題に適用できることは当業者には自明であろう。
本発明の一態様において、対象物体を記述するデータおよび複数のパラメータを含んでいると共にノイズ部を有する入力信号をソースから受信し、初期パラメータのグループを選択し、二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数(pdf)を推定し、二次尤度汎関数(QLF)を計算し、QLFに基づいてデータに対する初期パラメータの適合度を評価し、所定の条件が満たされるまで新たなパラメータのグループを選択して当該新たなパラメータのグループの適合度を評価することによりQLFを反復的に最適化することにより、コンピュータプロセッサ内に対象物体のモデルを構築する方法およびシステムを提供する。受容可能な適合が達成されたならば、最適化されたパラメータを用いて構築された対象物体のモデルの出力を表示することができる。
本発明の別の態様において、データ成分およびノイズ成分を有する入力信号から対象物体の再構築された画像を生成する方法は、コンピュータに、複数のパラメータを含む当該入力信号を画像ソースから受信し、複数のパラメータから初期パラメータのグループを選択し、二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数(pdf)を推定し、
Figure 2016534482
の形式(θはパラメータを表し、xは観測値を表し、f(x、θ)はpdf)で二次尤度汎関数(QLF)を計算し、当該QLFに基づいて当該データに対する当該初期パラメータの適合度を評価し、所定の条件が満たされるまで新たなパラメータのグループを選択して当該新たなパラメータのグループの適合度を評価することによりQLFを反復的に最適化して、当該最適化されたパラメータに基づいて当該対象物体の再構築された画像の表示を含む出力を生成する命令を実行させるステップを含んでいる。いくつかの実施形態において、入力信号は第1の平面画像データおよび第2の平面画像データを含み、QLFを選択、推定、計算、評価して反復的に最適化する各ステップは第1の平面画像データおよび第2の平面画像データの各々について実行され、出力を生成するステップは対象物体の3次元画像を表示するステップを含んでいる。
本発明の更に別の態様において 入力信号に含まれる対象物体を記述するデータをモデリングするシステムは、コンピュータ可読媒体と、当該コンピュータ可読媒体に接続されたパラメータ最適化プロセッサと、当該パラメータ最適化プロセッサに接続されていて前記パラメータ最適化プロセッサとの間で再構築されたモデルの電子表現を受送信すべく適合された通信インターフェースとを含み、当該コンピュータ可読媒体が、前記パラメータ最適化プロセッサにより実行されたならば前記パラメータ最適化プロセッサに以下の動作、すなわち被写体データを収集すべく構成されたソースから入力信号を受信し、当該被写体データに対応する初期パラメータのグループを生成し、二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数を推定し、
Figure 2016534482
の形式(θはパラメータを表し、xは観測値を表し、f(x、θ)はpdf)で二次尤度汎関数(QLF)を計算し、当該QLFに基づいて当該データに対する当該初期パラメータの適合度を評価し、所定の条件が満たされるまで新たなパラメータのグループを選択して当該新たなパラメータのグループの適合度を評価することによりQLFを反復的に最適化して、最適化されたパラメータを用いて構築された当該対象物体のモデルを含む出力を生成するステップを含む動作を実行させるソフトウェア命令を保存している。いくつかの実施形態において、当該データは重みwを含み、QLFは、
Figure 2016534482
の形式を有している。
本システムは、データおよび基底関数を用いてソース項を計算するソフトウェア命令を更に含んでいてよい。QLFは、基底関数を用いるグラム行列を計算して、グラム行列、パラメータおよびソース項を組み合せてQLFを求めることにより得られる。いくつかの実施形態において、入力信号は画像データであり、出力は、グラフィカルユーザーインターフェースに表示された対象物体の2次元、3次元または4次元表現を含んでいる。画像データは、X線、CT、放射断層撮影、SPECTおよびPETからなるグループから選択されてよく、対象物体は患者の身体部位である。画像データは、出力が3次元表現を含むように少なくとも2平面で取得されてよい。いくつかの実施形態において、画像データは、出力が4次元表現を含むように少なくとも2平面で取得され、更に時間を含んでいる。
本発明の更に別の態様において、統計的推定を用いて入力データからモデルを生成する方法を提供し、改良点として、モデルを生成するためのパラメータを最適化すべく対数尤度関数(LLF)を二次尤度汎関数(QLF)で代替することが含まれる。
図1は、本発明を実施可能な汎用計算環境のブロック図である。 図2は、従来の反復的パラメータ推定処理の模式的フロー図である。 図3は、本発明による反復的パラメータ推定処理の模式的フロー図である。 図4aは、黒い点で印を付けた7個の部位を有する円形領域の通常のボロノイ分割および小円で印を付けたボロノイセルの重心と、7個の部位を有する同じ領域のCVTとを示す。図4bは、黒い点で印を付けた7個の部位を有する円形領域の通常のボロノイ分割および小円で印を付けたボロノイセルの重心と、7個の部位を有する同じ領域のCVTとを示す。 図5は、本発明の少なくとも1個の実施形態による例示的な統計的推定システムのブロック図である。 図6は、強い周期的成分、若干の白色ノイズ、および定常的な背景からなるシミュレーションされた信号を示す。 図7は、単一のピークおよび若干の白色ノイズを示す図6の信号のパワースペクトル(ゼロ周波数は省略)を示す。 図8は、図6のソースからの「遠距離で検出された」200個の光子のヒストグラムである。 図9は、図8の200個の光子のヒストグラムの周波数(ゼロ周波数は省略)を有し、周期性を示さないパワースペクトルのプロット図である。 図10は、正しい周波数で明らかな周期性を示す、検出された200個の光子のQLFパラメータ化により決定された周波数(ゼロ周波数は省略)を有するパワースペクトルのプロット図である。 図11は、実時間定位生検の実行を支援すべくX線データをモデリングする本発明のQLF法の例示的アプリケーションを示すフロー図である。 図12は、血管造影プロシージャの実行中に実時間撮像用のX線データをモデリングする本発明のQLF法の例示的アプリケーションを示すフロー図である。 図13は、関心対象領域の3次元または4次元画像を生成すべくCTまたは放射断層撮影データをモデリングする本発明のQLF法の例示的アプリケーションを示すフロー図である。
本発明の態様について述べる前に、本発明を実装できる適当なコンピュータシステム環境100(図1)を簡単に説明することは有用であろう。コンピュータシステム環境100は、適当な計算環境の一例に過ぎず、本発明の利用範囲または機能に関して何らの限定を示唆するものではない。また、計算環境100が、例示的な動作環境100に示す構成要素の任意の1個または組み合わせに関して何らかの依存性または必要性を有しているものと解釈すべきでない。
本発明は、多数の他の汎用または専用コンピュータシステム環境または構成で動作する。本発明での利用に適した公知のコンピュータシステム、環境、および/または構成の例として、パーソナルコンピュータ、サーバコンピュータ、携帯またはラップトップ機器、マルチプロセッサシステム、マイクロプロセッサ利用システム、セットトップボックス、プログラマブル家電、ネットワークPC、ミニコンピュータ、メインフレームコンピュータ、電話システム、上述のシステムまたは装置の任意のものを含む分散計算環境等が含まれるがこれらに限定されない。
本発明は、コンピュータにより実行されるプログラムモジュール等、コンピュータ実行可能な命令の一般的な文脈で記述できる。一般に、プログラムモジュールには、特定のタスクを実行する、または特定の抽象的なデータ型を実装するルーティン、プログラム、オブジェクト、コンポーネント、データ構造等が含まれる。当業者は、本明細書の記述および/または図面を以下に述べる任意の形式の計算機可読媒体に保存できるコンピュータ実行可能な命令として実装することができる。
本発明はまた、通信ネットワークを介して接続されたリモート処理装置によりタスクが実行される分散型計算環境で実施されてもよい。分散計算環境において、プログラムモジュールは、メモリストレージ装置を含むローカルおよびリモートコンピュータ記憶媒体の両方に搭載することができる。
図1を参照するに、本発明を実装する例示的なシステムは、コンピュータ110の形式で汎用計算装置を含んでいる。コンピュータ110の構成要素は、処理部120、システムメモリ130、およびシステムメモリを含む各種のシステム要素を処理部120に接続するシステムバス121を含んでいてよいがこれらに限定されない。システムバス121は、メモリバスまたはメモリコントローラ、周辺バス、および各種バスアーキテクチャの任意のものを用いるローカルバスを含む数種類のバス構造の任意のものであってよい。非限定的な例として、このようなアーキテクチャは、業界標準アーキテクチャ(ISA)バス、マイクロチャネルアーキテクチャ(MCA)バス、拡張ISA(EISA)バス、ビデオ電子標準協会(VESA)ローカルバス、およびメザニンバスとしても知られる周辺装置相互接続(PCI)バスを含んでいる。
コンピュータ110は典型的に、各種の計算機可読媒体を含んでいる。計算機可読媒体は、コンピュータ110からアクセス可能な任意の入手可能な媒体であってよく、揮発性および不揮発性媒体、着脱可能および着脱不可能媒体のいずれも含まれる。非限定的な例として、計算機可読媒体は、コンピュータ記憶媒体および通信媒体を含んでいてよい。コンピュータ記憶媒体は、計算機可読命令、データ構造、プログラムモジュールその他のデータ等の情報を保存する任意の方法または技術で実装された揮発性および不揮発性の着脱可能および着脱不可能媒体を含んでいる。コンピュータ記憶媒体は、RAM、ROM、EEPROM、フラッシュメモリその他メモリ技術、CD−ROM、デジタル多用途ディスク(DVD)その他の光ディスクストレージ、磁気カセット、磁気テープ、磁気ディスクストレージその他の磁気ストレージ装置、または所望の情報を保存すべく利用可能であってコンピュータ110からアクセス可能な他の任意の媒体を含んでいるがこれらに限定されない。通信媒体は典型的に、計算機可読命令、データ構造、プログラムモジュールその他のデータを、搬送波その他の搬送機構等の変調されたデータ信号で表現し、且つ任意の情報配信媒体を含んでいる。「変調されたデータ信号」という用語は、信号内の情報を符号化すべく自身の1個以上の特徴が設定または変更された信号を意味する。非限定的な例として、通信媒体は、有線ネットワークまたは直接的有線接続等の有線媒体、および音響、RF、赤外線その他の無線媒体等の無線媒体を含んでいる。上述の任意の要素の組み合わせもまた、計算機可読媒体の範囲内に含めるべきである。
システムメモリ130は、読出し専用メモリ(ROM)131およびランダムアクセスメモリ(RAM)132等の揮発性および/または不揮発性メモリの形式でコンピュータ記憶媒体を含んでいる。起動中等にコンピュータ110内の要素間での情報転送を支援する基本ルーティンを含む基本入出力システム133(BIOS)は典型的にはROM131に保存されている。RAM132は典型的に、処理部120から直ちにアクセス可能な、および/または現在処理部120上で動作中であるデータおよび/またはプログラムモジュールを含んでいる。非限定的な例として、図1に、オペレーティングシステム134、アプリケーションプログラム135、他のプログラムモジュール136、およびプログラムデータ137を示す。
コンピュータ110はまた、他の着脱可能/着脱不可能な揮発性/不揮発性コンピュータ記憶媒体(「非一時的コンピュータ可読媒体」)を含んでいてよい。一例に過ぎないが、図1に、着脱不可能な不揮発性磁気媒体との間で読み出しまたは書き込みを行うハードディスク装置141、着脱可能な不揮発性磁気ディスク152との間で読み出しまたは書き込みを行う磁気ディスクドライブ151、およびCDROMその他の光媒体等の着脱可能な不揮発性光ディスク156との間で読み出しまたは書き込みを行う光ディスクドライブ155を示している。例示的な動作環境で利用可能な他の着脱可能/着脱不可能な揮発性/不揮発性コンピュータ記憶媒体は、磁気テープカセット、フラッシュメモリカード、デジタル多用途ディスク、デジタルビデオ・テープ、ソリッドステートRAM、ソリッドステートROM等を含むがこれらに限定されない。ハードディスク装置141は典型的にインターフェース140等の着脱不可能なメモリインターフェースを介してシステムバス121に接続され、磁気ディスクドライブ151および光ディスクドライブ155は典型的にインターフェース150等の着脱可能なメモリインターフェースを介してシステムバス121に接続されている。
上述の、および図1に示すドライブおよび付随するコンピュータ記憶媒体は、計算機可読命令、データ構造、プログラムモジュールおよびコンピュータ110用の他のデータを保存する。例えば図1において、ハードディスク装置141を、オペレーティングシステム144、アプリケーションプログラム145、他のプログラムモジュール146、およびプログラムデータ147を保存しているものとして示す。これらの構成要素は、オペレーティングシステム134、アプリケーションプログラム135、他のプログラムモジュール136、およびプログラムデータ137と同一であっても、または異なっていてもよいことに注意されたい。オペレーティングシステム144、アプリケーションプログラム145、他のプログラムモジュール146、およびプログラムデータ147は、少なくとも異なる複製であることを示すために本明細書では異なる番号を付与されている。
ユーザーは、キーボード162等の入力装置、マイクロホン163(電話を介して与えられた入力を表場合もある)、およびマウス、トラックボールまたはタッチパッド等のポインティング装置161を介してコンピュータ110に命令および情報を入力することができる。他の入力装置(図示せず)としてジョイスティック、ゲームパッド、衛星パラボラアンテナ、スキャナ等を含んでいてよい。これらおよび他の入力装置は多くの場合、システムバスに接続されたユーザー入力インターフェース160を介して処理部120に接続されているが、パラレルポート、ゲームポートまたは汎用シリアルバス(USB)等の他のインターフェースおよびバス構造により接続されていてもよい。モニタ191その他の種類の表示装置もまた、ビデオインターフェース190等のインターフェースを介してシステムバス121に接続されている。モニタに加え、コンピュータはまた、出力周辺装置インターフェース195を介して接続可能なスピーカー197およびプリンタ196等の他の周辺出力装置を含んでいてよい。
コンピュータ110は、リモートコンピュータ180等の1個以上のリモートコンピュータへの論理接続を用いてネットワーク化された環境で動作可能である。リモートコンピュータ180は、パーソナルコンピュータ、携帯機器、サーバ、ルータ、ネットワークPC、ピア装置その他の一般的なネットワークノードであってよく、典型的にはコンピュータ110との関連で上に述べた要素の多くまたは全てを含んでいる。図1に示す論理接続は、ローカルエリアネットワーク(LAN)171および広域ネットワーク(WAN)173を含んでいるが、他のネットワークをも含んでいてよい。そのようなネットワーク環境は、オフィス、企業内コンピュータネットワーク、イントラネットおよびインターネットで普及している。
LANネットワーク環境で用いる場合、コンピュータ110はネットワークインターフェースまたはアダプタ170を介してLAN171に接続される。WANネットワーク環境で用いる場合、コンピュータ110は典型的に、インターネット等のWAN173上で通信を確立するモデム172その他の手段を含んでいる。モデム172は、内蔵または外付けのいずれにせよ、ユーザー入力インターフェース160その他の適当な機構を介してシステムバス121に接続されていてよい。ネットワーク化された環境において、コンピュータ110またはその一部との関連で示すプログラムモジュールは、リモートメモリストレージ装置に保存されていてよい。非限定的な例として、図1に、リモートアプリケーションプログラム185をリモートコンピュータ180に常駐するものとして示している。図に示すネットワーク接続が例示的であって、コンピュータ間で通信リンクを確立する他の手段が用いてもよいことが理解されよう。
図5は、少なくとも1個の実施形態による、入力画像、信号その他のデータの正確なモデルを構築すべく利用可能な例示的な統計的推定システム500の簡略ブロック図である。システム500は、プロセッサ502、オペレーティングシステム504、メモリ506、および入出力インターフェース508を含んでいる。メモリ506はQLFエンハンサ510、およびパラメータデータ512およびQLFエンハンサ510を用いるパラメータの最適化に用いるデータを保存する区画を含んでいてよい。
動作時においてプロセッサ502は、メモリ506に保存されたパラメータ推定部510を実行することができる。QLFエンハンサ510は、プロセッサにより実行されたならば本開示に従いパラメータを推定する動作をプロセッサに実行させるソフトウェア命令を含んでいてよい。QLFエンハンサ510は、オペレーティングシステム504と協働的に動作可能であって、本明細書に記述するようにパラメータデータ512を利用して、入出力インターフェース508で表示すべく得られたモデルの出力、すなわち再構築されたデータを生成することができる。再構築されたデータの例として画像、プロット、グラフ等が含まれるがこれらに限定されない。具体的な説明例を以下の第12節に示す。
以下の詳細な記述は、二次尤度汎関数が統計的に整合性を有しているため統計的推定に統計的推定に適していることを最初に示すために編集されている。対数尤度関数が非漸近領域で本質的に利点を有しておらず、二次尤度汎関数が精度および計算効率の面で理想的な推定ツールとなることを示す。データ空間における関数の効率的な積分を含む、二次尤度汎関数を構築および最適化する明示的なプロシージャを提供する。
以下の記述では、読者の便宜を図って表1に掲載する多くの略語を用いる。
Figure 2016534482
本明細書で用いる「信号」とは、コンピュータプロセッサにアップロードされて解釈されデジタルフォーマットで構成されたデータを意味する。信号の例として、デジタル画像データ、通信伝送、工業測定、市場調査データ、サプライチェーン管理データ、在庫管理、および金融市場等が含まれるがこれらに限定されず、いずれも信号内にノイズ成分が含まれており、ノイズの除去はデータの解析結果を最適化するために望ましい。
本明細書で用いる「ノイズ」とは、測定/受信設備の限界、対象物体と信号検出器/受信器との干渉(例:大気、電磁気、機械的、その他)、および非強調信号に混入または導入されて、入力に含まれるデータを正確に反映する所望の品質の出力表現を生成する能力を損なう恐れがある他のあらゆる不確実性の結果、信号に含まれる望ましくない成分を意味する。
本明細書で用いる「対象」または代替的に「対象物体」とは、データにより記述される項目を指す。データが画像データである場合、対象は、関心対象の特徴、例えば身体の内臓、地質学的構造、天体等、検出された信号を用いて情報が得られる任意の構造である。データが市場調査データである場合、対象は、当該データに関する、または当該データが収集された母集団であってよい。
本発明の技術の詳細な記述を明快にすべく章に分割する。以下に各章の目次を示す。
1.二次尤度関数(QLF)
2.フィッシャー一致の実証
3.漸近的効率
4.線形ノンパラメトリック確率密度関数
5.重み付き二次尤度汎関数
6.不均一性補正
7.バックグラウンド
8.QLFの最適化
9.前処理部
10.非線形pdf
11.数値積分
12.例
i.疎な信号検知
ii.デジタルマンモグラフィ
iii.血管造影
iv.計算断層撮影
v.放射断層撮影
1.二次尤度汎関数(QLF)
LLFには漸近的効率の利点があるが、一般に最適化が困難である。対照的に、二次尤度汎関数(QLF)は、常に二次であるがLLFではない。QLFの利点は、漸近的効率が若干低下する代償として計算が便利な点である。LLFの場合と同様に、QLFもデータのビン化を必要とせず、次式で与えられる。
Figure 2016534482
ここに、θはデータのモデル構築に用いるパラメータを表し、xは、を表す。
ノイズがガウスpdfに従うデータのQLFとLLFとの差異を強調することが重要である。QLFはpdfにおいて二次形式のであるのに対し、LLFは対数形式である。また、ガウスpdfの対数の平均値が二次形式となり得る。(これらは、ベル状のガウスpdfの標準偏差、すなわち幅の二次形式ではない)。関心対象パラメータがガウス平均である場合、ガウスLLFはこれらのパラメータの二次関数であり、ML法は最小二乗法に帰着する。対照的に、ガウスpdfのQLFは平均の指数関数である。
一方、pdfがパラメータの一次形式である、すなわち振幅が求めるパラメータである(3章および4章を参照)基底関数の線形結合である場合、pdfが二次形式であるQLFはまた、パラメータの二次形式でもある。一方、線形pdfのLLFは振幅の対数関数である。
従って、線形pdfのガウスLLFおよびQLFは共に各々のパラメータの二次関数に帰着し、類似の技術を適用してパラメータを求めることができる。しかし、基本的なノイズパターンは全く異なっていて互いに排他的である。ガウスpdfを有する問題のパラメータ推定は、LLFを用いて求めることができ、且つそのように求めるべきであるのに対し、線形pdfはQLFを用いる方が扱いやすい。
ガウスLLFと線形QLFの差異を別の仕方で述べると、LLFは何らかの物理量を計量する際に用いられるのに対し、QLFは事象の発生率を測定するために用いられる。無論、事象をビン単位でカウントし、次いでポアソンLLF推定量、すなわち式(1)を用いてもよい。実際、カウント値が高ければポアソンpdfはガウス分布を近似し、従ってカウントを測定したい物理量と見なすことができるが、このアプローチは低いカウント値において顕著に問題となる。このときQLFが非ビン化データを用いる代替的な統計的推定法を提供する。
2.フィッシャー一致
LLFおよびQLFは共にpdfの同一の完全無制約な推定を与える。式(1)および(4)をfに関して変化させることで次式が各々得られる。
Figure 2016534482
確率変数δfは任意であるため、角括弧内の2項は等しくゼロでなければならず、両式は同一の解を与える。
Figure 2016534482
完全無制約pdfは、観測毎に変化する全てのランダムな統計的変動を含むため実用面で殆ど関心対象外である。従って次に区分的定数pdfを考察する。同様の導出により、両方の推定量を最小化することで再び同一の推定pdfが得られ、ビン内の周波数密度は定数pdfの領域に対応することを示す。
Figure 2016534482
ここに、NはビンK内の観測値の個数(カウント値)、Vはその量である。
上述の考察はまた、LLFおよびQLFがフィッシャー(1922)一致、すなわち、サンプルが無限に増加するに従い正しいpdfを与えることを示す。これは、両方の推定量が正当であることを意味するため、重要な結論である。両者の間の選択は、各々の利点と短所を加味するユーザーに委ねられる。
3.漸近的効率
MLの大きな利点はその漸近的効率にある。本章は、QLFの効率がどれくらいMLを下回るのかを評価する。
本発明は、二乗可積分な基底関数の組の線形結合として表現可能なノンパラメトリックpdfに焦点を合わせる。
Figure 2016534482
ここに、Bα(x)は基底関数であり、係数θ αはpdfを特徴付ける未知パラメータθαの複素共役であり、上式の第2項は簡潔なベクトル行列表記を用いる。
pdfは、実数且つ非負でなければならないが、pdfが実数且つ非負のままであるように制約されている限り、式(9)は複素係数を有する複素基底関数、例えばフーリエ変換を許す。実基底関数の場合、パラメータは実数である。
Figure 2016534482
最も簡単な場合は、単一の実非負基底関数であり、その唯一の未知数が非負パラメータである。
f(x)=θB(x)、B(x)≧0(単一実基底関数)(11)
1個の実基底関数しか存在しないため複素数およびベクトル表記は省略しており、以下で具体的には基底関数の正規化は単位元に設定される。
Figure 2016534482
LLF推定は、
Figure 2016534482
であって、Nは観測されたカウント値である一方、QLF推定は次式で与えられる。
Figure 2016534482
最初に、LLFおよびQLF推定は共に偏りが無い点に注意されたい。所定個数の観測値について基底関数B(x)にわたり式(13)および(14)の平均を求めて、式(12)の正規化を用いれば、いずれの場合も次式が得られる。
Figure 2016534482
観測値の個数が、平均λを有するランダムなポアソン変数であれば期待値に偏りが無い。
Figure 2016534482
次に分散を見るに、LLF推定の分散は単にポアソン分散である。
Figure 2016534482
QLFの場合、式(14)の和の項を二乗し、B(x)にわたり平均を求め、次いでNにわたり平均を求める必要がある。結果的に得られる分散は次式で与えられる。
Figure 2016534482
ここで式(18)が式(12)の単位積分だけでなく任意の正規化について正しくなるように一般的な正規化を再導入する。QLF推定量の効率は従って次式で与えられる。
Figure 2016534482
また式(19)から、効率が並進およびスケーリングに対して不変である点にも注意されたい。
ドット積不等式(コーシー、1821年;シュヴァルツ、1888年)を適用すれば、式(19)から、QLF効率が必然的に単位元を上界とし、当該上界に達するのはB(x)が自身の非ゼロ領域にわたり定数である場合に限られる。また、Bが三乗可積分でなければ効率がゼロまで減少し続けることが分かる。興味深い問題は、実際に関心対象の基底関数に対して効率が単位元からどこまで低下するかである。従って次に、一般的に用いられる基底関数を考える。(簡潔を旨として議論を一次元に限定するが、より高い次元に計算を拡張することは容易である)。
べき法則基底関数に対して次式が成り立つ。
Figure 2016534482
式(19)からQLF効率が得られる。
e=(3α+1)(α+1)/(2α+1)(べき法則QLF効率)(21)
予想されるように、効率は、定数基底関数α=0に対して最大値に達し、α→−1/3に近づくにつれてB(x)にわたる積分が発散する場合ゼロに低下する。より現実的な線形すなわちα=1の場合、効率はe=8/9≒0.89である。任意の三角基底関数についても2個の線形部分に分割にすることで分かるように同じ効率が得られる。次式のように非ゼロの切片を有する一般的な線形基底関数の場合、
Figure 2016534482
効率は、上述の限界の間、0.89≒8/9<e<1に存在することが示される。
関心対象の他の場合では次式のガウス基底関数
B(x)∝exp(−x)(ガウス基底関数)、(23)
に対して
Figure 2016534482
および、
次式の逆放物線
Figure 2016534482
に対してe=14/15≒0.93が含まれる。
最後に、2個の確率pおよびq=1−pを有する圏論的二項の場合を考える。式(19)は次式となり、
e=(p+q/(p+q)(二項QLF効率)、(25)
p=1、q=0、p=0、q=1、およびp=0.5、q=0.5に対してe=1であり、その中間で最小値e=8/9≒0.89に低下する。
従って、実用的な基底関数に対して、QLFの漸近的効率の損失は10〜15%の範囲にあり、これは以下に詳述するようにQLFの利点を前提として許容可能であると結論付けることができる。しかし、LLFおよびQLFの効率はいずれも非漸近領域では知られていないことに注意されたい。
4.線形ノンパラメトリックpdf
QLFの利点は、複数の二乗可積分基底関数(式)(9)の組み合わせである線形ノンパラメトリックpdfを考慮することにより直ちに分かる。
Figure 2016534482
θ αに関する各推定量の式(1)および(4)の勾配をゼロに設定することで
Figure 2016534482
が得られる。
3章で議論した単一の基底関数の場合を除いてLLF勾配の式はパラメータに関して極めて非線形式に見えるが、QLF勾配は線形である。
QLFおよびその勾配を簡潔なベクトル行列表記で表すことで次式が得られる。
Figure 2016534482
ここに、Qはグラム行列
Figure 2016534482
であり、bはソース項
Figure 2016534482
である。
後述する5章のQLFの最適化は、グラム行列Qが常に正半定値であり、基底関数が線形独立ならば正定値であるという事実により極めて容易になる。また、Qの要素が有限であるために必然的に基底関数が二乗可積分でなければならない点にも注意されたい。
式(28)におけるQLFの第2の利点は、Qが使用する基底関数系、すなわち式(30)だけに依存するのに対し、全てのデータの依存性がソース項b、すなわち式(31)に現れる点である。
5.重み付き二次尤度汎関数
上述のQLFの変型例は、データ点に重み付けするものであり、各事象iは位置xおよび既知の重みwにより特徴付けられる。例えば、観測される事象は既知だが異なるエネルギーを有する光子の検出であってよく、関心対象の分布関数は堆積エネルギーの分布関数であって光子の個数の分布関数ではない。別の例において、事象は異なる価格での販売であってよく、所望の分布関数はドル建ての販売高の分布関数であって販売事象の個数の分布関数ではない。このような場合、QLFは次式のように変更できる。
Figure 2016534482
式(32)の完全無制約pdfおよび区分的定数pdfの解は、式(7)および(8)を反映し、
Figure 2016534482
推定されたpdfの整合性を同様に示す。
3章の基底関数の線形結合の解析は上と同様に続く。グラム行列Qは不変であって、ソース項は次式のようになる。
Figure 2016534482
特に、Qはここでも、使用する基底関数系だけに依存するのに対し、重み付けを含む全てのデータの依存性はソース項bに現れる。
以下の議論において、bが式(35)により与えられるものと仮定し、重み付けの無い特別なケースは次式の通りである。
≡1(重み付けの無い特別なケース)(36)
6.不均一性の補正
恐らくはピクセル毎に急速に変化する、または「無効ピクセル」を含んでいるために、検出器全体で検出確率が均一でない場合、軽微な混乱が生じる。
この問題に対処する標準的な方法は、強い一定のソースについてカウント数U(x)を測定するものである。(慣習的に、U(x)は通常、単位元を平均または中央値として正規化される。)本発明の状況に当該較正を含める最も簡単な方法は、基底関数を次式で与えられる新たな基底関数で代替するものである。
Figure 2016534482
解析はこれまで通り続く。
7.背景
「背景」はxに対して緩慢に変化するpdfの部分であるのに対し、「信号」はより急速に変化し得る。背景は往々にして異なる処理に起因するため、推定されたモデルにおける信号と背景の差異を維持する意味がある。
いくつかのアプリケーションにおいて、背景は未知のパラメータ、例えば1個の定数およびxに対して緩慢に変化する項を有する数個の基底関数により特徴付けることができる。背景基底関数は次いで推定処理に含められ、最終的には単に背景としてラベル付けされて個別に報告される。
他のアプリケーションでは、背景は既知または他のデータから推定可能であって、推定パラメータの一部ではない個別の既知項A(x)としてpdfに追加すべきである。信号および背景用の組み合わせpdfは従って次式で与えられる。
Figure 2016534482
信号が基底関数すなわち式(9)の線形結合として記述される場合、pdfは次式となる。
Figure 2016534482
定数項を除けば、推定量は
Figure 2016534482
但し
Figure 2016534482
となる。
固定された外部背景を導入した結果、単に解空間のソース項bから背景項aを減算するだけで済む。但し、aは全てのxにわたる積分すなわち式(41)で与えられるのに対し、bは観測値xについての和すなわち式(31)である点に注意されたい。
8.最適化
式(4)の一般的なQLFの最適化は、pdfがどのようにパラメータ化されるかに依存する。主な関心対象である線形ノンパラメトリックpdfの場合、QLFは凸二次形式すなわち(式)(28)であって、Qおよびbは式(30)および(31)により与えられる。解きたい勾配方程式である式(29)は従ってθの一次方程式である。しかし、殆ど全ての場合において、図2の基本フロー図200に示すように、一次方程式は反復的に解く必要がある。本フロー図は当分野で公知の従来方式の反復的処理に従うため、「従来技術」とラベル付けされている。
フロー図200に簡単に示すように、従来方式反復的処理は、モデリングしたい初期パラメータ202およびデータ203を提供する入力201から始まる。初期パラメータ201を用いてデータモデル207を構築する。ステップ204でデータ203に対するデータモデル207の適合度をパラメータ推定部により計算し、次いでステップ205で適合度の良否を評価する。適合度が所定の範囲または閾値に収まっている場合、推定パラメータは表示(ステップ208)および/またはメモリに保存(ステップ209)されてよい。指定された適合度条件に達していない場合、所望の条件が満たされるまでステップ206でパラメータを更新し、更新されたパラメータに基づいて変更されたデータモデルを構築して、適合度を評価することにより本処理を反復する。
図2に示す基本タイプの反復的スキームは通常、パラメータの推定には効果的であるが、同時に以下を含む重大な問題点があり得る。
1.大規模な問題には数百万個のデータ点および数百万個のパラメータが含まれる場合がある。従ってB、BおよびQを数兆個の要素を有する完全な行列としてメモリに保存するのは現実的でない。疎な行列表現は大抵可能であるが、多くのシステムにおいて疎な行列でも大き過ぎる。そのような場合、B、BおよびQを演算子として扱う必要があり、各行列乗算が膨大なパラメータ化を用いて新たに計算される。
2.必然的に、B、BおよびQ、または任意の関連行列の膨大な再配置を必要とするいかなる方法も用いることができない。特に、逆グラム行列Q−1またはQのコレスキー分解の計算は論外である(例:Press他、2007年)。
3.多くの問題においてQは扱い難い、すなわち、最大と最小の固有値の比の絶対値が極めて大きい。これは、最適化においてソース項bの統計的ノイズが、制御されなければ大幅に強調され、解における顕著なアーチファクトが生じることを意味する。最適化方法は従って、良好な制約機構を有していなければならない。
上述の考察は、反復的な共役勾配(CG)法の利用を示唆する(Fox、Huskey&Wilkinson、1948年;Hestenes&Stiefel、1952年)。CGには、行列を明示的に計算される必要が一切なく、B、BおよびQを演算子として適用できる利点がある。QLF勾配すなわち式(29)のQθ項は2段階のステップ、すなわちパラメータ空間からデータ空間への前方射影、およびデータ空間からパラメータ空間への後方射影Bに分解される。
推定パラメータの個数がデータ点の個数と同等か上回る非漸近的問題では、ノイズが信号として認識されて解にアーチファクトが生じるのを回避すべく追加的な制約が必要である。制約方法(例:Puetter他による考察、2005年)を、他の統計推定量に適用する場合と同様にQLF推定量に適用できるが、追加的な統計推定量が必要な場合がある。
大多数のPIXON(登録商標)法の適用では、例えば、パラメータ推定部とは異なる適合度(GOF)統計量を用いて平滑化の程度を判定する。例外として、パラメータ推定部の勾配にPIXON(登録商標)平滑化を適用するが追加的なGOF統計量を必要としない新たな技術(米国特許第8,675,983号)がある。
CG法もまた、凸集合への射影法(Biemond、LagendijkおよびMersereau、1990年)により制約を課すものである。例えば、パラメータは多くの場合負でなくなければならない。
θ≧0(非負条件)(43)
上の条件は、各反復の後で単に負の成分をゼロに切り捨て、あたかも切り捨てが生じなかったように継続することで課せられる(Puetter他、2005年)。
9.前処理部
CG法は、前処理部を追加して、次式によりQLFの勾配を代替することにより大幅に加速できる(例:GolubおよびVan Loan、1996年:Saad、2003年)。
G=P−1∇L=P−1(Qθ−b)(前処理部を有する勾配)、(44)
従って、一次方程式
Pθ=c(前処理方程式)(45)
を容易に解くことができる。結果的に得られる方法は、前処理済共役勾配(PCG)法として知られる。
知的な問題依存型前処理部は、PCG法の収束速度を大幅に向上させることができる。フェイルセーフな一般的選択は、PをQすなわち次式で要素が与えられる対角行列の対角部に設定する単純な再スケーリングである。
Figure 2016534482
スケーリング前処理部が機能するのは半定値行列Qの対角支配に起因する。
Figure 2016534482
前処理により収束が加速する程度は、Qの非対角項が式(47)内の極限からの距離およびQの条件数に依存する。
しかし、式(46)は、各々の基底関数について個別に評価する必要があり、演算子としてB α次いでBαの適用を含むため、大規模な問題で実行的に計算するのは困難であろう。代替策(AhnおよびFessler、2003年)は、式(46)を次式で代替するものである。
Figure 2016534482
計算上の利点は、式(48)における和を単位パラメータの配列への単一の前方射影として計算ができる点である。
Figure 2016534482
10.非線形の場合
式(9)における解の線形ノンパラメトリック表現の簡便性は強調して余りある。にもかかわらず、基底関数を特徴付けるために追加的な非線形パラメータを導入することが有用な場合が多い。一例として幅が調節可能な基底関数が挙げられる。この場合非線形パラメータに関してQLFの偏導関数を追加して、非線形CG最適化(例:Press他、2007年)を求めることが必要である。
11.数値積分
上の10章においてグラム行列がシステムだけに依存し、且つ全てのデータ依存性がソース項bに現れると述べた。この差異の実用性に現れるのは、Qがデータ空間にわたる空間積分すなわち式(30)であるのに対し、bはデータ点にわたる和すなわち式(31)であることである。
固定されたメッシュにわたる空間積分は、空間解像度が各次元において相応の数のメッシュ点を必要とする(「次元の呪い」)ため、問題の次元が増えるにつれて計算量が膨大になる。この難点に対処する2通りの方法がある。一つの方法は、データ空間内で適応型の階層的直交格子を形成するものである。代替的な方法は、「自然な」無構造メッシュはデータ点自体の位置により与えられるものと認識することである。データ位置のボロノイ(1908年)分割によりデータ空間を分割することができ、その内部で(通常は有限な)データボリュームVがセルに分割され、空間内の各点が自身に最も近いデータ位置周辺のセルに割り当てられる。
={x∈V||x−x|<|x−x|、∀k≠i}、i=1、...、N(ボロノイセル)(50)
ボロノイ分割に関する議論についてはPress他(2007年)を参照されたい。ボロノイセルのボリュームVへの制約は、無境界ボロノイセルを切り捨てることにより実現される。
データ空間が分割されたならば、データ点の位置でpdfを計算するだけでよい。関数F(x)にわたる空間積分は、各データ点における被積分関数の値とそれらのボロノイセルのボリュームとの積にわたる和で近似される。
Figure 2016534482
ボロノイセルは確かに、データ点の位置の統計的変動に起因して変動するボリュームを有するが、積分された関数は多数のデータ点にわたり滑らかでなければならないため、ボロノイボリュームの変動により生じた余分の分散は、計算された積分に重大な影響を及ぼさない筈である。
別の可能な方法は、分割のノードが自身の重心でもある(Du、FaberおよびGunzburger、 1999年)重心ボロノイ分割(CVT)を構築するものである。pdfは次いでCVTのノードで計算されて、より滑らかなセルボリュームにわたり積分される。ソース項bに対して必要なデータ点の位置におけるpdf値は、内挿により計算することができる。CVTの決定には多少の計算労力が必要であるが、これを行うために並列化された確率的方法(Ju、DuおよびGunzburger、2002年)およびCVTエネルギー関数のC平滑度(Liu他、2009年)を利用するより新たなニュートン型最適化がある。
本明細書に示す全てのアプリケーションの中心に、反復的パラメータ推定におけるQLFの計算があり、その基本ステップを図2に示す。4〜11章で述べたように、当該計算の本質を図3のフロー図に示す。データと重み301、および基底関数302を用いてソース項304およびグラム行列305を計算し、次いでこれらをパラメータ303と組み合わせてQLF307を与える。グラム行列307は通常、特に大規模な問題においては明示的に計算されないが、基底関数302およびそれらの共役はステップ306で演算子として適用される。追加的なオプションは、11章に詳述するようにデータのボロノイ分割309に基づいて、データ空間にわたる積分をボロノイ積分308で近似するものである。このオプションのステップを、当該フロー図の破線接続で示す。4〜11章に詳述するいくつかの他のオプション的特徴は、明快さのため省略しているが、当業者には明らかな筈である。
図4a、4bは通常のボロノイ分割をCVTと比較しており、図4aは黒点で印を付けた7個の部位を有する円形領域の通常のボロノイ分割、および小さい(開いた)円で印を付けたボロノイセルの重心を示し、図4bは7個の部位を有する同じ領域のCVTを示す。
12.事例
以下に、統計的推定に本発明のアルゴリズムを適用した例を示す。
事例i:疎な信号検出
強い信号であっても、過度に疎にサンプリングされたならば検出が困難な場合がある。例えば、軍事信号は、少数の光子しか検出されていなくても遠距離から観測することができる。図6に、強い周期的信号に加えて非ゼロの定数背景を上回る若干の白色ノイズを発するソースのシミュレーションを示す。図7に、異なる周波数(ゼロ周波数は省略)の相対的な強度を示す図7の信号のパワースペクトル(Press他、2007年)を示す。
検出器が「ソースから非常に遠い」ため、比較的少数の光子、例えば200個の光子が検出される。図8の光子の到着時間のヒストグラムは若干の高周波挙動を示すが、図9のヒストグラムのパワースペクトルに周期性は検出されない。一方、QLF推定は図10に示すように、正しい周波数でパワースペクトルに周期的信号を明瞭に検出する。
事例ii:デジタルマンモグラフィ
デジタルマンモグラフィは基本的に従来のマンモグラフィと同じX線マンモグラフィシステムを用いるが、当該システムはフィルムカセットの代わりにデジタルレセプタおよびコンピュータを備えている。標準的マンモグラフィにおいて、画像はX線カセットを用いてフィルムに記録され、当該フィルムは放射線科医により光ボックスを用いて確認されてから施設のアーカイブのジャケットに保存される。デジタルマンモグラフィにおいて、胸部画像は特別な電子X線検出器を用いて撮像され、当該検出器により当該画像をデジタル画像に変換してコンピュータモニタで確認してデジタルアーカイブシステムに保存する。デジタルマンモグラフィにおいて、画像の拡大、向き、明度、およびコントラストを変更して放射線科医が特定の領域をより明確に視認できるよう支援することができる。
デジタル全視野マンモグラフィは乳癌のスクリーニングおよび診断を可能にするものであり、いずれは従来のマンモグラフィを代替するであろう。デジタルスポット視野マンモグラフィは、より高速且つより正確な定位生検を可能にするものであり、異なる二方向から胸部を撮像するため、より正確に生検の3次元誘導を行う。その結果、検査時間が短縮されると共に、患者が静止姿勢をとり続けなければならない時間が大幅に短縮されるため患者の快適性および利便性が大幅に向上する。デジタルスポット視野マンモグラフィにおいて画像はデジタル的に撮像されて直ちにシステムモニタに表示されるのに対し、従来の定位生検ではマンモグラムフィルムを露光、現像、次いで視認する必要があるため、胸部生検が完了までの時間が極めて長くなる。
デジタルマンモグラフィは、標準的フィルムマンモグラフィに比べて多くの利点をもたらす潜在的可能性がある。すなわち(a)密な胸部組織と疎な胸部組織とのコントラストが向上し、(b)画像取得が高速化(1分未満)され、(c)試験時間が短縮(フィルムベースのマンモグラフィの約半分)され、(d)画像の保存が容易になり、(e)医師が画像を操作することで乳癌の診断がより正確になり、(f)マンモグラムを反復する必要なしにフィルムの露出不足/過剰を修正可能であり、(g)他の医師による遠隔診断用に画像のデジタル送信が可能になる。
QLFは、PIXON(登録商標)法等のノイズ低減技術と協働して、収集したデジタルデータをよりうまく利用してマンモグラフィ画像の鮮明さを向上させることでこれらの目標の実現を支援することができる。これにより医師は、より少ないX線カウント数で同じ画質を得ることが可能になり、より少ない線量、より短い走査時間、またはこれら二つの組み合わせが可能なる。
具体的には、QLFは勾配が線形であるため、勾配が非線形式であるポアソンLLFよりも極めて高速なノイズ低減処理を行うことができる。これは特に、医師が生検針を誘導できるように画像を実時間で見る必要がある定位生検法で重要である。図11に、ステップ1102でQLFを用いて単一平面マンモグラムを強調し(ステップ1104)、その結果をステップ1106で医師に提示する例示的なフロー図を示す。第2の平面マンモグラムの誘導(ステップ1112)を用いて任意選択の定位生検を実行し、QLF(ステップ1114)を用いて当該マンモグラムを強調して、生検の実時間誘導を行うべくステップ1116で医師に提示する。重心ボロノイ分割は、特に低カウント領域で更なるノイズ低減をもたらす。
また、より良質な画像を用いて、潜在的に異常な領域へ医師を自動的に案内する機械学習ツールを開発することができる。このようなツールは、異常の領域を見逃がす(偽陰性)ことがないが、一旦証明されたならば、試験結果の読み込み効率を大幅に向上させ、読み込みの完了までに要する時間を短縮して、診断に対する医師の信頼度が増すことを保証すべく徹底的に試験されなければならない。
QLFは、より良質な画像の生成だけでなく、基底関数の振幅の形式で画像の信頼性が高い数値表現を与えることによっても上記処理を支援することができる。経験的に、スペクトル表現として知られるこれらの振幅は多くの場合、元のピクセル画像強度よりも有用な機械学習ツール用の入力データ(より優れた機械学習言語)であることが分かっている。
事例iii:血管造影法
血管造影法は、特に動脈、静脈および心室に注目して、身体の血管および器官の内部を視覚化すべく用いる医療画像技術である。これは従来、放射線不透過性造影剤を血管内に注入してX線検出器を用いて撮像することにより行われる。血管造影法の種類に応じて、心臓の左側および動脈系を観察する場合は大抵大腿動脈を、または心臓の右側および静脈系を観察する場合は頸静脈または大腿動脈を通して血管内に進入する。ガイドワイヤおよびカテーテルシステムを用いて、ある種の造影剤(X線を吸収することで現れる)を血液に添加してX線画像で視認可能にする。今日用いられている血管造影システムは一般に、定位生検と同様に、より良質な3次元位置調整を可能にすべく二つの異なる視点からの画像を与える。
撮像するX線画像は、静止画像または動画像のいずれでもよい。心臓以外の全ての構造について、画像は通常、デジタル減算血管造影法(DSA)と呼ばれる技術を用いて撮像される。この場合画像は通常、毎秒2〜3フレームで撮像されるため、介在する放射線科医が1個または複数の血管を通る血液流を評価することができる。DSAは、造影剤で満たされた血管だけを観察できるように骨および他の器官を減ずる。減算技術は患者が静止したままでいる必要があり、拍動および動いている心臓に用いることができないため、心臓の画像はDSAを使用せずに毎秒15〜30フレームで撮像される。これら両方の技術により、介在する放射線科医または心臓外科医が血流を阻害して痛みの原因になり得る血管内部の塞がりまたは狭まり(狭窄)を観察することができる。
従来の血管造影X線検出器は、蛍光透視法、すなわち蛍光媒体を用いてX線光子を吸収し、それらのエネルギーを画像増倍管により光学光子として再放射して電荷結合素子(CCD)ビデオカメラにより検出する処理を利用することにより、画像を記録してモニタで再生することができる。当該装置は、各X線光子から大量且つ可変的な数の光学光子が生じるため、厳密には光子カウント装置ではない。また、膨大な放射線量が生じるため、患者に影響を及ぼすだけでなく、臨床医が診療を実施できる時間も制限される。(防護服を着用しているが、完全に保護される訳でないため、週次および通算診療時間を制限する必要がある)。この放射線量を減らすべく、光子カウント装置として直接X線を検出する固体X線検出器を構築する努力を現在継続中である。導入予定の光子カウント装置においてQLFは特に、画質を向上させて放射線量を減らす有効なツールであろう。
再び、QLFはノイズ低減技術と協働して、将来的な光子カウント装置において、画質を向上させるとともに放射線量を減らす特に有効なツールとなろう。その処理は、一般に二つの撮像平面を有し、医師が血管を通して導線を誘導するのと同時に処理を実時間で実行する必要がある点を除いてマンモグラフィで概説したものと同様である。図12を参照するに、ステップ1202A、1202BにおいてX線データを2平面で収集する。各平面は、ステップ1204A、1204BにおいてQLF強調を用いて別々に処理される。処理された画像はステップ1206A、1206Bにおいて医師に提示され、血管造影法処置を実行している医師を正確に誘導する。臨床医は患者の近くにあって、防護服を着用しているにもかかわらず相当量の放射線に晒されるため、放射線量を低減する潜在能力は本出願において特に重要である。放射線量が大幅に低減することで、臨床医がより長い時間患者に接触できるようになるため、当該医療を実施する方法が変わる。
事例iv:コンピュータ断層撮影
コンピュータ断層撮影(CT)は、患者または試験対象に対して侵襲的外科治療を行うことなく当該患者または試験対象の内部構造を調べられることができる医療診断および測定方法並びに試験技術を提供する。この場合、調査対象の物体の多くの射影像が様々な角度から記録され、そこから当該物体の3次元画像を計算することができる。
観測された射影像(データ)を3次元画像に変換することにより断層撮影画像が生成される。例えば、X線CT撮影において、X線ビームが物体に向けられ、当該物体内で変動する構造によりビームは様々な程度に減衰する。物体の反対側で、減衰したビームが検出器により測定される。このような射影像が、物体周辺において様々な角度で生成される。次いでこれらの射影像を用いて、断層撮影画像再構築として知られる処理により3次元画像を計算する。
CT画像再構築は当初、および今日でも依然として多く、フィルタリングされた後方射影として知られる簡単な変換により実行されている。しかし、射影測定値はノイズが多く、且つ相対ノイズレベルは減衰量に依存する。骨および特に金属等の高密度材料を通る射影は、肉体、水、または他のより低密度の材料を通る射影よりも信号対ノイズ比が低い。検出される光子の個数の大幅且つ空間的な変動に対処することで、画質を向上させるべく統計画像再構築技術の開発に至った。
適切な統計モデリングによりノイズ画像を低減することができるため、患者へのX線放射線量を減少させることが可能になる。問題は、(恐らくは非線形な)物理モデルおよび統計的ノイズモデルに従い測定に最も良く適合する画像を見出すことである。再構築は反復的であって、多くのパスを必要とする。QLFは統計画像再構築、特に個々の光子を検出する現在開発中の装置を用いるもので特に重要な役割を果たし得る。最終的な形式として、データは検出された光子のリストに編成される。各光子は、自身のエネルギーと、検出された時間および位置と共に個々に記録される。これは個々の事象を扱うべく設計されているため、QLFが理想的に適している検出の微細な量子限界である。
図13に、例えば心臓、脳、または内臓あるいは身体構造の非侵襲的体積または機能評価を支援すべく3次元または4次元モデルの生成に用いるコンピュータ断層撮影または放射断層撮影データを強調するステップを示すフロー図を示す。処理フローは一般に、図2に示すものに続く。データ(デジタル化されたCTまたはET画像データ1301)および初期画像にシステムおよび構築パラメータを加えた値1302を用いて初期pdf1307を計算する。pdf1307は、光子リストモードデータ1303(各光子のエネルギー、時間および位置)に作用してQLF1304を計算する。QLF推定値のデータ1301、1303との適合度をステップ1305で評価して適合度の良さを求める。適合度が所定の範囲または閾値に収まる場合、推定パラメータを用いて、表示(ステップ1308)および/またはメモリに保存(ステップ1309)可能な多次元画像を生成する。指定された適合条件に達しない場合、所望の条件が満たされるまで、ステップ1306で画像パラメータを更新し、更新されたパラメータに基づいて、変更されたデータモデルのpdfをステップ1307で再計算し、適合度を評価(1305)することにより処理を反復して、結果を表示(1308)および/またはメモリに保存(1309)することができる。
時間情報が利用可能であることで、第4次元として時間を含む4次元CT撮影を可能にする潜在的能力を有している。これによりCTは、画像が時間と共に変化する動的なプロセス、例えば拍動および動いている心臓を測定することができる。制約因子は、各々の短い時間間隔での光子の個数であり、QLFはリストモードデータの非ビン化再構築によりこれを支援することができる。
事例v:放射断層撮影
放射断層撮影は、CTで用いるものと同様の断層撮影技術を用いて身体の機能工程の3次元画像を生成する核医療画像技術である。違いは、ガンマ線または陽電子を放射する放射性同位元素(放射性核種と呼ばれる)が患者の血流に注入される点である。ガンマ線を放射する放射性核種は単一の光子を放射し、その撮像方法は単一光子放射コンピュータ断層撮影(SPECT、または時折SPET)として知られる。対照的に、放射される陽電子は、体内の電子により消滅して反対方向に移動する2個の光子を形成し、これらは同時に検出される。当該撮像方法は、陽電子放射断層撮影(PET)として知られる。
殆どの時間、放射特性だけが関心対象であるマーカー放射性同位元素が特定のリガンドに取り付けられて放射性リガンドを生成し、これはある種の組織への化学結合特性により関心対象である。この合体により、リガンドと放射性同位元素の組み合わせ(放射性医薬品)を搬送して体内の関心対象部位と結合させることができ、次いでこれにより(同位元素の直接または間接的ガンマ線放射により)リガンド濃度を撮像することができる。
核スキャンは、CTまたは磁気共鳴撮像(MRI)走査と合わせて行われるケースが増えており、その組み合わせにより解剖学的および代謝的情報の両方(すなわち、何の構造か、および何を生化学的に行っているか)を与える。核撮像は解剖学的撮像と組み合せて最も有用であるため、最新の核スキャナを一体化された最高級の多検出器列CTスキャナ、またはより最近ではMRIと合わせて利用できる。2種類の走査の間に患者が姿勢を変えることなく2回の走査が同一セッション内で連続的に、または同時に実行できるため、2組の画像がより正確に登録され、核画像上の異常領域を、より完全にCTまたはMRI画像上の解剖学的知見と相関させることができる。これは、動いている器官または脳の外側でより一般的である、より高い解剖学的変動を有する構造の詳細な画像を表示するのに極めて有用である。より最近では、核撮像よりも高い解像度を有する解剖学的モダリティを用いて核画像再構築自体を誘導している。(従来は、CT画像の役割は、上述の解剖学的相関に加え、核減衰マップの提供に限定されていた)。
コンピュータ断層撮影(CT)の再構築に極めて類似した技術を3次元画像の生成に用いることは一般的だが、核撮像において収集されたデータセットの光子はCTよりもはるかに少ないため、再構築技術はより困難である。次いでQLF処理を用いて画質を向上させる、および/または患者に浴びせられる放射線量を減らすことができる。放射断層撮影装置の物理的構成はCTとは極めて異なるが、これは前方および後方射影演算子(8章参照)だけに影響を及ぼし、実際にCTおよび放射断層撮影と極めて良く類似している基本的再構築方法論には影響を及ぼさない。図13のフロー図は従って、システム構成およびパラメータは異なるものの、放射断層撮影にも同様に適用可能である。
上述の詳細な説明および例に基づいて、データの統計モデリングを行う本発明のQLF法の他の応用例は、当業者には容易に明らかになろう。従って、本発明の範囲は、本明細書に記述する具体例だけに限定されるものではない。
以下の参考文献のリストは、統計画像再構築の背景および全般的な最先端技術の理解に役立つ。一覧に示す参考文献は、当業者が本発明を完全に理解すると共に実施するために必要な程度に、当分野の技術レベルを示すと共に、本開示の一部として本明細書に引用している。
参考文献
References
1. Ahn S & Fessler JA 2003,“Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms”,IEEE Trans.Med.Imag.,22,613.
2. Biemond J,Lagendijk RL & Mersereau RM 1990,“Iterative Methods for Image Deblurring”,Proc.IEEE,78,856.
3. Cauchy A 1821,Oeuvres 2,III,373.
4. Dirac P 1958,Principles of Quantum Mechanics,4th ed.,(Oxford: Clarendon Press),p.58.
5. Du Q,Faber V & Gunzburger M 1999,“Centroidal Voronoi Tessellations:Applications and Algorithms”,SIAM Rev.,41,637.
6. Fisher RA 1912,“On an Absolute Criterion for Fitting Frequency Curves”,Messenger Math.,41,155.
7. Fisher RA 1922,“On the Mathematical Foundations of Theoretical Statistics”,Phil. Trans. Roy. Soc. A,222,309.
8. Fox L,Huskey HD & Wilkinson JH 1948,“Notes on the Solution of Algebraic Linear Simultaneous Equations”,Quart.J.Mech.Appl.Math.,1,149.
9. Gauss CF 1809,Theoria Motus Corporum Coelestium in sectionibus conicis solem ambientium,English translation,1963,Theory of Motion of the Heavenly Bodies(New York:Dover).
10. Golub GH & Van Loan CF 1996,Matrix Computations,3rd ed.(Baltimore:Johns Hopkins University Press).
11. Hestenes MR & Stiefel E 1952,“Methods of Conjugate Gradients for Solving Linear Systems”,J.Res.Nat.Bur.Standards,49,409.
12. Ju L,Du Q & Gunzburger M 2002,“Probabilistic Methods for Centroidal Voronoi Tessellations and their Parallel Implementations”,Parallel Computing,28,1477.
13. Liu Y,Wang W,Levy B,Sun F,Yan D−M,Lu L & Yang C 2009,“On Centroidal Voronoi Tessellation−Energy Smoothness and Fast Computation”,ACM Trans.Graphics 28,Article 101.
14.Pinya RK & Puetter RC 1993,“Bayesian Image Reconstruction:the Pixon and Optimal Image Modeling”Publ.Astron.Soc.Pacific 105,630.
15. Poisson SD 1837,Probabilite des Jugements en Matiere Criminelle et en Matiere Civile,Precedees des Regles Generales du Calcul des Probabilites (Paris: Bachelier),p.206.
16. Press WH,Teukolsky SA,Vetterling WT & Flannery BP 2007,Numerical Recipes,3rd ed.,(Cambridge:Cambridge University Press).
17. Puetter RC,Gosnell TR & Yahil A 2005,“Digital Image Reconstruction:Deblurring and Denoising”,Annu.Rev.Astron.Astrophys.,43,139.
18. Puetter RC & Yahil A 1999,“The Pixon Method of Image Reconstruction”,in Astronomical Data Analysis Software and Systems VIII,Mehringer DM,Plante RL & Roberts DA,eds.,ASP Conference Series,172,307.
19. Saad Y 2003,Iterative Methods for Sparse Linear Systems,2nd ed.(Philadelphia: SIAM).
20. Schwarz HA 1888,“Ueber ein Flaechen kleinsten Flaecheninhalts betreffendes Problem der Variationsrechnung”,Acta Societatis Scientiarum Fennicae 45,318.
21. Voronoi G 1908,“Nouvelles Applications des Parametres Continus a la Theorie des Formes Quadratiques”,Journal fuer die Reine und Angewandte Mathematik 133,97.

Claims (22)

  1. 対象物体のモデルを構築するコンピュータ実装方法であって、コンピュータに、
    対象物体を記述するデータおよび複数のパラメータを含んでいると共にノイズ部を有する入力信号をソースから受信し、
    前記複数から初期パラメータのグループを選択し、
    二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数(pdf)を推定し、
    Figure 2016534482
    の形式(θはパラメータを表し、xは観測値を表し、f(x、θ)はpdf)で二次尤度汎関数(QLF)を計算し、
    前記QLFに基づいて前記データに対する前記初期パラメータの適合度を評価し、
    所定の条件が満たされるまで新たなパラメータのグループを選択して前記新たなパラメータのグループの適合度を評価することにより前記QLFを反復的に最適化して、
    最適化されたパラメータを用いて構築された前記対象物体のモデルを含む出力を生成する命令を実行させることを特徴とする方法。
  2. 請求項1に記載のコンピュータ実装方法において、前記データが重みwを含み、前記QLFが、
    Figure 2016534482
    の形式を有していることを特徴とする方法。
  3. 請求項1に記載のコンピュータ実装方法において、前記データおよび基底関数を用いてソース項を計算するステップを更に含んでいることを特徴とする方法。
  4. 請求項3に記載のコンピュータ実装方法において、前記QLFが、
    前記基底関数を用いてグラム行列を計算し、
    前記グラム行列、前記パラメータ、および前記ソース項を組み合せて前記QLFを生成することにより得られることを特徴とする方法。
  5. 請求項1に記載のコンピュータ実装方法において、前記入力信号が画像データであり、前記出力が、グラフィカルユーザーインターフェースで表示される前記対象物体の2次元、3次元、または4次元表現を含んでいることを特徴とする方法。
  6. 請求項5に記載のコンピュータ実装方法において、前記画像データが、X線、CT、放射断層撮影、SPECT、およびPETからなるグループから選択され、前記対象物体が患者の身体部位であることを特徴とする方法。
  7. 請求項6に記載のコンピュータ実装方法において、前記画像データが少なくとも2個の平面で取得され、前記出力が3次元表現を含んでいることを特徴とする方法。
  8. 請求項6に記載のコンピュータ実装方法において、前記画像データが少なくとも2個の平面で取得されて更に時間を含み、前記出力が4次元表現を含んでいることを特徴とする方法。
  9. 入力信号に含まれる対象物体を記述するデータをモデリングするシステムであって、
    コンピュータ可読媒体と、
    前記コンピュータ可読媒体に接続されたパラメータ最適化プロセッサと、
    前記パラメータ最適化プロセッサに接続されていて前記パラメータ最適化プロセッサとの間で再構築されたモデルの電子表現を受送信すべく適合された通信インターフェースとを含み、前記コンピュータ可読媒体に、前記パラメータ最適化プロセッサにより実行されたならば前記パラメータ最適化プロセッサに以下の動作、すなわち、
    被写体データを収集すべく構成されたソースから入力信号を受信し、
    前記被写体データに対応する初期パラメータのグループを生成し、
    二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数を推定し、
    Figure 2016534482
    の形式(θはパラメータを表し、xは観測値を表し、f(x、θ)はpdf)で二次尤度汎関数(QLF)を計算し、
    前記QLFに基づいて前記データに対する当前記初期パラメータの適合度を評価し、
    所定の条件が満たされるまで新たなパラメータのグループを選択して前記新たなパラメータのグループの適合度を評価することにより前記QLFを反復的に最適化して、
    最適化されたパラメータを用いて構築された前記対象物体のモデルを含む出力を生成するステップを含む動作を実行させるソフトウェア命令を保存していることを特徴とするシステム。
  10. 請求項9に記載のシステムにおいて、前記データが重みwを含み、前記QLFが、
    Figure 2016534482
    の形式を有していることを特徴とするシステム。
  11. 請求項9に記載のシステムにおいて、前記データおよび基底関数を用いてソース項を計算するステップを更に含んでいることを特徴とするシステム。
  12. 請求項11に記載のシステムにおいて、前記QLFが、
    前記基底関数を用いてグラム行列を計算し、
    前記グラム行列、前記パラメータ、および前記ソース項を組み合せて前記QLFを生成することを特徴とするシステム。
  13. 請求項9に記載のシステムにおいて、前記入力信号が画像データであり、前記出力が、グラフィカルユーザーインターフェースで表示される前記対象物体の2次元、3次元、または4次元表現を含んでいることを特徴とするシステム。
  14. 請求項13に記載のシステムにおいて、前記画像データが、X線、CT、放射断層撮影、SPECT、PETからなるグループから選択され、前記対象物体が患者の身体部位であることを特徴とするシステム。
  15. 請求項14に記載のシステムにおいて、前記画像データが少なくとも2個の平面で取得され、前記出力が3次元表現を含んでいることを特徴とするシステム。
  16. 請求項13に記載のコンピュータ実装方法において、前記画像データが少なくとも2個の平面で取得されて更に時間を含み、前記出力が4次元表現を含んでいることを特徴とする方法。
  17. データ成分およびノイズ成分を有する入力信号から対象物体の再構築された画像を生成する方法であって、コンピュータに、
    複数のパラメータを含む前記入力信号を画像ソースから受信し、
    前記複数のパラメータから初期パラメータのグループを選択し、
    二乗可積分な基底関数の組の線形結合を含むノンパラメトリック確率分布関数(pdf)を推定し、
    Figure 2016534482
    の形式(θはパラメータを表し、xは観測値を表し、f(x、θ)はpdf)で二次尤度汎関数(QLF)を計算し、
    前記QLFに基づいて前記データに対する当該初期パラメータの適合度を評価し、
    所定の条件が満たされるまで新たなパラメータのグループを選択して前記新たなパラメータのグループの適合度を評価することにより前記QLFを反復的に最適化して、
    前記最適化されたパラメータに基づいて前記対象物体の再構築された画像の表示を含む出力を生成する命令を実行させることを特徴とする方法。
  18. 請求項17に記載の方法において、前記入力信号が第1の平面画像データおよび第2の平面画像データを含み、前記QLFを選択、推定、計算、推定して反復的に最適化する前記ステップが、前記第1の平面画像データおよび第2の平面画像データの各々について実行され、出力を生成するステップが前記対象物体の3次元画像を表示するステップを含んでいることを特徴とする方法。
  19. 請求項18に記載の方法において、前記データが重みwを含み、前記QLFが
    Figure 2016534482
    の形式を有していることを特徴とする方法。
  20. 請求項17に記載のシステムにおいて、前記データおよび基底関数を用いてソース項を計算するステップを更に含んでいることを特徴とするシステム。
  21. 請求項20に記載のシステムにおいて、前記QLFが、
    前記基底関数を用いてグラム行列を計算し、
    前記グラム行列、前記パラメータ、および前記ソース項を組み合せて前記QLFを生成することを特徴とするシステム。
  22. 統計的推定を用いて入力データからモデルを生成する改良された方法において、前記改良が、対数尤度関数(LLF)を二次尤度汎関数(QLF)で代替してモデルを生成するパラメータを最適化することを含んでいることを特徴とする方法。
JP2016543996A 2013-09-18 2014-09-18 二次尤度汎関数を用いてデータの統計モデリングを行う方法およびシステム Active JP6453892B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US201361879629P 2013-09-18 2013-09-18
US61/879,629 2013-09-18
PCT/US2014/056392 WO2015042314A1 (en) 2013-09-18 2014-09-18 Method and system for statistical modeling of data using a quadratic likelihood functional

Publications (2)

Publication Number Publication Date
JP2016534482A true JP2016534482A (ja) 2016-11-04
JP6453892B2 JP6453892B2 (ja) 2019-01-16

Family

ID=52668729

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2016543996A Active JP6453892B2 (ja) 2013-09-18 2014-09-18 二次尤度汎関数を用いてデータの統計モデリングを行う方法およびシステム

Country Status (6)

Country Link
US (1) US10068327B2 (ja)
EP (1) EP3047391B1 (ja)
JP (1) JP6453892B2 (ja)
CN (1) CN105556507B (ja)
HU (1) HUE063524T2 (ja)
WO (1) WO2015042314A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2020056769A (ja) * 2018-10-02 2020-04-09 キヤノンメディカルシステムズ株式会社 核医学診断装置、核医学画像再構成方法及び核医学画像再構成プログラム
JP2020516345A (ja) * 2017-04-05 2020-06-11 ゼネラル・エレクトリック・カンパニイ 深層学習に基づくトモグラフィ再構成

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150356056A1 (en) * 2014-06-09 2015-12-10 The Mathworks, Inc. Methods and systems for calculating joint statistical information
DE102015215938A1 (de) * 2015-08-20 2017-02-23 Siemens Healthcare Gmbh Verfahren zur lokalen Verbesserung der Bildqualität
US10338012B2 (en) * 2016-03-09 2019-07-02 Toshiba Medical Systems Corporation Photon counting detector and X-ray computed tomography (CT) apparatus
CN105787905B (zh) * 2016-03-24 2019-03-26 中国人民解放军信息工程大学 一种基于动态电流的锥束ct环状伪影校正方法
CN108268193B (zh) * 2017-01-04 2021-03-02 珠海金山办公软件有限公司 电子表格中函数对应的提示信息的显示方法及装置
EP3503016B1 (en) * 2017-12-19 2021-12-22 Nokia Technologies Oy Apparatus, method and computer program for processing a piecewise-smooth signal
CN108470209B (zh) * 2018-03-27 2021-06-04 北京工业大学 一种基于格拉姆矩阵正则化的卷积神经网可视化方法
GB201808878D0 (en) * 2018-05-31 2018-07-18 Imperial Innovations Ltd Optimisation system and method
EP3629335A1 (de) * 2018-09-26 2020-04-01 Siemens Healthcare GmbH Verfahren zum automatischen auswählen von bildgebungsparametern für bildgebungsverfahren
CN109830286B (zh) * 2019-02-13 2022-09-30 四川大学 基于非参数统计的脑功能磁共振编码能量成像方法
WO2021113728A1 (en) * 2019-12-05 2021-06-10 The Regents Of The University Of California Generating synthetic patient health data
CN114200328B (zh) * 2022-01-20 2022-09-30 中国科学技术大学 一种非高斯Lévy噪声下的锂离子电池SOC估计方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002067201A1 (en) * 2001-02-15 2002-08-29 The Regents Of The University Of Michigan Statistically reconstructing an x-ray computed tomography image with beam hardening corrector
US20080063135A1 (en) * 2006-09-08 2008-03-13 General Electric Company Method and system for generating a multi-spectral image of an object
JP2009536499A (ja) * 2006-05-05 2009-10-08 トムソン ライセンシング 2次元画像から3次元オブジェクトを再構成するシステム及び方法
US20090296984A1 (en) * 2006-05-04 2009-12-03 Yousef Wasef Nijim System and Method for Three-Dimensional Object Reconstruction from Two-Dimensional Images
JP2011152255A (ja) * 2010-01-27 2011-08-11 Hitachi Medical Corp 再構成演算装置、再構成演算方法、及びx線ct装置
WO2013008702A1 (ja) * 2011-07-08 2013-01-17 株式会社 日立メディコ 画像再構成装置及び画像再構成方法

Family Cites Families (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6353688B1 (en) 1993-06-08 2002-03-05 The Regents Of The University Of California Accelerated signal encoding and reconstruction using pixon method
US5912993A (en) 1993-06-08 1999-06-15 Regents Of The University Of Calif. Signal encoding and reconstruction using pixons
US7328182B1 (en) 1999-09-23 2008-02-05 Pixon, Llc System and method for prediction of behavior in financial systems
US6507633B1 (en) 2001-02-15 2003-01-14 The Regents Of The University Of Michigan Method for statistically reconstructing a polyenergetic X-ray computed tomography image and image reconstructor apparatus utilizing the method
US8549434B2 (en) * 2001-10-18 2013-10-01 Microsoft Corporation Method for graphical representation of a content collection
US6993204B1 (en) 2002-01-04 2006-01-31 Pixon Llc High speed signal enhancement using pixons
CA2685237A1 (en) 2007-04-25 2008-11-06 Pixon Imaging, Llc Image compression and decompression using the pixon method
WO2008147416A1 (en) * 2007-05-31 2008-12-04 General Electric Company Methods and systems to facilitate correcting gain fluctuations in image
US9647780B2 (en) * 2007-08-24 2017-05-09 Invention Science Fund I, Llc Individualizing a content presentation
US8090179B2 (en) * 2007-10-31 2012-01-03 Siemens Medical Solutions Usa, Inc. External pixon smoothing for tomographic image reconstruction technical field
US8014580B2 (en) 2007-10-31 2011-09-06 Siemens Medical Solution Usa, Inc. Determining a pixon map for image reconstruction
US8086011B2 (en) 2007-10-31 2011-12-27 Siemens Medical Solutions Usa, Inc. Reconstructing a tomographic image
US8577103B2 (en) * 2008-07-16 2013-11-05 Siemens Medical Solutions Usa, Inc. Multimodal image reconstruction
US8160340B2 (en) 2008-09-04 2012-04-17 Siemens Medical Solutions Usa, Inc. Reconstructing a tomographic image
RU2431196C1 (ru) * 2010-03-31 2011-10-10 Закрытое Акционерное Общество "Импульс" Способ определения уровня яркости в зоне интереса цифрового медицинского рентгеновского изображения
CN101894377B (zh) * 2010-06-07 2012-09-05 中国科学院计算技术研究所 三维标记点序列的跟踪方法及其系统
ES2697611T3 (es) 2011-04-15 2019-01-25 Siemens Medical Solutions Usa Inc Método para determinar un mapa de pixon en la reconstrucción de imagen iterativa
CN102833191A (zh) * 2011-06-13 2012-12-19 中兴通讯股份有限公司 一种信噪比估计方法与装置
US10242126B2 (en) * 2012-01-06 2019-03-26 Technoimaging, Llc Method of simultaneous imaging of different physical properties using joint inversion of multiple datasets
WO2013106512A1 (en) * 2012-01-10 2013-07-18 The Johns Hopkins University Methods and systems for tomographic reconstruction
US20140234009A1 (en) * 2013-02-19 2014-08-21 Elwha Llc Writing implement sensing a user's hand cleanliness
US20150066346A1 (en) * 2013-08-28 2015-03-05 Elwha LLC, a limited liability company of the State of Delaware Vehicle collision management system responsive to a situation of an occupant of an approaching vehicle

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002067201A1 (en) * 2001-02-15 2002-08-29 The Regents Of The University Of Michigan Statistically reconstructing an x-ray computed tomography image with beam hardening corrector
US20090296984A1 (en) * 2006-05-04 2009-12-03 Yousef Wasef Nijim System and Method for Three-Dimensional Object Reconstruction from Two-Dimensional Images
JP2009536499A (ja) * 2006-05-05 2009-10-08 トムソン ライセンシング 2次元画像から3次元オブジェクトを再構成するシステム及び方法
US20080063135A1 (en) * 2006-09-08 2008-03-13 General Electric Company Method and system for generating a multi-spectral image of an object
JP2008062035A (ja) * 2006-09-08 2008-03-21 General Electric Co <Ge> 対象の多重スペクトル画像を形成する方法及びシステム
JP2011152255A (ja) * 2010-01-27 2011-08-11 Hitachi Medical Corp 再構成演算装置、再構成演算方法、及びx線ct装置
WO2013008702A1 (ja) * 2011-07-08 2013-01-17 株式会社 日立メディコ 画像再構成装置及び画像再構成方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
HU, Y. ET.AL: ""Sparse reconstruction from a limited projection number of the coronary artery tree in X-ray rotatio", 2012 9TH IEEE INTERNATIONAL SYMPOSIUM ON BIOMEDICAL IMAGING(ISBI), JPN6018028367, 2 May 2012 (2012-05-02), pages 804 - 807, XP032199128, ISSN: 0003920764, DOI: 10.1109/ISBI.2012.6235670 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2020516345A (ja) * 2017-04-05 2020-06-11 ゼネラル・エレクトリック・カンパニイ 深層学習に基づくトモグラフィ再構成
JP7187476B2 (ja) 2017-04-05 2022-12-12 ゼネラル・エレクトリック・カンパニイ 深層学習に基づくトモグラフィ再構成
JP2020056769A (ja) * 2018-10-02 2020-04-09 キヤノンメディカルシステムズ株式会社 核医学診断装置、核医学画像再構成方法及び核医学画像再構成プログラム
JP7237621B2 (ja) 2018-10-02 2023-03-13 キヤノンメディカルシステムズ株式会社 核医学診断装置、核医学画像再構成方法及び核医学画像再構成プログラム

Also Published As

Publication number Publication date
CN105556507A (zh) 2016-05-04
HUE063524T2 (hu) 2024-01-28
EP3047391B1 (en) 2023-06-28
US10068327B2 (en) 2018-09-04
CN105556507B (zh) 2020-05-19
US20150081262A1 (en) 2015-03-19
JP6453892B2 (ja) 2019-01-16
EP3047391A4 (en) 2017-08-02
WO2015042314A1 (en) 2015-03-26
EP3047391A1 (en) 2016-07-27

Similar Documents

Publication Publication Date Title
JP6453892B2 (ja) 二次尤度汎関数を用いてデータの統計モデリングを行う方法およびシステム
Gong et al. Iterative PET image reconstruction using convolutional neural network representation
Kaplan et al. Full-dose PET image estimation from low-dose PET image using deep learning: a pilot study
Han et al. Algorithm-enabled low-dose micro-CT imaging
Geyer et al. State of the art: iterative CT reconstruction techniques
Zhang et al. Regularization strategies in statistical image reconstruction of low‐dose x‐ray CT: A review
US10282820B2 (en) Structure propagation restoration for spectral CT
US8055050B2 (en) Motion compensation in energy-sensitive computed tomography
JP6100772B2 (ja) 画像処理方法及びコンピューティング装置
CN111540025B (zh) 预测用于图像处理的图像
JP2016152916A (ja) X線コンピュータ断層撮像装置及び医用画像処理装置
CN105793894B (zh) 根据图像数据来进行骨骼分割
EP3230946B1 (en) Statistically weighted regularization in multi-contrast imaging
Li et al. Multienergy cone-beam computed tomography reconstruction with a spatial spectral nonlocal means algorithm
Roser et al. X-ray scatter estimation using deep splines
Gomi et al. Evaluation of digital tomosynthesis reconstruction algorithms used to reduce metal artifacts for arthroplasty: A phantom study
Qi et al. A quantitative study of motion estimation methods on 4D cardiac gated SPECT reconstruction
CN114503118A (zh) 通过将图像形成建模为一个或多个神经网络的图像重建
WO2016097981A1 (en) Penalized maximum likelihood material decomposition
Choi et al. Integration of 2D iteration and a 3D CNN-based model for multi-type artifact suppression in C-arm cone-beam CT
Park et al. Deep learning-based noise reduction algorithm using patch group technique in cadmium zinc telluride fusion imaging system: A Monte Carlo simulation study
Seyyedi et al. Incorporating a noise reduction technique into X-ray tensor tomography
Zubair et al. Enabling Predication of the Deep Learning Algorithms for Low-Dose CT Scan Image Denoising Models: A Systematic Literature Review
US11704795B2 (en) Quality-driven image processing
Whiteley Deep Learning in Positron Emission Tomography Image Reconstruction

Legal Events

Date Code Title Description
A711 Notification of change in applicant

Free format text: JAPANESE INTERMEDIATE CODE: A711

Effective date: 20170828

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20170828

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20170915

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20170919

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20180718

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20180724

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20181023

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20181024

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20181120

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20181120

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

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20181213

R150 Certificate of patent or registration of utility model

Ref document number: 6453892

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250