JP2003190148A - Method and device for operating index related to local hemokinesis - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To improve efficiency in the diagnosis of a map acquired by a CBS study. <P>SOLUTION: A first temporal concentration curve related to an artery within a specific site and a second temporal concentration curve related to a tissue within the specific site are prepared from a plurality of continuous images related to the specific site of an examinee injected with a contrast medium. A transmission function indicating the local hemokinesis of the tissue relative to respective arteries is calculated by curve fitting to minimize the residual of the second temporal concentration curve to the convolution of the transmission function and the first temporal concentration curve. An index related to the local hemokinesis relative to the respective arteries is calculated from the transmission function. The map of the indices corresponding to the arteries is prepared and the maps of the indices are synthesized to one map according to the residual relative to the respective first temporal concentration curves. <P>COPYRIGHT: (C)2003,JPO

Description

Detailed Description of the Invention

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to an index calculation method and calculation apparatus for local blood flow dynamics of brain tissue and the like.

[0002]

2. Description of the Related Art In an X-ray CT examination, morphological information can be obtained from a simple CT image, and dynamic information of blood flow around a lesion can be obtained as visual information by a dynamic scan by contrast CT. It is considered that in recent years, high-speed scanning by multi-slice CT has become possible, and the range of application of dynamic scanning of contrast CT is expanding more and more.

As one of the directions, there is a method called CBP study for calculating an index relating to blood flow dynamics of capillaries in brain tissue. In the CBP study, CBP that quantitatively shows the local blood flow dynamics in the tissue, that is, the dynamics of blood flow passing through capillaries in the local tissue,
Indexes such as CBV, MTT, and Err are obtained, and a map of these indexes is output.

[0004] CBP represents the blood volume per unit volume and unit time in the capillaries of brain tissue [ml / 100 ml / min], and CBV represents the blood volume per unit volume in brain tissue [ml].
/ 100 ml], MTT represents the average blood transit time [seconds] of the capillaries, and Err represents the sum of the residuals when the transfer function is approximated or the square root of the squared sum of the residuals.

Indexes CBP, CBV and MTT which quantitatively represent the blood flow dynamics of capillaries in these brain tissues.
Is expected to be useful information for evaluating the ischemic cerebral vascular disease pathology, the presence or absence of capillary expansion, and blood flow velocity, as well as information on the time elapsed since the onset of cerebral ischemic stroke. There is. For example, generally in ischemic cerebrovascular accidents, the blood pressure of the provided artery is lowered, and the blood flow velocity in the blood vessel is reduced. As a result, even if CBV is constant, MTT is extended and CBP is lowered. Also, in the hyperacute phase of cerebral infarction,
The blood flow CB is increased by expanding the capillaries and increasing the blood flow rate to compensate for the decrease in blood flow rate due to the decrease in blood pressure.
It has a function of trying to suppress the decrease of P (auto regulation). Therefore, by extending MTT,
Even if CBP is lowered, if CBV is increased, it is information that suggests the possibility of recanalization of capillaries.

In the CBP study, a contrast agent having no cerebrovascular 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 injected intravenously by the injector flows into the cerebral artery via the heart and lungs. Then, the contrast medium flows out from the cerebral artery to the cerebral vein through the capillaries in the brain tissue. At this time, the iodine contrast agent passes through the capillaries in the normal brain tissue without leaking out of the blood vessels. FIG. 1 schematically shows this state.

The state of passage of the contrast agent is 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 (capillary blood vessel) are obtained from the consecutive images. (T), time-density curve C of pixels on the cerebral vein
Each sss (t) is measured.

Here, in the CBP study, the concentration time curve Ca (t) of the cerebral artery and the concentration time curve Ci of the brain tissue are shown.
The ideal relationship that holds with (t) is used as the analytical model. When a contrast medium is injected from a blood vessel immediately before entering the brain tissue, the time-density curve within a unit volume (1 pixel) of the brain tissue has a vertical rising edge, maintains a constant value for a while, and then rises steeply. It goes down. This is approximated by a rectangular function (box-MTF method: box-Modulat
Ion Transfer Function met
hod).

That is, the time concentration curve Ca (t) of the cerebral artery
Is an input function and the time-concentration curve Ci (t) of the brain tissue is an output function, and the transfer function between the input function and the output function is approximated by a rectangular function. The transfer function represents the process by which the tracer passes through the capillaries.

The problems with the CBP study are as follows.
Since the indexes of CBP, CBV, MTT, and Err are calculated for each pixel (x, y, z), it is possible to form an image having the value as a pixel value, and map such an image. Call. For example, when R types of indexes are obtained, R maps can be constructed. The R maps thus created can be regarded as one map (vector value map) in which each pixel has a vector value. That is, it can be expressed as follows. V _k (x, y, z) = <P _{k, 1} (x, y, z), P _{k, 2} (x, y,
z), ..., P _{k, R} (x, y, z)> For example, in a CBP study, typically R = 4,
P _{k, 1} (x, y, z) is the value of CBP, P _{k, 2} (x, y, z) is the value of CBV, and P _{k, 3} (x, y, z) is the value of MTT. , P
_{k, 4} (x, y, z) can be configured to represent the value of the residual Err.

Such a vector value map V _k is created for each time concentration curve Ca (t) _k of the referenced cerebral artery. For example, if the time-concentration curve of the cerebral artery is obtained from the left and right middle cerebral arteries, the anterior cerebral artery, and the posterior cerebral artery, K = 6, and if the time-concentration curve of the cerebral artery is obtained from several places around the lesion, K = about 10-15.

Thus, the time concentration curve Ca of the cerebral artery
If the number K of (t) _k (k = 1, 2, ..., K) is large, the resulting vector value map V _k (k =
It is inconvenient to observe because of the large number of (1, 2, ..., K). That is, to observe as a normal gray scale image or color scale image, one map is composed of R images, and since there are K sets, a total of K × R images must be compared. . Furthermore, it is not always obvious which part is nourished by which artery, using anatomical knowledge,
Which map V _k (k = 1, 2, ...,
You have to decide whether to observe K). Especially in cases of cerebrovascular disease such as cerebral infarction, it is not always consistent with anatomical knowledge as to which artery depends on the tissue, and abnormal dominance is often seen. . Due to these problems, it is difficult to interpret the vector value map.

Further, in a dynamic CT image taken by multi-slice CT or volume CT, more arteries are observed. This is because the same artery can be observed in multiple slices. If the time-density curve of the cerebral artery is created for all the tomographic images of these arteries, the number becomes very large.

The CBP study also has the following problems. When bolus injection is performed in the elbow vein, the CT effect of CT observed in blood (tens of HU when not contrasted) rises up to several hundreds HU. However, in order to effectively analyze the cerebral blood flow, the contrast effect must be measured with an error within a few percent at most.
That is, the blood contrast effect (increase in CT value) is 20 to 4
It is necessary to be able to detect this even at 0 HU.

The volume ratio of capillaries in a unit volume of brain tissue is about 3 to 4% at most. 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.

In the CT image, the standard deviation of noise (sd)
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 the contrast effect of 0.5 HU, the X-ray dose must be increased about 10 to 100 times, which means that the exposure dose to the patient is significantly increased. Also, in dynamic CT, the same area is photographed several tens of times, so the exposure of the skin at the photographed area is several hundreds to several thousand times the usual exposure, resulting in inflammation, hair loss, necrosis, and the like. It is not realistic considering radiation damage such as carcinogenesis.

On the contrary, in the dynamic CT, the X-ray dose must be reduced as compared with the usual radiography. In general,
X-ray dose per scan is, for example, 1/2 to 1 / normal
It is reduced to about 10. This makes it possible to limit the exposure to X-rays to several times to 20 times that of normal CT imaging, which does not cause radiation damage. However, in such a CT image with reduced X-ray dose, sd is, for example, about 15 to 20 HU,
A contrast effect of about 0.5 to 1.5 HU cannot be detected at all.

Therefore, suppressing the noise component of the image is one of the important problems in the CBP study. Therefore, 1) increasing the slice thickness, 2) averaging neighboring pixels, and 3) passing the image through smoothing are generally possible measures. However, these have the following problems.

In order to "thicken the slice thickness", the slice thickness is set thick at the time of photographing, or several consecutive thin slice images are averaged to generate a thick slice image.
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 uniform brain tissue, and multiple tissues (white matter, gray matter, blood vessels, sulci, ventricles, and ventricles). The probability of representing an average CT value of () and the like increases, and the error of the value such as cerebral blood flow obtained as an analysis result increases.

Pixels including the influence of blood vessels cannot be normally analyzed. For this reason, increasing the slice thickness results in very poor quality results, including many inaccurate and unanalyzable pixels.

"Averaging neighboring pixels" sacrifices spatial resolution to some extent. For example, an average value of a square area (including n × n pixels) having one side of n pixels is obtained, and this is set as an average CT value of the entire square, and such a square is regarded as a pixel, and this is regarded as a pixel. "Panel bundling image" is constructed by laying out. If the original image is composed of 512 pixels on one side (including 512 × 512 pixels), and n = 2, the “pixel bundle image” is one side (512 pixels).
/ 2) is an image composed of pixels (including 256 × 256 pixels). According to this method, the noise is
It can be reduced in inverse proportion to n. further,
Since the number of pixels to be analyzed becomes 1 / (n × n) times, there is an advantage that the amount of calculation becomes small.

