JPWO2014030262A1 - 形状データ生成プログラム、形状データ生成方法及び形状データ生成装置 - Google Patents

形状データ生成プログラム、形状データ生成方法及び形状データ生成装置 Download PDF

Info

Publication number
JPWO2014030262A1
JPWO2014030262A1 JP2014531476A JP2014531476A JPWO2014030262A1 JP WO2014030262 A1 JPWO2014030262 A1 JP WO2014030262A1 JP 2014531476 A JP2014531476 A JP 2014531476A JP 2014531476 A JP2014531476 A JP 2014531476A JP WO2014030262 A1 JPWO2014030262 A1 JP WO2014030262A1
Authority
JP
Japan
Prior art keywords
voxel
target object
shape
luminance value
value
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
JP2014531476A
Other languages
English (en)
Other versions
JP5954846B2 (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.)
Fujitsu Ltd
University of Tokyo NUC
Original Assignee
Fujitsu Ltd
University of Tokyo NUC
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 Fujitsu Ltd, University of Tokyo NUC filed Critical Fujitsu Ltd
Application granted granted Critical
Publication of JP5954846B2 publication Critical patent/JP5954846B2/ja
Publication of JPWO2014030262A1 publication Critical patent/JPWO2014030262A1/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • G06T15/08Volume rendering
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/02Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03Computerised tomographs
    • A61B6/032Transmission computed tomography [CT]
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • G06T15/50Lighting effects
    • 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
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/143Segmentation; Edge detection involving probabilistic approaches, e.g. Markov random field [MRF] modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/149Segmentation; Edge detection involving deformable models, e.g. active contour models
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16HHEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/20ICT 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
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B2576/00Medical imaging apparatus involving image processing or analysis
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/46Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with special arrangements for interfacing with the operator or the patient
    • A61B6/461Displaying means of special interest
    • A61B6/466Displaying means of special interest adapted to display 3D data
    • 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/20Special algorithmic details
    • G06T2207/20076Probabilistic image processing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20081Training; Learning
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20092Interactive image processing based on input by user
    • G06T2207/20104Interactive definition of region of interest [ROI]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30048Heart; Cardiac

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Medical Informatics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Biomedical Technology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Pathology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Surgery (AREA)
  • Molecular Biology (AREA)
  • Optics & Photonics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Biophysics (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • Pulmonology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Quality & Reliability (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Epidemiology (AREA)
  • Primary Health Care (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本形状データ生成方法は、複数の断層画像のうち対象物体が占める領域が指定されている断層画像を用いて、対象物体を表す三次元のボクセルデータを生成し、対象物体が占める領域が指定されている断層画像から、指定されている領域に含まれるボクセルの輝度値を抽出し、抽出された当該輝度値を用いて、対象物体が占める領域に含まれるボクセルである確率を算出するための関数を生成し、複数の断層画像を包含するボクセル空間における各ボクセルについて、当該ボクセルの輝度値と関数とを用いて確率を算出し、対象物体の形状を表す三次元のボクセルデータとボクセル空間における各ボクセルについて算出された確率とを用いて、対象物体の形状を表す三次元のボクセルデータを再生成する処理を含む。

Description

形状データ生成プログラム、形状データ生成方法及び形状データ生成装置に関する。
医療分野においては、手術シミュレータ及び臓器シミュレータ等のシミュレータが、治療方針の決定、診断、術後予測、医薬品及び医療機器の開発等に利用されている。このようなシミュレータによるシミュレーションにおいては、臓器の三次元形状データを使用するが、臓器の三次元形状データを生成するのは容易ではないことが多い。臓器は生体内部に存在しているため目視や直接計測を行うことができない上に、臓器の形状は元来複雑であるためである。
従来、CT(Computed Tomography)画像及びMRI(Magnetic Resonance Imaging)画像等の断層画像の輝度値に基づき対象の領域を抽出することで、シミュレーションに利用する三次元形状データを生成する技術(例えばリジョングローイング法)が存在する。例えば、図1に示した心臓の断層画像のデータから、リジョングローイング法により、図2において実線21で囲まれた領域(右室流体の領域)及び実線22で囲まれた領域(左室流体の領域)を抽出し、右室流体及び左室流体の三次元形状データを生成したいとする。ところが、図1に示した心臓の断層画像データに対して実際にリジョングローイング法を実行すると、図3に示すような結果が得られる。図3に示す抽出結果においては、右室流体の一部の領域を、造影剤の影響により抽出できていない。また、左室流体において乳頭筋に該当する領域も抽出できていない。このように、輝度値のみに基づいた領域の抽出には限界がある。
また、輪郭を示す画素であるか否かをその特徴量に基づき評価する評価関数を求め、その評価関数を用いて、領域内の画素が輪郭を示す画素であるかを示す評価値を算出し、評価値に基づき輪郭を決定する技術が存在する。
また、誤差拡散法によって2値化された画像において、細かいエッジのつぶれを最小限にとどめ且つ階調性が高い鮮明な多値画像へ変換する技術が存在する。具体的には、入力された2値画像における注目画素の少なくともエッジ強度を検出し、検出結果に基づき画像の明るさの閾値を設定する。複数の多値化フィルタ内の黒画素及び白画素の割合と閾値とを比較することにより多値化フィルタの大きさを決定する。そして、決定された大きさの多値化フィルタ内の白画素及び黒画素の割合に基づき2値化画像データを多階調画像データに変換する。
しかし、上記のような技術を用いても、対象の領域を正確に抽出できないことがあるため、精度が高い三次元形状データを生成することはできない。
特開2007−307358号公報 特開平8ー191387号公報
A. Tremeau, N. Borel, A Region Growing and Merging Algorithm to Color Segmentation, Pattern recognition, Volume 30, Issue 7, July 1997, Pages 1191-1203
従って、1つの側面では、本発明の目的は、高精度の三次元形状データを生成するための技術を提供することである。
本発明に係る形状データ生成方法は、複数の断層画像のうち対象物体が占める領域が指定されている断層画像を用いて、対象物体を表す三次元のボクセルデータを生成し、対象物体が占める領域が指定されている断層画像から、指定されている領域に含まれるボクセルの輝度値を抽出し、抽出された当該輝度値を用いて、対象物体が占める領域に含まれるボクセルである確率を算出するための関数を生成し、複数の断層画像を包含するボクセル空間における各ボクセルについて、当該ボクセルの輝度値と関数とを用いて確率を算出し、対象物体の形状を表す三次元のボクセルデータとボクセル空間における各ボクセルについて算出された確率とを用いて、対象物体の形状を表す三次元のボクセルデータを再生成する処理を含む。
高精度の三次元形状データを生成できるようになる。
図1は、心臓の断層画像を示す図である。 図2は、心臓の断層画像にシミュレーションに利用する領域を重ねて表示した図である。 図3は、心臓の断層画像に対してリジョングローイング法を実行した結果を示す図である。 図4は、形状データ生成装置の機能ブロック図である。 図5は、形状データ生成部の機能ブロック図である。 図6は、確率マップ生成部の機能ブロック図である。 図7は、DICOMデータに含まれるヘッダ情報の一例を示す図である。 図8は、DICOMデータに含まれる輝度値情報の一例を示す図である。 図9は、ピクセルの値の記述方法について説明するための図である。 図10は、セグメント画像データ格納部に格納されているデータの一例を示す図である。 図11は、メインの処理フローを示す図である。 図12Aは、初期値生成処理の概要について説明するための図である。 図12Bは、初期値生成処理の処理フローを示す図である。 図13は、形状データ生成処理の処理フローを示す図である。 図14は、TPSWarpについて説明するための図である。 図15Aは、ボクセル化処理の処理フローを示す図である。 図15Bは、ボクセル化処理の途中で生成されるデータの一例を示す図である。 図15Cは、ボクセル化処理の途中で生成されるデータの一例を示す図である。 図16Aは、ボクセル化処理の処理フローを示す図である。 図16Bは、ボクセル化処理の途中で生成されるデータの一例を示す図である。 図17は、ボクセル化処理により生成されるデータを示す図である。 図18は、境界特定処理の処理フローを示す図である。 図19は、ぼかし処理により生成されるデータを示す図である。 図20は、確率マップ生成処理の処理フローを示す図である。 図21は、エロージョン処理について説明するための図である。 図22は、輝度値格納部に格納されているデータの一例を示す図である。 図23は、確率マップ生成処理の処理フローを示す図である。 図24は、内側領域についてのヒストグラム及び外側領域についてのヒストグラムを示す図である。 図25は、確率マップ生成処理の結果を特定のz座標について示した図である。 図26は、初期値生成処理の結果を特定のz座標について示した図である。 図27は、領域抽出処理の結果を特定のz座標について示した図である。 図28は、表示装置に表示される画像の一例を示す図である。 図29は、表示装置に表示される画像の一例を示す図である。 図30は、表示装置に表示される画像の一例を示す図である。 図31は、コンピュータの機能ブロック図である。
図4に、本実施の形態における形状データ生成装置の機能ブロック図を示す。形状データ生成装置1は、断層画像データ格納部10と、入力処理部11と、セグメント画像データ格納部12と、形状データ生成部13と、確率マップ生成部14と、マージ部15と、初期値格納部16と、領域抽出部17と、抽出結果格納部18と、出力部19と、標準形状データ格納部20とを含む。
入力処理部11は、断層画像データ格納部10に格納されている断層画像データに対する、シミュレーションに用いる領域の指定をユーザから受け付ける。そして、入力処理部11は、領域指定された断層画像データ(以下、セグメント画像データと呼ぶ)をセグメント画像データ格納部12に格納する。形状データ生成部13は、セグメント画像データ格納部12に格納されているデータ及び標準形状データ格納部20に格納されているデータを用いて、後述する形状データ生成処理を実行し、形状データを生成する。確率マップ生成部14は、セグメント画像データ格納部12に格納されているデータ及び断層画像データ格納部10に格納されているデータを用いて、後述する確率マップ生成処理を実行し、確率マップを生成する。マージ部15は、形状データ生成部13により生成された形状データ及び確率マップ生成部14により生成された確率マップを用いて初期値を生成し、初期値格納部16に格納する。領域抽出部17は、初期値格納部16に格納されている初期値を用いて領域抽出処理を実行し、処理結果を抽出結果格納部18に格納する。出力部19は、抽出結果格納部18に格納されているデータを表示装置(例えばディスプレイ)に表示させる。
図5に、形状データ生成部13の機能ブロック図を示す。形状データ生成部13は、変形処理部131と、第1形状データ格納部132と、ボクセル化処理部133と、ボクセル化データ格納部134と、境界特定部135と、特定結果格納部136と、ぼかし処理部137と、第2形状データ格納部138とを含む。
変形処理部131は、セグメント画像データ格納部12に格納されているデータを用いて、標準形状データ格納部20に格納されている標準形状のデータ(本実施の形態においては、左室流体(乳頭筋部分を含む)の標準形状のデータ)を変形し、変形後の形状データを第1形状データ格納部132に格納する。ボクセル化処理部133は、第1形状データ格納部132に格納されている形状データに対して、後述するボクセル化処理を実行し、処理結果をボクセル化データ格納部134に格納する。境界特定部135は、ボクセル化データ格納部134に格納されているデータを用いて左室流体と左室流体ではない部分との境界に該当するボクセルを特定する処理を実行し、処理結果を特定結果格納部136に格納する。ぼかし処理部137は、ボクセル化データ格納部134に格納されているデータ及び特定結果格納部136に格納されているデータを用いてぼかし処理を実行し、処理結果を第2形状データ格納部138に格納する。
図6に、確率マップ生成部14の機能ブロック図を示す。確率マップ生成部14は、正規化部141と、第1データ格納部142と、モルフォロジ処理部143、第2データ格納部144、輝度値抽出部145及び輝度値格納部149を含むヒストグラム生成部150と、ヒストグラムデータ格納部146と、マッピング部147と、確率マップ格納部148とを含む。
正規化部141は、断層画像データ格納部10に格納されている断層画像データの輝度値を正規化する処理を実行し、処理結果を第1データ格納部142に格納する。モルフォロジ処理部143は、セグメント画像データ格納部12に格納されているセグメント画像データのうち、領域指定された部分に対してエロージョン処理を実行し、処理結果を第2データ格納部144に格納する。輝度値抽出部145は、第1データ格納部142に格納されているデータ、第2データ格納部144に格納されているデータを用いて処理を行い、処理結果を輝度値格納部149に格納する。また、輝度値抽出部145は、輝度値格納部149に格納されているデータを用いてヒストグラムを生成し、ヒストグラムデータ格納部146に格納する。マッピング部147は、ヒストグラムデータ格納部146に格納されているデータを用いて、断層画像データ格納部10に格納されている断層画像データを包含するボクセル空間における各ボクセルに確率データをマッピングする処理を実行し、処理結果を確率マップ格納部148に格納する。
断層画像データ格納部10には、例えば図1に示すような断層画像データが複数枚分格納されている。断層画像データのフォーマットは、例えばDICOM(Digital Imaging and COmmunication in Medicine)である。DICOMのデータは、ヘッダ情報と輝度値の情報とを含む。ヘッダ情報には、例えば図7に示すように、番号毎に患者情報、医療用測定装置の種類及び断層画像の解像度等の情報が含まれる。図7において、「0.456」はピクセル1つ分の距離(ミリ)を表しており、「0.5」は断層画像間の距離(ミリ)を表している。また、輝度値情報には、例えば図8に示すように、R(赤)、G(緑)及び青(B)の各々についての値が含まれる。CTの場合には、CT値(単位はHounsfield Unit)と呼ばれる値が含まれ、最小値は−1000であり、最大値は3095である。図8における各ピクセルの値は、図9に示すように、断層画像の左上隅から右方向に記述され、右端に到達すると、次の行の左端から右方向へ同様に記述される。なお、CT及びMRI等によって取得された複数枚の断層画像を積層すると、三次元の臓器の形状を包含するデータを生成することができる。このデータは、ボリュームデータと呼ばれる。
図10に、セグメント画像データ格納部12に格納されているデータの一例を示す。図10の例は、図1に示した断層画像データに、シミュレーションに用いる部分であるとしてユーザによって指定された領域(本実施の形態においては、左室流体(乳頭筋部分を含む)の部分)に所定のラベル値が付与された画像データである。セグメント画像データ格納部12には、このようなセグメント画像データが複数枚分格納されている。ユーザは、断層画像データ格納部10に格納されている断層画像データの一部(例えば、30枚に1枚)に、手動で領域指定を行っていく。なお、領域指定はセグメンテーションと呼ばれることもある。
標準形状データ格納部20には、左室流体の標準的な形状(以下、標準形状と呼ぶ)のデータが格納されている。具体的には、形状の頂点の情報及び頂点の接続情報を含むSTL(Standard Triangulated Language)データが格納されており、形状は三角形メッシュからなる。但し、データのフォーマットはSTLデータのフォーマットに限られるわけではない。
次に、図11乃至図30を用いて、図4に示した形状データ生成装置1の動作について説明する。
まず、入力処理部11は、断層画像データ格納部10に格納されている断層画像データの一部(例えば、30枚に1枚)を表示装置に順に出力する(図11:ステップS1)。これにより、医師等であるユーザに領域指定を促す。ユーザは、各断層画像データに対し、目的の領域(本実施の形態においては、左室流体である領域)に対し、マウス等を操作してラベル値を付与する。これに応じ、入力処理部11は、領域指定を受け付け(ステップS3)、ユーザによる領域指定により生成されたセグメント画像データをセグメント画像データ格納部12に格納する。
そして、形状データ生成装置1は、初期値生成処理を実行する(ステップS5)。初期値生成処理については、図12A乃至図26を用いて説明する。
まず、図12Aを用いて、初期値生成処理の概要について説明する。初期値生成処理により生成される初期値は、ボクセル空間上のデータであり、後述する反応拡散フィルタの初期値になる。初期値は、形状データと確率マップとをマージすることにより生成される。形状データは、標準形状データをセグメント画像データにより特定される形状に合うように変形することにより得られる。確率マップは、セグメント画像データを用いて作成された輝度値のヒストグラムを用いて、確率データを断層画像データにおける各ボクセルにマッピングすることにより得られる。
以下では、初期値生成処理の処理フローについて説明する。
まず、形状データ生成部13は、形状データ生成処理を実行する(図12B:ステップS11)。形状データ生成処理については、図13乃至図19を用いて説明する。
変形処理部131は、標準形状データ格納部20から標準形状のデータを読み出すと共に、セグメント画像データ格納部12からセグメント画像データを読み出す。そして、変形処理部131は、標準形状のデータ及びセグメント画像データを含むランドマーク設定画面を表示装置に表示させる(図13:ステップS21)。
ユーザは、表示装置に表示されたランドマーク設定画面を見て、標準形状とセグメント画像データにより特定される形状(以下、目標形状と呼ぶ)との大まかな位置合わせを行う。具体的には、(1)標準形状における所定の部分にソースランドマークを設定する。(2)目標形状における、ソースランドマークを配置した位置に対応する部分にターゲットランドマークを設定する。なお、所定の部分とは、例えば左室流体における特徴的な部分である。
変形処理部131は、ソースランドマーク及びターゲットランドマークの設定を受け付け、ソースランドマーク及びターゲットランドマークのデータ(例えば三次元座標)をメインメモリ等の記憶装置に格納する(ステップS23)。
変形処理部131は、後述するTPSWarp(Thin Plate Spline Warp)等の手法により標準形状を変形する処理を実行する(ステップS25)。そして、変形処理部131は、処理結果の形状のデータを第1形状データ格納部132に格納する。なお、処理結果の形状データのフォーマットは、STLデータのフォーマットである。
ここで、TPSWarpについて説明する。図14に示すように、TPSWarpにおいては、ソースランドマーク及び当該ソースランドマークに対応するターゲットランドマークが与えられると、ソースランドマークがターゲットランドマークに重なるように変形を行う。なお、TPSWarpの詳細については、例えば「Fred L. Bookstein "Principal Warps: Thin-Plate Splines and the Decomposition of Deformation",IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE. VOL. 11, NO. 6, PP. 567-585, JUNE 1989」を参照のこと。
このようにすれば、標準形状を、ユーザによって指定された領域に合うように変形できるようになる。
そして、ボクセル化処理部133は、第1形状データ格納部132に格納されている形状データを用いてボクセル化処理を実行する(ステップS27)。ボクセル化処理については、図15A乃至図17を用いて説明する。
はじめに、ボクセル化処理を説明するにあたって前提となる事項を説明する。
をボクセル空間における座標(x,y,z)におけるボクセル値とする。初期値は−1とする。
をボクセル空間におけるボクセルとする。vNは、ボクセル空間に含まれるボクセルの数である。ボクセル空間において、x軸方向におけるボクセルの数をxN、y軸方向におけるボクセルの数をyN、z軸方向におけるボクセルの数をzNとする。従って、xN×yN×zN=vNである。
をボクセル空間の下限(左下)の座標とする。
をボクセル空間の上限(右上)の座標とする。
をボクセルのサイズ(xs×ys×zs)とする。
第1形状データ格納部132に格納されている形状データをSとし、形状データSを以下のように定める。
ここで、{pi=1,...,Nは形状データSにおける頂点を表し、ベクトル{mj=1,...,Mは形状データSにおける三角形メッシュを表す。Nは形状データSにおける頂点の数である。Mは形状データSにおける三角形メッシュの数である。
メッシュmに含まれる頂点を以下のように定める。
なお、ボクセル化処理によって生成されるボクセル空間は、断層画像データ格納部10に格納されている断層画像データによって特定されるボクセル空間と一致するボクセル空間である。すなわち、下限及び上限の座標、各軸方向のボクセル数並びにボクセルサイズが同じであるとする。
以下では、ボクセル化処理の処理フローについて説明する。
まず、ボクセル化処理部133は、三角形メッシュを識別するための変数であるjをj=1と設定する(図15A:ステップS31)。
ボクセル化処理部133は、ボクセル空間におけるx座標を表す変数であるrをr=rminと設定する(ステップS32)。rminは、以下のような値である。
ここで、ガウス記号中の分子の第1項は、xm(j),1((j)はjがmの下付きの文字であることを表す)の最小値を求める関数である。
ボクセル化処理部133は、ボクセル空間におけるy座標を表す変数であるsをs=sminと設定する(ステップS33)。sminは、以下のような値である。
ボクセル化処理部133は、ボクセル空間におけるz座標を表す変数であるtをt=tminと設定する(ステップS34)。tminは、以下のような値である。
ボクセル化処理部133は、(r,s,t)が形状データSの内部に有るか判断する(ステップS35)。ステップS35においては、以下の式を満たすか判断する。
ここで、
とする。
数式11において、((vs・v+L)−pm(j),1)((j)はjがmの下付きの文字であることを表す)は、始点をpm(j),1とし終点をボクセル(r,s,t)とするベクトルを表す。((pm(j),2−pm(j),1)×(pm(j),3−pm(j),1))は、三角形メッシュを含む平面を形状の内側から外側の方向へ貫く法線ベクトルを表す。前者のベクトルと法線ベクトルとの内積が0より小さいか否かを判断することにより、ボクセル(r,s,t)が形状データSの内部に有るか否かを判断することができる。但し、ステップS35の処理における判断の対象となるボクセルは、左室流体と左室流体でない部分との境界付近のボクセルに限られている。
(r,s,t)が形状データSの内部に有る場合(ステップS35:Yesルート)、ボクセル化処理部133は、V(r,s,t)=1と設定する(ステップS36)。一方、(r,s,t)が形状データSの内部に無い場合(ステップS35:Noルート)、ボクセル化処理部133は、V(r,s,t)=0と設定する(ステップS37)。
ボクセル化処理部133は、t≧tmaxであるか判断する(ステップS38)。tmaxは、以下のような値である。
ここで、ガウス記号中の分子の第1項は、zm(j),1((j)はjがmの下付きの文字であることを表す)の最大値を求める関数である。
t≧tmaxではない場合(ステップS38:Noルート)、ボクセル化処理部133は、tを1インクリメントし(ステップS39)、ステップS35の処理に戻る。一方、t≧tmaxである場合(ステップS38:Yesルート)、ボクセル化処理部133は、s≧smaxであるか判断する(ステップS40)。smaxは、以下のような値である。
s≧smaxではない場合(ステップS40:Noルート)、ボクセル化処理部133は、sを1インクリメントし(ステップS41)、ステップS34の処理に戻る。一方、s≧smaxである場合(ステップS40:Yesルート)、ボクセル化処理部133は、r≧rmaxであるか判断する(ステップS42)。rmaxは、以下のような値である。
r≧rmaxではない場合(ステップS42:Noルート)、ボクセル化処理部133は、rを1インクリメントし(ステップS43)、ステップS33の処理に戻る。一方、r≧rmaxである場合(ステップS42:Yesルート)、ボクセル化処理部133は、j≧Mであるか判断する(ステップS44)。j≧Mではない場合(ステップS44:Noルート)、ボクセル化処理部133は、jを1インクリメントし(ステップS45)、ステップS32の処理に戻る。一方、j≧Mである場合(ステップS44:Yesルート)、処理は端子Aを介して図16AのステップS51に移行する。
ここまでの処理によって生成されるボクセルデータの一例を図15Bに示す。図15Bは、特定のz座標におけるボクセル値を示している。図15Bにおいては、ボクセル値が1であるボクセルを白で示し、ボクセル値が0であるボクセルを灰色で示し、ボクセル値が−1であるボクセルを黒で示している。ここまでの処理では、左室流体と左室流体ではない部分との境界付近のボクセル以外のボクセルについては、初期値である−1が設定されたままである。
図15Cは、図15Bにおいて白線150で囲まれた領域を拡大した図である。図15Cに示すように、境界付近のボクセルのうち形状データSの内部にあるボクセルのボクセル値は1であり、境界付近のボクセルのうち形状データSの外部にあるボクセルのボクセル値は0である。
図16Aの説明に移行し、ボクセル化処理部133は、ボクセル空間に対して、シード点を(0,0,0)としてリジョングローイング法を実行し(図16A:ステップS51)、処理結果をメインメモリ等の記憶装置に格納する。リジョングローイング法は、シード点に連結する点を取り込むことで領域を拡張し、対象領域を抽出する方法である。リジョングローイング法については、例えば「A. Tremeau, N. Borel, A Region Growing and Merging Algorithm to Color Segmentation, Pattern recognition, Volume 30, Issue 7, July 1997, Pages 1191-1203」を参照のこと。なお、ステップS51においてはシード点を(0,0,0)としているが、左室流体ではない部分(すなわち、形状データSの外側の部分)のボクセルであればシード点として利用することができる。
ボクセル化処理部133は、x座標を識別するための変数xをx=1と設定する(ステップS52)。
ボクセル化処理部133は、y座標を識別するための変数yをy=1と設定する(ステップS53)。
ボクセル化処理部133は、z座標を識別するための変数zをz=1と設定する(ステップS54)。
ボクセル化処理部133は、ボクセル(x,y,z)がリジョングローイング法によって抽出されたか判断する(ステップS55)。抽出された場合(ステップS55:Yesルート)、ボクセル(x,y,z)は形状データSの外側に存在するので、ボクセル化処理部133は、V(x,y,z)=0と設定する(ステップS56)。抽出されていない場合(ステップS55:Noルート)、ボクセル化処理部133は、V(x,y,z)が−1であるか判断する(ステップS57)。
V(x,y,z)が−1である場合(ステップS57:Yesルート)、ボクセル(x,y,z)は形状データSの内側にあり且つ値が初期値のままなので、ボクセル化処理部133は、V(x,y,z)=1と設定する(ステップS58)。一方、V(x,y,z)が−1ではない場合(ステップS57:Noルート)、V(x,y,z)=0又は1と既に設定されているので、ステップS59の処理に移行する。
ボクセル化処理部133は、z≧zNであるか判断する(ステップS59)。z≧zNではない場合(ステップS59:Noルート)、ボクセル化処理部133は、zを1インクリメントし(ステップS60)、ステップS55の処理に戻る。
z≧zNである場合(ステップS59:Yesルート)、ボクセル化処理部133は、y≧yNであるか判断する(ステップS61)。y≧yNではない場合(ステップS61:Noルート)、ボクセル化処理部133は、yを1インクリメントし(ステップS62)、ステップS54の処理に戻る。
y≧yNである場合(ステップS61:Yesルート)、ボクセル化処理部133は、x≧xNであるか判断する(ステップS63)。x≧xNではない場合(ステップS63:Noルート)、ボクセル化処理部133は、xを1インクリメントし(ステップS64)、ステップS53の処理に戻る。一方、x≧xNである場合(ステップS63:Yesルート)、元の処理に戻る。
図16Bを用いて、ステップS56の処理について説明する。本実施の形態においては、シード点を(0,0,0)としている(すなわち、形状データSの外側の部分としている)ため、リジョングローイング法においては、形状データSの外側の部分が抽出される。よって、ステップS56において、リジョングローイング法によって抽出されたボクセルの値を0とすると、図16Bのようになる。図16Bは、図15Bと同様、特定のz座標におけるボクセル値を示しており、ボクセル値が1であるボクセルを白で示し、ボクセル値が0であるボクセルを灰色で示し、ボクセル値が−1であるボクセルを黒で示している。但し、ステップS37の処理においてV(x,y,z)=0と設定されたボクセルについては、灰色の上に黒の縞を付している。そして、図16Bにおいて形状データSの内側の黒色部分に該当するボクセルの値は、ステップS58の処理により1に置換される。従って、形状データSの内側のボクセル値が1となり、形状データSの外側のボクセル値が0となる。
以上のようなボクセル化処理を実行することにより、STLデータであった形状データを、ボクセル空間上のデータとする。ボクセル化処理を実行すると、例えば図17に示すようなデータを生成することができる。図17においては、ボクセル値が1である部分を明るい色(白)で示しており、ボクセル値が0である部分を黒で示している。この段階では、左室流体と左室流体ではない部分との境界のボクセル値は0又は1であるため、境界は滑らかではなく凸凹がある。
図13の説明に戻り、ボクセル化処理(ステップS27)が終了すると、境界特定部135は、境界特定処理を実行する(ステップS28)。境界特定処理については、図18を用いて説明する。
まず、境界特定部135は、ボクセルを識別するための変数iをi=1と設定する(図18:ステップS71)。
境界特定部135は、vは左室流体と左室流体ではない部分との境界に該当するか判断する(ステップS73)。ステップS73においては、以下の式を満たすか判断する。
が境界に該当する場合(ステップS73:Yesルート)、境界特定部135は、v=(x,y,z)を境界b={v}に追加し(ステップS75)、v=(x,y,z)の座標を特定結果格納部136に格納する。vが境界に該当しない場合(ステップS73:Noルート)、ステップS77の処理に移行する。
境界特定部135は、i≧vNであるか判断する(ステップS77)。i≧vNではない場合(ステップS77:Noルート)、境界特定部135は、iを1インクリメントし(ステップS79)、ステップS73の処理に戻る。一方、i≧vNである場合(ステップS77:Yesルート)、元の処理に戻る。
このように、境界に該当するボクセルを特定しておけば、ぼかし処理に利用できるようになる。
図13の説明に戻り、ぼかし処理部137は、ボクセル化データ格納部134に格納されているデータにより特定される形状の境界に対してぼかし処理を実行し、処理結果を第2形状データ格納部138に格納する(ステップS29)。そして元の処理に戻る。
ステップS29においては、特定結果格納部136にデータが格納されているボクセルのボクセル値の平均が0.5(すなわち、0と1との平均)になるように、ガウシアンフィルタを実行する。ガウシアンフィルタを実行すると、ボクセル値の分布は図19に示すようになる。図19において、ボクセル値は0から1までの間の値を取るようになっており、ボクセル値が大きいほど明るい色で示している。境界においては、ボクセル値は1から0までに徐々に変化するようになっている。
図12Bの説明に戻り、形状データ生成処理が終了すると、確率マップ生成部14は、確率マップ生成処理を実行する(ステップS13)。確率マップ生成処理については、図20乃至図25を用いて説明する。
まず、正規化部141は、断層画像データ格納部10に格納されている断層画像データの輝度値を正規化する処理を実行し、処理結果を第1データ格納部142に格納する(図20:ステップS81)。ステップS81においては、例えば、輝度値を−1000から1000の間の値に正規化する。但し、このような値の範囲に限られるわけではない。
ヒストグラム生成部150におけるモルフォロジ処理部143は、セグメント画像データ格納部12から未処理のセグメント画像データを1枚分特定する(ステップS83)。
モルフォロジ処理部143は、特定されたセグメント画像データに対して、モルフォロジ処理の一種であるエロージョン処理を実行し、処理結果を第2データ格納部144に格納する(ステップS85)。なお、モルフォロジ処理については、例えば「L. Vincent, Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms, IEEE Transactions on Image Processing, 2(2):176-201, 1993.」を参照のこと。
図21を用いて、ステップS85におけるエロージョン処理について説明する。図21において、曲線202で囲まれた領域(以下、内側領域と呼ぶ)はユーザにより指定された領域である。ステップS85におけるエロージョン処理においては、ユーザにより指定された領域を、曲線201で囲まれた領域によりマスクする。以下では、曲線201と曲線202との間の領域を外側領域と呼ぶ。
輝度値抽出部145は、ステップS83において特定されたセグメント画像データの内側領域に対応する領域を、第1データ格納部142に格納されている当該セグメント画像データに対応する断層画像データから特定する(ステップS87)。
輝度値抽出部145は、ステップS87において特定された領域に含まれるボクセルの輝度値を抽出し、輝度値格納部149に格納する(ステップS89)。
図22に、輝度値格納部149に格納されているデータの一例を示す。図22の例では、内側領域及び外側領域について、断層画像データにおけるボクセルの座標と、そのボクセルの輝度値とが格納されている。
輝度値抽出部145は、ステップS83において特定されたセグメント画像データの外側領域に対応する領域を、第1データ格納部142に格納されている当該セグメント画像データに対応する断層画像データから特定する(ステップS91)。
輝度値抽出部145は、ステップS91において特定された領域に含まれるボクセルの輝度値を抽出し、輝度値格納部149に格納する(ステップS93)。
輝度値抽出部145は、未処理のセグメント画像データがセグメント画像データ格納部12に有るか判断する(ステップS95)。未処理のセグメント画像データが有る場合(ステップS95:Yesルート)、次のセグメント画像データについて処理するため、ステップS83の処理に戻る。一方、未処理のセグメント画像データが無い場合(ステップS95:Noルート)、処理は端子Bを介して図23のステップS97に移行する。
図23の説明に移行し、輝度値抽出部145は、輝度値格納部149に格納されている輝度値を用いて、内側領域についてのヒストグラムHin(Ii,j,k)及び外側領域についてのヒストグラムHout(Ii,j,k)を生成する(ステップS97)。そして、輝度値抽出部145は、生成されたヒストグラムをヒストグラムデータ格納部146に格納する。
図24に、ステップS97において生成されるヒストグラムを示す。図24の縦軸は頻度であり、横軸は輝度値である。内側領域についてのヒストグラムHin(Ii,j,k)は、輝度値が300付近にピークがある。一方、外側領域についてのヒストグラムHout(Ii,j,k)は、輝度値が70付近にピークがある。
マッピング部147は、ヒストグラムデータ格納部146に格納されているヒストグラムを用いて関数を生成する(ステップS98)。また、マッピング部147は、断層画像データにおける各ボクセルに対し、生成された関数を用いて確率データをマッピングする(ステップS99)。そして、マッピング部147は、マッピングの結果を確率マップ格納部148に格納する。そして元の処理に戻る。ステップS99においては、以下の関数を用いてマッピングをする。
数式17によって求められる確率は、断層画像データにおけるボクセル(i,j,k)が左室流体に含まれる確率を表している。確率マップ生成処理によって生成されたデータの一例を図25に示す。図25は、特定のz座標における確率データを示している。図25においては、確率が1に近いほど明るい色になり、確率が0に近いほど暗い色になっている。なお、本実施の形態においては左室流体についての確率マップを生成しているため、図25において左室流体以外の領域についての確率は意味を持たない。
図12Bの説明に戻り、確率マップ生成処理が終了すると、マージ部15は、第2形状データ格納部138に格納されているデータ及び確率マップ格納部148に格納されている確率データをマージすることにより初期値を生成する(ステップS15)。そして、マージ部15は、生成された初期値を初期値格納部16に格納し、元の処理に戻る。ステップS15においては、以下の式に従って初期値を生成する。
ここで、λは0<λ<1であり、例えば0.5とする。D1i,j,kは、第2形状データ格納部138に格納されている、ボクセル(i,j,k)におけるボクセル値である。D2i,j,kは、確率マップ格納部148に格納されている、ボクセル(i,j,k)の確率データである。epsは十分に小さな値であり、例えば0.01とする。
図26に、ステップS15において生成された初期値の一例を示す。図26は、特定のz座標における初期値を示している。図26においては、値が1に近いほど明るい色になり、値が0に近いほど暗い色になっている。また、数式18においてepsの条件を利用しているため、左室流体以外であるとみなされる領域は除去されている(すなわち、値が0になっている)。
図11の説明に戻り、領域抽出部17は、初期値格納部16に格納されている初期値を用いて領域抽出処理を実行し(ステップS7)、処理結果を抽出結果格納部18に格納する。
ステップS7においては、以下の反応拡散方程式を利用して、初期値格納部16に格納されている初期値に対してフィルタをかけることにより、左室流体の領域を抽出している。
ここで、uはボクセル値である。αは反応係数であり、βは拡散係数である。
第1項は反応項であり、以下の関数を用いる。
反応項は、uを代表値に集める効果をもたらす。上記の反応拡散方程式の場合、uが0から0.5までの値である場合には、uは次第に0に近付く。uが0.5から1までの値である場合には、uは次第に1に近付く。
第2項は拡散項であり、ガウシアンフィルタと同様の効果(すなわち、平滑化の効果)をもたらす。
数式19を離散化すると、以下のようになる。
ここで、ui,j,kは、ボクセル(i,j,k)におけるボクセル値である。Dは以下のようになる。
また、外縁のボクセルの値uは、以下のようにする。
セグメント画像データにおいて指定された領域に含まれるボクセルの値uは、ノイマン境界条件として、以下のようにする。
セグメント画像データにおいて指定された領域に含まれないボクセルの値uは、ノイマン境界条件として、以下のようにする。
心室と心房の境界に該当するボクセルの値uは、ディリクレ境界条件として、以下のようにする。
以上のような条件のもとで領域抽出処理を実行する(すなわち、所定回数計算を実行する)と、例えば図27に示すようなデータが得られる。図27は、特定のz座標におけるボクセル値を示している。領域抽出処理により、ボクセル値は0又は1に近付くようになっている。そのため、図26に示した初期値と比較すると、境界付近に存在していたぼかされた部分は少なくなっている。
図11の説明に戻り、出力部19は、抽出結果格納部18に格納されているデータを用いて、抽出した領域を示すデータを生成し、表示装置に表示させる(ステップS9)。そして処理を終了する。
図28乃至図30は、表示装置に表示される画像の一例を示す図である。図28においては、図27と同じz座標における断層画像データと、図27に示した領域(白い部分)を示す太線280とが重ねて表示されている。図29においては、図28とは異なるz座標における断層画像データと、領域抽出処理によって抽出された領域を示す太線290とが重ねて表示されている。また、図30においては、図28及び図29とは異なるz座標における断層画像データと、領域抽出処理によって抽出された領域を示す太線300とが重ねて表示されている。
図28乃至図30に示すように、領域抽出処理によって抽出された領域には、乳頭筋部分が含まれるようになっている。また、造影剤の影響によるムラが無くなっている。すなわち、輝度値のみによって左室流体に該当する領域を抽出した場合と比較して、シミュレーションに用いる領域を適切に抽出することができるようになる。
以上本技術の一実施の形態を説明したが、本技術はこれに限定されるものではない。例えば、上で説明した形状データ生成装置1の機能ブロック図は必ずしも実際のプログラムモジュール構成に対応するものではない。
また、上で説明した処理フローは、処理結果が変わらなければ処理の順番を入れ替えることも可能である。さらに、並列に実行させるようにしても良い。
例えば、初期値生成処理においては形状データ生成処理を実行した後に確率マップ生成処理を実行しているが、これらの処理を並行で実行してもよい。
また、上では左室流体を例にして説明したが、その他の複雑な形状を有する臓器であっても本実施の形態の処理を適用することができる。
なお、上で述べた形状データ生成装置1は、コンピュータ装置であって、図31に示すように、メモリ2501とCPU(Central Processing Unit)2503とハードディスク・ドライブ(HDD:Hard Disk Drive)2505と表示装置2509に接続される表示制御部2507とリムーバブル・ディスク2511用のドライブ装置2513と入力装置2515とネットワークに接続するための通信制御部2517とがバス2519で接続されている。オペレーティング・システム(OS:Operating System)及び本実施例における処理を実施するためのアプリケーション・プログラムは、HDD2505に格納されており、CPU2503により実行される際にはHDD2505からメモリ2501に読み出される。CPU2503は、アプリケーション・プログラムの処理内容に応じて表示制御部2507、通信制御部2517、ドライブ装置2513を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ2501に格納されるが、HDD2505に格納されるようにしてもよい。本発明の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク2511に格納されて頒布され、ドライブ装置2513からHDD2505にインストールされる。インターネットなどのネットワーク及び通信制御部2517を経由して、HDD2505にインストールされる場合もある。このようなコンピュータ装置は、上で述べたCPU2503、メモリ2501などのハードウエアとOS及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
以上述べた本実施の形態をまとめると以下のようになる。
本実施の形態に係る形状データ生成方法は、(A)複数の断層画像のうち対象物体が占める領域が指定されている断層画像を用いて、対象物体を表す三次元のボクセルデータを生成し、(B)対象物体が占める領域が指定されている断層画像から、指定されている領域に含まれるボクセルの輝度値を抽出し、抽出された当該輝度値を用いて、対象物体が占める領域に含まれるボクセルである確率を算出するための関数を生成し、(C)複数の断層画像を包含するボクセル空間における各ボクセルについて、当該ボクセルの輝度値と関数とを用いて確率を算出し、(D)対象物体の形状を表す三次元のボクセルデータとボクセル空間における各ボクセルについて算出された確率とを用いて、対象物体の形状を表す三次元のボクセルデータを再生成する処理を含む。
上で述べたように、対象物体が占める領域に含まれるボクセルである確率を用いて三次元のボクセルデータを再生成すれば、高精度の三次元形状になる。
また、上で述べた輝度値を抽出する処理において、(b1)指定されている領域以外の領域に含まれるボクセルの輝度値を抽出し、関数を生成する処理において、(b2)指定されている領域に含まれるボクセルの輝度値についての第1の度数分布と、指定されている領域以外の領域に含まれるボクセルの輝度値についての第2の度数分布とを生成し、第1の度数分布の値を第1の度数分布の値と第2の度数分布の値との和で除した値を確率とする関数を生成してもよい。
このような関数とすれば、指定されている領域に含まれるボクセルの輝度値と同程度の輝度値を有するボクセルの確率が高くなる。
また、上で述べた輝度値を抽出する処理において、(b3)対象物体が占める領域が指定されている断層画像のうち指定されている領域に対しエロージョン処理を実行し、当該エロージョン処理により特定される領域のうち指定されている領域以外の領域から輝度値を抽出してもよい。
エロージョン処理により特定される領域は指定されている領域に肉付けされたような領域になるので、このようにすれば、指定されている領域以外の領域の輝度値を適切に抽出できるようになる。
また、上で述べた対象物体の形状を表す三次元のボクセルデータを再生成する処理において、(d1)対象物体の形状を表す三次元のボクセルデータに含まれるボクセルの輝度値と、当該ボクセルに対応するボクセル空間上のボクセルの確率との加重平均を算出してもよい。
このようにすれば、算出された確率が適切に反映された三次元のボクセルデータを生成できるようになる。
また、本形状データ生成方法が、(E)再生成された対象物体の形状を表す三次元のボクセルデータに含まれるボクセルの輝度値を、当該輝度値についての反応拡散方程式を解くことにより、対象物体の内部であることを示す第1の輝度値又は対象物体の外部であることを示す第2の輝度値のいずれかに近付けてもよい。
このようにすれば、対象物体と対象物体ではない部分との境界がぼやけている場合であっても、境界を明瞭にすることができるようになる。
また、上で述べた対象物体の形状を表す三次元のボクセルデータを生成する処理において、(a1)対象物体の標準的な形状を表す三角形メッシュデータを、指定されている領域に合うように変形し、(a2)変形後の三角形メッシュデータにより特定される形状の内部が第1の輝度値になり且つ外部が第2の輝度値となるように、対象物体の形状を表す三次元のボクセルデータを生成してもよい。
このように標準的な形状を利用すれば、最終的に生成される三次元のボクセルデータが現実的なものになる。
また、上で述べた対象物体の形状を表す三次元のボクセルデータを生成する処理において、(a3)対象物体の形状を表す三次元のボクセルデータに含まれるボクセルのうち、対象物体と対象物体以外の部分との境界に該当するボクセルを特定し、(a4)特定された境界に該当するボクセルの輝度値が第1の輝度値と第2の輝度値との平均値になるようにガウシアンフィルタを実行してもよい。
ボクセル値が第1の輝度値又は第2の輝度値のいずれかであるならば、その平均値が境界に該当すると考えるのが妥当であるからである。
また、上で述べた対象物体の形状を表す三次元のボクセルデータを生成する処理において、(a5)対象物体の形状を表す三次元のボクセルデータに含まれるボクセルのうち対象物体以外の部分のボクセルをシード点としてリジョングローイング法による処理を実行してもよい。
このようにすれば、対象物体以外の部分を抽出することができるので、対象物体の内部と外部とを分けることができるようになる。
また、上で述べた対象物体は心臓であってもよい。心臓のような複雑な臓器であっても、上で述べた処理であれば高精度の三次元形状になる。
なお、上記方法による処理をコンピュータに行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、CD−ROM、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。
まず、入力処理部11は、断層画像データ格納部10に格納されている断層画像データの一部(例えば、30枚に1枚)を表示装置に順に出力する(図11:ステップS1)。これにより、医師等であるユーザに領域指定を促す。ユーザは、各断層画像データにおける目的の領域(本実施の形態においては、左室流体である領域)に対し、マウス等を操作してラベル値を付与する。これに応じ、入力処理部11は、領域指定を受け付け(ステップS3)、ユーザによる領域指定により生成されたセグメント画像データをセグメント画像データ格納部12に格納する。
ここで、TPSWarpについて説明する。図14に示すように、TPSWarpにおいては、ソースランドマーク及び当該ソースランドマークに対応するターゲットランドマークが与えられると、ソースランドマークがターゲットランドマークに重なるように変形を行う。なお、TPSWarpの詳細については、例えば「Fred L. Bookstein, "Principal Warps: Thin-Plate Splines and the Decomposition of Deformations",IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 11, NO. 6, PP. 567-585, JUNE 1989」を参照のこと。
図15Cは、図15Bにおいて白線15で囲まれた領域を拡大した図である。図15Cに示すように、境界付近のボクセルのうち形状データSの内部にあるボクセルのボクセル値は1であり、境界付近のボクセルのうち形状データSの外部にあるボクセルのボクセル値は0である。

Claims (11)

  1. コンピュータに、
    複数の断層画像のうち対象物体が占める領域が指定されている断層画像を用いて、前記対象物体を表す三次元のボクセルデータを生成させ、
    前記対象物体が占める領域が指定されている断層画像から、指定されている領域に含まれるボクセルの輝度値を抽出させ、抽出された当該輝度値を用いて、前記対象物体が占める領域に含まれるボクセルである確率を算出するための関数を生成させ、
    前記複数の断層画像を包含するボクセル空間における各ボクセルについて、当該ボクセルの輝度値と前記関数とを用いて前記確率を算出させ、
    前記対象物体の形状を表す三次元のボクセルデータと前記ボクセル空間における各ボクセルについて算出された確率とを用いて、前記対象物体の形状を表す三次元のボクセルデータを再生成させる
    ことを特徴とする形状データ生成プログラム。
  2. 前記輝度値を抽出させる処理において、
    前記指定されている領域以外の領域に含まれるボクセルの輝度値を抽出させ、
    前記関数を生成させる処理において、
    前記指定されている領域に含まれるボクセルの輝度値についての第1の度数分布と、前記指定されている領域以外の領域に含まれるボクセルの輝度値についての第2の度数分布とを生成させ、前記第1の度数分布の値を前記第1の度数分布の値と前記第2の度数分布の値との和で除した値を前記確率とする関数を生成させる
    ことを特徴とする請求項1記載の形状データ生成プログラム。
  3. 前記輝度値を抽出させる処理において、
    前記対象物体が占める領域が指定されている断層画像のうち前記指定されている領域に対しエロージョン処理を実行させ、当該エロージョン処理により特定される領域のうち前記指定されている領域以外の領域から前記輝度値を抽出させる
    ことを特徴とする請求項2記載の形状データ生成プログラム。
  4. 前記対象物体の形状を表す三次元のボクセルデータを再生成させる処理において、
    前記対象物体の形状を表す三次元のボクセルデータに含まれるボクセルの輝度値と、当該ボクセルに対応する前記ボクセル空間上のボクセルの確率との加重平均を算出させる
    ことを特徴とする請求項1乃至3のいずれか1つ記載の形状データ生成プログラム。
  5. 前記コンピュータに、さらに、
    再生成された前記対象物体の形状を表す三次元のボクセルデータに含まれるボクセルの輝度値を、当該輝度値についての反応拡散方程式を解かせることにより、前記対象物体の内部であることを示す第1の輝度値又は前記対象物体の外部であることを示す第2の輝度値のいずれかに近付けさせる
    ことを特徴とする請求項1乃至4のいずれか1つ記載の形状データ生成プログラム。
  6. 前記対象物体の形状を表す三次元のボクセルデータを生成させる処理において、
    前記対象物体の標準的な形状を表す三角形メッシュデータを、前記指定されている領域に合うように変形させ、
    変形後の前記三角形メッシュデータにより特定される形状の内部が前記第1の輝度値になり且つ外部が前記第2の輝度値となるように、前記対象物体の形状を表す三次元のボクセルデータを生成させる
    ことを特徴とする請求項5記載の形状データ生成プログラム。
  7. 前記対象物体の形状を表す三次元のボクセルデータを生成させる処理において、
    前記対象物体の形状を表す三次元のボクセルデータに含まれるボクセルのうち、前記対象物体と前記対象物体以外の部分との境界に該当するボクセルを特定させ、
    特定された前記境界に該当するボクセルの輝度値が前記第1の輝度値と前記第2の輝度値との平均値になるようにガウシアンフィルタを実行させる
    ことを特徴とする請求項6記載の形状データ生成プログラム。
  8. 前記対象物体の形状を表す三次元のボクセルデータを生成させる処理において、
    前記対象物体の形状を表す三次元のボクセルデータに含まれるボクセルのうち前記対象物体以外の部分のボクセルをシード点としてリジョングローイング法による処理を実行させる
    ことを特徴とする請求項6記載の形状データ生成プログラム。
  9. 前記対象物体は心臓である
    ことを特徴とする請求項1乃至8のいずれか1つ記載の形状データ生成プログラム。
  10. コンピュータが、
    複数の断層画像のうち対象物体が占める領域が指定されている断層画像を用いて、前記対象物体を表す三次元のボクセルデータを生成し、
    前記対象物体が占める領域が指定されている断層画像から、指定されている領域に含まれるボクセルの輝度値を抽出し、抽出された当該輝度値を用いて、前記対象物体が占める領域に含まれるボクセルである確率を算出するための関数を生成し、
    前記複数の断層画像を包含するボクセル空間における各ボクセルについて、当該ボクセルの輝度値と前記関数とを用いて前記確率を算出し、
    前記対象物体の形状を表す三次元のボクセルデータと前記ボクセル空間における各ボクセルについて算出された確率とを用いて、前記対象物体の形状を表す三次元のボクセルデータを再生成する
    処理を実行する形状データ生成方法。
  11. 複数の断層画像のうち対象物体が占める領域が指定されている断層画像を用いて、前記対象物体を表す三次元のボクセルデータを生成する第1生成部と、
    前記対象物体が占める領域が指定されている断層画像から、指定されている領域に含まれるボクセルの輝度値を抽出し、抽出された当該輝度値を用いて、前記対象物体が占める領域に含まれるボクセルである確率を算出するための関数を生成し、前記複数の断層画像を包含するボクセル空間における各ボクセルについて、当該ボクセルの輝度値と前記関数とを用いて前記確率を算出する算出部と、
    前記対象物体の形状を表す三次元のボクセルデータと前記ボクセル空間における各ボクセルについて算出された確率とを用いて、前記対象物体の形状を表す三次元のボクセルデータを再生成する第2生成部と、
    を有する形状データ生成装置。
JP2014531476A 2012-08-24 2012-08-24 形状データ生成プログラム、形状データ生成方法及び形状データ生成装置 Active JP5954846B2 (ja)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2012/071502 WO2014030262A1 (ja) 2012-08-24 2012-08-24 形状データ生成プログラム、形状データ生成方法及び形状データ生成装置

Publications (2)

Publication Number Publication Date
JP5954846B2 JP5954846B2 (ja) 2016-07-20
JPWO2014030262A1 true JPWO2014030262A1 (ja) 2016-07-28

Family

ID=50149595

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2014531476A Active JP5954846B2 (ja) 2012-08-24 2012-08-24 形状データ生成プログラム、形状データ生成方法及び形状データ生成装置

Country Status (4)

Country Link
US (1) US9390549B2 (ja)
EP (1) EP2889001B1 (ja)
JP (1) JP5954846B2 (ja)
WO (1) WO2014030262A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TWI594206B (zh) 2016-01-08 2017-08-01 Nat Yang-Ming Univ Cardiac medical imaging single chamber mapping system and method
US10453200B2 (en) * 2016-11-02 2019-10-22 General Electric Company Automated segmentation using deep learned priors

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007144056A (ja) * 2005-11-30 2007-06-14 Konica Minolta Medical & Graphic Inc モデリング装置、領域抽出装置、モデリング方法及びプログラム
WO2011106440A1 (en) * 2010-02-23 2011-09-01 Loma Linda University Medical Center Method of analyzing a medical image

Family Cites Families (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08191387A (ja) 1995-01-12 1996-07-23 Fuji Xerox Co Ltd 画像処理装置
US6008813A (en) * 1997-08-01 1999-12-28 Mitsubishi Electric Information Technology Center America, Inc. (Ita) Real-time PC based volume rendering system
US6396898B1 (en) * 1999-12-24 2002-05-28 Kabushiki Kaisha Toshiba Radiation detector and x-ray CT apparatus
WO2002073536A2 (en) 2001-03-09 2002-09-19 Koninklijke Philips Electronics N.V. Image segmentation
DE10111661A1 (de) 2001-03-09 2002-09-12 Philips Corp Intellectual Pty Verfahren zum Segmentieren einer in einem Objekt enthaltenen dreidimensionalen Struktur, insbesondere für die medizinische Bildanalyse
JP3932360B2 (ja) 2003-03-04 2007-06-20 独立行政法人産業技術総合研究所 ランドマーク抽出装置およびランドマーク抽出方法
IL164030A0 (en) * 2003-09-12 2005-12-18 Revital Pery Shechter Photoacoustic analyzer of a region of interest in a human body
JP2007098028A (ja) 2005-10-07 2007-04-19 Konica Minolta Medical & Graphic Inc モデリング装置、モデリング方法、領域抽出装置、およびプログラム
JP4720478B2 (ja) 2005-12-15 2011-07-13 コニカミノルタエムジー株式会社 モデリング装置、領域抽出装置およびプログラム
JP4999163B2 (ja) 2006-04-17 2012-08-15 富士フイルム株式会社 画像処理方法および装置ならびにプログラム
JP2008009549A (ja) 2006-06-27 2008-01-17 Yamaguchi Univ 画像処理方法、画像処理装置及び画像処理プログラム
JP4191753B2 (ja) * 2006-07-12 2008-12-03 ザイオソフト株式会社 画像処理方法
JP4260177B2 (ja) * 2006-09-27 2009-04-30 ザイオソフト株式会社 画像処理方法
US8098918B2 (en) * 2007-09-21 2012-01-17 Siemens Corporation Method and system for measuring left ventricle volume
JP5610129B2 (ja) * 2010-03-26 2014-10-22 富士通株式会社 3次元テンプレート変形方法、装置及びプログラム
JP5709216B2 (ja) 2011-05-26 2015-04-30 富士通株式会社 画像処理プログラム、方法及び装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007144056A (ja) * 2005-11-30 2007-06-14 Konica Minolta Medical & Graphic Inc モデリング装置、領域抽出装置、モデリング方法及びプログラム
WO2011106440A1 (en) * 2010-02-23 2011-09-01 Loma Linda University Medical Center Method of analyzing a medical image

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
JPN6016019979; 佐竹 功次 koji SATAKE: '領域拡張法と確率アトラスを用いた3次元非造影腹部CT画像からの肝臓領域抽出法 Liver Segmentation Met' 電子情報通信学会技術研究報告 Vol.108 No.47 IEICE Technical Report 第108巻, 社団法人電子情報通信学会 The Institute of Electro *
JPN6016019982; 坂本 佳隆 Yoshitaka SAKAMOTO: '局所濃度値特徴と多クラス識別器を用いた腹部3次元CT像からの肝血管腫検出 Detection of the liver hem' 電子情報通信学会技術研究報告 Vol.111 No.389 IEICE Technical Report 第111巻, 社団法人電子情報通信学会 The Institute of Electro *
JPN6016019985; 北川 輝彦 Teruhiko KITAGAWA: '体幹部非造影X線CT像における肝臓アトラスの構築とその肝臓自動抽出法への応用 Generation of a Probab' 電子情報通信学会論文誌 (J91-D) 第7号 THE IEICE TRANSACTIONS ON INFORMATION AND SYSTEMS (J 第J91-D巻, 社団法人電子情報通信学会 THE INSTITUTE OF ELECTRO *
JPN6016019987; 木本 達也 Tatsuya KIMOTO: '被検者固有アトラスを用いた3次元造影CT像からの膵臓領域抽出 Pancreas segmentation from three dimen' 電子情報通信学会技術研究報告 Vol.108 No.385 IEICE Technical Report 第108巻, 社団法人電子情報通信学会 The Institute of Electro *

Also Published As

Publication number Publication date
JP5954846B2 (ja) 2016-07-20
US20150199840A1 (en) 2015-07-16
US9390549B2 (en) 2016-07-12
EP2889001B1 (en) 2017-09-13
WO2014030262A1 (ja) 2014-02-27
EP2889001A1 (en) 2015-07-01
EP2889001A4 (en) 2016-02-17

Similar Documents

Publication Publication Date Title
EP2916738B1 (en) Lung, lobe, and fissure imaging systems and methods
JP2021529051A (ja) 3dオブジェクトの正準ポーズの自動化判定、および深層学習を使った3dオブジェクトの重ね合わせ
US8483462B2 (en) Object centric data reformation with application to rib visualization
JP2017174039A (ja) 画像分類装置、方法およびプログラム
WO2008091583A2 (en) Image-based extraction for vascular trees
JP2004520923A (ja) デジタル画像をセグメント化する方法
JP6458166B2 (ja) 医用画像処理方法及び装置及びシステム及びプログラム
WO2018097880A1 (en) Systems and methods for an integrated system for visualizing, simulating, modifying and 3d printing 3d objects
Jaffar et al. Anisotropic diffusion based brain MRI segmentation and 3D reconstruction
JP2012245085A (ja) 画像処理プログラム、方法及び装置
WO2015150320A1 (en) Segmentation of tubular organ structures
CN111476764B (zh) 一种运动模糊ct图像三维重建的方法
JP5954846B2 (ja) 形状データ生成プログラム、形状データ生成方法及び形状データ生成装置
CN109313803B (zh) 一种用于映射对象的身体的至少部分的图像中的结构的至少部分的方法和装置
CN111918611B (zh) 胸部x线图像的异常显示控制方法、记录介质及装置
JP2005511177A (ja) 体の構造の孤立した視覚化されたものを形成する方法及び装置
CN113889238B (zh) 一种图像识别方法、装置、电子设备及存储介质
EP3989172A1 (en) Method for use in generating a computer-based visualization of 3d medical image data
Boskamp et al. Geometrical and structural analysis of vessel systems in 3D medical image datasets
JP4571378B2 (ja) 画像処理方法および装置並びにプログラム
Abdolali et al. Mandibular canal segmentation using 3D Active Appearance Models and shape context registration
CN113962957A (zh) 医学图像处理方法、骨骼图像处理方法、装置、设备
Alom et al. Automatic slice growing method based 3D reconstruction of liver with its vessels
Reska et al. Fast 3D segmentation of hepatic images combining region and boundary criteria
JP2022505428A (ja) アルゴリズムによるセグメンテーションの正確性の予測

Legal Events

Date Code Title Description
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: 20160607

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160609

R150 Certificate of patent or registration of utility model

Ref document number: 5954846

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150