(57)【特許請求の範囲】
【請求項１】
被験者の体表の測定点から電位および磁場の少なくとも一方を計測する計測手段と、
前記計測手段により測定された前記被験者の電位および磁場の少なくとも一方に基づいて、当該被験者の心臓内の電気現象を表す情報としての、当該心臓内の興奮到達時刻の分布情報と当該心臓の心筋の電気活動の大きさを表す心筋電気活動量の分布情報とを同時に推定する分布情報推定手段と、
前記興奮到達時刻の分布情報と前記心筋電気活動量の分布情報とを表示する表示手段と、を備えたことを特徴とする心臓内電気現象の診断装置。
【請求項２】
請求項１記載の心臓内電気現象の診断装置において、
前記分布情報推定手段は、前記心臓内での前記興奮到達時刻の分布と前記心筋電気活動量の分布とを有限個の複数のパラメータで表現する表現手段と、このパラメータにより表現された前記興奮到達時刻の分布と前記心筋電気活動量の分布とに基づいて前記電位および磁場の少なくとも一方の各測定点での時間波形を計算する波形計算手段と、前記測定電位および磁場の少なくとも一方の波形と前記計算電位および磁場の少なくとも一方の時間波形とが一致するように非線型最適化法を用いて前記複数のパラメータを推定するパラメータ推定手段と、を備えることを特徴とした心臓内電気現象の診断装置。
【請求項３】
請求項２記載の心臓内電気現象の診断装置において、
前記分布情報推定手段は、前記被験者の体内の３次元画像データから前記心臓のモデルを設定する心臓モデル設定手段を備えたことを特徴とする心臓内電気現象の診断装置。
【請求項４】
請求項３記載の心臓内電気現象の診断装置において、
前記表現手段は、前記心臓内での前記興奮到達時刻の分布と前記心筋電気活動量の分布とを前記心臓モデルの代表点の数よりも少ない数で表現する手段であり、前記分布情報推定手段は、この少ない数のパラメータに基づき前記興奮到達時刻の分布と前記心筋電気活動量の分布とを算出する算出手段を備えることを特徴とした心臓内電気現象の診断装置。
【請求項５】
請求項３記載の心臓内電気現象の診断装置において、
前記表現手段は、前記心臓モデルの各代表点での前記興奮到達時刻、前記心臓内での前記心筋電気活動量で表された前記興奮到達時刻の分布、および前記心筋電気活動量の分布を、前記心臓モデルの代表点よりも少ない数のパラメータの線形関数で表現する手段であり、前記分布情報推定手段は、この線形関数で表現された少ない数のパラメータに基づき前記興奮到達時刻の分布と前記心筋電気活動量の分布とを算出する算出手段を備えることを特徴とした心臓内電気現象の診断装置。
【請求項６】
請求項５記載の心臓内電気現象の診断装置において、
前記線形関数はフーリエ変換または逆フーリエ変換による関数であることを特徴とする心臓内電気現象の診断装置。
【請求項７】
請求項６記載の心臓内電気現象の診断装置において、
前記線形関数は行列演算による関数であることを特徴とする心臓内電気現象の診断装置。
【請求項８】
請求項７記載の心臓内電気現象の診断装置において、
前記行列演算の行列は、測定した又は仮想した興奮到達時刻の共分散行列から求めた行列であることを特徴とする心臓内電気現象の診断装置。
【請求項９】
請求項３記載の心臓内電気現象の診断装置において、
前記心臓モデル設定手段は、前記心臓モデルとして心室の形状を多数の多面体の組み合わせにより設定する手段であることを特徴とする心臓内電気現象の診断装置。
【請求項１０】
請求項９記載の心臓内電気現象の診断装置において、
前記多面体は、６面体、５面体、もしくは４面体、または、それらの多面体の内の少なくとも２種類の組み合わせであることを特徴とする心臓内電気現象の診断装置。
【請求項１１】
請求項１記載の心臓内電気現象の診断装置において、
医用モダリティから得た前記被検者の体内の３次元画像データを入力し転送する画像転送手段を備え、前記分布情報推定手段は、この画像転送手段から転送されてきた前記３次元画像データから前記被験者の心臓のモデルを設定する心臓モデル設定手段を備えたことを特徴とする心臓内電気現象の診断装置。
【請求項１２】
請求項３または１０記載の心臓内電気現象の診断装置において、
前記心臓モデル設定手段は、予め与えられた既定形状の心室モデルを移動、回転、および／または変形させ、非線形最適化法により前記３次元画像データによる心室表面をフィットさせることで、前記被験者の心室モデルの心室の位置、向き、形状を設定する手段であることを特徴とする心臓内電気現象の診断装置。
【請求項１３】
請求項１２記載の心臓内電気現象の診断装置において、
前記心臓モデル設定手段は、前記既定形状の心室モデルを対話的に移動、回転、および／または変形させる手段を備えたことを特徴とする心臓内電気現象の診断装置。
【請求項１４】
請求項１２または１３記載の心臓内電気現象の診断装置において、
前記心臓モデル設定手段は、前記３次元画像データと前記被験者の前記心室モデルとを重ねて表示する手段を備えたことを特徴とする心臓内電気現象の診断装置。
【請求項１５】
請求項３または１１記載の心臓内電気現象の診断装置において、
前記心臓モデル設定手段は、前記心臓モデルとしての心室形状を複数種類の既定形状の画像を構成するユニットを配列して表現する手段を備えたことを特徴とする心臓内電気現象の診断装置。
【請求項１６】
請求項１５記載の心臓内電気現象の診断装置において、
前記既定形状は６面体であることを特徴とする心臓内電気現象の診断装置。
【請求項１７】
請求項１記載の心臓内電気現象の診断装置において、
前記分布情報推定手段は、１種類または複数種類の医用モダリティから得られた１種類または複数種類の３次元画像データから生体モデルを作成する生体モデル作成手段を備え、この生体モデル作成手段は、前記３次元画像データから３次元の画像領域内の生体の各領域の組織種別を分別する組織種別分別手段を有することを特徴とする心臓内電気現象の診断装置。
【請求項１８】
請求項１７記載の心臓内電気現象の診断装置において、
前記生体モデル作成手段は、前記組織種別分別手段の分別結果に基づき前記領域毎の導電率を設定する導電率設定手段を有することを特徴とした心臓内電気現象の診断装置。
【請求項１９】
請求項１８記載の心臓内電気現象の診断装置において、
前記生体モデル作成手段は、前記３次元の画像領域内の生体を多数のボクセルにより構成するボクセル構成手段を備え、前記分布情報推定手段は、差分法を用いて前記多数のボクセルの電位波形および磁場波形の少なくとも一方を計算する計算手段を有することを特徴とした心臓内電気現象の診断装置。
【請求項２０】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記興奮到達時刻の分布情報を等時線、色の濃淡、あるいは色相の相違、または、それらの要素の少なくとも２つの組み合わせにより表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２１】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記心筋電気活動量の分布情報を等高線、色の濃淡、あるいは色相の相違、または、それらの要素の少なくとも２つの組み合わせにより表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２２】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記興奮到達時刻の分布情報を等時線で表示し、かつ、前記心筋電気活動量の分布情報を色の濃淡あるいは色相の相違で表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２３】
請求項２２記載の心臓内電気現象の診断装置において、
前記表示手段は、前記興奮到達時刻の分布情報と前記心筋電気活動量の分布情報とを重畳して表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２４】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記心臓の心臓表面の展開図上に、前記興奮到達時刻の分布情報、前記心筋電気活動量の分布情報、またはその両方の分布情報を表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２５】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記心臓の１つまたは複数の断面図上に、前記興奮到達時刻の分布情報、前記心筋電気活動量の分布情報、またはその両方の分布情報を表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２６】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記心臓を任意の方向から透視したときの透視投影図上に、前記興奮到達時刻の分布情報、前記心筋電気活動量の分布情報、またはその両方の分布情報を表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２７】
請求項１記載の心臓内電気現象の診断装置において、
前記表示手段は、前記心臓の任意方向の投影図と心臓表面の展開図との一方または両方に、前記興奮到達時刻の分布情報、前記心筋電気活動量の分布情報、またはその両方の分布情報を表示する手段であることを特徴とした心臓内電気現象の診断装置。
【請求項２８】
請求項１記載の心臓内電気現象の診断装置において、
前記計測手段は、前記磁場を計測するマルチチャンネル型ＳＱＵＩＤ磁束計および前記電位を計測する心電計装置のいずれかまたは両方であることを特徴とした心臓内電気現象の診断装置。
【請求項２９】
請求項１乃至２８のいずれか１項に記載の心臓内電気現象の診断装置において、
前記心筋電気活動量は、前記心筋の活動電位の立ち上がりの時相における活動電位の振幅値、前記心筋の導電率の分布、それら活動電位の振幅値および導電率の分布の積である電流ダイポール密度の絶対値、またはそれらの量に関連する量であることを特徴とした心臓内電気現象の診断装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、被検者の心臓内の電気現象を診断する装置に係り、とくに、心臓から発生する電位または磁界を計測して心臓内の電気現象を推定する心電図逆問題や心磁図逆問題の新規な解法を実施する診断装置に関する。心臓内電気現象は、本発明では具体的には、心筋の興奮到達時刻、心筋電気活動量としての活動電位振幅、心室内の導電率分布、電流ダイポール密度などで表される心臓興奮伝播過程として捉えられる。
【０００２】
【従来の技術】
従来技術１
この種の従来技術の１つとして、例えば、本発明者が既に提案している特開平８−２８９８７７号に記載の診断装置が知られている。この診断装置は、図１０に処理の流れを示すように、計算機内にソフトウエア的に心臓モデルを構築し、このモデルを基にダイクラスト法やウォーシャルフロイト法を用いて心臓の興奮伝播を効率良くシミュレートし、このシミュレーション結果から心電図や心磁図を導出する手法を提案している。
【０００３】
この手法を用いる場合、心電図または心磁図を計算するため、心室の形状、心室の位置、心室の向き、１つ又は複数の初期興奮位置とそれらの位置の興奮時刻、興奮伝播速度の心室内での分布、導電率の心室内での分布などの値を心室モデルにパラメータとして与える必要がある。
【０００４】
さらに、特許公開８−２８９８７７号において、本発明者は上述した興奮伝播シミュレーシヨン方法に基づく逆問題解法を提案している。具体的には、測定した心電図または心磁図と上記シミュレーション方法で計算した心電図や心磁図との差がなるべく小さくなるように、心室モデルに与えたバラメー夕のうちのいくつか（例えば、初期興奮位置とそれらの位置での興奮時刻）のパラメータを修正し、その差が一番小さくなったときの値（例えば、初期興奮位置とそれらの位置での興奮時刻）を推定結果として出力する。この際、選択した以外のパラメータ（例えば、初期興奮位置および興奮時刻以外のパラメータであり、興奮伝播速度分布や活動電位振幅分布）には適当な値を与え、推定のときには、これらのパラメータは固定値として扱うようにしている。
【０００５】
従来技術２
また別の従来技術として、体表面の多チャンネルの心電図のＱＲＳ波形から心室表面での興奮到達時刻分布を推定する方法が、文献「H.Roozen et al.:Compu-ting the activation sequence at the ventricular heart surface from body surface potentials, Medical & Biological Engineering & Computing, May , 1987」で報告されている。具体的には、この論文には、心筋内の導電率が一様かつ等方的であり、しかも活動電位振幅が一定であるという仮定の下で、心筋内の興奮到達時刻分布と各チャンネルに誘導される心電図ＱＲＳ波形との関係を記述する体積積分方程式（積分領域は心室筋内）を境界積分方程式に変換し、この方程式を解くことで、心電図のＱＲＳ波形の計測結果から心筋境表面の興奮到達時刻の分布を解析するという手法が示されている。
【０００６】
【発明が解決しようとする課題】
しかしながら、上述した従来技術には以下のような種々の解決しなければならない課題があった。
【０００７】
従来技術１の場合
従来技術１による心臓内電気現象の診断装置の場合、心磁図と推定するパラメータとが非線形関係にあるため、このパラメータを推定するには非常に多数の反復演算を行う必要があり、しかも、その反復演算毎に興奮伝播シミュレーシヨンを行う必要がある。興奮伝播シミユレーシヨンには多くの計算時間が必要であるため、結果として、１回の導出に要する時間が膨大になってしまうという第１の問題があった。
【０００８】
また、従来技術１による心臓内電気現象の推定法については、初期興奮部位と初期興奮時刻を変数として推定し、興奮伝播時刻や活動電位振幅分布は固定値としたときの解析結果が電気学会論文誌Ａ，Ｖｏｌ．１１６−Α，Ｎｏ．８，１９９６，６９８−７０４ページに報告されている。
【０００９】
この論文誌報告のものは、興奮伝播時刻や活動電位振幅分布を固定値として扱っているので、実際に沿わないことがある。例えば、従来技術１による心臓内電気現象の推定法を心筋梗塞などの患者に実際に適用することを想定すると、梗塞部位では興奮伝導の遅延があり、しかも活動電位の振幅は低下していると考えられる。また、初期興奮部位や興奮時刻は個人差が大きいから、興奮伝播速度分布、活動電位振幅分布、導電率分布、初期興奮部位、および初期興奮時刻を何らかの固定値で扱うということは無理がある。すなわち、これらのパラメータを全て変数として推定する必要がある。
【００１０】
しかし、その一方で、これら全てのパラメータを変数とすると、推定パラメータの数が非常に多くなり、パラメータの探索空間の大きさが膨大化し、安定に推定することはきわめて困難になる。従って、従来技術１を、心筋梗塞などを持つ患者に実際に適用し、その異常を診断することは極めて困難であるという第２の問題があった。
【００１１】
一方、心筋梗塞などの基礎疾患を持たない、不整脈などの患者に、従来技術１による心臓内電気現象の推定法を適用する場合、興奮伝播速度分布や活動電位振幅分布、導電率分布は既定値として与えておき、電気学会論文誌Α，Ｖｏｌ．１１６−Α，Ｎｏ．８，１９９６，６９８−７０４ページ記載のように初期興奮部位と初期興奮時刻のみを推定することが考えられる。
【００１２】
この場合も、望ましくは、正確に解析を行うためには心筋中の興奮伝播速度、活動電位振幅、導電率の心臓内の分布を詳細に与える必要がある。しかし、現状ではそのような詳細な情報は知られていないため、大まかな情報を与えざるを得ず、興奮伝播を高精度にシミュレートすることができない。結果として、推定結果が事前に与えた興奮伝播速度、活動電位振幅、導電率の曖昧な値に大きく影響され、解析結果が不安定になるという第３の問題があった。
【００１３】
また、不整脈の患者は現実には心筋梗塞や虚血などの基礎疾患を持っていることか多い。このため、心臓内電気現象の推定装置がこのように基礎疾患を持たない不整脈患者だけに適用できるものの場合、適用できる患者が限定されてしまうために汎用性が低いという第４の問題があった。
【００１４】
これに対して、従来技術２は、興奮伝播シミュレーションを行う必要ないため、解析に膨大な時間がかかるという問題は緩和されるものの、活動電位振幅分布、導電率分布が一様であると仮定しており、本方法の基礎方程式である体積積分方程式を境界積分方程式に変換する際、この仮定は必須のものであるため、心筋梗塞などの疾患を持つ患者には適用することができず、従来技術１と同様に適用患者が限定され、汎用性が低くなるという第５の問題があった。
【００１５】
さらに従来技術１および２に共通して、心室モデル胸部導電率分布モデルに関する問題がある。
【００１６】
これを詳述すると、従来技術１および２において、興奮到達時刻分布などの推定精度を向上させるためには被験者個々の心室形状に適合させた心室モデルを用いて正確に解析を行うことが必要である。従来法では、被験者のＭＲΙ画像から心臓のサイズと位置を手作業で測定し、代表的な心室の形状を持った心室モデルを測定値に合わせて拡大、縮小して設置していた。
【００１７】
つまり、心臓のサイズを自動的に検出することは難しいことから、被験者個々の心室形状をＭＲＩやＸ−ＣＴの３次元画像から心臓の形状を完全に自動的に抽出することは困難であり、多数の患者全てに対し日常的に患者個々に適合したモデルを作成するのは実質的に不可能である。
【００１８】
また心磁図逆問題では一般に、磁場計算の際も被験者個々の胸部形状および導電率分布に適合させた胸部導電率分布モデルを用いる必要があることが指摘されている。しかしながら、従来は、個々の被験者の導電率分布を考慮するために境界要素法や有限要素法を用いているため、心室や胸部の組織境界を３角形メッシュで構成、または胸部領域を４面体などのメッシュで構成する必要がある。その際、組織境界をメッシュに適合させる必要があるが、そのような処理の完全自動化は困難であるため、被験者個々の心室形状や胸部の導電率分布をＭＲＩやＸ−ＣＴの３次元画像から完全に自動的に抽出することができず、結果として、多数の患者全てに対し日常的に胸部導電率分布モデルを作成することは実質的に不可能である。
【００１９】
これらの状況から、患者に適合した心室形状や胸部導電率分布を考慮した推定法が実施されておらず、興奮到達時刻分布を推定精度が低く、実際の臨床には供し難いという第６の問題があった。
【００２０】
本発明の１つの目的は、上述した従来技術が有する種々の未解決の問題に鑑みてなされたもので、心臓内の興奮伝播速度の心臓内の分布に関する詳細な情報を事前に必要とせずに、心臓内の興奮伝播過程（心筋の興奮到達時刻、活動電位振幅など）を従来よりも高速に、安定に、且つ精度良く推定できるようにする、ことである。
【００２１】
本発明の別の目的は、心臓内の興奮伝播速度の心臓内の分布に関する詳細な情報を事前に必要とせずに、心臓内の興奮伝播過程を従来よりも高速に、安定に、且つ精度良く推定するとともに、心筋梗塞などの異常を持つ患者に対しても適用可能にし、汎用性を持たせる、ことである。
【００２２】
本発明のさらに別の目的は、患者個々の心室形状および胸部導電率分布を患者毎に、容易に、かつ自動的に抽出できるようにし、また、多くの患者を簡単に診断できるようにして実用性を高める、ことである。
【００２３】
本発明のさらに別の目的は、心臓内の興奮伝播速度の心臓内の分布に関する詳細な情報を事前に必要とせずに、心臓内の興奮伝播過程を従来よりも高速に、安定に、且つ精度良く推定し、また、心筋梗塞などの異常を持つ患者に対しても適用可能にし、汎用性を持たせ、さらに、患者個々の心室形状および胸部導電率分布を患者毎に、容易に、かつ自動的に抽出できるようにし、また、多くの患者を簡単に診断できるようにして実用性を高める、ことである。
【００２４】
【課題を解決するための手段】
ここで、心筋電気活動量は、具体的には、心筋の活動電位の立ち上がりの時相における活動電位の振幅値、心筋の導電率の分布、それら両者の積である電流ダイポール密度の絶対値、またはそれらの量に関連する量であり、活動電位の瞬時値のような各時刻での値ではなく、心筋各部位における各心拍の電気活動の大きさを表す量である。
【００２５】
上記種々の目的を達成するため、本発明に基づく心臓内電気現象の診断装置は、被験者の体表の測定点から電位および磁場の少なくとも一方を計測する計測手段と、前記計測手段により測定された前記被験者の電位および磁場の少なくとも一方に基づいて、当該被験者の心臓内の電気現象を表す情報としての、当該心臓内の興奮到達時刻の分布情報と当該心臓の心筋の電気活動の大きさを表す心筋電気活動量の分布情報とを同時に推定する分布情報推定手段と、前記興奮到達時刻の分布情報と前記心筋電気活動量の分布情報とを表示する表示手段と、を備えたことを特徴とする。
【００２７】
好適には、前記分布情報推定手段は、前記心臓内での前記興奮到達時刻の分布と前記心筋電気活動量の分布とを有限個の複数のパラメータで表現する表現手段と、このパラメータにより表現された前記興奮到達時刻の分布と前記心筋電気活動量の分布とに基づいて前記電位および磁場の少なくとも一方の各測定点での時間波形を計算する波形計算手段と、前記測定電位および磁場の少なくとも一方の波形と前記計算電位および磁場の少なくとも一方の時間波形とが一致するように非線型最適化法を用いて前記複数のパラメータを推定するパラメータ推定手段と、を備える。この場合、例えば、分布情報推定手段は、前記被験者の体内の３次元画像データから前記心臓のモデルを設定する心臓モデル設定手段を備える。また、前記表現手段は、前記心臓内での前記興奮到達時刻の分布と前記心筋電気活動量の分布とを前記心臓モデルの代表点の数よりも少ない数で表現する手段であり、前記分布情報推定手段は、この少ない数のパラメータに基づき前記興奮到達時刻の分布と前記心筋電気活動量の分布とを算出する算出手段を備える、ことも好適である。
【００２８】
さらに好適には、前記表現手段は、前記心臓モデルの各代表点での前記興奮到達時刻、前記心臓内での前記心筋電気活動量で表された前記興奮到達時刻の分布、および前記心筋電気活動量の分布を、前記心臓モデルの代表点よりも少ない数のパラメータの線形関数で表現する手段であり、前記分布情報推定手段は、この線形関数で表現された少ない数のパラメータに基づき前記興奮到達時刻の分布と前記心筋電気活動量の分布とを算出する算出手段を備える。このとき、例えば、前記線形関数はフーリエ変換または逆フーリエ変換による関数である。さらに、例えば、この線形関数は行列演算による関数である。また一態様として、この行列演算の行列は、測定した又は仮想した興奮到達時刻の共分散行列から求めた行列である。
【００２９】
さらに好適には、前記心臓モデル設定手段は、前記心臓モデルとして心室の形状を多数の多面体の組み合わせにより設定する手段である。この多面体は、例えば、６面体、５面体、もしくは４面体、または、それらの多面体の内の少なくとも２種類の組み合わせであることを特徴とする。
【００３０】
さらに好適な一態様として、医用モダリティから得た前記被検者の体内の３次元画像データを入力し転送する画像転送手段を備え、前記分布情報推定手段は、この画像転送手段から転送されてきた前記３次元画像データから前記被験者の心臓のモデルを設定する心臓モデル設定手段を備える。
【００３１】
さらに別の好適な一態様として、前記心臓モデル設定手段は、予め与えられた既定形状の心室モデルを移動、回転、および／または変形させ、非線形最適化法により前記３次元画像データによる心室表面をフィットさせることで、前記被験者の心室モデルの心室の位置、向き、形状を設定する手段であることを特徴とする。この心臓モデル設定手段は、前記既定形状の心室モデルを対話的に移動、回転、および／または変形させる手段を備えることが望ましい。
【００３２】
また、心臓モデル設定手段は、前記３次元画像データと前記被験者の前記心室モデルとを重ねて表示する手段を備える、ことも望ましい。
【００３３】
さらに別の好適な一態様として、前記心臓モデル設定手段は、前記心臓モデルとしての心室形状を複数種類の既定形状の画像を構成するユニットを配列して表現する手段を備える。この既定形状は例えば６面体である。
【００３４】
さらに別の好適な一態様として、前記分布情報推定手段は、１種類または複数種類の医用モダリティから得られた１種類または複数種類の３次元画像データから生体モデルを作成する生体モデル作成手段を備え、この生体モデル作成手段は、前記３次元画像データから３次元の画像領域内の生体の各領域の組織種別を分別する組織種別分別手段を有することを特徴とする。この場合、前記生体モデル作成手段は、前記組織種別分別手段の分別結果に基づき前記領域毎の導電率を設定する導電率設定手段を有することもできる。また、この場合、前記生体モデル作成手段は、前記３次元の画像領域内の生体を多数のボクセルにより構成するボクセル構成手段を備え、前記分布情報推定手段は、差分法を用いて前記多数のボクセルの電位波形および磁場波形の少なくとも一方を計算する計算手段を有することもできる。
【００３６】
これにより、本発明では例えば、心臓内の興奮到達時刻分布と心筋電気活動量分布が同時に推定される。また好適には、興奮到達時刻分布や心筋電気活動量分布を生成する比較的少数のパラメータが推定対象となる。心筋電気活動量分布の代わりに、導電率分布や、心筋電気活動量と導電率の積（電流ダイポール密度）の分布を推定してもよい。さらに好適には、心室形状、胸部形状が多数のボクセルにより構成され、差分法を用いて電位計算・磁場計算が実行される。
【００３７】
本発明のそのほかの特徴および効果は、以下に説明する発明の実施の形態および添付の図面により明らかになる。
【００３８】
【発明の実施の形態】
以下、本発明の実施の形態を添付図面を参照して説明する。
【００３９】
第１の実施の形態
本発明の第１の実施形態を図１〜図４、図９に基づき説明する。
【００４０】
図１に、この実施形態に係る心臓内電気現象の診断装置の構成の概要を示す。この診断装置は、ＳＱＵＩＤ（超伝導量子干渉素子）磁束計１１、画像収集装置１２、画像転送装置１３、およびコンピュータ装置１４を備える。
【００４１】
ＳＱＵＩＤ磁束計１１は、ＳＱＵＩＤセンサ１１ａ、駆動回路１１ｂ、および信号収集回路１１ｃを備え、高感度な磁気センサであるマルチチャネルＳＱＵＩＤ磁束計に構成されている。このＳＱＵＩＤ磁束計１１は、ＳＱＵＩＤセンサを用いて胸面上の数点（測定点）での心臓磁場波形を計測し、この波形信号を増幅／フィルタ処理し、Α／Ｄ変換によりデジタルデータに変換する。このデジタル量の磁場波形データはコンピュータ装置１４に供給される。
【００４２】
画像収集装置１２は例えばＭＲＩ（磁気共鳴イメージング）装置、Ｘ線ＣＴ  （Ｘ−ＣＴ）、超音波診断装置などの３次元画像データを取得するモダリティの１つまたは複数で構成され、患者の体内の画像データ（生画像データ）を収集し、画像転送装置１３に供給する。画像転送装置１３は例えばＣＰＵおよびメモリを備えて構成され、送られてきた画像データを保管するとともに、コンピュータ装置１４からの指令に応じて画像データをコンピュータ装置１４に供給する。画像転送装置１３は好ましくは、ネットワーク経由または記憶媒体経由で画像データをコンピュータ装置１４に供給するように構成され、また必要に応じて画像フォーマットの変換や３次元画像データへの再構成も行うように構成されている。
【００４３】
コンピュータ装置１４は、心臓内の興奮伝播過程を推定するための解析および演算を行うもので、ＣＰＵを要部とする演算処理器１４ａ、メモリ１４ｂ、入力器１４ｃ、および表示器１４ｄを備える。メモリ１４ｂには、かかる解析／演算用のプログラムが予め内蔵されており、演算処理器１４ａがこのプログラムを実行する。この実行の一例を図２のフローチャートに示す。
【００４４】
このフローチャートを説明する。まず、演算処理器１４ａは、ＳＱＵＩＤ磁束計１１から供給された磁場波形データを入力し（ステップＳ１）、この波形データの信号処理を実行する（ステップＳ２）。この信号処理では、雑音除去などのためにデジタルフィルタや心拍毎の加算平均処理などが行われる。そして、心磁図のＱＲＳ群に相当する時間が検出され、各チャンネルのＱＲＳ群に相当する磁場波形が抽出される。時刻の基準点の取り方は任意で良いが、例えば、ＱＲＳ群の開始点付近を時刻＝０としたり、最大振幅の時刻を０としたりすればよい。
【００４５】
この後、演算処理器１４ａは、画像転送装置１３から複数枚の断層画像データ（３次元データ）を入力し（ステップＳ３）、この画像データを基に心臓形状の測定および設定を行う（ステップＳ４）。これにより、例えばＭＲＩ装置の収集に拠る３次元画像データから心室筋の領域が検出され、心室形状モデルが構築される。
【００４６】
図３に、構築する心室モデルの基礎になるモデルを示す。この例では心室の形状が多数の６面体の要素（ユニット）の組み合わせで表現されている。この心室モデルは、長軸方向の長さ、短軸直交方向それぞれの左心室外径、左心室壁厚、右心室壁厚、中隔と右心室壁の間の隙間の長さなどのパラメータを、被験者個々に合わせて入力器１４ｃを介して操作し、変形させることが出来る。長軸方向の長さ以外のパラメータは、心尖部、心基部、およびその間の任意のレベル毎に指定することができる。
【００４７】
この心室モデルを被験者個々の心室形状にフィットさせるフィット処理は、ここでは例えば手作業による対話的手法で、あるいは非線型最適化手法を用いて自動的に行われる。いずれの手法においても、被験者のＭＲＩ、Ｘ−ＣＴ、超音波診断装置などから収集された３次元画像の心室の表面と心室モデルの表面とがなるべく一致するように、心室モデルを変形する各パラメータと心室モデルの被験者体内の位置、向きを表現するパラメータとが変更される。
【００４８】
手作業による対話的手法を用いて心室モデルの形状を３次元画像の心室形状にフィットさせる場合、３次元画像の任意方向の投影図または任意位置での断面図と心室モデルの形状とを同時に表示器１４ｄに表示し、操作者が表示画面を確認しながらパラメータを変更する態様が好ましい。さらに、その変更したパラメータにしたがって心室モデルの位置、向き、形状が変化し、その変化の様子が表示画面上に連動して表示器１４ｄに表示させる構成も望ましい。
【００４９】
非線形最適化手法を用いて自動的にフィットさせる場合も、フィッティング中の心室位置、向き、形状や、フィットさせた最終的な心室位置、向き、形状が３次元画像の心室形状と共に表示器１４ｄに表示させる構成が望ましい。さらに、必要に応じて手作業で心室の形状を変化させ、その結果を初期値としてさらにフィッティング作業を自動的処理で再開できるように構成することが望ましい。
【００５０】
コンピュータ装置１４の演算処理器１４ａは、この後、設定した心臓形状のデータと測定（入力）した被験者のＱＲＳ期間の心磁図のデータとに基づいて、心臓の心室内の興奮伝播過程の推定値としての興奮到達時刻分布および活動電位振幅分布の解析を行う（ステップＳ５〜Ｓ１２）。
【００５１】
興奮到達時刻分布および活動電位振幅分布を解析するアルゴリズムとしては、各種の非線形最適化アルゴリズムが使用できる。これらのアルゴリズムとしては、準ニュートン法、共役勾配法、コンジュゲートディレクション法、シミュレーテッドアニーリング法、遺伝的アルゴリズム法など、多様なものがある。これらの内、典型的なアルゴリズムには、興奮到達時刻分布・活動電位振幅分布の初期値設定処理（図２、ステップＳ５）、心臓内電流源分布の計算処理（同図、ステップＳ６）、磁場計算処理（図２、ステップＳ７）、誤差算出処理（図２、ステップＳ１０）、および、興奮到達時刻分布・活動電位振幅分布の変更処理（図２、ステップＳ９）が含まれる。これらの処理ステップのそれぞれをどのような順序で実行するは、それぞれの非線形最適化アルゴリズムにおいて周知であり、採用するアルゴリズムにしたがってその順序が決められる。図２には、準ニュートン法など、多くのアルゴリズムに適用される順序を模式的に示している。
【００５２】
以下に、上記処理ステップＳ５〜Ｓ１０の処理内容を説明する。
【００５３】
ステップＳ５に係る興奮到達時刻分布・活動電位振幅分布の初期値設定処理では、心内膜・外膜興奮到達時刻の初期値設定および活動電位振幅の初期値設定が行われる。
【００５４】
心内膜・外膜興奮到達時刻の初期値設定の場合、図３の心室モデルの心内膜側、あるいは心外膜側の全ての頂点に対して興奮到達時刻の初期値を設定する。これは予め定められた範囲内、例えばＱＲＳ開始時刻を基準として０から８０msの間でランダムに決定してもよいし、正常心臓の典型的な興奮開始時刻を設定してもよい。
【００５５】
活動電位振幅も同様に心室モデルの心内膜側、あるいは心外膜側の全ての頂点に対して活動電位振幅の初期値を設定する。その範囲は、ランダムでもよく、あるいは、心筋電気活動量として具体的には活動電位振幅の量を採用する場合には、正常の活動電位振幅である約−９０mVを全ての頂点に与えてもよい。さらに、図３の心室モデルにおいて心外膜と心内膜以外の点にも代表点を設け、心内膜、心外膜の頂点とそれ以外の点を含む代表点全ての興奮到達時刻と活動電位振幅の値を初期値として設定してもよい。
【００５６】
ステップＳ６に係る心臓内電流原分布の計算処理では、心臓モデルの各６面体の内部の興奮到達時刻や活動電位振幅を補間演算し、それら補間によって求めた興奮到達時刻と活動電位振幅の値を用いてＱＲＳ期間中のいくつかの時刻の興奮波面の形状と興奮波面上の電流源分布が計算される。補間方法としては有限要素法で頻繁に用いられる例えば次式で表される補間多項式を用いればよい。なお、この補間方法は興奮到達時刻分布を例に説明するが、活動電位振幅に対しても同様の方法がとられる。
【００５７】
【数１】
0003660781.tif 000002
【外１】
0003660781.tif 000003
【００５８】
【数２】
0003660781.tif 000004
【外２】
0003660781.tif 000005
【００５９】
さらにステップＳ７における磁場計算処理では、計算された多数の電流ダイポールの位置・向き・大きさに基づき、胸面上に設置された各ＳＱＵＩＤセンサで計測されるべき磁場を、ＱＲＳ期間中の各時刻毎に計算する。この計算は、胴体の導電率分布を無限一様媒質または半無限平板と仮定してビオサバールの法則に基づいて行ってもよいし、同心球状導電率球と仮定してサーバスの式を用いて行ってもよい。あるいは、境界要素法や有限要素法を用いて心筋、血液、肺、脂肪層、皮膚などの実形状を考慮した磁場計算を行うこともできる。
【００６０】
このステップＳ６およびＳ７の処理は、必要回数だけ反復される（ステップＳ８）。
【００６１】
さらにステップＳ９の興奮到達時刻分布・活動電位振幅分布の変更処理では、心臓モデルに設定されている興奮到達時刻分布と活動電位振幅分布が最終的に、ステップＳ１，Ｓ２で計測した心磁図とステップＳ７で計算した心磁図の差が小さくなるように、ステップＳ８で反復計算した値を基に変更される。具体的にどのように変更するかは各種非線形最適化アルゴリズムにより異なる。
【００６２】
さらにステップＳ１０の誤差算出処理では、ステップＳ７で計算したＱＲＳ期間の心磁図と、ステップＳ１，Ｓ２の処理を介して得たＱＲＳ期間の心磁図とがどれ位異なるのかを示す誤差ｅが計算される。例えば、チャネルｉが測定したｊ番目の時刻の心磁図をＢij、同様に、計算により求めた心磁図をＣijとしたとき、
【数３】
0003660781.tif 000006
に式により誤差ｅを求める。ｂはＢijの２乗平均平方根、ｃはＣijの２乗平均平方根である。
【００６３】
そして、この興奮到達時刻分布・活動電位振幅分布の解析処理においては、上述した各処理ステップＳ６〜Ｓ１０が、採用した非線形最適化アルゴリズムで要請されるように反復的に実行される（ステップＳ１１）。これにより、誤差ｅがなるべく小さくなるような興奮到達時刻と活動電位振幅の心室内の分布が求められる。
【００６４】
このようにして推定した興奮到達時刻分布および活動電位振幅分布のデータは、その後、表示器１４ｄに表示される（ステップＳ１２）。
【００６５】
この表示結果の例を図４に示す。この例は、心内膜・外膜の展開図の上に興奮到達時刻の等時線図を表示したものである。等時線図は同時刻に興奮を開始する点を線で結ぶ表示方法である。また、図のように、色の濃淡あるいは色相の違いにより活動電位振幅の分布を重ねて表示している。
【００６６】
この例は心室を長軸方向に沿って切り込みを入れ、切り開いた展開図で示しているが、心外膜、左心室内膜、右心室内膜をそれぞれ円状に展開した図の上に等時線図や活動電位振幅分布を表示してもよい。また、心室の断面上に興奮到達時刻分布と活動電位振幅分布を表示してもよい。また、心臓モデルの透視投影図上に重ねて表示してもよい。さらに、等時線図と活動電位振幅を重ねるものではなく、別々に表示してもよい。さらに興奮開始時刻や活動電位振幅の数値を図中に示しようにしてもよい。
【００６７】
さらに、推定した興奮到達時刻分布および活動電位振幅分布のデータの表示例を図１１〜図１４に示す。これらの表示を行う処理は演算処理器１４ａにより実行される。
【００６８】
図１１は、心筋の２枚の同一長軸断面図それぞれの上に興奮到達時刻分布の等時線図と活動電位振幅分布の濃淡画像を個別に重ね、それらを同時に表示する例である。同図に示すように、左心室心尖部に活動電位振幅が低下している領域が存在することから、この部位に何らかの心筋異常が在ることを容易に判別することが可能になる。つまり、興奮到達時刻分布および活動電位振幅分布の双方を同時に表示していることから、興奮伝導異常を伴う各種の不整脈、活動電位の異常を伴う心筋虚血や心筋梗塞など、様々な心疾患の診断表示に優れている。また、従来の電流源密度分布表示に比べ、各種の心疾患を判別するのに必要な画像枚数も大幅に低減でき、読影に要する手間も低減できる。
【００６９】
図１２は、心筋の同一の長軸断面図上に興奮到達時刻分布の等時線図と活動電位振幅分布の濃淡画像とを重ねて表示する例である。これにより、図１１で説明した利点に加え、１枚の画像で済むので、画面上の表示スペースが節約できる。このため、相異なる複数枚の心室断面上の解析結果（興奮到達時刻分布および活動電位振幅分布）を１枚の表示画面や紙に並べて提示できる。一方、興奮到達時刻分布と活動電位振幅が重ねられて表示されているので、例えば、活動電位振幅の異常領域内の興奮到達時刻や、異常領域中の興奮波面形状の対応関係の読取りを容易化させることができ、診断能向上や診断能率向上に寄与可能になる。なお、図１３には図１２と同様の表示を心室の短軸断面上に行った例を示す。
【００７０】
さらに図１４（ａ）、（ｂ）は、心室形状モデルの透視投影図上に、心室表面での興奮到達時刻分布を重ねて表示する例に関する。同図（ｂ）は透視投影図の視線方向の指定方法を説明する図である。同図（ｂ）の画面左下および画面左上の断面は前述したと同様に表示するものであり（図１１〜１３参照）、画面右は任意方向からの透視投影図あるいは平行投影図上に興奮到達時刻分布または活動電位振幅分布、またはその両方を表示するものである。画面左下の矢印は視線方向を示す。この矢印は操作者が入力器１４ｃから断面上で対話的に動かすことができる。この矢印をさらに画面左上の断面上に移動させることもできる。操作者は、心室の回転と矢印の移動により、視線方向を任意方向に設定できる。同図  （ａ）は同図（ｂ）の画面右側で表示される一例である、投影図上に興奮到達時刻分布を重ねて表示する例を具体的に示している。このように表示することで、上述した各種の利点に加え、心筋虚血を観察する場合、心外膜上の虚血領域が容易に見つけ易くなるという利点もある。
【００７１】
なお、第１の実施形態の上述した処理において、ステップＳ１，Ｓ２の処理が信号処理手段を、ステップＳ３，Ｓ４の処理が心臓形状設定手段を、ステップＳ５の処理が初期値設定手段を、ステップＳ６，Ｓ８，Ｓ１１の処理が電流原分布計算手段を、ステップＳ７，Ｓ８，Ｓ１１の処理が磁場計算手段を、ステップＳ９，Ｓ１１の処理が変更手段を、ステップＳ１０，Ｓ１１の処理が誤差算出手段を、ステップＳ１２の処理が表示指令手段にそれぞれ対応している。この各手段の機能的な繋がり具合を図９に示す。
【００７２】
以上のように、本実施形態およびその変形例によれば、活動電位振幅分布と興奮到達時刻分布を同時に推定できる。心室の電気的活動は、心拍が一定のレート以下であるどのような正常及び異常の心室に対しても、活動電位振幅分布と興奮到達時刻分布によって近似的に表現することができる。例えば、心筋虚血あるいは心筋梗塞は、虚血および梗塞部位の興奮伝播が正常部位よりも遅れ、活動電位振幅が小さいというように表現することが可能である。さらに、心室頻拍などの不整脈は、心室細動時のように心拍数が極端に小さくならなければ正常と異なる興奮到達時刻分布として表現することができる。したがって、本実施形態による心臓内電気現象の診断装置は、心筋梗塞や不整脈などの様々な疾患を持つ患者に対して適用でき、臨床の場での汎用性が非常に高いという利点がある。
【００７３】
また、本実施形態の診断装置の場合、従来法のように計算時間がかかる興奮伝播シミュレーシヨンの過程が反復中に含まれていないという特徴がある。上述した診断装置の解析の反復処理中に含まれる演算の内、計算量が最も多いのは磁場計算処理であるが、これは四則演算の繰り返しの多い線形演算であるので、べクトル化、並列化、バイプライン処理化、および／またはキャッシュメモリなどの採用により現在使用されている計算機で高速に演算することができる。したがって、従来法よりも解析に要する時間を大幅に短縮することができる。
【００７４】
さらに、従来法のように、推定演算中に興奮伝播速度の心筋内の詳細な情報を事前に必要とする構成ではないので、不正確な興奮伝播速度分布を推定に使用することに因る推定結果の不安定性を回避することができる。
【００７５】
第２の実施の形態
本発明の第２の実施形態を図５に基づき説明する。なお、この実施形態に係る心臓内電気現象の診断装置の構成は上述したものと同一または同等である。
【００７６】
この実施形態は、心臓形状の表現方法の他の例を説明するものである。
【００７７】
図５に本実施形態で用いる心室モデルを示す。第１の実施形態の図３との違いは、図３のものが６面体を用いて心臓形状を表現するのに対して、図５のものは５面体（三角柱状の形状）を用いて心室形状を表現する点にある。
【００７８】
この比較を詳述する。第１の実施形態では心臓を６面体で分割しているので、心尖部の６面体の大きさよりも心基部の６面体の大きさの方が小さくなってしまう。このため、計算量が同一の場合、均一な大きさの６面体を用いる場合に比べて、磁場計算の精度が劣る。しかし、この第２の実施形態はその点も改善できるものであり、心臓を５面体で分割することにより、心臓全体を均一な大きさに分割できる。したがって、同一の計算量の場合でも、高い磁場計算の精度を維持できる。
【００７９】
具体的には、前述した図２のステップＳ３、Ｓ４の処理において、この５面体による心室モデルが使用され、この心室モデルが被検者個々の心臓に形状に合わせたモデルが設定される。この設定した５面体による心室モデルを使って、図２のステップＳ６の処理にて、第１の実施形態のときと同様に電流原分布が計算される。
【００８０】
したがって、前述したように、この５面体による心室モデルは６面体のものよりも、同じ計算量で比較して、磁場計算の精度、すなわち推定解析結果の精度をさらに向上させることができるという効果がある。
【００８１】
なお、前述した図３では６面体により心室形状を表現し、上記図５では５面体により心室形状を表現したが、４面体など他の多面体により心室形状を設定するように、図２のステップＳ４に係るソフトウエアを構成してもよい。さらに、それら数種類の多面体を適宜に組み合わせて心室形状を設定するように構成してもよい。
【００８２】
第３の実施の形態
本発明の第３の実施形態を図６に基づき説明する。
【００８３】
本実施形態は、興奮到達時刻や活動電位振幅の分布を、心内膜・外膜上の全ての頂点での値を用いて表現するのではなく、もっと少数のパラメータで表現することを目的とする。なお、以下の説明では、興奮到達時刻を例示するが、活動電位振幅に対しても同様の方法で処理することができる。
【００８４】
【外３】
0003660781.tif 000007
【数４】
0003660781.tif 000008
と表される。
【００８５】
【外４】
0003660781.tif 000009
【数５】
0003660781.tif 000010
のように表される。Ｐ1 ，Ｐ2 はＰの部分行列で、Ｐ1 はＮ×Ｍ，Ｐ2 はＮ×  （Ｎ−Ｍ）である。また、ＭはＮより小さい。この行列Ｐ1 を用いて
【数６】
0003660781.tif 000011
【外５】
0003660781.tif 000012
【数７】
0003660781.tif 000013
によりその近似値を得ることはできる。
【００８６】
本実施形態に係る心臓内電気現象の診断装置は上述した原理に基づくものであり、そのハード的構成は前述した第１の実施形態のものと同一または同等に形成されている。この構成において、コンピュータ装置１４の演算処理器１４ａは図６に示す一連の処理を行うようになっている。
【００８７】
この図６の処理は、前述した図２のものと比較して、ステップＳ５とＳ６との間に新しくステップＳ５ａ，Ｓ５ｂが挿入され、ステップＳ９がステップＳ９ａに置換され、ステップＳ１１で反復判断がなされたときの戻るフロー上の位置はステッププＳ５ａとＳ５ｂの間になっていることである。
【００８８】
【外６】
0003660781.tif 000014
ステップＳ５ａでは、演算処理器１４ａによって、変数ｘ，ｙの初期値が設定される。具体的には、前述したステップＳ５で設定した興奮到達時刻と活動電位振幅の初期値を元に式（７）を用いて変数ｘ，ｙの初期値が設定される。
【００８９】
次いで、ステップＳ５ｂでは、興奮到達時刻および活動電位振幅が計算される。すなわち、式（６）を用いて変数ｘ，ｙから心室内の各代表点での興奮到達時刻と活動電位振幅が算出される。
【００９０】
次いで、ステップＳ６，Ｓ７の処理を介して（ステップＳ５ｂ，Ｓ６，Ｓ７の反復処理を含む）ステップＳ９ａの変数ｘ，ｙの変更処理が実施される。この変更処理では、第１の実施形態における心興奮到達時刻分布・活動電位振幅分布の変更処理と同様に変数ｘ，ｙが変更される。
【００９１】
この後、さらにステップＳ１０、Ｓ１１の処理が同様に実行される。
【００９２】
このため、本実施形態の解析法によれば、第１、第２の実施形態のときより少ない数の変数を推定すればよいので、推定が安定し、しかも，より短時間に推定することができる。
【００９３】
しかし、この第３の実施形態の解析法では、Ｍにあまり小さな値を用いると、心内膜・外膜興奮到達時刻分布の自由度か小さくなるため、推定した心内膜・外膜興奮到達時刻分布の空間分解能が低下することがある。また、推定した心内膜・外膜興奮到達時刻に基づいて計算された心磁図か実測の心磁図と一致し難くなり、推定誤差が大きくなることがある。
【００９４】
これを避けるためには、非線型最適化アルゴリズム実行中の反復初期においては、Ｍに小さい値を用い、収束か進むにつれてＭを大きくしてＮに近づけるように構成することが望ましい。このとき、Ｍを増やすことに伴って新たに増える変数には初期値０を与えておけばよい。このように構成することで、Ｍを小さくとったために興奮到達時刻分布の推定値の推定誤差が大きくなったり、空間分解能が低下したりする問題を回避ことができる。
【００９５】
また、興奮到達時刻分布のパラメータ数Ｍと活動電位振幅分布のパラメータ数は互いに異なる値を用いてもよい。これにより、例えば、興奮到達時刻分布の空間分解能が活動電位振幅分布の空間分解能より低くても済む場合、興奮到達時刻分布パラメータ数を少なくすることができ、推定の安定性をさらに向上させることができる。
【００９６】
なお、上述した図６の処理において、初めから変数ｘ，ｙに限定した処理を行う場合、ステップＳ５の初期値設定処理は省略することもできる。
【００９７】
第４の実施の形態
本発明の第４の実施形態を図７、図８に基づき説明する。この実施形態による心臓内電気現象の診断装置は、多数のボクセルにより構成された心室モデルを使用するものである。ハード的構成は第１の実施形態のものと同一または同等に形成されている。
【００９８】
本実施形態では、心室興奮到達時刻分布を以下のように表現する。活動電位振幅分布も同様に表現可能であるので説明を省略する。
【００９９】
まず、各ボクセルのｘ，ｙ，ｚ方向の各方向の配列の順番ｉ，ｊ，ｋにより各ボクセルを表す。ｉ，ｊ，ｋ方向のボクセル数をＮｉ，Ｎｊ，Ｎｋとする。ボクセル（ｉ，ｊ，ｋ）の興奮到達時刻をｔ（ｉ，ｊ，ｋ）と表す。ｔ（ｉ，ｊ，ｋ）の３次元離散的フーリエ変換をＴ（ｌ，ｍ，ｎ）と表す。Ｔ（ｌ，ｍ，ｎ）は複素数である。ｔ（ｉ，ｊ，ｋ）は実数なのでｌ，ｍ，ｎの個数は、それぞれ、Ｎｉ，Ｎｊ，Ｎｋの半分である。
【０１００】
ここで小さいｌ，ｍ，ｎの成分以外はＴ（ｌ，ｍ，ｎ）＝０であると仮定する。これはｔ（ｉ，ｊ，ｋ）に高周波成分が含まれないと仮定することに相当する。心室に相当しないボクセルの興奮到達時刻も表現されることになるが、そのような値は無視することにする。この手法を用いると、ｌ，ｍ，ｎの小さい成分のみ、すなわち少数のパラメータにより心室内の興奮到達時刻ｔ（ｉ，ｊ，ｋ）を表現することができる。
【０１０１】
図７および図８に、コンピュータ装置１４の演算処理器１４ａによって実行される、ボクセル状の心室モデルと、本実施形態特有の興奮到達時刻分布・活動電位振幅分布の表現法を用いた興奮到達時刻分布・活動電位振幅分布の解析の概要をフローチャートで示す。
【０１０２】
演算処理器１４ａは、かかる解析を行うために、図７に示すステップＳ２１〜Ｓ３３までの処理を行う。この内、ステップＳ２１〜Ｓ２３までの処理は前記図２のステップＳ１〜Ｓ３までの処理と同様に実行される。その後、ステップＳ２４に移行すると、ボクセル状の心室モデルが構成される。この構成に関わる、より詳細な処理の流れを図８に示す。
【０１０３】
図８の処理を説明する。演算処理器１４ａがここで行う処理は、以下に詳述する位置合わせ処理（ステップＳ２４−１）、画像処理（ステップＳ２４−２）、組織種別分別処理（ステップＳ２４−３）、ボクセル再作成処理（ステップＳ２４−４）、および導電率分布設定処理（ステップＳ２４−５）である。
【０１０４】
まず、位置合わせ処理（ステップＳ２４−１）では、ＭＲＩ，Ｘ−ＣＴで測定した胸部の各種３次元画像とＳＱＵＩＤセンサの各チャネルのそれぞれの座標系の位置関係を求める。各３次元画像間の位置合わせには、各種画像上に複数のマーカを予め写し込んでおき、それぞれのマーカが物理的に同じ位置にあるとして位置関係を求めることが好適な一例である。ＳＱＵＩＤセンサと各画像間の位置合わせにも上記マーカと同じ位置から磁場を発生させ、その計測値からマーカの位置を求め、３次元画像上のマーカと同じ位置にあるとして両者の位置関係を求める。
【０１０５】
３次元画像間の他の位置合わせ法として、各画像に撮影された被験者の胸部の表皮面を一致させるように位置合わせを行ってもよい。この場合、ＳＱＵＩＤセンサとの位置合わせは、胸面上に付けた磁場発生器の位置をＳＱＵＩＤ磁束計を用いて計測し、それらの位置と３次元画像土の胸部の表皮面とを一致させるように位置合わせを行えばよい。
【０１０６】
画像処理（ステップＳ２４−２）では、この後の組織種別分別処理を行い易くするための前処理が行われる。例えば、画像の低階調化、雑音の除去、平滑化、画像サイズの変換、微分処理、エッジ検出などが行われる。
【０１０７】
なお、組織種別データーベースＴＤＢが設けられている。このデータベースはＭＲＩのＴ１，Ｔ２強調像、Ｘ−ＣＴ画像の各ボクセルの信号強度と組織種別の間の関連を記述している。この記述において、心室筋、心房筋、血管、弁、その他筋組織など同じような信号強度を示すものは、同じカテゴリー（筋肉など）として分類されている。
【０１０８】
さらに、組織種別分別処理（ステップＳ２４−３）により、信号処理を施した各種画像の信号強度から組織種別データーベースを参照して、組織種別をボクセル毎に各カテゴリに分別される。このようにして得られたカテゴリの画像を元に、心室など各組織の形状や位置などの既知の情報を用いて、各ボクセルが、例えば、心室筋、心房筋、心臓弁、心臓以外の筋組織、血液、肺、脂肪、皮膚、空気、外気に分類される。とくに、同じカテゴリに含まれる心室筋、心房筋、心臓弁、心臓以外の筋組織の間の分類には、筋組織の厚さ、外気からの距離を用いて分類が行われる。
【０１０９】
またボクセル再作成処理（ステップＳ２４−４）では、無駄な外気領域や首、下腹部を除去し、胸部だけが画像に収まるように画像サイズが変換され、なるべく少数のボクセルで胴体部分の導電率分布や心室形状を表現できるようにボクセルサイズの変換が行われる。このとき、全てのボクセルを同じサイズに処理するのではなく、所定の最小単位の大きさのボクセルの整数倍の大きさとなるようにしても処理してもよい。
【０１１０】
最後に、導電率分布設定処理（ステップＳ２４−５）では、分類した組織種別の代表的な導電率の値が記述された導電率データーベースＣＤＢを参照することにより、各ボクセルの導電率が設定される。
【０１１１】
次いで、興奮到達時刻分布・活動電位振幅分布の初期値設定処理が行われる  （ステップＳ２５）。この処理では、具体的には、正常心臓の典型的な興奮到達時刻分布の３次元離散的フーリエ変換の小さいｌ，ｍ，ｎ以外の値を０と置いたＴ（ｌ，ｍ，ｎ）を変数Ｔの初期値として設定する。あるいは、Ｔ（ｌ，ｍ，ｎ）の小さいｌ，ｍ，ｎの成分をランダムに設定し、それ以外の成分を０と置いて、これを変数Ｔの初期値として設定するようにしてもよい。なお、活動電位振幅分布に対しても同様の初期設定処理が行われる。
【０１１２】
次いで、興奮到達時刻分布・活動電位振幅分布の算出処理が行われる（ステップＳ２６）。具体的には、Ｔ（ｌ，ｍ，ｎ）が３次元離散的逆フーリエ変換され、興奮到達時刻の分布が求められる。活動電位振幅分布に対しても同様の算出処理を行われる。
【０１１３】
【外７】
0003660781.tif 000015
【数８】
0003660781.tif 000016
のように求める。これにより、
【数９】
0003660781.tif 000017
と表される。ここで、σ（ｉ，ｊ，ｋ）はボクセル（ｉ，ｊ，ｋ）での細胞内導電率、Ａｘ（ｉ，ｊ，ｋ），Αｙ（ｉ，ｊ，ｋ），Αｚ（ｉ，ｊ，ｋ）はボクセル（ｉ，ｊ，ｋ）のｉ，ｊ，ｋ方向にそれぞれ直角な面でのボクセルの断面積、Φ（ｉ，ｊ，ｋ）はボクセル（ｉ，ｊ，ｋ）での活動電位振幅である。
【０１１４】
次いで、磁場計算処理が実行される（ステップＳ２８）。この磁場計算処理は、第１の実施形態のときと同様に、ビオサバールの法則やサーバスの式などを用いて、ＳＱＵＩＤセンサで計測される磁場の大きさが計算される。なお、心臓形状測定・設定処理（ステップＳ２４）で設定した組織の導電率分布を利用し、差分法を用いて磁場を計算するようにしてもよい。差分法は格子状の代表点を用いて偏微分方程式を解く手法であり、境界要素法や有限要素法に比べて、アルゴリズムが簡単であるという特徴を有している。
【０１１５】
ステップＳ２６〜Ｓ２８の処理は必要回数だけ反復して実行される（ステップＳ２９）。
【０１１６】
次いで、興奮到達時刻分布・活動電位振幅分布の変更処理、および誤差算出処理が、第１の実施形態のときと基本的に同様に必要回数だけ繰り返して実行される（ステップＳ３０〜３２）。さらに、興奮到達時刻分布・活動電位振幅分布の解析結果の表示処理が、第１の実施形態のときと基本的に同様に実行される（ステップＳ３３）。
【０１１７】
従来法では、推定精度を向上させるためには、心臓や肺、胴体などの実形状を表現することが必要で、このために境界要素法や有限要素法を使用していた。例えば境界要素法では３角形の組み合わせで心臓や肺、体表の表面を構成する。このようなモデルを作成するためにはＭＲＩ、Ｘ−ＣΤからの３次元画像における組織の境界を抽出し、境界面を近似する三角形メッシュを作成する必要がある。しかし、この作成を全て自動的に行うのは困難であるのが実状であり、検査を行う全ての患者に対して患者個々に合わせた組織境界の３角形メッシュモデルを作成するのは実際上不可能であった。有限要素法の場合、３次元領域を４面体や６面体のメッシュで分割するが、４面体や６面体の境界を組織の境界と一致させる必要があり、境界要素法と同様の困難な状況にあった。
【０１１８】
これに対して、本実施形態では、心臓や胸部の形状あるいは導電率分布を表現するため、三角形メッシュや４面体、６面体メッシュなどを組織境界に一致させたモデルを作成する必要が無く、格子状に規則的に並んだ６面体それぞれの組織を分類するのみでよい。このため、ＭＲΙやＸ−ＣＴの３次元画像から容易に心室・胴体モデルを作成できるので、検査を行う全ての患者に対して患者個々に合わせた心室モデルを容易に作成し、その心室モデルを用いて興奮到達時刻分布と活動電位振幅分布を同時に解析できる。その上、患者個々の導電率分布モデルを容易に作成し、電位計算および磁場計算に利用することができるので、興奮到達時刻と活動電位振幅の分布の解析を従来法よりも高精度に行うことができる。とくに、この実施形態の解析は、前述した各実施形態のものと比べても、解析の手間が少なくて済むから、多数の患者に対して容易に且つスループット良く、興奮到達時刻と活動電位振幅の分布の解析を行うことができ、実用性の面で高い価値を有する診断装置が提供される。
【０１１９】
その他の実施の形態
上述した各実施形態では、心臓から発生する磁場の計測値から心臓内の興奮到達時刻の分布を測定する方法を具体的に説明したが、心電計装置（図示せず）で得られる心電図データからも同様の方法に沿って心臓内の興奮到達時刻を推定することができる。
【０１２０】
また、上述した各実施形態における解析対象は、その一例として、興奮到達時刻と活動電位振幅の分布であるとして説明したが（活動電位振幅は心筋電気活動量の一態様として採用されたもの）、これも種々の変形が可能である。例えば、活動電位振幅の代わりに、心筋電気活動量として、心室内の導電率の分布、活動電位振幅と導電率の積である電流ダイポール密度の絶対値の心室内の分布、または、それらの量に関連または変形した量を推定対象としてもよい。その場合の実施方法およびその効果は興奮到達時刻分布と活動電位振幅分布を解析対象としたときと同等のものが得られる。
【０１２１】
【発明の効果】
以上説明したように本発明によれば、被験者の体表の測定点から測定した電位および磁場の少なくとも一方に基づいて、その被験者の心臓内の興奮到達時刻の分布情報とその心臓の心筋の電気活動の大きさを表す心筋電気活動量の分布情報とを同時に推定し、興奮到達時刻の分布情報と心筋電気活動量の分布情報とを表示するようにしたため、従来法のように膨大な演算時間が掛かる興奮伝播シミュレーションを行う必要がないことなどに因って、従来法よりも著しく高速に、高精度に、かつ安定して心臓内の電気生理活動を推定でき、心筋梗塞などの疾患を持つ患者など、不特定多数の患者の診断に容易に供することができる実用性および汎用性に非常に優れた心臓内電気現象の診断装置およびその診断方法を提供することができる。
【０１２２】
また、本発明の別の側面によれば、患者個々の心室形状および胸部導電率分布を患者毎に、容易に、かつ自動的に抽出できるようにし、これにより、多くの患者を簡単に診断でき、実用性を一層高めることができる。
【図面の簡単な説明】
【図１】本発明の実施形態に係る心臓内電気現象の診断装置の概略構成を示すブロック図。
【図２】第１の実施形態において実行される興奮到達時刻および活動電位振幅の分布情報の推定解析を示す概略フローチャート。
【図３】心室モデル（心室形状）の例を示す図。
【図４】興奮到達時刻および活動電位振幅の分布の表示法の例を示す図。
【図５】第２の実施形態における心室モデル（心室形状）のほかの例を示す図。
【図６】第３の実施形態において実行される興奮到達時刻および活動電位振幅の分布情報の推定解析を示す概略フローチャート。
【図７】第４の実施形態において実行される興奮到達時刻および活動電位振幅の分布情報の推定解析を示す概略フローチャート。
【図８】第４の実施形態での心臓モデル設定の処理の概要を示すフローチャート。
【図９】第１の実施形態に係る図２の処理における推定解析部分を機能的に表す図。
【図１０】従来の推定法の一例における演算の流れを示す説明図。
【図１１】解析結果である興奮到達時刻分布および活動電位振幅分布の表示の別の例を示す図。
【図１２】解析結果である興奮到達時刻分布および活動電位振幅分布の表示のさらに別の例を示す図。
【図１３】解析結果である興奮到達時刻分布および活動電位振幅分布の表示のさらに別の例を示す図。
【図１４】解析結果である興奮到達時刻分布および活動電位振幅分布の表示のさらに別の例を示す図。
【符号の説明】
１１  ＳＱＵＩＤ磁束計
１２  画像収集装置
１３  画像転送装置
１４  コンピュータ装置
１４ａ  演算処理器
１４ｂ  メモリ
１４ｃ  入力器
１４ｄ  表示器
【図１】
0003660781.tif 000018
【図２】
0003660781.tif 000019
【図３】
0003660781.tif 000020
【図４】
0003660781.tif 000021
【図５】
0003660781.tif 000022
【図６】
0003660781.tif 000023
【図７】
0003660781.tif 000024
【図８】
0003660781.tif 000025
【図９】
0003660781.tif 000026
【図１０】
0003660781.tif 000027
【図１１】
0003660781.tif 000028
【図１２】
0003660781.tif 000029
【図１３】
0003660781.tif 000030
【図１４】
0003660781.tif 000031