However, when n is increased, the spatial resolution is lowered, and the partial volume effect is generated accordingly, that is, one pixel does not represent a uniform brain tissue, and a plurality of tissues (white matter, gray The probability of representing an average CT value (quality, blood vessel, sulci, ventricles, etc.) increases, and the error in the values such as cerebral blood flow obtained as an analysis result increases. In particular, a pixel including the influence of blood vessels cannot be normally analyzed. For this reason, if n is increased, only a very poor quality result is obtained, which has a low spatial resolution, is inaccurate, and includes many pixels that cannot be analyzed. Therefore, practically, the limit is about n = 2 to 4, and this alone cannot provide a sufficient noise suppressing effect.

Further, smoothing of the image, that is, one CT
If a method of smoothing by applying a two-dimensional spatial filter to each image is used, the spatial resolution is significantly impaired in exchange for a sufficient noise suppressing effect. In particular, a pixel adjacent to a location where a thick blood vessel (artery / vein) exists is affected by the contrast effect produced in the thick blood vessel, and the time-density curve of these pixels becomes incorrect. Therefore, only a slight degree of smoothing has to be done. Here, what is important for performing a very slight smoothing is to make the size of the image filter extremely small, for example, to set it to about 3 × 3. In order to obtain the maximum image noise suppressing effect using the 3 × 3 smoothing filter, the upper limit is to reduce the noise sd to 1/3, and it is impossible to suppress the noise more than that. Is. Therefore, a sufficient noise suppressing effect cannot be obtained.

On the other hand, temporal smoothing, that is, the time density curve obtained for each pixel is regarded as a curve, and this is set to 1
When the method of smoothing with a dimensional filter is used, the time resolution is significantly impaired when trying to obtain a sufficient noise suppression effect. Originally, the dynamic CT in the CBP study obtains a high time resolution by capturing an image in a short sampling period, and a slight and fast change in the time-density curve (particularly, how smoothing effect caused by physiological structure is Since the purpose is to precisely measure (whether it has occurred), temporal smoothing is not appropriate at all.

[0025]

The object of the present invention is to provide CB
It is to improve the diagnostic efficiency of maps obtained by P-study. Another object of the present invention is to suppress a decrease in spatial and temporal resolution and suppress noise, thereby achieving CBP.
To improve the analysis accuracy of the study.

[0026]

According to a first aspect of the present invention, a first time concentration curve of an artery in a specific site is obtained from a plurality of consecutive images of a specific site of a subject injected with a contrast medium. , A second time-density curve for the tissue in the specific site is created, and a transfer function that represents the local hemodynamics of the tissue for each of the arteries is calculated using the second convolution of the transfer function and the first time-density curve. Calculation is performed by curve fitting so that the residual of the 2-hour concentration curve is minimized, an index relating to the local hemodynamics for each of the arteries is calculated from the transfer function, and a map of the index corresponding to the artery is created. However, the map of the index is combined into one map according to the residual for each of the first time concentration curves. A second aspect of the present invention relates to a first time-density curve for an artery in the specific site and a tissue for the specific site from a plurality of consecutive images of the specific site of a subject injected with a contrast agent. By creating a 2-hour concentration curve and selecting one of the first time-concentration curves that best fits each of the second time-concentration curves, it is most likely that the local hemodynamics of each of the tissues are dependent. A high one of the arteries is specified, and a map for distinguishing a subordinate region of the arteries is created based on the one of the arteries specified for each tissue. According to a third aspect of the present invention, in order to reduce noise from a plurality of consecutive images of a specific site of a subject injected with a contrast agent, the plurality of images are weighted according to the similarity between pixels. From the plurality of noise-reduced images, a first time-density curve for an artery in the specific site and a second time-density curve for tissue in the specific site.
A time-density curve is created, and a transfer function representing the local blood flow dynamics of the tissue with respect to each of the arteries is defined as the first transfer function and the first transfer function.
Calculating by curve fitting such that the residual of the second time-concentration curve with respect to convolution with the time-concentration curve is minimized, and calculating an index relating to the local hemodynamics for each of the arteries from the transfer function. is there.

[0027]

BEST MODE FOR CARRYING OUT THE INVENTION The present invention will be described below with reference to the accompanying drawings in preferred embodiments. A feature of this embodiment is that by combining many index maps generated by the CBP study into one sheet, the diagnostic efficiency of the CBP study is improved, and noise is reduced by using a coherent filter. It is to improve the accuracy of the index by satisfying both suppression of deterioration of spatial and temporal resolution.

The present embodiment relates to a method and apparatus for calculating an index representing local blood flow dynamics from a plurality of temporally consecutive images of a specific part of a subject, and a plurality of target objects The modality that generates the image of is not limited to a specific device.
-Ray computer tomography system (X-ray CT system), single photon emission tomography system (SPE)
CT), positron emission tomography device (P
ET) or a magnetic resonance imaging apparatus (MRI). Note that, here, an X-ray CT apparatus will be described as an example.

(Apparatus Configuration) FIG. 2 shows the configuration of the X-ray CT apparatus according to this 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 (D
AS; Data Acquisition System
m) and. The X-ray tube 101 and the X-ray detector 102 are mounted at positions facing each other with a subject P sandwiched in a rotating ring (not shown) that rotates at high speed and continuously.

The computer device 20 is the image processing device 3.
0, an image display unit 107, and an input unit 109. The image processing apparatus 30 has a control unit 108 as a center, and a preprocessing unit 104 that converts raw data output from the data collection unit 103 into projection data through a correction process and the like, a memory unit 105 that stores projection data, and projection data. To C
An image reconstruction unit 106 that reconstructs T image data, a storage device 10M that stores CT image data, a coherent filter processing unit 110 that executes coherent filter processing on CT image data, and a CT image that has undergone coherent filter processing. The CBP study processing unit 120 executes CBP study processing using data.

The coherent filter processing section 110 comprises a variance value estimation section 111, a weighting function calculation section 112, and a pixel value calculation section (coherent filter section) 113. These variance value estimation unit 111 and weighting function calculation unit 11
2. The function of the pixel value calculation unit 113 will be described in the detailed description of the coherent filter processing described later.

The CBP study processing unit 120 includes an ROI setting support unit 121, a time density curve creation unit 122, a cerebral artery time density curve correction unit 123, an MTF processing unit 124, an index calculation unit 125, a map creation unit 126, and a map synthesis unit. It is composed of 127.

The ROI setting support unit 121 is information (AT map, PT map, TT map for cerebral artery ROI) for supporting the work of setting the region of interest ROI for the cerebral artery and cerebral vein on the CT image. Etc.) are created and provided.

The cerebral artery ROI is, for example, the anterior cerebral artery (ACA), the middle cerebral artery (MCA), the posterior cerebral artery (PC).
Target area A) is set individually in the left and right brain regions. Therefore, in this example, six cerebral artery ROIs are set, three on the left and three on the right. Further, another time concentration curve Csss (t) is used to correct the time concentration curve Ca (t) of the cerebral artery. This time concentration curve Cs
ss (t) is created for a ROI set on a blood vessel that is thick enough that there are pixels that do not include partial volumes. The ROI of Csss (t) is, for example,
It is set in the uppermost sagittal sinus of the cerebral blood vessels.

The time-density curve creating section 122 includes the storage device 1.
A time-density curve for a cerebral artery, a cerebral vein, and a cerebral tissue (capillary blood vessel) is created from the dynamic CT image data (a plurality of temporally consecutive image data) stored in 0M. The time-density curve Ca (t) of the cerebral arteries is individually created for the set six cerebral arteries ROIs, for example. The time-concentration curve Csss (t) of the cerebral vein is created for the cerebral vein ROI set in the superior sagittal sinus.
Further, the time-density curve Ci (t) of the brain tissue is created for every pixel on all the pixels on the brain tissue.

The cerebral artery time-density curve correction unit 123 uses the time-density curve Ca (t) of the cerebral artery as the time-density curve Csss (t of the superior sagittal sinus in order to remove the influence of noise and the partial volume effect. ). This correction method will be described later. MTF processing unit 124
Is a box-based on the corrected time-concentration curve Ca (t) of the cerebral artery and the time-concentration curve Ci (t) of the brain tissue.
By the MTF method, the transfer function MTF is calculated for each pixel for all the pixels in the brain tissue region.

The index calculator 125 calculates the index (CBP, CBV, MTT, Err) representing the blood flow dynamics of the brain tissue from the calculated transfer function MTF for each pixel for all the pixels in the brain tissue region. . The map creation unit 126 assigns the maps of the calculated indexes to
It is generated for each cerebral artery (ACA, MCA, PCA). The map shows the index type (=
4) x number of cerebral arteries (three of ACA, MCA, PCA) =
12 types are created. In multi-slice, the number of types of maps is increased. A map synthesizing unit 127 is provided to reduce the enormous number of maps by the synthesizing process and improve the diagnostic efficiency.

The following is the coherent filtering process and the CB.
The P study process will be described in order. The principle of the coherent filter will be briefly described. On the basis of weighted average of pixels in the neighborhood such as 3 × 3, and using the weighted average value as the value of the local center pixel, the weight of each peripheral pixel is centered. It is characterized in that it is changed according to the degree of similarity between a pixel and a peripheral pixel. The similarity referred to here is an index indicating the degree of possibility that the tissues are close to anatomical, specifically brain tissues (capillaries) under the control of the same cerebral artery, between pixels, By giving a high weight to pixels with a high degree of similarity and conversely, giving a low weight to pixels with a low degree of similarity close to zero, while suppressing noise, it is possible to suppress a decrease in spatial resolution. It is possible. In the following, the similarity is
Appropriately translated into fitness or risk.

In the present embodiment, a plurality of CT images (dynamic CT) obtained by continuously capturing the brain of a subject injected (intravenous) with a contrast medium having no cerebrovascular permeability, for example, an iodine contrast medium (dynamic CT). Image) is used to calculate the degree of similarity by comparing the time density curves of each pixel. Therefore, the certainty of the similarity depends 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.

(Coherent Filter) (General Description of Coherent Filter) (Pixel Value v (x)) Generally, a digital image acquired through an image pickup means such as a camera or a CT scanner is composed of a plurality of pixels (pixels). Configured (or the image can be thought of as a collection of such pixels). In the following description, the position of the pixel is represented as a vector x (that is, a vector of coordinate values), and the value of the pixel x (for example, a numerical value representing grayscale, the CT value HU) is represented as a K-dimensional vector. In the case of a two-dimensional image, the pixel x is a two-dimensional vector indicating the coordinate value (x, y) indicating the position on the image.
The “pixel value v (x)” defined for a certain pixel x is expressed as v (x) = (v ₁ (x), v ₂ (x), ..., V _K (x)) (1) . Note that v _{1 on} the right side of the equation (1)
Each of (x), v ₂ (x), ..., V _K (x) is hereinafter referred to as a “scalar value” for the pixel x.

For example, if the image is a "color image"
The brightness of each pixel is of three primary colors (red, green, blue).
Since it has (scalar value), the pixel of each of these pixels
Consider the value v (x) to be a vector whose dimension is K = 3.
Can be obtained. That is, each term on the right side of the above equation (1)
The letters are associated with, for example, "red", "green", and "blue".
It Further, for example, the image is composed of K still images.
Each pixel in the nth image is a moving image and is a scalar value.
v_nIn the case of having (x), on K still images,
Pixel value (scan value) of pixel x at the same common point (same coordinate)
Error vector), which is a K-dimensional vector value v
_n(X) = (v₁(X), v_Two(X), ..., v _K(X)
Is a pixel value as a vector value described below.

(Similarity (fitness or danger rate) and weight) Consider a set N (x) of appropriate peripheral pixels for the pixel x (this set N (x) may include the pixel x). ). Next, consider the weight w (p (x, y)) of the peripheral pixel y which is an element of N (x) with respect to the central pixel x. This weight w (p
(X, y)) has the following properties.

(Similarity p (x, y)) First, the weight w (p
The meaning of the function p (x, y) that influences the value of (x, y) will be described. This p (x, y) is a means for quantifying the “similarity” in the present embodiment, and generally speaking,
To what extent the central pixel x and the peripheral pixel yεN (x) are similar in some sense (for example, the statistics observed between the pixel values v (x) and v (y) of the both pixels x and y). The specific numerical value indicating the degree of difference) is given.

More specifically, for example, when p (x, y) shows a large value, "there is no statistically significant difference (that is, the similarity is high)" between pixel x and pixel y. , P (x, y) is determined to have a high possibility of being similar, and when p (x, y) has a small value, there is a statistically significant difference between the pixel x and the pixel y (that is, the similarity is low. ) ”, Is determined.

By the way, pixel values v (x) and v (y)
(Or scalar values v ₁ (x), ..., V _K (x) and v
₁ (y), ..., V _K (y)) always contains noise. For example, considering the case where an image is acquired by a CCD image sensor, for each pixel that constitutes it,
There are noises and the like due to dark current in the element and irregular fluctuations of the amount of light incident from the outside.

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 external world), the noise is actually obtained. It may not have the same value on the observed image. Conversely speaking, assuming that the noise is removed from the pixel x and the pixel y, both of which reflect 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.

Therefore, in consideration of the above-mentioned characteristics of noise,
Regarding the above p (x, y), using the concept of the "null hypothesis" well known in the statistical test method, this p (x, y)
Regarding y), it can be specifically stated as follows. That is, the null hypothesis H "pixel x and pixel y have the same pixel value when their respective noises are removed". In other words, "v (x) = v (y),""Exclude caused difference" (that is, when such a proposition holds, consider that "the similarity between both pixels x and y is high (the degree of matching is large)") and the function p ( x, y) can be configured as a risk rate (or a significance level) when rejecting this hypothesis H (in this case,
p (x, y) is defined as a function whose range is [0,1] (p (x, y) ε [0,1])).

Therefore, it can be said that the above hypothesis H is likely to be satisfied when the risk rate p (x, y) is large, that is, when the rejection is erroneous, and when it is small, that is, the rejection is rejected. It can be said that there is a high possibility that the hypothesis H will not be satisfied if the risk of being false is small. (Although it is a well-known matter in statistical tests, even if the hypothesis H is not rejected, it is It does not mean that it is “true.” In this case, it just means that the proposition given by Hypothesis H cannot be denied.)

