JP2008100121A - Index arithmetic method related to blood flow movement state of capillary blood vessel in cerebral tissue, apparatus therefor, and storage medium - Google Patents

Index arithmetic method related to blood flow movement state of capillary blood vessel in cerebral tissue, apparatus therefor, and storage medium Download PDF

Info

Publication number
JP2008100121A
JP2008100121A JP2008010951A JP2008010951A JP2008100121A JP 2008100121 A JP2008100121 A JP 2008100121A JP 2008010951 A JP2008010951 A JP 2008010951A JP 2008010951 A JP2008010951 A JP 2008010951A JP 2008100121 A JP2008100121 A JP 2008100121A
Authority
JP
Japan
Prior art keywords
pixel
pixels
time
curve
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
JP2008010951A
Other languages
Japanese (ja)
Other versions
JP4714228B2 (en
Inventor
Kyojiro Nanbu
恭二郎 南部
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.)
Toshiba Corp
Original Assignee
Toshiba Corp
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 Toshiba Corp filed Critical Toshiba Corp
Priority to JP2008010951A priority Critical patent/JP4714228B2/en
Publication of JP2008100121A publication Critical patent/JP2008100121A/en
Application granted granted Critical
Publication of JP4714228B2 publication Critical patent/JP4714228B2/en
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

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/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/507Apparatus 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 determination of haemodynamic parameters, e.g. perfusion CT

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Engineering & Computer Science (AREA)
  • Optics & Photonics (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Dentistry (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Physics & Mathematics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

<P>PROBLEM TO BE SOLVED: To enhance the analytical precision of CBP study by supporting a decrease in space and time resolving powers to suppress noise. <P>SOLUTION: This index arithmetic apparatus is equipped with a coherent filter processing part 110 wherein the similarity between pixels constituting the plurality of respective CT images continuously acquired by setting the brain of a subject into which a contrast medium having no cerebral blood vessel permeability is injected as a photographing target is judged on the basis of a time density curve of respective pixels over plurality of the CT images and the weighted average of the pixels is locally taken at every image of the basis of the weight corresponding to the similarity and a CBP study processing part 120 wherein the time density curve of cerebral artery pixels and the time density curve of cerebral tissue pixels are formed from a plurality of CT images subjected to weighted average and the transmission function of the time density curve of cerebral tissue pixels to the time density curve of cerebral artery pixels is calculated and CBP, CBV, MTT, and Err related to the blood flow kinetics of blood capillaries in cerebral tissue are operated on the basis of the transmission function. <P>COPYRIGHT: (C)2008,JPO&INPIT

Description

本発明は、脳組織内毛細血管の血流動態に関するインデックスの演算方法及び演算装置に関する。   The present invention relates to an index calculation method and calculation apparatus related to blood flow dynamics of capillaries in brain tissue.

従来からX線CT検査では単純CT像から形態情報を、造影CTによるダイナミックスキャンで病巣の周りの血流の動態情報をそれぞれ視覚情報として得ることができる。近年、マルチスライスCTによる高速スキャンが可能になり益々造影CTのダイナミックスキャンの活用範囲が拡大していくものと考えられる。   Conventionally, in X-ray CT examinations, morphological information can be obtained from simple CT images, and dynamic information on blood flow around a lesion can be obtained as visual information by dynamic scanning using contrast CT. In recent years, high-speed scanning by multi-slice CT has become possible, and it is considered that the range of utilization of dynamic scanning by contrast CT will be expanded.

その一つの方向性として、脳組織内毛細血管の血流動態に関するインデックスを演算するためのCBPスタディと呼ばれる方法がある。CBPスタディでは、”毛細血管を通過する血流”の動態を定量的に表すCBP、CBV、MTT、Err等のインデックスを求め、またこれらインデックスのマップを出力する。   As one of the directions, there is a method called a CBP study for calculating an index related to blood flow dynamics of capillaries in brain tissue. In the CBP study, indexes such as CBP, CBV, MTT, and Err that quantitatively indicate the dynamics of “blood flow through capillaries” are obtained, and a map of these indexes is output.

CBPは、脳組織内の単位体積及び単位時間あたりの血流量[ml/100ml/min]を表し、CBVは、脳組織内の単位体積あたりの血液量[ml/100ml]、MTTは毛細血管の血液平均通過時間[秒]を表している。   CBP represents the blood volume per unit volume and unit time in the brain tissue [ml / 100 ml / min], CBV represents the blood volume per unit volume in the brain tissue [ml / 100 ml], and MTT represents the capillary blood vessel. It represents the average blood passage time [seconds].

これら脳組織中の毛細血管の血流動態を定量的に表しているインデックスCBP、CBV、MTTは、脳虚血が発症してからの経過時間情報ととともに、虚血性脳血管障害の病体鑑別、毛細血管の拡大の有無、血流速などの評価のための有益な情報として期待されている。例えば、一般に虚血性の脳血管障害では、提供する動脈の血圧が低下し、その血管内の血流速の低下が見られる。その結果、CBVは一定でも、MTTが延長し、CBPは低下する。また、脳梗塞超急性期では、血圧低下による血流速の低下を補うために、毛細血管を拡張させ、血流速を増加させることにより、血流量CBPの低下を抑制しようとする働き(オートレギュレーション)がある。従って、MTTが延長することにより、CBPが低下しても、CBVが増加していれば、毛細血管の再開通の可能性を示唆する情報となる。   Indexes CBP, CBV, and MTT, which quantitatively represent the blood flow dynamics of the capillaries in these brain tissues, together with the elapsed time information from the onset of cerebral ischemia, are used to identify the pathophysiology of ischemic cerebrovascular disorders, It is expected to be useful information for evaluating the presence or absence of capillaries and blood flow rate. For example, generally in an ischemic cerebrovascular disorder, the blood pressure of the provided artery is reduced, and the blood flow velocity in the blood vessel is reduced. As a result, even if CBV is constant, MTT is prolonged and CBP is lowered. In the cerebral infarction hyperacute phase, in order to compensate for the decrease in blood flow rate due to the decrease in blood pressure, the capillaries are expanded and the blood flow rate is increased to suppress the decrease in blood flow CBP (auto Regulation). Therefore, even if CBP decreases due to the extension of MTT, if CBV increases, it is information suggesting the possibility of reopening of capillaries.

CBPスタディではトレーサーとして脳血管透過性を持たない造影剤、例えばヨード造影剤が使用される。ヨード造影剤は例えばインジェクターにより肘静脈から注入される。インジェクターにより静注されたヨード造影剤は、心臓、肺を経由して、脳動脈から流れ込む。そして、脳動脈から、脳組織内の毛細血管を経て、脳静脈へと流れ出ていく。このとき、ヨード造影剤は正常な脳組織内の毛細血管では造影剤は血管外へ漏れ出ることなく通過する。図1はこの様子を模式的に示している。   In the CBP study, a contrast agent having no cerebral vascular permeability, such as an iodine contrast agent, is used as a tracer. The iodine contrast agent is injected from the cubital vein by an injector, for example. The iodine contrast agent intravenously injected by the injector flows from the cerebral artery via the heart and lungs. Then, it flows out from the cerebral artery to the cerebral vein through the capillaries in the brain tissue. At this time, the iodine contrast medium passes through the capillary blood vessels in normal brain tissue without leaking out of the blood vessels. FIG. 1 schematically shows this state.

造影剤の通過の様子をダイナミックCTで撮影して、その連続画像から、脳動脈の画素の時間濃度曲線Ca(t)、毛細血管を含む脳組織の画素の時間濃度曲線Ci(t)、脳静脈の画素の時間濃度曲線Csss(t)をそれぞれ測定する。   The state of passage of the contrast agent is photographed by dynamic CT, and from the continuous image, the temporal density curve Ca (t) of the cerebral artery pixel, the temporal density curve Ci (t) of the brain tissue pixel including the capillary, the brain The time density curve Csss (t) of each vein pixel is measured.

ここで、CBPスタディでは、造影剤の血中濃度について脳組織に近い脳血管の血中濃度Ca(t)と毛細血管の血中濃度Ci(t)との間で成り立つ理想的な関係を解析モデルとしており、つまり脳組織に入る直前の血管から造影剤を注入した場合、毛細血管を含む脳組織単位体積(1画素)内の時間濃度曲線は立ち上がりが垂直で、しばらくは一定の値を維持し、その後、急勾配で立ち下がる形になる。これを矩形関数で近似する(box−MTF法:box−Modulation Transfer Function method)。   Here, in the CBP study, an ideal relationship established between the blood concentration Ca (t) of the cerebral blood vessels close to the brain tissue and the blood concentration Ci (t) of the capillaries is analyzed for the blood concentration of the contrast agent. This is a model. In other words, when a contrast medium is injected from a blood vessel just before entering the brain tissue, the time density curve in the brain tissue unit volume (1 pixel) including the capillaries has a vertical rise and remains constant for a while. After that, it will fall down steeply. This is approximated by a rectangular function (box-MTF method: box-Modulation Transfer Function method).

脳動脈血中時間濃度曲線Ca(t)を入力関数、脳組織の時間濃度曲線Ci(t)を出力関数として、トレーサーが毛細血管を通過する過程の特徴を、矩形の伝達関数として求めることができる。   Using the cerebral artery blood time concentration curve Ca (t) as an input function and the brain tissue time concentration curve Ci (t) as an output function, the characteristics of the process of the tracer passing through the capillary can be obtained as a rectangular transfer function. .

このようなCBPスタディにおける第1の問題は次の通りである。肘静脈に、ボーラスインジェクションを行った場合、CTで観察される造影効果は、血液のCT値(造影されない時数十HU)が最大数百HUに上昇する。しかし、脳血流を有効に解析するためには造影効果を高々数パーセント以内の誤差で計測できなくてはならない。すなわち、血液の造影効果(CT値の上昇)が20〜40HU程度であってもこれを検出できる必要がある。   The first problem in such a CBP study is as follows. When bolus injection is performed on the cubital vein, the contrast effect observed by CT increases the CT value of blood (several tens of HU when not contrasted) to a maximum of several hundred HU. However, in order to effectively analyze cerebral blood flow, the contrast effect must be measurable with an error of no more than a few percent. That is, even if the blood contrast effect (CT value increase) is about 20 to 40 HU, it needs to be detected.

単位体積の脳組織中に占める毛細血管の体積比率は高々3〜4パーセント程度である。従って、血液のCT値が20〜40HU上昇した場合、脳組織の平均CT値は、0.5〜1.5HU程度上昇するに過ぎない。   The volume ratio of capillaries in a unit volume of brain tissue is at most about 3 to 4 percent. Therefore, when the CT value of blood increases by 20 to 40 HU, the average CT value of brain tissue only increases by about 0.5 to 1.5 HU.

CT画像ではノイズの標準偏差(sd)は、照射X線量の平方根に反比例し、典型的な照射条件においてsdは例えば5〜10HU程度である。従って、0.5HUの造影効果を検出するためには、X線量を10〜100倍程度増やさねばならず、これは患者の被曝線量が著しく大きくなることを意味する。さらに、ダイナミックCTにおいては同一箇所を数十回に渡って撮影するのであるから、撮影箇所に於ける皮膚の被曝は通常の数百〜数千倍に至ることになり、炎症・脱毛・壊死・発癌等の放射線障害が生じる恐れがある。   In CT images, the standard deviation (sd) of noise is inversely proportional to the square root of the irradiation X-ray dose, and sd is, for example, about 5 to 10 HU under typical irradiation conditions. Therefore, in order to detect a contrast effect of 0.5 HU, the X-ray dose must be increased by about 10 to 100 times, which means that the patient exposure dose is significantly increased. Furthermore, in the dynamic CT, the same part is photographed several tens of times, and therefore, the skin exposure at the photographing part is several hundred to several thousand times the normal, and inflammation, hair loss, necrosis, Radiation damage such as carcinogenesis may occur.

むしろダイナミックCTにおいてはX線量を通常の撮影よりも減らさなくてはならない。一般に、1スキャン当たりのX線量を例えば通常の1/2〜1/10程度に減じることが行われる。これによって、通常の1回のCT撮影に比べて数倍〜20倍程度のX線被曝に留めることができ、これは放射線障害を生じない程度である。しかし、このようなX線量を低減したCT画像において、sdは例えば15〜20HU程度であり、0.5〜1.5HU程度の造影効果は到底検出できない。   Rather, in dynamic CT, the X-ray dose must be reduced compared to normal imaging. In general, the X-ray dose per scan is reduced to, for example, about 1/2 to 1/10 of the usual amount. As a result, the X-ray exposure can be reduced to several to 20 times that of normal one CT imaging, which is a level that does not cause radiation damage. However, in such a CT image with a reduced X-ray dose, sd is, for example, about 15 to 20 HU, and a contrast effect of about 0.5 to 1.5 HU cannot be detected.

そこで、画像のノイズ成分を抑制することが、CBPスタディでは最も重要な課題の1つである。そのため、1)スライス厚を厚くする、2)画素束ね、3)画像の平滑化、が一般的に取りうる方策である。しかしこれらには以下のような問題点がある。   Therefore, suppressing the noise component of the image is one of the most important issues in the CBP study. For this reason, 1) increasing the slice thickness, 2) bundling pixels, and 3) smoothing the image are generally available measures. However, these have the following problems.

“スライス厚を厚くする”ために、撮影時にスライス厚を厚く設定するか、連続する薄いスライスの画像数枚を平均して厚いスライスの画像を生成する。スライス厚に比例して画素当たりのX線量が増えるため、画像ノイズのsdは、スライス厚の平方根に反比例して小さくなる。しかしながら、スライス厚を厚くすることによって、パーシャルボリューム効果が生じ、すなわち1個の画素が、一様な脳組織を表して居らず、複数の組織(白質、灰白質、血管、脳溝、脳室など)の平均的なCT値を表すことになる確率が大きくなり、解析結果として得られる脳血流量等の値の誤差が大きくなる。   In order to “increase the slice thickness”, the slice thickness is set to be thick at the time of shooting, or an image of a thick slice is generated by averaging several consecutive thin slice images. Since the X-ray dose per pixel increases in proportion to the slice thickness, the image noise sd decreases in inverse proportion to the square root of the slice thickness. However, by increasing the slice thickness, a partial volume effect occurs, that is, one pixel does not represent a uniform brain tissue, and multiple tissues (white matter, gray matter, blood vessels, sulci, ventricles) Etc.), and the error in the value of the cerebral blood flow obtained as an analysis result increases.

特に血管の影響を含む画素は、正常な解析が不可能である。このためスライス厚を厚くすると、不正確で、しかも解析不可能な画素を沢山含む非常に品質の悪い結果しか得られなくなる。   In particular, normal analysis is not possible for pixels including the influence of blood vessels. For this reason, when the slice thickness is increased, only a very poor result including many pixels that are inaccurate and cannot be analyzed can be obtained.

“画素束ね”では、空間解像度が或る程度犠牲になる。例えば一辺がn個の画素からなる正方形の領域(n×n個の画素を含む)の平均値を求め、これをその正方形全体の平均CT値とし、このような正方形を画素とみなし、これを敷き詰めて「画素束ね画像」を構成する。もとの画像が例えば一辺512個の画素からなる(512×512個の画素を含む)とし、n=2とすれば、「画素束ね画像」は一辺(512/2)個の画素から構成される(256×256個の画素を含む)画像となる。この方法によれば、ノイズは、nに反比例して減少させることが可能である。さらに、解析対象となる画素の数が1/(n×n)倍になるため、計算量も小さくなるという利点がある。   “Pixel bundling” sacrifices some spatial resolution. For example, an average value of a square area (including n × n pixels) having n pixels on one side is obtained, and this is set as an average CT value of the entire square. Such a square is regarded as a pixel, It spreads out to form a “pixel bundled image”. For example, if the original image is composed of 512 pixels on one side (including 512 × 512 pixels) and n = 2, the “pixel bundle image” is composed of pixels on one side (512/2). (Including 256 × 256 pixels). According to this method, noise can be reduced in inverse proportion to n. Furthermore, since the number of pixels to be analyzed becomes 1 / (n × n) times, there is an advantage that the amount of calculation is also reduced.

しかしながら、nを大きくすると、空間解像度が低下し、それに伴ってパーシャルボリューム効果が生じ、すなわち1個の画素が、一様な脳組織を表して居らず、複数の組織(白質、灰白質、血管、脳溝、脳室など)の平均的なCT値を表すことになる確率が大きくなり、解析結果として得られる脳血流量等の値の誤差が大きくなる。特に血管の影響を含む画素は、正常な解析が不可能である。このため、nを大きくすると、空間解像度が低く、不正確で、しかも解析不可能な画素を沢山含む非常に品質の悪い結果しか得られなくなる。このため、実用上は、n=2〜4程度が限界であり、これだけでは十分なノイズ抑制効果が得られない。   However, when n is increased, the spatial resolution is reduced, and a partial volume effect is caused accordingly. That is, one pixel does not represent a uniform brain tissue, and a plurality of tissues (white matter, gray matter, blood vessels) , Cerebral sulcus, ventricle, etc.) is likely to represent an average CT value, and an error in values such as cerebral blood flow obtained as an analysis result increases. In particular, normal analysis is not possible for pixels including the influence of blood vessels. For this reason, if n is increased, only a very poor result including a low spatial resolution, inaccurate and unanalyzable pixels can be obtained. For this reason, in practical use, n = 2 to 4 is the limit, and a sufficient noise suppression effect cannot be obtained only by this.

また、画像の平滑化、すなわち1枚のCT画像ごとに、2次元の空間フィルタを作用させて平滑化を行う方法を用いると、十分なノイズ抑制効果と引き換えに、空間解像度が著しく損なわれる。特に、太い血管(動脈・静脈)が存在する箇所に近接している画素には、太い血管において生じた造影効果の影響が及ぶことになり、これらの画素の時間濃度曲線は正しくなくなってしまう。従ってごく軽度の平滑化を行うに留めねばならない。ここで、ごく軽度の平滑化を行うに際して重要なのは、画像フィルタのサイズをごく小さくする事、例えば、3×3程度に設定することである。3×3の平滑化フィルタを用いて最大の画像ノイズ抑制効果を得ようとすると、その上限は、ノイズsdを1/3に低減することであり、それ以上にノイズを抑制するのは不可能である。従って十分なノイズ抑制効果は得られない。   Further, if a method of smoothing an image, that is, smoothing by applying a two-dimensional spatial filter for each CT image, the spatial resolution is significantly impaired in exchange for a sufficient noise suppression effect. In particular, a pixel close to a location where a thick blood vessel (artery / vein) exists is affected by the contrast effect generated in the thick blood vessel, and the time density curve of these pixels becomes incorrect. Therefore, only a slight smoothing must be performed. Here, it is important to make the size of the image filter very small, for example, to set it to about 3 × 3 when performing very smooth smoothing. When trying to obtain the maximum image noise suppression effect using a 3 × 3 smoothing filter, the upper limit is to reduce the noise sd to 1/3, and it is impossible to suppress the noise further. It is. Therefore, a sufficient noise suppression effect cannot be obtained.

一方、時間的平滑化、すなわち各画素について得られた時間濃度曲線を曲線とみなして、これを1次元フィルターで平滑化する手法を用いると、十分なノイズ抑制効果を得ようとすると時間分解能を著しく損なう。元来、CBPスタディでダイナミックCTを行うのは短いサンプリング周期で撮影を行うことによって高い時間分解能を得て、時間濃度曲線の僅かな変化(特に生理学的構造に起因する平滑化効果がどの程度生じているか)を精密に計測することが目的であるから、時間的平滑化は全く適当でない。   On the other hand, temporal smoothing, that is, using a technique in which the time density curve obtained for each pixel is regarded as a curve and smoothed with a one-dimensional filter, the time resolution is reduced to obtain a sufficient noise suppression effect. Significantly damaged. Originally, dynamic CT in CBP studies is performed with a short sampling period to obtain a high temporal resolution, and a slight change in the time density curve (especially how much smoothing effect is caused by physiological structure) Therefore, temporal smoothing is not appropriate at all.

本発明の目的は、空間及び時間分解能の低下を抑えて、ノイズを抑制することにより、CBPスタディの解析精度を向上することにある。   An object of the present invention is to improve the analysis accuracy of a CBP study by suppressing noise by suppressing a decrease in spatial and temporal resolution.

本発明は、脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のCT画像各々を構成する画素間の類似度を、前記複数のCT画像にわたる各画素の時間濃度曲線に基づいて判定し、前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均し、前記加重平均された複数のCT画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成し、前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数を計算し、前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算する。   In the present invention, the similarity between pixels constituting each of a plurality of CT images obtained by continuously capturing the brain of a subject into which a contrast medium having no cerebral vascular permeability has been injected as a subject of imaging is represented by the plurality of CT images. Based on the time density curve of each pixel, weighted average of pixels for each image with weights according to the similarity, and time density curves of cerebral artery pixels from the plurality of weighted averaged CT images And a time density curve of the brain tissue pixel, calculate a transfer function of the time density curve of the brain tissue pixel with respect to the time density curve of the cerebral artery pixel, and based on the transfer function, capillary blood vessels in the brain tissue Multiple types of indices related to blood flow dynamics are calculated.

本発明によれば、コヒーレントフィルタ又はコヒーレントレグレッションを併用することにより、空間及び時間分解能の低下を抑えて、ノイズを抑制し、それによりCBPスタディの解析精度を向上することができる。   According to the present invention, by using a coherent filter or coherent regression together, it is possible to suppress a decrease in space and time resolution and suppress noise, thereby improving the analysis accuracy of a CBP study.

以下、図面を参照して本発明を好ましい実施形態により説明する。
本実施形態の特徴としては、コヒーレントフィルタを用いて空間及び時間分解能の低下を抑えながら、ノイズを効果的に抑制することにより、脳組織の毛細血管の血流動態を定量的に表すCBP等のインデックスを高精度に計算することにある。
Hereinafter, the present invention will be described by way of preferred embodiments with reference to the drawings.
As a feature of this embodiment, a coherent filter is used to effectively suppress noise while suppressing a decrease in spatial and temporal resolution, and thus, such as CBP that quantitatively represents blood flow dynamics of capillaries in brain tissue. It is to calculate the index with high accuracy.

(装置構成)
図2には、本実施形態に係るX線CT装置の構成を示している。X線CT装置は、ガントリ部10とコンピュータ装置20とから構成される。ガントリ部10は、X線管101、高電圧発生装置101a、X線検出器102、データ収集部103(DAS;Data Aquisition System)とを有する。X線管101とX線検出器102とは、高速で且つ連続的に回転する図示しない回転リングに被検体Pを挟んで互いに対向する位置に搭載される。
(Device configuration)
FIG. 2 shows the configuration of the X-ray CT apparatus according to the present embodiment. The X-ray CT apparatus includes a gantry unit 10 and a computer device 20. The gantry unit 10 includes an X-ray tube 101, a high voltage generator 101a, an X-ray detector 102, and a data acquisition unit 103 (DAS; Data Acquisition System). The X-ray tube 101 and the X-ray detector 102 are mounted at positions facing each other with the subject P sandwiched between rotating rings (not shown) that rotate continuously at high speed.

