JP2016529054A - 構造力学上の逆問題に対する直接的な解法を有限要素法に基づいて実現するための画像処理方法 - Google Patents
構造力学上の逆問題に対する直接的な解法を有限要素法に基づいて実現するための画像処理方法 Download PDFInfo
- Publication number
- JP2016529054A JP2016529054A JP2016539616A JP2016539616A JP2016529054A JP 2016529054 A JP2016529054 A JP 2016529054A JP 2016539616 A JP2016539616 A JP 2016539616A JP 2016539616 A JP2016539616 A JP 2016539616A JP 2016529054 A JP2016529054 A JP 2016529054A
- Authority
- JP
- Japan
- Prior art keywords
- component
- matrix
- image
- node
- target object
- 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.)
- Pending
Links
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/08—Detecting organic movements or changes, e.g. tumours, cysts, swellings
- A61B8/0891—Detecting organic movements or changes, e.g. tumours, cysts, swellings for diagnosis of blood vessels
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/12—Diagnosis using ultrasonic, sonic or infrasonic waves in body cavities or body tracts, e.g. by using catheters
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/48—Diagnostic techniques
- A61B8/485—Diagnostic techniques involving measuring strain or elastic properties
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5215—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
- G16H50/20—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for computer-aided diagnosis, e.g. based on medical expert systems
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/08—Detecting organic movements or changes, e.g. tumours, cysts, swellings
- A61B8/0833—Detecting organic movements or changes, e.g. tumours, cysts, swellings involving detecting or locating foreign bodies or organic structures
- A61B8/085—Detecting organic movements or changes, e.g. tumours, cysts, swellings involving detecting or locating foreign bodies or organic structures for locating body or organic structures, e.g. tumours, calculi, blood vessels, nodules
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10132—Ultrasound image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20112—Image segmentation details
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Public Health (AREA)
- Radiology & Medical Imaging (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Theoretical Computer Science (AREA)
- Biophysics (AREA)
- Veterinary Medicine (AREA)
- Animal Behavior & Ethology (AREA)
- Surgery (AREA)
- Molecular Biology (AREA)
- Heart & Thoracic Surgery (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Vascular Medicine (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Quality & Reliability (AREA)
- Databases & Information Systems (AREA)
- Primary Health Care (AREA)
- Epidemiology (AREA)
- Data Mining & Analysis (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本発明に係る画像処理方法は、当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する工程(100)と、当該変形状態イメージをセグメント分割する工程(200)と、当該変形状態イメージを表す格子点配列構造を定義する工程(300)と、当該格子点配列構造内の各節点“i”に対して、対象物体の弾力性特性を表す未知の節点変数の対を割り当て、対象物体を構成する2つ以上の構成物質間の境界線上に局所化された一つ以上の節点に対して強化処理を行う工程(400)と、格子点配列構造内における構成要素セルの各々について構成要素マトリクスを算出する工程(500)と、全体構造を記述する一のマトリクスを取得するために、複数の構成要素マトリクスを組み合わせる工程(600)と、全体構造のマトリクスから対象物体の弾性率分布を表す画像を算出する工程(700)を含む。【選択図】図2
Description
本発明は、一般的には、構造力学における画像処理技術の分野と関係し、より具体的には、線形弾性解析における逆問題を解くための画像処理技術と関係する。本発明に係る非限定的な実施形態は、アテローム性の動脈硬化プラークが破裂する危険性の高さを評価することを目的として、アテローム性動脈硬化に起因する病変に特有の弾力性を画像で再現するために応用することが可能である。しかしながら、本発明の応用用途が上記のような特性の応用用途に限定されないことは、当該技術分野における当業者にとっては明らかに自明である。
<1>アテローム性の動脈硬化プラークが破裂することによりもたらされる結果について
アテローム又はアテローム性動脈硬化とは、大規模の血管又は平均的規模の血管(大動脈、冠状動脈、大脳動脈、下腰動脈など)の内膜の転位に対応し、脂質、複合炭水化物、血液、血液製剤、脂肪組織、石灰質沈着物およびその他の無機物などが血管内壁上で区域性集積することによって引き起こされる。
アテローム又はアテローム性動脈硬化とは、大規模の血管又は平均的規模の血管(大動脈、冠状動脈、大脳動脈、下腰動脈など)の内膜の転位に対応し、脂質、複合炭水化物、血液、血液製剤、脂肪組織、石灰質沈着物およびその他の無機物などが血管内壁上で区域性集積することによって引き起こされる。
この血管病変は、一般的には、数十年かけてゆっくりと進行するので、状態が安定し、患者にとって差し迫った危険性を示さないものであるかも知れないが、これが悪化すると、動脈硬化プラークの破裂に結び付きかねない不安定な状態となり、その場合、数日中に致死性の又は病的な心臓血管障害または脳卒中を引き起こす(CVA)可能性がある。
実際、動脈硬化プラークの破裂は、当該破裂したプラークの内容物を血管内で運ばれる血液に接触させることとなり、その結果として、血栓が形成され得る。当該血栓の形成は、血管の作用部位における血流を阻害する。また、血栓は、血管内壁から剥がれ、血流によって運ばれ、最も深刻な場合には、血管内腔を完全に詰まらせて、病変部位よりも先にある領域への血液供給を止めて、当該領域における虚血状態を生じさせる。
生体組織の力学的な特性(特にそれらの弾力性)を調べることは、医学的診断における根本的な関心事であり、特に、アテローム性の動脈硬化プラークが破裂する危険性を評価する際の根本的な関心事である。
<2>危険なアテローム性動脈硬化プラークを診断により見つけるための繊維性被膜の厚み測定について
アテローム性動脈硬化プラークが破裂する危険性を評価することを可能にするために、血管内を映像化するための幾つかの技術が開発されてきた。これらは、具体的には、
−血管内超音波検査法(IVUS:Intravascular Echography)
−光干渉断層撮影法(OCT:Optical Coherence Tomography)
−磁気共鳴映像法(IV−MRI:Intravascular Magnetic Resonance Imaging)
アテローム性動脈硬化プラークが破裂する危険性を評価することを可能にするために、血管内を映像化するための幾つかの技術が開発されてきた。これらは、具体的には、
−血管内超音波検査法(IVUS:Intravascular Echography)
−光干渉断層撮影法(OCT:Optical Coherence Tomography)
−磁気共鳴映像法(IV−MRI:Intravascular Magnetic Resonance Imaging)
これら複数の画像化技法は、アテローム性動脈硬化プラークの映像の取得およびアテローム性動脈硬化プラークにおける繊維性被膜の厚み測定を実現する可能性を与える。しかしながら、アテローム性動脈硬化プラークが破裂する危険性を評価するために、繊維性被膜の厚みを知るだけでは、アテローム性動脈硬化プラークの安定性を示す情報としては充分ではない。
<3>プラークが破裂する危険性を検出するために繊維性被膜に作用する力の強さを評価する処理について
従来からの知見によれば、繊維性被膜に作用する最大力の大きさ(PCS:ピーク荷重力)はアテローム性動脈硬化プラークの脆弱性を生体力学的に評価するための良好な指標となるという経験則を確立することが可能である。
従来からの知見によれば、繊維性被膜に作用する最大力の大きさ(PCS:ピーク荷重力)はアテローム性動脈硬化プラークの脆弱性を生体力学的に評価するための良好な指標となるという経験則を確立することが可能である。
しかしながら、繊維性被膜に作用する最大力の大きさ(PCS:ピーク荷重力)を生体内部において数値化することは、依然として解決困難な技術的課題である。何故なら、繊維性被膜に作用する力学的荷重力は、アテローム性動脈硬化プラークの形態的側面に依存するのみならず、当該プラークの構成成分(例えば、当該プラークに含まれている内包物など)の力学的性質にも依存して決まるからである。
しかし、アテローム性動脈硬化プラークの変形状態の分布を抽出するための幾つかの手法が開発済みではあるものの、アテローム性動脈硬化プラークの幾何学的構造が複雑であることに起因して、上記変形状態の分布に関する知見から当該プラークの力学的性質を直接的に推定することは実現困難である。
今まで20年間にわたって、新たな医療画像診断技法が研究開発され、これは超音波による弾性率画像計測方式と呼ばれる。上述した弾性率画像計測方式は、触診による診断と全く同じ原理に基づいて、荷重力が作用している状態の下での生体媒質の弾力性挙動を局所的に調べる。この検診法は、上記荷重力を作用させる前後において取得された又は上記荷重力を複数の異なる強さで作用させた際に取得される無線周波数帯域の超音波信号に基づいている。
<4>弾性率画像計測方式の分野においてすでに知られている従来技術について
(4.1)弾性率画像計測方式の概念について
超音波式の弾性率画像計測法は、柔らかい生体組織が変形した際に、一連の後方散乱超音波検査画像から成る時系列を処理するための技術的手段を介して当該柔らかい生体組織の弾性率を推定する演算処理から構成される。これらの画像を画像処理することにより、この生体組織の硬度や弾性率の分布を視覚的に表す画像を決定することが可能である。
(4.1)弾性率画像計測方式の概念について
超音波式の弾性率画像計測法は、柔らかい生体組織が変形した際に、一連の後方散乱超音波検査画像から成る時系列を処理するための技術的手段を介して当該柔らかい生体組織の弾性率を推定する演算処理から構成される。これらの画像を画像処理することにより、この生体組織の硬度や弾性率の分布を視覚的に表す画像を決定することが可能である。
(4.2)弾性率画像計測方式の分野において解くべき直接問題と逆問題について
標準的な構造解析計算においては、直接問題を解くための処理は、生体組織の弾性率の分布を知ることにより、当該生体組織内に作用する荷重力を予測する演算処理から構成される。より具体的には、以下に列挙する情報を知ることにより、生体組織内の任意の点位置における変位ベクトル場を決定することができるので、生体組織に加わる荷重力の空間的な分布を逆算により導き出すことが可能となり、その結果、媒質の弾性挙動についての法則性を見出すことが可能となる。
−生体組織の構造特性(すなわち、任意の点位置における生体組織の幾何学的構造)。
−生体組織を構成する複数の異なる構成物質の硬度や弾性率。
−生体組織に作用する荷重力(すなわち、生体組織に外力として加えられる力場)。
標準的な構造解析計算においては、直接問題を解くための処理は、生体組織の弾性率の分布を知ることにより、当該生体組織内に作用する荷重力を予測する演算処理から構成される。より具体的には、以下に列挙する情報を知ることにより、生体組織内の任意の点位置における変位ベクトル場を決定することができるので、生体組織に加わる荷重力の空間的な分布を逆算により導き出すことが可能となり、その結果、媒質の弾性挙動についての法則性を見出すことが可能となる。
−生体組織の構造特性(すなわち、任意の点位置における生体組織の幾何学的構造)。
−生体組織を構成する複数の異なる構成物質の硬度や弾性率。
−生体組織に作用する荷重力(すなわち、生体組織に外力として加えられる力場)。
上述した計算方法とは逆に、逆問題を解くための処理は、複数の異なる圧縮条件の下で放射された超音波信号により計測された変位ベクトル場が与えられた場合に、生体組織の弾性率を決定する演算処理から構成される。より具体的には、以下に列挙する情報を知ることにより、生体組織を構成する複数の異なる構成物質の硬度や弾性率を決定することが可能となる。
−生体組織の構造特性(これは、超音波検査画像の画像処理方式を使用して得られる)。
−生体組織に対して徐々に強めながら加えられる荷重力(これは、生体組織に対して加えられる外力を測定することにより得られる)。
−生体組織の任意の点位置における変形状態の分布すなわち変位ベクトル場(これは、複数の異なる強さの荷重力が作用している複数の時点でそれぞれ取得された超音波信号から超音波検査画像を生成し、当該超音波検査画像を表す全ての点について各点同士を対比することにより得られる)。
−生体組織の構造特性(これは、超音波検査画像の画像処理方式を使用して得られる)。
−生体組織に対して徐々に強めながら加えられる荷重力(これは、生体組織に対して加えられる外力を測定することにより得られる)。
−生体組織の任意の点位置における変形状態の分布すなわち変位ベクトル場(これは、複数の異なる強さの荷重力が作用している複数の時点でそれぞれ取得された超音波信号から超音波検査画像を生成し、当該超音波検査画像を表す全ての点について各点同士を対比することにより得られる)。
(4.3)生体組織の硬度や弾性率を決定するための既知の方法について
血管内を画像化する技術により、アテローム性動脈硬化に起因する病変部位内部における変形状態の分布場の推定結果を取得し、当該変形状態の分布場の推定結果から血管壁の弾性率分布を表すマップ情報を推定するために、有限要素法に基づく幾つかの方式がすでに開発されている。
血管内を画像化する技術により、アテローム性動脈硬化に起因する病変部位内部における変形状態の分布場の推定結果を取得し、当該変形状態の分布場の推定結果から血管壁の弾性率分布を表すマップ情報を推定するために、有限要素法に基づく幾つかの方式がすでに開発されている。
構造力学解析の分野においては、有限要素法が実行する処理内容は、解析対象となる対象物体の格子点配列イメージを生成する演算処理から構成される。格子点配列イメージを生成する演算処理により、解析対象物体を空間的に離散化した形で表現することが可能となる。この格子点配列イメージは、構成要素である複数のセルから構成され、これら複数のセルは、互いに隣接しあう格子構造(各々が三角形、多角形、正方形などの形状を有する格子構造)とも呼ばれ、各々の格子構造は、解析対象物体を形成する各部の輪郭線を定める辺と節点とから構成される(例えば、米国特許出願US2010/0160778号公報などを参照されたい)。
(下記の非特許文献1記載の従来技術について)
以下の非特許文献1記載の技術開発成果は、有限要素法に基づいて生体組織の弾性率の分布を画像で再現するための直接的な手法である。この方法は、生体組織の力学的性質についての複数の前提条件を仮定することにより、生体組織の変形状態の分布場とヤング率で表された弾性率との間の直接的な関係性を取得するものである。具体的には、この方法では、アテローム性動脈硬化プラークの力学的性質は、構成要素となる各々のセルについて一定であるとの前提を置いている。この方法はロバスト性の高い方法であるけれども、非常に非均質な構造特性を有するアテローム性動脈硬化プラークの検査期間中において、非常に数多くの構成要素セルについて演算を行わなくてはならないことにより、計算に要する時間が長くなってしまう。
以下の非特許文献1記載の技術開発成果は、有限要素法に基づいて生体組織の弾性率の分布を画像で再現するための直接的な手法である。この方法は、生体組織の力学的性質についての複数の前提条件を仮定することにより、生体組織の変形状態の分布場とヤング率で表された弾性率との間の直接的な関係性を取得するものである。具体的には、この方法では、アテローム性動脈硬化プラークの力学的性質は、構成要素となる各々のセルについて一定であるとの前提を置いている。この方法はロバスト性の高い方法であるけれども、非常に非均質な構造特性を有するアテローム性動脈硬化プラークの検査期間中において、非常に数多くの構成要素セルについて演算を行わなくてはならないことにより、計算に要する時間が長くなってしまう。
(下記の非特許文献2記載の従来技術について)
非特許文献1に関して上述した技術的限界を乗り越えるために、以下の非特許文献2では、(構造力学解析分野において解くべきあらゆるタイプの逆問題について)弾性率分布のマップ情報を画像で再現する手法が提案されており、当該手法では、有限要素法で各セルに設定される節点変数として荷重力と変位量の分布が考慮されているだけでなく、ヤング率とポアソン係数がさらに考慮されている。このようなアプローチはNMP(Nodal Mechanical Properties)と名付けられている。
非特許文献1に関して上述した技術的限界を乗り越えるために、以下の非特許文献2では、(構造力学解析分野において解くべきあらゆるタイプの逆問題について)弾性率分布のマップ情報を画像で再現する手法が提案されており、当該手法では、有限要素法で各セルに設定される節点変数として荷重力と変位量の分布が考慮されているだけでなく、ヤング率とポアソン係数がさらに考慮されている。このようなアプローチはNMP(Nodal Mechanical Properties)と名付けられている。
従って、構成要素となる一のセルが与えられた場合、非特許文献2記載の手法は、当該セルを構成する節点において弾性率変数を考慮した上で、多項式を含む内挿補間関数を使用して内挿補間処理を行うことで、構成要素となるこれらのセル内部の弾性率を推定する仕組みを提案している。
この手法の欠点は、アテローム性動脈硬化プラーク内における構成物質の不連続性を考慮に入れた演算処理を実現可能な仕組みを持たないことである。実際に、この手法では、多項式を含む内挿補間関数を用いているために、内挿補間される一連の変数間に連続性があることが必須の条件となってしまっている。従って、この手法で用いられている空間形状関数を使用したとしても、アテローム性動脈硬化プラーク中に内包される複数の異なる物質間における不連続性が著しい場合には、そのような非均質な構造特性を有するアテローム性動脈硬化プラークの特徴抽出を行うことは実現困難である。
(下記の非特許文献3記載の従来技術について)
以下の非特許文献3記載の技術開発成果は、生体組織の弾性率分布を画像で再現するために反復的に実行可能な方法であり、この方法は、IMOD(Imaging MODulography)法と名付けられている。反復実行により弾性率分布を再現するこの方法は、一般的には以下の処理工程から構成される。
−アテローム性動脈硬化プラークを構成する複数の異なる構成部分や内包物に対して弾性率の値を割り当てる処理工程。
−割り当てられた弾性率の値から生体組織における変形状態の分布場を推定する処理工程。
−生体組織の変形に関し、計算により算出された変形状態の分布場を血管内の画像化処理により実測された変形状態の分布場と対比する処理工程。
−上記計算により算出された変形状態の分布場と上記実測された変形状態の分布場との間の差異が所定の閾値を下回るまで、上述したより工程を繰り返し実行する。
以下の非特許文献3記載の技術開発成果は、生体組織の弾性率分布を画像で再現するために反復的に実行可能な方法であり、この方法は、IMOD(Imaging MODulography)法と名付けられている。反復実行により弾性率分布を再現するこの方法は、一般的には以下の処理工程から構成される。
−アテローム性動脈硬化プラークを構成する複数の異なる構成部分や内包物に対して弾性率の値を割り当てる処理工程。
−割り当てられた弾性率の値から生体組織における変形状態の分布場を推定する処理工程。
−生体組織の変形に関し、計算により算出された変形状態の分布場を血管内の画像化処理により実測された変形状態の分布場と対比する処理工程。
−上記計算により算出された変形状態の分布場と上記実測された変形状態の分布場との間の差異が所定の閾値を下回るまで、上述したより工程を繰り返し実行する。
非特許文献3記載のIMOD法は、生体組織の構造的非均質性を検出するための評価基準を用いて事前の調整を行う処理工程を備えることにより、アテローム性動脈硬化プラークを構成する全ての構成部分の輪郭線を自動的に識別することが可能である。このような評価基準を用いる処理は、パラメータ化されたモデルを使用して実行され、当該パラメータ化されたモデルは、1セグメント毎に処理分岐によって制御される。
しかしながら、IMOD法のアプローチは高い効率性とロバスト性を有するにもかかわらず、この手法は、現在時刻における弾性率分布の再現処理をリアルタイムに行うことが困難である。実際、反復実行により状態分布を再現する他の全ての方法と同様に、この手法は、数分間のオーダーに及ぶ非常に長い計算所要時間を必要とする。何故なら、アテローム性動脈硬化プラークにおけるヤング率の局所的な分布を視覚的に表すことより血管壁の弾力性について生成されるマップ情報を取得するためには、上述した計算処理工程(弾性率割り当て、変形状態の分布場の推定、対比)を非常に多くの回数にわたって繰り返し実行しなくてはならないからである。
Zhu Y他著:「ヤング率を画像で再現するための有限要素法解析に基づくアプローチ」、IEEE、医用画像技術に関する研究会議事録2003年版、第22巻、第890ページ〜第901ページ。
Oberai AA他著:「アドジョイント法を使用して弾性画像の逆問題を解くための解法」、逆問題2003年版、第19巻、第297ページ〜第313ページ。
Le Floc’h S他著:「張力測定を使用しながら、セグメント化により促進される最適化手順に基づいて、脆弱なアテローム性動脈硬化プラークの弾力性を画像で再現する技術(理論的な枠組みとしての研究)」、IEEE、医用画像技術に関する研究会議事録2009年版、第28巻、第1126ページ〜第1137ページ。
Nicolas Moes著:「格子点配列構造の再構成を行わずにヒビ割れ構造の成長過程を扱うための有限要素法」、工学分野における数値解法に関する国際学会誌、第46巻、第131ページ〜第150ページ、1999年9月10日。
MM Doyley著:「モデル・ベース弾性画像生成: 逆弾性問題へのアプローチ法に関する調査」、医学と生物学のための物理学に関する国際雑誌、2012年、第57号、R35−73
<本発明の目的>
本発明は、非特許文献1〜3に関して上述した技術的欠点を克服することを可能にする弾性率分布の再現方法を提案することを目的とする。具体的には、本発明の目的は、対象物体(例えば、アテローム性動脈硬化プラーク)に荷重力を加える前後において、または複数の異なる強さの荷重力を段階的に加えた際に取得される無線周波数帯域の超音波信号から対象物体の弾性率画像をリアルタイムに画像として再現する方法を提案することである。
本発明は、非特許文献1〜3に関して上述した技術的欠点を克服することを可能にする弾性率分布の再現方法を提案することを目的とする。具体的には、本発明の目的は、対象物体(例えば、アテローム性動脈硬化プラーク)に荷重力を加える前後において、または複数の異なる強さの荷重力を段階的に加えた際に取得される無線周波数帯域の超音波信号から対象物体の弾性率画像をリアルタイムに画像として再現する方法を提案することである。
<本発明の概略の開示>
この目的を達成するために、本発明に係る画像処理方法は、対象物体を形成する複数の構成物質に応じて、当該対象物体の弾性率画像を形成することを可能にし、当該画像処理方法は、以下の処理工程を備えることを特徴とする。
−当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する工程。
−有限要素法を適用することにより、当該変形状態イメージの格子点配列構造を定義する工程であって、当該格子点配列構造は、複数の構成要素セル(格子構造)によって構成され、各々のセルは、同一の構成物質を他の構成物質から区別する境界線を定め、少なくとも3つの節点を含み、節点の各々は、格子点配列構造内で隣接しあう一つ以上のセルに属している、ことを特徴とする工程。
−当該格子点配列構造内の各節点“i”に対して、未知の節点変数の対を割り当てる工程であって、当該未知の節点変数の対は、決定すべき弾力性特性を表す変数の対であり、対を成す各変数は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi)に対応する。
−複数の構成物質の間を区切る境界線上に位置し、当該境界線を挟んで隣接する少なくとも2つの構成要素セルにとって共通の節点である少なくとも一つの節点を検出する工程。
−当該少なくとも一つの共通の節点を強化処理するための強化処理工程であって、当該強化処理工程は、当該共通の節点に対して割り当てられた未知の節点変数の対を、強化処理用の未知変数を含む2つ以上の対で置き換える処理段階を含み、当該検出された共通の節点における強化処理用の未知変数の各対は、それぞれ対応する構成要素セルに割り当てられ、上記それぞれ対応する構成要素セルは、複数の構成物質の間を区切る境界線を挟んで隣接する少なくとも2つの構成要素セルの中から選択される、ことを特徴とする強化処理工程。
−当該格子点配列構造を構成する構成要素セルの各々について構成要素マトリクスをそれぞれ算出することによって、複数の構成要素マトリクスを取得する工程であって、一の構成要素セルについての当該構成要素マトリクスは、当該構成要素セルに割り当てられた未知の節点変数の対と、当該構成要素セルに割り当てられた強化処理用の未知変数の複数の対とを考慮して算出される、ことを特徴とする工程。
−当該複数の構成要素マトリクスから選択された一以上の構成要素マトリクスを組み合わせることにより、一の全体構造マトリクスを取得する工程。
−当該変形状態イメージと当該全体構造マトリクスから当該対象物体の弾性率分布画像を算出する工程。
この目的を達成するために、本発明に係る画像処理方法は、対象物体を形成する複数の構成物質に応じて、当該対象物体の弾性率画像を形成することを可能にし、当該画像処理方法は、以下の処理工程を備えることを特徴とする。
−当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する工程。
−有限要素法を適用することにより、当該変形状態イメージの格子点配列構造を定義する工程であって、当該格子点配列構造は、複数の構成要素セル(格子構造)によって構成され、各々のセルは、同一の構成物質を他の構成物質から区別する境界線を定め、少なくとも3つの節点を含み、節点の各々は、格子点配列構造内で隣接しあう一つ以上のセルに属している、ことを特徴とする工程。
−当該格子点配列構造内の各節点“i”に対して、未知の節点変数の対を割り当てる工程であって、当該未知の節点変数の対は、決定すべき弾力性特性を表す変数の対であり、対を成す各変数は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi)に対応する。
−複数の構成物質の間を区切る境界線上に位置し、当該境界線を挟んで隣接する少なくとも2つの構成要素セルにとって共通の節点である少なくとも一つの節点を検出する工程。
−当該少なくとも一つの共通の節点を強化処理するための強化処理工程であって、当該強化処理工程は、当該共通の節点に対して割り当てられた未知の節点変数の対を、強化処理用の未知変数を含む2つ以上の対で置き換える処理段階を含み、当該検出された共通の節点における強化処理用の未知変数の各対は、それぞれ対応する構成要素セルに割り当てられ、上記それぞれ対応する構成要素セルは、複数の構成物質の間を区切る境界線を挟んで隣接する少なくとも2つの構成要素セルの中から選択される、ことを特徴とする強化処理工程。
−当該格子点配列構造を構成する構成要素セルの各々について構成要素マトリクスをそれぞれ算出することによって、複数の構成要素マトリクスを取得する工程であって、一の構成要素セルについての当該構成要素マトリクスは、当該構成要素セルに割り当てられた未知の節点変数の対と、当該構成要素セルに割り当てられた強化処理用の未知変数の複数の対とを考慮して算出される、ことを特徴とする工程。
−当該複数の構成要素マトリクスから選択された一以上の構成要素マトリクスを組み合わせることにより、一の全体構造マトリクスを取得する工程。
−当該変形状態イメージと当該全体構造マトリクスから当該対象物体の弾性率分布画像を算出する工程。
このようにして、本発明は、線形弾性構造解析における逆問題を解くために、有限要素法による新規で直接的な定式化手法を提案している。この定式化手法は、新規な内挿補間関数を新たに採用することにより、解析対象となる対象物体の内部における構成物質間の不連続性を考慮に入れて計算を行うことを実現可能にする仕組みに基づいている。この新規な手法は、線形弾性構造解析分野における有限要素法を使用した構造解析計算を行うための任意の形態のソフトウェアに応用可能である点で技術的に優れている。
本発明に係る画像処理方法を実施するための実施形態のうち、本発明の権利範囲を限定しない好適な実施形態は、上記格子点配列構造を定義する工程に先立って変形状態イメージをセグメント分割する工程をさらに含み、それにより、この実施形態では、複数の分割区域から構成されるセグメント化された変形状態イメージが取得され、当該複数の分割区域は、対象物体を構成する各部分領域のうち、互いに異なる構成物質に含まれると考えられる部分領域を表現している。また、上記格子点配列構造を定義する工程は、セグメント分割された変形状態イメージの各々に対して実行され、当該格子点配列構造を構成する複数の構成要素セルの寸法と配置は、セグメント化された変形状態イメージを構成するそれぞれの分割区域の寸法と配置に従って決定される。
−対象物体の弾性率分布画像を算出する工程は、以下のマトリクス方程式
(数1)
{R}=([Q’]T[Q’])−1[Q’]T{F’}
を解く演算処理を備え、
ここで、[Q’]は、簡約化された全体構造マトリクスを表し、[Q’]Tは、[Q’]の転置行列であり、{R}は、対象物体の複数の節点位置における弾性率分布を示す分布場のマトリクスを表し、{F’}は、これらの節点位置に作用する荷重力の分布場の簡約化されたマトリクスを表している。
(数1)
{R}=([Q’]T[Q’])−1[Q’]T{F’}
を解く演算処理を備え、
ここで、[Q’]は、簡約化された全体構造マトリクスを表し、[Q’]Tは、[Q’]の転置行列であり、{R}は、対象物体の複数の節点位置における弾性率分布を示す分布場のマトリクスを表し、{F’}は、これらの節点位置に作用する荷重力の分布場の簡約化されたマトリクスを表している。
−検出された共通の節点を強化処理するための強化処理工程は、当該検出された共通の節点に割り当てられた未知の節点変数の一対を強化された未知変数から成るn個の対で置き換える処理を含み、nは、当該共通の節点が位置する境界線に隣接する複数の異なる構成物質の個数に相当する、ことを特徴とする。
また、本発明に係る画像処理システムは、対象物体を形成する複数の構成物質に応じて、当該対象物体の弾性率画像を形成することを可能にし、当該画像処理システムは、以下を備えることを特徴とする。
−当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する受信装置;
−一連の処理工程を実行するようにプログラミングされた演算処理回路であって、当該一連の処理工程は、
有限要素法を適用することにより、当該変形状態イメージの格子点配列構造を定義する工程であって、当該格子点配列構造は、複数の構成要素セル(格子構造)によって構成され、各々のセルは、同一の構成物質を他の構成物質から区別する境界線を定め、少なくとも3つの節点を含み、節点の各々は、格子点配列構造内で隣接しあう一つ以上のセルに属している、ことを特徴とする工程;
当該格子点配列構造内の各節点“i”に対して、未知の節点変数の対を割り当てる工程であって、当該未知の節点変数の対は、決定すべき弾力性特性を表す変数の対であり、対を成す各変数は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi)に対応する、ことを特徴とする工程;
複数の構成物質の間を区切る境界線上に位置し、当該境界線を挟んで隣接する少なくとも2つの構成要素セルにとって共通の節点である少なくとも一つの節点を検出する工程;
当該少なくとも一つの共通の節点を強化処理するための強化処理工程であって、当該強化処理工程は、当該共通の節点に対して割り当てられた未知の節点変数の対を、強化処理用の未知変数を含む2つ以上の対で置き換える処理段階を含み、当該検出された共通の節点における強化処理用の未知変数の各対は、それぞれ対応する構成要素セルに割り当てられ、上記それぞれ対応する構成要素セルは、複数の構成物質の間を区切る境界線を挟んで隣接する少なくとも2つの構成要素セルの中から選択される、ことを特徴とする強化処理工程;
当該格子点配列構造を構成する構成要素セルの各々について構成要素マトリクスをそれぞれ算出することによって、複数の構成要素マトリクスを取得する工程であって、一の構成要素セルについての当該構成要素マトリクスは、当該構成要素セルに割り当てられた未知の節点変数の対と、当該構成要素セルに割り当てられた強化処理用の未知変数の複数の対とを考慮して算出される、ことを特徴とする工程;
当該複数の構成要素マトリクスから選択された一以上の構成要素マトリクスを組み合わせることにより、一の全体構造マトリクスを取得する工程;および
当該変形状態イメージと当該全体構造マトリクスから当該対象物体の弾性率分布画像を算出する工程、を含むことを特徴とする演算処理回路。
−当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する受信装置;
−一連の処理工程を実行するようにプログラミングされた演算処理回路であって、当該一連の処理工程は、
有限要素法を適用することにより、当該変形状態イメージの格子点配列構造を定義する工程であって、当該格子点配列構造は、複数の構成要素セル(格子構造)によって構成され、各々のセルは、同一の構成物質を他の構成物質から区別する境界線を定め、少なくとも3つの節点を含み、節点の各々は、格子点配列構造内で隣接しあう一つ以上のセルに属している、ことを特徴とする工程;
当該格子点配列構造内の各節点“i”に対して、未知の節点変数の対を割り当てる工程であって、当該未知の節点変数の対は、決定すべき弾力性特性を表す変数の対であり、対を成す各変数は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi)に対応する、ことを特徴とする工程;
複数の構成物質の間を区切る境界線上に位置し、当該境界線を挟んで隣接する少なくとも2つの構成要素セルにとって共通の節点である少なくとも一つの節点を検出する工程;
当該少なくとも一つの共通の節点を強化処理するための強化処理工程であって、当該強化処理工程は、当該共通の節点に対して割り当てられた未知の節点変数の対を、強化処理用の未知変数を含む2つ以上の対で置き換える処理段階を含み、当該検出された共通の節点における強化処理用の未知変数の各対は、それぞれ対応する構成要素セルに割り当てられ、上記それぞれ対応する構成要素セルは、複数の構成物質の間を区切る境界線を挟んで隣接する少なくとも2つの構成要素セルの中から選択される、ことを特徴とする強化処理工程;
当該格子点配列構造を構成する構成要素セルの各々について構成要素マトリクスをそれぞれ算出することによって、複数の構成要素マトリクスを取得する工程であって、一の構成要素セルについての当該構成要素マトリクスは、当該構成要素セルに割り当てられた未知の節点変数の対と、当該構成要素セルに割り当てられた強化処理用の未知変数の複数の対とを考慮して算出される、ことを特徴とする工程;
当該複数の構成要素マトリクスから選択された一以上の構成要素マトリクスを組み合わせることにより、一の全体構造マトリクスを取得する工程;および
当該変形状態イメージと当該全体構造マトリクスから当該対象物体の弾性率分布画像を算出する工程、を含むことを特徴とする演算処理回路。
本発明に係る画像処理方法を実施するための実施形態のうち、本発明の権利範囲を限定しない好適な実施形態では、当該画像処理システムは、算出された弾性率分布画像を画面表示するために、ディスプレイ装置をさらに備え、当該演算処理回路は、上述した複数の実施形態に係る画像処理方法を実施することができるように構成されている。本発明はさらに、コンピュータによって読み出し可能なデータ記憶媒体上に記録されるプログラム・コードを含んだコンピュータ・プログラムとも関係し、当該プログラム・コードは、当該コンピュータ・プログラムが実行されるために当該コンピュータに読み込まれた際に、上述した複数の実施形態に係る画像処理方法を当該コンピュータに実行させることを特徴とする。
上記に加え、本発明に係るその他の技術的特徴や技術的優位性は、本発明の実施形態に関する下記の詳細な説明から明らかとなる。ここで、本発明の実施形態に関する下記の詳細な説明は、本発明の権利範囲を限定することを意図しない単に例示的なものである。また、本発明の実施形態に関する下記の詳細な説明を読むに当たっては、以下において簡単に説明する添付図面を参照されたい。
<1>はじめに
上記「背景技術」の欄で述べたように、非特許文献2は、構造力学解析における逆問題を解くために、有限要素法によって弾性率分布を画像で再現する方法を開示している。対象物体の様々な箇所での異なる弾性率を表現する画像を再現するための非特許文献2記載の方法は、節点変数として荷重力と変位量のみならず、ヤング率とポアソン係数(2つのラメ乗数により同様の数量を数値化することも可能)をさらに考慮している点が独創的である。
上記「背景技術」の欄で述べたように、非特許文献2は、構造力学解析における逆問題を解くために、有限要素法によって弾性率分布を画像で再現する方法を開示している。対象物体の様々な箇所での異なる弾性率を表現する画像を再現するための非特許文献2記載の方法は、節点変数として荷重力と変位量のみならず、ヤング率とポアソン係数(2つのラメ乗数により同様の数量を数値化することも可能)をさらに考慮している点が独創的である。
しかしながら、この方法の欠点は、連続的な空間形状関数を使用している点に関係している。この事は、構成物質間の不連続性が著しい対象物体について特徴抽出を行うことを実現困難にしている。さらに、この弾性率画像再現方法は、反復実行により解を算出する方法であるため、対象物体の複数の箇所での異なる弾性率を表す画像をリアルタイムで再現することが実現困難である。従って、本発明は、構成物質同士が不連続性である少なくとも2つの部分領域によって構成される対象物体の複数の箇所における異なる弾性率を表す画像をリアルタイムに再現することを可能にする画像処理方法を得ることを目的として創出された。
<2>本発明に係る画像処理方法を応用可能な分野について
以下、本発明に係る幾つかの例示的な実施形態について詳しく説明する。具体的には、この方法は、少なくとも一つの内包物と少なくとも一つの内腔部分を含む対象物体における弾性率分布のマップ情報を決定する用途に適用可能であり、当該対象物体は、複数の異なる構成物質によって充分な非圧縮性を有する物体として形成されているものである。
以下、本発明に係る幾つかの例示的な実施形態について詳しく説明する。具体的には、この方法は、少なくとも一つの内包物と少なくとも一つの内腔部分を含む対象物体における弾性率分布のマップ情報を決定する用途に適用可能であり、当該対象物体は、複数の異なる構成物質によって充分な非圧縮性を有する物体として形成されているものである。
以下においては、対象物体内に含まれる複数の異なる内包物の弾性率に関連する情報の断片をリアルタイムに提供する技術に基づいて本発明の幾つかの実施形態を説明する。より具体的には、本発明の幾つかの実施形態は、内腔部分を取り巻く非均質で変形可能な構造部分におけるヤング率の分布を表す弾性率画像を生成するための方法に基づいて説明されるものである。例えば、本発明の幾つかの実施形態は、アテローム性の動脈硬化プラークを有する血管の診断を適用事例として参照する形で説明される。
当該技術分野における当業者にとって明らかなように、本発明に従う画像処理方法は、内腔部分を有するか否かにかかわらず、心臓、膵臓、肝臓、前立腺、胸部などのようなその他の種別の対象物体に適用することが可能である。また、本発明に従う画像処理方法は、内腔部分を有するか否かにかかわらず、変形状態を数値することが可能で工業生産される任意の種別の人工的構造物にも適用することが可能である。同時に、本発明に係る画像処理方法を実行した結果として、ヤング率Eとポアソン係数νの分布を表すマップ情報(2つのラメ乗数λとμにより同様の数量を数値化し、当該数値化した値の分布を表す同様のマップ情報を得ることもできる)を画像で再現することが可能となる。
<3>本発明に係る画像処理方法の基本原理について
以下において後述する画像処理方法は有限要素法(FEM:Finite Element Method)に基づいて実現され、これにより、検査対象となる対象物体の格子点配列構造を取得することが可能となる。この格子点配列構造は、複数のセルによって構成され、各セルの輪郭線は、セル自体の力学的特性と共に、各構成物質の外縁を定めている。
以下において後述する画像処理方法は有限要素法(FEM:Finite Element Method)に基づいて実現され、これにより、検査対象となる対象物体の格子点配列構造を取得することが可能となる。この格子点配列構造は、複数のセルによって構成され、各セルの輪郭線は、セル自体の力学的特性と共に、各構成物質の外縁を定めている。
この画像処理方法は、対象物体の各構成物質についての節点力学特性(NMP:Nodal Mechanical Properties)に基づくアプローチ法を採用している。上述した格子点配列構造内における複数のセル間の不連続性をモデル化するために、本発明に係る画像処理方法では、強化型有限要素法(XFEM:Extended Finite Element Method)をNMPアプローチ法に対して適用可能となるように適合させたNMP−XFEM法が採用されている。
今日におけるXFEM法は、有限要素モデルにおけるヒビ割れ構造を適切に処理し、変位量の不連続性の問題に対処するための最も信頼性の高い手法である(上述した非特許文献4を参照)。ヒビ割れ構造や穿孔部分のような特異点の存在は、有限要素法(FEM)における計算収束性を著しく劣化させるので、良好な収束解を得るために、特異点近傍領域において格子点配列構造の粒度を大幅に精細化するだけでは十分ではない。このような問題を解決するために複数の異なるアプローチ法が提案されてきた。これらのアプローチ法の多くは、ヒビ割れ構造を有する対象物体のヒビ割れ構造部分における力学的挙動を表現可能な形状を変位量の関数に取り入れる仕組みに依拠している。このようなアプローチ法を実現するために用いられるXFEM法では、ヒビ割れ構造部分に位置する節点の強化処理(すなわち、節点の2重化処理)の仕組みが提案されている。このようなXFEM法をNMPアプローチ法に対して適用可能とするように適合させることは新規な着想であり、検査対象となる対象物体の内部における複数の構成物質間の構造的不連続の問題に対処することを可能にする。
本発明に従ってアテローム性動脈硬化プラークの弾性率分布を直接的に画像で再現するこの方法は、血管内超音波検査法(IVUS)を用いて患者の生体内部で撮影された6箇所の動脈病変部位に対して適用された。
<4>本発明に係る画像処理方法と関連した基本原理
図1は、例示的なアテローム性動脈硬化プラークを示す図である。この動脈硬化プラークは、動脈の断面構造の一部として示されている。図1に示す動脈の断面構造は、内側を血流が流れる内腔部分の外周縁を定める内壁部分10と血管壁の外表面を定める外壁部分20を備えている。この動脈硬化プラークはさらに、区域性集積した脂質によって充填された内部空間(壊死性コア)30を備えている。また、この動脈硬化プラークはさらに、石灰化した領域のようなその他の内包物領域40を備えている。その結果、上記の動脈硬化プラークは、一般的には、それぞれ異なる物質で構成された一つ以上の内包物(脂質、カルシウム等)を備える。
図1は、例示的なアテローム性動脈硬化プラークを示す図である。この動脈硬化プラークは、動脈の断面構造の一部として示されている。図1に示す動脈の断面構造は、内側を血流が流れる内腔部分の外周縁を定める内壁部分10と血管壁の外表面を定める外壁部分20を備えている。この動脈硬化プラークはさらに、区域性集積した脂質によって充填された内部空間(壊死性コア)30を備えている。また、この動脈硬化プラークはさらに、石灰化した領域のようなその他の内包物領域40を備えている。その結果、上記の動脈硬化プラークは、一般的には、それぞれ異なる物質で構成された一つ以上の内包物(脂質、カルシウム等)を備える。
図2を参照すると、本発明に係る画像処理方法は、以下の処理工程を含んでいる。
−当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する工程(図2の工程100)。
−任意付加的に、当該変形状態イメージをセグメント分割する工程(図2の工程200)。
−当該変形状態イメージを表す格子点配列構造を定義する工程(図2の工程300)。
−当該格子点配列構造内の各節点“i”に対して、対象物体の弾力性特性を表す未知の節点変数の対(ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi))を割り当てる工程。
−対象物体を構成する2つ以上の構成物質間の境界線上に局所化された一つ以上の節点を検出する工程。
−上記検出された節点に対して強化処理を行う工程(図2の工程400)。
−格子点配列構造内における構成要素セルの各々について構成要素マトリクスを算出する工程(図2の工程500)。
−全体構造を記述する一のマトリクスを取得するために、複数の構成要素マトリクスを組み合わせる工程(図2の工程600)。
−上記のように得られた全体構造のマトリクスから対象物体の弾性率分布を表す画像を算出する工程(図2の工程700)。
−当該対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する工程(図2の工程100)。
−任意付加的に、当該変形状態イメージをセグメント分割する工程(図2の工程200)。
−当該変形状態イメージを表す格子点配列構造を定義する工程(図2の工程300)。
−当該格子点配列構造内の各節点“i”に対して、対象物体の弾力性特性を表す未知の節点変数の対(ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi))を割り当てる工程。
−対象物体を構成する2つ以上の構成物質間の境界線上に局所化された一つ以上の節点を検出する工程。
−上記検出された節点に対して強化処理を行う工程(図2の工程400)。
−格子点配列構造内における構成要素セルの各々について構成要素マトリクスを算出する工程(図2の工程500)。
−全体構造を記述する一のマトリクスを取得するために、複数の構成要素マトリクスを組み合わせる工程(図2の工程600)。
−上記のように得られた全体構造のマトリクスから対象物体の弾性率分布を表す画像を算出する工程(図2の工程700)。
(4.1)検査対象となる対象物体を形成する複数の点位置における変形状態/変位量の分布を示すマッピング画像を受信する工程について
上述したように、本発明に係る画像処理方法は、対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変形状態の分布場を視覚的に示す変形状態イメージ(以下、「変形状態エラストグラム」と呼ぶ)を受信する工程100を備えている。
上述したように、本発明に係る画像処理方法は、対象物体内の圧力変化に応じた当該対象物体内の複数の点位置における変形状態の分布場を視覚的に示す変形状態イメージ(以下、「変形状態エラストグラム」と呼ぶ)を受信する工程100を備えている。
変形状態エラストグラムは、検査対象となる対象物体が圧縮された結果生じる内部的な変形状態を表現するものであり、ここでは、対象物体は血管組織であるので、変形状態エラストグラムは、血圧に応じて生じる血管組織の内部的な変形状態を表現する。変形状態エラストグラムは、当該技術分野における当業者にとって周知である任意の手法によって取得されてもよい。
例えば、変形状態エラストグラムは、以下のような手法によって取得されてもよい。まず、患者の動脈のうち、検査を行うことを所望する動脈部位に超音波が照射される。続いて、心拍の働きによって血管組織が拡張と収縮を繰り返すのと並行して、動脈の検査対象部位を撮影した一連の超音波検査画像の系列が取得される。最後に、変形状態エラストグラムは、一連の超音波検査画像から成る時系列の動力学的挙動を分析した結果として取得される。
具体的には、上述した動力学的挙動の分析は、以下の処理工程を備えていてもよい。
−動脈の検査対象部位を心拍周期中の2つの異なる時点で撮影した2つの異なる超音波検査画像を対比する工程(ここで、当該2つの異なる超音波検査画像は、動脈の検査対象部位に2つの異なる荷重力がそれぞれ作用している状態でそれぞれ撮影されたものである)。
−上述した2つの異なる超音波検査画像の間において複数の点または画素の各位置を対比することにより、変位量の分布を表すマップ情報を決定する工程。
−上述した2つの異なる超音波検査画像をそれぞれ取得している期間中にそれぞれ測定された動脈の荷重圧の差分を求めることにより、当該2つの異なる超音波検査画像の間における動脈の圧力差を決定する。
−動脈の検査対象部位を心拍周期中の2つの異なる時点で撮影した2つの異なる超音波検査画像を対比する工程(ここで、当該2つの異なる超音波検査画像は、動脈の検査対象部位に2つの異なる荷重力がそれぞれ作用している状態でそれぞれ撮影されたものである)。
−上述した2つの異なる超音波検査画像の間において複数の点または画素の各位置を対比することにより、変位量の分布を表すマップ情報を決定する工程。
−上述した2つの異なる超音波検査画像をそれぞれ取得している期間中にそれぞれ測定された動脈の荷重圧の差分を求めることにより、当該2つの異なる超音波検査画像の間における動脈の圧力差を決定する。
このようにして、変形状態イメージ、すなわち変形状態エラストグラムは、検査対象となる動脈部位について得られる一方で、上記のように決定された圧力差は、この変形状態イメージを得ることを可能にする。これら両方の情報断片(すなわち、検査対象となる対象物体の変形状態/変位量の分布を示すマッピング情報と上記変形状態/変位量のマッピング情報と関係付けられた圧力差)は、対象物体内の内腔部分の内壁の弾性率分布を決定するために実行される後続の処理工程に適用すべき情報として有益である。
(4.2)セグメント分割処理について
変形状態イメージをセグメント分割する工程は、輪郭線によって区切られた複数の部分領域から構成されるセグメント化された変形状態イメージを取得することを可能にする。これら複数の部分領域の各々は、検査対象となる対象物体の各エリアを表しており、より具体的には、周囲を取り囲む他のエリアとは異なる構成物質によって構成されていると想定される対象物体中のエリアである。例えば、第1の部分領域は、図1に示すアテローム性動脈硬化プラーク内において、脂質が堆積した内包物部分の輪郭を定め、第2の部分領域は、石灰質が堆積した内包物部分の輪郭を定める。具体的には、上記セグメント分割する処理工程は、図1に示すアテローム性動脈硬化プラーク内において、例えば、以下に列挙するような、複数の異なる内包物部分の外縁を定める輪郭線を検出することを可能にする。
−動脈の検査対象部位における内側の輪郭線であって、内側を血流が流れる内腔部分を定義する血管内壁10を表す輪郭線。
−動脈の検査対象部位における外側の輪郭線であって、動脈の外表面を定義する血管外壁20を表す輪郭線。
−図1に示すアテローム性動脈硬化プラーク内において、様々な内包物領域30および40の外周縁を定める輪郭線。
変形状態イメージをセグメント分割する工程は、輪郭線によって区切られた複数の部分領域から構成されるセグメント化された変形状態イメージを取得することを可能にする。これら複数の部分領域の各々は、検査対象となる対象物体の各エリアを表しており、より具体的には、周囲を取り囲む他のエリアとは異なる構成物質によって構成されていると想定される対象物体中のエリアである。例えば、第1の部分領域は、図1に示すアテローム性動脈硬化プラーク内において、脂質が堆積した内包物部分の輪郭を定め、第2の部分領域は、石灰質が堆積した内包物部分の輪郭を定める。具体的には、上記セグメント分割する処理工程は、図1に示すアテローム性動脈硬化プラーク内において、例えば、以下に列挙するような、複数の異なる内包物部分の外縁を定める輪郭線を検出することを可能にする。
−動脈の検査対象部位における内側の輪郭線であって、内側を血流が流れる内腔部分を定義する血管内壁10を表す輪郭線。
−動脈の検査対象部位における外側の輪郭線であって、動脈の外表面を定義する血管外壁20を表す輪郭線。
−図1に示すアテローム性動脈硬化プラーク内において、様々な内包物領域30および40の外周縁を定める輪郭線。
上述したセグメント分割工程は、変形状態イメージ内に含まれるデータ量を大幅に減少させ、本発明に係る画像処理方法における後続する処理工程にとって無関係な情報を除去することを可能にすると同時に、変形状態イメージが表す重要な構造特性に関する情報をそのまま維持することを可能にする。このセグメント分割工程は、重要な境界部分を境に複数のセグメントに分割する技法を用いて実現されてもよいし、当該技術分野における当業者にとって周知であるその他のセグメント分割技法を用いて実現されてもよい。
(4.3)変形状態イメージを表す格子点配列構造の定義について
変形状態イメージを表す格子点配列構造を定義する工程は、有限要素法の実施に該当する。このような格子点配列構造を定義する工程は、検査対象物体を空間的に離散化することを可能にする。変形状態イメージを表す格子点配列構造を定義する工程は、(必要に応じてセグメント分割された)変形状態イメージを複数の構成要素セルに対応する立方体の集まりとして組み上げる処理手順から構成される。複数の構成要素セルの各々は、モデル化された媒体を自身に固有の力学的特性に応じて区切るものである。複数の構成要素セルの各々は、は、複数の節点と複数の辺から成る。一つの節点は、互いに隣接する複数の構成要素セルに属している。さらに、一つの辺は、互いに隣接する2つの構成要素セルにのみ属している。
変形状態イメージを表す格子点配列構造を定義する工程は、有限要素法の実施に該当する。このような格子点配列構造を定義する工程は、検査対象物体を空間的に離散化することを可能にする。変形状態イメージを表す格子点配列構造を定義する工程は、(必要に応じてセグメント分割された)変形状態イメージを複数の構成要素セルに対応する立方体の集まりとして組み上げる処理手順から構成される。複数の構成要素セルの各々は、モデル化された媒体を自身に固有の力学的特性に応じて区切るものである。複数の構成要素セルの各々は、は、複数の節点と複数の辺から成る。一つの節点は、互いに隣接する複数の構成要素セルに属している。さらに、一つの辺は、互いに隣接する2つの構成要素セルにのみ属している。
(4.4)構成要素セルの各節点に対して未知変数の対を割り当てる処理工程について
繰り返しの説明になるが、本発明に係る画像処理方法は、検査対象物体を形成している複数の構成物質に応じて、検査対象物体の弾性率分布画像を決定することを可能にする。従って、検査対象物体を形成している様々な構成物質の弾性率は、変形状態イメージ内の任意の点位置において未知である。
繰り返しの説明になるが、本発明に係る画像処理方法は、検査対象物体を形成している複数の構成物質に応じて、検査対象物体の弾性率分布画像を決定することを可能にする。従って、検査対象物体を形成している様々な構成物質の弾性率は、変形状態イメージ内の任意の点位置において未知である。
本発明によれば、検査対象物体を構成する全ての構成要素セルの弾性率を検査するのではなく、構成要素セルの各節点“i”に対して、決定すべき弾性特性を表す一対の未知変数(λi,μi)(ラメ乗数のペアとして知られている)が割り当てられる。ここで、各節点“i”は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数から成るペア(Ei,νi)のいずれか一方を割り当てられるようにしてもよい。
様々な構成要素セルの複数の節点におけるこれらの未知変数の対と関係付けられる値が決定されれば、多項式に基づく内挿補間関数の助けを借りて構成要素セル内部における弾性特性を算出することが可能となる。より具体的には、一の構成要素セルが与えられた場合、当該構成要素セルの複数の節点において一対の未知変数の値が算出されれば、当該構成要素セルの各節点における一対の未知変数の値から内挿補間処理により、当該構成要素セル内部の任意の点位置における一対の未知変数の値を算出することができるようになる。各構成要素セルの複数の節点に対して未知変数の対を割り当てることは、構成要素セル内部における力学的特性の様々なバリエーションをモデル化可能にすると共に、弾性率分布を示すマッピング情報を画像としてより精細に再現することを可能にする。
(4.5)検出処理工程と強化処理工程について
構成要素セルの複数の節点に対して未知変数の対が一旦割り当てられたならば、本発明に係る画像処理方法は、対象物体を構成する複数の異なる構成物質がそれぞれ位置するエリアを表す2つ以上の部分領域同士の間における境界線上に局所化された共通の節点を検出する工程を実施する。続いて、当該検出された共通の節点に対して強化処理が実行される。上記のように検出された共通の節点の各々について実行される強化処理は、当該検出された共通の節点に対して割り当てられた未知の節点変数の対を強化処理された未知変数から成るn個の対で置き換える処理手順を含み、nは、複数の異なる構成物質のエリアを表す部分領域の個数に相当する。
構成要素セルの複数の節点に対して未知変数の対が一旦割り当てられたならば、本発明に係る画像処理方法は、対象物体を構成する複数の異なる構成物質がそれぞれ位置するエリアを表す2つ以上の部分領域同士の間における境界線上に局所化された共通の節点を検出する工程を実施する。続いて、当該検出された共通の節点に対して強化処理が実行される。上記のように検出された共通の節点の各々について実行される強化処理は、当該検出された共通の節点に対して割り当てられた未知の節点変数の対を強化処理された未知変数から成るn個の対で置き換える処理手順を含み、nは、複数の異なる構成物質のエリアを表す部分領域の個数に相当する。
例えば、検出された共通の節点が2つの部分領域の境界線上に位置しているならば、当該共通の節点に割り当て済みである2つの未知変数(変数の対)は、4つの未知変数(2×2個の未知変数、すなわち2対の未知変数)で置き換えられることとなる。検出された共通の節点が3つの部分領域の境界線上に位置しているならば、共通の節点に割り当て済みである2つの未知変数(変数の対)は、6つの未知変数(2×3個の未知変数、すなわち3対の未知変数)で置き換えられることとなる。検出された共通の節点が4つの部分領域の境界線上に位置しているならば、当該共通の節点に割り当て済みである2つの未知変数(変数の対)は、8つの未知変数(2×4個の未知変数、すなわち4対の未知変数)で置き換えられ、以下同様である。図3Aと図3Bは、共通の節点が2個の部分領域の境界線上および3個の部分領域の境界線上に位置している場合の強化処理の事例について示している。
上記のように検出された共通の節点に対して割り当てられた強化処理済みの未知変数から成る各対は、それぞれ対応する構成要素セルと関係付けられ、それぞれ対応する構成要素セルは、複数の異なる構成物質間の境界を区切る形で隣接する複数の構成要素セルの中から選択された構成要素セルである。
例えば、2つの異なる構成物質M1およびM2のエリアを表す2つの部分領域同士の間の境界線上に位置する共通の節点N2が検出された場合を考えると、第1の部分領域M1の輪郭の一部を定める構成要素セルの節点N10に対して強化処理済みの未知変数から成る第1の対が割り当てられ、第2の部分領域M2の輪郭の一部を定める構成要素セルの節点N11に対して強化処理済みの未知変数から成る第2の対が割り当てられることとなる。
3つの異なる構成物質M1、M2およびM3のエリアを表す3つの部分領域同士の間の境界線上に位置する共通の節点N5が検出された場合を考えると、第1の部分領域M1の輪郭の一部を定める構成要素セルの節点N14に対して強化処理済みの未知変数から成る第1の対が割り当てられ、第2の部分領域M2の輪郭の一部を定める構成要素セルの節点N15に対して強化処理済みの未知変数から成る第2の対が割り当てられ、第3の部分領域M3の輪郭の一部を定める構成要素セルの節点N16に対して強化処理済みの未知変数から成る第3の対が割り当てられることとなる。
(XFEM法と同様の基本原理に基づいてはいるが、構成要素セルの節点力学特性をモデルとして簡略化されている)強化処理工程は、検査対象となる対象物体を形成している複数の構成物質(この例においては、アテローム性動脈硬化プラークの複数の内包物)間における構造的不連続性を適切に処理することを可能にする。
(4.6)構成要素セルの構成要素マトリクスを算出し、これらを組み合わせる工程について
幾つかの部分領域の境界線上に位置する共通の節点に対する強化処理が一旦行われたならば、本発明に係る画像処理方法は、構成要素セルの各々について構成要素マトリクスを算出する工程を実施する。
幾つかの部分領域の境界線上に位置する共通の節点に対する強化処理が一旦行われたならば、本発明に係る画像処理方法は、構成要素セルの各々について構成要素マトリクスを算出する工程を実施する。
この構成要素マトリクスは、複数の節点の座標位置、複数の節点における変位量(変形状態)の分布場および構成物質の節点定数と節点変位について使用される形状関数を考慮に入れて算出される。本明細書の末尾の追補の章では、3つの節点を有する三角形形状の構成要素セルと4つの節点を有する矩形形状の構成要素セルを用いて、構成要素マトリクスを算出するための2つの例について説明している。
全ての構成要素セルについて構成要素マトリクスが算出されたならば、これらの構成要素マトリクスは、有限要素法において用いられる標準的かつ当業者にとって周知な組み合わせ方法に従って組み合わせ処理がされる。このようにして、全体構造を記述する一のマトリクスが得られる。
この全体構造マトリクスおよび複数の節点に加えられる荷重力の分布場は、格子点配列構造内の複数の異なるセルに属する複数の節点における弾性率を算出するのに用いられる。その結果として、アテローム性の動脈硬化プラークを形成している複数の異なる構成物質についての弾性率の分布を示すマップ情報が得られる。続いて、この弾性率分布を表すマップ情報は、コンピュータの画面のような表示手段によって画面表示され(図2の工程800)、それによって利用者が診断を実施することを可能にする(より具体的には、破裂する可能性のあるアテローム性動脈硬化プラークを利用者が識別することを可能にする)。
技術的に優れている点として、本発明に従って上述した画像処理方法は、プロセッサを備えるコンピュータ上でアプリケーション・ソフトウェアとして実行可能とするために、プログラム・コードの形で実装することが可能であり、その際、上述したコンピュータのプロセッサは、本発明に従って上述した一連の処理工程を実行するように構成されることとなる。
この画像処理方法は、例えば、IVUSなどの手法を用いて利用者が取得した超音波信号からアテローム性動脈硬化プラークの弾性率分布を示すマップ画像をコンピュータのディスプレイ装置上にリアルタイムで画面表示することを可能にする。以下、本発明に係る画像処理方法の理論的な側面についてより詳しく説明する。
<5>本発明に係る画像処理方法と関連する理論
(5.1)有限要素モデルで表された構造を算出するための新規な予測定式化について
まず、構成物質の構造的不連続性に対して行われる強化処理の概念を導入するために、以下の説明においては、以下のような単純な動脈硬化プラークを題材として考慮する。この動脈硬化プラークは、3つの別々の構成物質から構成され、各々が4つの節点を含む4つの構成要素セルから成る格子点配列構造を用いて格子点配列としてモデル化されている(この事例は、図3Aおよび図3Bにおいて図示されている)。
(5.1)有限要素モデルで表された構造を算出するための新規な予測定式化について
まず、構成物質の構造的不連続性に対して行われる強化処理の概念を導入するために、以下の説明においては、以下のような単純な動脈硬化プラークを題材として考慮する。この動脈硬化プラークは、3つの別々の構成物質から構成され、各々が4つの節点を含む4つの構成要素セルから成る格子点配列構造を用いて格子点配列としてモデル化されている(この事例は、図3Aおよび図3Bにおいて図示されている)。
ガレルキン有限要素法に基づいて、この構造モデルに対して適用される弾性方程式をマトリクス形式で表すと以下の式のようになる。
(数2)
[K(λ、μ)]{U}={F} 式(1)
ここで、{F}は、複数の節点に作用する荷重力を表し、{U}は、複数の節点の変位量を表す変位ベクトルであり、λとμは、それぞれの構成要素セルに対応する構成物質の弾性特性と関連付けられたラメ乗数であり、[K]は、全体構造に関する対称行列である剛性マトリクスを表す。
(数2)
[K(λ、μ)]{U}={F} 式(1)
ここで、{F}は、複数の節点に作用する荷重力を表し、{U}は、複数の節点の変位量を表す変位ベクトルであり、λとμは、それぞれの構成要素セルに対応する構成物質の弾性特性と関連付けられたラメ乗数であり、[K]は、全体構造に関する対称行列である剛性マトリクスを表す。
また、図3Aおよび図3Bに示すような規則的な格子点配列構造と関連付けられた複数の節点の変位量を有限要素モデルで近似すると、以下の式のようになる。
ここで、uとvは、複数の節点の変位量を表す変位量ベクトルの2つの方向成分であり、uiとviは、上記2つの方向成分について節点“i”における節点変位量を表し、
は、節点“i”と関連付けられた変位量に対して適用される形状関数であり、
は、節点の位置を表す位置ベクトルである。
ここで、uとvは、複数の節点の変位量を表す変位量ベクトルの2つの方向成分であり、uiとviは、上記2つの方向成分について節点“i”における節点変位量を表し、
上述した技法における核心的な技術的課題は、強化処理の対象とすべき節点を如何にして適切に選択するかという点と、そのような節点に適用すべき強化処理用の形状関数の形をどのようにするかという点である。
以下に述べる節点強化法は、アテローム性動脈硬化プラークを構成する複数の異なる構成物質同士の間の境界に位置する節点の自由度を高める追加的な節点を新たに設けることにより、アテローム性動脈硬化プラークを構成する複数の異なる構成物質に固有の構造的不連続性による制約条件を満足することを可能にする。この節点強化法によれば、本発明に係る画像処理方法において、ある節点が複数の構成物質間の境界線上に位置しているならば、当該節点に対して、強化処理が行われなくてはならないと見なされる。
図3Aと図3Bは、動脈硬化プラークの単純な一例に対して上述した規則を適用したケースを図示している。この例において、節点N2、N4、N5およびN6に対して強化処理が行われている。従って、図3Aと図3Bにおける格子点配列構造と関連付けられた構成物質に関する2つの定数(すなわち、構成要素セルの2つのラメ乗数λとμ)を有限要素モデルで近似すると以下の式のようになる。
上記の式(3a)と(3b)において、
ここで、φiは、節点“i”と関連する力学的特性についての形状関数であり、ei(1≦i≦4)は、構成要素セルの識別番号であり、Iは、強化処理がされていない節点の全てを要素として含む集合(つまり、I={節点N1、N3、N7、N8およびN9})であり、Jは、強化処理がされた節点の位置における関連する節点を要素として含む集合(つまり、J={節点N10,…,N18})であり、λkとμk(k∈I,J)は、節点物性ベクトル{λ}と{μ}のそれぞれに含まれるベクトル要素である。
形状関数φi、
および、
形状関数
の両者が相違する点は、構成物質の弾性特性に関して使用されるか、それとも変位量の分布場に関して使用されるかという点である。例えば、力学的特性に関する3つの節点と変位量に関する6つの節点について構成要素セルの同一の幾何学的形状が考慮されるようにしてもよい。しかしながら、当該技術分野における当業者にとって明らかなように、物質特性と変位量の分布場の両者について全く同一の形状関数を用いるようにしてもよい。
および、
形状関数
(5.2)逆問題として定義された有限要素モデルを使用して構造を算出するための新たな定式化について
予測型有限要素法(PMN−XFEM法)に基づく上記の方程式は、節点変位量(つまり、節点における変位量)と関連する内挿補間関数および構成物質中の節点の力学的特性と関連する内挿補間関数を考慮に入れて書き直される。
予測型有限要素法(PMN−XFEM法)に基づく上記の方程式は、節点変位量(つまり、節点における変位量)と関連する内挿補間関数および構成物質中の節点の力学的特性と関連する内挿補間関数を考慮に入れて書き直される。
続いて、予測型有限要素法(PMN−XFEM法)に基づくこれらの同じ方程式を再度定式化し直すことによって、節点の力学的特性を弾性解析の逆問題における未知の数量を表すベクトル{R}として分離する。このようにして、変位量を表すマトリクスQが抽出され、節点における力学的特性を表すベクトル{R}を決定するために使用される。その結果、有限要素法に基づいて定義された逆問題における構造を算出するための新規な手法について、以下の式(4)に示すような複数の線形平衡状態方程式から成る組が導かれる。
変形の態様が平面変形である(つまり、動脈硬化プラークの全ての点位置における変位量が同一断面内に分布しているような場合(本明細書の末尾の「追補」の章を参照))と想定した場合、3つの節点から成る三角形状の格子構造および4つの節点から成る矩形形状の格子構造について変位量を表す構成要素マトリクスの表現[Qe]が与えられる。変位量を表す構成要素マトリクス[Qe]は、正方行列であるけれども、対称行列ではない。変位量を表す全体構造マトリクス[Q]を得るために、複数の異なる構成要素セルの変位量に関する構成要素マトリクス[Qe]は、当業者にとって周知な方法に従って組み合わせられ、このように構成要素マトリクス[Qe]を組み合せた結果から、全体構造を表す剛性マトリクス[K]が推定される。
(5.3)数値解法について
変位量を表す全体構造マトリクス[Q]は、矩形行列であるのに対して剛性マトリクス[K]は、正方行列である。より具体的には、2つのマトリクス[Q]および[K]は共に同じ数の行を有しているが列の個数が以下のように異なる。つまり、変位量を表す全体構造マトリクス[Q]に含まれる列の数は、強化処理がされた節点位置において新たに追加された節点の個数を2倍した数と、マトリクス[Q]自身に含まれる行の個数と、を加算した合計値に等しい。他方、節点変位量と関係した未知変数の代わりに、節点変形状態に関連した未知変数を使用して同様の系を導くようにしてもよい。
変位量を表す全体構造マトリクス[Q]は、矩形行列であるのに対して剛性マトリクス[K]は、正方行列である。より具体的には、2つのマトリクス[Q]および[K]は共に同じ数の行を有しているが列の個数が以下のように異なる。つまり、変位量を表す全体構造マトリクス[Q]に含まれる列の数は、強化処理がされた節点位置において新たに追加された節点の個数を2倍した数と、マトリクス[Q]自身に含まれる行の個数と、を加算した合計値に等しい。他方、節点変位量と関係した未知変数の代わりに、節点変形状態に関連した未知変数を使用して同様の系を導くようにしてもよい。
ところで、逆問題について考えるならば、変位量の分布場を視覚的に表す変形状態イメージから全ての節点についての変位量(ukおよびvk)が導かれれば、マトリクス[Q]は既知となる。しかしながら、これだけでは、変位量が加えられた節点位置において節点に加わる荷重力が未知であるので、逆問題を解くための条件付けは不充分である。その結果、全体構造を表す荷重力ベクトル{F}においては、幾つかの要素の値が未知である。
全体構造を表す有限要素系におけるこれら未知の荷重力成分を包含する方程式を縮約することにより、以下の式(5)で表すように、複数の線型方程式から成る簡約化された系が導かれる。
(数20)
[Q’]{R}={F’} 式(5)
ここで、{F’}は、簡約化された節点荷重力ベクトルを表し、[Q’]は、簡約化された変位量マトリクスを表す。
(数20)
[Q’]{R}={F’} 式(5)
ここで、{F’}は、簡約化された節点荷重力ベクトルを表し、[Q’]は、簡約化された変位量マトリクスを表す。
矩形行列に基づくこの線形系から解を求めるために、当業者にとって周知な数値解析アルゴリズムが使用される。例えば、上記式(5)で表される方程式の左辺と右辺にある項に対して簡約化された変位量マトリクスを乗算して行列積を求めると、正方行列に基づく線形系として以下の式が導かれる。
(数21)
[Q’]T[Q’]{R}=[Q’]T{F’} 式(6)
続いて、上記のように導かれた線形系を解くことにより、複数の節点位置における力学的特性が以下のように決定される。
(数22)
{R}=([Q’]T[Q’])−1[Q’]T{F’} 式(7)
(数21)
[Q’]T[Q’]{R}=[Q’]T{F’} 式(6)
続いて、上記のように導かれた線形系を解くことにより、複数の節点位置における力学的特性が以下のように決定される。
(数22)
{R}=([Q’]T[Q’])−1[Q’]T{F’} 式(7)
以上のように、本発明に従って上述した画像処理方法は、利用者がアテローム性の動脈硬化プラークが破裂する危険性を予測できるようにするために活用可能な情報の断片をリアルタイムに提供することを実現可能にする。
<6>追補
(6.1)力学的特性と変位量の両者を求めるために3つの節点から成る三角形状を有する同一の構成要素セルが題材として考慮される場合に、変位量を表す構成要素マトリクス[Qe]を表すための数学的表現について
変形の態様が平面変形である(つまり、動脈硬化プラークの全ての点位置における変位量が構造の同一平面内に分布しているような場合)と想定した場合、この例では、3つの節点(節点1、2および3)を有する三角形状の構成要素セルについて、変位量を表す構成要素マトリクス[Qe]として、正方行列ではあるが対称行列ではない行列を表す数学的表現が以下のように与えられる。ここで、(xi,yi)および(ui,vi)は、節点“i”(1≦i≦3)の座標位置および変位量をそれぞれ表すものとする。
(6.1)力学的特性と変位量の両者を求めるために3つの節点から成る三角形状を有する同一の構成要素セルが題材として考慮される場合に、変位量を表す構成要素マトリクス[Qe]を表すための数学的表現について
変形の態様が平面変形である(つまり、動脈硬化プラークの全ての点位置における変位量が構造の同一平面内に分布しているような場合)と想定した場合、この例では、3つの節点(節点1、2および3)を有する三角形状の構成要素セルについて、変位量を表す構成要素マトリクス[Qe]として、正方行列ではあるが対称行列ではない行列を表す数学的表現が以下のように与えられる。ここで、(xi,yi)および(ui,vi)は、節点“i”(1≦i≦3)の座標位置および変位量をそれぞれ表すものとする。
正方行列ではあるが対称行列ではなく、変位量を表す6行×6列の構成要素マトリクス[Qe]に含まれる36個の行列要素
は、xij,yij(ここで、xij=xj−xiかつyij=yj−yiであり、1≦i,j≧3)の関数として以下の式のように与えられ、以下の式においてAは三角形状の構成要素セルの表面積(すなわち、A=(x21y31−y21x31)/2)である。
構成要素マトリクス[Qe]の第1行目は以下のように与えられる。
構成要素マトリクス[Qe]の第2行目は以下のように与えられる。
構成要素マトリクス[Qe]の第3行目は以下のように与えられる。
構成要素マトリクス[Qe]の第4行目は以下のように与えられる。
構成要素マトリクス[Qe]の第5行目は以下のように与えられる。
構成要素マトリクス[Qe]の第6行目は以下のように与えられる。
(6.2)力学的特性と変位量の両者を求めるために4つの節点から成る矩形形状を有する同一の構成要素セルが題材として考慮される場合に、変位量を表す構成要素マトリクス[Qe]を表すための数学的表現について
変形の態様が平面変形である(つまり、動脈硬化プラークの全ての点位置における変位量が構造の同一平面内に分布しているような場合)と想定した場合、この例では、4つの節点(それぞれ位置座標(−a,−b)、(a,−b)、(a,b)および(−a,b)を有する4つの節点1、2、3および4)を有する矩形形状の構成要素セルについて、変位量を表す構成要素マトリクス[Qe]として、正方行列ではあるが対称行列ではない行列を表す数学的表現が以下のように与えられる。ここで、(ui,vi)は、節点“i”(1≦i≦4)の変位量をそれぞれ表すものとする。
変形の態様が平面変形である(つまり、動脈硬化プラークの全ての点位置における変位量が構造の同一平面内に分布しているような場合)と想定した場合、この例では、4つの節点(それぞれ位置座標(−a,−b)、(a,−b)、(a,b)および(−a,b)を有する4つの節点1、2、3および4)を有する矩形形状の構成要素セルについて、変位量を表す構成要素マトリクス[Qe]として、正方行列ではあるが対称行列ではない行列を表す数学的表現が以下のように与えられる。ここで、(ui,vi)は、節点“i”(1≦i≦4)の変位量をそれぞれ表すものとする。
正方行列ではあるが対称行列ではなく、変位量を表す8行×8列の構成要素マトリクス[Qe]に含まれる64個の行列要素
は、以下の式のように与えられる。
構成要素マトリクス[Qe]の第1行目は以下のように与えられる。
構成要素マトリクス[Qe]の第2行目は以下のように与えられる。
構成要素マトリクス[Qe]の第3行目は以下のように与えられる。
構成要素マトリクス[Qe]の第4行目は以下のように与えられる。
構成要素マトリクス[Qe]の第5行目は以下のように与えられる。
構成要素マトリクス[Qe]の第6行目は以下のように与えられる。
構成要素マトリクス[Qe]の第7行目は以下のように与えられる。
構成要素マトリクス[Qe]の第8行目は以下のように与えられる。
Claims (9)
- 対象物体を形成する複数の構成物質に応じて、前記対象物体の弾性率画像を形成することを可能にする画像処理方法であって、
前記対象物体内の圧力変化に応じた前記対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する工程と、
有限要素法を適用することにより、前記変形状態イメージの格子点配列構造を定義する工程であって、前記格子点配列構造は、複数の構成要素セル(格子構造)によって構成され、各々の前記構成要素セルは、同一の構成物質を他の構成物質から区別する境界線を定め、少なくとも3つの節点を含み、前記節点の各々は、前記格子点配列構造内で隣接しあう一つ以上の前記構成要素セルに属している、ことを特徴とする工程と、
前記格子点配列構造内の各節点“i”に対して、未知の節点変数の対を割り当てる工程であって、前記未知の節点変数の対は、決定すべき弾力性特性を表す変数の対であり、対を成す変数は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi)に対応する、ことを特徴とする工程と、
複数の構成物質の間を区切る境界線上に位置し、前記境界線を挟んで隣接する少なくとも2つの前記構成要素セルにとって共通の節点である少なくとも一つの節点を検出する工程と、
前記少なくとも一つの共通の節点を強化処理するための強化処理工程であって、前記強化処理工程は、前記共通の節点に対して割り当てられた前記未知の節点変数の対を、強化処理用の未知変数を含む2つ以上の対で置き換える処理段階を含み、前記検出された共通の節点における強化処理用の未知変数の各対は、それぞれ対応する前記構成要素セルに割り当てられ、それぞれ対応する前記構成要素セルは、複数の構成物質の間を区切る境界線を挟んで隣接する少なくとも2つの前記構成要素セルの中から選択される、ことを特徴とする強化処理工程と、
前記格子点配列構造を構成する前記構成要素セルの各々について構成要素マトリクスをそれぞれ算出することによって、複数の構成要素マトリクスを取得する工程であって、一の前記構成要素セルについての前記構成要素マトリクスは、前記構成要素セルに割り当てられた前記未知の節点変数の対と、前記構成要素セルに割り当てられた強化処理用の未知変数の複数の対とを考慮して算出される、ことを特徴とする工程と、
前記複数の構成要素マトリクスから選択された一以上の構成要素マトリクスを組み合わせることにより、一の全体構造マトリクスを取得する工程と、
前記変形状態イメージと前記全体構造マトリクスから前記対象物体の弾性率分布画像を算出する工程と、を備える、画像処理方法。 - 前記格子点配列構造を定義する工程に先立って、
前記変形状態イメージをセグメント分割する工程をさらに含み、
前記セグメント分割する工程においては、複数の分割区域から構成されるセグメント化された変形状態イメージが取得され、前記複数の分割区域は、対象物体を構成する各部分領域のうち、互いに異なる構成物質に含まれると考えられる部分領域を表現している、
ことを特徴とする請求項1記載の画像処理方法。 - 前記格子点配列構造を定義する工程は、セグメント分割された前記変形状態イメージの各々に対して実行され、前記格子点配列構造を構成する複数の構成要素セルの寸法と配置は、セグメント化された前記変形状態イメージを構成するそれぞれの分割区域の寸法と配置に従って決定される、
ことを特徴とする請求項2記載の画像処理方法。 - 対象物体の弾性率分布画像を算出する工程は、以下のマトリクス方程式
(数1)
{R}=([Q’]T[Q’])−1[Q’]T{F’}
を解く演算処理を備え、
ここで、[Q’]は、簡約化された全体構造マトリクスを表し、[Q’]Tは、[Q’]の転置行列であり、{R}は、対象物体の複数の節点位置における弾性率分布を示す分布場のマトリクスを表し、{F’}は、これらの節点位置に作用する荷重力の分布場の簡約化されたマトリクスを表す、
ことを特徴とする、請求項1乃至請求項3のいずれか一項に記載の画像処理方法。 - 前記検出された共通の節点を強化処理するための前記強化処理工程は、前記検出された共通の節点に割り当てられた前記未知の節点変数の一対を前記強化された未知変数から成るn個の対で置き換える処理を含み、nは、前記共通の節点が位置する境界線に隣接する複数の異なる構成物質の個数に相当する、
ことを特徴とする、請求項1乃至請求項4のいずれか一項に記載の画像処理方法。 - 対象物体を形成する複数の構成物質に応じて、前記対象物体の弾性率画像を形成することを可能にする画像処理システムであって、
−前記対象物体内の圧力変化に応じた前記対象物体内の複数の点位置における変位量または変形状態の分布場を視覚的に示す変形状態イメージを受信する受信装置と、
−一連の処理工程を実行するようにプログラミングされた演算処理回路と、を備え、
前記一連の処理工程は、
有限要素法を適用することにより、前記変形状態イメージの格子点配列構造を定義する工程であって、前記格子点配列構造は、複数の構成要素セル(格子構造)によって構成され、各々の前記構成要素セルは、同一の構成物質を他の構成物質から区別する境界線を定め、少なくとも3つの節点を含み、節点の各々は、格子点配列構造内で隣接しあう一つ以上の前記構成要素セルに属している、ことを特徴とする工程;
前記格子点配列構造内の各節点“i”に対して、未知の節点変数の対を割り当てる工程であって、前記未知の節点変数の対は、決定すべき弾力性特性を表す変数の対であり、対を成す変数は、ラメ乗数のペア(λi,μi)またはヤング率とポアソン係数のペア(Ei,νi)に対応する、ことを特徴とする工程;
複数の構成物質の間を区切る境界線上に位置し、前記境界線を挟んで隣接する少なくとも2つの前記構成要素セルにとって共通の節点である少なくとも一つの節点を検出する工程;
前記少なくとも一つの共通の節点を強化処理するための強化処理工程であって、前記強化処理工程は、前記共通の節点に対して割り当てられた前記未知の節点変数の対を、強化処理用の未知変数を含む2つ以上の対で置き換える処理段階を含み、前記検出された共通の節点における強化処理用の未知変数の各対は、それぞれ対応する構成要素セルに割り当てられ、それぞれ対応する前記構成要素セルは、複数の構成物質の間を区切る境界線を挟んで隣接する少なくとも2つの前記構成要素セルの中から選択される、ことを特徴とする強化処理工程;
前記格子点配列構造を構成する前記構成要素セルの各々について構成要素マトリクスをそれぞれ算出することによって、複数の構成要素マトリクスを取得する工程であって、一の前記構成要素セルについての前記構成要素マトリクスは、前記構成要素セルに割り当てられた未知の節点変数の対と、前記構成要素セルに割り当てられた強化処理用の未知変数の複数の対とを考慮して算出される、ことを特徴とする工程;
前記複数の構成要素マトリクスから選択された一以上の構成要素マトリクスを組み合わせることにより、一の全体構造マトリクスを取得する工程;および
前記変形状態イメージと前記全体構造マトリクスから前記対象物体の弾性率分布画像を算出する工程、
を含むことを特徴とする画像処理システム。 - 算出された前記弾性率分布画像を画面表示するために、ディスプレイ装置をさらに備えることを特徴とする、請求項6記載の画像処理システム。
- 前記演算処理回路は、請求項2乃至請求項5のいずれか一項に記載された画像処理方法を実施することができるように構成されている、
ことを特徴とする、請求項6または請求項7に記載の画像処理システム。 - コンピュータによって読み出し可能なデータ記憶媒体上に記録されるプログラム・コードを含んだコンピュータ・プログラムであって、
前記プログラム・コードは、前記コンピュータ・プログラムが実行のために前記コンピュータに読み込まれた際に、請求項1乃至請求項5のいずれか一項に記載された画像処理方法を前記コンピュータに実行させることを特徴とする、コンピュータ・プログラム。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
FR1358425 | 2013-09-03 | ||
FR1358425A FR3010218B1 (fr) | 2013-09-03 | 2013-09-03 | Procede de traitement d'image base sur la technique des elements finis pour la resolution directe de problemes inverses en mecanique des structures. |
PCT/FR2014/052172 WO2015033058A1 (fr) | 2013-09-03 | 2014-09-03 | Procédé de traitement d'image basé sur la technique des éléments finis pour la résolution directe de problèmes inverses en mécanique des structures |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2016529054A true JP2016529054A (ja) | 2016-09-23 |
Family
ID=49949791
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2016539616A Pending JP2016529054A (ja) | 2013-09-03 | 2014-09-03 | 構造力学上の逆問題に対する直接的な解法を有限要素法に基づいて実現するための画像処理方法 |
Country Status (5)
Country | Link |
---|---|
US (1) | US10007983B2 (ja) |
EP (1) | EP3041417A1 (ja) |
JP (1) | JP2016529054A (ja) |
FR (1) | FR3010218B1 (ja) |
WO (1) | WO2015033058A1 (ja) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
FR3002672B1 (fr) * | 2013-02-22 | 2016-10-07 | Univ Joseph Fourier - Grenoble 1 | Procede de generation d'une image d'elasticite |
JP6125281B2 (ja) * | 2013-03-06 | 2017-05-10 | 東芝メディカルシステムズ株式会社 | 医用画像診断装置、医用画像処理装置及び制御プログラム |
JP6707249B2 (ja) * | 2015-03-13 | 2020-06-10 | 学校法人 名城大学 | マイクロ断層可視化方法およびシステム |
EP3282949B1 (en) * | 2015-04-13 | 2019-06-05 | Case Western Reserve University | Dual energy x-ray coronary calcium grading |
WO2018133098A1 (zh) | 2017-01-23 | 2018-07-26 | 上海联影医疗科技有限公司 | 血管壁应力应变状态获取方法及系统 |
EP3828824A1 (en) * | 2019-11-28 | 2021-06-02 | Dassault Systèmes | Polyline contributor in civil engineering |
CN114241156A (zh) * | 2021-12-15 | 2022-03-25 | 上海交通大学医学院附属第九人民医院 | 用于仿真软组织变形的装置、仿真系统 |
CN116487038B (zh) * | 2023-06-25 | 2023-08-18 | 四川大学华西医院 | 轻度认知障碍向阿尔茨海默发展的预测系统和存储介质 |
Family Cites Families (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP0908137A1 (en) | 1997-10-06 | 1999-04-14 | Technologiestichting STW | A method and apparatus for making an image of a lumen or other body cavity and its surrounding tissue |
US6277074B1 (en) * | 1998-10-02 | 2001-08-21 | University Of Kansas Medical Center | Method and apparatus for motion estimation within biological tissue |
US20040009459A1 (en) * | 2002-05-06 | 2004-01-15 | Anderson James H. | Simulation system for medical procedures |
JP4263943B2 (ja) | 2003-05-07 | 2009-05-13 | テルモ株式会社 | 超音波診断装置 |
US7318804B2 (en) * | 2003-12-09 | 2008-01-15 | The Regents Of The University Of Michigan | Methods and systems for measuring mechanical property of a vascular wall and method and system for determining health of a vascular structure |
ES2379468T3 (es) * | 2004-08-24 | 2012-04-26 | The General Hospital Corporation | Procedimiento, sistema y configuración de software para determinar el módulo de elasticidad |
WO2009064715A1 (en) * | 2007-11-14 | 2009-05-22 | Auckland Uniservices Limited | Method for multi-scale meshing of branching biological structures |
US8394026B2 (en) * | 2008-11-03 | 2013-03-12 | University Of British Columbia | Method and apparatus for determining viscoelastic parameters in tissue |
FR2938957B1 (fr) * | 2008-11-21 | 2011-01-21 | Univ Joseph Fourier Grenoble I | Procede de traitement d'image pour l'estimation d'un risque de rupture de plaque d'atherome |
WO2012162413A2 (en) * | 2011-05-23 | 2012-11-29 | University Of Pittsburgh-Of The Commonwealth System Of Higher Education | Methods and apparatuses for measuring tissue stiffness changes using ultrasound elasticity imaging |
-
2013
- 2013-09-03 FR FR1358425A patent/FR3010218B1/fr not_active Expired - Fee Related
-
2014
- 2014-09-03 US US14/915,690 patent/US10007983B2/en not_active Expired - Fee Related
- 2014-09-03 JP JP2016539616A patent/JP2016529054A/ja active Pending
- 2014-09-03 WO PCT/FR2014/052172 patent/WO2015033058A1/fr active Application Filing
- 2014-09-03 EP EP14777693.4A patent/EP3041417A1/fr not_active Withdrawn
Also Published As
Publication number | Publication date |
---|---|
FR3010218B1 (fr) | 2016-12-30 |
WO2015033058A1 (fr) | 2015-03-12 |
US20160196645A1 (en) | 2016-07-07 |
EP3041417A1 (fr) | 2016-07-13 |
US10007983B2 (en) | 2018-06-26 |
FR3010218A1 (fr) | 2015-03-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP2016529054A (ja) | 構造力学上の逆問題に対する直接的な解法を有限要素法に基づいて実現するための画像処理方法 | |
Li et al. | Stress analysis of carotid plaque rupture based on in vivo high resolution MRI | |
Ohayon et al. | Necrotic core thickness and positive arterial remodeling index: emergent biomechanical factors for evaluating the risk of plaque rupture | |
Xiong et al. | Simulation of blood flow in deformable vessels using subject‐specific geometry and spatially varying wall properties | |
van Rietbergen | Micro-FE analyses of bone: state of the art | |
Urban et al. | Measurement of viscoelastic properties of in vivo swine myocardium using lamb wave dispersion ultrasound vibrometry (LDUV) | |
JP4884528B2 (ja) | 腔部の画像を評定する方法、装置ならびにコンピュータプログラム製品 | |
Creane et al. | Finite element modelling of diseased carotid bifurcations generated from in vivo computerised tomographic angiography | |
Mei et al. | Estimating the non-homogeneous elastic modulus distribution from surface deformations | |
Leischik et al. | Intraobserver and interobserver reproducibility for radial, circumferential and longitudinal strain echocardiography | |
Mei et al. | Regularizing biomechanical maps for partially known material properties | |
US8660326B2 (en) | Image processing method for estimating a risk of atheroma plaque breakage | |
Auricchio et al. | A clinically applicable stochastic approach for noninvasive estimation of aortic stiffness using computed tomography data | |
Bracco et al. | Fast strain mapping in abdominal aortic aneurysm wall reveals heterogeneous patterns | |
WO2009081599A1 (ja) | 画像処理装置、画像処理プログラム、記憶媒体及び超音波診断装置 | |
Morita et al. | Elastic modulus estimation based on local displacement observation of elastic body | |
KR101530352B1 (ko) | 물질특성에 기반한 전산유체역학 모델링 및 분석 방법 | |
Poloni et al. | Velocity vector comparison between vector flow imaging and computational fluid dynamics in the carotid bifurcation | |
Aglyamov et al. | Model‐based reconstructive elasticity imaging using ultrasound | |
Peters et al. | Digital image elasto-tomography: combinatorial and hybrid optimization algorithms for shape-based elastic property reconstruction | |
JP6367837B2 (ja) | 弾性率分布画像の生成方法及びコンピュータ・プログラム | |
Veress et al. | The direct incorporation of perfusion defect information to define ischemia and infarction in a finite element model of the left ventricle | |
Hoskins et al. | Patient specific modelling | |
Benitez Mendieta et al. | Imaging-Based Patient-Specific Biomechanical Evaluation of Atherosclerosis and Aneurysm: A Comparison Between Structural-Only, Fluid-Only and Fluid–Structure Interaction Analysis | |
Berthomier et al. | Deep Venous Thrombosis: Database creation and image preprocessing |