(Weight w (p (x, y))) Now, the weight w
(P (x, y)) is, as is clear from the way it is expressed,
A function of the risk factor p (x, y) as described above (more generally, a fitness function (if the fitness is ρ (x, y), w
(It can be configured to be ρ (x, y)), and in order to obtain this weight w (p (x, y)), the risk rate p obtained for each combination of x and y
Generally speaking, the weighting function w acting on (x, y) has a function of embodying the above-mentioned “rejection”. Specifically, when the risk rate 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 in the opposite case, it is small. The value is adjusted to take a positive value (or “0”) (the specific form of the weighting 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 are likely to satisfy the proposition shown in the above hypothesis H, and takes a small value in the opposite case. As an example, in particular, it may be configured such that there are only two possible values of w that are “0” or a constant value that is not “0”.

The above-mentioned hypothesis H and the risk rate p
Summarizing the relationship between (x, y) and weight w (p (x, y)), when the null hypothesis H is highly likely to be correct, the similarity p
Becomes higher, and the weight w given to the pixel is increased, while
When the null hypothesis H is unlikely to be correct, the similarity p is also low, and the weight w given to the pixel is low. By changing the degree of contribution (weight) to the weighted average value according to the degree of similarity in this way, it is possible to effectively suppress noise while suppressing a decrease in resolution. Also, the weighting function w
(T) can be more generally referred to as a “non-negative monotonically increasing function defined by tε [0,1]”, and the w
The property to be satisfied by (t) should be at least such a property.

(Coherent Filter Processing) From the above description, the “coherent filter” is derived as follows. That is, first, for a certain pixel x forming an image, the weight w (p (x, y)) is calculated for all the pixels y that are elements of the set N (x). Next, using the plurality of weights w (p (x, y)), the pixel x
A new scalar value v ′ _k (x) that constitutes the above is calculated by the following equation (2). That is,

[Equation 1] However, k = 1, 2, ..., K. Then, using v ′ _k (x) obtained by this equation, the pixel value (new pixel value) v ′ (x) of the pixel x after conversion is v ′ (x) = (v ′ ₁ (X), v ′ ₂ (x), ..., V ′ _K (x)) (3).

Here, the pixel value v (y) = (v ₁ (y), v ₂ (y), ..., V represented by the equation (2) above.
Let _K (y)) (including the case where y = x) be v ′.
(X) = (v ′ ₁ (x), v ′ ₂ (x), ..., V ′
_The filter that transforms into _K (x)) is in the form of a "coherent filter". As is clear from the expression, this represents the weighted average value of the scalar values v _k (y) that form the pixel value.

Such processing has the following results.
Let me down. That is, the pixel value v of the pixel x from which noise is removed
It is certain that it will take the same pixel value as ´ (x) (= above
The probability of satisfying the proposition of hypothesis H is high) The contribution of pixel y
Weighted average value v ' _k(X)
It will represent Koutor. In addition, the number of such pixels y is 10
If there is a sufficient number, the pixel value v ′ (x) is
I have mentioned above without deviating from the true value that it should have
Noise is suppressed by the effect of averaging, but it is also decomposed.
Noh will have a value that does not decrease.

Note that the risk rate p (x, y) is small, and therefore the null hypothesis H is "rejected", and the weight w (p (x, y
Even if y)) becomes small, as is clear from the above description, this is not always completely “rejected”. This depends on the specific form of the weighting function w, which will be described later, but the risk factor p (x,
Even when y) is close to “0” (= 0%), w (p
(X, y)) ≠ 0 (provided that p (x, y) is a smaller positive value than when p (x, y) is close to “1”) (note that p (x, y) = 1) The case of is when v (x) = v (y), as will be described later.).

That is, it is possible to construct so as to allow a small contribution rather than a complete rejection (in such a case, w (p (x, y))).
If = 0, it is synonymous with performing a complete rejection.

Such processing can generally be said as follows. That is, when there are a plurality of pixels x forming a certain image, the matching degree p (x, y) between this pixel x and a certain arbitrary pixel y (in the above, yεN (x)) is quantified. If the goodness of fit is large, a large contribution of the pixel y is recognized in the weighted averaging process using the pixel value v (y), and if the goodness of fit is small, only a small contribution is recognized. By doing so, it can be said that the image processing method effectively suppresses the noise of the pixel x. In other words, when the pixel x and the pixel y are “similar things”, the pixel y is more contributed to the averaging process, and when the pixel “is not similar”, the pixel y is ignored or hardly ignored. May be paraphrased (the weight is zero or its approximate value).

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 decrease in spatial resolution. Further, not only for the purpose of noise suppression, but also in the field of pattern recognition, for example, the weighting function or the coherent filter can exhibit an excellent effect by making it into a suitable concrete form.

The above-mentioned "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, continuous rotation C).
In the T-apparatus, repetitive photographing by continuous rotation is often performed. ) And obtain projection data one after another, and perform a reconstruction process one after another based on the projection data to obtain a series of images in time series (in this case,
The image display on the image display unit 107 is controlled by, for example, a counter (not shown) or the like so as to be performed after a predetermined time from the scan start point or the end point related to the projection data acquisition that is the source of the image).

Therefore, the image thus obtained and displayed is a so-called moving image, which is composed of a plurality of time-series still images like a movie or the like. Note that such an imaging method typically involves injecting a contrast agent into the subject P,
It is used for observing and analyzing the change over time to analyze the pathological condition of other lesions such as stenosis and occlusion in blood vessels. Further, a method of performing CT imaging of the same site only twice before and after administration of a contrast agent can be considered as dynamic CT imaging in a broad sense.

In the prior art, during the "dynamic CT" imaging as described above, for example, some change in the subject P (for example, a change in the concentration of the contrast agent, respiratory motion, etc.) is generally performed while performing K times of imaging. However, in order to suppress image noise without spoiling the spatial resolution, there is no choice but to perform smoothing in the time direction. As a result, the adverse effect of impairing the time resolution was unavoidable.

However, the image obtained by the dynamic CT imaging is a moving image as described above, and is taken for the purpose of observing the temporal changes in detail, so that the time resolution is impaired. That is not a favorable situation.

If a coherent filter is used, the following dynamic coherent filter processing that can suppress noise from all K still images (a plurality of images) without impairing the temporal resolution is performed. can do.

First, regarding the pixel x defined for the K still images which are moving images obtained as described above,
As described above, v (x) = (v ₁ (x), v ₂ (x), ..., V _K (x)) (1) is reconfigured as the pixel value v (x). it can. Here, subscripts 1, 2, ..., K in each term on the right side are serial numbers of the K still images.

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

[0065]

[Equation 2]

However, yεN (x), and this set N (x) may be arbitrarily set for each pixel x (= any reference may be set). However, in reality, the pixel x and the pixel y located far from the pixel x can satisfy the hypothesis “v (x) = v (y). However, the difference due to noise of both pixels is excluded”. Since it can be said that the property is generally low, limiting the set N (x) on the basis of a set of pixels close to x has a practical significance such as an improvement in calculation speed.

Therefore, here, as an example, the set N (x) is a set of pixels included in a rectangular area around the pixel x. More specifically, as the set N (x), for example, when all the pixels forming one still image of interest are 128 × 128 pixels, 3 × with the pixel x as the center is set. It may be an area for 3 pixels, or in the case of 512 × 512 pixels, an area for 13 × 13 pixels around the pixel x may be used.

Further, σ _k in the above equation (4) is the standard deviation of noise estimated assuming that each pixel of the k-th still image has a certain degree common to all the pixels. On the other hand, C is an adjustable parameter that determines the degree of action when the weight w1 (p (x, y)) is substituted into the above equation (4).

Below, these σ_kAnd the explanation about C
Do in order. First, σ in equation (4)_kExplain about
(In the following, the variance σ_k ^TwoExplained as). This σ_k
^TwoIs, as described above, for each pixel on the k-th still image.
It is the variance of the noise component that the scalar value has. And
Also, the variance σ in the above equation (4)_k ^TwoIs the kth image
Variance σ that is a constant value for the scalar value of each pixel of_k ^TwoTo
Estimated assuming that it contains the noise
is there. In general, such assumptions are as follows:
Against the backdrop, there is sufficient justification.

When the size of the subject P, the structure of the X-ray tube 101, the X-ray detector 102, the reconstruction unit 106, and the like are constant and the energy of the irradiation X-rays is constant, the noise of the CT image is increased. Is determined by the irradiation X-ray dose, that is, the product of the tube current and the irradiation time in the X-ray tube 101 which is proportional to this (so-called tube current time product (mA · s)).

On the other hand, noise in CT images is additive,
It is also known to follow a Gaussian distribution. That is,
An arbitrary scalar that makes up the pixel value v (x) of a pixel x
Value v _nThe truth about (x) (n = 1, 2, ..., K)
The value (value from which the contribution of noise is removed) is v_n ⁰(X)
Then, the difference value v_n(X) -v_n ⁰(X) is roughly
Mean 0, variance σ_k ^TwoGaussian distribution of
X-ray dose or tube current time product mA · s and noise variance σ_k
^TwoAnd are almost in inverse proportion to each other. ).

The variance σ _k ² is the position of the pixel x itself (as described above, for example, each coordinate value x = (x,
y)), however, in the normal X-ray CT apparatus 100, the X-ray tube 101 and the X-ray detector 102 have an X-ray detector.
Ignore this because it has a physical X-ray filter (for example, a so-called "wedge" or "X-ray filter" made of copper foil, metal block, etc.) that adjusts the radiation dose. can do. This is because the wedge is irradiated so that the same amount of X-ray dose can be detected by any X-ray detector 102 by utilizing the fact that the subject P is made of a substance having almost the same density as water. It has the effect of attenuating part of the X-rays. Therefore, such a wedge results in the effect of making the noise variance σ _k ^{2 a} substantially constant value that hardly depends on the position of the pixel x (by the way, this wedge is
Generally, it is installed for the purpose of effectively utilizing the dynamic range of the X-ray detector 102).

From the above, on the K still images acquired by the dynamic CT imaging, the variance σ _k ^{2 is} calculated for all the pixels on the kth still image.
It is reasonable to assume that is almost constant. Of course, it can be easily considered that the present embodiment can be extended to the case where the variance is different for each pixel.

Next, in order to specifically calculate the equation (4), what numerical value should be applied as the variance σ _k ² becomes a problem. This is a problem because the shape of the noise distribution can usually be assumed (Gaussian distribution above), but the specific value of the variance σ _k ² is often unknown.

Further, in general, the irradiation dose (X-ray tube current × irradiation time (mA · s)) may be changed for each photographing, and in this case, the variance σ _k ² is different.

Now, the kth image (k = 1, 2, ...,
In K), the variance of the noise of the scalar value of each pixel is σ _k ^2, and the irradiation dose used to capture the kth image is R
_{When k} , σ _k ² is proportional to R _k . Therefore, if σk _O ² can be specified for at least one k = k ₀ , for other k,

[Equation 3] Allows σ _k ² to be accurately estimated.

In the present embodiment, at least one
For k, σ _k ^TwoOf the concrete numerical value of
An estimate can be made.

Of the K times of photography, there is almost no change in the subject P.
N times (1 <N ≦ K) that can be assumed that
The variance σ is measured by using the image._k ^TwoExpected value for
E [σ_k ^TwoIs effective. Simple explanation below
In order to obtain
And therefore σ for k = 1, 2, ... N_k ^TwoIs one
Constant (σ^TwoWrite). In these N images
Some pixel x_fPixel value v (x_f) Each ska
Error value v₁(X_f), V_Two(X_f), ..., v _K(X_f)
As described above, the noise included in the^TwoMoth
These mean values are
Equation (6) below,

[Equation 4] , The expected value E [σ ² ] for the true variance σ ² is

[Equation 5] Can be asked as Then, as described above, the expected value E [σ ² ] of the variance can be considered to be valid for all the pixels x on all the K still images, and is used as a substitute for the true variance σ ^2. , It is a value that guarantees a certain degree of certainty. Therefore, in the actual calculation of the above formula (4), this E [σ ² ] is expressed by σ of the formula (4).
Substitute in ² .

More specifically, such E [σ ² ] may be obtained by an actual measurement value based on, for example, the first and second still images in K still images (above). This is equivalent to setting N = 2 in the expressions (6) and (7). As for the pixel x _{f used} for the actual calculation of the above equations (6) and (7), for example, only an appropriate pixel x _f excluding a portion where air or bone is imaged is selected (plurality). If selected, the average of all obtained E [σ ² ] is taken) and the like. In addition, in general, it is better to take measures to suppress the influence of the movement of the subject P.

Irradiation dose in capturing these N images
Σ is not constant_k ^TwoIs R_kProportional to
Correctly using that_k ^TwoIt is easy to estimate
I can think of it.

Next, the parameter C in the equation (4) will be described. First, in the equation (4), the concept of the risk rate p (x, y) described in the above general form is included as follows. That is, the expression in the radical in the numerator on the right-hand side of equation (4) corresponds to the χ ² value that is said to follow the so-called χ ² distribution, and this is divided by (2σ) ² to obtain the entire parentheses. The value that is placed on the shoulder of e is the risk rate p1 (x, y) itself. That is,