コンピュータ装置20は、画像処理装置30と、画像表示部107と、入力部109とから構成される。画像処理装置30は、制御部108を中枢として、データ収集部103から出力される生データを補正処理等を経て投影データに変換する前処理部104、投影データを記憶するメモリ部105、投影データからCT画像データを再構成する画像再構成部106、CT画像データを保管する記憶装置10M、CT画像データに対してコヒーレントフィルタ処理を実行するコヒーレントフィルタ処理部110、及びコヒーレントフィルタ処理を受けたCT画像データを使ってCBPスタディ処理を実行するCBPスタディ処理部120とから構成される。   The computer device 20 includes an image processing device 30, an image display unit 107, and an input unit 109. The image processing apparatus 30 has a control unit 108 as a center, a pre-processing unit 104 that converts raw data output from the data collection unit 103 into projection data through correction processing and the like, a memory unit 105 that stores projection data, and projection data Image reconstructing unit 106 for reconstructing CT image data from, storage device 10M for storing CT image data, coherent filter processing unit 110 for performing coherent filter processing on CT image data, and CT that has undergone coherent filter processing And a CBP study processing unit 120 that executes CBP study processing using image data.

コヒーレントフィルタ処理部110は、分散値推定部111、重み関数演算部112、画素値演算部(コヒーレントフィルタ部)113とから構成される。これら分散値推定部111、重み関数演算部112、画素値演算部113の機能については後述するコヒーレントフィルタ処理の詳細説明の中で説明する。   The coherent filter processing unit 110 includes a variance value estimation unit 111, a weight function calculation unit 112, and a pixel value calculation unit (coherent filter unit) 113. The functions of the variance value estimation unit 111, the weight function calculation unit 112, and the pixel value calculation unit 113 will be described in the detailed description of coherent filter processing described later.

CBPスタディ処理部120は、ROI設定支援部121、時間濃度曲線作成部122、脳動脈時間濃度曲線補正部123、MTF処理部124、インデックス計算部125、マップ作成部126、マップ合成部127から構成される。   The CBP study processing unit 120 includes an ROI setting support unit 121, a time concentration curve creation unit 122, a cerebral artery time concentration curve correction unit 123, an MTF processing unit 124, an index calculation unit 125, a map creation unit 126, and a map composition unit 127. Is done.

ROI設定支援部121は、CT画像上に脳動脈や脳静脈に対して関心領域ROIを設定する作業を支援するための情報(脳動脈ROIのためのATマップ、PTマップ、TTマップ等)を作成し提供する。   The ROI setting support unit 121 provides information (AT map, PT map, TT map, etc. for the cerebral artery ROI) for supporting the operation of setting the region of interest ROI for the cerebral artery and cerebral vein on the CT image. Create and provide.

なお、脳動脈ROIは、前大脳動脈(ACA)、中大脳動脈(MCA)、後大脳動脈(PCA)を対象として、左脳、右脳それぞれの領域に個別に設定される。従って、左右に3個ずつ、合計6個の脳動脈ROIが設定される。また、脳動脈の時間濃度曲線Ca(t)を補正するために、他の時間濃度曲線Csss(t)が利用される。この時間濃度曲線Csss(t)は、パーシャルボリュームを含まない画素が存在するのに充分に太い血管上に設定されたROIに関して作成される。例えば、脳血管の中で最も太い上矢状静脈洞に設定することが好ましい。   The cerebral artery ROI is individually set in the left brain region and the right brain region for the anterior cerebral artery (ACA), middle cerebral artery (MCA), and posterior cerebral artery (PCA). Accordingly, a total of six cerebral arteries ROI are set, three on each side. Further, another time concentration curve Csss (t) is used to correct the time concentration curve Ca (t) of the cerebral artery. This time density curve Csss (t) is created for an ROI set on a blood vessel that is sufficiently thick so that there are pixels that do not contain a partial volume. For example, it is preferable to set the thickest upper sagittal sinus in the cerebral blood vessel.

時間濃度曲線作成部122は、記憶装置10Mに記憶されているダイナミックCT画像データ(時間的に連続した複数枚の画像データ)から脳動脈、脳静脈及び脳組織(毛細血管)に関する時間濃度曲線を作成する。なお、脳動脈の時間濃度曲線Ca(t)は、設定された6つの脳動脈ROIに関して個々に作成される。脳静脈の時間濃度曲線Csss(t)は、上矢状静脈洞に設定された脳静脈ROIに関して作成される。また、脳組織の時間濃度曲線Ci(t)は、脳組織上の全画素を対象として画素ごとに作成される。   The time concentration curve creation unit 122 generates time concentration curves related to cerebral arteries, cerebral veins, and brain tissues (capillaries) from dynamic CT image data (a plurality of temporally continuous image data) stored in the storage device 10M. create. The cerebral artery time concentration curve Ca (t) is individually created for the set six cerebral arteries ROI. The cerebral venous time concentration curve Csss (t) is created for the cerebral venous ROI set in the upper sagittal sinus. Further, the time density curve Ci (t) of the brain tissue is created for each pixel for all the pixels on the brain tissue.

脳動脈時間濃度曲線補正部123は、ノイズやパーシャルボリューム効果の影響を除去するために、脳動脈の時間濃度曲線Ca(t)を、上矢状静脈洞の時間濃度曲線Csss(t)に基づいて補正する。この補正方法については後述する。MTF処理部124は、補正された脳動脈の時間濃度曲線Ca(t)と、脳組織の時間濃度曲線Ci(t)とに基づいて、box−MTF法により、伝達関数MTFを、脳組織上の全画素を対象として画素ごとに計算する。   The cerebral artery time concentration curve correction unit 123 uses the cerebral artery time concentration curve Ca (t) based on the upper sagittal sinus time concentration curve Csss (t) in order to remove the influence of noise and the partial volume effect. To correct. This correction method will be described later. The MTF processing unit 124 calculates the transfer function MTF on the brain tissue by the box-MTF method based on the corrected cerebral artery time concentration curve Ca (t) and the brain tissue time concentration curve Ci (t). The calculation is performed for each pixel for all the pixels.

インデックス計算部125は、計算された伝達関数MTFから脳組織の血流動態を表すインデックス(CBP、CBV、MTT、Err)を、脳組織上の全画素を対象として画素ごと計算する。マップ作成部126は、計算されたインデックス各々のマップを、脳動脈(ACA,MCA,PCA)ごとに生成する。このマップは、各スライスに関して、インデックスの種類(=4)×脳動脈の種類の数(ACA,MCA,PCAの3つ)=12種類作成される。マルチスライスでは、そのスライス数倍の種類のマップが作成される。この膨大な枚数のマップの枚数を合成処理により減らして診断効率を向上させるためにマップ合成部127が設けられている。   The index calculation unit 125 calculates an index (CBP, CBV, MTT, Err) representing the blood flow dynamics of the brain tissue from the calculated transfer function MTF for every pixel on the brain tissue. The map creation unit 126 generates a map of each calculated index for each cerebral artery (ACA, MCA, PCA). This map is created for each slice, index type (= 4) × number of types of cerebral arteries (3 of ACA, MCA, and PCA) = 12 types. In multi-slice, the number of types of maps is created that is the number of slices. A map composition unit 127 is provided in order to reduce the enormous number of maps by the composition process and improve the diagnosis efficiency.

以下に、コヒーレントフィルタ処理とCBPスタディ処理について順番に説明する。   Hereinafter, the coherent filter process and the CBP study process will be described in order.

コヒーレントフィルタの原理について簡単に説明すると、近傍の3×3等の局所内画素を加重平均し、その加重平均値を局所中心画素の値とすることを基本として、周辺画素各々の重みを中心画素と周辺画素との間の類似度に従って変えることを特徴としたものである。ここで言う類似度とは、画素間で、解剖学的に近い組織、具体的には同じ脳動脈の支配下にある毛細血管どうしである可能性の度合いを示す指標であり、この類似度が高い画素に対しては高い重みを与え、逆に類似度が低い画素に対してはゼロに近い低い重みを与えることにより、ノイズ抑制を果たしながらも、空間分解能の低下を抑制することを可能としている。ここで類似度の計算が重要になるが、本実施形態では、脳血管透過性を持たない造影剤、例えばヨード造影剤を注入(静注)された被検体の脳を撮影対象として連続的に取得した複数のCT画像(ダイナミックCT画像)を用いて、各画素の時間濃度曲線の比較により類似度を計算する。そのため類似度の確からしさは、サンプリング周波数、つまり単位時間あたりの画像枚数と、サンプリング数、つまり全画像枚数とに依存して決まる。そこでスキャン間隔を例えば0.5秒に短縮することが効果的である。   The principle of the coherent filter will be briefly described. Based on the weighted average of neighboring local pixels such as 3 × 3 and the weighted average value as the value of the local central pixel, the weight of each peripheral pixel is determined as the central pixel. It is characterized in that it is changed according to the similarity between the pixel and the peripheral pixel. The similarity referred to here is an index indicating the degree of possibility that the pixels are close to the anatomical structure, specifically, the capillaries under the control of the same cerebral artery between the pixels. By giving a high weight to high pixels and conversely giving a low weight close to zero to pixels with low similarity, it is possible to suppress a decrease in spatial resolution while suppressing noise. Yes. Here, the calculation of the similarity is important. In this embodiment, a contrast medium having no cerebral vascular permeability, for example, the brain of a subject into which an iodinated contrast medium is injected (intravenous injection) is continuously taken as an imaging target. Using the acquired plurality of CT images (dynamic CT images), the similarity is calculated by comparing the time density curves of each pixel. Therefore, the probability of similarity is determined depending on the sampling frequency, that is, the number of images per unit time, and the number of samplings, that is, the total number of images. Therefore, it is effective to shorten the scan interval to 0.5 seconds, for example.

(コヒーレントフィルタ)
(コヒーレントフィルタの一般的説明)
(画素値v(x))
一般に、カメラ等の撮像手段を介して取得されたデジタル画像は、複数の画素(pixel)から構成されている(あるいは、当該画像をそのような画素の集合として考えることができる。)。以下の説明では、当該画素の位置をベクトルx(すなわち座標値のベクトル)として表し、画素xが有する値(例えば濃淡を表わす数値、CT値HU)をK次元ベクトルとして表す。2次元画像の場合、画素xとは画像上における位置を表す座標値(x、y)を示す2次元ベクトルである。ある画素xについて定義される「画素値v(x)」を、
v(x)=(v(x),v(x),…,v(x)) … (1)
と表記する。なお、この(1)式の右辺における、v(x),v(x),…,v(x)それぞれを、以下では、画素xについての「スカラー値」と呼ぶことにする。
(Coherent filter)
(General description of coherent filter)
(Pixel value v (x))
In general, a digital image acquired via an imaging means such as a camera is composed of a plurality of pixels (or the image can be considered as a set of such pixels). In the following description, the position of the pixel is expressed as a vector x (that is, a vector of coordinate values), and a value (for example, a numerical value indicating shading, CT value HU) of the pixel x is expressed as a K-dimensional vector. In the case of a two-dimensional image, the pixel x is a two-dimensional vector indicating coordinate values (x, y) representing a position on the image. The “pixel value v (x)” defined for a pixel x is
v (x) = (v 1 (x), v 2 (x),..., v K (x)) (1)
Is written. Note that each of v 1 (x), v 2 (x),..., V K (x) on the right side of the equation (1) is hereinafter referred to as a “scalar value” for the pixel x.

