JP7169514B2 - 診断支援プログラム - Google Patents

診断支援プログラム Download PDF

Info

Publication number
JP7169514B2
JP7169514B2 JP2022066013A JP2022066013A JP7169514B2 JP 7169514 B2 JP7169514 B2 JP 7169514B2 JP 2022066013 A JP2022066013 A JP 2022066013A JP 2022066013 A JP2022066013 A JP 2022066013A JP 7169514 B2 JP7169514 B2 JP 7169514B2
Authority
JP
Japan
Prior art keywords
image
lung
diaphragm
frequency
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.)
Active
Application number
JP2022066013A
Other languages
English (en)
Other versions
JP2022095871A (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.)
Mediott Co., Ltd.
Paramevia Pte Ltd
Original Assignee
Mediott Co., Ltd.
Radwisp Pte Ltd
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 Mediott Co., Ltd., Radwisp Pte Ltd filed Critical Mediott Co., Ltd.
Publication of JP2022095871A publication Critical patent/JP2022095871A/ja
Application granted granted Critical
Publication of JP7169514B2 publication Critical patent/JP7169514B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis 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
    • A61B6/5217Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data extracting a diagnostic or physiological parameter from medical diagnostic data
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • A61B6/5264Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to motion
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/0033Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/0033Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
    • A61B5/004Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/05Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/055Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves  involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/72Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235Details of waveform analysis
    • A61B5/7253Details of waveform analysis characterised by using transforms
    • A61B5/7257Details of waveform analysis characterised by using transforms using Fourier transforms
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/72Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7271Specific aspects of physiological measurement analysis
    • A61B5/7275Determining trends in physiological measurement data; Predicting development of a medical condition based on physiological measurements, e.g. determining a risk factor
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03Computed tomography [CT]
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5288Devices using data or image processing specially adapted for radiation diagnosis involving retrospective matching to a physiological signal
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/54Control of apparatus or devices for radiation diagnosis
    • A61B6/541Control of apparatus or devices for radiation diagnosis involving acquisition triggered by a physiological signal
    • 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
    • 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
    • G16H30/00ICT specially adapted for the handling or processing of medical images
    • G16H30/40ICT specially adapted for the handling or processing of medical images for processing medical images, e.g. editing
    • 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
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03Computed tomography [CT]
    • A61B6/032Transmission computed tomography [CT]
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/503Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of the heart
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10116X-ray image
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20021Dividing image into blocks, subimages or windows
    • 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/20048Transform domain processing
    • G06T2207/20056Discrete and fast Fourier transform, [DFT, FFT]
    • 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
    • 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/30061Lung
    • 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/30101Blood vessel; Artery; Vein; Vascular
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V2201/00Indexing scheme relating to image or video recognition or understanding
    • G06V2201/03Recognition of patterns in medical or anatomical images
    • G06V2201/031Recognition of patterns in medical or anatomical images of internal organs

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Public Health (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biophysics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Physiology (AREA)
  • Optics & Photonics (AREA)
  • Epidemiology (AREA)
  • Artificial Intelligence (AREA)
  • Psychiatry (AREA)
  • Signal Processing (AREA)
  • Primary Health Care (AREA)
  • Quality & Reliability (AREA)
  • Mathematical Physics (AREA)
  • Databases & Information Systems (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)
  • Measuring And Recording Apparatus For Diagnosis (AREA)
  • Image Processing (AREA)

Description

本発明は、人体の画像を解析し、解析結果を表示する技術に関する。
医師が胸部の動態画像によって肺の診断を行なう際、被写体が自然呼吸の状態で撮影された時系列の胸部動態画像の観察が重要である。生理学的なデータを取得しやすいスパイロメータ、RI(Radio Isotope)検査、形態的なデータが得られる単純X線写真、CT(Computed Tomography)などが肺機能を評価するための手法として知られている。しかし、生理学的なデータと形態的なデータの両者を効率良く取得することは容易ではない。
近年、FPD(Flat panel detector)等の半導体イメージセンサを利用して、人体の胸部の動態画像を撮影し、診断に用いる方法が試みられている。例えば、非特許文献1には、動態画像を構成する複数のフレーム画像の間で、信号値の差を示す差分画像を生成し、その差分画像から各信号値の最大値を求めて表示する技術が開示されている。
また、特許文献1には、人体の胸部の動態を示す複数のフレーム画像のそれぞれのフレーム画像から肺野領域を抽出し、その肺野領域を複数の小領域に分割し、複数のフレーム画像の間において、分割された小領域を互いに対応付けて解析する技術が開示されている。この技術によれば、分割された小領域の動きを示す特徴量が表示される。
特許第5874636号明細書
しかしながら、非特許文献1に記載されている技術のように、単に動態画像の画素ごとのフレーム間差分値の最大値を表示するだけでは、医師が病態を把握することは容易ではない。また、特許文献1に記載されている技術のように、特徴量を表示するだけでは、やはり病態の把握には十分ではない。このため、呼吸や肺血管の状態に即した画像を表示することが望ましい。すなわち、被写体である人体の呼吸状態および血管動態全体を把握し、呼吸、心臓、肺門部の血管または血流の波形若しくは周波数、または画像の変化傾向に基づいて、実際の動きを示す画像を表示することが望ましい。
本発明は、このような事情に鑑みてなされたものであり、呼気または吸気の全部または一部を含む呼吸要素ごとに形状が変化する領域の動きを表示することが可能な診断支援プログラムを提供することを目的とする。より具体的には、計測しようとしている新たな対象のデータに対し、既に取得している波の形およびHzに対する一致率やその他の不一致率を数値化することによって、診断の補助となる数値を算出し、さらに、これらの数値を画像化することにより、診断の補助となる画像を生成することを目的とする。
(1)上記の目的を達成するために、本願は、以下のような手段を講じた。すなわち、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の少なくとも一つの周波数を特定する処理と、前記特定した呼吸要素の少なくとも一つの周波数に基づいて、肺野を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記呼吸要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
(2)また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルのうち、ノイズの周波数を含み、前記フレーム画像から得られる呼吸要素の周波数以外の周波数、または入力された周波数若しくは周波数帯域に対応するスペクトルを含む一定の帯域内のスペクトルを、フィルタを用いて抽出する処理をさらに含むことを特徴とする。
(3)また、本発明の一態様に係る診断支援プログラムは、前記呼吸要素の周波数および前記各フレーム画像に基づいて、前記フレーム間の画像を生成する処理と、をさらに含むことを特徴とする。
(4)また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数を特定する処理と、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の少なくとも一つの周波数を特定する処理と、前記特定した呼吸要素の少なくとも一つの周波数に基づいて、肺野を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記心血管拍要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
(5)また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数を特定する処理と、肺野を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記心血管拍要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
(6)また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルのうち、ノイズの周波数を含み、前記フレーム画像から得られる心血管拍要素の周波数以外の周波数、または入力された周波数若しくは周波数帯域に対応するスペクトルを含む一定の帯域内のスペクトルを、フィルタを用いて抽出する処理をさらに含むことを特徴とする。
(7)また、本発明の一態様に係る診断支援プログラムは、前記特定した心血管拍要素の周波数および前記各フレーム画像に基づいて、前記フレーム間の画像を生成する処理をさらに含むことを特徴とする。
(8)また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、被写体の血管拍から抽出される血管拍要素の少なくとも一つの周波数を特定する処理と、前記各フレーム画像について設定された解析範囲を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画像の変化を計算する処理と、前記各フレーム画像における各ブロック領域の画像の変化をフーリエ変換する処理と、前記フーリエ変換後に得られるスペクトルのうち、前記血管拍要素の少なくとも一つの周波数に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する処理と、前記一定の帯域から抽出したスペクトルに対して逆フーリエ変換する処理と、前記逆フーリエ変換後の各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
(9)また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルのうち、ノイズの周波数を含み、前記フレーム画像から得られる血管拍要素の周波数以外の周波数、または入力された周波数若しくは周波数帯域に対応するスペクトルを含む一定の帯域内のスペクトルを、フィルタを用いて抽出する処理をさらに含むことを特徴とする。
(10)また、本発明の一態様に係る診断支援プログラムは、前記特定した心血管拍要素の周波数および前記各フレーム画像に基づいて、前記フレーム間の画像を生成する処理をさらに含むことを特徴とする。
(11)また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の少なくとも一つの周波数を特定する処理と、前記特定した呼吸要素の少なくとも一つの周波数に基づいて、肺野および横隔膜を検出する処理と、前記検出した肺野を複数のブロック領域に分割し、前記各フレーム画像におけるブロック領域の画素の変化率を算出する処理と、前記ブロック領域の画素の変化率と、呼吸と連動する動的部位の変化率との比の値である同調率を用いて、前記同調率が予め定められた一定の範囲内にあるブロック領域のみを抽出する処理と、前記抽出したブロック領域のみを含む各画像をディスプレイに表示する処理と、をコンピュータに実行させることを特徴とする。
(12)また、本発明の一態様に係る診断支援プログラムは、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数または血管拍から抽出される血管拍要素の少なくとも一つの周波数を特定する処理をさらに含むことを特徴とする。
(13)また、本発明の一態様に係る診断支援プログラムは、前記同調率の対数の値が、0を含む一定の範囲として定められることを特徴とする。
(14)また、本発明の一態様に係る診断支援プログラムは、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲線(Bezier curve)を用いて、他のフレームにおける肺野を検出する処理をさらに含むことを特徴とする。
(15)また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野内に内部制御点を選定し、前記肺野内の内部制御点を通る曲線または直線によって前記肺野を分割することを特徴とする。
(16)また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野の外延およびその近傍における制御点の間隔を相対的に大きくし、前記検出した肺野内における部位毎の膨張比率に応じて、前記内部制御点の間隔を相対的に小さくすることを特徴とする。
(17)また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野において、制御点の間隔を、人体に対して頭尾方向に進むに従って相対的に大きくし、または特定のベクトル方向に従って相対的に大きくすることを特徴とする。
(18)また、本発明の一態様に係る診断支援プログラムは、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲面(Bezier surface)を用いて、他のフレームにおける肺野を検出する処理をさらに含むことを特徴とする。
(19)また、本発明の一態様に係る診断支援プログラムは、特定のフレームにおいて予め定められた解析範囲上に、少なくとも一つ以上のベジエ曲線(Bezier curve)を用いて、他のフレームにおいて前記解析範囲に対応する範囲を検出する処理をさらに含むことを特徴とする。
(20)また、本発明の一態様に係る診断支援プログラムは、少なくとも一つ以上のベジエ曲線(Bezier curve)を用いて、少なくとも肺野、血管または心臓を描画する処理をさらに含むことを特徴とする。
(21)また、本発明の一態様に係る診断支援プログラムは、人体の画像を解析し、解析結果を表示する診断支援プログラムであって、前記画像を格納するデータベースから複数のフレーム画像を取得する処理と、前記取得したすべてのフレーム画像について、ベジエ曲線を用いて解析範囲を特定する処理と、前記解析範囲内のインテンシティ(intensity)の変化に基づいて解析対象を検出する処理と、をコンピュータに実行させることを特徴とする。
(22)また、本発明の一態様に係る診断支援プログラムは、前記検出した解析対象の辺縁の特徴を算出する処理をさらに含むことを特徴とする。
(23)また、本発明の一態様に係る診断支援プログラムは、連続する各画像について、インテンシティ(intensity)の差分を算出することで横隔膜を検出し、前記検出した横隔膜または呼吸と連動する動的部位の位置または形状を示す指標を表示することを特徴とする。
(24)また、本発明の一態様に係る診断支援プログラムは、インテンシティ(intensity)の閾値を変化させることで、横隔膜以外の部位によって遮られていない横隔膜を表示し、横隔膜の全体形状を補間することを特徴とする。
(25)また、本発明の一態様に係る診断支援プログラムは、前記検出した横隔膜の位置若しくは形状、または呼吸と連動する動的部位の位置若しくは形状から、前記呼吸要素の少なくとも一つの周波数を計算する処理と、をさらに含むことを特徴とする。
(26)また、本発明の一態様に係る診断支援プログラムは、前記検出した肺野を空間的に正規化し、またはリコンストラクション(reconstruction)を利用して時間的に正規化する処理をさらに含むことを特徴とする。
(27)また、本発明の一態様に係る診断支援プログラムは、前記呼吸要素の少なくとも一つの周波数の位相を変化させ、または呼吸要素の波形を円滑化させることで、呼吸要素を補正することを特徴とする。
(28)また、本発明の一態様に係る診断支援プログラムは、解析範囲内のいずれかの部位の波形を特定し、前記特定した波形の周波数の構成要素を抽出し、前記波形の周波数の構成要素に対応する画像を出力することを特徴とする。
(29)また、本発明の一態様に係る診断支援プログラムは、解析範囲のデンシティ(density)を検出し、デンシティが相対的に大きく変化する箇所を除去することを特徴とする。
(30)また、本発明の一態様に係る診断支援プログラムは、前記フーリエ変換後に得られるスペクトルから、臓器特有の周期的な変化のスペクトル構成比に基づいて、逆フーリエ変換を行なう際の少なくとも一つの周波数を選択する処理をさらに含むことを特徴とする。
(31)また、本発明の一態様に係る診断支援プログラムは、前記呼吸要素の少なくとも一つの周波数に応じて、X線の照射間隔を調整するよう、X線撮影装置を制御することを特徴とする。
(32)また、本発明の一態様に係る診断支援プログラムは、前記逆フーリエ変換後に、振幅値が相対的に大きいブロックのみを抽出して表示することを特徴とする。
(33)また、本発明の一態様に係る診断支援プログラムは、前記肺野を同定した後、横隔膜または胸郭を特定し、横隔膜または胸郭の変化量を算出し、前記変化量から変化率を算出する処理をさらに含むことを特徴とする。
(34)また、本発明の一態様に係る診断支援プログラムは、特定のスペクトルに係数を乗算する処理をさらに含み、前記係数が乗算された特定にスペクトルに基づいて強調表示を行なうことを特徴とする。
(35)また、本発明の一態様に係る診断支援プログラムは、画像を格納するデータベースから複数のフレーム画像を取得した後、呼吸要素の周波数または波形を特定するために、解析対象となる部位にデジタルフィルタを施すことを特徴とする。
(36)また、本発明の一態様に係る診断支援プログラムは、前記各フレーム画像の特定領域の画素に基づいて、呼気または吸気の全部または一部を含む呼吸要素の複数の周波数を特定し、前記呼吸要素の複数の周波数のそれぞれに対応する各画像をディスプレイに表示することを特徴とする。
(37)また、本発明の一態様に係る診断支援プログラムは、ある一枚以上のフレーム画像の特定の範囲について、ある一定の値に集簇する画像を選択し、ディスプレイに表示することを特徴とする。
本発明の一態様によれば、呼気または吸気の全部または一部を含む呼吸要素ごとに形状が変化する領域の動きを表示することが可能となる。
本実施形態に係る診断支援システムの概略構成を示す図である。 肺領域の分割手法の一例を示す図である。 肺の形態が時間の経過によって変化する様子を示す図である。 肺の形態が時間の経過によって変化する様子を示す図である。 特定ブロックの「intensity」変化と、それをフーリエ解析した結果を示す図である。 心拍に近い周波数成分を抜き出したフーリエ変換結果と、これを逆フーリエ変換して心拍に近い周波数成分の「intensity」変化を示す図である。 フーリエ変換後に得られたスペクトルのうち、ある一定の帯域を抽出する例を示す図である。 肺の変化率を模式的に表した図である。 肺野領域のパターン画像の例を示す図である。 肺野領域のパターン画像の例を示す図である。 肺野領域のパターン画像の例を示す図である。 肺野領域のパターン画像の例を示す図である。 肺野の輪郭を、ベジエ曲線および直線の両方を用いて描画した例を示す図であり、肺野が最大の状態を示す。 肺野の輪郭を、ベジエ曲線および直線の両方を用いて描画した例を示す図であり、肺野が最小の状態を示す。 前と次のフレームの間において、肺野の画像の前後を重ねた図である。 図4Aの2枚の元画像の差分を取った結果、「ギャップの強いline」が生じた状態を示す図である。 図4Bにおいて、画像の上下方向の各位置における「intensity」値の合計「density」の差分値を示す図である。 曲線回帰を行ない、横隔膜の相対的位置を近似した結果を示す図である。 本実施形態に係る呼吸機能解析の概要を示すフローチャートである。 ディスプレイに表示される画像の一例を示す図である。 ディスプレイに表示される画像の一例を示す図である。 本実施形態に係る肺血流解析の概要を示すフローチャートである。 本実施形態に係るその他の血流解析の概要を示すフローチャートである。 フーリエ変換後に得られたスペクトルのうち、ある一定のスペクトルに係数を掛ける例を示す図である。 肺野を、ベジエ曲線を用いて描画した例を示す図である。 肺野を、ベジエ曲線を用いて分割した例を示す図である。 肺野を、ベジエ曲線を用いて分割した例を示す図である。 大動脈血流量の波形と心室容量の波形とを対比した一例を示す図である。 肺と肺の近辺の画素値の一例を示す図である。 人体の血管の概略構成を模式化した図である。
まず、本発明の基本的な概念について説明する。本発明では、人体における呼吸や血管、肺野の面積及び体積、その他の生体運動において、一定の周期で反復するように捉えられる動きに対し、全体若しくはある部分的な範囲について、時間軸における一定の反復若しくは一定の運き(ルーチーン)を、波として捉え、計測する。波の計測結果については、(ア)波の形態自体、若しくは、(イ)波の間隔(周波数:Hz)を用いる。この2つの概念をまとめて「ベースデータ」と呼称する。
同時期に同じようにリンクされる波が存在し得る。例えば、呼吸であれば、以下の近似を概念することができる。
(ある大雑把な範囲の「density」変化の平均)≒(胸郭の変化)≒(横隔膜の動き)≒(肺機能検査)≒(胸腹呼吸センサ)
上記「(ア)波の形態自体」について、「波形同調性」という概念を用い、これに基づいて画像を表示する(Wave form tunable imaging)。また、上記「(イ)波の間隔(周波数:Hz)」について、「周波数同調性」という概念を用い、これに基づいて画像を表示する(Frequency tunable imaging)。
例えば、心臓の場合、図13に示す「大動脈血流量の波形と心室容量の波形とを対比した一例」のように、大動脈血流量のピークと、心室容積のピークや波形は一致しない。しかし、図13中、時刻t1からt2、時刻t2からt3、時刻t3からt4…といったように等間隔の時間的幅を1サイクルと定めると、大動脈血流量の1サイクルおよび心室容量の1サイクルが何度も繰り返されることとなり、各波形は周波数が同調していると言える。この波形について着目すると、図13に示すような実測値から1サイクルを特定し、モデル波形を利用することによって、波の形(Wave form)を予想することができる。すなわち、「ベースデータとしての波形」の作り方として、実測しても良いし、周波数(サイクル)から生成しても良いし、モデル波形を利用しても良いし、個人間の波形を平均化して利用しても良い。心臓などの周波数を有する臓器のサイクル(周期)が分かれば、波の形(Wave form)が予想できるため、大動脈血流量や心室容量などの波形を把握し、この波形に基づいて、臓器の動的な画像を表示することが可能となる。
なお、呼吸、心臓、肺門などの「density」の変化を取得する際、他の要素が混ざらないよう、予めデジタルフィルタをかけても良い。
また、本発明では、「呼吸要素」という概念を用いる。「呼吸要素」とは、呼気または吸気の全部または一部を含む。例えば、「1呼吸」を、「1呼気」と「1吸気」に分けて考えることもできるし、「1呼気または1吸気」の「0%、10%、20%、30%、40%、50%、60%、70%、80%、90%、100%」のいずれかに限定して考えることもできる。さらに、各呼気の一定の割合のみ、例えば、呼気の10%のみを抽出して評価することも可能である。これらのいずれかのデータ、またはこれらを組み合わせたデータを用いて、より精度の高い画像の抽出を可能とする。この際、何度も相互的に計算することもある。
このような考え方は、「呼吸要素」のみならず、「心血管要素」についても同様に適用することが可能である。
ここで、ベースデータを作る際、単一若しくは複数のモダリティから得られる特徴量(例えば、ある一定範囲の「density」、「volumetry」より構成される変化量、胸郭の動き、横隔膜の動き、「spirometry」、胸腹呼吸センサの2つ以上)、または、同じ呼吸周回などの複数回の波形測定により、お互いの成分抽出を補足し合い、精度を高めていく。これにより、アーティファクトの軽減、ライン(line)などのある一定の予想をもとに精度を高めることが可能となる。ここで、「density」とは、「密度」と訳されるが、画像においては、特定の領域における画素の「吸収値」を意味する。例えば、CTでは、空気は「-1000」、骨は「1000」、水は「0」として用いられている。
また、お互いの成分抽出による波の軸、幅、範囲およびHzの揺らぎ、幅を推定する。すなわち、複数回の重ね合わせによって、Hzの軸設定が平均化、分散によって軸、幅、範囲、Hzの最適レンジ(range)が計算される。その際に他の行動のHz(ノイズ)が抽出され、その波があればそれが入らない程度も相対的に計測していく場合がある。すなわち、波形要素の全体のうち、一部の波形のみを抽出する場合がある。
本明細書では、「density」と「intensity」とを区別して用いる。「density」は、上述したように、吸収値を意味し、XPやXP動画の元画において、空気の透過性が高く、透過性の高い部分を白であることを数値にし、空気を「-1000」、水を「0」、骨を「1000」として表示するものとする。一方、「intensity」は、「density」から相対的に変化したもの、例えば、normalizedして濃度の幅、信号の程度に“変換”して表示したものとする。すなわち、「intensity」は、画像において、明暗や強調度などの相対的な値である。XP画像の吸収値を直接扱っている間は、「density」または「densityの変化(Δdensity)」として表す。そして、これを画像表現上の都合で、上記のような変換を行なって、「intensity」として表す。例えば、0から255の256階調にカラー表示する場合は、「intensity」となる。このような用語の区別は、XPやCTの場合にあてはまる。
一方、MRIの場合は、空気を「-1000」、水を「0」、骨を「1000」と定めようとしても、MRIの画素値、測定機械の種類、測定時の人の体調、体つき、測定時間によって、値が非常に変化してしまうという事情があり、また、T1強調像などMRIの信号の採り方についても、その施設、測定機械の種類によって様々であり、一定ではない。このため、MRIの場合は、XPやCTの場合のような「density」定義ができない。このため、MRIでは、最初に描出する段階から、相対値を取り扱うこととしており、最初から「intensity」として表現する。そして、その処理する信号も「intensity」である。
以上により、ベースデータを得ることが可能となる。上記ベースデータに対し、計測したい新たな対象について、上記ベースデータの波形、波のHzのある一定の幅、範囲で抽出する。例えば、呼吸抽出のみや、血管抽出程度の幅、範囲、波形要素で抽出する。なお、この波形、Hzの幅に関しては、他の機能における波形要素、ノイズなどの「artifact」、他の同調性があると思われる他の「modality」の波形、複数回行なう再現性など用いて、相対的に、また、統計をもとに総合的に判断される。そこに調整、経験が必要となる(機械学習を適用することも可能である)。これは、幅、範囲を広げると他の機能の要素が入り始める一方、狭すぎると機能自体の要素がそぎ落とされてしまうので、そのレンジに関しては、調整が必要となるからである。例えば、複数回のデータがあると、レンジ、Hzと測定一致幅などが規定しやすい。
[同調一致率について]
本明細書では、画像変化の傾向を、同調一致率として説明する。例えば、肺野を検出し、複数のブロック領域に分割し、各フレーム画像におけるブロック領域の「平均density(画素値x)」を算出する。そして、「平均density(画素値x)」の最小値から最大値の変化幅(0%~100%)に対する各フレーム画像におけるブロック領域の平均画素値の割合(x’)を算出する。一方、横隔膜の最小位置から最大位置の変化幅(0%~100%)に対する各フレーム画像の横隔膜の変化(y)の割合(y’)との比の値(x’/y’)を用いて、比の値(x’/y’)が予め定められた一定の範囲内にあるブロック領域のみを抽出する。
ここで、y’=x’若しくはy=ax(aは横隔膜の振幅の数値や「density」の数値の係数)となる場合は、完全一致である。しかし、完全一致の場合のみが有意義な値であるわけではなく、ある一定の幅を持った値を抽出すべきである。そこで、本発明の一態様では、対数(log)を用いて、一定の幅を以下のように定める。すなわち、y=xの割合(%)で計算すると、同調の完全一致は「log y’/x’=0」である。さらに、同調一致率の範囲が狭いもの(数式的に狭い)範囲を抽出する場合は、例えば、0に近い範囲で「log y’/x’=-0.05~+0.05」と定め、同調一致率の範囲が広いもの(数式的に広い)範囲のものであれば例えば0に近い範囲で「log y’/x’=-0.5~+0.5」と定める。すなわち、同調率の対数の値が、0を含む一定の範囲として定められる。この範囲が狭いものであればあるほど、また、その範囲内で一致する数値が高いほど、一致率が高いといえる。
画素のピクセルごとにこの比の値を求めて個数をカウントすると、健康な人の場合は、完全一致の場合をピークとした正規分布が得られる。これに対し、疾患を有する人の場合は、この比の値の分布が崩れることとなる。なお、上記のように、対数を用いて幅を定める手法は、あくまでも一例であり、本発明はこれに限定されない。すなわち、本発明は、(ある大雑把な範囲の「density」の変化)≒(胸郭の変化)≒(横隔膜の動き)≒(肺機能検査)≒(胸腹呼吸センサの動き)≒(肺野の面積及び体積)として、“画像抽出”を行なうものであり、対数を用いる手法以外の手法も適用可能である。このような手法により、同調性画像を表示することが可能となる。
血管の場合は、一連の心臓の収縮(y)に呼応して生じる一連の「density」の変化(x(肺門部における一波形)において、そのままの形でわずかな時間の遅れ(位相の変化)が存在するため、y=a’(x-t)と表される(すなわち、Y≒X)。完全一致の場合は、t=0であるため、y=x、またはy=a’xである。横隔膜の場合と同様に、同調一致率の範囲が狭いもの(数式的に狭い)範囲を抽出する場合は、例えば、0に近い範囲で「log y’/x’=-0.05~+0.05」と定め、同調一致率の範囲が広いもの(数式的に広い)範囲のであれば例えば0に近い範囲で「log y’/x’=-0.5~+0.5」と定める。この範囲が狭いものであればあるほど、また、その範囲内で一致する数値が高いほど、一致率が高いといえる。
その他の血管の場合は、上記の「心臓に呼応する部分」が除外され、肺門からプロットした中枢側の「density」を用いる。末梢の血管の場合も同様に取り扱うことができる。
さらに、循環器についても本発明を適用することができ、例えば、心臓の「density」の変化が、肺門部~末梢肺野への血流の「density」の変化に直接関連し、一連の心臓の「density」の変化や肺門部の「density」の変化は、一種の変換を受けてそのまま伝播される。それは、心臓の「density」の変化と肺門部の「density」の変化の関係より若干の位相の差を得て生じると考えられる。また、肺門部などの「density」の変化が、そのまま肺野の血流への「density」の変化に関連するので、そのままの率で反映したもの(Y≒Xの一致率の関係)で同調性を表現することも可能である。また、頸部血管系や、胸部、腹部、骨盤、四肢などの大血管系においても、同様に、近傍の中枢の心臓血管でプロットした「density」の変化が、直接関連し、またはわずかな位相を伴って関連していると考えられる。そして、その「density」が、背景に応じて変動し、伝搬するときには「density」の変化の様が伝わるとして、同調一致率として考察することが可能となる。
ここで、1枚の画像の変化量と1枚の画像の変化率のそれぞれにおいて、「吸気量合計≒呼気量合計」とすることができる。そこで、周囲空気の透過性との差から相対的に数値を出す場合、肺野の「density」から変化量を1としたときの相対的な値(Standard Differential Signal Density/Intensity)として表示しようとすると、(1)画像1枚ごとの差の画像で、1枚ごとに1とした時の画像(通常想定)、(2)1枚ごとの差の画像で「density(変化量や変化率)」を足した吸気全体もしくは呼気全体、もしくは吸気呼気の絶対値を1としたときの割合、さらに、(3)複数回の撮影におけるそれぞれの呼吸時(10%の時を数回select)における「density」総量を1としてその割合、について、それぞれ、変化量、変化率の描出を行なうことができる。
また、MRなどの3Dの場合であるが、吸気全体の「intensity(MRの場合)」や「density」(CTの場合)を合計した値(その際はそれを1としたとき)、その「intensity」や「density」の差は吸気(安静時や努力呼吸でも)の「peak flow volume deta」に換算でき、その値をその「intensity」や「density」の割合を出すことによって、少なくともMRIやCTなどでの「3D×time」の計算する場合に、各肺野の部分における実測呼吸量、呼吸率を換算することができる。同様に、1回心拍出量を入力することで、肺野の「flow」における「capillary phase」における分布が、肺血流末梢量の分布、容量に換算する推定値を提示することも可能である。
すなわち、(画像1枚あたりの吸気変化量)×(吸気の枚数全部)≒(画像1枚当たりの呼気変化量)×(呼気の枚数全部)≒(その時の吸気呼吸:自然呼吸もしくは努力呼吸のvolume)≒(その時の呼気呼吸:自然呼吸もしくは努力呼吸のvolume)≒(その時の自然呼吸もしくは努力呼吸の「volume」の吸気もしくは呼気の変化量)が成り立つ。10%や20%の1枚の変化量だけを取り出す場合は、(すべての枚数)×(その時間の変化量)を計算することによって、推定値を算出することが可能である。
その抽出変化量を可視化して、画像に描出する。これが、以下に説明する呼吸機能解析、血管解析である。そして、胸郭や横隔膜の変化率を可視化する。その際に、再度結果に対するアーティファクトを除外し、新しいデータ抽出波形や最初のベースとなるデータ波形、その他のモダリティなどの波形、周囲、複数回の波形から抽出して、機能の抽出を行なう場合もある。アーティファクトを除外する手法は、後述する。
また、上記抽出したもの以外から抽出した変化成分を除外したものでも特徴量を把握する場合がある。例えば、腹部腸管の動きを把握する際、腹部から呼吸の影響と血管の影響を除外して、腹部腸管の動きの抽出を図る。
以下、本発明の実施形態について図面を参照して説明する。図1Aは、本実施形態に係る診断支援システムの概略構成を示す図である。この診断支援システムは、コンピュータに診断支援プログラムを実行させることにより特定の機能を発揮する。基本モジュール1は、呼吸機能解析部3、肺血流解析部5、その他の血流解析部7、フーリエ解析部9、波形解析部10および視覚化・数値化部11から構成されている。基本モジュール1は、入力インタフェース13を介してデータベース15から画像データを取得する。データベース15には、例えば、DICOM(Digital Imaging and COmmunication in Medicine)による画像が格納されている。基本モジュール1から出力された画像信号は、出力インタフェース17を介してディスプレイ19に表示される。次に、本実施形態に係る基本モジュールの機能について説明する。
[呼吸要素の周期解析]
本実施形態では、以下の指標に基づいて呼吸要素の周期を解析する。「呼吸要素」とは、上述したように、呼気または吸気の全部または一部を含む概念である。すなわち、肺野内のある一定領域における「density」/「intensity」、横隔膜の動き、胸郭の動きの少なくとも一つを用いて呼吸要素の少なくとも一つの周波数を解析する。この「呼吸要素の少なくとも一つの周波数」において、呼吸要素が示す周波数スペクトルは一つ以上であり、一定の帯域幅を有する場合を含む概念である。肺野をブロックの集合体と考え、各ブロックから複数の周波数が抽出されることから、本実施形態では、これらを周波数群として処理する。なお、前述したように、ベースデータは、「波の形態自体」、および、「波の間隔(周波数:Hz)」の両方の概念を有するため、波の形態として処理することも可能である。また、X線(その他CT、MRIなどの複数種類のモダリティ)の透過性が高い部位で測定されるある一定のvolume 「density」/「intensity」で構成される範囲、スパイログラムなどの他の測定方法から得られるデータや外部入力情報を用いても良い。
なお、一呼吸毎の解析結果を比較し、複数のデータから傾向を解析して、データの確度を高めることもできる。
また、呼吸要素の少なくとも一つの周波数の位相を変化させ、または呼吸要素の波形を円滑化させることで、呼吸要素を補正することも可能である。この場合、(胸郭、その他横隔膜の動き)≒(胸郭の動き)≒(density)≒(精密肺機能)≒(胸郭センサ)などの動きを用いてその波に位相を合わせる。また、肺野平均の「density」を追跡し、最後の変化は波の形態として波の二乗法などの近似を行ない、波の同定を行なう。ここで、胸部の「density」などの場合、最も変化する値が肺の「density」であるため、画面全体の「density」の評価により、肺の「density」の変化を評価する場合もある。波のプロットをする場合、実際動いている場合と、計測値で位相のずれが生じる場合がある。その場合は位相差を最大、最小値の位置、波の形態全体などで位相を補正する場合がある。
[波形解析]
呼吸要素の波形から、波形の周波数の構成要素を算出することができる。これにより、上述した「波形同調性画像」を取得する。具体的には、解析範囲内のいずれかの部位の波形を特定し、前記特定した波形の周波数の構成要素を抽出し、前記波形の周波数の構成要素に対応する画像を出力する。
[心血管拍解析および血管拍解析]
本実施形態では、以下の指標に基づいて心血管拍解析および血管拍を解析する。すなわち、心電図や脈拍計等の他のモダリティの計測結果、または肺輪郭から心臓・肺門位置・主要血管を特定し、各部位の「density」/「intensity」の変化を用いて血管拍を解析する。また、マニュアルで画像上にプロットし、対象部位の「density」/「intensity」の変化を解析しても良い。そして、心拍または血管拍から得られる心血管拍要素の少なくとも一つの周波数(波形)を特定する。なお、一拍毎の解析結果を比較し、複数のデータから傾向を解析して、データの確度を高めることが望ましい。また、各部位の「density」/「intensity」の抽出は、複数回実施したり、一定の範囲に対して行なうことで精度を高めることが可能となる。また、心血管拍周波数もしくは周波数帯を入力する方法もある。
[肺野同定]
データベース(DICOM)から画像を抽出し、上記の呼吸要素の周期解析結果を用いて、肺輪郭を自動検出する。この肺輪郭の自動検出については、従来から知られている技術を用いることができる。例えば、特開昭63-240832号公報、または特開平2-250180号公報に開示されている技術を用いることが可能である。次に、肺野を複数のブロック領域に分けて、各ブロック領域の変化を計算する。ここで、撮影速度に応じてブロック領域の大きさを定めても良い。撮影速度が遅い場合は、あるフレーム画像の次のフレーム画像で対応する部位が特定しにくくなるため、ブロック領域を大きくする。一方、撮影速度が速い場合は、単位時間当たりのフレーム画像数が多いため、ブロック領域が小さくても追従することが可能となる。また、呼吸要素の周期のうちどのタイミングを選ぶかに応じて、ブロック領域の大きさを計算しても良い。ここで、肺野領域のずれを補正することが必要になる場合がある。その際には、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、また、肺輪郭の相対位置を把握し、その動きに基づいて相対的に評価する。なお、ブロック領域が小さすぎると、画像のちらつきが発生する場合がある。これを防止するため、ブロック領域は一定の大きさを有する必要がある。
上記自動検出した肺野領域に少なくとも一つのベジエ曲線を用いて、肺野を点および制御点の座標として表すことができる。そして、少なくとも一つのベジエ曲線を用いて囲まれた閉じた曲線、いわゆる「単純閉曲線」を、複数用いることによって、肺野を表すことも可能である。同様に、一つまたは複数の単純閉曲線を用いて解析対象を表すことも可能である。
各フレームの肺野は、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲線(Bezier curve)を用いて、他のフレームにおける肺野を検出することもできる。例えば、最大と最小の2つの肺野を検出し、その2つの肺野を用いて、その他のフレームの肺野を算出する方法が挙げられる。ここでは、その他のフレームに「変化率」という変数を定義する。「変化率」は、肺野の大きさ、すなわち呼吸の状態を表現する値であり、横隔膜の位置や画像全体の「intensity」平均値などから算出できる。スパイログラフィーなどの外部装置の計測データを用いて算出したり、モデル化された変化率を用いたりすることも可能である。このように、「変化率」という変数を任意に定めることができるため、例えば、肺野が一定の割合(10%, 20%, 30%…)で変化していると仮定して算出することもできる。このように定義した変化率は、誤差を含んでいる場合があるため、誤差の自動・手動除去を行なった結果、または最小二乗法等を用いて近似を行なった結果などを使用して以後の処理を行なう場合もある。最大肺野と最小肺野まで、線形に変形すると仮定し、それぞれのフレームの変化率を用いて、線形変換などの技法を用いてそれぞれのフレームにおける肺野を算出する。
また、連続するフレームの任意の範囲において、上述の処理を適用することが可能である。例えば、呼吸において、肺野は極大と極小への変化を繰り返すものであるが、実際の測定においては、極大時の形状が常に一定とは限らない。例えば、極大から極小、極小から極大の各範囲において、上述の処理を適用することによって、最大と最小の2つの肺野を定義して計算するよりも、精度高く肺野を算出することができると期待される。なお、ここでは、具体例として、極大と極小を用いて説明したが、本発明はこれに限定されるわけではなく、「任意の範囲」であるため、例えば、呼吸の途中、0%と30%、30%と100%といった位置で行なうことも可能である。
また、精度は低下するものの、1つの肺野から各フレームの肺野を算出することも可能である。例えば、胸郭の形状等から類推することで、肺野の変化ベクトルを規定する。具体的には、ベジエ曲線の制御点のそれぞれに変化ベクトルを規定する手法を採るが、本発明は、これに限定されるわけではない。そして、検出した1つの肺野と変化ベクトル、それぞれのフレームにおける変化率を使用して、それぞれのフレームにおける肺野を算出する。この算出結果に対して自動、または手動で補正を行なうことでさらに精度を高めることができる。また、3Dにおいても、上述の手法は有効である。すなわち、3Dの場合であっても、特定のフレームにおいて検出した肺野上の少なくとも一つ以上のベジエ曲面(Bezier surface)を用いて、他のフレームにおける肺野を検出する処理を実行することも可能である。これにより、フレーム間の肺野の画像を得ることが可能となる。
図6Cは、呼吸要素の周期を示すグラフである。図6Cの画像中、白い垂直線が示されているが、これは、呼吸要素の周期中、現時点の位置を示す直線(指標)であり、図6Bに示す肺の動画の動きに応じて、呼吸要素の周期中の現在位置を示すように動く。呼吸要素の周期の現在位置を示すことによって、肺の動きの周期における現在位置を明確に把握することが可能となる。なお、本発明においては、呼吸要素の周期をグラフで表すだけではなく、血流の「density」、胸郭、横隔膜など、肺の動きと連動するものについては、すべてグラフ化することが可能である。
また、被写体が「息を止めた場合」は、呼吸要素の周波数を特定できない場合がある。この場合は、被写体の心拍または血管拍から抽出される心血管拍要素の少なくとも一つの周波数を用いて、後述するフーリエ解析を行なう。この場合、心臓、横隔膜、または呼吸と連動する動的部位の動き方に応じて、後述するブロック領域の分割の仕方を変えるようにしても良い。
[辺縁の検出とその評価]
本発明は肺の辺縁を検出し、その辺縁を評価することが可能である。例えば、前述の方法で肺野を算出した後、辺縁の位置および形状を、改めて精度高く検出することができる。算出された肺野内の任意の位置に点をプロットして、そこから四方八方に線を延伸し、各線において画素値の変化を評価する。例えば、図14に示すように、肺を切断する線分Sに沿って画素値を算出すると、辺縁で画素が大きく変動することがわかるが、その変動の絶対値は異なる。例えば、左側の辺縁と右側の辺縁検出時の閾値を調整することで、辺縁検出の精度が高まる。また、領域ごとの画素値変動の特性を利用することもできる。図14に示すように、S2領域とS3領域の縁の差分が小さかったとしても、画素値の変動の分散からS2領域とS3領域の縁を特定することができる。ここでは分散に着目したが、本発明はこれに限定されるわけではない。
さらに、同様の考え方によって、肺以外の臓器、血管、腫瘍などの解析範囲の辺縁を検出することが可能になる。例えば、血管中に造影剤が存在する場合、血管の内部は明確に可視化できるが、血管の外側や厚さを明確に算出することは容易でない。本実施形態では辺縁を正確に検出できるため、解析範囲内にある血管の形状、特徴を算出することができる。これにより、従来は把握が容易でなかった血管の厚さや外周を定量的にも把握し、診断に用いることが可能となる。
[ブロック領域の作成]
肺野を複数のブロック領域に分ける手法について説明する。図1Bは、肺野を肺門から放射状に分割する手法を示す図である。肺は、肺尖側よりも横隔膜側の方が大きく動くため、横隔膜側に近いほど粗く分割した点をプロットするようにしても良い。なお、図1Bにおいて、縦方向の線(点線)を追加的に描画し、複数の矩形(正方形)のブロック領域に分けてもよい。これにより、肺の動作をより正確に表すことが可能となる。なお、「肺の縦方向に点をプロットし、肺を横断的に分ける手法」、「肺の横方向に点をプロットし、肺を縦断的に分ける手法」、「肺尖部における接線と横隔膜における接線を引き、その接線が交わる点を中心点として定め、その点を含む直線(例えば、鉛直線)からある一定の角度で引いた線分で肺を分割する手法」、「肺を肺尖(または肺門)から横隔膜端部を結ぶ直線と直交する複数の平面で切断する手法」などの手法で、肺野を分割することも可能である。なお、これらの手法は、三次元立体画像にも適用可能である。3Dの場合は、複数の曲面または平面で囲まれた空間として、各臓器を捉える。臓器をさらに細かく分けることもできる。例えば、右肺の3Dモデルを考えた場合、上葉、中葉、下葉に分けて取り扱うこともできる。
肺野領域は、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、肺輪郭の相対的位置を把握し、その動きに基づいて相対的に評価すべきものである。このため、本願発明では、肺輪郭を自動検出した後、肺輪郭によって特定される領域を複数のブロック領域に分割し、各ブロック領域に含まれる画像の変化の値(画素値)を平均化する。例えば、図10に示すように、ベジエ曲線上で対向する肺の辺縁上に点をプロットし、それらを接続した上で、その中間の点を通る曲線を用いることも可能である。その結果、図1Cに示されているように、肺の形態が時間の経過によって変化しても、注目する領域の経時的変化を追跡することが可能となる。一方、図1Dは、解析対象となる臓器(この場合は肺)の形態を考慮せずに、ブロック領域に分割した場合の経時的変化を示す図である。肺野領域とは、上述した通り、胸郭の動き、横隔膜の動き、肺野全体の血管の位置関係を同定し、肺輪郭の相対的位置を把握し、その動きに基づいて相対的に評価すべきものであるが、図1Dに示すように、肺野領域を特定することなくブロック領域に分割すると、肺の経時的変化により、注目領域が肺野領域から外れ、意味のない画像となってしまう。特に、横隔膜の動きは肺野を縮める動作が強いため、横隔膜や全体的な数値だけを補正するのではなく、胸郭成分やその他の複数の要素を取り込んで、肺野領域を補正することが好ましい。また、呼吸要素周波数もしくは周波数帯を入力する方法もある。3Dについても、同様に領域分割の計算ができる。
さらに、図11に示すように、肺野Aにおいて、ベジエ曲線を用い、検出した肺野内に内部制御点を選定し、肺野内の内部制御点を通る曲線または直線によって肺野を分割することも可能である。すなわち、肺野の枠だけではなく、肺野領域の内部にも制御点を設け、これらの制御点を用いて肺野領域(A)を分割する。この場合、図12に示すように、検出した肺野の外延およびその近傍における制御点の間隔を相対的に大きくし(1)、検出した肺野内における部位毎の膨張比率に応じて、内部制御点の間隔を相対的に小さく(2)するようにしても良い。また、肺野Aにおいて、制御点の間隔を、人体に対して頭尾方向に進むに従って相対的に大きくしても良く、または、特定のベクトル方向に従って相対的に大きくしても良い。このベクトルの定め方は任意であるが、例えば、肺尖から肺野の反対側へ向かう方向に定めても良いし、図1Bに示すように、肺門から肺野の反対側へ向かう方向に定めても良い。また、肺の構造に応じた方向にベクトルを定めることも可能である。このように、肺野の分割の仕方を「不等分割」とすることによって、領域毎の特徴を考慮した画像の表示が可能となる。例えば、肺野の外周は動きが大きく、ずれが大きくなるため、ブロックを大きくする一方、肺野の内部は動きが小さく、ずれも小さくなるため、ブロックを小さくして細かくする。また、例えば、肺野の横隔膜側は動きが大きく、ずれが大きくなるため、ブロックを大きくする一方、肺野の頭部側は動きが小さく、ずれも小さくなるため、ブロックを小さくして細かくする。これにより、表示の精度を高めることが可能となる。この手法は、肺野に限定されるわけではなく、呼吸と連動する動的部位などにも適用することが可能である。このような手法は、肺を肺葉ごとに3次元的に分割する場合にも適用できる。また、横隔膜の下側の部位、例えば、心臓やその他の臓器をベジエ曲線で囲って表示する場合に用いることも可能である。この場合も、心臓やその他の臓器の構造に応じた方向にベクトルを定めて、領域を不等分割することも可能である。
次に、アーティファクトを排除して画像データを補間する。すなわち、解析範囲内に骨などが含まれるとノイズとして表れてしまうため、ノイズカットフィルタを用いてノイズを除去することが望ましい。X線画像においては、通例では、空気を-1000とし、骨を1000としているため、透過性が高い部分は画素値が低く、黒く表示され、透過性の低い部分は、画素値が高く、白く表示される。例えば、画素値を256階調で表す場合、黒は0で白は255となる。肺野領域内で、血管や骨が存在しない位置の周辺は、X線が透過しやすいため、X線画像の画素値が低くなり、X線画像は黒くなる。一方、血管や骨が存在する位置は、X線が透過し難いため、X線画像の画素値が高くなり、X線画像は白くなる。その他CT、MRIにおいても同様のことが言える。ここで、上記の呼吸要素の周期解析の結果から、一呼吸あたりの波形に基づいて、同一位相の値を用いてデータを補間し、アーティファクトを排除することが可能となる。また、「座標が異なること」、「画素値が極端に変動すること」、「周波数やdensityが異常に高くなること」を検出した場合に、それらに対してカットオフを行ない、残りの得られた画像に対して、例えば、最小二乗法等を用いて連続的で滑らかな波の形を同定することによって、横隔膜のHz計算、肺野の調節に使用できるようにしても良い。また、画像を重ね合わせる場合、(1)前後で片方の画像を取得した取得比較画像をその座標のまま重ねる場合と、(2)前後片方の画像をbaseに取得した後、画像を相対的に拡張してその相対的位置情報をbaseに重ねる方法とがある。以上のような手法によって、肺野の形態を修正したり、ブロック領域の画像の変化を修正したりすることが可能となる。その際、再度、結果に対するアーティファクト(artifact)を除外し、新しいデータ抽出波形や最初のベースとなるデータ波形、その他のモダリティなどの波形、周囲、複数回の波形から抽出して、機能の抽出を行なう。その際、回数は一回でも複数回でも良い。
ここで、時間軸における「reconstruction」について説明する。例えば、15f/sの吸気時間が2秒の場合、30+1枚の画像が得られる。その場合、10%ずつの「reconstruction」は、単純に3枚ずつ重ね合わせれば実施できる。その際、例えば、0.1秒が10%で、その画像が0.07秒と0.12秒の写真しか取得していない場合は、0.1秒の「reconstruction」が必要となる。その場合、10%前後の画像の中間の値(両者の平均)値を与えて「reconstruction」を行なう。また、時間軸で捉え、その時間の割合で係数を変えても良い。例えば、時間軸の差があって、0.1秒の撮影の値がなく、0.07秒と0.12秒の撮影時間があるときは、「(その0.07秒の値)×2/5+(0.12秒の値)×3/5」と計算し直して、「reconstruction」を行なうことができる。さらに、呼吸の平均や横隔膜の係数の変化量からその秒における変化位置関係を認識し、その値を係数にして数字の割合を求める。なお、「Maximum Differential Intensity Projection」の0~100%を含み、10%から20%の「reconstruction」や、10%から40%の「reconstruction」など、厚みを持たせて計算することが望ましい。このように、撮影していない部分についても、1呼吸割合での「reconstruction」を行なうことが可能である。なお、本発明は、呼吸のみならず、血流、胸郭の動き、横隔膜、その他これらと連動する一連の動きに対しても同様に「reconstruction」を行なうことが可能である。また、ブロックごと、またはピクセルごとに「reconstruction」をすることも可能である。なお、「Maximum Differential Intensity Projection」の0~100%を含み、10%から20%の「reconstruction」や、10%から40%の「reconstruction」など、厚みを持たせて計算することが望ましい。
また、上述した手法で肺野を検出し、この検出した肺野を正規化しても良い。すなわち、検出した肺野を空間的に正規化し、またはリコンストラクション(reconstruction)を利用して時間的に正規化する。人体の相違によって肺野の大きさや形状が異なるが、これを正規化することで一定の領域内に表示することができる。
[横隔膜および胸郭]
上記のように肺野を同定すると、横隔膜の動きや胸郭についても把握することが可能となる。すなわち、認識した横隔膜のXp上(2D画像)の横隔膜による曲線や胸郭の曲線を細かな座標の集まりとして計算し、その平均や曲線の局部における下方への変化率や変化量、また横隔膜を曲線として「curve fitting」してその変形率を数値化することにより画像からの機能評価の位置づけを行なうことができる。また、横隔膜面以外の胸部で描いた辺縁の曲線についても、同様に細かな座標の集まりとして計算し、その平均や曲線の変化率を数値化することにより画像からの機能評価を行なうことができる。上記の2つの変化率、変化を、相対的・相互連動として評価し、変化率が異なる(同じように連動して動かない部位など)を数値化、画像化することによりmovementの機能評価を行なう。
ここで、「横隔膜および胸郭評価方法」について説明する。まず、横隔膜について、身体の軸(いわゆる正中線)に直交する左右水平線を軸にその動きを表示する。次に、横隔膜のラインを基線に平坦化する。すなわち、横隔膜のラインを水平の直線に合わせる。そして、横隔膜の動きを評価する。さらに、横隔膜のラインを伸ばして平坦化し、曲線の直交する動きを評価する。次に、胸郭外側では、肺尖から横隔膜胸郭角を結んだ線を基線として(軸として)動きを評価する。胸郭のラインを基線として平坦化し、すなわち、胸郭のラインを「肺尖-肋横隔膜角」との直線に合わせて、動きを評価する。胸郭のラインを基線に伸ばして平坦化して曲線の直交する動きを評価する。そして、上記胸郭、横隔膜ラインの曲率や曲率半径を評価する。そして、上記の変化を「変化量」として算出し、この変化量を微分して“変化率”として評価する。
図6Bおよび図6Cは、ディスプレイに表示される画像の一例を示す図である。図6Bでは、左肺の動きが動画として表示される。図6Bの画像中、白い水平線が示されているが、これは、横隔膜の位置を示す直線(指標)であり、動画が再生されると、横隔膜の動きに追従して上下に動くこととなる。このように、横隔膜を検出し、検出した横隔膜の位置を示す指標、すなわち、横隔膜の位置を示す白い水平線を示すことによって、医師による画像診断を行なうことができるようになる。また、横隔膜の一部のみならず、肺野-横隔膜のラインの認識を使用し、すべての点を認識して、左右、外内側などの横隔膜の一領域、また、横隔膜全体の診断を行なうことが可能となる。同様に、横隔膜のみならず、呼吸と連動する動的部位、例えば、胸郭などの動きも同様に、接線位置などの直線や肺野認識による胸郭の直線によって、胸郭の動きを判断することが可能となる。このように、辺縁は動くものという仮定で、連続画像で差分を取ることによって、辺縁を検出することが可能となる。例えば、腫瘍は堅く、その周りは柔らかいことが多い。このため、腫瘍はあまり動かず、その周りは活発に動くこととなるため、差分を取ることによって、腫瘍の辺縁を検出することができる。
また、MRIやCTなどの3D画像においても、横隔膜の面を一つの座標や立体的な曲面として捉え、その座標や曲面を細かな座標の集まり(横隔膜の辺縁の輪郭、平面および座標の集合群)として計算し、その平均や曲面の局部における下方への変化率や変化量、また横隔膜を曲面として「curve fitting」してその変形率を数値化することにより画像からの機能評価の位置づけを行なうことができる。また、横隔膜面以外の胸部で描いた辺縁の曲面についても、同様に細かな座標の集まりとして計算し、その平均や曲面の変化率を数値化することにより、画像からの機能評価を行なうことができる。上記2つの変化率、変化を相対的、相互連動として評価し、変化率が異なる(同じように連動して動かない部位など)を数値化、画像化することによりmovementの機能評価を行なう。
[フーリエ解析]
上記のように解析した呼吸要素の周期および血管拍周期に基づいて、各ブロック領域の「density」/「intensity」の値や、また、その変化量について、フーリエ解析を実施する。図2Aは、特定ブロックの「intensity」変化と、それをフーリエ解析した結果を示す図である。図2Bは、心拍に近い周波数成分を抜き出したフーリエ変換結果と、これを逆フーリエ変換して心拍に近い周波数成分の「intensity」変化を示す図である。例えば、特定ブロックの「intensity」変化をフーリエ変換(フーリエ解析)すると、図2Aに示すような結果が得られる。そして、図2Aに示した周波数成分から、心拍に近い周波数成分を抜き出すと、図2Bの紙面に対して右側に示すような結果が得られる。これを逆フーリエ変換することによって、図2Bの紙面に対して左側に示すように、心拍の変化に同調した「intensity」変化を得ることができる。
図9に示すように、特定のスペクトルに係数をかけて重みを付けることも可能である。例えば、波形同調性を実現するために、この手法を用いることが可能である。すなわち、フーリエ逆変換を行なう時の周波数の選び方として、複数の周波数を選択して、かつその比率を掛けた後、フーリエ逆変換を行なう。例えば、抽出する帯域中、最も周波数の高いスペクトルを強調表示したい場合に、そのスペクトル強度を2倍することが可能となる。この場合、周波数の連続性は無くても構わない。とびとびに存在するスペクトルを選択することが可能である。
また、心臓の「density」の位置を、左肺(内臓逆位などの場合右芯の場合もある)の形態(肺野抽出の形態から左肺のくぼんだ部位の領域)および椎体、横隔膜の位置から類推することが可能となる。この場合、心臓のROIをとって「density」の抽出を行なう。この抽出をする際には、呼吸、血流の相対的なスペクトル値でおおまかな領域を用いて類推する。また、予め、心血管拍で生じるHz帯(心拍40~150Hz、≒0.67Hz~2.5Hz)などを用いて「filtering」を行なうことによって、呼吸やその他の「artifact」による周波数を除去する場合もある。また、心臓の位置も呼吸状況に応じて変わるので、胸郭の位置が変わるにつれて胸郭の形態値から相対的に心臓の位置を変更し、より正確な心血管拍の抽出や肺門、大血管などの抽出を行なうことがある。さらに、横隔膜の動きと同様、規則的に動く心臓の輪郭に基づいて、周波数を計算する方法がある。
ここで、周波数成分からなるスペクトルに対して逆フーリエ変換を行なう際に、呼吸や血流の「density」から特定される周波数要素(呼吸周波数、心血管拍周波数)と、スペクトルの帯域(BPF:band pass filterを用いても良い)とを両方を加味し、または、そのどちらかの要素に基づいて逆フーリエ変換を行なうようにしても良い。また、上記フーリエ変換後に得られるスペクトルから、臓器特有の周期的な変化のスペクトル構成比に基づいて、逆フーリエ変換を行なう際の少なくとも一つの周波数を選択してもよい。さらに、フーリエ変換後に得られる複数の周波数の構成割合によって、特定の臓器や解析対象となる領域の波形を特定することも可能である(波形同調性画像の作成)。
なお、フーリエ変換を実行する際には、短時間で計算ができるように、AR法(Autoregressive Moving average model)を用いることが可能である。AR法では、自己回帰移動平均モデルにおいて、ユールウォーカー方程式(Yule-walker equiation)やカルマンフィルタを用いる方法があり、そこで導きだされるユールウォーカー推定値(Yule-walker estimates)、PARCOR法、最小二乗法を用いて、計算を補足することができる。これにより、より早く、リアルタイムに近い画像を取得したり、計算の補助やアーティファクト(artifact)の補正を行なったりすることが可能となる。このようなフーリエ解析により、各ブロック領域における画像の性質を抜き出して表示することが可能となる。
また、このフーリエ解析の際に、「デジタルフィルタ」を用いる手法を採ることも可能である。すなわち、元の波形にフーリエ変換を行なって、各スペクトルのパラメータを取得し、元の波に演算処理を施す「デジタルフィルタ」を用いる。この場合は、逆フーリエ変換を行わず、デジタルフィルタを用いる。
ここで、各フレーム画像における各ブロック領域の画像の変化をフーリエ変換し、フーリエ変換後に得られるスペクトルのうち、呼吸要素の周期に対応するスペクトルを含む一定の帯域内のスペクトルを抽出することができる。図2Cは、フーリエ変換後に得られたスペクトルのうち、ある一定の帯域を抽出する例を示す図である。合成波のスペクトルの周波数fは、合成元となる各周波数f(呼吸成分)、f(血流成分)との間に、「1/f=1/f+1/f」という関係が成り立っており、スペクトルを抽出する際に、以下の方法を採ることが可能である。
(1)血流のスペクトル比率が高い部分を抽出する。
(2)呼吸/血流に対応するスペクトルのピークとその近辺の複数の合成波のピークの中間で区切り、スペクトルを抽出する。
(3)呼吸/血流に対応するスペクトルのピークとその近辺の複数の合成波のスペクトルの谷の部分で区切り、スペクトルを抽出する。
上述したように、本発明では、固定的なBPFを用いているわけではなく、呼吸要素の周期に対応するスペクトルを含む一定の帯域内のスペクトルを抽出する。さらに、本願発明では、フーリエ変換後に得られるスペクトルのうち、フレーム画像から得られる呼吸要素以外の周波数(例えば、また、各部位の「density」/「intensity」、心拍または血管拍から得られる心拍要素)、またはオペレータによって外部から入力された周波数に対応するスペクトル(例えば、スペクトルモデル)を含む一定の帯域内のスペクトルを抽出することも可能である。
ここで、合成波のスペクトルの成分は2つの成分(呼吸、血流)のみであれば50%+50%となり、3つの成分の場合は、1/3ずつの配分となる。このため、呼吸成分のスペクトルが何%、血流成分のスペクトルが何%と、スペクトルの成分およびその高さからある程度、合成波のスペクトルを計算することができる。その割合(%)が高いところでスペクトルを抽出することが可能である。すなわち、血流成分/呼吸成分と合成波成分との割合を計算し、血流成分/呼吸成分の高いスペクトル値を計算して抽出する。なお、横隔膜の同定等において、呼吸や心臓血管の周波数を取得したdataから、Hz(周波数)が比較的一定になる部位、すなわち、Hzの変化が少ない領域に対応するスペクトルやその重ね合わせのみを抽出する場合もある。また、スペクトルの帯域を定める場合、横隔膜の同定等をする際に、Hzの変化が生じたrangeおよびその周囲の領域でスペクトルの帯域を定める場合もある。波形の構成要素を加味することもある。
なお、逆フーリエ変換する際のスペクトルについては、「単純モデル化した周波数および周波数帯から高い部位(一つもしくは複数)を用いて抽出する場合(シュミレーション主義)」と、「実際の周波数または周波数帯に基づいて、スペクトル値に応じて周波数の高い部位または周波数の低い部位を抽出する場合(現場主義)」とを選択することが可能である。また、心臓の周波数がA、肺の周波数がBである場合、周波数帯域全体からAを減算することによって、Bが得られる。また、フーリエ変換から取得されるスペクトルについては、周波数軸上の一か所のみならず複数の箇所を抽出することも可能である。
以上により、呼吸要素の周期や血管拍周期に完全一致する場合のみならず、考慮した方が良いスペクトルも抽出することができ、画像診断に寄与することが可能となる。なお、「呼吸」や「心拍」は、特定の周波数帯域に含まれることが知られている。このため、呼吸の場合は、例えば、「0~0.5Hz(呼吸数0~30回/分)」、循環器の場合は、例えば、「0.6~2.5(心拍/脈拍数36~150回/分)Hz」というフィルタを用いて、予めこのフィルタで呼吸周波数や循環器の周波数を特定していくことも可能である。これにより、周波数同調性画像を表示することが可能となる。これは、心臓の「density」変化を取得する際に、呼吸(肺)の「density」変化を拾ってしまったり、肺の「density」変化を取得する際に、心臓の「density」変化を拾ってしまったりする場合があるからである。
[視覚化・数値化]
上記のように解析した結果を、視覚化および数値化する。視覚化および数値化をする際に、本明細書では、「モデル化した肺」を定義する。肺を動画像で表示する際、位置関係が動いてしまうため、相対的判断が容易ではない。このため、位置関係のずれを、空間的に統一化・平均化する。例えば、肺の形状を扇形などの図形に当てはめ、形を整えた状態で表示する。そして、リコンストラクションの概念を用いて時間的に統一化する。例えば、「複数の呼吸のうち、20%の肺の状況」を抽出し、それを「一呼吸の20%の肺の状況」として定めることができる。このように、空間的、時間的に統一化した肺を「モデル化した肺」とする。これにより、異なる患者同士を比較したり、一人の患者の現在と過去とを比較したりする際に、相対的判断が容易となる。
例えば、standard uptakeとして、計測された肺野全域の「density」/「intensity」から平均値を1として相対的/対数的に値を表示することがある。また、血流の方向だけを採用するため、特定方向への変化を切り出すことがある。これにより、意味のある方法のデータだけを取り出すことが可能となる。肺野同定結果を用いて、解析範囲の変化に追従して疑似カラー化を行なう。すなわち、フェーズに合わせた特定の形(最小、最大、平均、中央値)に沿って、各個人(被写体)の解析結果を相対的な領域に当てはめる。
また、複数の解析結果を比較できる特定の形状・フェーズに変形させる。さらに、モデル化した肺を作成する際、上記呼吸要素の周期解析の結果を用いて、肺野内の相対的な位置関係を計算する。なお、モデル化した肺は、複数の患者の胸郭ライン、「density」、横隔膜などを総合的に平均化したラインを用いて作成する。モデル化した肺の作成の際、肺血流の場合は、肺門から肺端部にかけて放射状に距離を測ることができる。また、呼吸の場合は、胸郭や横隔膜の動きに応じて補正する必要がある。さらに、肺尖からの距離を考慮して複合的に計算しても良い。
また、逆フーリエ変換後に、振幅値が相対的に大きいブロックのみを抽出して表示しても良い。すなわち、ブロック毎にフーリエ解析を行なう場合、逆フーリエ変換の後、波の振幅が大きいブロックと、波の振幅が小さいブロックが存在する。そこで、振幅が相対的に大きいブロックのみを抽出し、視覚化することも有効である。また、逆フーリエ変換後、各数値の実部と虚部とをそれぞれ使い分けることができる。例えば、実部のみから画像を再構成したり、虚部のみから画像を再構成したり、実部および虚部の絶対値から画像を再構成することが可能である。
モデル化した肺に対して、フーリエ解析をしても良い。すなわち、呼吸数回の画像を合わせたり、フーリエ解析や相対位置把握をしたりする際にもモデル化した肺を用いることが可能である。モデル化した肺を用いることによって、取得した複数のフレームをモデル化した肺に当てはめたり、血管の場合、心拍(例えば、肺門部から得られる心拍など)に応じて計算したモデル化した肺に当てはめたりすることによって、フーリエ解析を行なう場合の相対位置を一定とすることが可能となる。ベースとなる呼吸の状態を取得する際に、モデル化した肺を使うことによって、安定的に計算結果を得ることも可能となる。また、肺をモデル化することで、空間の相違を固定化することができ、肺の動きを見やすくさせることが可能となる。
画像化において、相対評価の標識方法は、以下の通りである。すなわち、画像を相対的に白黒、カラーマッピングで標識する。差分によって得られた「density」/「intensity」の数%前後の値をカットし、その上下残りを相対的に表示することがある。もしくは、得られた差分の前後数%前後の値はとびぬけた値となる場合があるため、これを「artifact」として除去し、残りの部分を相対的に表示することがある。0~255諧調などの方法の他、0~100%の値として表示することもある。
なお、ピクセルをある程度曖昧に表示し、ぼやけた状態にして全体を表示することも可能である。特に、肺血管の場合、高い信号値の間に低い信号の値が混在するが、高い信号値のみを大雑把に把握することができれば、全体として曖昧であっても構わない。例えば、血流の場合は、閾値以上の信号を抜き出し、呼吸の場合は閾値以上の信号を抜き出さなかったりしてもよい。具体的には、次の表の数字を1ピクセルとして真ん中の数値を取得する場合、真ん中の数値が占める割合を取得し、1ピクセル内で平均化すると、隣接する画素との間で滑らかに表現することができる。この手法は、ブロックごとの平均インテンシティ算出の際にも用いることが可能である。
Figure 0007169514000001
この手法は、肺野のみならず、任意の解析範囲のデンシティ(density)を検出し、デンシティが相対的に大きく変化する箇所を除去する際にも適用することが可能である。また、予め設定した閾値を、大きく超える点をカットオフする。また、肋骨の形態認識、例えば、突然出現する高/低信号ラインを認識し、除去する。また、位相から同様に、突然出現する信号、例えば、reconstructionのphaseが15%~20%前後でartifactが認められる患者の特徴など、通常の波の変化と異なる突然の信号を除去する場合がある。なお、最初にベースデータを取る際、(横隔膜)≒(胸郭)≒(胸郭の動き)≒(スパイロメータ)≒(肺野)、fieldの(density)≒(volumetry)などの計算において、位相が異なる場合があり、その位相を実際認識できる形態(XPの輪郭)に当てはめる場合がある。
モデル化した肺が作成できると、上述したように、同調性、一致率、不一致率を数値化して提示することが可能となる(周波数同調性画像または波長同調性画像の表示)。これにより、正常な状態からの逸脱が表示できる。本実施形態によれば、フーリエ解析を実行することにより、新しい病気の可能性の発見、普通の自分との比較、手と足との比較や、反対側の手および足との比較が可能となる。さらに、足の動かし方、嚥下などでどこがおかしいのかを同調性の数値化で把握することが可能となる。また、病気の状態の人が一定時間経過後に変化しているかどうかを判断したり、変化している場合は、変化の前後を比較することが可能となる。また、肺野を、末梢からの距離を一定として放射状に見やすいような形態(円形~類円形)とすることで、内層~中層、外層などの評価をしやすくすることもでき、また、「血管の末梢優位」か「中層優位」かに応じて表現することも可能である。
なお、視覚化の際に、フーリエ変換後の画像とフーリエ変換前の画像を切り替えて表示したり、両者を一つの画面に並べて表示したりすることも可能である。
図2Dに示すように、モデル化した肺を100としたときに、当該人体において、何パーセントの相違があるのかどうかを把握し、変化率を表示することが可能となる。なお、肺全体のみならず、肺の一部分でも相違を把握することが可能である。特に、上述したように、横隔膜の動きのみを特定すると共に、横隔膜以外の肺野の形状を固定して、横隔膜の動きを表示すると共に、同調一致率や変化率を表示することが可能である。さらに、肺野の全部を固定して、同調一致率や変化率を表示することも可能である。なお、「Variation分類」をすることによって、標準血流の特定も可能である。すなわち、呼吸要素の周期を特定し、血管の相対的な位置関係を計算し、被写体の血流動態を標準血流として特定することが可能となる。
また、パターンマッチングの手法を用いて肺の検出をしても良い。図2E~図2Hは、肺野領域のパターン画像の例を示す図である。図2E~図2Hに示すように、肺の形を、パターン分類をしておき、これらのなかで近いものを抽出するようにしても良い。この手法により、対象の画像が、片肺を表すのか、両肺を表すのかを特定することができる。また、右肺であるのか、左肺であるのかを特定することもできる。パターン数は、限定しないが、4~5パターンを持っておくことが想定される。なお、このように、肺野の形態(形状)のみで右肺、左肺、両肺の認識をする方法もある。さらに、椎体・縦隔による太い帯状の“透過性低下部位”を認識し、その帯状の透過性低下部位との位置関係と、肺野の“透過性亢進部位”との位置関係から左右もしくは両肺を認識する方法を採ることも可能である。また、図2Hに示すように、横隔膜の下側の領域についてもこの方法を適用することができる。これにより、横隔膜の下側の部位や、心臓を認識することも可能である。
さらに、空気は、透過性が一番高く、肺野よりも透過性が高い部位であることから、空気も考慮して計算することが望ましい。すなわち、画面上の空気の位置によって、以下のように判断することができる。
(画面の右上の空気の領域)>(画面の左上の空気の領域)である場合は、左肺と認識する。これは、肩周りは人体外の空気の領域が撮影上広くなるからである。
(画面の左上の空気の領域)>(画面の右上の空気の領域)である場合は、右肺と認識する。これも、上記と同様に、肩周りは人体外の空気の領域が撮影上広くなるからである。
次に、(画面の右上の空気の領域)≒(画面の左上の空気の領域)である場合は、両肺と認識する。これは、空気の領域が左右同程度であるからである。
なお、横隔膜下に腸管の空気が入ることがあり、その際は認識されなくなる場合がある。このため、肺野の中心部から、縦隔側、心臓側、横隔膜側など、最初に大まかな肺野およびその周りの透過性低下部位を認識し、そのラインを肺野の淵を認識することもできる。この手法は、例えば
「https://jp.mathworks.com/help/images/examples/block-processing-large-images_ja_JP.html」に開示されている技術を用いることも可能である。
これにより、ある患者と他の患者との比較や数値化が可能となる。また、正常肺または正常血管と、典型的な異常肺機能または異常血流との比較や数値化が可能となる。さらに、ある患者の異なる時間における肺機能や肺血流の相対評価として、モデル化した肺および標準血流を使用することが可能となる。このようなモデル化した肺および標準血流は、種々のタイプの典型患者、健康な人の典型例を集合させ、モデル化した肺および標準血流とし、ある患者に形態的に当てはめて評価する際の指標として用いることが可能である。
[肺野の描画]
一般的に、肺野には透過性の低い肋骨が含まれるため、「density」のみを指標として肺の輪郭を機械的に同定することは難しい。そこで、本明細書では、ベジエ曲線および直線の組み合わせを用いて肺野の輪郭を仮に描画し、合致性が高くなるように、肺輪郭を調整する手法を採用する。
例えば、左肺の輪郭を4本のベジェ曲線と1本の直線で表現すると、肺輪郭上の5点と、制御点4点を求めることで、肺輪郭を描画することが可能になる。点の位置をずらして、複数の肺輪郭を描画し、“輪郭内の「density」の合計値が最大になる”、“輪郭線の内側と外側の数ピクセルの「density」合計の差分が最大になる”等の条件を用いて合致性を評価することで、肺輪郭を精度高く検出することができるようになる。実際には、比較的エッジが検出しやすい肺の上部の輪郭や、後述の方法で検出した横隔膜の位置から数点の位置を同定することも可能であり、上述のシミュレーションの試行回数を抑えることができる。古典的な二値化による輪郭抽出により、外縁に近い点を抽出し、最小二乗法等を利用して、ベジェ曲線の制御点位置を調整することも可能である。
図3Aおよび図3Bは、肺野の輪郭を、ベジエ曲線および直線の両方を用いて描画した例を示す図である。図3Aは肺の面積が最大となる場合(極大輪郭)を示し、図3Bは肺の面積が最小となる場合(極小輪郭)を示す。各図において、「cp1~cp5」は制御点を示し、「p1~p5」は、ベジエ曲線上または直線上の点を示す。このように、極大輪郭と極小輪郭が把握できると、途中の輪郭を計算によって求めることが可能となる。例えば、呼気の10%、20%…の状態を表示することが可能となる。このように、本実施形態によれば、少なくとも一つ以上のベジエ曲線(Bezier curve)を用いて、少なくとも肺野、血管または心臓を描画することが可能となる。なお、以上の手法は、肺に限定されるわけではなく、「臓器の検出」として、他の臓器にも適用可能である。また、例えば、特定のフレームにおいて予め定められた解析範囲(腫瘍、脳の視床下部、基底核、内包の境界など)上に、少なくとも一つ以上のベジエ曲線(Bezier curve)を用いて、他のフレームにおいて解析範囲に対応する範囲を検出する処理を実行することも可能である。
また、平面的な画像のみならず、立体的な画像(3D画像)についても適用可能である。曲面の方程式を定義し、その制御点を設定することによって、複数の曲面で囲まれた範囲を解析対象とすることが可能となる。
[横隔膜または呼吸と連動する動的部位の動きの検出]
連続撮影された画像において、横隔膜または呼吸と連動する動的部位の動きを検出することが可能である。連続撮影された画像において、任意の間隔で画像を選択し、画像間の差分を計算すると、特にコントラストの大きい領域について、差分が大きくなる。この差分を適切に可視化することによって、動きのあった領域を検出することができる。可視化の際には、閾値によるノイズ除去や、最小二乗法などを活用したカーブフィッティング等で差分の絶対値が大きいエリアの連続性を強調することもできる。
肺野においては横隔膜や心臓が接するラインのコントラストが際立っており、図4Aに示すように、2枚の肺画像において差分を取り、一定の閾値を設定して差分を可視化すると、図4Bに示すように、横隔膜や心臓が接するラインを可視化することができる。
[横隔膜の動きの推定]
本手法では、対象画像間において、横隔膜が動いている場合は横隔膜位置を検出可能であるが、横隔膜の動きが緩やかになる箇所の検出は困難となる。すなわち、呼気吸気が切り替わるタイミングや、呼吸を止めている間、撮影の開始直後や終了直前では検出が難しい。本手法においては、任意の補完方法を用いて、横隔膜の動きを推定する。
前述の方法を用いて、図4Bに示すように横隔膜ラインを可視化した後、縦1024pxの画像を縦8pxごとに128個の長方形に分割し、各長方形領域に含まれる信号値を合計し、図4Cに示すように棒グラフ化した。複数のピークのうち、点線の矩形で示す、最も下の座標に来るピークが横隔膜の位置を示していると期待される。通常の立位XP画像において、横隔膜は曲線として表示されるが、この座標を横隔膜の位置と近似する。
本手法で全画像について横隔膜位置を検出すると、図5に示すように「ピーク位置」が検出された。この検出した値に対して補正を行なうことで、横隔膜の動きを推定する。まず、差分が一定値よりも大きい場合は外れ値とみなして除外する(図5中の細い実線)。外れ値を除外したデータを、任意のクラスターに分割し、それぞれのクラスターに対して4次の曲線回帰を行ない、結果を繋ぎ合わせた(図5中の太い実線)。本解析では回帰分析を行なったが、本発明はこれに限定されるわけではなく、スプライン補間など任意の補完方法を用いる事が可能である。
[動的部位検出の精緻化]
動的部位のコントラストはラインに沿って一様でない場合がある。その場合はノイズ除去に使用する閾値を変更して、複数回検出処理を行なうことによって、動的部位の形状をより正確に検出することができる。例えば、左肺において、横隔膜のラインのコントラストは、人体内部にいくに従って弱くなる傾向がある。図4Bにおいては、横隔膜の右半分しか検出できていない。このとき、ノイズ除去に利用した閾値の設定を変えることによって、横隔膜の左半分の残りの部分を検出することもできる。この処理を複数回繰り返すことによって、横隔膜全体の形状を検出することが可能となる。本手法によって、横隔膜の位置だけでなく、形状について線や面の変化率や変化量を数値化することも可能となり、新たな診断に役立てることができる。
このように検出された横隔膜の位置または形状を診断に利用することが可能となる。すなわち、本願発明では、横隔膜の座標をグラフ化し、上述したように計算された曲線(局面)、若しくは直線を用いて、胸郭や横隔膜の座標の計算をし、また、心拍や血管拍、肺野の「density」などを、周期に対応した位置、座標としてグラフ化したりすることが可能である。このような手法は、呼吸と連動する動的部位についても適用可能である。
このような手法により、吸気、呼気でのHzだけでなく、横隔膜または呼吸と連動する動的部位の周波数(Hz)が変化した場合、その変化に応じた周波数帯域で計測できるようになる。そして、BPF(band pass filter)のスペクトル抽出の際に、一定の範囲において、呼吸それぞれの状態に応じてBBFを据えること、呼吸それぞれの「reconstruction phase」でBPFの位置の軸が変動し、最適な状態が生じえること、それを合わせた変動性のBPFを作成することが可能となる。これにより、呼吸が遅くなったり、止めたとき(Hz=0)のように、呼吸のリズムの変動があっても、それに応じた画像を提供することが可能となる。
また、呼吸要素が呼気または吸気の全体に占める割合に基づいて、呼気または吸気の全体の周波数を計算するようにしても良い。なお、横隔膜の検出において、複数回施行し、信号または波形が安定しているものを選択するようにしても良い。以上により、検出した横隔膜の位置若しくは形状、または呼吸と連動する動的部位の位置若しくは形状、から、呼吸要素の少なくとも一つの周波数を計算することが可能となる。横隔膜または動的部位の位置または形状が把握できると、呼吸要素の周波数を把握することが可能となる。この手法によれば、波形の一部を区切ったとしても、その後の波形を追跡することができる。このため、呼吸要素の周波数が途中で変わっても、元々の呼吸要素を追従することが可能である。また、心臓の拍動などが突然変わることがあるが、心血管についても同様に適用することが可能となる。次に、本実施形態に係る各モジュールの動作について説明する。
[呼吸機能解析]
まず、呼吸機能解析について説明する。図6Aは、本実施形態に係る呼吸機能解析の概要を示すフローチャートである。基本モジュール1がデータベース15からDICOMの画像を抽出する(ステップS1)。ここでは、少なくとも、一呼吸周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像において、少なくとも肺野内のある一定領域における密度(density/intensity)を用いて、呼吸要素の周期を特定する(ステップS2)。なお、特定した呼吸周期やこの呼吸周期から特定される波形については、以下の各ステップで用いることが可能である。
呼吸要素の周期の特定は、さらに、横隔膜の動き、胸郭の動きを用いることも可能である。また、X線の透過性が高い部位で測定される、ある一定のvolume、「density」/「intensity」で構成される範囲、スパイログラムなどの他の測定方法から得られるデータを用いても良い。なお、予め各臓器(ここでは肺)が有する周波数を特定しておき、その特定した周波数に対応する「density」/「intensity」を抽出しても良い。
次に、図6Aにおいて、肺野を自動検出する(ステップS3)。肺輪郭は連続的に変化するため、最大の形状と最小の形状が検出できれば、間の形状は計算により補間が可能である。ステップS2において特定した呼吸要素の周期に基づいて、各フレーム画像を補間することによって、各フレーム画像における肺輪郭を特定する。また、図2E~図2Hで示したようなパターンマッチングを行なって、肺野を検出しても良い。なお、検出した肺野について、カットオフによるノイズ除去を行なっても良い。次に、検出した肺野を複数のブロック領域に分割する(ステップS4)。そして、各フレーム画像における各ブロック領域の変化を計算する(ステップS5)。ここでは、各ブロック領域内での変化の値を平均化し、1つのデータとして表現する。
なお、各ブロック領域内での変化の値について、カットオフによるノイズ除去を行なっても良い。次に、各ブロック領域の「density」/「intensity」の値や、また、その変化量について、上記呼吸要素の周期に基づいて、フーリエ解析または同調一致率の解析を実施する(ステップS6)。
次に、フーリエ解析または同調一致率の解析により得られた結果について、ノイズ除去を行なう(ステップS7)。ここでは、上述したようなカットオフや、アーティファクト(artifact)の除去を行なうことができる。以上のステップS5からステップS7の動作を1回以上行ない、完了するかどうかを判断する(ステップS8)。ここで、ディスプレイで表示される特徴量については、合成波や他の波の混在により純度の高い要素、例えば、呼吸要素や血流要素、その他の要素の周波数同調性画像が一度のスペクトル抽出では表示できない場合もある。その際は、ディスプレイで表示される特徴量を画素値として、再度、複数回ディスプレイに至るすべてもしくは一部の解析をし直す場合がある。この作業によってさらに要素、例えば、呼吸要素や血流要素の同調性や一致性に関し純度の高い画像を取得することが可能となる。この操作については、オペレータがディスプレイの画像を視認しながら手動で行なっても良いし、出力結果からスペクトルを抽出してその分布割合を計算し直すことを自動的に行なっても良い。さらに、計算後においても、場合に応じて、ノイズカット処理、最小二乗法による穴埋め(補間)、周囲の「density」を用いた補正を行なっても良い。
ステップS8において、完了しない場合は、ステップS5に遷移し、完了する場合は、フーリエ解析または同調一致率の解析により得られた結果を、疑似カラー画像としてディスプレイに表示する(ステップS9)。なお、白黒画像を表示しても良い。このように、複数のサイクルを繰り返すことによって、データの確度を高める場合もある。これにより、所望の動画を表示することが可能となる。また、ディスプレイに表示された画像を修正することで、所望の動画を得るようにしても良い。
本実施形態では、所望の周波数または周波数帯域を計算により算出するが、実際の画像として見ると、必ずしも良い画像が表示できるとは限らない。このため、以下の手法を採る場合もある。
(1)いくつかの周波数帯を複数提示し、人的に選択する方法
(2)いくつかの周波数帯を複数提示し、AI技術によりパターン認識でよい画像を抽出する方法
(3)HISTGRAMの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、HISTGRAMの傾向、形態から選択しても良い。
[肺血流解析]
次に、肺血流解析について説明する。図7は、本実施形態に係る肺血流解析の概要を示すフローチャートである。基本モジュール1がデータベース15からDICOMの画像を抽出する(ステップT1)。ここでは、少なくとも、一心拍周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像に基づいて、血管拍周期を特定する(ステップT2)。なお、特定した血管拍周期やこの血管拍周期から特定される波形については、以下の各ステップで用いることが可能である。血管拍周期は、上述したように、例えば、心電図や脈拍計等の他のモダリティの計測結果、心臓・肺門・主要血管など、任意の部位の「density」/「intensity」の変化を用いて血管拍を解析する。なお、予め各臓器(ここでは肺血流)が有する周波数を特定しておき、その特定した周波数に対応する「density」/「intensity」を抽出しても良い。
次に、図7において、上述した方法で呼吸要素の周期を特定し(ステップT3)、その呼吸要素の周期を用いて肺野を自動検出する(ステップT4)。肺の輪郭の自動検出では、フレーム画像毎に、ばらつきが生ずることもあるが、ステップT3において特定した呼吸要素の周期に基づいて、各フレーム画像を補間することによって、各フレーム画像における肺輪郭を特定する。また、図2E~図2Hで示したようなパターンマッチングを行なって、肺野を検出しても良い。なお、検出した肺野について、カットオフによるノイズ除去を行なっても良い。次に、検出した肺野を複数のブロック領域に分割する(ステップT5)。そして、各フレーム画像における各ブロック領域の変化を計算する(ステップT6)。ここでは、各ブロック領域内での変化の値を平均化し、1つのデータとして表現する。なお、各ブロック領域内での変化の値について、カットオフによるノイズ除去を行なっても良い。次に、各ブロック領域の「density」/「intensity」の値や、また、その変化量について、上記血管拍周期に基づいて、フーリエ解析または同調一致率の解析を実施する(ステップT7)。
次に、フーリエ解析または同調一致率の解析により得られた結果について、ノイズ除去を行なう(ステップT8)。ここでは、上述したようなカットオフや、アーティファクト(artifact)の除去を行なうことができる。以上のステップT6からステップT8の動作を1回以上行ない、完了するかどうかを判断する(ステップT9)。ここで、ディスプレイで表示される特徴量については、合成波や他の波の混在により純度の高い要素、例えば、呼吸要素や血流要素、その他の要素の周波数同調性画像が一度のスペクトル抽出では表示できない場合もある。その際は、ディスプレイで表示される特徴量を画素値として、再度、複数回ディスプレイに至るすべてもしくは一部の解析をし直す場合がある。この作業によってさらに要素、例えば、呼吸要素や血流要素の同調性や一致性に関し純度の高い画像を取得することが可能となる。この操作については、オペレータがディスプレイの画像を視認しながら手動で行なっても良いし、出力結果からスペクトルを抽出してその分布割合を計算し直すことを自動的に行なっても良い。さらに、計算後においても、場合に応じて、ノイズカット処理、最小二乗法による穴埋め(補間)、周囲の「density」を用いた補正を行なっても良い。
ステップT9において、完了しない場合は、ステップT6に遷移し、完了する場合は、フーリエ解析または同調一致率の解析により得られた結果を、疑似カラー画像としてディスプレイに表示する(ステップT10)。なお、白黒画像を表示しても良い。これにより、データの確度を高めることが可能となる。また、ディスプレイに表示された画像を修正することで、所望の動画を得るようにしても良い。
本実施形態では、所望の周波数または周波数帯域を計算により算出するが、実際の画像として見ると、必ずしも良い画像が表示できるとは限らない。このため、以下の手法を採る場合もある。
(1)いくつかの周波数帯を複数提示し、人的に選択する方法
(2)いくつかの周波数帯を複数提示し、AI技術によりパターン認識でよい画像を抽出する方法
(3)HISTGRAMの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、HISTGRAMの傾向、形態から選択しても良い。
[その他の血流解析]
次に、その他の血流解析について説明する。本発明の一態様は、図15に示すように、心臓、大動脈、肺血管、上腕動脈、頸部血管などの血流解析についても適用可能である。さらに、図示しない腹部血管や、末梢の血管などについても、同様に血流解析が可能である。図8は、本実施形態に係るその他の血流解析の概要を示すフローチャートである。基本モジュール1がデータベース15からDICOMの画像を抽出する(ステップR1)。ここでは、少なくとも、一心拍周期内に含まれる複数のフレーム画像を取得する。次に、取得した各フレーム画像に基づいて、血管拍周期を特定する(ステップR2)。なお、特定した血管拍周期やこの血管拍周期から特定される波形については、以下の各ステップで用いることが可能である。血管拍周期は、上述したように、例えば、心電図や脈拍計等の他のモダリティの計測結果、心臓・肺門・主要血管など、任意の部位の「density」/「intensity」の変化を用いて血管拍を解析する。なお、予め各臓器(例えば、主要血管)が有する周波数を特定しておき、その特定した周波数に対応する「density」/「intensity」を抽出しても良い。
次に、解析範囲を設定し(ステップR3)、設定した解析範囲を複数のブロック領域に分割する(ステップR4)。そして、各ブロック領域内での変化の値を平均化し、1つのデータとして表現する。なお、各ブロック領域内での変化の値について、カットオフによるノイズ除去を行なっても良い。次に、各ブロック領域の「density」/「intensity」の値や、また、その変化量について、上記血管拍周期に基づいて、フーリエ解析または同調一致率の解析を実施する(ステップR5)。
次に、フーリエ解析または同調一致率の解析により得られた結果について、ノイズ除去を行なう(ステップR6)。ここでは、上述したようなカットオフや、アーティファクト(artifact)の除去を行なうことができる。以上のステップR5からステップR6の動作を1回以上行ない、完了するかどうかを判断する(ステップR7)。ここで、ディスプレイで表示される特徴量については、合成波や他の波の混在により純度の高い要素、例えば、呼吸要素や血流要素、その他の要素の周波数同調性画像が一度のスペクトル抽出では表示できない場合もある。その際は、ディスプレイで表示される特徴量を画素値として、再度、複数回ディスプレイに至るすべてもしくは一部の解析をし直す場合がある。この作業によってさらに要素、例えば、呼吸要素や血流要素の同調性や一致性に関し純度の高い画像を取得することが可能となる。この操作については、オペレータがディスプレイの画像を視認しながら手動で行なっても良いし、出力結果からスペクトルを抽出してその分布割合を計算し直すことを自動的に行なっても良い。さらに、計算後においても、場合に応じて、ノイズカット処理、最小二乗法による穴埋め(補間)、周囲の「density」を用いた補正を行なっても良い。
ステップR7において、完了しない場合は、ステップR5に遷移し、完了する場合は、フーリエ解析または同調一致率の解析により得られた結果を、疑似カラー画像としてディスプレイに表示する(ステップR8)。なお、白黒画像を表示しても良い。これにより、データの確度を高めることが可能となる。また、ディスプレイに表示された画像を修正することで、所望の動画を得るようにしても良い。
本実施形態では、所望の周波数または周波数帯域を計算により算出するが、実際の画像として見ると、必ずしも良い画像が表示できるとは限らない。このため、以下の手法を採る場合もある。
(1)いくつかの周波数帯を複数提示し、人的に選択する方法
(2)いくつかの周波数帯を複数提示し、AI技術によりパターン認識でよい画像を抽出する方法
(3)HISTGRAMの傾向、形態から選択する。すなわち、結果の信号における「Histgram」の中心部の値が高くなる傾向があり、また、動きに応じて「histgram」の値が変動するため、HISTGRAMの傾向、形態から選択しても良い。
なお、3Dで解析した場合、呼吸量、心拍出量、中枢の血流量を別の装置で測定することによって、相対的な値であるフーリエ解析結果から各ブロック領域の呼吸量、心拍出量、中枢の血流量を算出することが可能となる。すなわち、呼吸機能解析の場合は、呼吸量から肺換気量の推定が可能となり、肺血流解析の場合は、心(肺血管)拍出量から肺血流量の推定が可能となり、その他の血流量解析の場合は、中枢側の血流量(割合)から描出される分岐血管における推定血流量(割合)の推定が可能となる。
また、上記のように、取得したdatabaseは、すべてを計算できればより精度の高い判断が可能であるが、コンピュータ解析を実行する上では時間を要する場合がある。そのため、任意の枚数(例えば、特定のphase)だけを抜き出し、計算を行なうようにすることもできる。これにより、解析時間を短縮することができ、さらに、呼吸の初めに観られるようなイレギュラーな個所をカットすることが可能となる。また、解析結果を表示する際に、任意の範囲を表示することができる。例えば、「呼気/吸気」の変わり目から「吸気/呼気」の変わり目の範囲を表示することによって、繰り返し再生をする際に、いわゆる「エンドレス再生」が可能となり、医師による診断をし易くすることが可能となる。
以上説明したように、本実施形態によれば、X線動画装置で人体の画像を評価することが可能となる。デジタルデータを取得できれば、現存施設装置で概ね良好に計算可能であり、導入費用が安価で済む。例えば、Flat panel detectorを用いたX線動画装置において、被写体の検査を簡便に済ませることが可能となる。また、肺血流については、肺血栓塞栓症のスクリーニングが可能となる。例えば、Flat panel detectorを用いたX線動画装置においては、CTを行なう前に本実施形態に係る診断支援プログラムを実行することによって、無駄な検査を除外することができる。また、検査が簡便であるため、緊急性の高い疾患を早期に発見し、優先的に対応することが可能となる。なお、現時点における撮影方法では、CT、MRIなどの他のmodalityでは、いくつかの課題があるが、これを解決できれば、各領域の細かい診断が可能となる。
また、各種の血管、例えば、頸部血流狭小化のスクリーニングにも適用可能であり、また、大血管の血流評価やスクリーニングにも適用可能である。また、肺呼吸データについては、肺の部分機能検査として有効となり、肺機能検査として使用することが可能となる。また、COPD、肺気腫などの疾患の同定も可能となる。さらに、術前、術後の性状の把握にも適用可能である。さらに、呼吸要素の周期および血流周期をフーリエ解析し、腹部のX線画像において、呼吸の波形および血流の波形を除去することで、残りの生体運動の異変、例えば、腸管イレウスなどが観察可能となる。
なお、最初に取得した画像が、ある程度高精細である場合は、画素数が多いため、計算時間に時間がかかることもある。その場合、一定のピクセル数に画像を減らしてから計算しても良い。例えば、「4096×4096」のピクセルを実際には「1024×1024」にしてから計算することで計算時間を抑えることが可能である。
[その他]
なお、X線画像を撮影する際に、例えば、AR法(Autoregressive Moving average model)などの予測アルゴリズムを用いることができる。呼吸要素の少なくとも一つの周波数が特定できると、この周波数に応じて、X線の照射間隔を調整するよう、X線撮影装置を制御することも可能である。例えば、呼吸要素の周波数が小さい場合(周期が長い場合)、X線撮影の回数を減らすことができる。これにより、人体の被ばく量を減らすことが可能となる。なお、頻呼吸や頻脈など、呼吸要素や心血管拍要素の周波数が大きい場合(周期が短い場合)は、照射頻度を高めて最適な画像作成を行なっても良い。
また、DICOMデータの保存形式であるが、圧縮をすると画像の質が低下する場合があることから、非圧縮にて保存することが望ましい。また、データの圧縮形式に応じて、計算方法を変えるようにしても良い。
1 基本モジュール
3 呼吸機能解析部
5 肺血流解析部
7 その他の血流解析部
9 フーリエ解析部
10 波形解析部
11 視覚化・数値化部
13 入力インタフェース
15 データベース
17 出力インタフェース
19 ディスプレイ

Claims (4)

  1. 人体の画像を解析し、解析結果を表示する診断支援プログラムであって、
    複数のフレーム画像を取得する処理と、
    前記各フレーム画像内の画素値を相対値または対数値として算出する処理と、
    前記相対値若しくは対数値として算出されたいずれかのフレーム画像と前記相対値若しくは対数値として算出され任意の間隔で選択された他のフレーム画像との比を算出する処理と、
    前記算出された比を可視化する処理と、
    を備え
    肺野の濃度変化または血管の位置関係に基づいて、一つ以上のフレーム画像の肺野領域の画素値を算出することを特徴とする診断支援プログラム。
  2. 人体の画像を解析し、解析結果を表示する診断支援プログラムであって、
    複数のフレーム画像を取得する処理と、
    前記各フレーム画像内の画素値を相対値または対数値として算出する処理と、
    前記相対値若しくは対数値として算出されたいずれかのフレーム画像と前記相対値若しくは対数値として算出され任意の間隔で選択された他のフレーム画像との比を算出する処理と、
    前記算出された比を可視化する処理と、
    を備え、
    アーティファクトを除外したフレーム画像の画素値を算出することを特徴とする診断支援プログラム。
  3. 心室容量の波形に基づいて、臓器の動的な画像を表示することを特徴とする請求項1または請求項2に記載の診断支援プログラム。
  4. 肺野を複数のブロック領域に分割し、各フレーム画像におけるブロック領域の平均画素値を算出することを特徴とする請求項1から請求項3のいずれかに記載の診断支援プログラム。
JP2022066013A 2018-01-05 2022-04-13 診断支援プログラム Active JP7169514B2 (ja)

Applications Claiming Priority (6)

Application Number Priority Date Filing Date Title
JP2018000855 2018-01-05
JP2018000855 2018-01-05
JP2018043511 2018-03-09
JP2018043511 2018-03-09
PCT/JP2019/000102 WO2019135412A1 (ja) 2018-01-05 2019-01-07 診断支援プログラム
JP2019563986A JP7462898B2 (ja) 2018-01-05 2019-01-07 診断支援プログラム

Related Parent Applications (1)

Application Number Title Priority Date Filing Date
JP2019563986A Division JP7462898B2 (ja) 2018-01-05 2019-01-07 診断支援プログラム

Publications (2)

Publication Number Publication Date
JP2022095871A JP2022095871A (ja) 2022-06-28
JP7169514B2 true JP7169514B2 (ja) 2022-11-11

Family

ID=67144237

Family Applications (4)

Application Number Title Priority Date Filing Date
JP2019563986A Active JP7462898B2 (ja) 2018-01-05 2019-01-07 診断支援プログラム
JP2022066013A Active JP7169514B2 (ja) 2018-01-05 2022-04-13 診断支援プログラム
JP2022141459A Active JP7310048B2 (ja) 2018-01-05 2022-09-06 診断支援プログラム
JP2023126624A Pending JP2023153937A (ja) 2018-01-05 2023-08-02 診断支援プログラム

Family Applications Before (1)

Application Number Title Priority Date Filing Date
JP2019563986A Active JP7462898B2 (ja) 2018-01-05 2019-01-07 診断支援プログラム

Family Applications After (2)

Application Number Title Priority Date Filing Date
JP2022141459A Active JP7310048B2 (ja) 2018-01-05 2022-09-06 診断支援プログラム
JP2023126624A Pending JP2023153937A (ja) 2018-01-05 2023-08-02 診断支援プログラム

Country Status (13)

Country Link
US (1) US20210052228A1 (ja)
EP (1) EP3735906A4 (ja)
JP (4) JP7462898B2 (ja)
KR (1) KR20200106050A (ja)
CN (1) CN111867471A (ja)
AU (1) AU2019205878A1 (ja)
BR (1) BR112020013421A2 (ja)
CA (1) CA3087702A1 (ja)
MX (1) MX2020006984A (ja)
PH (1) PH12020551042A1 (ja)
SG (1) SG11202006059WA (ja)
WO (1) WO2019135412A1 (ja)
ZA (1) ZA202004575B (ja)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111867471A (zh) * 2018-01-05 2020-10-30 帕拉梅维亚私人有限公司 诊断辅助程序
JP7047574B2 (ja) * 2018-04-26 2022-04-05 コニカミノルタ株式会社 動態画像解析装置、動態画像解析システム、動態画像解析プログラム及び動態画像解析方法
JP2020141841A (ja) * 2019-03-06 2020-09-10 コニカミノルタ株式会社 動態解析装置及びプログラム
TWI684433B (zh) * 2019-04-19 2020-02-11 鉅怡智慧股份有限公司 生物影像處理方法以及生理資訊檢測裝置
TW202116253A (zh) * 2019-07-04 2021-05-01 新加坡商派拉梅維爾私人有限公司 診斷支援程式
CN114173665A (zh) * 2019-08-01 2022-03-11 国立大学法人大阪大学 血管壁厚度推测方法、血管壁厚度推测装置及血管壁厚度推测系统
JP7435242B2 (ja) * 2020-05-15 2024-02-21 コニカミノルタ株式会社 動態画像解析装置、動態画像解析方法及びプログラム
CN111879724B (zh) * 2020-08-05 2021-05-04 中国工程物理研究院流体物理研究所 基于近红外光谱成像的人皮面具识别方法及系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005151099A (ja) 2003-11-14 2005-06-09 Canon Inc 放射線画像処理装置及び処理方法
JP2008219654A (ja) 2007-03-06 2008-09-18 Canon Inc 画像処理装置及び画像処理方法
WO2012026145A1 (ja) 2010-08-27 2012-03-01 コニカミノルタエムジー株式会社 診断支援システム及びプログラム
JP2012110397A (ja) 2010-11-22 2012-06-14 Konica Minolta Medical & Graphic Inc 動態診断支援情報生成システム
JP2015043880A (ja) 2013-08-28 2015-03-12 コニカミノルタ株式会社 胸部診断支援システム
US20160256059A1 (en) 2015-03-04 2016-09-08 Dmytro Volkov Method and system for processing of medical images for generating a prognosis of cardiac function
JP2017113344A (ja) 2015-12-25 2017-06-29 コニカミノルタ株式会社 動態撮影装置及び動態撮影システム
JP2017124325A (ja) 2008-01-15 2017-07-20 コニカミノルタ株式会社 動態画像解析装置

Family Cites Families (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH07100060B2 (ja) * 1987-03-30 1995-11-01 コニカ株式会社 肺野部輪郭線を求める画像計測方法
JP3222273B2 (ja) * 1993-07-09 2001-10-22 株式会社日立製作所 核磁気共鳴診断装置における動画像の画質改善方法
US6577752B2 (en) 2001-06-15 2003-06-10 Arch Development Corporation Automated method and system for the delineation of the chest wall in computed tomography scans for the assessment of pleural disease
US20040260188A1 (en) * 2003-06-17 2004-12-23 The General Hospital Corporation Automated auscultation system
JP5402056B2 (ja) * 2009-02-16 2014-01-29 コニカミノルタ株式会社 画像処理装置、画像処理方法、およびプログラム
JP5446443B2 (ja) * 2009-05-15 2014-03-19 日産自動車株式会社 心拍数測定装置および心拍数測定方法
JP5804340B2 (ja) * 2010-06-10 2015-11-04 株式会社島津製作所 放射線画像領域抽出装置、放射線画像領域抽出プログラム、放射線撮影装置および放射線画像領域抽出方法
JP5556413B2 (ja) 2010-06-18 2014-07-23 コニカミノルタ株式会社 動態画像処理装置及びプログラム
FI20115053A0 (fi) * 2011-01-19 2011-01-19 Delfin Technologies Oy Menetelmä ja järjestelmä kardiovaskulaaristen sykeaaltojen visualisoimiseksi
JP5962237B2 (ja) * 2012-06-11 2016-08-03 コニカミノルタ株式会社 胸部診断支援情報生成方法
JP5907185B2 (ja) * 2014-01-27 2016-04-26 コニカミノルタ株式会社 動態画像診断支援システム
US20160189394A1 (en) * 2014-12-30 2016-06-30 Huazhong University Of Science And Technology Method for iteratively extracting motion parameters from angiography images
JP6348865B2 (ja) * 2015-03-30 2018-06-27 株式会社リガク Ct画像処理装置および方法
JP6638729B2 (ja) * 2015-07-22 2020-01-29 コニカミノルタ株式会社 コンソール及び動態画像撮影診断システム
JP2017131310A (ja) * 2016-01-26 2017-08-03 株式会社日立製作所 X線透視撮影装置、撮影方法
JP6701880B2 (ja) * 2016-03-30 2020-05-27 コニカミノルタ株式会社 動態解析装置、動態解析システム、動態解析方法及びプログラム
US9947093B2 (en) * 2016-05-03 2018-04-17 Konica Minolta, Inc. Dynamic analysis apparatus and dynamic analysis system
EP3488782B1 (en) * 2016-07-19 2021-09-29 Paramevia Pte. Ltd. Diagnosis support program
CN111867471A (zh) * 2018-01-05 2020-10-30 帕拉梅维亚私人有限公司 诊断辅助程序

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005151099A (ja) 2003-11-14 2005-06-09 Canon Inc 放射線画像処理装置及び処理方法
JP2008219654A (ja) 2007-03-06 2008-09-18 Canon Inc 画像処理装置及び画像処理方法
JP2017124325A (ja) 2008-01-15 2017-07-20 コニカミノルタ株式会社 動態画像解析装置
WO2012026145A1 (ja) 2010-08-27 2012-03-01 コニカミノルタエムジー株式会社 診断支援システム及びプログラム
JP2012110397A (ja) 2010-11-22 2012-06-14 Konica Minolta Medical & Graphic Inc 動態診断支援情報生成システム
JP2015043880A (ja) 2013-08-28 2015-03-12 コニカミノルタ株式会社 胸部診断支援システム
US20160256059A1 (en) 2015-03-04 2016-09-08 Dmytro Volkov Method and system for processing of medical images for generating a prognosis of cardiac function
JP2017113344A (ja) 2015-12-25 2017-06-29 コニカミノルタ株式会社 動態撮影装置及び動態撮影システム

Also Published As

Publication number Publication date
PH12020551042A1 (en) 2021-09-06
JP7310048B2 (ja) 2023-07-19
KR20200106050A (ko) 2020-09-10
BR112020013421A2 (pt) 2020-12-01
EP3735906A4 (en) 2021-11-03
JP7462898B2 (ja) 2024-04-08
CN111867471A (zh) 2020-10-30
JP2023153937A (ja) 2023-10-18
CA3087702A1 (en) 2019-07-11
US20210052228A1 (en) 2021-02-25
MX2020006984A (es) 2020-10-05
WO2019135412A1 (ja) 2019-07-11
JP2022095871A (ja) 2022-06-28
ZA202004575B (en) 2022-12-21
SG11202006059WA (en) 2020-07-29
JP2022172305A (ja) 2022-11-15
EP3735906A1 (en) 2020-11-11
JPWO2019135412A1 (ja) 2021-02-25
AU2019205878A1 (en) 2020-08-06

Similar Documents

Publication Publication Date Title
JP7169514B2 (ja) 診断支援プログラム
US20210233243A1 (en) Diagnostic support program
US11189025B2 (en) Dynamic image analysis apparatus, dynamic image analysis method, and recording medium
WO2023063318A1 (ja) 診断支援プログラム
TWI828661B (zh) 診斷支援系統
EP3995081A1 (en) Diagnosis assisting program
JP2020171475A (ja) 動態画像解析装置、動態画像解析方法及びプログラム
OA20419A (en) Diagnostic support program
EA040692B1 (ru) Программа поддержки диагностики

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20220428

A871 Explanation of circumstances concerning accelerated examination

Free format text: JAPANESE INTERMEDIATE CODE: A871

Effective date: 20220428

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20220518

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20220607

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20220808

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20220823

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20220830

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

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20221011

R150 Certificate of patent or registration of utility model

Ref document number: 7169514

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150