[Equation 6] Then, the above equation (4) is related to p1 (x, y) represented by this equation (8),

[Equation 7] It is nothing but In addition, A is a constant and p1 is (0 to
It is standardized to have the value of 1).

After all, in the equation (4), the risk rate p (x, y) described in the general form as described above is not explicitly displayed, but the weight w1 (p (x, y) ), As described above, is exactly the risk rate (= p1 (x,
y)) can be regarded as a function (equation (9)), that is, a function of fitness (however, the risk rate and the fitness are, as described above, the other one, the more the other increases). There is an increasing relationship).

Then, 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 rate p1 (x, y). That is, when C is increased, p1 (x, y)
W1 (p (x, y)) approaches 0 with only slightly decreasing. Further, if C is made small, such a hypersensitive reaction can be suppressed. It should be noted that C may be specifically about 1 to 10, and preferably C = 3.

In the present embodiment, the similarity determination between the central pixel x and the peripheral pixel y, in other words, the rejection determination of the null hypothesis H described above regarding both pixels x and y is clear from the above description. Thus, the risk rate p1 (x, y)
It is determined based on the so-called χ-square test method (statistical test method).

As can be seen from the expression (4) above, in the present invention, after calculating the risk rate p (x, y) for each combination of x and y, the weight w (p
It is not always necessary to take the procedure of obtaining (x, y)), and it is also possible to directly calculate (w ° p) as a composite function without specifically obtaining the risk rate p (x, y). Good.

As described above, the variance σ^TwoAnd estimate
(For example, E [σ in equation (7)^Two]) And parameter C
By appropriately determining (for example, C = 3), (4)
A set N (x) defined for a pixel x using the formula
(As described above, for example, a 3 × 3 image centered on the pixel x
For all pixels y included in the prime area, etc.,
A specific weight w1 (p (x, y)) can be obtained.
It After that, in w (p (x, y)) in the above formula (2),
Instead, by using this w1 (p (x, y))
To perform specific numerical calculations for coherent filters.
It is possible to And as a result, time resolution
That is, noise is strongly suppressed without impairing the spatial resolution.
Restricted pixel value v ′ (x) = (v ′₁(X), v '
_Two(X), ..., v ' _K(X)) (= equation (3))
To obtain such K still or moving images
You can

FIGS. 3 (a) to 3 (c) show such image processing for easy conceptual understanding. That is, first, in FIG.
..., in a still image of K sheets, for a certain pixel x, a rectangular area N of 3 × 3 pixels centered on the pixel x
^{3 × 3} (x) is assumed. This rectangular area N
If the pixel at the left corner of ^{3 × 3} (x) is y ₁ , then
This pixel y ₁ has a pixel value v (y ₁ ).

Then, the pixel value v (y₁)
Scalar value v₁(Y₁), V_Two(Y ₁), ..., v_K(Y
₁) And the scalar value v at the pixel value v (x)₁(X),
v_Two(X), ..., v_KWith each of (x),
The weight w1 (p (x, y₁)) Is calculated
(FIG. 3 (b)). Also, the rectangular area N
^3x3The remaining pixel y of (x)_Two, ..., y₈Also for
Finally, as shown in FIG. 3B, w1 (p (x,
y₁)), ..., w1 (p (x, y₈)) And w1 (p
(X, x)) is obtained. In this case, it is more dangerous than equation (8)
The rate p (x, x) is "1" and therefore the weight w1
(P (x, x)) is also “1” according to the equation (9) (= maximum)
It is heavily weighted).

Next, the weight w1 thus obtained is obtained.
(P (x, y₁)), ..., w1 (p (x, y₈)), W
1 (p (x, x)) is the kth image of the corresponding pixel
The scalar value v at_k(Y₁), V_k(Y_Two),…,
v_k(Y₈), V_k(X) is multiplied by each to get the sum.
(Corresponding to the numerator in the above formula (2))
Rectangular area N^3x3Sum of weights w1 for (x)
(Similarly, it corresponds to the denominator of equation (2).)
Noise for pixel x in the kth image is
Suppressed scalar value v '_kCan determine (x)
(FIG. 3 (c)). Also, all of k = 1, 2, ..., K
Of the same weight w1 (p (x, y ₁)), ...,
w1 (p (x, y₈)), W1 (p (x, x))
, The scalar value v ′ with suppressed noise_kFind (x)
Image in which noise in the pixel x is suppressed by
Elementary value v '_k(X) = (v '₁(X), v '_Two(X),
…, V '_K(X)) is obtained. For every pixel x
If the above calculation is repeated, K images with suppressed noise will be displayed.
The image is obtained.

In the image composed of the pixel values v '(x) calculated by the coherent filter in this way, the random noise seen in the original image is sufficiently suppressed.

The processing described above may be performed in accordance with the flowcharts shown in FIGS. 4A and 4B, and the calculation / image In order to realize the display and the like on the actual X-ray CT apparatus 100, for example, as shown in FIG.
The image processing unit 110 including the weight calculation unit 112, the weight calculation unit 112, and the pixel value calculation unit 113 may be provided and implemented.

Of these, the weight calculator 112 directly weights the weight w1 from the pixel values v (x) and v (y) according to the procedure described above.
The configuration is such that (p (x, y)) is obtained. Therefore, the calculation unit 112 is a device that directly calculates the weight without specifically calculating the value of the risk rate p1 (x, y). It should be noted that, instead of the configuration as described above, a risk rate calculation unit (fitness quantification unit) that specifically calculates the value of the risk ratio p1 (x, y) and the weight w1 (p (x,
The configuration may be a two-step procedure called a weight calculation section for obtaining y)). In any case, the weight calculation unit 112
The variance σ ² estimated by the variance value estimation unit 111, and v
Weight w1 (p (x, y)) using (x) and v (y)
To calculate.

Further, the pixel value calculation unit 113 determines that the pixel value v
The pixel value v ′ (x) is calculated using (x) and v (y) and the weight w1 (p (x, y)) numerically calculated by the weight calculator 112. That is, the calculation unit 113
Performs the process of suppressing the noise of the original image, that is, actually applies the coherent filter (hereinafter, this is referred to as “coherent filter is applied”).

In the dynamic coherent filter processing as described above, when a coherent filter is applied to a moving image composed of K still images, the processing in the image processing section 110 is performed once for all still images. After being configured, these may be stored in the storage device 10M and a coherent filter may be applied to them later as post-processing, but the present embodiment is not limited to such a form, and is described above. The process of applying a coherent filter may be performed in real time in the sequence of continuous scan, continuous projection data acquisition, continuous reconstruction, and continuous display (hereinafter, this is referred to as “real-time coherent filter process”).

Real-time coherent filtering
In the preferred embodiment of, a new image is taken.
Each time it is reconfigured, the following processing is performed. At first
The latest image (image number) from the obtained image (image number 1)
Up to M), the image numbers M, M-1, ..., M-K + 1
Of the same point (same coordinate) on K still images with
Pixel values (scalar values) of pixel x are arranged and K-dimensional
Toll value v (x) = (v _M(X), v_M-1(X), ...,
v_{M-K + 1}(X)). Thus, above
Exactly the same as "dynamic coherent filtering"
You can apply a coherent filter like this. However
However, the pixel value calculation unit 113 actually calculates the pixel value v ′ (x)
Instead of calculating all elements, the latest image (image number
Scalar value v corresponding to No. M)_MCalculate only ´ (x)
It As a result, the calculation speed is improved, so real time
The latest image with suppressed noise can be displayed.

As another preferred embodiment of this "real-time coherent filter process", when the first K images are obtained, a coherent filter is applied in exactly the same manner as described above, v ₁ ′ (x) ,. , V _K ′ (x) is obtained, and thereafter, the K-dimensional vector value is set to the image number M, M.
V using K still images with -1, ..., M-K + 1
_{(X) = (v M (} x), v M-1 '(x), ..., v
_{M−K + 1} ′ (x)), and the above real-time coherent filter processing may be applied thereto. It should be noted that when these real-time coherent filter processes are performed, the pixel value vector v
It is convenient if the dimension K of (x) can be changed at any time by manual setting or automatic setting.

Thus, with the coherent filter,
A CBP study was performed using a CT image in which only noise was effectively suppressed without reducing spatial and temporal resolution, and blood flow passing through local blood flow dynamics in a tissue, that is, capillaries in a local tissue was performed. The flow dynamics are analyzed quantitatively, and the indices (CBP, CBV,
By obtaining MTT and Err), improvement in accuracy and reliability can be expected.

As described above, the CBP study process is performed on the image from which the noise is removed while suppressing the deterioration of the resolution. (CBP Study) As described above, in the CBP study, the indexes of CBP, CBV, MTT, and Err that quantitatively represent the dynamics of “blood flow through capillaries” in brain tissue are obtained, and the space of these indexes is also calculated. Output a map showing the statistical distribution.

CBP: blood volume per unit volume and unit time in capillaries of brain tissue [ml / 100 ml / min] CBV: blood volume per unit volume in brain tissue [ml / 10]
0 ml] MTT: Mean blood transit time of capillaries [seconds] Err: Deviation of the actual measurement value from the analytical model Residual index Note that the low residual Err may be controlled by the reference cerebral artery. Means high, and conversely the residual Er
A high r means that it is unlikely to be dominated by the reference cerebral artery.

In the CBP study, a contrast agent having no cerebrovascular permeability, for example, an iodine contrast agent is used as a tracer. The iodine contrast medium rapidly injected from the cubital vein by the injector flows from the cerebral artery via the heart and lungs. Then, it flows out from the cerebral artery through the capillaries in the brain tissue to the cerebral vein. At this time, a contrast agent having no cerebrovascular permeability, for example, an iodine contrast agent, passes through the capillaries in the normal brain tissue without leaking out of the blood vessel.

The state of passage of the contrast agent is continuously photographed by dynamic CT, and the time-concentration curve Ca (t) of the pixels on the cerebral artery and the time of the pixels on the brain tissue including capillaries are taken from the consecutive images. The concentration curve Ci (t) and the time concentration curve Csss (t) of pixels on the cerebral vein are measured.

In the CBP study, the blood concentration of the contrast agent, the time curve Ca of the blood concentration in the cerebral blood vessels close to the brain tissue,
An ideal relationship that holds between (t) and the time curve Ci (t) of the blood concentration in the capillaries is adopted as an analytical model. That is, when a contrast agent is injected from a blood vessel immediately before entering the brain tissue, the time-density curve in the brain tissue unit volume (1 pixel) including the capillary rises vertically, maintains a constant value, and has a slight gradient. Fall down. This can be approximated by a rectangular function (box-MTF method: box-
Modulation Transfer Funct
ion method).

The transfer function represented by a rectangular function is used to describe the characteristics of the process of passing through the capillaries using the time concentration curve Ca (t) in the cerebral artery blood as an input function and the time concentration curve Ci (t) of the brain tissue as an output function. Can be asked as

(Specific Procedure) FIGS. 5 and 6 show a typical procedure of the CBP study according to this embodiment.
First, a bolus injection (contrast agent is administered all at once) to a blood vessel such as an elbow vein, and dynamic CT immediately after or immediately before (imaging repeatedly at the same location)
I do. As the most typical procedure, when performing bolus injection into the elbow vein, it takes about 20 to 40 seconds.
For example, photographing is repeated at intervals of 0.5 to 2 seconds. The CT value of each j-th pixel (x, y) of the N CT images obtained by the dynamic CT is v (x, y, j). This is nothing more than the sampling of the time-density curve (smooth curve) f (t, x, y) at this pixel (x, y).

First, as preprocessing, in step S1, C
Pixels that are clearly determined to be tissues other than brain tissue are excluded from the analysis target from each T image. That is, the range of CT values of brain tissue (eg CT
Pixels showing values that do not fall within the range of 10 to 60 HU) are pixels corresponding to air, bone, fat, etc., and are not related to the quantification of cerebral blood flow, so these can be ignored. This analysis range is
As a default, it is set to 10 to 60 HU, but it can be arbitrarily set through the input unit 109.

As a preprocess, the contrast effect is initialized in step S2. In order to obtain the contrast effect (increase in CT value) at each pixel, for each pixel (x, y), an image before the contrast agent reaches the tissue corresponding to that pixel (generally multiple images are obtained) Serial number 1,
When expressed as 2, ... K, the temporal average value is