例えば、画像が「カラー画像」であるとき、各画素が、それぞれ三原色(赤,緑,青)の明るさ(スカラー値)を有することから、これら各画素の画素値v(x)は、その次元がK=3のベクトルであると考えることができる(上記(1)式の右辺各項で、その添え字が例えば「赤」,「緑」及び「青」である場合を想定されたい。また例えば、画像がK枚の静止画像から構成される動画像であって、第n番目の画像の各画素はスカラー値v(x)を持つという場合には、K枚の静止画像上、共通の同一点(同一座標)の画素xの持つ画素値(スカラー値)を並べて構成される、K次元ベクトル値v(x)=(v(x),v(x),…,v(x))が以下で述べるベクトル値としての画素値である。 For example, when the image is a “color image”, each pixel has brightness (scalar value) of three primary colors (red, green, and blue), and therefore the pixel value v (x) of each pixel is It can be considered that the dimension is a vector of K = 3 (assuming the case where the subscripts are “red”, “green”, and “blue” in each term on the right side of the above equation (1). For example, if the image is a moving image composed of K still images and each pixel of the nth image has a scalar value v n (x), K-dimensional vector values v n (x) = (v 1 (x), v 2 (x),..., Arranged by arranging pixel values (scalar values) of pixels x at the same common point (same coordinates). v K (x)) is a pixel value as a vector value described below.

(類似度(適合度ないし危険率)p(x,y)と重みw(p(x,y)))
上記画素xに対して、適当な画素の集合N(x)を考える(この集合N(x)は画素xを含んでよい。)。次に、N(x)の要素であるそれぞれの画素yと、前記画素xとの間で、重みw(p(x,y))を考える。この重みw(p(x,y))は、次に記す性質を有する。
(Similarity (fitness or risk factor) p (x, y) and weight w (p (x, y)))
Consider an appropriate pixel set N (x) for the pixel x (this set N (x) may include the pixel x). Next, a weight w (p (x, y)) is considered between each pixel y which is an element of N (x) and the pixel x. The weight w (p (x, y)) has the following properties.

(類似度p(x,y))
まず、重みw(p(x,y))の値を左右する関数p(x,y)の意味について述べる。このp(x,y)は、本実施形態にいう「類似度」を定量化する手段であり、一般的にいえば、画素xと画素y∈N(x)とが、何らかの意味でどの程度類似しているか(例えば、両画素x及びyの上記画素値v(x)及びv(y)間に認められる統計的差異の程度)、を示す具体的数値を与える。
(Similarity p (x, y))
First, the meaning of the function p (x, y) that affects the value of the weight w (p (x, y)) will be described. This p (x, y) is a means for quantifying the “similarity” referred to in the present embodiment, and generally speaking, to what extent the pixel x and the pixel y∈N (x) are in some sense. A specific numerical value indicating whether they are similar (for example, the degree of statistical difference recognized between the pixel values v (x) and v (y) of both pixels x and y) is given.

より具体的には例えば、p(x,y)が小さな値を与えるときには、画素xと画素yとが、その画素値v(x)及びv(y)間に「統計的に有意な差がなく(つまり類似度が高い)」、類似である可能性が高いと判断され、p(x,y)が大きな値を与えるときには、「統計的に有意な差があり(つまり類似度が低い)」、の如く判断されるということである。   More specifically, for example, when p (x, y) gives a small value, the pixel x and the pixel y have a “statistically significant difference between the pixel values v (x) and v (y)”. Is not (that is, the degree of similarity is high) ", and when p (x, y) gives a large value when it is determined that there is a high possibility of being similar," there is a statistically significant difference (that is, the degree of similarity is low). "Is judged as follows.

ところで、画素値v(x)及びv(y)(ないしスカラー値v(x),…,v(x)及びv(y),…,v(y))には、必ずノイズが含まれていると考えなければならない。例えば、画像がCCD撮像素子により取得された場合を考えると、それを構成する各画素については、素子内の暗電流や外界から入射する光量の不規則変動に起因するノイズ等が存在する。 By the way, the pixel values v (x) and v (y) (or scalar values v 1 (x),..., V K (x) and v 1 (y),..., V K (y)) are always noise. Must be considered to be included. For example, when considering a case where an image is acquired by a CCD image pickup device, each pixel constituting the image includes noise or the like due to dark current in the device or irregular fluctuation of the amount of light incident from the outside.

このようなノイズは、一般に、全画素についてまちまちな値をとるため、画素xと画素yとが、仮に(外界における)同一物体を反映したものである場合であっても、実際に観測される画像上では、同一の値を持たないことがある。このことを逆にいえば、いずれも同一物体を反映した画素xと画素yにおいて、それぞれのノイズを除去した状況を仮に想定すれば、これらは該同一物体を表象するものとして画像上に表示され(=そのように認識され)るし、また、両者は本来同一の(あるいはごく近い)画素値を有する。   Since such noise generally takes various values for all pixels, even if the pixel x and the pixel y reflect the same object (in the outside world), they are actually observed. On the image, it may not have the same value. In other words, if it is assumed that the respective noises are removed from the pixel x and the pixel y reflecting the same object, they are displayed on the image as representing the same object. (= Recognized as such) and both have essentially the same (or very close) pixel values.

そこで、上述したノイズの性質を踏まえ、上記のp(x,y)に関し、統計的検定法でよく知られている「帰無仮説」の概念を用いると、このp(x,y)については、具体的に次のように言うことができる。すなわち、帰無仮説H「画素xと画素yとはそれぞれのノイズを除去した場合に同一の画素値を有する」言いかえれば「v(x)=v(y)、ただし、両画素のノイズに起因する差異を除く」を立てる(つまり、このような命題が成立する場合、「両画素x及びyとの間の類似度が高い(適合度が大きい)」と考える。)と、関数p(x,y)は、この仮説Hを棄却する場合の危険率(あるいは、有意水準)であるということができる(この場合、p(x,y)は、その値域が[0,1]であるような関数として定義される(p(x,y)∈[0,1])。)。   Therefore, based on the above-mentioned noise characteristics, regarding the above p (x, y), using the concept of the “null hypothesis” that is well known in the statistical test method, this p (x, y) Specifically, we can say as follows. In other words, the null hypothesis H “pixel x and pixel y have the same pixel value when the respective noises are removed”, in other words, “v (x) = v (y)” "Excluding the difference caused by" (that is, when such a proposition is established, it is considered that "the similarity between both pixels x and y is high (the degree of matching is high)"). x, y) can be said to be a risk factor (or significance level) when rejecting this hypothesis H (in this case, p (x, y) has a value range of [0, 1]. (P (x, y) ε [0, 1])).

したがって、危険率p(x,y)が大きい場合、すなわち棄却が誤りである危険性が大きい場合には上記仮説Hを満たす可能性が高いといえ、逆に小さい場合、すなわち棄却が誤りである危険性が小さい場合には仮説Hを満たさない可能性が高いということができる(なお、統計的検定における周知事項ではあるが、仮説Hが「棄却」されないといっても、それが「真」であることを意味するわけではない。この場合、仮説Hが示す命題が、否定し得ないことを意味するに過ぎない。)。   Therefore, when the risk factor p (x, y) is large, that is, when there is a high risk of rejection being incorrect, it can be said that the above hypothesis H is likely to be satisfied, and conversely, when it is small, that is, rejection is incorrect. If the risk is small, it can be said that there is a high possibility that the hypothesis H is not satisfied (note that although it is a well-known matter in the statistical test, the hypothesis H is not “rejected”, it is “true”. (In this case, it simply means that the proposition indicated by hypothesis H cannot be denied.)

(重みw(p(x,y)))
さて、重みw(p(x,y))は、その表され方から明らかな通り、上記したような危険率p(x,y)の関数(より一般には、適合度の関数(適合度をρ(x,y)とすれば、w(ρ(x,y))となるように構成できる)であり、また、この重みw(p(x,y))を求めるため、x及びyの組み合わせそれぞれについて求められた危険率p(x,y)に作用させる重み関数wは、一般的にいうと、上記「棄却」を具現化する作用を有するものである。具体的には、危険率p(x,y)が大きい場合には重み関数wの値、すなわち重みw(p(x,y))が大きな正の値をとり、その逆の場合には小さな正の値(又は“0”)をとる、等というように調整されている(重み関数wの具体的形式については後述する。)。つまり、重みw(p(x,y))は、画素xと画素yとが、上記仮説Hに示される命題を満たすらしい場合には、大きい値をとり、その逆の場合には小さい値をとる。一例として特に、wのとりうる値が”0”かまたは”0”でない一定値の2通りしかないように構成してもよい。
(Weight w (p (x, y)))
The weight w (p (x, y)) is a function of the risk factor p (x, y) as described above (more generally, a function of fitness (the fitness is If ρ (x, y), it can be configured to be w (ρ (x, y)), and in order to obtain this weight w (p (x, y)), x and y Generally speaking, the weighting function w to be applied to the risk factor p (x, y) obtained for each combination has the effect of embodying the “rejection.” Specifically, the risk factor. When p (x, y) is large, the value of the weighting function w, that is, the weight w (p (x, y)) takes a large positive value, and vice versa, a small positive value (or “0”). ")", Etc. (The specific form of the weight function w will be described later.) That is, the weight w (p (X, y)) takes a large value when the pixel x and the pixel y seem to satisfy the proposition shown in the above hypothesis H, and takes a small value in the opposite case. It may be configured so that there are only two possible values of w, “0” or a fixed value other than “0”.

なお、以上までに述べた仮説H、危険率p(x,y)、重みw(p(x,y))間の関係をまとめると、帰無仮説Hが正しい可能性が高いとき、類似度pも高くなり、その画素に与える重みwを高くし、一方、帰無仮説Hが正しい可能性が低いとき、類似度pも低くなり、その画素に与える重みwを低くする。このように加重平均値への寄与度(重み)を類似度に応じて変えることにより、分解能の低下を抑えながら、ノイズを効果的に抑制することが可能となる。また、重み関数w(t)は、より一般に、「t∈[0,1]で定義される非負の単調増加関数」ということができ、また、該w(t)の満たすべき性質は、少なくともそのようであればよい。   In addition, when the relationship between the hypothesis H, the risk factor p (x, y), and the weight w (p (x, y)) described above is summarized, when the null hypothesis H is highly likely to be correct, When p is also increased and the weight w given to the pixel is increased, on the other hand, when the possibility that the null hypothesis H is correct is low, the similarity p is also lowered and the weight w given to the pixel is lowered. Thus, by changing the contribution (weight) to the weighted average value according to the similarity, it is possible to effectively suppress noise while suppressing a decrease in resolution. Further, the weight function w (t) can be more generally referred to as “a non-negative monotonically increasing function defined by t∈ [0, 1]”, and the property to be satisfied by w (t) is at least That is all that is necessary.

(コヒーレントフィルタ処理)
以上までの説明により、「コヒーレントフィルタ」は次のように導かれる。すなわちまず、画像を構成するある画素xに対し、集合N(x)の要素たる画素yのすべてについて上記した重みw(p(x,y))を計算する。次に、これら複数の重みw(p(x,y))を用いて、当該画素xを構成する新たなスカラー値v´(x)を、以下の(2)式で計算する。すなわち、

Figure 2008100121
(Coherent filter processing)
From the above description, the “coherent filter” is derived as follows. That is, first, the weight w (p (x, y)) described above is calculated for all the pixels y that are elements of the set N (x) for a certain pixel x constituting the image. Next, a new scalar value v ′ k (x) constituting the pixel x is calculated by the following equation (2) using the plurality of weights w (p (x, y)). That is,
Figure 2008100121

ただし、k=1,2,…,Kである。そして、この式で求められたv´(x)を用いて、当該画素xの変換後の画素値(新たな画素値)v´(x)を、
v´(x)=(v´(x),v´(x),…,v´(x)) … (3)
として構成する。
However, k = 1, 2,..., K. Then, using v ′ k (x) obtained by this equation, the pixel value (new pixel value) v ′ (x) after the conversion of the pixel x is
v ′ (x) = (v ′ 1 (x), v ′ 2 (x),..., v ′ K (x)) (3)
Configure as.

ここに、上記(1)式で表される、画素値v(y)=(v(y),v(y),…,v(y))(y=xである場合を含む。)を、v´(x)=(v´(x),v´(x),…,v´(x))に変換するフィルタが、「コヒーレントフィルタ」の形式である。これはその表式から明らかな通り、画素値を構成するスカラー値v(y)の重み付け平均値を表している。 Here, the case where the pixel value v (y) = (v 1 (y), v 2 (y),..., V K (y)) (y = x) expressed by the above equation (1) is included. .) Is converted to v ′ (x) = (v ′ 1 (x), v ′ 2 (x),..., V ′ K (x)) is a “coherent filter” format. As is apparent from the expression, this represents a weighted average value of the scalar values v k (y) constituting the pixel values.

このような処理は、以下のような結果をもたらす。すなわち、画素値v´(x)は、画素xとノイズを除いて同一の画素値をとることが確からしい(=上記仮説Hの命題を満たす可能性が高い)画素yを重視した重み付け平均値v´(x)から構成されたベクトルを表すこととなる。また、このような画素yが十分な数存在するならば、画素値v´(x)は、画素xが本来有すべきその真値から外れることなく、上記したような平均化の作用によりノイズのみを抑制した値を有することとなる。 Such processing yields the following results. That is, the pixel value v ′ (x) is likely to be the same pixel value as the pixel x except for noise (= the weighted average value that places importance on the pixel y that is likely to satisfy the proposition of the above hypothesis H). It represents a vector composed of v ′ k (x). Further, if there are a sufficient number of such pixels y, the pixel value v ′ (x) is not deviated from its true value that the pixel x should originally have, and noise is caused by the above-described averaging operation. It has the value which suppressed only.

なお、危険率p(x,y)が小さく、したがって、帰無仮説Hが「棄却」され、重みw(p(x,y))が小さくなるような場合であっても、上記記述からもわかる通り、必ずしもこれを完全に「棄却」するとは限らない。このようなことは、後述する重み関数wの具体的形式に依存するところであるが、危険率p(x,y)が“0”(=0%)に近いような場合でも、w(p(x,y))≠0(ただし、p(x,y)が“1”に近い場合に比べて、より小さな正の値ではある。)としてよい(なお、p(x,y)=1である場合とは、後述するように、v(x)=v(y)のときである。)。   Even if the risk factor p (x, y) is small and therefore the null hypothesis H is “rejected” and the weight w (p (x, y)) is small, the above description also indicates As you can see, this is not necessarily completely “rejected”. Although this depends on the specific form of the weight function w described later, even when the risk factor p (x, y) is close to “0” (= 0%), w (p ( x, y)) ≠ 0 (where p (x, y) is a smaller positive value than when p (x, y) is close to “1”) (p (x, y) = 1) (Some cases are when v (x) = v (y), as will be described later.)

すなわち、完全な棄却ということではなく、小さな寄与は認めてよいということである(なおこのような場合に、w(p(x,y))=0とするのであれば、完全な棄却を行うのと同義である。   In other words, it is not a complete rejection but a small contribution may be accepted (in this case, if w (p (x, y)) = 0, a complete rejection is performed. Is synonymous with

このような処理は、一般的に次のように言える。すなわち、ある画像を構成する(複数の)画素xが存在するとき、この画素xとある任意の画素y(上記ではy∈N(x)とされた。)との適合度を定量化し(上記では、p(x,y)に基づいていた。)、該適合度が大きい場合には、画素値v(y)を利用した重み付き平均化処理において、当該画素yについて大きな寄与を認め、適合度が小さい場合には小さな寄与しか認めないようにすることで、当該画素xのノイズを有効に抑制する画像処理方法である、といえる。いわば、画素xと画素yとが「似たもの同士」のときには、該画素yを前記平均化処理に、より貢献させ、「似ていないもの同士」のときには、該画素yを殆ど又は全く無視する(重みをゼロ又はその近似値)、と言い換えてもよい。   Such processing can be generally described as follows. That is, when there are (a plurality of) pixels x constituting an image, the degree of matching between this pixel x and an arbitrary pixel y (in the above, yεN (x)) is quantified (above Is based on p (x, y).) When the degree of matching is large, a large contribution is recognized for the pixel y in the weighted averaging process using the pixel value v (y). It can be said that the image processing method effectively suppresses noise of the pixel x by allowing only a small contribution when the degree is small. In other words, when the pixel x and the pixel y are “similar”, the pixel y contributes more to the averaging process. When the pixel x and the pixel y are “similar”, the pixel y is almost or completely ignored. In other words, the weight may be zero or an approximate value thereof.

このような処理を画像全体に施すことにより、画像のぼけ、つまり空間的分解能の低下を殆ど生じることなく、極めて高いノイズ抑制効果を発揮することができる。また、ノイズ抑制という用途に限定せず、例えばパターン認識の分野においても、重み関数、あるいはコヒーレントフィルタを好適な具体的形式にすることによって、優れた効果を発揮することができる。   By applying such processing to the entire image, an extremely high noise suppression effect can be exhibited with almost no blurring of the image, that is, a reduction in spatial resolution. Further, the present invention is not limited to the use of noise suppression. For example, also in the field of pattern recognition, an excellent effect can be exhibited by making a weighting function or a coherent filter suitable in a specific form.

ここで上記した「ダイナミックCT」撮影とは、上記X線管101及びX線検出器102が被検体Pの同一部位を反復撮影(反復スキャン、連続回転型CT装置では、連続回転による反復撮影がしばしば行われる。)して次々に投影データを取得するとともに、該投影データに基づいて次々に再構成処理を行って時系列的な一連の画像を得る撮影方式のことをいう(この場合、画像表示部107における画像表示は、例えば図示しないカウンタ等によって、その画像の元となった投影データ収集に係るスキャン開始点又は終点から一定時間後に行われるように制御される。)。   Here, “dynamic CT” imaging means that the X-ray tube 101 and the X-ray detector 102 repeatedly image the same part of the subject P (repetitive scanning, in a continuous rotation CT apparatus, repeated imaging by continuous rotation is performed). This is an imaging method in which projection data is acquired one after another, and reconstruction processing is performed one after another based on the projection data to obtain a series of time-series images (in this case, images The image display on the display unit 107 is controlled so as to be performed after a predetermined time from the scan start point or end point related to the collection of projection data that is the source of the image, for example, by a counter (not shown).

したがって、このように取得・表示される画像は、映画等と同様に時系列的な複数枚の静止画像からなる、いわゆる動画像となる。なお、このような撮影方式は、典型的には、被検体Pに対し造影剤を注入し、その経時変化を観察・解析して、例えば血管における狭窄や閉塞等その他病変部の病態を分析するために用いられる。また、造影剤投与の前後2回だけに限り同一部位のCT撮影を行う方式も、広義のダイナミックCT撮影と考えることができる。   Therefore, the image acquired / displayed in this way is a so-called moving image composed of a plurality of time-series still images as in a movie or the like. Note that such an imaging method typically injects a contrast medium into the subject P, observes and analyzes changes over time, and analyzes other pathological conditions such as stenosis and occlusion in blood vessels, for example. Used for. A method of performing CT imaging of the same part only twice before and after contrast medium administration can also be considered as dynamic CT imaging in a broad sense.

さて、従来においては、上記のような「ダイナミックCT」撮影時、例えばK回の撮影を実施する間に被検体Pに何らかの変化(例えば、造影剤の濃度変化や呼吸動等が一般的に考えられる)があった場合、空間解像度を損なわず画像ノイズを抑制するためには、時間方向の平滑化を行うほかなかった。その結果、時間分解能が損なわれるという弊害は避け得なかった。   Conventionally, at the time of “dynamic CT” imaging as described above, for example, some change in the subject P (for example, a change in the concentration of a contrast medium, respiratory movement, etc.) is generally considered while K imaging is performed. In order to suppress the image noise without impairing the spatial resolution, smoothing in the time direction was necessary. As a result, the adverse effect that the time resolution is impaired cannot be avoided.

ところが、ダイナミックCT撮影により取得される画像は、上述したように、動画像であって時間的変化を仔細に観察する目的で行うものであるから、その分解能が損なわれるというのは、本来、好ましい状況とは言えない。   However, since the image acquired by dynamic CT imaging is a moving image as described above and is performed for the purpose of closely observing temporal changes, it is originally preferable that the resolution is impaired. It's not a situation.

コヒーレントフィルタを利用すれば、分解能を損ねず、K枚の静止画像のすべて(複数枚の画像)につきそのノイズを抑制することが可能な、次のようなダイナミック・コヒーレントフィルタ処理を実施することができる。   If a coherent filter is used, the following dynamic coherent filter processing that can suppress noise for all (a plurality of images) of K still images without degrading resolution can be performed. it can.

まず、上記のようにして得られた動画像たるK枚の静止画像につき定義される画素xについては、既に述べたように、画素値v(x)として、
v(x)=(v(x),v(x),…,v(x)) … (1再掲)
を構成することができる。ここで右辺各項における添え字1,2,…,Kは、K枚の各静止画像の通し番号である。
First, for the pixel x defined for the K still images, which are the moving images obtained as described above, as described above, as the pixel value v (x),
v (x) = (v 1 (x), v 2 (x),..., v K (x)) (1 reprint)
Can be configured. Here, subscripts 1, 2,..., K in each term on the right side are serial numbers of K still images.

次に、この場合における重み関数w1の具体的形式を、例えば次の(4)式により与える。

Figure 2008100121
Next, a specific form of the weight function w1 in this case is given by the following equation (4), for example.
Figure 2008100121

ただし、y∈N(x)であって、かつ、この集合N(x)は、画素xにつき任意に設定してよい(=どのような基準によって設定してもよい。)。しかし実際上は、画素xと該画素xから遠く離れた位置にある画素yとが仮説「v(x)=v(y)。ただし、両画素のノイズに起因する差異を除く」を満たす可能性は一般に低いといえるから、集合N(x)をxに近接している画素の集合という基準で限定することは、演算速度向上等の実用的な意義がある。   However, yεN (x) and the set N (x) may be arbitrarily set for each pixel x (= may be set according to any criterion). However, in practice, the pixel x and the pixel y far away from the pixel x can satisfy the hypothesis “v (x) = v (y). However, excluding differences caused by noise between the two pixels”. Therefore, limiting the set N (x) on the basis of a set of pixels close to x has practical significance such as an increase in calculation speed.

したがってここでは、その一例として、集合N(x)を、当該画素xを中心としたその周囲の矩形状エリアに含まれる画素の集合、とする。より具体的に、集合N(x)としては、例えば、いま注目している静止画像一枚を構成する全画素が128×128画素であるような場合に、前記画素xを中心とした3×3画素分のエリアとしたり、また、512×512画素であるような場合に、当該画素xを中心とした13×13画素分のエリア等としてもよい。   Therefore, here, as an example, the set N (x) is a set of pixels included in the surrounding rectangular area centered on the pixel x. More specifically, as the set N (x), for example, when all pixels constituting one still image of interest are 128 × 128 pixels, 3 × It may be an area for 3 pixels, or in the case of 512 × 512 pixels, it may be an area for 13 × 13 pixels centered on the pixel x.

また、上記(4)式におけるσは、k枚目の静止画像の各画素が、そのどれにも共通な一定の程度で有するものと仮定して推定されたノイズの標準偏差であり、一方Cは、重みw1(p(x,y))が、上記(4)式に代入された場合における作用の程度を決定調節可能なするパラメータである。 In addition, σ k in the above equation (4) is a noise standard deviation estimated on the assumption that each pixel of the k-th still image has a certain degree common to all the pixels. C is a parameter that can determine and adjust the degree of action when the weight w1 (p (x, y)) is substituted into the above equation (4).

以下、これらσ及びCについての説明を順に行う。
まず、(4)式におけるσについて説明する(以下では、分散σ として説明する。)。このσ は、上述したように、k枚目の静止画像上の各画素のスカラー値が有するノイズ成分の分散である。そしてまた、上記(4)式における分散σ は、k枚目の画像の各画素のスカラー値について一定値たる分散σ を持つノイズを含んでいるものと仮定して推定したものである。一般に、このような仮定は、次に記すようなことを背景として、十分な正当性を持つ。
Hereinafter, these σ k and C will be described in order.
First, σ k in the equation (4) will be described (hereinafter, described as variance σ k 2 ). As described above, σ k 2 is the variance of the noise component of the scalar value of each pixel on the kth still image. The variance σ k 2 in the above equation (4) is estimated on the assumption that the scalar value of each pixel of the k-th image includes noise having a constant variance σ k 2. is there. In general, such assumptions are sufficiently valid against the background described below.

被検体Pの大きさ、X線管101及びX線検出器102、再構成部106等の構造が一定で、かつ、照射X線のエネルギを一定にした状態では、CT画像のノイズは、照射X線量、すなわちこれと比例関係にあるX線管101における管電流と照射時間との積(いわゆる管電流時間積(mA・s))によって決定される。   When the size of the subject P, the structure of the X-ray tube 101 and the X-ray detector 102, the reconstruction unit 106, etc. are constant and the energy of the irradiation X-ray is constant, the noise of the CT image is irradiated. It is determined by the X-ray dose, that is, the product of the tube current and the irradiation time in the X-ray tube 101 that is proportional to the X-ray dose (so-called tube current time product (mA · s)).

一方、CT画像のノイズは加法的であり、概ねガウス分布に従うことも知られている。すなわち、ある画素xの画素値v(x)を構成する任意のスカラー値v(x)(n=1,2,…,K)について、その真値(ノイズの寄与分を除去した値)をv (x)とすると、これらの差の値v(x)−v (x)は、概ね平均0、分散σ のガウス分布に従う(なお、照射X線量ないし管電流時間積m・Asとノイズの分散σ とは、概ね反比例関係にある。)。 On the other hand, it is also known that CT image noise is additive and generally follows a Gaussian distribution. That is, for an arbitrary scalar value v n (x) (n = 1, 2,..., K) constituting the pixel value v (x) of a certain pixel x, the true value (value obtained by removing the noise contribution). V n 0 (x), these difference values v n (x) −v n 0 (x) generally follow a Gaussian distribution with an average of 0 and a variance σ k 2 (note that the irradiation X-ray dose or the tube current The time product m · As and the noise variance σ k 2 are generally in inverse proportion.)

また、この分散σ は、画素xの位置そのもの(上で述べたように、例えば各座標値x=(x,y))にも依存するが、通常のX線CT装置100においては、X線管101及びX線検出器102の間に、X線照射量を調節する物理的なX線フィルタ(例えば銅箔や金属塊等により構成された、いわゆる「ウェッジ」あるいは「X線フィルタ」と呼称されるもの)を備えているため、これを無視することができる。なぜならばウェッジは、被検体Pが水とほぼ同じ密度を持つ物質から構成されていることを利用して、どのX線検出器102においても同程度のX線量が検出されるよう、照射されるX線量の一部減弱する作用を有するものであり、従ってこのようなウェッジによれば、結果的に、ノイズの分散σ を画素xの位置に殆ど依らない概ね一定値にする効果を生じるからである(ちなみに、このウェッジは、一般に、X線検出器102のダイナミックレンジを有効に利用することを本来の目的として設置されるものである。)。 The variance σ k 2 also depends on the position of the pixel x itself (as described above, for example, each coordinate value x = (x, y)), but in the normal X-ray CT apparatus 100, Between the X-ray tube 101 and the X-ray detector 102, a physical X-ray filter (for example, a so-called “wedge” or “X-ray filter” composed of a copper foil, a metal lump or the like) that adjusts the X-ray irradiation amount. So that it can be ignored. This is because the wedge is irradiated so that the same amount of X-rays can be detected by any X-ray detector 102 by using the fact that the subject P is made of a substance having the same density as water. The X-ray dose has a function of partially reducing the X-ray dose. Therefore, according to such a wedge, as a result, the noise variance σ k 2 is brought to an almost constant value almost independent of the position of the pixel x. (Incidentally, this wedge is generally installed for the purpose of effectively using the dynamic range of the X-ray detector 102).

以上のことから、ダイナミックCT撮影により取得されたK枚の静止画像上においては、k枚目の静止画像上におけるすべての画素について、分散σ がほぼ一定であると推定することは妥当である。むろん、画素ごとに分散が異なる場合について本実施形態を拡張することも容易に推考できる。 From the above, on the K still images acquired by dynamic CT imaging, it is reasonable to estimate that the variance σ k 2 is almost constant for all pixels on the kth still image. is there. Of course, it can be easily inferred that the present embodiment can be extended in the case where the variance is different for each pixel.

さて次に、上記(2)式を具体的に演算するためには、その分散σ として、どのような数値をあてるか、が問題となる。このようなことが問題となるのは、通常、ノイズの分布の形は想定できても(上記ではガウス分布)、分散σ の具体値は不明であることが多いからである。 Now, in order to specifically calculate the above equation (2), what value is assigned as the variance σ k 2 becomes a problem. This is a problem because usually the shape of the noise distribution can be assumed (Gaussian distribution in the above), but the specific value of the variance σ k 2 is often unknown.

更に、一般的に、毎回の撮影毎に照射線量(X線管電流×照射時間(mAs))を変更して撮影を行ってもよい。   Furthermore, in general, imaging may be performed by changing the irradiation dose (X-ray tube current × irradiation time (mAs)) for each imaging.

さて、k枚目の画像(k=1,2,…,K)に於いて各画素のスカラー値が持つノイズの分散をσ とし、k枚目の画像の撮影に用いた照射線量をRとするとき、σ はRに比例する。従って少なくともひとつのk=kについてσkO が指定できれば、他のkに関しても、

Figure 2008100121
In the k-th image (k = 1, 2,..., K), the noise dispersion of the scalar value of each pixel is σ k 2, and the irradiation dose used for photographing the k-th image is When R k , σ k 2 is proportional to R k . Therefore, if σk O 2 can be specified for at least one k = k 0 ,
Figure 2008100121

によってσ を正確に推定することができる。 Can accurately estimate σ k 2 .

本実施形態(このような事情が当てはまる)に於いては少なくともひとつのkについて、以下のような方法でσ の具体的数値の推定を行うことができる。 In the present embodiment (this situation applies), a specific numerical value of σ k 2 can be estimated for at least one k by the following method.

K回の撮影のうち、被検体Pに殆ど変化がなかったと仮定することのできるN回(1<N≦K)の画像を用いて、実測により、分散σ に対する期待値E[σ ]を求める方法が有効である。以下説明を簡単にするために、これらN枚の画像における照射線量は同じであり、従ってk=1,2,…Nに関してσ は一定(σと書く)と仮定する。これらN枚の画像における、ある画素xの画素値v(x)を構成する各スカラー値v(x),v(x),…,v(x)が含むノイズは、上述したように平均0、分散σのガウス分布に従うと予想されるから、これらの平均値を以下の(6)式、

Figure 2008100121
Of the K times of imaging, the expected value E [σ k for the variance σ k 2 is measured by using N times (1 <N ≦ K) images that can be assumed to have little change in the subject P. 2 ] is effective. For simplicity of explanation below, it is assumed that the irradiation doses in these N images are the same, and therefore σ k 2 is constant (denoted as σ 2 ) for k = 1, 2 ,. Noise included in each of the scalar values v 1 (x f ), v 2 (x f ),..., V K (x f ) constituting the pixel value v (x f ) of a certain pixel x f in these N images. Is expected to follow a Gaussian distribution with mean 0 and variance σ 2 as described above, so these average values are expressed by the following equation (6):
Figure 2008100121

を用いると、真の分散σに対する期待値E[σ]を、

Figure 2008100121
Is the expected value E [σ 2 ] for the true variance σ 2 ,
Figure 2008100121

として求めることができる。そして、この分散の期待値E[σ]は、上述した通り、K枚すべての静止画像上の全画素xにつき妥当するものと考えることができ、真の分散σの代用として用いるのに、一定程度以上確からしさが保証された値である。したがって、上記(4)式の実際の演算においては、このE[σ]を(4)式のσに代入すればよい。 Can be obtained as The expected value E [σ 2 ] of the variance can be considered to be appropriate for all pixels x on all K still images as described above, and can be used as a substitute for the true variance σ 2. This is a value for which the certainty is guaranteed over a certain level. Therefore, in the actual calculation of the above equation (4), this E [σ 2 ] may be substituted for σ 2 of equation (4).

なお、このようなE[σ]は、より具体的には、K枚の静止画像中、例えば1枚目と2枚目の静止画像に基づく実測値により求めてもよい(上記(6)及び(7)式で言えば、N=2とすることに該当する。)。また、上記(6)及び(7)式の実際の演算に供される画素xについては、例えば、空気や骨が撮像されている部分を除いた適当な画素xのみを選定する(複数選定した場合は得られるE[σ]すべての平均をとる)等といった工夫を施してもよい。さらに、その他一般的には、被検体Pの動きによる影響を抑える工夫等を施すと尚よい。 More specifically, such E [σ 2 ] may be obtained from measured values based on, for example, the first and second still images in K still images (the above (6)). And, in terms of equation (7), this corresponds to N = 2). Further, for the pixel x f to be subjected to the actual operation of the (6) and (7), for example, air or bone are selected only appropriate pixel x f excluding the portion being imaged (s If selected, an average of all the obtained E [σ 2 ] may be taken. Furthermore, in general, it is better to devise measures to suppress the influence of the movement of the subject P.

これらN枚の画像の撮影において照射線量が一定でない場合においても、σ がRに比例することを利用して正しくσ を推定することは容易に推考できるであろう。 Even when the irradiation dose is not constant in taking these N images, it can be easily estimated that σ k 2 is correctly estimated by using the fact that σ k 2 is proportional to R k .

さて次に、上記(4)式におけるパラメータCについての説明を行う。まず、(4)式においては、上記一般的形態で述べた危険率p(x,y)の考え方が、以下のようにして含まれている。すなわち、(4)式の右辺分子における根号内の表式は、いわゆるχ二乗分布に従うとされる当該χ値に一致するものであり、これを(2σ)で除し、括弧の全体をeの肩に置いた値は、危険率p1(x,y)そのものである。つまり、

Figure 2008100121
Next, the parameter C in the above equation (4) will be described. First, the expression (4) includes the concept of the risk factor p (x, y) described in the above general form as follows. That is, the expression in the root sign in the numerator on the right side of the equation (4) matches the χ 2 value that follows the so-called χ square distribution, which is divided by (2σ) 2 and the whole parenthesis Is the risk factor p1 (x, y) itself. That means
Figure 2008100121

そして、上記(4)式は、この(8)式のように表されるp1(x,y)に関し、

Figure 2008100121
And the above equation (4) relates to p1 (x, y) expressed as this equation (8).
Figure 2008100121

としたものに他ならない。尚、Aは定数でp1が(0〜1)の値になるように規格化されたものである。 It is nothing else. A is a constant and is standardized so that p1 has a value of (0 to 1).

結局、(4)式においては、上記したような一般的形態で述べた危険率p(x,y)が陽には表示されてはいないが、重みw1(p(x,y))の実態は、上述したように、まさしく危険率(=p1(x,y))の関数であると見ることができ((9)式)、すなわち「適合度の関数」である(ただし、危険率と適合度とは、上述したように、一方が増えれば他方も増加する関係にある)。   Eventually, in formula (4), the risk factor p (x, y) described in the general form as described above is not explicitly displayed, but the actual condition of the weight w1 (p (x, y)). Can be regarded as a function of the risk factor (= p1 (x, y)) as described above (equation (9)), that is, “function of fitness” (however, the risk factor and As described above, the degree of fitness has a relationship in which when one increases, the other increases.

そして、上記(9)式からわかるように、パラメータCは、重みw1(p(x,y))が、危険率p1(x,y)にどの程度敏感に反応するかを決める効果がある。つまり、Cを大きくすると、p1(x,y)がわずかに小さくなるだけで、w1(p(x,y))は0に近づく。また、Cを小さくするとそのような過敏な反応を抑制することができる。なお、Cとして、具体的には1乃至10程度とすればよく、好適にはC=3とするとよい。   As can be seen from the above equation (9), the parameter C has an effect of determining how sensitively the weight w1 (p (x, y)) reacts to the risk factor p1 (x, y). That is, when C is increased, p1 (x, y) is slightly decreased, and w1 (p (x, y)) approaches 0. Moreover, when C is made small, such a sensitive reaction can be suppressed. Specifically, C may be about 1 to 10 and preferably C = 3.

本実施形態においては、両画素x及びyに関する類似判定、言い換えると、両画素x及びyに関する上述した帰無仮説Hの棄却の判定は、上述したことから明らかなように、上記危険率p1(x,y)に基づいて、いわゆるχ二乗検定法(統計的検定法)によって決定されている。   In this embodiment, the similarity determination regarding both the pixels x and y, in other words, the determination of rejection of the above-described null hypothesis H regarding both the pixels x and y is, as is clear from the above, the risk factor p1 ( x, y) based on the so-called chi-square test method (statistical test method).

また、上記(4)式の表式からわかるように、本発明においては、危険率p(x,y)をx,yの組み合わせそれぞれについて計算した後、重みw(p(x,y))を求めるといった手順を踏む必要は必ずしもなく、危険率p(x,y)を具体的に求めずに、合成関数としての(wOp)を、直接計算する構成としてもよい。   In addition, as can be seen from the expression (4), in the present invention, after calculating the risk factor p (x, y) for each combination of x and y, the weight w (p (x, y)) It is not always necessary to go through the procedure of obtaining, and it may be configured to directly calculate (wOp) as a composite function without specifically obtaining the risk factor p (x, y).

以上述べたように、分散σの推定をし(例えば、(7)式のE[σ])、かつ、パラメータCを適当に決める(例えば、C=3)ことにより、(4)式を用いて、ある画素xにつき定義される集合N(x)(上述したように、例えば画素xを中心とした3×3画素分のエリア等)に含まれるすべての画素yについて、具体的な重みw1(p(x,y))を求めることができる。後は、上記(2)式におけるw(p(x,y))に代えて、このw1(p(x,y))を用いることにより、コヒーレントフィルタの具体的な数値演算を実施することが可能となる。そしてその結果、時間分解能は勿論のこと、空間分解能をも損なわずに、ノイズを強く抑制した画素値v´(x)=(v´(x),v´(x),…,v´(x))(=(3)式)、すなわちそのようなK枚の静止画像ないし動画像を、得ることができる。 As described above, the variance σ 2 is estimated (for example, E [σ 2 ] in the equation (7)), and the parameter C is appropriately determined (for example, C = 3). For all pixels y included in the set N (x) defined for a certain pixel x (for example, an area corresponding to 3 × 3 pixels centered on the pixel x, for example) The weight w1 (p (x, y)) can be obtained. After that, by using this w1 (p (x, y)) instead of w (p (x, y)) in the above equation (2), it is possible to carry out a specific numerical calculation of the coherent filter. It becomes possible. As a result, the pixel values v ′ (x) = (v ′ 1 (x), v ′ 2 (x),..., V, in which noise is strongly suppressed without impairing the spatial resolution as well as the temporal resolution. ' K (x)) (= Equation (3)), that is, such K still images or moving images can be obtained.

このような画像処理を、概念的に把握しやすいよう図示したものが、図3である。すなわちまず、図3(a)においては、1,2,…,K枚ある静止画像において、ある画素xにつき、該画素xを中心とした3×3画素分の矩形状エリアN3×3(x)が想定されている。この矩形状エリアN3×3(x)の左角隅における画素を、yとすれば、この画素yは、図3に併せて示すように、画素値v(y)を有している。 FIG. 3 illustrates such image processing so that it can be conceptually easily understood. That is, in FIG. 3A, in a still image of 1, 2,..., K still images, for a certain pixel x, a rectangular area N 3 × 3 (3 × 3 pixels centered on the pixel x). x) is assumed. The pixels in the left corner corner of the rectangular area N 3 × 3 (x), if y 1, the pixel y 1, as also shown in FIG. 3, a pixel value v (y 1) ing.

そして、この画素値v(y)を構成するスカラー値v(y),v(y),…,v(y)と画素値v(x)におけるスカラー値v(x),v(x),…,v(x)とのそれぞれにより、上記(4)式によって重みw1(p(x,y))が計算される(図3(b))。また、矩形状エリアN3×3(x)の残る画素y,…,yについても同様で、結局図3(b)に示すように、w1(p(x,y)),…,w1(p(x,y))及び、w1(p(x,x))が得られる。(この場合、(8)式より危険率p(x,x)は、“1”であり、したがって重みw1(p(x,x))も、(9)式より“1”である(=最大の重み付けがされている))。 Then, the scalar value v 1 constituting the pixel value v (y 1) (y 1 ), v 2 (y 1), ..., v K (y 1) and the scalar value of the pixel values v (x) v 1 ( x), v 2 (x),..., v K (x), the weight w1 (p (x, y 1 )) is calculated by the above equation (4) (FIG. 3B). The same applies to the remaining pixels y 2 ,..., Y 8 in the rectangular area N 3 × 3 (x), and as shown in FIG. 3B, w1 (p (x, y 1 )),. , W1 (p (x, y 8 )) and w1 (p (x, x)). (In this case, the risk factor p (x, x) is “1” from equation (8), and therefore the weight w1 (p (x, x)) is also “1” from equation (9) (= The maximum weight has been)).

次に、このようにして得られた、重みw1(p(x,y)),…,w1(p(x,y)),w1(p(x,x))を、対応する画素の、k枚目の画像におけるスカラー値v(y),v(y),…,v(y),v(x)にそれぞれ乗算して総和を取り(上記(2)式における分子に該当する。)、これを矩形状エリアN3×3(x)に関する重みw1の総和(同じく(2)式の分母に該当する。)により除せば、当該k枚目の画像における画素xについての、ノイズが抑制されたスカラー値v´(x)を求めることができる(図3(c))。また、k=1,2,…,Kのすべての画像につき、同じ重みw1(p(x,y)),…,w1(p(x,y)),w1(p(x,x))を用いて、ノイズが抑制されたスカラー値v´(x)を求めることによって、画素xにおけるノイズが抑制された画素値v´(x)=(v´(x),v´(x),…,v´(x))が得られる。すべての画素xにつき、上記演算を繰り返せば、ノイズを抑制したK枚の画像が得られる。 Next, the weights w1 (p (x, y 1 )),..., W1 (p (x, y 8 )), w1 (p (x, x)) obtained in this way are used as the corresponding pixels. Are multiplied by scalar values v k (y 1 ), v k (y 2 ),..., V k (y 8 ), v k (x) in the k-th image , respectively, to obtain the sum ((2 ), And this is divided by the sum of the weights w1 related to the rectangular area N 3 × 3 (x) (also corresponds to the denominator of equation (2)), the kth A scalar value v ′ k (x) in which noise is suppressed can be obtained for the pixel x in the image (FIG. 3C). Further, for all images of k = 1, 2,..., K, the same weights w1 (p (x, y 1 )),..., W1 (p (x, y 8 )), w1 (p (x, x) )) To obtain a scalar value v ′ k (x) in which noise is suppressed, thereby obtaining a pixel value v ′ k (x) = (v ′ 1 (x), v in which noise in the pixel x is suppressed. ' 2 (x), ..., v' K (x)) is obtained. If the above calculation is repeated for all the pixels x, K images with reduced noise can be obtained.

このようにしてコヒーレントフィルタで算出された画素値v´(x)で構成される画像では、オリジナル画像で見られたランダムなノイズが、十分に抑制される。   In the image composed of the pixel values v ′ (x) calculated by the coherent filter in this way, random noise seen in the original image is sufficiently suppressed.

なお、以上までに述べた各処理は、例えば図4に示すようなフローチャートに則ってこれを行えばよく、また、当該各処理に係る演算・画像表示等を実際のX線CT装置100上で実現するためには、例えば、図2に示すように、分散値推定部111、重み演算部112及び画素値演算部113により構成される画像処理部110を設けて、これを実施すればよい。   The above-described processes may be performed in accordance with a flowchart as shown in FIG. 4, for example, and calculations and image display related to the processes are performed on the actual X-ray CT apparatus 100. In order to realize this, for example, as shown in FIG. 2, an image processing unit 110 configured by a variance value estimation unit 111, a weight calculation unit 112, and a pixel value calculation unit 113 may be provided and implemented.

このうち重み演算部112は、上述した手順通り、画素値v(x)及びv(y)から直接重みw1(p(x,y))を求める構成となっている。したがって当該演算部112は、危険率p1(x,y)の値を具体的に求めることなく(すなわち、「危険率演算部(本発明にいう「適合度定量化部」を内蔵し)、重みを直接に求める装置である。なお、上記したような構成ではなく、具体的に危険率p1(x,y)の値を求める「危険率演算部(適合度定量化部)」と、その出力に基づいて重みw1(p(x,y))を求める「重み演算部」という、二段の手順を踏む構成としてもよい。いずれにせよ、重み演算部112は、分散値推定部111により推定された分散σと、v(x)及びv(y)を用いて重みw1(p(x,y))を算出する。 Among these, the weight calculation unit 112 is configured to obtain the weight w1 (p (x, y)) directly from the pixel values v (x) and v (y) as described above. Therefore, the calculation unit 112 does not specifically calculate the value of the risk factor p1 (x, y) (that is, the “risk rate calculation unit (incorporating the“ fitness quantification unit ”according to the present invention) and weights). It should be noted that the “risk rate calculation unit (fitness quantification unit)” that specifically calculates the value of the risk factor p1 (x, y) and the output thereof are not the above-described configuration. The weight calculation unit 112 may be configured to take a two-step procedure called “weight calculation unit” that calculates the weight w1 (p (x, y)) based on the above. The weight w1 (p (x, y)) is calculated using the variance σ 2 and v (x) and v (y).

また、画素値演算部113は、画素値v(x)及びv(y)、並びに重み演算部112により数値演算された重みw1(p(x,y))を使って、画素値v´(x)を演算する。すなわち当該演算部113は、元となる画像のノイズを抑制する処理、すなわちコヒーレントフィルタの適用を実際に行う(以下、これを「コヒーレントフィルタをかける」と表現する。)。   In addition, the pixel value calculation unit 113 uses the pixel values v (x) and v (y) and the weight w1 (p (x, y)) numerically calculated by the weight calculation unit 112 to generate a pixel value v ′ ( x) is calculated. That is, the calculation unit 113 actually performs processing for suppressing noise of the original image, that is, application of a coherent filter (hereinafter, this is expressed as “coherent filter is applied”).

上記のようなダイナミック・コヒーレントフィルタ処理においてK枚の静止画像から構成される動画像に、コヒーレントフィルタをかける場合には、上記画像処理部110における処理は、一旦すべての静止画像を再構成した後、これらを上記記憶装置10Mに蓄え、後処理として後にこれらに対してコヒーレントフィルタをかけるようにしてもよいが、本実施形態はこのような形態に限定されるものではなく、上述した連続スキャン、連続投影データ収集、連続再構成及び連続表示という流れの中で、コヒーレントフィルタをかける処理をリアルタイムに実施する(以下、これを「リアルタイム・コヒーレントフィルタ処理」と呼ぶ。)のでもよい。   In the case of applying a coherent filter to a moving image composed of K still images in the dynamic coherent filter processing as described above, the processing in the image processing unit 110 is performed after once reconstructing all still images. These may be stored in the storage device 10M and subjected to a coherent filter later as post-processing. However, the present embodiment is not limited to such a form. In the flow of continuous projection data collection, continuous reconstruction, and continuous display, a process for applying a coherent filter may be performed in real time (hereinafter referred to as “real time coherent filter process”).

リアルタイム・コヒーレントフィルタ処理の好ましい実施形態においては、新しい画像が撮影され再構成されるたびに、以下のような処理を行う。最初に得られた画像(画像番号1)から最新の画像(画像番号M)までのうち、画像番号M,M−1,…,M−K+1を持つK枚の静止画像上、共通の同一点(同一座標)の画素xの持つ画素値(スカラー値)を並べてK次元ベクトル値v(x)=(v(x),vM−1(x),…,vM−K+1(x))を構成する。こうして、上記の「ダイナミック・コヒーレントフィルタ処理」と全く同様にコヒーレントフィルタをかけることができる。ただし、画素値演算部113は実際には画素値v´(x)の全ての要素を計算するのではなく、最新の画像(画像番号M)に対応するスカラー値v´(x)だけを計算する。この結果、計算速度が向上するので、リアルタイムでノイズが抑制された最新の画像を表示できる。 In a preferred embodiment of real-time coherent filtering, each time a new image is taken and reconstructed, the following processing is performed. Among the first obtained image (image number 1) to the latest image (image number M), common identical points on K still images having image numbers M, M-1,..., M-K + 1. The pixel values (scalar values) of the pixel x at (same coordinates) are arranged and the K-dimensional vector value v (x) = (v M (x), v M−1 (x),..., V M−K + 1 (x) ). In this way, a coherent filter can be applied in exactly the same manner as the “dynamic coherent filter processing” described above. However, the pixel value calculation unit 113 does not actually calculate all the elements of the pixel value v ′ (x), but only the scalar value v M ′ (x) corresponding to the latest image (image number M). calculate. As a result, since the calculation speed is improved, the latest image in which noise is suppressed can be displayed in real time.

この「リアルタイム・コヒーレントフィルタ処理」の別の好ましい実施例として、最初のK枚の画像が得られた時点で、上記と全く同様にコヒーレントフィルタをかけてv´(x),…,v´(x)を求めておき、以後は、K次元ベクトル値を画像番号M,M−1,…,M−K+1を持つK枚の静止画像を用いてv(x)=(v(x),vM−1´(x),…,vM−K+1´(x))によって構成し、これに対して上記のリアルタイム・コヒーレントフィルタ処理を適用するように構成してもよい。なお、これらのリアルタイム・コヒーレントフィルタ処理の際に画素値ベクトルv(x)の次元Kを、マニュアル設定、あるいは自動設定によって、随時変更できるように構成しておくと便利である。 As another preferred embodiment of this “real-time coherent filtering”, when the first K images are obtained, a coherent filter is applied in the same manner as described above, and v 1 ′ (x) ,. ′ (X) is obtained, and thereafter, using K still images having K-dimensional vector values having image numbers M, M−1,..., M−K + 1, v (x) = (v M (x ), V M−1 ′ (x),..., V M−K + 1 ′ (x)), and the real-time coherent filter processing described above may be applied to this. It is convenient that the dimension K of the pixel value vector v (x) can be changed at any time by manual setting or automatic setting during the real-time coherent filter processing.

このようにコヒーレントフィルタにより、空間及び時間分解能を低下させることなく、ノイズだけを効果的に抑制したCT画像を使ってCBPスタディを実行し、脳組織(毛細血管)の血液の動態を定量的に解析し、そのインデックス(CBP、CBV、MTT、Err)を求めることにより、その精度及び信頼性の向上が期待できる。   In this way, CBP study is performed using CT images that effectively suppress only noise without reducing spatial and temporal resolution by using a coherent filter, and the dynamics of blood in brain tissue (capillary blood vessels) are quantitatively analyzed. By analyzing and obtaining the indexes (CBP, CBV, MTT, Err), improvement in accuracy and reliability can be expected.

以下に、CBPスタディについて説明する。
(CBPスタディ)
(原理)
CBPスタディでは、脳組織内の”毛細血管を通過する血流”の動態を定量的に表すCBP、CBV、MTT、Errのインデックスを求め、またこれらインデックスのマップを出力する。
Below, the CBP study will be described.
(CBP study)
(principle)
In the CBP study, indexes of CBP, CBV, MTT, and Err that quantitatively indicate the dynamics of “blood flow through capillaries” in brain tissue are obtained, and a map of these indexes is output.

CBP:脳組織の毛細血管内の単位体積及び単位時間あたりの血流量[ml/100ml/min]
CBV:脳組織内の単位体積あたりの血液量[ml/100ml]
MTT:毛細血管の血液平均通過時間[秒]
Err:解析モデルからの実測値のずれ残差の指標。この指標の多少により、脳動脈の支配組織と非支配組織との識別等の解析が可能となる。
CBP: Unit volume and blood flow per unit time in capillaries of brain tissue [ml / 100 ml / min]
CBV: Blood volume per unit volume in brain tissue [ml / 100ml]
MTT: Average blood passage time in capillaries [sec]
Err: An index of the residual deviation of the actual measurement value from the analysis model. Depending on the degree of this index, it is possible to analyze the discriminating between the dominating tissue and the non-dominating tissue of the cerebral artery.

CBPスタディでは、トレーサーとして脳血管透過性を持たない造影剤、たとえばヨード造影剤が使用される。ヨード造影剤は例えばインジェクターにより肘静脈から注入される。インジェクターにより静注されたヨード造影剤は、心臓、肺を経由して、脳動脈から流れ込む。そして、脳動脈から、脳組織内の毛細血管を経て、脳静脈へと流れ出ていく。このとき、脳血管透過性を持たない造影剤、たとえばヨード造影剤は正常な脳組織内の毛細血管では造影剤は血管外へ漏れ出ることなく通過する。   In the CBP study, a contrast agent having no cerebral vascular permeability, such as an iodine contrast agent, is used as a tracer. The iodine contrast agent is injected from the cubital vein by an injector, for example. The iodine contrast agent intravenously injected by the injector flows from the cerebral artery via the heart and lungs. Then, it flows out from the cerebral artery to the cerebral vein through the capillaries in the brain tissue. At this time, a contrast agent having no cerebral vascular permeability, for example, an iodinated contrast agent, passes through the capillary in normal brain tissue without leaking out of the blood vessel.

造影剤の通過の様子をダイナミックCTで連続的に撮影して、その連続画像から、脳動脈上の画素の時間濃度曲線Ca(t)、毛細血管を含む脳組織上の画素の時間濃度曲線Ci(t)、脳静脈上の画素の時間濃度曲線Csss(t)をそれぞれ測定する。   The state of passage of the contrast agent is continuously photographed by dynamic CT, and the time density curve Ca (t) of the pixel on the cerebral artery and the time density curve Ci of the pixel on the brain tissue including the capillary are obtained from the continuous image. (T) The time density curve Csss (t) of the pixel on the cerebral vein is measured.

CBPスタディでは、造影剤の血中濃度について脳組織に近い脳血管の血中濃度の時間曲線Ca(t)と、毛細血管の血中濃度の時間曲線Ci(t)との間で成り立つ理想的な関係を解析モデルとしており、つまり脳組織に入る直前の血管から造影剤を注入した場合、毛細血管を含む脳組織単位体積(1画素)内の時間濃度曲線は立ち上がりが垂直で、若干の勾配を持って立ち下がる形になる。これは、矩形関数で近似することができる(box−MTF法:box−Modulation Transfer Function method)。   In the CBP study, the blood concentration of the contrast agent is ideal between the time curve Ca (t) of the blood concentration of the cerebral blood vessels close to the brain tissue and the time curve Ci (t) of the blood concentration of the capillary blood vessels. When the contrast medium is injected from the blood vessel immediately before entering the brain tissue, the time density curve in the brain tissue unit volume (1 pixel) including the capillary blood vessel has a vertical rise and a slight gradient. It will be in the form of falling with. This can be approximated by a rectangular function (box-MTF method: box-Modulation Transfer Function method).

脳動脈血中時間濃度曲線Ca(t)を入力関数、脳組織の時間濃度曲線Ci(t)を出力関数として、毛細血管を通過する過程の特徴を、矩形関数で表される伝達関数として求めることができる。   Using the cerebral artery blood time concentration curve Ca (t) as the input function and the brain tissue time concentration curve Ci (t) as the output function, the characteristics of the process of passing through the capillaries are obtained as a transfer function represented by a rectangular function. Can do.

(具体的な手順)
図5、図6には、本実施形態によるCBPスタディの典型的手順を示している。
まず、肘静脈等の血管にボーラスインジェクション(造影剤を一気に投与する)を行い、その直後あるいは直前からダイナミックCT(同じ箇所を反復して撮影する)を行う。最も典型的な手技として、肘静脈へボーラスインジェクションを行った場合、概ね20〜40間に渡って、例えば0.5〜2秒間隔で撮影を繰り返す。ダイナミックCTで得たN枚のCT画像のうちのj枚目の各ピクセル(x、y)のCT値をv(x、y、j)とする。これはこの画素(x、y)における時間濃度曲線(滑らかな曲線である)f(t、x、y)をサンプリングしたものに他ならない。
(Specific steps)
5 and 6 show a typical procedure of the CBP study according to the present embodiment.
First, bolus injection (contrast medium is administered all at once) is performed on a blood vessel such as an elbow vein, and dynamic CT (images of the same part are repeatedly taken) is performed immediately after or immediately before. As a most typical procedure, when bolus injection is performed on the elbow vein, imaging is repeated at intervals of, for example, 0.5 to 2 seconds over a period of approximately 20 to 40 seconds. The CT value of each j-th pixel (x, y) in N CT images obtained by dynamic CT is represented by v (x, y, j). This is nothing but a sample of the time density curve (smooth curve) f (t, x, y) at this pixel (x, y).

まず、前処理として、ステップS1で、CT画像各々から、明らかに脳組織以外の組織であることが判別される画素を、解析対象から除外する。すなわち、脳組織のCT値として考えられる範囲(例えばCT値10〜60HU)に入らない値を示す画素は、空気や骨、脂肪などに対応する画素であり、脳血流の定量とは関係ないのでこれらは無視して良い。この解析範囲は、デフォルトとして、10〜60HUに設定されるが、入力部109を介して任意に設定可能である。   First, as preprocessing, in step S1, pixels that are clearly identified as tissues other than brain tissue are excluded from the analysis target from each CT image. That is, a pixel indicating a value that does not fall within a range that can be considered as a CT value of brain tissue (for example, a CT value of 10 to 60 HU) is a pixel corresponding to air, bone, fat, or the like, and is not related to cerebral blood flow quantification. So these can be ignored. This analysis range is set to 10 to 60 HU as a default, but can be arbitrarily set via the input unit 109.

また、前処理として、ステップS2で、造影効果の初期化が行われる。各画素に於ける造影効果(CT値の上昇)を得るためには、各画素(x,y)について、その画素に対応する組織に造影剤が到達する以前の画像(一般に複数枚得られる)を、通し番号1,2,…Kで表すと、その時間的平均値は、

Figure 2008100121
Further, as preprocessing, the contrast effect is initialized in step S2. In order to obtain a contrast effect (an increase in CT value) in each pixel, for each pixel (x, y), an image before the contrast agent reaches the tissue corresponding to that pixel (generally a plurality of images are obtained). Is represented by serial numbers 1, 2,... K, the temporal average value is
Figure 2008100121

を求め、この値をb(x,y)とする。そして、j=K+1、K+2、…,Nの各画像の画素値v(x,y,j)について、
q(x,y,j)=v(x,y,j)−b(x,y)
j<Kについて
q(x,y,j)=0
とすればよい。処理を簡単にするためには、どの画素に関しても同じKを採用しても良い。こうして得られたq(x,y,j)は、造影効果の(滑らかな)時間濃度曲線q(t,x,y)をt=t1,t2,…tNにおいてサンプリングしたものに他ならないと考えることができる。このq(t,x,y)を用いて脳血流の定量解析を行う。
And this value is b (x, y). For pixel values v (x, y, j) of each image of j = K + 1, K + 2,.
q (x, y, j) = v (x, y, j) -b (x, y)
About j <K
q (x, y, j) = 0
And it is sufficient. In order to simplify the processing, the same K may be adopted for any pixel. The q (x, y, j) obtained in this way is considered to be nothing but the sampled (smooth) time density curve q (t, x, y) of the contrast effect at t = t1, t2,. be able to. Quantitative analysis of cerebral blood flow is performed using this q (t, x, y).

定量解析にあたってはまず、右脳エリアと左脳エリアをCT画像上で分離する必要がある。上述したようにCBPスタディでは、毛細血管の血流動態の様子を、脳動脈の時間濃度曲線Ca(t)に対する脳組織の時間濃度曲線Ci(t)の伝達関数MTFとして求めるものであり、従って、解析対象の脳組織が、Ca(t)の脳動脈の支配下にがあることが前提とされる。少なくとも左脳と右脳とでそれぞれ別々の脳動脈の時間濃度曲線Ca(t)を使って個別に解析する、つまり左脳の脳動脈の時間濃度曲線Ca(t)は同じ左脳の脳組織の解析にだけ使用し、動揺に、右脳の脳動脈の時間濃度曲線Ca(t)は同じ右脳の脳組織の解析にだけ使用することが必要とされる。   In the quantitative analysis, it is necessary to first separate the right brain area and the left brain area on the CT image. As described above, in the CBP study, the blood flow dynamic state of the capillary is obtained as the transfer function MTF of the brain tissue time concentration curve Ci (t) with respect to the cerebral artery time concentration curve Ca (t). It is assumed that the brain tissue to be analyzed is under the control of Ca (t) cerebral artery. At least the left brain and the right brain are separately analyzed using the time concentration curves Ca (t) of the separate cerebral arteries. That is, the time concentration curve Ca (t) of the left cerebral artery is used only for analyzing the brain tissue of the same left brain. In use, the right brain cerebral artery time concentration curve Ca (t) needs to be used only for analysis of the same right brain brain tissue.

脳を左脳エリアと右脳エリアとを分割するために、図7に示すように、CT画像上に分割線が画面上に図形として重ねて表示される(S3)。分割線が最初は画像中央に表示されるように構成しても良い。操作者は、画像を参照して、分割線を移動し、また分割線を構成する複数の構成点を移動して任意に屈曲させることにより、左右エリアを分割する。   In order to divide the brain into the left brain area and the right brain area, as shown in FIG. 7, a dividing line is displayed on the screen as a figure superimposed on the CT image (S3). The dividing line may be initially displayed at the center of the image. The operator refers to the image, moves the dividing line, and moves the plurality of constituent points constituting the dividing line to arbitrarily bend, thereby dividing the left and right areas.

次に、CT画像上で脳動脈上に脳動脈ROIを設定するのであるが、この設定の精度を向上し且つ容易にするためにROI設定支援部121により支援マップが作成され、CT画像とは別に、またはCT画像に重ねて表示される(S4)。支援マップとしては、例えば、AT(アピアランスタイム)マップ、PT(ピークタイム)マップ、TT(トランジットタイム)マップがあげられる。各画素について、図8に示すように、造影前の任意の時刻(例えばデータ収集開始時刻)T0から造影剤濃度がピークpeakの数パーセント(例えば1パーセント)に達する時刻までの時間AT、時刻T0から造影剤濃度がピークに達した時刻までの時間(ピークタイム)PT、又は造影剤の移動時間を例えば半値幅で表すTTが計算され、マップとして生成され表示される。デフォルトでは、これらATマップ、PTマップ、TTマップの全種類が生成され表示されるようになっているが、任意の1種類、又は任意の2種類を操作者が選択することが可能である。   Next, the cerebral artery ROI is set on the cerebral artery on the CT image. In order to improve and facilitate the accuracy of this setting, a support map is created by the ROI setting support unit 121. Separately or overlaid on the CT image (S4). Examples of the support map include an AT (appearance time) map, a PT (peak time) map, and a TT (transit time) map. For each pixel, as shown in FIG. 8, time AT, time T0 from an arbitrary time before contrast (for example, data collection start time) T0 to a time when the contrast agent concentration reaches several percent (for example, 1%) of the peak peak The time (peak time) PT from the time until the contrast agent concentration reaches the peak, or TT that represents the moving time of the contrast agent, for example, by a half-value width, is calculated and generated and displayed as a map. By default, all types of the AT map, PT map, and TT map are generated and displayed, but the operator can select any one type or any two types.

脳動脈のこれら数値は他と比較して高いある程度固有値で現れる傾向にあるので、その値を中心としたウインドウ内が抽出して表示されるよう設定されたカラーマップ(ルックアップテーブル)を通してカラー表示させることで、脳動脈の場所を容易に識別し、脳動脈ROIを正確に設定することが可能となる(S5)。典型的には、左右脳エリアそれぞれに、脳動脈ROIは、前大脳動脈(ACA)、中大脳動脈(MCA)、後大脳動脈(PCA)の3箇所ずつ設定される。   These values of cerebral arteries tend to appear as eigenvalues to some extent higher than others, so color display through a color map (look-up table) set to extract and display the window centered on that value By doing so, the location of the cerebral artery can be easily identified and the cerebral artery ROI can be accurately set (S5). Typically, in each of the left and right brain areas, three cerebral arteries ROI are set: an anterior cerebral artery (ACA), a middle cerebral artery (MCA), and a posterior cerebral artery (PCA).

なお、マルチスライス、撮影の場合、例えば隣接する4枚のスライスを解析対象とする場合、図12に示すように、その各スライスで個々に脳動脈ROIを設定することは作業負担が大きいばかりで、解析を行う上では必要のない作業である。従って、ある任意の1スライスで設定した脳動脈ROIを他のスライスにも共用する。または、後述するコヒーレントレグレッション法を使って全スライスで共通に用いることのできる脳動脈の時間濃度曲線Ca(t)を作成するようにしてもよい。   In the case of multi-slice imaging, for example, when four adjacent slices are to be analyzed, setting the cerebral artery ROI individually for each slice as shown in FIG. This is not necessary for analysis. Therefore, the cerebral artery ROI set in one arbitrary slice is also shared with other slices. Alternatively, a cerebral artery time concentration curve Ca (t) that can be commonly used in all slices may be created using a coherent regression method described later.

次に、設定された脳動脈ROI各々に関して時間濃度曲線Ca(t)が、ダイナミックCTによる連続画像データから時間濃度曲線作成部122により作成される(S6)。   Next, a time density curve Ca (t) is created for each set cerebral artery ROI by the time density curve creating unit 122 from continuous image data by dynamic CT (S6).

ここで、画素に比べて、脳動脈は非常に細いものが多く、しかも一般にCTの撮影スライスに対して直交していないため、画像上のどの画素も正確に動脈血のCT値を表しておらず、1画素が脳動脈と他の組織との混在により構成され、そのパーシャルボリューム効果のためにそれよりも低い造影効果しか示さないことがほとんどである。また、動脈をパーシャルボリュームとして含むこれらの画素において、画像ノイズが大きい。特に脳梗塞を生じている部位の動脈等においては、造影効果が比較的小さいために、ノイズの影響は甚大である。画像ノイズに関しては上述したコヒーレントフィルタにより抑制されているものの、パーシャルボリューム効果の影響は依然として残存している。   Here, many cerebral arteries are much thinner than pixels, and since they are generally not orthogonal to the CT imaging slice, none of the pixels on the image accurately represents the CT value of arterial blood. One pixel is composed of a mixture of cerebral arteries and other tissues, and because of its partial volume effect, it usually exhibits a lower contrast effect. Further, image noise is large in these pixels including the artery as a partial volume. In particular, in an artery or the like at a site where cerebral infarction occurs, the effect of noise is enormous because the contrast effect is relatively small. Although the image noise is suppressed by the above-described coherent filter, the influence of the partial volume effect still remains.

この問題は、脳動脈の時間濃度曲線は、単一のスライス画像において計測するのではなく、その動脈を含む立体内の画素を用いて後述するコヒーレントレグレッション法を適用することで抑圧することが可能である。従って、上述したコヒーレントフィルタ法に代えて、この段階でコヒーレントレグレッション法を適用するようにしてもよい。   This problem can be suppressed by applying a coherent regression method to be described later using pixels in a solid including the artery rather than measuring the time density curve of the cerebral artery in a single slice image. Is possible. Therefore, instead of the above-described coherent filter method, a coherent regression method may be applied at this stage.

また、この方式によれば、動脈ごとにそれに対応する唯一の脳動脈の時間濃度曲線が得られ、従って、撮影範囲内にある全スライス中の任意の部位の解析に利用することができ、これによって、特定の動脈について、その脳動脈の時間濃度曲線が最も明瞭に得られるスライスを選んで、その脳動脈の時間濃度曲線を全スライスに適用することができ、脳動脈の時間濃度曲線の数を減らすことができる。   In addition, according to this method, for each artery, a cerebral artery time density curve corresponding to it can be obtained. Therefore, it can be used for analysis of any part in all slices within the imaging range. Allows you to select the slice for which the time density curve of the cerebral artery is most clearly obtained for a particular artery, and apply the time density curve of that cerebral artery to all slices. Can be reduced.

(コヒーレントレグレッション法)
上記時間濃度曲線の作成では、パーシャルボリューム効果、ランダムノイズの影響を除去することが重要である。まず、「時間濃度曲線」とは、上記ダイナミックCT画像中の特定の部位における像の濃度値の経時的変化を表す曲線である。ことに、上記医用画像診断装置においては、人体組織等における血流動態や代謝機能等の詳細を調べる事を目的として、人体の特定組織内の造影剤濃度等の経時的変化を時間濃度曲線として計測することが行われている。また、天体観測等においては、特定の天体の光度変化等を解析する目的で、時間濃度曲線が用いられる。より形式的に明示すると、すなわち、時間濃度曲線とは、時刻tにおけるある部位の濃度値をdとするとき、対の列{<t,d>(k=1,2,・・・,K)}として表現される。また、時間濃度曲線の多くの用途においては、必ずしもdの絶対的な値が必要なのではなく、むしろ最初の画像1を基準とする増分(d−d)だけが得られれば十分である。さらにそのような用途のうちの多くでは、単に(d−d)に比例するデータ A(d−d)(ここにAは未知の比例係数)だけが得られれば十分である。この場合には、従って、対の列{<t,A(d−d)>(k=1,2,・・・,K)}が、求める時間濃度曲線である。
(Coherent regression method)
In creating the time concentration curve, it is important to remove the effects of the partial volume effect and random noise. First, the “time density curve” is a curve representing a change with time of an image density value at a specific portion in the dynamic CT image. In particular, in the above-mentioned medical image diagnostic apparatus, the temporal change of the contrast medium concentration in a specific tissue of a human body is used as a time concentration curve for the purpose of examining details such as blood flow dynamics and metabolic functions in the human tissue. Measuring is done. In astronomical observation or the like, a time concentration curve is used for the purpose of analyzing a change in luminous intensity of a specific celestial object. And more formally specified, i.e., the time-density curve, when the density value of a site at time t k and d k, the pair row {<t k, d k> (k = 1,2, · .., K)}. Also, for many applications of time density curves, it is not always necessary to have an absolute value of d k , but rather it is sufficient if only the increment (d k −d 1 ) relative to the first image 1 is obtained. is there. Moreover many of such applications, (here A unknown proportional coefficient) simply (d k -d 1) proportional to the data A (d k -d 1) is sufficient as long obtained only. In this case, therefore, the pair of columns {<t k , A (d k −d 1 )> (k = 1, 2,..., K)} is the time density curve to be obtained.

このような時間濃度曲線を求めるためには、原理的には、上記ダイナミックCT画像を構成する各画像k(k=1,2,・・・,K)における、該時間濃度曲線を測定しようとする部位に含まれる画素xのスカラー値v(x)を用いて、対の列{<t,v(x)>}あるいは、{<t,A(v(x)−v(x))>}を構成すればよい。 In order to obtain such a time density curve, in principle, an attempt is made to measure the time density curve in each image k (k = 1, 2,..., K) constituting the dynamic CT image. Using the scalar value v k (x) of the pixel x included in the region to be paired, {<t k , v k (x)>} or {<t k , A (v k (x) −v 1 (x)) >>.

しかし、実用においては、上記医用画像診断装置等によって撮影されたダイナミックCT画像にランダムなノイズが含まれているために、本来測定しようとする時間濃度曲線を正確に求められないという問題がある。   However, in practice, there is a problem that the time density curve to be originally measured cannot be obtained accurately because random noise is included in the dynamic CT image taken by the medical image diagnostic apparatus or the like.

さらに、実用においては、これらのダイナミックCT画像においては、いわゆる「パーシャルボリューム効果」が生じる。パーシャルボリューム効果とは、すなわち、被検体内の微小な物体の像は、画像上では少数個の画素によって表現されるが、これら少数個の画素には、被検体内の隣接する物体の像も影響を与えるため、これら少数個の画素の画素値は(本来計測しようとする濃度値の変動に比例するものの)比較的小さな変動しか示さない、という現象である。言い換えれば、これら少数個の画素の画素値は僅かな信号しか含まない。従って、パーシャルボリューム効果が生じている場合には、どの画素xを取っても対の列{<t,v(x)>(k=1,2,・・・,K)}は非常に信号レベルが低く、本来計測しようとしているのではない組織における濃度値の変化の影響を受け、さらにランダムなノイズが存在するために、本来測定しようとする時間濃度曲線{<t,d>}を正確に求められないという問題がある。 Furthermore, in practical use, a so-called “partial volume effect” occurs in these dynamic CT images. The partial volume effect means that an image of a minute object in a subject is represented by a small number of pixels on the image, and an image of an adjacent object in the subject is also included in these few pixels. This is a phenomenon in which the pixel values of these small number of pixels exhibit relatively small fluctuations (although they are proportional to the fluctuations of the density value to be measured). In other words, the pixel values of these few pixels contain only a few signals. Therefore, when the partial volume effect is generated, the pair of columns {<t k , v k (x)> (k = 1, 2,..., K)} is extremely high regardless of which pixel x is taken. Since the signal level is low and affected by the change of the concentration value in the tissue that is not originally measured, and there is random noise, the time concentration curve to be originally measured {<t k , d k There is a problem that it is not possible to accurately obtain>

そこで、従来は、ランダムなノイズを抑制するために、時間的又は空間的平滑化を用いていたが、時間平均を行うと時間分解能が損なわれてしまい、また、空間平均を行うと、本来の該測定しようとする部位以外の部位の濃度の経時変化が計測値に混入するという問題点があった。このような問題点を解決し、より正確な時間濃度曲線を得るために、コヒーレントフィルタを採用する。   Therefore, in the past, temporal or spatial smoothing has been used to suppress random noise, but temporal resolution is lost when temporal averaging is performed, and when spatial averaging is performed, There has been a problem that a change with time in the concentration of a portion other than the portion to be measured is mixed into the measured value. In order to solve such a problem and obtain a more accurate time density curve, a coherent filter is employed.

まず、本実施形態のコヒーレントフィルタにおいて用いるべき、帰無仮説について説明する。計測しようとする部位における真の時間濃度曲線を{<t,d>(k=1,2,・・・,K)}であると仮定するとき、その一次変換である{<t,A(d−d)>(k=1,2,・・・,K)}(だたしAは未知の係数)を計測することを目的とする場合において、計測しようとする部位に概ね相当する画素の集合Rを設定する。この集合Rの要素である任意の画素x∈Rについて、条件Q:「もし、この画素xが上記真の時間濃度曲線を良く反映し、しかも他の部位の経時的濃度変化の影響をほとんど受けていない」のであれば、(ベクトル値としての)画素値v(x)=(v(x),v(x),...,v(x))について、パーシャルボリューム効果およびランダムノイズの影響を考慮することによって、
(x)=p(x) d+q(x)+γ(x) …(11)
(k=1,2,・・・,K)
が成り立つと仮定することができる。ここに、p(x)およびq(x)は、画素xごとに異なるが画像番号k(すなわち撮影時刻t)によっては変化しない未知の係数であり、パーシャルボリューム効果をモデル化したものである。またγ(x)はランダムなノイズをモデル化したものであって、画素xごとに、しかも画像番号kごとに値が異なるが、その期待値は0であり、またその統計分布は画素xにも画像番号kにも依存しない。
First, the null hypothesis to be used in the coherent filter of this embodiment will be described. When it is assumed that the true time concentration curve at the site to be measured is {<t k , d k > (k = 1, 2,..., K)}, {<t k , A (d k −d 1 )> (k = 1, 2,..., K)} (where A is an unknown coefficient), the part to be measured A set R of pixels substantially corresponding to is set. For any pixel x∈R that is an element of the set R, the condition Q: “If this pixel x well reflects the true time density curve, it is almost affected by changes in density over time in other parts. If not, the partial volume effect and random for pixel values v (x) = (v 1 (x), v 2 (x),..., V K (x)) (as vector values) By considering the effects of noise,
v k (x) = p (x) d k + q (x) + γ k (x) (11)
(K = 1, 2,..., K)
Can be assumed to hold. Here, p (x) and q (x) are unknown coefficients that are different for each pixel x but do not change depending on the image number k (that is, the photographing time t k ), and model the partial volume effect. . Γ k (x) is a model of random noise, and has a value different for each pixel x and for each image number k, but its expected value is 0, and its statistical distribution is the pixel x. And image number k.

以上の仮定によれば、該集合Rの要素である任意の2個の画素x,yに関して、もし「画素x,yが共に(上記の)条件Qを満たす。」という命題が成り立つのであれば、次式の関係が成り立つことが証明できる。   According to the above assumption, for any two pixels x and y that are elements of the set R, if the proposition “both pixels x and y satisfy the condition Q (above)” holds. It can be proved that the relationship of the following equation holds.

(x)=a(y)+a+ξ (k=1,2,・・・,K)
… (12)
ここに、aおよびaは、画素の組x,yごとに異なるが画像番号k(すなわち撮影時刻t)によっては変化しない未知の係数である。またξはランダムなノイズであって、画素の組x,yごとに、しかも画像番号kごとに値が異なるが、その期待値は0である。
v k (x) = a 1 v k (y) + a 2 + ξ k (k = 1, 2,..., K)
(12)
Here, a 1 and a 2 are unknown coefficients that differ for each pixel set x and y but do not change depending on the image number k (ie, the photographing time t k ). Ξ k is random noise, and the value is different for each pixel set x, y and for each image number k, but the expected value is zero.

(12)式は以下のようにして導かれる。すなわち、xにyを代入して得られる式
(y)= p(y) d+q(y)+γ(y)… (13)
を変形すると、

Figure 2008100121
Equation (12) is derived as follows. That is, an expression obtained by substituting y for x
v k (y) = p (y) d k + q (y) + γ k (y) (13)
Transforming
Figure 2008100121

とおくことによって、(12)式が導かれる。ここで、(16)式のaとaは、パーシャルボリューム効果を現すパラメータであり、また(16)式のξは、ランダムなノイズを表す。 By doing so, the equation (12) is derived. Here, a 1 and a 2 in the equation (16) are parameters representing the partial volume effect, and ξ k in the equation (16) represents random noise.

以上から、「画素x,yが共に条件Qを満たす。」という命題は、帰無仮説H´「v(x)=a(y)+a+ξ (k=1,…,K)である。」と等価であることが示された。 From the above, the proposition that “pixels x and y both satisfy condition Q” is the null hypothesis H 0 ′ “v k (x) = a 1 v k (y) + a 2 + ξ k (k = 1,... , K) ”.

次に帰無仮説H´「v(x)=a(y)+a+ξ (k=1,…,K)である。」を、実質的に等価であり、かつ実際に検定できる形式の命題に変換する方法について述べる。この帰無仮説を改めて数学的に厳密な表現で述べると、帰無仮説H´「ある定数aおよびaが存在して、ξ=v(x)−a(y)−a(k=1,…,K)は平均0、分散(σh(a))の正規分布に従う。」となる。ここに係数h(a)は、
h(a)=1+a … (17)
である。((17)式はaとξの定義である(16)式、および、ランダム変数に関する分散の持つ一般的な性質から直ちに導かれる。)また、上記の分散σの値は、簡便かつ実用上十分正確に推定できる。
Next, the null hypothesis H 0 ′ “v k (x) = a 1 v k (y) + a 2 + ξ k (k = 1,..., K)” is substantially equivalent and actually. Describes how to convert the proposition to a form that can be tested. When this null hypothesis is described in mathematically exact expression, the null hypothesis H 0 ′ “There are certain constants a 1 and a 2 and ξ k = v k (x) −a 1 v k (y ) −a 2 (k = 1,..., K) follows a normal distribution with an average of 0 and variance (σ 2 h (a 1 )). Here, the coefficient h (a 1 ) is
h (a 1 ) = 1 + a 1 2 (17)
It is. (Equation (17) is immediately derived from equation (16), which is the definition of a 1 and ξ k , and the general nature of the variance with respect to random variables.) The value of the variance σ 2 is simple And it can be estimated sufficiently accurately for practical use.

以上から、もし、上記の定数aおよびaを決定することができれば、上記の帰無仮説H´を検定することが可能である。そして実際上は、これらの定数の最適な推定値a およびa が得られれば十分である。 From the above, if the above constants a 1 and a 2 can be determined, it is possible to test the null hypothesis H 0 ′. In practice, it is sufficient to obtain optimal estimates a 1 to a 2 to these constants.

このような、定数aおよびaの最適な推定値の算出には、公知の当てはめ法(fitting)がそのまま利用できる。そこで、以下では、そのような当てはめ法の典型的な具体例として、線形最小二乗法を用いる場合における概要を説明する。線形最小二乗法を本実施例に適用するには、単に、上記の帰無仮説のξの二乗和をS(a)として、すなわち

Figure 2008100121
A known fitting method (fitting) can be used as it is for calculating the optimum estimated values of the constants a 1 and a 2 . Therefore, an outline in the case of using the linear least square method will be described below as a typical example of such a fitting method. In order to apply the linear least square method to this embodiment, the sum of squares of ξ k of the above null hypothesis is simply set as S (a), that is,
Figure 2008100121

を定義する。S(a)の値は定数ベクトルa=(a,a)、すなわち上記の定数aおよびa、の値に依存する。このS(a)が最小の値を取るような定数ベクトルaを算出すれば、定数aおよびaに関する、不偏推定の意味での最適な推定値a およびa が得られる。なお、線形最小二乗法の具体的な計算方法としては、様々な公知の方法を利用することができ、しかも、これら公知の計算方法はいずれも非常に簡単であり、必要な計算時間はごく僅かである。 Define The value of S (a) depends on the constant vector a = (a 1 , a 2 ), that is, the values of the above constants a 1 and a 2 . Be calculated constant vector a as the S (a) takes a minimum value, to the constant a 1 and a 2, the optimum estimated values a 1 ~ and a 2 ~ of the sense of unbiased estimate can be obtained. In addition, as a specific calculation method of the linear least square method, various known methods can be used, and these known calculation methods are all very simple and require a very short calculation time. It is.

このようにして、上記の定数a,aの最適な推定値a ,a を算出した結果、次式で定義される残差
(x,y)=v(x)−a (y)−a …(19)
を具体的に計算することができる。従って、この残差r を用いて、上記の帰無仮説H´を、実質的に等価な帰無仮説H”「r (x,y)(k=1,…,K)は平均0、分散(1+(a )σの正規分布に従う。」と言い換えることができる。これは、実際に検定の計算を実行可能な具体的命題である。
Thus, as a result of calculating the optimal estimated values a 1 to a 2 to the constants a 1 and a 2 , residuals defined by the following equations are calculated.
r k to (x, y) = v k (x) −a 1 to v k (y) −a 2 to (19)
Can be calculated specifically. Therefore, using this residual r k ˜ , the above-described null hypothesis H 0 ′ is substantially equivalent to the null hypothesis H 0 ″ “r k ˜ (x, y) (k = 1,..., K ) mean 0, variance (1+ (a 1 ~) 2 ) σ 2 of normally distributed. "and can be translated. This is a concrete proposition that can actually perform the calculation of the test.

なお、さらに、ベクトルによる表現

Figure 2008100121
In addition, expression by vector
Figure 2008100121

(ただし、ベクトルa及びξは画素の組x,yに依存する。)を導入し、また、次式
f(a~,v(y))=a v(y)+a … (21)
で定義されるベクトル値関数fを用いて帰無仮説H´を言い換えると、帰無仮説H”は「v(x)=f(a~,v(y))+ξ、(ただし、ξは平均0、分散(1+(a )σの正規分布に従う。)」となり、これは上述した帰無仮説Hと全く同じ形式である。すなわち、本実施形態は上述したコヒーレントフィルタの一変形例であることは明らかである。なお、ここで、上記f(a~,v(y))とは、すなわち、画素yの画素値v(y)に対して、パーシャルボリューム効果を現すパラメータaを最適に調節して、画素xの画素値v(x)と最も高い適合度を持つように変換したものを意味する。
(However, the vectors a and ξ depend on the pixel set x, y.) And the following expression f (a ~ , v (y)) = a 1 ~ v (y) + a 2 ~ ... ( 21)
In other words the null hypothesis H 0 'using defined as a vector value function f, the null hypothesis H 0 "is" v (x) = f (a ~, v (y)) + ξ, ( provided that, xi] mean 0, variance (1+ (a 1 ~) 2 ) follows the normal distribution of the sigma 2.) ", and this is exactly the same format as the null hypothesis H 0 described above. That is, it is clear that this embodiment is a modification of the above-described coherent filter. Here, the f (a ~, v (y )) and is, i.e., the pixel value v (y) of the pixel y, optimally adjusted to the parameter a indicating the partial volume effect, a pixel x The pixel value v (x) is converted so as to have the highest matching degree.

次に、本実施形態において、上記の帰無仮説H”を用いて、コヒーレントフィルタによって時間濃度曲線を求める方法について説明する。計測しようとする部位に概ね相当する画素の集合Rについて、この集合Rに含まれるあるひとつの画素x∈Rについて、集合Rの要素である全ての画素y∈Rに対して、以下の計算を行う。すなわち、上記の方法を用いて実際に残差r (x,y)(k=1,…,K)を算出し、次に、上記の帰無仮説H”「r (x,y)(k=1,…,K)は平均0、分散(1+(a )σの正規分布に従う。」を棄却する場合の危険率p(x,y)ないし重みw(p(x,y))を具体的に計算する。そして、重み付き平均v´(x)を下式(22)によって計算し、画素xにおける時間濃度曲線{<t,v´(x)−v´(x)>(k=1,2,・・・,K)}を構成する。

Figure 2008100121
Next, in the present embodiment, a method for obtaining a time density curve by a coherent filter using the above-described null hypothesis H 0 ″ will be described. For a set R of pixels roughly corresponding to a region to be measured, this set For one pixel xεR included in R, the following calculation is performed for all the pixels yεR that are elements of the set R. That is, using the above-described method, the residual r k ~ (X, y) (k = 1,..., K) is calculated, and then the above-described null hypothesis H 0 ″ “r k ˜ (x, y) (k = 1,..., K) is an average of 0 , According to the normal distribution of variance (1+ (a 1 to ) 2 ) σ 2 ”, the risk factor p (x, y) or weight w (p (x, y)) when specifically rejecting is calculated. Then, the weighted average v ′ k (x) is calculated by the following equation (22), and the time density curve {<t k , v ′ k (x) −v ′ 1 (x)> (k = 1) at the pixel x. , 2,..., K)}.
Figure 2008100121

こうして得られた時間濃度曲線は、画素xにおける真の時間濃度曲線{<t,d>}の一次変換である{<t,A(d−d)>}(だたしAは未知の係数)を近似している計測値であり、しかも、重み付き平均の効果によって、ランダムなノイズが抑制されている。また、他の画素yの画素値ベクトルに対しては、式から明らかなように、パーシャルボリューム効果の影響を補正したものが用いられている。さらに、本実施形態はコヒーレントフィルタの共通の特徴である、「時間平均を全く使用せず、また空間平均を画素xとの適合度に基づく重みを使って計算する」という性質を有する。従って、本実施形態によって、時間分解能を損なわず、パーシャルボリューム効果の影響を抑制し、しかもランダムなノイズが抑制された時間濃度曲線を得ることができる。なお、このようにして時間濃度曲線を求める方式を、特に「コヒーレントレグレッション法」と称す。 The time density curve thus obtained is {<t k , A (d k −d 1 )>} (which is a linear transformation of the true time density curve {<t k , d k >} at the pixel x. A is a measurement value approximating an unknown coefficient), and random noise is suppressed by the effect of the weighted average. In addition, as is apparent from the equation, a pixel value vector of another pixel y corrected for the influence of the partial volume effect is used. Furthermore, the present embodiment has a common feature of the coherent filter, that is, “the time average is not used at all, and the spatial average is calculated using a weight based on the degree of fitness with the pixel x”. Therefore, according to the present embodiment, it is possible to obtain a time density curve in which the influence of the partial volume effect is suppressed and random noise is suppressed without impairing the time resolution. The method for obtaining the time density curve in this way is particularly referred to as “coherent regression method”.

次に、具体的に、医療用のX線CTにおけるダイナミックCT撮影等で得られたダイナミックCT画像における、時間濃度曲線の臨床的利用の一例を説明する。この応用例では、造影剤を血管に急速に注入しながら、ダイナミックCT等の撮影を行い、人体組織中に存在する動脈の像の濃度変化を時間濃度曲線として計測することによって、当該組織における血流動態を診断しようとするものである。   Next, an example of clinical use of a time density curve in a dynamic CT image obtained by dynamic CT imaging or the like in medical X-ray CT will be specifically described. In this application example, imaging such as dynamic CT is performed while rapidly injecting a contrast medium into a blood vessel, and a change in the density of an image of an artery present in a human tissue is measured as a time density curve, whereby blood in the tissue is measured. It is intended to diagnose the flow dynamics.

この応用例において、多くの場合、人体組織中の動脈は一般に非常に細いために、CTによる断層画像上に現れる動脈の像は、パーシャルボリューム効果を生じる。さらに、像にはランダムなノイズが含まれていることは言うまでもない。このため、従来の方法では、動脈に関する十分に正確な時間濃度曲線を得ることは困難であり、強いて計測を行えば、動脈に関する真の時間濃度曲線<t,D>の一次変換である<t,A(D−D)>(ここにDは動脈の像に相当する一群の画素の、時刻tにおける(スカラー値である)画素値を表す。また、k=1,2,・・・,K)をある程度近似する測定値<t,(v(x)−v(x))>しか得られなかった。この測定値はランダムなノイズを含む。また、パーシャルボリューム効果の影響のために、係数Aは未知のままである。 In this application example, in many cases, arteries in human tissue are generally very thin, and thus an arterial image that appears on a tomographic image by CT produces a partial volume effect. Furthermore, it goes without saying that the image contains random noise. For this reason, in the conventional method, it is difficult to obtain a sufficiently accurate time concentration curve related to the artery, and if measurement is performed forcibly, it is a linear transformation of the true time concentration curve <t k , D k > related to the artery. <T k , A (D k −D 1 )> (where D k represents a pixel value (which is a scalar value) of a group of pixels corresponding to an artery image at time t k , and k = 1. , 2,..., K) to some extent, only measured values <t k , (v k (x) −v 1 (x))> were obtained. This measurement includes random noise. Also, the coefficient A remains unknown due to the effect of the partial volume effect.

そこで、<t,A(D−D)>を十分に近似する測定値<t,(v´(x)−v´(x))>(k=1,2,・・・,K)を得ることができる。一方、同じ断層画像上で観察できる静脈の中には、相当に太いものが存在し、従ってそれらの静脈に関しては、従来の方法で、時間濃度曲線の十分に良い近似値<t,(J−J)>(k=1,2,・・・,K)を得ることができる。ここにJは静脈の像に相当する一群の画素の、時刻tにおける画素値を表す。 Therefore, measured values <t k , (v ′ k (x) −v ′ 1 (x))> (k = 1, 2,...) Sufficiently approximating <t k , A (D k −D 1 )>. .., K) can be obtained. On the other hand, some veins that can be observed on the same tomographic image are considerably thick. Therefore, for those veins, a sufficiently good approximation value <t k , (J k− J 1 )> (k = 1, 2,..., K) can be obtained. Here, J k represents a pixel value at a time t k of a group of pixels corresponding to a vein image.

ところで、血液循環に関する時間濃度曲線においては、命題S:「もし、時刻tにおける血中の造影剤濃度が0であるならば、どの血管dに関する時間濃度曲線<t,(d-d)>も、その曲線下面積(AUC:Area Under CurVe)が一致する」という性質が成り立つことが知られている。ここで言う曲線下面積とは、時間濃度曲線<t,(d-d)>の時間tに関する積分を意味する。 By the way, in the time concentration curve relating to blood circulation, the proposition S: “If the contrast agent concentration in the blood at time t 1 is 0, the time concentration curve relating to which blood vessel d <t k , (d k -d 1 )> is also known to have the property that the area under the curve (AUC) matches. The area under the curve here means the integration of the time density curve <t k , (d k -d 1 )> with respect to time t.

従って、ある血管dに関する時間濃度曲線<t,(d-d)>の曲線下面積AUC(d)は、例えば次式によって近似的に計算することができる。

Figure 2008100121
Therefore, the area under the curve AUC (d) of the time concentration curve <t k , (d k −d 1 )> relating to a certain blood vessel d can be approximately calculated by the following equation, for example.
Figure 2008100121

従って、静脈に関して従来の方法で得られた時間濃度曲線{<t,(J-J)>}に関する曲線下面積AUC(J)を(22)式を用いて計算することができる。(dにJを代入すればよい。)また、動脈に関して、仮に、時間濃度曲線{<t,(D−D)>}が知られていれば、曲線下面積AUC(D)を(18)式を用いて同様に計算することができ、しかも上記命題Sに従って
AUC(D)≒AUC(J) … (24)
が成り立つ筈である。しかし、実際には、時間濃度曲線<t,(D−D)>は未知であるため、AUC(D)は計算できない。
Therefore, the area under the curve AUC (J) for the time concentration curve {<t k , (J k −J 1 )>} obtained by the conventional method for veins can be calculated using the equation (22). (J can be substituted for d.) Also, regarding the artery, if the time concentration curve {<t k , (D k −D 1 )>} is known, the area under the curve AUC (D) is It can be calculated in the same way using equation (18), and AUC (D) ≈AUC (J) according to the proposition S above (24)
Should be true. However, in practice, the time concentration curve <t k , (D k −D 1 )> is unknown, so AUC (D) cannot be calculated.

一方、本発明に係る方式で得られた時間濃度曲線<t,(v´(x)−v´(x))>は、<t,A(D−D)>を近似するものであり、後者は未知の係数Aを含んでいる。このため、{<t,(v´(x)−v´(x))>}から(23)式を用いて具体的に計算できる曲線下面積AUC(v´)は、AUC(D)のちょうどA倍でなくてはならない。すなわち、
AUC(v´)≒A AUC(D) … (25)
である。すなわち、(24)式と(25)式から、
A≒AUC(v´)/AUC(J) … (26)
という関係が成り立つ。(26)式の右辺は(23)式を用いて具体的に計算できるため、未知であった係数Aの値が具体的に決定できる。そこで、この係数Aの値を用いて時間濃度曲線<t,(v´(x)−v´(x))/A>を構成すれば、これは、動脈の時間濃度曲線<t,(D−D)>を近似するものに他ならない。このように、曲線下面積を用いて、未知であった比例係数Aの値を決定した時間濃度曲線を構成する方法を「AUC法」と呼ぶ。
On the other hand, the time concentration curve <t k , (v ′ k (x) −v ′ 1 (x))> obtained by the method according to the present invention represents <t k , A (D k −D 1 )>. The latter contains an unknown coefficient A. For this reason, the area under the curve AUC (v ′) that can be specifically calculated using the equation (23) from {<t k , (v ′ k (x) −v ′ 1 (x))>} is AUC ( It must be exactly A times D). That is,
AUC (v ′) ≈A AUC (D) (25)
It is. That is, from the equations (24) and (25),
A≈AUC (v ′) / AUC (J) (26)
This relationship holds. Since the right side of equation (26) can be specifically calculated using equation (23), the unknown value of coefficient A can be specifically determined. Therefore, if the value of the coefficient A is used to construct a time concentration curve <t k , (v ′ k (x) −v ′ 1 (x)) / A>, this is equivalent to the arterial time concentration curve <t It is nothing but an approximation of k , (D k −D 1 )>. The method of constructing the time concentration curve in which the unknown value of the proportionality coefficient A is determined using the area under the curve in this way is referred to as “AUC method”.

以上から、ダイナミックCT撮影等で得られたダイナミックCT画像における、時間濃度曲線の臨床的利用において、上記コヒーレントレグレッション法に、さらに上記AUC法を組み合わせることによって、従来の方法では計測が困難あるいは不可能であった、細い動脈の時間濃度曲線に関しても、パーシャルボリューム効果およびランダムなノイズの影響を排除し、しかも、未知の比例係数Aを含まない測定値が得られる。   From the above, in the clinical use of the time density curve in dynamic CT images obtained by dynamic CT imaging or the like, it is difficult or impossible to measure with the conventional method by combining the AUC method with the coherent regression method. Even for a fine arterial time-concentration curve, which is possible, a measurement value that eliminates the influence of the partial volume effect and random noise and does not include the unknown proportionality factor A is obtained.

なお、もちろんAUC法は、単独で従来の方法で計測された動脈に関する時間濃度曲線<t,(v´(x)−v´(x))>に対しても適用でき、(ランダムなノイズやパーシャルボリューム効果の影響は排除できないものの、)未知であった比例係数Aの値を決定した時間濃度曲線を構成できる。 Of course, the AUC method can also be applied to a time-concentration curve <t k , (v ′ k (x) −v ′ 1 (x))> relating to an artery measured by a conventional method alone (randomly). Although the influence of noise and the partial volume effect cannot be excluded, it is possible to construct a time density curve in which the value of the proportionality coefficient A that has been unknown is determined.

(上矢状静脈洞の時間濃度曲線を使った脳動脈の時間濃度曲線の補正(パーシャルボリューム効果の影響を抑圧))
パーシャルボリュームの影響を抑圧するために、このコヒーレントレグレッションに代えて又は併用して、上矢状静脈洞の時間濃度曲線Csss(t)を使って脳動脈の時間濃度曲線Ca(t)を補正するようにしてもよい。
(Correction of cerebral artery time concentration curve using upper sagittal sinus time concentration curve (suppression of partial volume effect))
In order to suppress the influence of partial volume, the time concentration curve Ca (t) of the cerebral artery is corrected using the time concentration curve Csss (t) of the superior sagittal sinus instead of or in combination with this coherent regression. You may make it do.

まず、図5のステップS7において、図11に示すように、CT画像上で上矢状静脈洞を取り囲むように大き目に上矢状静脈洞ROIを設定する。上矢状静脈洞は動脈に比べて大きく、また位置も固定的であるので、上矢状静脈洞ROI設定は容易である。この大き目の上矢状静脈洞ROIには、複数の画素が含まれている。   First, in step S7 of FIG. 5, as shown in FIG. 11, a large upper sagittal sinus ROI is set so as to surround the upper sagittal sinus on the CT image. Since the upper sagittal sinus is larger than the artery and the position is fixed, setting the upper sagittal sinus ROI is easy. The large upper sagittal sinus ROI includes a plurality of pixels.

次に、上矢状静脈洞ROIの中の全ての画素がその全域に渡って上矢状静脈洞に含まれるように、上矢状静脈洞ROIを縮小処理する(S8)。縮小処理としては、例えば、まず、上矢状静脈洞ROI内の画素各々について、しきい値処理(二値化)を実行し、ROI内の二値マップ(“0”,“1”)を作成する。当該しきい値は、上矢状静脈洞の像を、その周辺の組織や骨の像と分離する値に設定される。“1”は上矢状静脈洞の像上の画素であることを表し、“0”は周辺の組織や骨の像上の画素であることを表している。この二値マップの各画素(中心画素)を、その近傍4または8画素の値に従って、置換する。中心画素が“1”であって、且つ、近傍4または8画素全てが“1”である場合のみ、中心画素の値を“1”に維持する。つまり、中心画素が“0”である場合はもちろん、たとえ“1”であったとしても、近傍4または8画素のなかの1つでも“0”を示しているときには、当該中心画素の値を“0”に置換する。従って、上矢状静脈洞ROIは、上矢状静脈洞の像の外形よりも、少なくとも1画素分縮小される。それにより縮小処理を受けた上矢状静脈洞ROIの中の全ての画素は、上矢状静脈洞像上の画素であるという条件を高い確度で実現され得る。   Next, the upper sagittal sinus ROI is reduced so that all the pixels in the upper sagittal sinus ROI are included in the upper sagittal sinus over the entire area (S8). As the reduction processing, for example, first, threshold processing (binarization) is performed for each pixel in the upper sagittal sinus ROI, and a binary map (“0”, “1”) in the ROI is obtained. create. The threshold value is set to a value that separates the image of the superior sagittal sinus from the images of surrounding tissues and bones. “1” represents a pixel on the image of the upper sagittal sinus, and “0” represents a pixel on an image of surrounding tissue or bone. Each pixel (center pixel) of this binary map is replaced according to the value of its neighboring 4 or 8 pixels. Only when the central pixel is “1” and all the neighboring 4 or 8 pixels are “1”, the value of the central pixel is maintained at “1”. That is, when the central pixel is “0”, even if it is “1”, if any one of the neighboring 4 or 8 pixels indicates “0”, the value of the central pixel is changed. Replace with “0”. Accordingly, the upper sagittal sinus ROI is reduced by at least one pixel from the outline of the image of the upper sagittal sinus. As a result, the condition that all the pixels in the upper sagittal sinus ROI subjected to the reduction process are pixels on the upper sagittal sinus sinus image can be realized with high accuracy.

また、この手法に代えて、時間濃度曲線の曲線下面積AUCを使って上矢状静脈洞ROIを修正するようにしてもよい。この場合、大き目のROIを探索範囲として、その中の画素各々について時間濃度曲線の曲線下面積AUCを計算する。造影効果により上矢状静脈洞像上の画素の曲線下面積AUCは、周辺画素のそれに比べて明らかに高位を示す。従って、この曲線下面積AUCに対してしきい値処理を実行することにより、ROIの中から上矢状静脈洞像上の画素だけを選別することができる。   Instead of this method, the upper sagittal sinus ROI may be corrected using the area AUC under the time density curve. In this case, the area under the curve AUC of the time density curve is calculated for each of the pixels in the large ROI as a search range. The area under the curve AUC of the pixel on the upper sagittal sinus image is clearly higher than that of the peripheral pixel due to the contrast effect. Therefore, by executing threshold processing on the area AUC under the curve, only pixels on the upper sagittal sinus image can be selected from the ROI.

このようにしていずれかの手法又は両手法を併用してアンド条件によりピックアップされた上矢状静脈洞像上であるという確度の高い複数の画素に対して、各画素の時間濃度曲線が平均化され、上矢状静脈洞の時間濃度曲線Csss(t)が作成される(S9)。   In this way, the time density curve of each pixel is averaged for a plurality of highly accurate pixels that are on the upper sagittal sinus image picked up by the AND condition using either method or both methods in combination. Then, a time concentration curve Csss (t) of the superior sagittal sinus is created (S9).

ここで、ヨード造影剤は関門非通過であるので、原理的に、ヨード濃度は脳動脈と脳静脈とで変化しない、つまり、上矢状静脈洞の時間濃度曲線Csss(t)の曲線下面積AUCは、S6で作成した脳動脈の時間濃度曲線Ca(t)の曲線下面積AUCに略等価になる。従って、図10に示すように、上矢状静脈洞の時間濃度曲線Csss(t)の曲線下面積AUCsssに対して、S6で作成した脳動脈の時間濃度曲線Ca(t)の曲線下面積AUCaが略等価になるように、S6で作成した脳動脈の時間濃度曲線Ca(t)の各時刻値に、AUC(sss/AUCa)を乗算することで補正する(S10)。   Here, since the iodine contrast medium does not pass through the barrier, in principle, the iodine concentration does not change between the cerebral artery and the cerebral vein, that is, the area under the curve of the time density curve Csss (t) of the upper sagittal sinus. AUC is substantially equivalent to the area under the curve AUC of the time concentration curve Ca (t) of the cerebral artery created in S6. Therefore, as shown in FIG. 10, the area under the curve AUCa of the time concentration curve Ca (t) of the cerebral artery created in S6 with respect to the area under the curve AUCss of the time concentration curve Csss (t) of the upper sagittal sinus. Are corrected by multiplying each time value of the time concentration curve Ca (t) of the cerebral artery created in S6 by AUC (sss / AUCa) (S10).

次に、以上のようにノイズ及びパーシャルボリューム効果が抑圧された図13(a)に示す脳動脈の時間濃度曲線Ca(t)を使って、脳組織(毛細血管)の血流動態の様子を定量化する。そのためにまず、脳組織上の各画素について、図13(b)に示す時間濃度曲線Ci(t)が作成される(S11)。   Next, using the time concentration curve Ca (t) of the cerebral artery shown in FIG. 13A in which the noise and the partial volume effect are suppressed as described above, the state of the blood flow dynamics of the brain tissue (capillary blood vessel) is shown. Quantify. For this purpose, first, a time density curve Ci (t) shown in FIG. 13B is created for each pixel on the brain tissue (S11).

次に、S12に示すように、左右エリアで別々の脳動脈時間濃度曲線Ca(t)を使って、画素ごとに、脳動脈時間濃度曲線Ca(t)を入力関数、脳組織の時間濃度曲線Ci(t)を出力関数として、トレーサーが毛細血管を通過する過程の特徴を、伝達関数MTFとして求める。つまり、左エリアの脳組織の時間濃度曲線Ci(t)に対しては同じ左エリアの脳動脈の時間濃度曲線Ca(t)を使い、また右エリアの脳組織の時間濃度曲線Ci(t)に対しては同じ右エリアの脳動脈の時間濃度曲線Ca(t)を使って、伝達関数MTFを求める。また、上述したように脳動脈時間濃度曲線Ca(t)は、ACA,MCA,PCAごとに作成されるので、各Ca(t)ごとに伝達関数MTFの計算が繰り返される。   Next, as shown in S12, the cerebral artery time density curve Ca (t) is used for each pixel by using separate cerebral artery time density curves Ca (t) in the left and right areas, and the brain tissue time density curve is used as an input function. Using Ci (t) as an output function, the characteristic of the process of the tracer passing through the capillary is obtained as the transfer function MTF. That is, for the time density curve Ci (t) of the brain tissue in the left area, the time density curve Ca (t) of the same cerebral artery in the left area is used, and the time density curve Ci (t) of the brain tissue in the right area. , The transfer function MTF is obtained using the time concentration curve Ca (t) of the cerebral artery in the same right area. Moreover, since the cerebral artery time concentration curve Ca (t) is created for each ACA, MCA, and PCA as described above, the calculation of the transfer function MTF is repeated for each Ca (t).

ここでは、図14に示すように、脳動脈の時間濃度曲線Ca(t)と、毛細血管の時間濃度曲線Ci(t)との間で成り立つ理想的な関係を解析モデルとして用いる。つまり脳組織に入る直前の血管から造影剤(トレーサー)を瞬間的に注入した場合、毛細血管を含む脳組織単位体積(1画素)内の時間濃度曲線は立ち上がりが垂直で、若干の勾配を持って立ち下がる形になり、これが組織のインパルス応答、すなわち伝達関数を表している。そこで伝達関数を、矩形関数で近似する。これをbox−MTF法と呼ぶ。   Here, as shown in FIG. 14, an ideal relationship established between the time concentration curve Ca (t) of the cerebral artery and the time concentration curve Ci (t) of the capillary is used as the analysis model. In other words, when a contrast medium (tracer) is injected instantaneously from a blood vessel just before entering the brain tissue, the time density curve in the brain tissue unit volume (1 pixel) including the capillary blood vessels has a slight rise in the vertical rise. This represents the impulse response of the tissue, that is, the transfer function. Therefore, the transfer function is approximated by a rectangular function. This is called a box-MTF method.

図15にbox−MTF法の原理を示している。脳動脈の時間濃度曲線Ca(t)に対して、矩形関数で表される伝達関数box−MTFをコンボリュートし、それにより得られたCi´(t)と、S11で作成した実測Ci(t)との残差を評価し、その残差の二乗和を減少させるように伝達関数box−MTFを修正する。このルーチンを繰り返すことで、残差を最小化させる。   FIG. 15 shows the principle of the box-MTF method. The transfer function box-MTF represented by a rectangular function is convoluted with the time concentration curve Ca (t) of the cerebral artery, and Ci ′ (t) obtained by the convolution and the measured Ci (t) created in S11. ) And the transfer function box-MTF is corrected so as to reduce the sum of squares of the residual. By repeating this routine, the residual is minimized.

残差を最小化する伝達関数box−MTFに基づいて、図14に示すように、CBP、CBV、MTTが計算され(S13)、またS12で最小化した残差の二乗和がそのままErrとして取り込まれる。厳密には、
CBP=CBP
CBV=(1−Ht)/(1−b*Ht)*CBV'
MTT=(1−Ht)/(1−b*Ht)*MTT'
で補正される。ここで、Htは大血管のヘマトクリット値であり、b*Htは末梢血管のヘマトクリット値(一般的にはbは0.7程度)である。
Based on the transfer function box-MTF that minimizes the residual, as shown in FIG. 14, CBP, CBV, and MTT are calculated (S13), and the square sum of the residual minimized in S12 is directly taken as Err. It is. Strictly speaking,
CBP = CBP
CBV = (1-Ht) / (1-b * Ht) * CBV ′
MTT = (1-Ht) / (1-b * Ht) * MTT ′
It is corrected by. Here, Ht is the hematocrit value of the large blood vessel, and b * Ht is the hematocrit value of the peripheral blood vessel (generally, b is about 0.7).

こうしてbox−MTF法により得られた伝達関数から計算されたCBP、CBV、MTT、Errに対して、図16に示すように、個々に出力範囲(適正範囲)が設定される(S14)。それぞれ対応する出力範囲内の値を持つ画素については、その値が維持され、出力範囲外の値を持つ画素にはその値から、例えば表示上で黒に対応する値に置換される(S15)。   As shown in FIG. 16, output ranges (appropriate ranges) are individually set for CBP, CBV, MTT, and Err calculated from the transfer function thus obtained by the box-MTF method (S14). For pixels having a value within the corresponding output range, the value is maintained, and for pixels having a value outside the output range, the value is replaced with a value corresponding to black on the display, for example (S15). .

次に、図17に示すように、出力適正化を受けたCBP、CBV、MTT、Errからそれぞれのマップが生成される(S16)。CBP、CBV、MTT、Errのインデックスは、前大脳動脈ACA,中大脳動脈MCA,後大脳動脈PCAごと、しかも各スライスに対して個々に計算される。従って、マップは、1スライスでも、図18に示すように、4×3=12枚になる。設定する脳動脈数を増加すれば、その増加数の4倍でマップは増える。このような多くのマップを総合的に評価することは現実的ではない。そこで、マップ枚数を減少させるために、マップを合成する(S17)。   Next, as shown in FIG. 17, respective maps are generated from CBP, CBV, MTT, and Err that have undergone output optimization (S16). The CBP, CBV, MTT, and Err indices are calculated individually for each anterior cerebral artery ACA, middle cerebral artery MCA, posterior cerebral artery PCA, and for each slice. Therefore, the map is 4 × 3 = 12 even in one slice as shown in FIG. If the number of cerebral arteries to be set is increased, the number of maps increases by four times the increase. It is not realistic to evaluate many such maps comprehensively. Therefore, to reduce the number of maps, the maps are synthesized (S17).

合成方法としては、前大脳動脈ACAのCBPマップと、中大脳動脈MCAのCBPマップと、後大脳動脈PCAのCBPマップとを、前大脳動脈ACAの残差Err、中大脳動脈MCAの残差Err、後大脳動脈PCAの残差Errに基づいて合成する。例えば前大脳動脈ACAの時間濃度曲線Ca(t)とその支配下にある脳組織の時間濃度曲線Ci(t)とから伝達関数MTFを求めた場合、その残差Errは比較的少なく、逆に、支配下にない脳組織の時間濃度曲線Ci(t)から伝達関数MTFを求めた場合、その残差Errは比較的多くなる。つまり、残差Errは各脳動脈の支配可能性を表している。   As a synthesis method, a CBP map of the anterior cerebral artery ACA, a CBP map of the middle cerebral artery MCA, and a CBP map of the posterior cerebral artery PCA are obtained. Synthesize based on residual Err of posterior cerebral artery PCA. For example, when the transfer function MTF is obtained from the time concentration curve Ca (t) of the anterior cerebral artery ACA and the time concentration curve Ci (t) of the brain tissue under its control, the residual Err is relatively small, and conversely When the transfer function MTF is obtained from the time concentration curve Ci (t) of the brain tissue that is not under control, the residual Err is relatively large. That is, the residual Err represents the dominance possibility of each cerebral artery.

従って、各画素ごとに、前大脳動脈ACAのCBP値と、中大脳動脈MCAのCBP値と、後大脳動脈PCAのCBP値との中から、最も残差Errの低い値に対応するCBP値をその画素の値として選択する。こうして合成されたマップは、前大脳動脈ACA,中大脳動脈MCA,後大脳動脈PCAの支配下にある可能性の高い脳組織のCBP値から構成される。他のインデックスCBV、MTTのマップ合成についても同様である。   Therefore, for each pixel, a CBP value corresponding to the lowest residual Err is selected from the CBP value of the anterior cerebral artery ACA, the CBP value of the middle cerebral artery MCA, and the CBP value of the posterior cerebral artery PCA. Select as the value of that pixel. The map synthesized in this way is composed of CBP values of brain tissue that are likely to be under the control of the anterior cerebral artery ACA, middle cerebral artery MCA, and posterior cerebral artery PCA. The same applies to the map synthesis of other indexes CBV and MTT.

また、最も残差Errが低い値を示す脳動脈が、前大脳動脈ACAと中大脳動脈MCAと後大脳動脈PCAとのいずれであるかを選択し、選択された前大脳動脈ACA、中大脳動脈MCA又は後大脳動脈PCAに応じて相違する色情報を与えることにより、前大脳動脈ACAの支配領域と、中大脳動脈MCAの支配領域と、後大脳動脈PCAの支配領域とを色分けにより区別された支配マップを作成するようにしても良い。   The cerebral artery having the lowest residual Err is selected from the anterior cerebral artery ACA, the middle cerebral artery MCA, and the posterior cerebral artery PCA, and the selected anterior cerebral artery ACA and middle cerebral artery are selected. By giving different color information depending on the MCA or the posterior cerebral artery PCA, the dominant region of the anterior cerebral artery ACA, the dominant region of the middle cerebral artery MCA, and the dominant region of the posterior cerebral artery PCA were distinguished by color coding A dominance map may be created.

ここで、マップ合成について以下に詳細に説明する。動脈に対応する位置にある画素から得られる時間濃度曲線は動脈血中造影剤濃度を反映しており、これに上記したコヒーレントレグレッション法等を適用して正確な動脈血中造影剤濃度の時間濃度曲線を得ることができる。このような脳動脈の時間濃度曲線は動脈ごとに作成可能であり、それぞれ血行状態によって違いがある。特に脳血管障害等を起こしている症例に於いてはこの違いが著しい場合がある。例えばK箇所の動脈で得た脳動脈の時間濃度曲線をA(t) (k=1,2,・・・,K)で表すことにする。 Here, the map composition will be described in detail below. The time density curve obtained from the pixel located at the position corresponding to the artery reflects the arterial blood contrast medium concentration, and the above-mentioned coherent regression method is applied to the time density curve of the correct arterial blood contrast medium concentration. Can be obtained. Such a time-concentration curve of the cerebral artery can be created for each artery, and there are differences depending on the blood circulation state. In particular, this difference may be significant in cases that have cerebrovascular disorders. For example, the time concentration curve of the cerebral artery obtained from K arteries is represented by A k (t) (k = 1, 2,..., K).

ある組織の時間濃度曲線を、その組織を栄養している(支配している)動脈の脳動脈の時間濃度曲線と比較することによって、当該組織に於ける微小循環(毛細血管系の構造、機能)を反映するCBP等のインデックスを得ることが出来る。これらのインデックスは各部位(x,y,z)ごとに算出されるので、その値を画素値とする画像を構成することができ、このような画像が、インデックスマップである。例えばR種類(典型的には上述したようにCBP,CBV,MTT,Errの4種)のインデックスが得られる場合、R枚のマップが構成できる。このようにして作成されたR枚のマップは、各画素がベクトル値をとる1枚のマップ(ベクトル値マップ)と見なすことができる。すなわち、
(x,y,z) = <Pk,1(x,y,z), Pk,2(x,y,z), ... , Pk,R(x,y,z)>
となる。
By comparing the time-concentration curve of a tissue with the time-concentration curve of the cerebral artery of the artery that is feeding (dominating) the tissue, the microcirculation (structure and function of the capillary system in the tissue) ) Reflecting CBP or the like can be obtained. Since these indexes are calculated for each part (x, y, z), an image having the values as pixel values can be configured, and such an image is an index map. For example, when R types of indexes (typically four types of CBP, CBV, MTT, and Err as described above) are obtained, R maps can be configured. The R maps created in this way can be regarded as one map (vector value map) in which each pixel takes a vector value. That is,
V k (x, y, z) = <P k, 1 (x, y, z), P k, 2 (x, y, z),..., P k, R (x, y, z) >
It becomes.

例えばCBPスタディでは、典型的には上述したようにR=4とし、Pk,1(x,y)はCBPの値を、Pk,2(x,y,z)はCBVの値を、Pk,3(x,y,z)はMTTの値を、Pk,4(x,y,z)は残差Errの値を表すように構成できる。 For example, in a CBP study, typically R = 4 as described above, P k, 1 (x, y) is the CBP value, P k, 2 (x, y, z) is the CBV value, P k, 3 (x, y, z) can be configured to represent the value of MTT, and P k, 4 (x, y, z) can represent the value of residual Err.

部位(x,y,z)の内で、解析の対象となる臓器に対応していないことが明らかであるようなものは初めから解析の対象外とし、Pk,r(x,y)には解析対象外を示す特殊な値を代入するとよい(上記ステップS14,S15)。そのような値として、負で絶対値が大きい値を用いると便利である。或いは、ベクトルV(x,y,z)に追加されるべきさらにもう一個の要素として、
k,R+1(x,y,z)=((x,y,z)が解析対象外であれば0, さもなくば1)
というマップを作っても良い。このようなマップは「マスク」と呼ばれる。
Among the parts (x, y, z), those that are clearly not corresponding to the organ to be analyzed are excluded from the analysis from the beginning, and P k, r (x, y) May be substituted with a special value indicating that it is not subject to analysis (steps S14 and S15). As such a value, it is convenient to use a negative value having a large absolute value. Alternatively, as yet another element to be added to the vector V k (x, y, z)
P k, R + 1 (x, y, z) = (0 if (x, y, z) is not subject to analysis, 1 otherwise)
You can make a map. Such a map is called a “mask”.

このようなベクトル値マップVは、参照した脳動脈の時間濃度曲線Aごとに作成される。例えば、左右の中大脳動脈、前大脳動脈、後大脳動脈から脳動脈の時間濃度曲線を得たとするとK=6、さらにおよび病変部周辺にある動脈数カ所から脳動脈の時間濃度曲線を得たとすると、K=10〜15程度になる。 Such a vector value map V k is created for each referenced cerebral artery time density curve A k . For example, if the time concentration curve of the cerebral artery is obtained from the left and right middle cerebral artery, the anterior cerebral artery, and the posterior cerebral artery, then K = 6, and further, the time concentration curve of the cerebral artery is obtained from several arteries around the lesion. K = about 10-15.

このうち、右半球にある動脈から得られた脳動脈の時間濃度曲線は右半球に属する部位(x,y,z)の解析にだけ、また左半球にある動脈から得られた脳動脈の時間濃度曲線は左半球に属する部位(x,y,z)の解析にだけ、用いられるべきである。そこで、右半球と左半球の境目(正中線)を直線、曲線もしくは折れ線、或いは平面、曲面等として操作者が指定し、それぞれの半球ごとにマップを作るように構成するのが望ましい。しかしそれでもなお、片側の半球ごとにK=3〜10程度の数の脳動脈の時間濃度曲線が存在しうる。   Of these, the cerebral artery time-concentration curve obtained from the artery in the right hemisphere is only for the analysis of the part (x, y, z) belonging to the right hemisphere, and the time of the cerebral artery obtained from the artery in the left hemisphere. The concentration curve should be used only for the analysis of the part (x, y, z) belonging to the left hemisphere. Therefore, it is desirable that the operator designates the boundary (midline) between the right hemisphere and the left hemisphere as a straight line, a curved line or a polygonal line, a plane, a curved surface, or the like, and creates a map for each hemisphere. However, there may still be as many cerebral artery time concentration curves as K = 3-10 for each hemisphere on one side.

このように脳動脈の時間濃度曲線A (k=1,2,・・・,K)の数Kが大きい場合に、結果として得られるベクトル値マップV(k=1,2,・・・,K)の枚数が多いために、観察するのに不便である。すなわち通常のグレースケール画像あるいはカラースケール画像として観察しようとすれば、一つのマップがR枚の画像から構成され、これがK組あるのだから、合計K×R枚の画像を比較しなくてはならない。さらに、どの部位がどの動脈によって栄養されているのかは必ずしも自明でなく、解剖学的知識を用いて、各部位ごとにどのマップV(k=1,2,・・・,K)を観察するべきかを判断しなくてはならない。特に脳梗塞等の脳血管障害を生じている症例においては、組織を支配しているのがどの動脈かは、必ずしも解剖学的知識とは一致せず、異常な支配がしばしば見られる。これらの問題によって、ベクトル値マップの読影が難しいという問題点がある。 Thus, when the number K of the cerebral artery time density curve A k (k = 1, 2,..., K) is large, the resulting vector value map V k (k = 1, 2,... ., K) is inconvenient to observe because of the large number of sheets. That is, if a normal grayscale image or color scale image is to be observed, one map is composed of R images, and there are K sets, so a total of K × R images must be compared. . Further, it is not always obvious which part is nourished by which artery, and which map V k (k = 1, 2,..., K) is observed for each part using anatomical knowledge. You must decide what to do. In particular, in cases where cerebrovascular disorders such as cerebral infarction have occurred, the arteries that control the tissue do not necessarily match the anatomical knowledge, and abnormal control is often observed. Due to these problems, there is a problem that it is difficult to interpret the vector value map.

この問題を解決するために、マップ合成を行う。つまり、残差マップを利用して、K個のベクトル値マップV(k=1,2,・・・,K)を一つのベクトル値マップVに集約する。例えば、Pk,R(x,y,z)が残差マップである場合に、
V(x,y,z) = V(x,y,z) ただしk はk=1,2,・・・,Kのうちで|Pk,R(x,y,z)|が最小であるようなkとする。
To solve this problem, map composition is performed. That is, the K vector value maps V k (k = 1, 2,..., K) are aggregated into one vector value map V using the residual map. For example, if P k, R (x, y, z) is a residual map,
V (x, y, z) = V k (x, y, z) where k is the minimum of | P k, R (x, y, z) | Let k be such that

また、各部位においてk=1,2,・・・,Kのうちどれが採用されたかを示すためのマップ
(x,y,z)=(k=1,2,・・・,Kのうちで|Pk,R(x,y,z)|が最小であるようなk)
を追加することもできる。
In addition, a map for indicating which of k = 1, 2,...
P 0 (x, y, z) = (k = 1,2,..., K, where | P k, R (x, y, z) | is the smallest)
Can also be added.

この方式によれば、すなわち通常のグレースケール画像あるいはカラースケール画像として観察しようとするときR枚ないしR+1枚の画像を観察すればよい。   According to this method, that is, when an image is to be observed as a normal grayscale image or color scale image, R or R + 1 images may be observed.

この方式によれば、本来脳動脈の時間濃度曲線Aを使って算出されるべき部位(x,y,z)において誤ってAを使った算出結果が用いられる可能性がある。しかしながらこのような誤りが生じるには、V(x,y,z)の定義から明らかなように、|Pk,R(x,y,z)|<|Pj,R(x,y,z)|
となることが必要であり、このような関係は、AとAが極めて類似している場合にしか生じない。このため、部位(x,y,z)においてはV(x,y,z)とVj(x,y,z)は元々類似していると考えられ、この誤りによって、V(x,y,z)の解釈に誤りが生じる可能性はほとんどない。
According to this method, there is a possibility that a calculation result using Ak is erroneously used at a site (x, y, z) that should be calculated using the time density curve A j of the cerebral artery. However, in order to make such an error, as is apparent from the definition of V (x, y, z), | P k, R (x, y, z) | <| P j, R (x, y, z) |
This relationship only occurs when A j and A k are very similar. For this reason, V k (x, y, z) and V j (x, y, z) are considered to be similar to each other at the site (x, y, z). There is almost no possibility of errors in the interpretation of y, z).

実際にこの方法を適用すると、AとAが極めて類似している場合にだけ、概ね一様であると思われる組織内において、部位ごとにP(x,y,z)=kであったりP(x,y,z)=jであったりすることが起こり、そのときPk,r(x,y,z)≒Pj,r(x,y,z) (r=1,2,...,R)であって、どちらを採用しても結果はほとんど違いがないことが観察される。 When this method is actually applied, P 0 (x, y, z) = k for each site in the tissue that seems to be almost uniform only when A j and A k are very similar. Or P 0 (x, y, z) = j, and P k, r (x, y, z) ≈P j, r (x, y, z) (r = 1) , 2,..., R), and it is observed that the result is almost the same no matter which one is employed.

逆に、特定の動脈に支配されている組織であって、それに対応する脳動脈の時間濃度曲線Aが他のカーブと似ていない場合には、本方式を用いることによって、当該組織中の部位(x,y,z)においてはほぼ確実に、しかも自動的にV(x,y,z)が選択される。従って、上記P(x,y,z)を観察することによって、解剖学的知識なしに、どの組織がどの動脈の支配を受けているかを観察することができる。 Conversely, a tissue that are dominated by a specific artery, in a case where time-density curve A k of the cerebral artery corresponding to it is not similar to other curve, by using this method, of the tissue In the region (x, y, z), V k (x, y, z) is selected almost certainly and automatically. Therefore, by observing the above P 0 (x, y, z), it is possible to observe which tissue is controlled by which artery without anatomical knowledge.

ここで、図6に戻る。図20(a)、図20(b)、図20(c)、図20(d)に示すように、以上のように合成された、又は各脳動脈で単独のCBPマップ、CBVマップ、MTTマップ、Errマップに対して、複数画素を含む関心領域ROIを設定し(S18)、そのROI内の画素値(CBP値,CBV値、MTT値、Err値)の平均値(CBP平均値、CBV平均値、MTT平均値、Err平均値)を計算し(S19)、その平均値を診断材料とすることがある。この平均化に際して、上記ステップS14でCBP、CBV、MTT、Err各々に対して適正範囲を設定し、その範囲内の値を維持し、その範囲から外れた値は、例えば黒色表現に対応した最小値に置換したので、この置換した値を含めて平均化すると、その平均値には当然にして誤差が含まれてしまう。そのためこの平均化処理にあたっては、適正範囲内の値だけを選択して、または置換した値を除外して、平均化処理をすることが必要である。   Returning now to FIG. As shown in FIGS. 20A, 20B, 20C, and 20D, the CBP map, the CBV map, and the MTT synthesized as described above or each cerebral artery alone. A region of interest ROI including a plurality of pixels is set for the map and the Err map (S18), and an average value (CBP average value, CBV) of pixel values (CBP value, CBV value, MTT value, Err value) in the ROI is set. An average value, an MTT average value, and an Err average value) are calculated (S19), and the average value may be used as a diagnostic material. At the time of this averaging, an appropriate range is set for each of CBP, CBV, MTT, and Err in step S14, and a value within the range is maintained. A value outside the range is, for example, a minimum corresponding to black expression. Since the value is replaced with a value, if the average including the replaced value is averaged, the average value naturally includes an error. Therefore, in this averaging process, it is necessary to perform the averaging process by selecting only values within the appropriate range or excluding the replaced values.

また、この平均化のための関心領域ROIの設定にあたっては、CBPマップ、CBVマップ、MTTマップ、Errマップのいずれかのマップ上で当該関心領域ROIを設定すれば、そのROIが他のマップにも共通で用いられるようになっており、、ROI設定作業の簡素化を図り、また、同じROIに関する平均値(CBP平均値、CBV平均値、MTT平均値、Err平均値)の計算を可能としている。   Further, in setting the region of interest ROI for averaging, if the region of interest ROI is set on any of the CBP map, CBV map, MTT map, and Err map, the ROI becomes another map. Are also used in common, simplifying ROI setting work, and enabling calculation of average values (CBP average value, CBV average value, MTT average value, Err average value) related to the same ROI Yes.

以上のように本実施形態によれば、コヒーレントフィルタ又はコヒーレントレグレッションを併用することにより、空間及び時間分解能の低下を抑えて、ノイズを抑制し、それによりCBPスタディの解析精度を向上することができる。   As described above, according to the present embodiment, by using a coherent filter or coherent regression together, it is possible to suppress a reduction in spatial and temporal resolution and suppress noise, thereby improving analysis accuracy of a CBP study. it can.

なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。   Note that the present invention is not limited to the above-described embodiment as it is, and can be embodied by modifying the constituent elements without departing from the scope of the invention in the implementation stage. In addition, various inventions can be formed by appropriately combining a plurality of constituent elements disclosed in the embodiment. For example, some components may be deleted from all the components shown in the embodiment. Furthermore, constituent elements over different embodiments may be appropriately combined.

CBPスタディの原理説明図。Explanatory drawing of the principle of CBP study. 本発明の実施形態に係る脳組織内毛細血管の血流動態に関するインデックス演算装置の構成を示すブロック図。The block diagram which shows the structure of the index calculating apparatus regarding the blood flow dynamics of the capillary in the brain tissue which concerns on embodiment of this invention. 本実施形態のコヒーレントフィルタによる画像処理の説明図。Explanatory drawing of the image process by the coherent filter of this embodiment. 本実施形態におけるコヒーレントフィルタによるノイズ抑制処理の流れを示すフローチャート。The flowchart which shows the flow of the noise suppression process by the coherent filter in this embodiment. 本実施形態におけるインデックス演算処理全体の前半部のフローチャート。The flowchart of the first half part of the whole index calculation process in this embodiment. 本実施形態におけるインデックス演算処理全体の後半部のフローチャート。The flowchart of the second half part of the whole index calculation process in this embodiment. 図5のステップS3の分割線の一例を示す中間調画像。FIG. 6 is a halftone image showing an example of a dividing line in step S <b> 3 of FIG. 5. 図5のステップS4のAT,PT,TTを示す図。FIG. 6 is a diagram showing AT, PT, and TT in step S4 of FIG. 5. 図5のステップS4のATマップ,PTマップ,TTマップを示す図。FIG. 6 is a diagram showing an AT map, a PT map, and a TT map in step S4 of FIG. 図5のステップS10の脳動脈の時間濃度曲線の補正に関する補足図。FIG. 6 is a supplementary diagram regarding correction of the cerebral artery time density curve in step S10 of FIG. 5; 図5のステップS7で設定される上矢状静脈洞ROIを示す中間調画像。FIG. 6 is a halftone image showing the upper sagittal sinus ROI set in step S <b> 7 of FIG. 5. 図5のステップS6でスライス間で共通される脳動脈ROIを示す中間調画像。FIG. 6 is a halftone image showing a cerebral artery ROI shared between slices in step S <b> 6 of FIG. 5. 図5のステップS10、S11で作成された脳動脈の時間濃度曲線Ca(t)と脳組織の時間濃度曲線Ci(t)との一例を示す図。The figure which shows an example of the time density | concentration curve Ca (t) of the cerebral artery created in step S10, S11 of FIG. 5, and the time density | concentration curve Ci (t) of brain tissue. 図6のステップS12のbox−MTF法の原理説明図。FIG. 7 is a diagram illustrating the principle of the box-MTF method in step S12 of FIG. 図6のステップS12のbox−MTF処理の説明図。Explanatory drawing of the box-MTF process of step S12 of FIG. 図6のステップS14の各インデックスの出力範囲設定画面の一例を示す図。The figure which shows an example of the output range setting screen of each index of step S14 of FIG. 図6のステップS16で作成されたCBP,CBV,MTT,Errの各マップの一例を示す中間調画像。FIG. 7 is a halftone image showing an example of each map of CBP, CBV, MTT, and Err created in step S <b> 16 of FIG. 6. 図6のステップS16で脳動脈ごとに作成されたCBPマップ,CBVマップ,MTTマップ,Errマップを一覧で表示した中間調画像。A halftone image displaying a list of CBP maps, CBV maps, MTT maps, and Err maps created for each cerebral artery in step S16 of FIG. 図6のステップS17のマップ合成法を説明するための中間調画像。FIG. 7 is a halftone image for explaining the map composition method in step S <b> 17 of FIG. 6. 図6のステップS19で計算された平均値の表示例を示す中間調画像。The halftone image which shows the example of a display of the average value calculated by step S19 of FIG.

符号の説明Explanation of symbols

10…ガントリ部、
20…コンピュータ装置、
101…X線管、
101a…高電圧発生装置、
102…X線検出器、
103…データ収集部、
30…画像処理装置、
107…画像表示部、
109…入力部、
108…制御部、
104…前処理部、
105…メモリ部、
106…画像再構成部、
110…コヒーレントフィルタ処理部、
120…CBPスタディ処理部、
111…分散値推定部、
112…重み関数演算部、
113…画素値演算部(コヒーレントフィルタ部)、
121…CBPスタディ処理部、
122…ROI設定支援部、
122…時間濃度曲線作成部、
123…脳動脈時間濃度曲線補正部、
124…MTF処理部、
125…インデックス計算部、
126…マップ作成部、
127…マップ合成部。
10 ... Gantry part,
20 ... Computer device,
101 ... X-ray tube,
101a ... high voltage generator,
102 ... X-ray detector,
103 ... data collection unit,
30. Image processing device,
107: Image display unit,
109 ... input section,
108 ... control unit,
104 ... Pre-processing unit,
105: Memory unit,
106: Image reconstruction unit,
110: Coherent filter processing unit,
120 ... CBP study processing unit,
111... Variance value estimation unit,
112 ... Weight function calculation unit,
113 ... Pixel value calculation unit (coherent filter unit),
121... CBP study processing unit,
122 ... ROI setting support part,
122... Time concentration curve creation unit,
123 ... cerebral artery time concentration curve correction unit,
124 ... MTF processing unit,
125: Index calculation unit,
126 ... map creation part,
127: Map composition unit.

Claims (6)

脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のCT画像各々を構成する画素間の類似度を、前記複数のCT画像にわたる各画素の時間濃度曲線に基づいて判定し、
前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均し、
前記加重平均された複数のCT画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成し、
前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数を計算し、
前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算することを特徴とする脳組織内毛細血管の血流動態に関するインデックス演算方法。
Similarity between pixels constituting each of a plurality of CT images obtained continuously by imaging the brain of a subject injected with a contrast agent having no cerebral vascular permeability is determined for each pixel across the plurality of CT images. Judgment based on time concentration curve,
A weighted average of the pixels in each region with a weight according to the similarity,
Generating a time density curve of cerebral artery pixels and a time density curve of brain tissue pixels from the plurality of weighted averaged CT images;
Calculating the transfer function of the time density curve of the brain tissue pixel to the time density curve of the cerebral artery pixel;
An index calculation method for blood flow dynamics of capillaries in brain tissue, wherein a plurality of types of indexes related to blood flow dynamics of capillaries in the brain tissue are calculated based on the transfer function.
時間軸に対する組織の濃度変化のインパルス応答が矩形関数で表されると仮定する解析モデルを用い、前記脳動脈画素の時間濃度曲線と前記伝達関数のコンボリューションにより得られる時間濃度曲線と、前記脳組織画素の時間濃度曲線との間の残差の指標が最小化するようにカーブフィッティングにより前記伝達関数を決定し、得られた伝達関数から前記インデックスを算出することを特徴とする請求項1記載のインデックス演算方法。 Using an analysis model that assumes that the impulse response of the tissue concentration change with respect to the time axis is represented by a rectangular function, the time concentration curve obtained by convolution of the cerebral artery pixel time concentration and the transfer function, and the brain The index is calculated from the obtained transfer function by determining the transfer function by curve fitting so that an index of the residual between the tissue pixel and the time density curve is minimized. Index calculation method. 前記インデックスとして、前記得られた伝達関数から、
CBP:前記脳組織内の単位体積及び単位時間あたりの血流量
CBV:脳組織内の単位体積あたりの血液量
MTT:毛細血管の血液平均通過時間
Err:前記最小化された残差の二乗和
とが計算されることを特徴とする請求項2記載のインデックス演算方法。
As the index, from the obtained transfer function,
CBP: blood volume per unit volume and unit time in the brain tissue
CBV: blood volume per unit volume in brain tissue
MTT: Average blood passage time in capillaries
Err: sum of squares of the minimized residuals
The index calculation method according to claim 2, wherein and are calculated.
前記類似度が比較的高い場合には比較的高い重み係数が与えられ、前記類似度が比較的低い場合には比較的低い重み係数が与えられることを特徴とする請求項1記載のインデックス演算方法。 2. The index calculation method according to claim 1, wherein when the similarity is relatively high, a relatively high weighting factor is given, and when the similarity is relatively low, a relatively low weighting factor is given. . 脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のCT画像各々を構成する画素間の類似度を、前記複数のCT画像にわたる各画素の時間濃度曲線に基づいて判定する手段と、
前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均する手段と、
前記加重平均された複数のCT画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成する手段と、
前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数をカーブフィッティングにより計算する手段と、
前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算する手段とを具備することを特徴とする脳組織内毛細血管の血流動態に関するインデックス演算装置。
Similarity between pixels constituting each of a plurality of CT images obtained continuously by imaging the brain of a subject injected with a contrast agent having no cerebral vascular permeability is determined for each pixel across the plurality of CT images. Means for determining based on the time concentration curve;
Means for performing a weighted average of pixels for each image with a weight according to the similarity,
Means for generating a time density curve of cerebral artery pixels and a time density curve of brain tissue pixels from the plurality of weighted averaged CT images;
Means for calculating a transfer function of the time density curve of the brain tissue pixel with respect to the time density curve of the cerebral artery pixel by curve fitting;
An index calculating apparatus for blood flow dynamics of capillaries in brain tissue, comprising means for calculating a plurality of types of indexes related to blood flow dynamics of capillaries in the brain tissue based on the transfer function.
脳組織内毛細血管の血流動態に関するインデックスをコンピュータに演算させるためのプログラムコードを記録したコンピュータ読み取り可能な記憶媒体であって、前記プログラムコードは、
脳血管透過性を持たない造影剤を注入された被検体の脳を撮影対象として連続的に取得した複数のCT画像各々を構成する画素間の類似度を、前記複数のCT画像にわたる各画素の時間濃度曲線に基づいて判定する手段と、
前記類似度に応じた重みにより画像ごとに画素を局所内で加重平均する手段と、
前記加重平均された複数のCT画像から脳動脈画素の時間濃度曲線と脳組織画素の時間濃度曲線とを生成する手段と、
前記脳動脈画素の時間濃度曲線に対する脳組織画素の時間濃度曲線の伝達関数をカーブフィッティングにより計算する手段と、
前記伝達関数に基づいて、前記脳組織内の毛細血管の血流動態に関する複数種類のインデックスを演算する手段とを具備することを特徴とする記憶媒体。
A computer-readable storage medium storing a program code for causing a computer to calculate an index relating to blood flow dynamics of capillaries in brain tissue, the program code comprising:
Similarity between pixels constituting each of a plurality of CT images obtained continuously by imaging the brain of a subject injected with a contrast agent having no cerebral vascular permeability is determined for each pixel across the plurality of CT images. Means for determining based on the time concentration curve;
Means for performing a weighted average of pixels for each image with a weight according to the similarity,
Means for generating a time density curve of cerebral artery pixels and a time density curve of brain tissue pixels from the plurality of weighted averaged CT images;
Means for calculating a transfer function of the time density curve of the brain tissue pixel with respect to the time density curve of the cerebral artery pixel by curve fitting;
A storage medium comprising: means for calculating a plurality of types of indexes related to blood flow dynamics of capillaries in the brain tissue based on the transfer function.
JP2008010951A 2008-01-21 2008-01-21 Index calculation method, apparatus and storage medium for blood flow dynamics of capillaries in brain tissue Expired - Lifetime JP4714228B2 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2008010951A JP4714228B2 (en) 2008-01-21 2008-01-21 Index calculation method, apparatus and storage medium for blood flow dynamics of capillaries in brain tissue

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2008010951A JP4714228B2 (en) 2008-01-21 2008-01-21 Index calculation method, apparatus and storage medium for blood flow dynamics of capillaries in brain tissue

Related Parent Applications (1)

Application Number Title Priority Date Filing Date
JP2001318344A Division JP4216496B2 (en) 2001-10-16 2001-10-16 Index calculation method, apparatus and program code for blood flow dynamics of capillaries in brain tissue

Publications (2)

Publication Number Publication Date
JP2008100121A true JP2008100121A (en) 2008-05-01
JP4714228B2 JP4714228B2 (en) 2011-06-29

Family

ID=39434858

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2008010951A Expired - Lifetime JP4714228B2 (en) 2008-01-21 2008-01-21 Index calculation method, apparatus and storage medium for blood flow dynamics of capillaries in brain tissue

Country Status (1)

Country Link
JP (1) JP4714228B2 (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010068958A (en) * 2008-09-17 2010-04-02 Toshiba Corp X-ray ct apparatus, medical image processing apparatus, and medical image processing program
JP2011525408A (en) * 2008-06-24 2011-09-22 メドラッド インコーポレーテッド Identification of target area and extraction of time value curve in imaging method
EP2623457A1 (en) 2012-02-02 2013-08-07 VTU Holding GmbH Use of an ionic liquid for storing hydrogen
US8724869B2 (en) 2010-09-01 2014-05-13 Kabushiki Kaisha Toshiba Medical image processing apparatus
JP7357060B2 (en) 2019-02-06 2023-10-05 ウィリアム イー バトラー, Improved methods of angiography

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000050109A (en) * 1998-05-29 2000-02-18 Stmicroelectronics Inc Nonlinear image filter for removing noise
WO2000057777A1 (en) * 1999-03-26 2000-10-05 Oestergaard Leif Method for determining haemodynamic indices by use of tomographic data
JP2001061157A (en) * 1999-06-14 2001-03-06 Nikon Corp Image processing method, machine-readable recording medium with image processing program recorded therein and image processor

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000050109A (en) * 1998-05-29 2000-02-18 Stmicroelectronics Inc Nonlinear image filter for removing noise
WO2000057777A1 (en) * 1999-03-26 2000-10-05 Oestergaard Leif Method for determining haemodynamic indices by use of tomographic data
JP2001061157A (en) * 1999-06-14 2001-03-06 Nikon Corp Image processing method, machine-readable recording medium with image processing program recorded therein and image processor

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011525408A (en) * 2008-06-24 2011-09-22 メドラッド インコーポレーテッド Identification of target area and extraction of time value curve in imaging method
JP2010068958A (en) * 2008-09-17 2010-04-02 Toshiba Corp X-ray ct apparatus, medical image processing apparatus, and medical image processing program
US9561011B2 (en) 2008-09-17 2017-02-07 Toshiba Medical Systems Corporation X-ray CT apparatus, medical image processing apparatus and medical image processing method
US8724869B2 (en) 2010-09-01 2014-05-13 Kabushiki Kaisha Toshiba Medical image processing apparatus
EP2623457A1 (en) 2012-02-02 2013-08-07 VTU Holding GmbH Use of an ionic liquid for storing hydrogen
JP7357060B2 (en) 2019-02-06 2023-10-05 ウィリアム イー バトラー, Improved methods of angiography

Also Published As

Publication number Publication date
JP4714228B2 (en) 2011-06-29

Similar Documents

Publication Publication Date Title
JP4216496B2 (en) Index calculation method, apparatus and program code for blood flow dynamics of capillaries in brain tissue
JP4363833B2 (en) Method and apparatus for calculating an index relating to local hemodynamics
US7826885B2 (en) Method and apparatus for calculating index concerning local blood flow circulations
JP4170767B2 (en) Image processing device
US7853309B2 (en) X-ray CT apparatus and myocardial perfusion image generating system
JP4354484B2 (en) Imaging internal structures
US7873141B2 (en) X-ray tomographic imaging apparatus
JP2018538076A (en) Perfusion digital subtraction angiography
US9070181B2 (en) System and method for extracting features of interest from an image
US20050277830A1 (en) X-ray CT apparatus and myocardial perfusion image generating system
US8452067B2 (en) Method for generating bone mask
EP2875781A1 (en) Apparatus and method for processing a medical image of a body lumen
JP2011244875A (en) Image processor, image display device, program, and x-ray ct apparatus
BR112020013421A2 (en) diagnostic support program
JP4714228B2 (en) Index calculation method, apparatus and storage medium for blood flow dynamics of capillaries in brain tissue
JP4302180B2 (en) Method and apparatus for calculating an index relating to local hemodynamics
JP4864909B2 (en) Image processing device
JP2020512133A (en) Spectral imaging FFR
KR102399792B1 (en) PRE-PROCESSING APPARATUS BASED ON AI(Artificial Intelligence) USING HOUNSFIELD UNIT(HU) NORMALIZATION AND DENOISING, AND METHOD
KR100709806B1 (en) Image processing device
KR20190002960A (en) Method for task based CT reconstruction protocol optimization using body-recognition
CN110730977A (en) Low dose imaging method and apparatus
RU2815435C2 (en) Method of contouring abnormal zones on results of percussion tomographic scanning
RU2644928C1 (en) Method for determination of frames relating to circulation phases borders, during angiographic study (versions)
EP3995081A1 (en) Diagnosis assisting program

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20080220

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20101207

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20110207

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20110301

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20110325

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20140401

Year of fee payment: 3

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313114

Free format text: JAPANESE INTERMEDIATE CODE: R313117

Free format text: JAPANESE INTERMEDIATE CODE: R313111

R371 Transfer withdrawn

Free format text: JAPANESE INTERMEDIATE CODE: R371

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313114

Free format text: JAPANESE INTERMEDIATE CODE: R313117

Free format text: JAPANESE INTERMEDIATE CODE: R313113

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

S533 Written request for registration of change of name

Free format text: JAPANESE INTERMEDIATE CODE: R313533

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350