[Equation 8] Is obtained, and this value is set as b (x, y). And j = K
+1, K + 2, ..., N pixel values v (x, y,
For j), q (x, y, j) = v (x, y, j) −b (x, y) For j <K, q (x, y, j) = 0. To simplify the processing, the same K may be adopted for any pixel. Q thus obtained
(X, y, j) is a smooth continuous curve time-density curve q
It can be considered that nothing but (t, x, y) is sampled at t = t1, t2, ... tN.
Quantitative analysis of cerebral blood flow is performed using this q (t, x, y).

In the quantitative analysis, first, the right brain area and the left brain area can be separated on the CT image. As described above, in the CBP study, the state of blood flow dynamics of capillaries is obtained as the transfer function MTF of the time concentration curve Ci (t) of the brain tissue with respect to the time concentration curve Ca (t) of the cerebral artery. If the brain tissue to be analyzed is not under the control of the cerebral artery of the reference curve Ca (t), the calculation is useless. At least the left and right brains are separately analyzed by using different time-concentration curves Ca (t) of the cerebral arteries, that is, the time-concentration curve Ca (t) of the cerebral arteries of the left cerebrum is only used to analyze the brain tissue of the same left cerebrum. Similarly, using the time-concentration curve Ca (t) of the cerebral artery of the right brain only for analyzing the brain tissue of the same right brain has the effect of reducing unnecessary calculation.

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 superposed as a figure on the screen and displayed on the CT image (S3). The dividing line may be initially displayed in the center of the image. The operator divides the left and right areas by referring to the image, moving the dividing line, and moving and arbitrarily bending a plurality of constituent points forming the dividing line.

By dividing the brain into the left-brain area and the right-brain area and limiting the analysis range to each area in this way, the number of analysis processing steps can be reduced. That is, the time-density curve Ca (t) of the cerebral arteries of the left brain (left ACA, left MCA, left PCA) is used only for the analysis of the left brain area (transfer function optimization processing), and similarly, the cerebral arteries of the right brain (right ACA ,
The time concentration curve Ca (t) of the right MCA and right PCA) is used only for the analysis of the brain tissue of the same right brain. In order to reduce the number of analysis processing steps, the left brain area and the right brain area may be further divided into areas, respectively, and the analysis processing may be limited to a narrower area.

A geometrical drawing method, for example, Voronoi drawing method is adopted for the region division. As is well known, the Voronoi projection method is a method that is often used in fields such as optimal placement of facilities such as hospitals, stores, and fire stations, and it is possible to use multiple points (corresponding to stores, mother points) arranged on a plane. It is characterized in that the plane is divided into a plurality of influence zones according to the distance.

As shown in FIG. 8, in this embodiment, the Voronoi projection method is applied individually to the left brain area and the right brain area. The left brain area is divided into a power area of the left ACA, a power area of the left MCA, and a power area of the left PCA with the left ACA, the left MCA, and the left PCA as three mother points. Left ACA,
A Voronoi point is set at the center of a circle passing through three generating points corresponding to the left MCA and the left PCA. Centering on the Voronoi point, the vertical bisector of the two generating points of the left ACA and the left MCA, the vertical bisector of the two generating points of the left MCA and the left PCA, and the two dividing points of the left ACA and the left PCA The vertical bisector of the generating point is connected. These vertical bisectors divide the left brain area into three areas of influence. Similarly, the right brain area has the right ACA, the right MCA, and the right PCA as three generating points, and the right A
It is divided into a power area of CA, a power area of right MCA, and a power area of right PCA.

Transfer function MTF of the time-density curve Ci (t) of the brain tissue with respect to the time-density curve Ca (t) of the left ACA
Is limited to the area of influence of the left ACA and is calculated for each pixel. Similarly, left MCA, left PCA, right ACA, right MCA, right P
The transfer function MTF is calculated for each pixel by confining it to each influence area with respect to CA.

As described above, by dividing the left brain area and the right brain area into a plurality of influence areas and confining the analysis range to each influence area, the number of analysis processing steps can be further reduced.

Next, on the CT image, the cerebral artery RO is placed on the cerebral artery.
I is set, but in order to improve and facilitate the accuracy of this setting, a support map is created by the ROI setting support unit 121 and displayed separately from the CT image or overlaid on the CT image (S4). ). As a support map, see Figure 9.
As shown in FIGS. 9A to 9C, for example, an AT (appearance time) map, a PT (peak time) map, and a TT (transit time) map can be cited.
For each pixel, as shown in FIG. 10, a time AT from an arbitrary time (for example, a data collection start time) before contrast to a time at which the concentration of the contrast agent reaches a few percent (for example, 1 percent) of the peak peak, a time T0. To peak time of contrast agent concentration (peak time)
PT, or TT that represents the moving time of the contrast agent, for example, in half width
Is calculated and generated and displayed as a map. By default, all types of these AT map, PT map, and TT map are generated and displayed, but the operator can select any one type or any two types.

Since these numerical values tend to appear at higher values in the cerebral arteries than in other tissues, color display is performed through a color lookup table that is set to display only pixels having values centered on those values. By doing so, it becomes possible to easily identify the location of the cerebral artery and accurately set the cerebral artery ROI (S5). Typically, in each of the left and right cerebral areas, the cerebral artery ROI is the anterior cerebral artery (AC
A), middle cerebral artery (MCA), and posterior cerebral artery (PCA) are set at three points each.

Note that in the case of multi-slice imaging, for example, when four adjacent slices are to be analyzed, FIG.
As shown in Fig. 1, the cerebral artery ROI is individually measured in each slice.
Setting is just a heavy work load, and is unnecessary work for analysis. Therefore, any one
The cerebral artery ROI set in the slice is also used in other slices. Alternatively, the time-density curve Ca (t) of the cerebral artery that can be commonly used for all slices may be created by using the coherent regression method described later.

Next, the time-density curve Ca (t) for each of the set cerebral artery ROIs is created by the time-density curve creating unit 122 from the continuous image data by dynamic CT (S6).

Here, many cerebral arteries are very thin as compared with the pixel size, and since they are generally not orthogonal to the CT imaging slice, every pixel on the image accurately represents the CT value of arterial blood. Not shown, one pixel is composed of a mixture of cerebral arteries and other tissues, and due to the partial volume effect, it shows almost no contrast effect lower than that. Further, image noise is large in these pixels including the artery as a partial volume. In particular, in the artery or the like where cerebral infarction occurs, the contrast effect is relatively small, so that the influence of noise is great. Although the image noise is suppressed by the above-mentioned coherent filter, the influence of the partial volume effect still remains.

This problem is suppressed by applying the coherent regression method, which will be described later, by using pixels in a solid including the artery, rather than measuring the time-density curve of the cerebral artery, in a single slice image. It is possible to Therefore, instead of the coherent filter method described above, the coherent regression method may be applied at this stage.

Further, according to this method, the time-density curve of the only slice image corresponding to each artery can be obtained,
Therefore, it can be used for analysis of any part in all slices within the imaging range, and by this, for a specific artery, the slice in which the time-concentration curve of the cerebral artery is obtained most clearly is selected. The cerebral artery time-density curve can be applied to all slices, and the number of cerebral artery time-density curves can be reduced.

(Coherent Regression Method) In creating the time density curve, it is important to remove the effects of partial volume effect and random noise. First, the “time density curve” is a curve that represents the change over time of the CT value (density value) at a specific site in the dynamic CT image. In particular, in the medical image diagnostic apparatus, for the purpose of investigating the details of blood flow dynamics and metabolic functions in human tissues, etc., a time-dependent concentration curve of changes over time in contrast agent concentration in specific tissues of the human body It is being measured. Further, in astronomical observation and the like, a time-density curve is used for the purpose of analyzing a change in luminous intensity of a specific celestial body. More formally, that is, the time concentration curve is the concentration value of a certain part at time t _k d _k
, The pair of columns {<t _k , d _k > (k = 1, 2, ...
.., K)}. Also, in many applications of time-density curves, the absolute value of d _k is not necessarily required, but rather only the increment (d _k −d ₁ ) relative to the first image 1 is needed. is there. Moreover many of such applications, (here A unknown constant) simply _(d k -d ₁₎ proportional to the data A _(d k _-d 1) is sufficient as long obtained only. In this case, therefore, the sequence of pairs {<t _k , A (d _k −d ₁ )> (k = 1,
2, ..., K)} is the time concentration curve to be obtained.

In order to obtain such a time density curve, in principle, the time density curve in each image k (k = 1, 2, ..., K) forming the above dynamic CT image is calculated. Using the scalar value v _k (x) of the pixel x included in the region to be measured, a pair of columns {<< t _k , v
_k (x)>} or {<t _k , A (v _k (x) -v
₁ (x) >>}.

However, in practical use, since the dynamic CT image taken by the medical image diagnostic apparatus or the like contains random noise, the time-density curve originally intended to be measured cannot be accurately obtained. There is.

Furthermore, in practical use, a so-called "partial volume effect" occurs in these dynamic CT images. What is the partial volume effect, that is, an image of a minute object in the subject is represented by a small number of pixels on the image, and these few pixels also include an image of an adjacent object in the subject. This is a phenomenon in which the pixel values of these small number of pixels show only a relatively small fluctuation (although proportional to the fluctuation of the density value to be measured) because of the influence. In other words, the pixel values of these few pixels contain very little signal.
Therefore, when the partial volume effect occurs, no matter which pixel x is taken, a pair of columns {<t _k , v _k (x)
> (K = 1, 2, ..., K)} has a very low signal level and is affected by the change in the concentration value in the tissue that is not originally intended to be measured, and further random noise exists. In addition, there is a problem that the time concentration curve {<t _k , d _k >} that is originally intended to be measured cannot be accurately obtained.

Therefore, in the past, in order to suppress random noise, temporal or spatial smoothing was used.
When time averaging is performed, the temporal resolution is impaired, and when spatial averaging is performed, there is a problem in that the change over time in the concentration of a portion other than the original 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 adopted.

First, the coherent filter of this embodiment
Describe the null hypothesis that should be used in. measurement
The true time-density curve at the part to be
t_k, D _k> (K = 1, 2, ..., K)}
The primary conversion is {<< t_k, A (d_k−
d₁)> (K = 1,2, ..., K)} (Dawashi A is not
In the case of measuring the knowledge coefficient),
Set up a set R of pixels that roughly corresponds to the part to be measured.
Set. For any pixel x ∈ R that is an element of this set R
Then, the condition Q: "If this pixel x is the true time density described above.
The curve is well reflected, and the change in concentration over time of other parts
If it is almost unaffected ", the vector value
Pixel value as v (x) = (v₁(X), v
_Two(X) ,. ． ． , V_K(X)), partial
Consider the effects of volume effects and random noise
By v_k(X) = p (x) d_k+ Q (x) + γ_k(X) (11) Assume that (k = 1, 2, ..., K) holds
You can Where p (x) and q (x) are the pixel x
Image number k (that is, shooting time t_kPhase
This is an unknown coefficient that does not change depending on
It is a model of the volume effect. Also γ_k
(X) is a model of random noise
The value is different for each pixel x and also for each image number k.
However, its expected value is 0, and its statistical distribution is pixel x
Neither does it depend on the image number k.

According to the above assumption, regarding any two pixels x and y which are elements of the set R, if “pixel x, y
Both satisfy the above condition Q. If the proposition “” holds, we can prove that the following relation holds.

V _k (x) = a ₁ v _k (y) + a ₂ + ξ _k (k = 1, 2, ..., K) (12) Here, a ₁ and a ₂ are a set of pixels. It is an unknown coefficient that differs for each x and y but does not change depending on the image number k (that is, the shooting time t _k ). Further, ξ _k is random noise, and for each pixel set x, y, the image number k
The value is different for each, but the expected value is 0.

Equation (12) is derived as follows.
That is, when the expression v _k (y) = p (y) d _k + q (y) + γ _k (y) ... (13) obtained by substituting y for x is transformed,

[Equation 9] Equation (12) is derived by setting here,
The a ₁ and a ₂ in the equation (16) are parameters expressing the partial volume effect, and ξ _k in the equation (16) represents random noise.

From the above, "Pixels x and y both satisfy the condition Q.
Add Is the null hypothesis H ₀´ “v_k(X) =
a₁v_k(Y) + a_Two+ Ξ_k  (K = 1, ..., K)
is there. It was shown to be equivalent to.

Next, the null hypothesis H ₀ '"v _k (x) = a ₁ v
_k (y) + a ₂ + ξ _k (k = 1, ..., K). Is converted into a proposition in a form that is substantially equivalent and can be actually tested. This null hypothesis can be described in a mathematically strict expression again.
₀ ′ 'There are certain constants a ₁ and a ₂ and ξ _k = v _k
(X) −a ₁ v _k (y) −a ₂ (k = 1, ..., K) follows a normal distribution with mean 0 and variance (σ ² h (a ₁ )). It will be. Here, the coefficient h (a ₁ ) is h (a ₁ ) = 1 + a ₁ ² (17). Equation (17) is the definition of a ₁ and ξ _k (16)
It is immediately derived from the equation and the general nature of the variance for random variables. Further, the value of the variance σ ² can be estimated easily and practically and sufficiently accurately.

From the above, if the above constants a ₁ and a
_{If 2} can be determined, it is possible to test the above null hypothesis H ₀ ′. Then, in practice, it is sufficient ^to obtain the optimum estimated values a ₁ ^to and a ₂ ^to these constants.

For the calculation of the optimum estimated values of the constants a ₁ and a ₂ as described above, a known fitting method can be used as it is. Therefore, in the following, as a typical specific example of such a fitting method, an outline in the case of using the linear least squares method will be described. In order to apply the linear least squares method to this embodiment, simply sum the square of ξ _k of the above null hypothesis as S (a), that is,

[Equation 10] Is defined. The value of S (a) is a constant vector a = (a ₁ ,
a ₂ ), ie the values of the constants a ₁ and a ₂ above. Be calculated constant vector a as the S (a) takes a minimum value, to the constant a ₁ and a _2, the optimum estimated values a ₁ ^~ and a ₂ ^~ of the sense of unbiased estimate can be obtained. As a concrete calculation method of the linear least squares method, various known methods can be used, and all of these known calculation methods are very simple, and the required calculation time is very short. Is.

[0134] Thus, the optimum estimated values _a ¹ ~ of the constants ^{_{a 1, a 2, a 2}} ~ result of calculating the residual _r ^k ~ defined by the following equation (x, y) = v _k a ^{_{(x) -a 1 ~ v k}} (y) -a 2 ~ ... (19) may be specifically calculated. Therefore, this residual r
_{Using k} ^~ , the null hypothesis H ₀ ′ described above is converted into a substantially equivalent null hypothesis H ₀ ″ “r _k ^~ (x, y) (k = 1, ...,
K) mean 0, follows a normal distribution of variance ^{_{(1+ (a 1 ~) 2}} ) σ 2. Can be paraphrased. This is a concrete proposition that can actually perform the calculation of the test.

Furthermore, a vector representation

[Equation 11] (However, the vectors a and ξ depend on the set of pixels x and y.) And the following equation f (a ^~ , v (y)) = a ₁ ^~ v (y) + a ₂ ^~ ... ( using vector value function f defined in 21) the null hypothesis H ₀ '
In other words, the null hypothesis H ₀ ″ is “v (x) = f (a
^~ , V (y)) + ξ, where ξ is mean 0 and variance (1
+ (A ₁ ^~ ) ² ) σ ² follows a normal distribution. ) ”,
This is exactly the same form as the null hypothesis H ₀ described above. That is, it is obvious that the present embodiment is a modification of the above-mentioned coherent filter. It should be noted that here, f (a ¹ , v (y)) means that the parameter a representing the partial volume effect is optimally adjusted with respect to the pixel value v (y) of the pixel y, and the pixel x Pixel value v
(X) means the one converted so as to have the highest conformity.

Next, in the present embodiment, a method for obtaining a time density curve by a coherent filter using the above null hypothesis H ₀ ″ will be described. Regarding a set R of pixels that roughly corresponds to a region to be measured , For one pixel xεR included in this set R, the following calculation is performed for all pixels yεR that are elements of the set R. That is, the residual error is actually obtained using the above method. r _k
^~ (X, y) (k = 1, ..., K) are calculated, and then the above null hypothesis H ₀ ″ “r _k ^~ (x, y) (k = 1, ..., K).
K) mean 0, follows a normal distribution of variance ^{_{(1+ (a 1 ~) 2}} ) σ 2. Specifically, the risk rate p (x, y) or the weight w (p (x, y)) in the case of 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 at the pixel x is calculated.
(X) -v ' ₁ (x)> (k = 1, 2, ..., K)}
Make up.

[0137]

[Equation 12]

The time-density curve thus obtained is the pixel x
Of the true time-density curve {<t _k , d _k >}, which is a linear transformation of {<t _k , A (d _k −d ₁ )>} (where A is an unknown coefficient), is approximated. It is a measured value, and random noise is suppressed by the effect of weighted averaging. Further, for the pixel value vectors of the other pixels y, the ones in which the influence of the partial volume effect is corrected are used, as is clear from the equation. Further, the present embodiment has the property that "a temporal average is not used at all, and a spatial average is calculated using a weight based on the degree of matching with the pixel x", which is a common feature of coherent filters. 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 without impairing the time resolution and random noise is suppressed. The method of obtaining the time-density curve in this way is particularly called a "coherent regression method".

Next, specifically, a dynamic CT obtained by dynamic CT imaging or the like in medical X-ray CT
An example of clinical use of the time-density curve in an image will be described. In this application example, while a contrast agent is rapidly injected into a blood vessel, imaging such as dynamic CT is performed, and a change in the concentration of an image of an artery existing in a human body tissue is measured as a time-density curve to obtain a blood concentration in the tissue. It is intended to diagnose the flow dynamics.

In this application example, in many cases, the human body
Since arteries in the weave are generally very thin, CT
The image of the artery appearing in the image is a partial volume effect.
Cause In addition, the image contains random noise
Needless to say For this reason, the traditional method
To obtain a sufficiently accurate time-density curve for the artery
Is difficult, and if you make a strong measurement, the true
Time concentration curve <t_k, D _k<T which is a linear transformation of_k，
A (D_k-D₁)> (D here)_kCorresponds to the image of the artery
Time t of a group of pixels_kAt (which is a scalar value)
Represents a pixel value. Also, there are k = 1, 2, ..., K)
Measured value approximating to some extent <t_k, (V_k(X) -v
₁Only (x))> was obtained. This measurement is random
Includes unwanted noise. Also, of the partial volume effect
Due to the effect, the coefficient A remains unknown.

[0141] _{Therefore, <t k, A (D} k -D 1)> measurements to sufficiently approximate the _{<t k, (v'k (} x) -v'
₁ (x) >> (k = 1, 2, ..., K) can be obtained. On the other hand, in the vein can be observed on the same tomography image, there are things considerably thicker, therefore For those veins, in a conventional manner, sufficiently good approximation of the time-density curve <t _{k, (J k -J 1)> (k =} 1,2, ···,
K) can be obtained. Here, J _k represents the pixel value of a group of pixels corresponding to a vein image at time t _k .

By the way, in the time-concentration curve concerning blood circulation, the proposition S: "If the contrast agent concentration in blood at time t ₁ is 0, the time-concentration curve <t _k , (d _k −d ₁ )>, the area under the curve (AUC: Area Under Curve) matches ”. The area under the curve referred to here, the time-density curve _{<t k, (d k -d} 1)> means the integration with respect to time t of.

Therefore, the area under the curve AUC (d) of the time-density curve <t _k , (d _k -d ₁ )> for a certain blood vessel d
Can be approximately calculated by the following equation.

[0144]

[Equation 13]

Therefore, the area under the curve AUC (J) for the time-concentration curve {<t _k , (J _k -J ₁ )>} obtained by the conventional method for veins should be calculated using the equation (22). You can (J may be substituted for d.) As for the artery, the time-density curve {<t _k , (D _k −
If D ₁ )>} is known, the area under the curve AUC (D) can be calculated in the same manner using the equation (18), and according to the above proposition S, AUC (D) ≈AUC (J) (24) should hold. However, in practice, since the time-density curve _{<t k, (D k -D} 1)> is unknown, AUC
(D) cannot be calculated.

[0146] On the other hand, the time-density curve obtained in a manner according to this embodiment _{<t k, (v'k (} x) -v' 1 (x))>
Is intended to approximate the _{<t k, A (D k} -D 1)>,
The latter contains an unknown coefficient A. Therefore, {<<
t _k, from _{(v'k (x) -v' 1} (x))>} (2
Area under the curve AUC that can be specifically calculated using the equation 3)
(V ') must be exactly A times AUC (D). That is, AUC (v ′) ≈A · AUC (D) (25). That is, the relationship of A≈AUC (v ′) / AUC (J) ... (26) is established from the expressions (24) and (25). Since the right side of equation (26) can be specifically calculated using equation (23), the unknown coefficient A
The value of can be specifically determined. Therefore, the time-density curve _{<t k} by using the value of the coefficient _{A, (v'k (x)} -v'
By configuring the _{1 (x)) / A>} , which is time-density curve of the artery _<t k, nothing but to approximate the (D _k -D _1)>. As described above, a method of forming a time-concentration curve in which the unknown value of the constant A is determined by using the area under the curve is referred to as “AUC method”.

From the above, in clinical use of the time-density curve in a dynamic CT image obtained by dynamic CT imaging or the like, by combining the coherent regression method with the AUC method,
Even with respect to the time-density curve of a thin artery, which is difficult or impossible to measure by the conventional method, the effect of partial volume effect and random noise can be eliminated, and a measured value that does not include the unknown constant A can be obtained. .

Of course, the AUC method alone is the time-concentration curve <t _k ,
(V ′ _k (x) −v ′ ₁ (x))> is also applicable, and the unknown value of the constant A is determined (although the effects of random noise and the partial volume effect cannot be excluded). A time-density curve can be constructed.

(Correction of time-density curve of cerebral artery using time-density curve of superior sagittal sinus (suppression of influence of partial volume effect)) In order to suppress the influence of partial volume, this coherent regression is used instead. Time curve of upper sagittal sinus Csss
The time concentration curve Ca (t) of the cerebral artery may be corrected using (t).

First, in step S7 of FIG.
As shown in 2, a large upper sagittal sinus ROI is set so as to surround the upper sagittal sinus on the CT image. Since the superior sagittal sinus is larger than the artery and its position is relatively fixed, it is easy to set the superior sagittal sinus ROI. The larger sagittal sinus ROI contains a plurality of pixels.

Next, all the pixels in the superior sagittal sinus ROI are included in the superior sagittal sinus sinus over the entire area.
The upper sagittal sinus ROI is reduced (S8). As the reduction processing, for example, first, threshold processing (binarization) is executed for each pixel in the upper sagittal sinus sinus ROI, and RO
A binary map (“0”, “1”) in I is created. The threshold value is set to a value that separates the image of the upper sagittal sinus from the images of tissues and bones around it. “1” represents a pixel on the image of the upper sagittal sinus, and “0” represents a pixel on the image of the 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 of the neighboring 4 or 8 pixels are “1”,
The value of the central pixel is maintained at "1". That is, not only when the central pixel is “0”, but even when the central pixel is “1”, even if one of the neighboring 4 or 8 pixels indicates “0”, the value of the central pixel is Replace with "0". Therefore, the superior sagittal sinus ROI is reduced by at least one pixel from the outline of the image of the superior sagittal sinus. Thereby, the condition that all the pixels in the superior sagittal sinus sinus ROI subjected to the reduction process are the pixels on the superior sagittal sinus sinus image can be realized with high accuracy.

Further, instead of this method, the upper sagittal sinus sinus ROI may be corrected by using the area AUC under the curve of the time concentration curve. In this case, the large ROI is set as the search range, and the area AUC under the curve of the time density curve is calculated for each pixel therein. Due to the contrast effect, the area AUC under the curve of the pixel on the upper sagittal sinus image clearly shows a higher value than that of the peripheral pixels. Therefore, by performing threshold processing on this area AUC under the curve, R
Only pixels on the upper sagittal sinus image can be selected from the OI.

In this way, with respect to a plurality of pixels having a high certainty that they are on the superior sagittal sinus sinus image picked up by the AND condition by using either one or both methods, the time density curve of each pixel is obtained. Are averaged and a time-density curve Csss (t) of the superior sagittal sinus is created (S9).

Here, the iodine contrast agent is the blood-brain barrier.
Since it does not pass through the brain barrier, the iodine concentration does not change in principle between the cerebral artery and the cerebral vein, that is, the area under the curve AU of the time concentration curve Csss (t) of the upper sagittal sinus.
C is the time-density curve Ca (t) of the cerebral artery created in S6
Is substantially equivalent to the area AUC under the curve. Therefore, as shown in FIG. 13, the time-density curve Csss of the superior sagittal sinus
For the area under the curve AUCsss of (t), the area under the curve AU of the time concentration curve Ca (t) of the cerebral artery created in S6
AUC (sss / is added to each time value of the time concentration curve Ca (t) of the cerebral artery created in S6 so that Ca is substantially equivalent.
It is corrected by multiplying (AUCa) (S10).

Next, using the time-concentration curve Ca (t) of the cerebral artery shown in FIG. 14 (a) in which the noise and the partial volume effect have been suppressed as described above, the blood flow dynamics of the brain tissue (capillaries) is used. To quantify the situation. For that, first
A time-density curve Ci (t) shown in FIG. 14B is created for each pixel on the brain tissue (S11).

Next, as shown in S12, the cerebral artery time-density curve Ca (t) is used for each pixel by using the cerebral artery time-density curve Ca (t) separately in the left and right areas. Using the time-density curve Ci (t) as an output function, the characteristics of the process in which the tracer passes through the capillaries are calculated using the transfer function MT
Calculate as F. That is, the same time concentration curve Ca (t) of the cerebral artery in the left area is used for the time concentration curve Ci (t) of the brain tissue in the left area, and the time concentration curve Ci (t) of the brain tissue in the right area is used. For, the transfer function MTF is obtained using the time concentration curve Ca (t) of the same cerebral artery in the right area. Further, as described above, the cerebral artery time concentration curve Ca
Since (t) is created for each of ACA, MCA, and PCA, the calculation of the transfer function MTF is repeated for each Ca (t).

Here, as shown in FIG. 15, the time concentration curve Ca (t) of the cerebral artery and the time concentration curve C of the capillaries are shown.
The box-MTF method is applied using an ideal relationship that holds with i (t) as an analytical model.

FIG. 16 shows the principle of the box-MTF method. The residual difference between the time-concentration curve Ca (t) of the cerebral artery and the convolution Ci ′ (t) of the transfer function box-MTF represented by a rectangular function and the measured Ci (t) created in S11 is evaluated. , The transfer function box-MTF is modified so as to reduce the sum of squares of the residuals. This procedure (procedu
The residual is minimized by repeating the re) routine.

Transfer function box-MTF that minimizes residuals
Based on CBP, CBV, M
TT is calculated (S13), and the residual sum of squares minimized in S12 is output as Err. Strictly speaking, CBP = CBP CBV = (1-Ht) / (1-b * Ht) * CBV 'MTT = (1-Ht) / (1-b * Ht) * MTT' is corrected. Here, Ht is a hematocrit value of a large blood vessel, and b * Ht is a hematocrit value of a peripheral blood vessel (generally, b is about 0.7).

The residual is given by y _i (x) -f (t _i , x). y _i (x) represents the scalar value of voxel x at time t _i , and corresponds to the time-density curve of brain tissue. f
(t _i , x) represents the scalar value at the time t _i of the model fitted to the vector pixel value of the voxel x, and corresponds to the convolution of the transfer function and the time-density curve of the cerebral artery. Err is the square root of the sum of squares of the residuals when the transfer function is approximated, and is calculated as shown in the following expression, for example.

[0161]

[Equation 14]

Here, S is a constant, for example, S = N-
p. p is the degree of freedom, that is, the model f to be approximated
Represents the number of parameters included in. w (t_i) Is the time
Tick t_iDetermines how much the residuals at contribute to Err
It is a weighting factor. For example, w (t _i) Does not depend on i
Then, it may be a fixed value such as 1 w (t_i) = Αe^{-ti Two / β} As | t_iThe larger the |, the smaller the weight w becomes.
So that the residual error becomes Err over time.
May be less likely to contribute to.

Thus, CBP, CBV, MTT, Er calculated from the transfer function obtained by the box-MTF method.
As shown in FIG. 17, an output range (appropriate range) is individually set for r (S14). Pixels having values within the corresponding output range are maintained, and pixels having values outside the output range are replaced with that value, for example, a value corresponding to black on the display (S15). .

Next, as shown in FIGS. 18 (a) to 18 (d), CBP, CBV, MT subjected to output optimization.
Each map is generated from T and Err (S1
6). The indexes of CBP, CBV, MTT, and Err are calculated for each of the anterior cerebral artery ACA, the middle cerebral artery MCA, the posterior cerebral artery PCA, and individually for each slice. Therefore, even with one slice, the map is 4 × 3 = 12, as shown in FIG. If the number of cerebral arteries to be set is increased, the map will increase by 4 times the increased number. It is not realistic to evaluate such many maps comprehensively. Therefore, in order to reduce the number of maps, the maps are combined (S17).

As shown in FIG. 20, the synthesizing method is as follows.
CBP map of anterior cerebral artery ACA and middle cerebral artery MCA
CBP map and posterior cerebral artery PCA CBP map, the residual Err of the anterior cerebral artery ACA, the middle cerebral artery MCA
And the residual Err of the posterior cerebral artery PCA. For example, the time concentration curve C of the anterior cerebral artery ACA
a (t) and the time-density curve Ci of the brain tissue under its control
When the transfer function MTF is calculated from (t), its residual E
rr 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 becomes relatively large. That is, the residual Err
Indicates the controllability of each cerebral artery.

Therefore, for each pixel, the anterior cerebral artery ACA
The CBP value corresponding to the value with the lowest residual error Err is selected as the value of the pixel from the CBP value of, the CBP value of the middle cerebral artery MCA, and the CBP value of the posterior cerebral artery PCA. The map thus synthesized is the anterior cerebral artery ACA,
It is composed of CBP values of brain tissues that are likely to be under the control of middle cerebral artery MCA and posterior cerebral artery PCA. The same applies to map composition of other indexes CBV and MTT.

Map synthesis will be described in detail below. The time-concentration curve obtained from pixels at the position corresponding to the artery reflects the contrast agent concentration in arterial blood, and by applying the above-mentioned coherent regression method, etc. Can be obtained. Such a time-concentration curve of a cerebral artery can be created for each artery, and there are differences depending on the blood circulation state. This difference may be significant especially in cases where cerebrovascular accidents occur. For example, the time-density curve of the cerebral artery obtained from K arteries is represented by A _k (t) (k = 1, 2, ...
・, K).

By comparing the time-density curve of a tissue with the time-density curve of the cerebral artery of the artery feeding (dominant) the tissue, the microcirculation (capillary system It is possible to obtain an index such as CBP that reflects the structure and function). Since these indexes are calculated for each part (x, y, z), it is possible to form an image having the values as pixel values, and such an image is an index map. For example, when indexes of R types (typically, four types of CBP, CBV, MTT, and Err as described above) are obtained, R maps can be configured. The R maps thus created can be regarded as one map (vector value map) in which each pixel has 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)>.

For example, in a CBP study, typically, R = 4 as described above, P _{k, 1} (x, y) is the value of CBP, and P _{k, 2} (x, y, z) is the value of CBV. The value is set to P _{k, 3} (x, y,
z) can be configured to represent the value of MTT, and P _{k, 4} (x, y, z) can be configured to represent the value of residual Err.

Among the parts (x, y, z), those which are clearly not corresponding to the organ to be analyzed are excluded from the analysis from the beginning, and P _{k, r} (x, It is advisable to substitute a special value indicating that the analysis target is not included in y) (step S above).
14, S15). It is convenient to use a negative value having a large absolute value as such a value. Alternatively, the vector V
_As yet another element to be added to _k (x, y, z), if P _{k, R + 1} (x, y, z) = ((x, y, z) is not an analysis target, 0 Alternatively, you can make a map called 1). Such a map is called a "mask".

Such a vector value map V _k is created for each time concentration curve A _k of the referenced cerebral artery. For example, suppose that the temporal concentration curve of the cerebral artery is obtained from the left and right middle cerebral arteries, the anterior cerebral artery, and the posterior cerebral artery. , K = about 10 to 15.

Of these, the time-density curve of the cerebral artery obtained from the artery in the right hemisphere shows the region (x, y, z) belonging to the right hemisphere.
The time-concentration curve of the cerebral artery obtained from the artery in the left hemisphere should be used only for the analysis of the region (x, y, z) belonging to the left hemisphere. Therefore, it is desirable that the operator designates the boundary (median line) between the right hemisphere and the left hemisphere as a straight line, a curved line or a polygonal line, or a flat surface or a curved surface, and creates a map for each hemisphere. However, nevertheless, there may be as many as K = 3 to 10 temporal concentration curves of cerebral arteries for each hemisphere.

In this way, the time-density curve A _k of the cerebral artery
When the number K of (k = 1, 2, ..., K) is large, the resulting vector value map V _k (k = 1, 2,
(..., K) is large, which is inconvenient to observe. That is, to observe as a normal gray scale image or color scale image, one map is composed of R images, and since there are K sets, a total of K × R images must be compared. . Furthermore, it is not always obvious which part is fed by which artery, and which map V _k (k = 1, 2, ..., K) is observed for each part using anatomical knowledge. You have to decide what to do. Particularly in cases where cerebrovascular disorders such as cerebral infarction occur, which arteries control the tissue does not always correspond to anatomical knowledge, and abnormal control is often seen. Due to these problems, it is difficult to interpret the vector value map.

In order to solve this problem, map synthesis 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, P _{k, R} (x, y, z)
Is a residual map, V (x, y, z) = V _k (x,
y, z) where k is | P among k = 1, 2, ..., K
_{Let k, R} (x, y, z) | be the minimum k.

Further, k = 1, 2, ...
, Map P ₀ (x, y, z) = (k = 1,2, ..., K for indicating which of K has been adopted | P of K
It is also possible to add _{k, R} (x, y, z) | such that k) is the smallest.

According to this method, that is, when observing as a normal gray scale image or a color scale image, R to R + 1 images may be observed.

According to this method, there is a possibility that the calculation result using A _k may be erroneously used at the site (x, y, z) which should be calculated using the time-density curve A _j of the cerebral artery. . However, if such an error occurs, V (x,
As is clear from the definition of y, z), | P _{k, R} (x, y, z)
It is necessary that | <| P _{j, R} (x, y, z) |, and such a relationship occurs only when A _j and A _k are extremely similar. Therefore, at the part (x, y, z), V
_{It is considered that k} (x, y, z) and V _j (x, y, z) are originally similar to each other, and it is almost possible that this error causes an error in the interpretation of V (x, y, z). Absent.

When this method is actually applied, A _j and A _k
P ₀ (x, y, z) = k for each part in the tissue that is considered to be almost uniform only when the values are very similar.
Or P ₀ (x, y, z) = j occurs, and then P _{k, r} (x, y, z) ≈P _{j, r} (x, y, z) (r =
It is observed that the results are almost the same regardless of which method is adopted.

On the contrary, when the tissue is controlled by a specific artery and the corresponding time-density curve A _k of the cerebral artery is not similar to other curves, the present method is used. At the site (x, y, z) in the tissue, V _k (x, y, z) is almost certainly and automatically selected. Therefore, by observing P ₀ (x, y, z), it is possible to observe which tissue is controlled by which artery without anatomical knowledge.

Now, let us return to FIG. As shown in FIGS. 21 (a) to 21 (d), a CBP map, a CBV map, an MT synthesized as described above or a single cerebral artery
A region of interest ROI including a plurality of pixels is set for the T map and Err map (S18), and the average value (CBP average value, CBP value, CBP value, CBV value, MTT value, Err value) of the pixel values in the ROI is set. CBV average value, MTT average value, Err
An average value is calculated (S19), and the average value may be used as a diagnostic material. In this averaging, the above step S
In 14, the proper range is set for each of CBP, CBV, MTT, and Err, the value within the range is maintained, and the value outside the range is replaced with the minimum value corresponding to the black expression. When the average including the replaced values is included, the average value naturally includes an error. Therefore, in this averaging process, it is necessary to select only values within the appropriate range or exclude replaced values and perform the averaging process.

The region of interest ROI for this averaging
When setting, CBP map, CBV map, M
If the region of interest ROI is set on any one of the TT map and the Err map, the ROI can be commonly used for other maps as well, which simplifies the ROI setting work. , Average values for the same ROI (CBP average value, CBV average value, MTT average value, Err
It is possible to calculate the average value).

Here, as described above, the minimum residual Err of the time-density curve of a certain tissue with respect to the time-density curve of a certain cerebral artery is the extent to which the cerebral artery controls the tissue, that is, the brain. It represents the extent to which the arteries are responsible for the blood supply to the tissue. In other words, it indicates how dependent the tissue is from the cerebral artery, that is, how much the cerebral tissue is supplied with blood flow from the cerebral artery. The cerebral artery corresponding to the small residual Err indicates that the pixel has a high dominance to the brain tissue, and the cerebral artery corresponding to the large residual Err indicates that the pixel has a low dominance to the brain tissue. It represents. Therefore, a map in which a cerebral artery having a high controllability for each pixel is distinguished by a label from the residual Err, that is, a region having a high controllability of the anterior cerebral artery ACA and a region having a high controllability of the middle cerebral artery MCA, It is possible to generate a control map that distinguishes the region of the posterior cerebral artery PCA that is highly likely to be controlled.

As shown in FIG. 22, the residual Err of the anterior cerebral artery ACA, the residual Err of the middle cerebral artery MCA, and the residual Err of the posterior cerebral artery PCA are shown for each pixel for each of the left and right areas of the brain. Compare. It indicates that the cerebral artery (ACA, MCA or PCA) showing the smallest residual error is most likely to control the brain tissue of the pixel.
For each pixel, it is most likely to dominate, that is, the residual Err
The cerebral artery that has the smallest value is identified. A label corresponding to the specified cerebral artery is given to each pixel.

FIG. 23 shows an example of the generated control map. This control map is displayed by distinguishing the labels by color and shade. In addition, by filtering the index map with an arbitrary label, the map shown in FIG.
As shown in (b) and FIG. 24 (c), it is possible to extract from the index map a region having a high controllability for each cerebral artery (ACA, MCA, or PCA).

As described above, according to the present embodiment, the coherent filter or the coherent regression can suppress the deterioration of the spatial and temporal resolution and suppress the noise, thereby improving the analysis accuracy of the CBP study. .

(Modification) The present invention is not limited to the above-described embodiment, and can be variously modified and implemented at the stage of implementation without departing from the scope of the invention.
Further, the above embodiment includes various stages,
Various inventions can be extracted by appropriately combining a plurality of disclosed constituent elements. For example, some constituent elements may be deleted from all the constituent elements shown in the embodiment.

[0187]

According to the present invention, the diagnostic efficiency of the map obtained by the CBP study can be improved.

[Brief description of drawings]

FIG. 1 is an explanatory diagram of the principle of a CBP study.

FIG. 2 is a block diagram showing the configuration of an index calculation device for blood flow dynamics of capillaries in brain tissue according to the embodiment of the present invention.

FIG. 3 is an explanatory diagram of image processing by the coherent filter of the present embodiment.

FIG. 4 is a flowchart showing the flow of noise suppression processing by the coherent filter in the present embodiment.

FIG. 5 is a flowchart of the first half of the entire index calculation process according to this embodiment.

FIG. 6 is a flowchart of the latter half of the entire index calculation process according to this embodiment.

FIG. 7 is a diagram showing an example of a dividing line in step S3 of FIG.

FIG. 8 is a diagram showing the influence area of each cerebral artery on the Voronoi diagram used for reducing the processing man-hour of step S12 of FIG.

FIG. 9 is a diagram showing an AT map, a PT map, and a TT map in step S4 of FIG.

10 is a diagram showing AT, PT, and TT in step S4 of FIG.

FIG. 11 is a diagram showing a cerebral artery ROI common between slices in step S6 of FIG.

FIG. 12 is a diagram showing the superior sagittal sinus ROI set in step S7 of FIG.

FIG. 13 is a supplementary diagram regarding correction of a time-density curve of a cerebral artery in step S10 of FIG.

14 is a diagram showing an example of a time-density curve Ca (t) of a cerebral artery and a time-density curve Ci (t) of a brain tissue created in steps S10 and S11 of FIG.

15 is an explanatory view of the principle of the box-MTF method of step S12 of FIG.

16 is an explanatory diagram of a box-MTF process of step S12 of FIG.

17 is a diagram showing an example of an output range setting screen for each index in step S14 of FIG.

FIG. 18 is a CBP created in step S16 of FIG.
The figure which shows an example of each map of CBV, MTT, and Err.

FIG. 19 is a CBP map, CBV map, MTT map, Er created for each cerebral artery in step S16 of FIG. 6;
The figure which displayed r map by the list.

FIG. 20 is a diagram for explaining the map synthesizing method of step S17 of FIG. 6;

FIG. 21 is a view showing a display example of the average value calculated in step S19 of FIG.

FIG. 22 is a diagram showing a cerebral artery (ACA, MC) that is highly likely to be controlled for each pixel (local tissue) in step S17 of FIG. 6;
The figure which shows the generation method of the map (control map) which distinguished A, PCA) by the label.

FIG. 23 is a diagram showing an example of a control map generated by the generating method of FIG. 22.

FIG. 24 shows a CBP map filtered using a dominance map.

[Explanation of symbols]

10 ... Gantry section, 20 ... Computer device, 30 ... Image processing device, 101 ... X-ray tube, 101a ... High voltage generator, 102 ... X-ray detector, 103 ... a data collection unit, 104 ... Preprocessing unit, 105 ... memory part, 106 ... an image reconstruction unit, 10M ... storage device, 107 ... image display unit, 108 ... control unit, 109 ... Input section, 110 ... Coherent filter processing unit, 111 ... Variance value estimation unit, 112 ... Weighting function calculation unit, 113 ... Pixel value calculation unit (coherent filter unit), 120 ... CBP study processing unit, 121 ... ROI setting support unit, 122 ... Time concentration curve creation unit, 123 ... Cerebral artery time-density curve correction unit, 124 ... MTF processing unit, 125 ... Index calculator, 126 ... Map creation section, 127 ... Map synthesis unit.

   ─────────────────────────────────────────────────── ─── Continued front page    F-term (reference) 4C093 AA22 AA24 BA09 DA02 DA04                       FD03 FD05 FD12 FF28 FG01                       FG14 FG16

Claims (23)

[Claims] 1. A first time-density curve for an artery in the specific site and a second time-density for a tissue in the specific site from a plurality of consecutive images of a specific site of a subject injected with a contrast agent. A curve is created so that the transfer function representing the local hemodynamics of the tissue for each of the arteries is minimized by the residual of the second time concentration curve with respect to the convolution of the transfer function and the first time concentration curve. Is calculated by curve fitting, the index relating to the local hemodynamics for each of the arteries is calculated from the transfer function, the map of the index corresponding to the artery is created, and the map of the index is calculated at the first time. An index calculation method characterized by combining into one map according to the residual for each concentration curve. 2. The index calculation method according to claim 1, wherein the specific part is a brain. 3. As the index, at least one of a blood volume per unit volume and time of brain tissue, a blood volume per unit volume in brain tissue, and an average transit time of blood flow through brain tissue is calculated. The index calculation method according to claim 2, which is characterized in that. 4. The transfer function is calculated separately for a left area and a right area of the brain.
The index calculation method described. 5. The index calculation method according to claim 2, wherein the transfer function is separately calculated in a power range of the artery obtained by dividing the brain region by a geometric method. 6. The index calculation method according to claim 5, wherein the geometric method is a Voronoi diagram. 7. The index calculation method according to claim 2, wherein the first time-density curve for the artery is corrected based on the time-density curve for the vein created from the image. 8. The first time concentration curve is corrected so that the area under the curve of the time concentration curve of the vein is substantially equivalent to the area under the curve of the time concentration curve of the first time concentration curve. The index calculation method according to claim 7, wherein: 9. The index calculation method according to claim 7, wherein the time density curve of the vein is created from an average density of a plurality of pixels in the area of the vein. 10. The characteristic value of a time-density curve of each of a plurality of pixels forming the image is calculated to support the setting of the artery, and a characteristic value map is displayed. Index calculation method. 11. The characteristic value is an appearance time,
11. The index calculation method according to claim 10, which is a peak time or a transit time. 12. The value of the index is replaced with a predetermined value corresponding to black, its approximate color or a color not used for image interpretation when the value of the index is out of a predetermined range. The index calculation method described. 13. A first time-density curve for an artery in the specific site and a second time-density for a tissue in the specific site from a plurality of consecutive images of a specific site of a subject injected with a contrast agent. By creating a curve and selecting one of the first time concentration curves that best fits each of the second time concentration curves,
Identifying one of the arteries that is most likely to depend on the local hemodynamics of each of the tissues, and creating a map that distinguishes the subordinate region of the artery based on the one of the arteries identified for each of the tissues. A method characterized by. 14. The method according to claim 13, wherein the map is displayed in such a manner that the dependent regions of the artery are color-coded. 15. The index map relating to the index representing the local hemodynamics corresponding to the artery is synthesized based on the map.
Method 3 16. Filtering the plurality of images using a weight according to the degree of similarity between pixels in order to reduce noise from a plurality of consecutive images of a specific part of a subject injected with a contrast agent. , A first time-density curve for an artery in the specific site and a second time-density curve for tissue in the specific site are created from the plurality of noise-reduced images, and local blood of the tissue for each of the arteries is generated. The transfer function representing the flow dynamics is calculated by curve fitting so that the residual of the second time concentration curve with respect to the convolution of the transfer function and the first time concentration curve is minimized, and from the transfer function, the artery is calculated. An index calculation method comprising calculating an index relating to the local hemodynamics for each. 17. The method according to claim 16, wherein the similarity between the pixels is determined based on a difference in a time density curve between the pixels. 18. The method according to claim 16, wherein the similarity between the pixels is determined based on a difference in a time density curve between the pixels and a variance expected value of the noise. 19. The method according to claim 16, wherein the pixel value of each of the pixels forming the image is replaced with a weighted average of the surrounding pixels and the weight. 20. The index calculation method according to claim 16, wherein the transfer function is a rectangular function. 21. A first time-density curve of an artery in the specific site and a second time-density of a tissue in the specific site from a plurality of consecutive images of a specific site of a subject injected with a contrast agent. A time-concentration curve creating unit for creating a curve, and a transfer function representing the local hemodynamics of the tissue for each of the arteries, A transfer function calculation unit that calculates by curve fitting so as to minimize the difference, an index calculation unit that calculates an index relating to the local hemodynamics for each of the arteries from the transfer function, and an index of the index corresponding to the artery. A map creation unit for creating a map, and a map of the index according to the residuals for each artery, Index calculation apparatus characterized by comprising a map synthesis section which synthesizes the flop. 22. A first time-density curve of an artery in the specific site and a second time-density of a tissue in the specific site from a plurality of consecutive images of a specific site of a subject injected with a contrast agent. By selecting a time-density curve creating unit that creates a curve and one of the first time-density curves that best fits each of the second time-density curves,
A unit that identifies one of the arteries that is most likely to depend on each of the tissues, and a map creation unit that creates a map that distinguishes the dependent regions of the arteries based on one of the arteries identified for each of the tissues; An apparatus comprising: 23. In order to reduce noise from a plurality of consecutive images of a specific part of a subject injected with a contrast agent, the plurality of images are filtered using a weight corresponding to a similarity between pixels. A filter section, and a time-density curve creating section that creates a first time-density curve for an artery in the specific site and a second time-density curve for tissue in the specific site from the plurality of noise-reduced images, A transfer function representing the local hemodynamics of the tissue for each of the arteries is calculated by curve fitting so as to minimize the residual of the second time concentration curve with respect to the convolution of the transfer function and the first time concentration curve. And a local hemodynamics for calculating the local hemodynamics for each of the arteries from the transfer function. And an index calculation unit that calculates an index for the index calculation device.