JP2009050726A - Method and apparatus for calculating index for local blood flow kinetics - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To improve diagnostic efficiency of a map obtained by the CBP study. <P>SOLUTION: First time-density curves on the arteries within a specific region and second time-density curves on the tissues within the specific region are created from a plurality of continuous images on the specific region of a subject to which a contrast agent is injected. By selecting one of the first time-density curves most fitted to each second time-density curve, one artery with the highest possibility of the local blood flow of each tissue being subordinate to is specified, and a map discriminating subordinate areas to the arteries is created based on the artery specified for each tissue. <P>COPYRIGHT: (C)2009,JPO&INPIT

Description

  The present invention relates to an index calculation method and calculation apparatus related to local blood flow dynamics such as brain tissue.

  In the X-ray CT examination, morphological information can be obtained as visual information from a simple CT image, and dynamic information of blood flow around a lesion can be obtained as a dynamic scan by contrast CT. In recent years, high-speed scanning by multi-slice CT has become possible, and it is considered that the range of utilization of dynamic scanning by contrast CT will be expanded.

  As one of the directions, there is a method called a CBP study for calculating an index related to blood flow dynamics of capillaries in brain tissue. In the CBP study, indices such as CBP, CBV, MTT, and Err that quantitatively represent the local blood flow dynamics in the tissue, that is, the blood flow dynamics through the capillaries in the local tissue are obtained. Output the map.

  CBP represents the unit volume in the capillary of the brain tissue and the blood flow volume per unit time [ml / 100 ml / min], CBP represents the blood volume per unit volume in the brain tissue [ml / 100 ml], and MTT represents The blood average transit time [seconds] of the capillary blood vessel is represented, and Err represents the sum of the residuals or the square root of the squared sum of the residuals when approximating the transfer function.

  Indexes CBP, CBV, and MTT that quantitatively represent the blood flow dynamics of the capillaries in these brain tissues, together with information on the elapsed time since the onset of cerebral ischemic stroke, are used for pathological differentiation of ischemic cerebrovascular disorders. It is expected as useful information for the evaluation of the presence or absence of enlargement of capillaries, blood flow velocity, and the like. For example, generally in an ischemic cerebrovascular disorder, the blood pressure of the provided artery is reduced, and the blood flow velocity in the blood vessel is reduced. As a result, even if CBV is constant, MTT is prolonged and CBP is lowered. In the cerebral infarction hyperacute phase, in order to compensate for the decrease in blood flow rate due to the decrease in blood pressure, the capillaries are expanded and the blood flow rate is increased to suppress the decrease in blood flow CBP (auto Regulation). Therefore, even if CBP decreases due to the extension of MTT, if CBV increases, it is information suggesting the possibility of reopening of capillaries.

  In the CBP study, a contrast agent having no cerebral vascular permeability, such as an iodine contrast agent, is used as a tracer. The iodine contrast agent is injected from the cubital vein by an injector, for example. The iodine contrast medium intravenously injected 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 medium passes through the capillary blood vessels in 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 from the continuous image, the time density curve Ca (t) of the pixel on the cerebral artery and the time density curve Ci (t) of the pixel on the brain tissue (capillary blood vessel) are obtained. The time density curve Csss (t) of the pixel on the cerebral vein is measured.

  Here, in the CBP study, an ideal relationship is established between the concentration time curve Ca (t) of the cerebral artery and the concentration time curve Ci (t) of the brain tissue as an analysis model. When contrast medium is injected from a blood vessel just before entering the brain tissue, the time density curve in the unit volume (1 pixel) of the brain tissue is vertical, 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-Modulation Transfer Function method).

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

The problems of the CBP study are as follows.
Since each index of CBP, CBV, MTT, and Err is calculated for each pixel (x, y, z), an image having the value as a pixel value can be configured, and such an image can be mapped. Call it. For example, when R types of indexes are obtained, R maps can be configured. The R maps created in this way can be regarded as one map (vector value map) in which each pixel takes a vector value. That is, 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 CBP value, P _{k, 2} (x, y, z) is the CBV value, P _{k, 3} (x, y, z) can be configured to represent the value of MTT, and P _{k, 4} (x, y, z) can be configured to represent the value of residual Err.

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

Thus, when the number K of the cerebral artery time concentration curve Ca (t) _k (k = 1, 2,..., K) is large, the resulting vector value map V _k (k = 1, 2). ,..., K) are inconvenient to observe because of the large number of sheets. That is, if a normal grayscale image or color scale image is to be observed, one map is composed of R images, and there are K sets, so a total of K × R images must be compared. . Further, it is not always obvious which part is nourished by which artery, and which map V _k (k = 1, 2,..., K) is observed for each part using anatomical knowledge. You must decide what to do. Especially in cases with cerebrovascular disorders such as cerebral infarction, which arteries depend on the tissue does not necessarily match the anatomical knowledge, and abnormal control is often seen . Due to these problems, there is a problem that it is difficult to interpret the vector value map.

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

  The CBP study also has the following problems. When bolus injection is performed on the cubital vein, the contrast effect observed by CT increases the CT value of blood (several tens of HU when not contrasted) to a maximum of several hundred HU. However, in order to effectively analyze cerebral blood flow, the contrast effect must be measurable with an error of no more than a few percent. That is, even if the blood contrast effect (CT value increase) is about 20 to 40 HU, it needs to be detected.

  The volume ratio of capillaries in a unit volume of brain tissue is at most about 3 to 4 percent. Therefore, when the CT value of blood increases by 20 to 40 HU, the average CT value of brain tissue only increases by about 0.5 to 1.5 HU.

  In CT images, the standard deviation (sd) of noise is inversely proportional to the square root of the irradiation X-ray dose, and sd is, for example, about 5 to 10 HU under typical irradiation conditions. Therefore, in order to detect a contrast effect of 0.5 HU, the X-ray dose must be increased by about 10 to 100 times, which means that the patient exposure dose is significantly increased. In addition, in dynamic CT, the same part is photographed several tens of times, so that the skin exposure at the photographing part is several hundred to several thousand times the usual, and inflammation, hair loss, necrosis, Considering radiation damage such as carcinogenesis, it is not realistic.

  Rather, in dynamic CT, the X-ray dose must be reduced compared to normal imaging. In general, the X-ray dose per scan is reduced to, for example, about 1/2 to 1/10 of the usual amount. As a result, the X-ray exposure can be reduced to several to 20 times that of normal one CT imaging, which is a level that does not cause radiation damage. However, in such a CT image with a reduced X-ray dose, sd is, for example, about 15 to 20 HU, and a contrast effect of about 0.5 to 1.5 HU cannot be detected.

  Therefore, suppressing the noise component of the image is one of the important issues in the CBP study. Therefore, 1) increase the slice thickness, 2) average neighboring pixels, and 3) pass the image through a smoothing process. However, these have the following problems.

  In order to “increase the slice thickness”, the slice thickness is set to be thick at the time of shooting, or an image of a thick slice is generated by averaging several consecutive thin slice images. Since the X-ray dose per pixel increases in proportion to the slice thickness, the image noise sd decreases in inverse proportion to the square root of the slice thickness. However, by increasing the slice thickness, a partial volume effect occurs, that is, one pixel does not represent a uniform brain tissue, and a plurality of tissues (white matter, gray matter, blood vessel, sulcus, ventricle) Etc.), and the error in the value of the cerebral blood flow obtained as an analysis result increases.

  In particular, normal analysis is not possible for pixels including the influence of blood vessels. For this reason, when the slice thickness is increased, only a very poor result including many pixels that are inaccurate and cannot be analyzed can be obtained.

  “Average neighboring pixels” sacrifices some spatial resolution. For example, an average value of a square area (including n × n pixels) having n pixels on one side is obtained, and this is set as an average CT value of the entire square. Such a square is regarded as a pixel, It spreads out to form a “pixel bundled image”. For example, if the original image is composed of 512 pixels on one side (including 512 × 512 pixels) and n = 2, the “pixel bundle image” is composed of pixels on one side (512/2). (Including 256 × 256 pixels). According to this method, noise can be reduced in inverse proportion to n. Furthermore, since the number of pixels to be analyzed becomes 1 / (n × n) times, there is an advantage that the amount of calculation is also reduced.

  However, when n is increased, the spatial resolution is reduced, and a partial volume effect is caused accordingly. That is, one pixel does not represent a uniform brain tissue, and a plurality of tissues (white matter, gray matter, blood vessels) , Cerebral sulcus, ventricle, etc.) is likely to represent an average CT value, and an error in values such as cerebral blood flow obtained as an analysis result increases. In particular, normal analysis is not possible for pixels including the influence of blood vessels. For this reason, if n is increased, only a very poor result including a low spatial resolution, inaccurate and unanalyzable pixels can be obtained. For this reason, in practical use, n = 2 to 4 is the limit, and a sufficient noise suppression effect cannot be obtained only by this.

  Further, if a method of smoothing an image, that is, smoothing by applying a two-dimensional spatial filter for each CT image, the spatial resolution is significantly impaired in exchange for a sufficient noise suppression effect. In particular, a pixel close to a location where a thick blood vessel (artery / vein) exists is affected by the contrast effect generated in the thick blood vessel, and the time density curve of these pixels becomes incorrect. Therefore, only a slight smoothing must be performed. Here, it is important to make the size of the image filter very small, for example, to set it to about 3 × 3 when performing very smooth smoothing. When trying to obtain the maximum image noise suppression effect using a 3 × 3 smoothing filter, the upper limit is to reduce the noise sd to 1/3, and it is impossible to suppress the noise further. It is. Therefore, a sufficient noise suppression effect cannot be obtained.

  On the other hand, temporal smoothing, that is, using a technique in which the time density curve obtained for each pixel is regarded as a curve and smoothed with a one-dimensional filter, the time resolution is reduced to obtain a sufficient noise suppression effect. Significantly damaged. Originally, dynamic CT in CBP studies is performed with a short sampling period to obtain a high time resolution, and a slight and fast change in the time density curve (especially how smoothing effect is due to physiological structure). Since the purpose is to accurately measure whether or not it has occurred, temporal smoothing is not at all appropriate.

  An object of the present invention is to improve the diagnostic efficiency of a map obtained by a CBP study.

According to a first aspect of the present invention, a first time concentration curve related to an artery in the specific site and a tissue related to the tissue in the specific site are obtained from a plurality of continuous images related to the specific site of the 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 to each of the second-time concentration curves, the local blood flow dynamics of each of the tissues are most likely dependent. A method is provided for identifying one of the arteries that is higher and creating a map that distinguishes the dependent region of the artery based on one of the arteries identified for each tissue.
According to a second aspect of the present invention, there is provided a first temporal concentration curve related to an artery in the specific site and a tissue related to the tissue in the specific site from a plurality of continuous images related to the specific site of the subject injected with a contrast agent. There is a possibility that each of the tissues may be subordinated by selecting a time concentration curve creating unit that creates a 2 hour concentration curve and one of the first time concentration curves that best fits to each of the second time concentration curves. An apparatus comprising: means for identifying one of the highest arteries; and a map creating unit for creating a map for distinguishing a dependent region of the artery based on one of the arteries identified for each tissue. I will provide a.

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

Hereinafter, the present invention will be described by way of preferred embodiments with reference to the drawings.
As a feature of the present embodiment, it is possible to improve the diagnostic efficiency of the CBP study by synthesizing many index maps generated by the CBP study into one sheet, and further reduce noise by using a coherent filter. The aim is to improve the accuracy of the index by achieving both reduction in space and temporal resolution.

  Note that the present embodiment relates to a method and apparatus for calculating an index representing local blood flow dynamics from a plurality of temporally continuous images related to a specific part of a subject. The generated modality is not limited to a specific apparatus. For example, X-ray computer tomography apparatus (X-ray CT apparatus), single photon emission tomography apparatus (SPECT), positron emission tomography apparatus (PET), magnetic Any of resonance imaging apparatus (MRI) may be used. Here, an X-ray CT apparatus will be described as an example.

(Device configuration)
FIG. 2 shows the configuration of the X-ray CT apparatus according to the present embodiment. The X-ray CT apparatus includes a gantry unit 10 and a computer device 20. The gantry unit 10 includes an X-ray tube 101, a high voltage generator 101a, an X-ray detector 102, and a data acquisition unit 103 (DAS; Data Acquisition System). The X-ray tube 101 and the X-ray detector 102 are mounted at positions facing each other with the subject P sandwiched between rotating rings (not shown) that rotate continuously at high speed.

  The computer device 20 includes an image processing device 30, an image display unit 107, and an input unit 109. The image processing apparatus 30 has a control unit 108 as a center, a pre-processing unit 104 that converts raw data output from the data collection unit 103 into projection data through correction processing and the like, a memory unit 105 that stores projection data, and projection data Image reconstructing unit 106 for reconstructing CT image data from, storage device 10M for storing CT image data, coherent filter processing unit 110 for performing coherent filter processing on CT image data, and CT that has undergone coherent filter processing And a CBP study processing unit 120 that executes CBP study processing using image data.

  The coherent filter processing unit 110 includes a variance value estimation unit 111, a weight function calculation unit 112, and a pixel value calculation unit (coherent filter unit) 113. The functions of the variance value estimation unit 111, the weight function calculation unit 112, and the pixel value calculation unit 113 will be described in the detailed description of coherent filter processing described later.

  The CBP study processing unit 120 includes an ROI setting support unit 121, a time concentration curve creation unit 122, a cerebral artery time concentration curve correction unit 123, an MTF processing unit 124, an index calculation unit 125, a map creation unit 126, and a map composition unit 127. Is done.

  The ROI setting support unit 121 provides information (AT map, PT map, TT map, etc. for the cerebral artery ROI) for supporting the operation of setting the region of interest ROI for the cerebral artery and cerebral vein on the CT image. Create and provide.

  The cerebral artery ROI is individually set in each of the left and right brain regions, for example, for the anterior cerebral artery (ACA), middle cerebral artery (MCA), and posterior cerebral artery (PCA). Accordingly, in this example, a total of six cerebral arteries ROI are set, three on each side. Further, another time concentration curve Csss (t) is used to correct the time concentration curve Ca (t) of the cerebral artery. This time density curve Csss (t) is created for an ROI set on a blood vessel that is sufficiently thick so that there are pixels that do not contain a partial volume. The ROI of Csss (t) is set to, for example, the thickest superior sagittal sinus in the cerebral blood vessel.

  The time concentration curve creation unit 122 generates time concentration curves related to cerebral arteries, cerebral veins, and brain tissues (capillaries) from dynamic CT image data (a plurality of temporally continuous image data) stored in the storage device 10M. create. The cerebral artery time concentration curve Ca (t) is individually created for, for example, six cerebral arteries ROI that have been set. The cerebral venous time concentration curve Csss (t) is created for the cerebral venous ROI set in the upper sagittal sinus. Further, the time density curve Ci (t) of the brain tissue is created for each pixel for all the pixels on the brain tissue.

  The cerebral artery time concentration curve correction unit 123 uses the cerebral artery time concentration curve Ca (t) based on the upper sagittal sinus time concentration curve Csss (t) in order to remove the influence of noise and the partial volume effect. To correct. This correction method will be described later. The MTF processing unit 124 converts the transfer function MTF into the brain tissue region by the box-MTF method based on the corrected time concentration curve Ca (t) of the cerebral artery and the time concentration curve Ci (t) of the brain tissue. The calculation is performed for each pixel with respect to all the pixels.

  The index calculation unit 125 calculates an index (CBP, CBV, MTT, Err) representing the blood flow dynamics of the brain tissue from the calculated transfer function MTF for every pixel in the brain tissue region. The map creation unit 126 generates a map of each calculated index for each cerebral artery (ACA, MCA, PCA). For each slice, 12 types of index (= 4) × the number of cerebral arteries (3 of ACA, MCA, and PCA) = 12 types are created. In multi-slice, the number of types of maps is created that is the number of slices. A map composition unit 127 is provided in order to reduce the enormous number of maps by the composition process and improve the diagnosis efficiency.

Hereinafter, the coherent filter process and the CBP study process will be described in order.
The principle of the coherent filter will be briefly described. The weights of neighboring local pixels such as 3 × 3 are weighted and the weighted average value is set as the value of the local central pixel. It is characterized by changing according to the similarity between the pixel and the peripheral pixel. The similarity referred to here is an index indicating the degree of possibility that the tissues are close to anatomical, specifically, the brain tissues (capillaries) under the control of the same cerebral artery, By giving a high weight to pixels with a high degree of similarity and conversely giving a low weight close to zero to a pixel with a low degree of similarity, while suppressing noise, the reduction in spatial resolution is suppressed. Making it possible. In the following description, the degree of similarity is appropriately rephrased as fitness or risk factor.

  In the present embodiment, a plurality of CT images (dynamic CT images) obtained continuously by imaging a subject's brain into which a contrast agent having no cerebral vascular permeability, for example, an iodine contrast agent is injected (intravenous injection), are obtained. Used to calculate the similarity by comparing the time density curves of each pixel. Therefore, the probability of similarity is determined depending on the sampling frequency, that is, the number of images per unit time, and the number of samplings, that is, the total number of images. Therefore, it is effective to shorten the scan interval to 0.5 seconds, for example.

(Coherent filter)
(General description of coherent filter)
(Pixel value v (x))
In general, a digital image acquired via an imaging means such as a camera or a CT scanner is composed of a plurality of pixels (or the image can be considered as a set of such pixels). . In the following description, the position of the pixel is expressed as a vector x (that is, a vector of coordinate values), and a value (for example, a numerical value indicating shading, CT value HU) of the pixel x is expressed as a K-dimensional vector. In the case of a two-dimensional image, the pixel x is a two-dimensional vector indicating coordinate values (x, y) representing a position on the image. The “pixel value v (x)” defined for a pixel x is
v (x) = (v ₁ (x), v ₂ (x),..., v _K (x)) (1)
Is written. Note that each of v ₁ (x), v ₂ (x),..., V _K (x) on the right side of the equation (1) is hereinafter referred to as a “scalar value” for the pixel x.

For example, when the image is a “color image”, each pixel has brightness (scalar value) of three primary colors (red, green, and blue), and therefore the pixel value v (x) of each pixel is It can be considered that the dimension is a vector of K = 3. That is, the subscripts on the right side of the above equation (1) are associated with, for example, “red”, “green”, and “blue”. For example, if the image is a moving image composed of K still images and each pixel of the nth image has a scalar value v _n (x), K-dimensional vector values v _n (x) = (v ₁ (x), v ₂ (x),..., Arranged by arranging pixel values (scalar values) of pixels x at the same common point (same coordinates). v _K (x) is a pixel value as a vector value described below.

(Similarity (fitness or risk factor) and weight)
Consider an appropriate set of neighboring pixels N (x) 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 that is an element of N (x) with respect to the center pixel x. The weight w (p (x, y)) has the following properties.

(Similarity p (x, y))
First, the meaning of the function p (x, y) that affects the value of the weight w (p (x, y)) will be described. This p (x, y) is a means for quantifying the “similarity” referred to in the present embodiment. Generally speaking, the central pixel x and the peripheral pixel y∈N (x) are in some sense. A specific numerical value indicating how similar (for example, the degree of statistical difference recognized between the pixel values v (x) and v (y) of both pixels x and y) is given.

  More specifically, for example, when p (x, y) shows a large value, “no statistically significant difference (that is, the degree of similarity is high)” between the pixel x and the pixel y. When it is determined that there is a high possibility and p (x, y) indicates a small value, “there is a statistically significant difference (that is, the similarity is low)” between the pixel x and the pixel y, It is that it is judged as follows.

By the way, the pixel values v (x) and v (y) (or scalar values v ₁ (x),..., V _K (x) and v ₁ (y),..., V _K (y)) are always noise. It is included. For example, when considering a case where an image is acquired by a CCD image pickup device, each pixel constituting the image includes noise or the like due to dark current in the device or irregular fluctuation of the amount of light incident from the outside.

  Since such noise generally takes various values for all pixels, even if the pixel x and the pixel y reflect the same object (in the outside world), they are actually observed. On the image, it may not have the same value. In other words, if it is assumed that the respective noises are removed from the pixel x and the pixel y reflecting the same object, they are displayed on the image as representing the same object. (= Recognized as such) and both have essentially the same (or very close) pixel values.

  Therefore, based on the above-mentioned noise characteristics, regarding the above p (x, y), using the concept of the “null hypothesis” that is well known in the statistical test method, this p (x, y) Specifically, we can say as follows. In other words, the null hypothesis H “pixel x and pixel y have the same pixel value when the respective noises are removed”, in other words, “v (x) = v (y)” "Excluding the difference caused by" (that is, when such a proposition is established, it is considered that "the similarity between both pixels x and y is high (the degree of matching is high)"). x, y) can be configured as a risk factor (or significance level) when rejecting this hypothesis H (in this case, p (x, y) is a function whose range is [0, 1]. (P (x, y) ε [0,1])).

  Therefore, when the risk factor p (x, y) is large, that is, when there is a high risk of rejection being incorrect, it can be said that the above hypothesis H is likely to be satisfied, and conversely, when it is small, that is, rejection is incorrect. If the risk is small, it can be said that there is a high possibility that the hypothesis H is not satisfied (note that although it is a well-known matter in the statistical test, the hypothesis H is not “rejected”, it is “true”. (In this case, it simply means that the proposition indicated by hypothesis H cannot be denied).

(Weight w (p (x, y)))
The weight w (p (x, y)) is a function of the risk factor p (x, y) as described above (more generally, a function of fitness (the fitness is If ρ (x, y), it can be configured to be w (ρ (x, y)), and in order to obtain this weight w (p (x, y)), x and y Generally speaking, the weighting function w to be applied to the risk factor p (x, y) obtained for each combination has the effect of embodying the “rejection.” Specifically, the risk factor. When p (x, y) is large, the value of the weighting function w, that is, the weight w (p (x, y)) takes a large positive value, and vice versa, a small positive value (or “0”). ")" (The specific form of the weighting function w will be described later), that is, the weight w (p ( x, y)) takes a large value if the pixel x and the pixel y seem to satisfy the proposition shown in the above hypothesis H, and takes a small value in the opposite case. It is also possible to configure so that there are only two possible values of “0” or a fixed value other than “0”.

  In addition, when the relationship between the hypothesis H, the risk factor p (x, y), and the weight w (p (x, y)) described above is summarized, when the null hypothesis H is highly likely to be correct, When p is also increased and the weight w given to the pixel is increased, on the other hand, when the possibility that the null hypothesis H is correct is low, the similarity p is also lowered and the weight w given to the pixel is lowered. Thus, by changing the contribution (weight) to the weighted average value according to the similarity, it is possible to effectively suppress noise while suppressing a decrease in resolution. Further, the weight function w (t) can be more generally referred to as “a non-negative monotonically increasing function defined by t∈ [0, 1]”, and the property to be satisfied by w (t) is at least That is all that is necessary.

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

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

Here, the case where the pixel value v (y) = (v ₁ (y), v ₂ (y),..., V _K (y)) (y = x) represented by the above equation (2) is included. .) Is converted to v ′ (x) = (v ′ ₁ (x), v ′ ₂ (x),..., V ′ _K (x)) is a “coherent filter” format. As is apparent from the expression, this represents a weighted average value of the scalar values v _k (y) constituting the pixel values.

Such processing yields the following results. That is, the contribution of the pixel y that is likely to take the same pixel value as the pixel value v ′ (x) of the pixel x from which noise is removed (highly likely to satisfy the proposition of the above hypothesis H) is increased. It represents a vector composed of the weighted average value v ′ _k (x). Further, if there are a sufficient number of such pixels y, the pixel value v ′ (x) is not deviated from its true value that the pixel x should originally have, and noise is caused by the above-described averaging operation. Although the resolution is suppressed, the resolution does not decrease.

  Even if the risk factor p (x, y) is small and therefore the null hypothesis H is “rejected” and the weight w (p (x, y)) is small, the above description also indicates As you can see, this is not necessarily completely “rejected”. Although this depends on the specific form of the weight function w described later, even when the risk factor p (x, y) is close to “0” (= 0%), w (p ( x, y)) ≠ 0 (where p (x, y) is a smaller positive value than when p (x, y) is close to “1”) (p (x, y) = 1) (Some cases are when v (x) = v (y), as will be described later.)

  That is, it is not a complete rejection, but may be configured to allow a small contribution (in this case, if w (p (x, y)) = 0, It is synonymous with complete rejection.

  Such processing can be generally described as follows. That is, when there are a plurality of pixels x constituting an image, the degree of matching p (x, y) between this pixel x and a given pixel y (in the above, yεN (x)) is quantified. In the weighted averaging process using the pixel value v (y), when the fitness is large, a large contribution of the pixel y is recognized, and when the fitness is small, only a small contribution is recognized. By doing so, it can be said that this is an image processing method for effectively suppressing noise of the pixel x. In other words, when the pixel x and the pixel y are “similar”, the pixel y contributes more to the averaging process. When the pixel x and the pixel y are “similar”, the pixel y is completely or almost ignored. In other words, the weight may be zero or an approximate value thereof.

  By applying such processing to the entire image, an extremely high noise suppression effect can be exhibited with almost no blurring of the image, that is, a reduction in spatial resolution. Further, the present invention is not limited to the use of noise suppression. For example, also in the field of pattern recognition, an excellent effect can be exhibited by making a weighting function or a coherent filter suitable in a specific form.

  Here, “dynamic CT” imaging means that the X-ray tube 101 and the X-ray detector 102 repeatedly image the same part of the subject P (repetitive scanning, in a continuous rotation CT apparatus, repeated imaging by continuous rotation is performed). This is an imaging method in which projection data is acquired one after another, and reconstruction processing is performed one after another based on the projection data to obtain a series of time-series images (in this case, images The image display on the display unit 107 is controlled to be performed after a predetermined time from the scan start point or end point related to the collection of projection data that is the source of the image, for example, by a counter (not shown).

  Therefore, the image acquired / displayed in this way is a so-called moving image composed of a plurality of time-series still images as in a movie or the like. Note that such an imaging method typically injects a contrast medium into the subject P, observes and analyzes changes over time, and analyzes other pathological conditions such as stenosis and occlusion in blood vessels, for example. Used for. A method of performing CT imaging of the same part only twice before and after contrast medium administration can also be considered as dynamic CT imaging in a broad sense.

  Conventionally, at the time of “dynamic CT” imaging as described above, for example, some change in the subject P (for example, a change in the concentration of a contrast medium, respiratory movement, etc.) is generally considered while K imaging is performed. In order to suppress the image noise without impairing the spatial resolution, smoothing in the time direction was necessary. As a result, the adverse effect that the time resolution is impaired cannot be avoided.

  However, the image acquired by dynamic CT imaging is a moving image as described above, and is taken for the purpose of closely observing temporal changes, so that the time resolution is impaired. Originally, this is not a preferable situation.

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

First, for the pixel x defined for the K still images, which are the moving images obtained as described above, as described above, as the pixel value v (x),
v (x) = (v ₁ (x), v ₂ (x),..., v _K (x)) (1 reprint)
Can be configured. Here, subscripts 1, 2,..., K in each term on the right side are serial numbers of K still images.

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

  However, yεN (x) and the set N (x) may be arbitrarily set for each pixel x (= may be set according to any criterion). However, in practice, the pixel x and the pixel y far away from the pixel x can satisfy the hypothesis “v (x) = v (y). However, excluding differences caused by noise between the two pixels”. Therefore, limiting the set N (x) on the basis of a set of pixels close to x has practical significance such as an increase in calculation speed.

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

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

Hereinafter, these σ _k and C will be described in order.
First, σ _{k in the} equation (4) will be described (hereinafter, described as variance σ _k ² ). As described above, σ _k ² is the variance of the noise component of the scalar value of each pixel on the kth still image. The variance σ _k ² in the above equation (4) is estimated on the assumption that the scalar value of each pixel of the k-th image includes noise having a constant variance σ _k ^2. is there. In general, such assumptions are sufficiently valid against the background described below.

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

On the other hand, it is also known that CT image noise is additive and generally follows a Gaussian distribution. That is, for an arbitrary scalar value v _n (x) (n = 1, 2,..., K) constituting the pixel value v (x) of a certain pixel x, the true value (value obtained by removing the noise contribution). V _n ⁰ (x), these difference values v _n (x) −v _n ⁰ (x) generally follow a Gaussian distribution with an average of 0 and a variance σ _k ² (note that the irradiation X-ray dose or the tube current The time product mA · s and the noise variance σ _k ² are generally in inverse proportion.

The variance σ _k ² also depends on the position of the pixel x itself (as described above, for example, each coordinate value x = (x, y)), but in the normal X-ray CT apparatus 100, Between the X-ray tube 101 and the X-ray detector 102, a physical X-ray filter (for example, a so-called “wedge” or “X-ray filter” composed of a copper foil, a metal lump or the like) that adjusts the X-ray irradiation amount. So that it can be ignored. This is because the wedge is irradiated so that the same amount of X-rays can be detected by any X-ray detector 102 by using the fact that the subject P is made of a substance having the same density as water. It has the effect of attenuating a part of X-rays. Therefore, such a wedge results in the effect that the noise variance σ _k ^{2 is made} to be a substantially constant value almost independent of the position of the pixel x. The purpose is to effectively use the dynamic range of the detector 102).

From the above, on the K still images acquired by dynamic CT imaging, it is reasonable to estimate that the variance σ _k ² is almost constant for all pixels on the kth still image. is there. Of course, it can be easily inferred that the present embodiment can be extended in the case where the variance is different for each pixel.

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

Furthermore, in general, imaging may be performed by changing the irradiation dose (X-ray tube current × irradiation time (mA · s)) for each imaging, and in this case, the variance σ _k ² is different for each image. .

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

Can accurately estimate σ _k ² .

In the present embodiment, a specific numerical value of σ _k ² can be estimated for at least one k by the following method.

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

Is the expected value E [σ ² ] for the true variance σ ² ,

Can be obtained as The expected value E [σ ² ] of the variance can be considered to be appropriate for all pixels x on all K still images as described above, and can be used as a substitute for the true variance σ ^2. This is a value for which the certainty is guaranteed over a certain level. Therefore, in the actual calculation of the above equation (4), this E [σ ² ] may be substituted for σ ² of equation (4).

More specifically, such E [σ ² ] may be obtained from measured values based on, for example, the first and second still images in K still images (the above (6)). And, in terms of equation (7), this corresponds to N = 2). Further, for the pixel x _f to be subjected to the actual operation of the (6) and (7), for example, air or bone are selected only appropriate pixel x _f excluding the portion being imaged (s If selected, an average of all the obtained E [σ ² ] may be taken. Furthermore, in general, it is better to devise measures to suppress the influence of the movement of the subject P.

Even when the irradiation dose is not constant in taking these N images, it can be easily estimated that σ _k ² is correctly estimated by using the fact that σ _k ² is proportional to R _k .

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

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

It is nothing else. A is a constant and is standardized so that p1 has a value of (0 to 1).

  Eventually, in formula (4), the risk factor p (x, y) described in the general form as described above is not explicitly displayed, but the actual condition of the weight w1 (p (x, y)). Can be regarded as a function of the risk factor (= p1 (x, y)) as described above (equation (9)), that is, “function of fitness” (however, the risk factor and As described above, the degree of fitness has a relationship in which when one increases, the other increases.

  As can be seen from the above equation (9), the parameter C has an effect of determining how sensitively the weight w1 (p (x, y)) reacts to the risk factor p1 (x, y). That is, when C is increased, p1 (x, y) is slightly decreased, and w1 (p (x, y)) approaches 0. Moreover, when C is made small, such a sensitive reaction can be suppressed. Specifically, C may be about 1 to 10 and preferably C = 3.

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

In addition, as can be seen from the expression (4), in the present invention, after calculating the risk factor p (x, y) for each combination of x and y, the weight w (p (x, y)) It is not always necessary to go through the procedure of obtaining, and the configuration may be such that (w ⁰ p) as a composite function is directly calculated without specifically obtaining the risk factor p (x, y).

As described above, the variance σ ² is estimated (for example, E [σ ² ] in the equation (7)), and the parameter C is appropriately determined (for example, C = 3). For all pixels y included in the set N (x) defined for a certain pixel x (for example, an area corresponding to 3 × 3 pixels centered on the pixel x, for example) The weight w1 (p (x, y)) can be obtained. After that, by using this w1 (p (x, y)) instead of w (p (x, y)) in the above equation (2), it is possible to carry out a specific numerical calculation of the coherent filter. It becomes possible. As a result, the pixel values v ′ (x) = (v ′ ₁ (x), v ′ ₂ (x),..., V, in which noise is strongly suppressed without impairing the spatial resolution as well as the temporal resolution. ' _K (x)) (= Equation (3)), that is, such K still images or moving images can be obtained.

FIG. 3A to FIG. 3C illustrate such image processing so that it can be conceptually easily understood. That is, in FIG. 3A, in a still image of 1, 2,..., K still images, for a certain pixel x, a rectangular area N ^{3 × 3} (3 × 3 pixels centered on the pixel x). x) is assumed. The pixels in the left corner corner of the rectangular area N ^{3 × 3} (x), if _{y 1,} the pixel _{y 1} has a pixel value v _{(y 1).}

Then, the scalar value _v 1 constituting the pixel value _{v (y 1) (y 1} ), v 2 (y 1), ..., v K (y 1) and the scalar value of the pixel values v (x) _v 1 ( x), v ₂ (x),..., v _K (x), the weight w1 (p (x, y ₁ )) is calculated by the above equation (4) (FIG. 3B). The same applies to the remaining pixels y ₂ ,..., Y ₈ in the rectangular area N ^{3 × 3} (x), and as shown in FIG. 3B, w1 (p (x, y ₁ )),. , W1 (p (x, y ₈ )) and w1 (p (x, x)). In this case, the risk factor p (x, x) is “1” from the equation (8), and therefore the weight w1 (p (x, x)) is also “1” from the equation (9) (= maximum). Is weighted).

Next, the weights w1 (p (x, y ₁ )),..., W1 (p (x, y ₈ )), w1 (p (x, x)) obtained in this way are used as the corresponding pixels. Are multiplied by scalar values v _k (y ₁ ), v _k (y ₂ ),..., V _k (y ₈ ), v _k (x) in the _{k-th image} , respectively, to obtain the sum ((2 ), And this is divided by the sum of the weights w1 related to the rectangular area N ^{3 × 3} (x) (also corresponds to the denominator of equation (2)), the kth A scalar value v ′ _k (x) in which noise is suppressed can be obtained for the pixel x in the image (FIG. 3C). Further, for all images of k = 1, 2,..., K, the same weights w1 (p (x, y ₁ )),..., W1 (p (x, y ₈ )), w1 (p (x, x) )) To obtain a scalar value v ′ _k (x) in which noise is suppressed, thereby obtaining a pixel value v ′ _k (x) = (v ′ ₁ (x), v in which noise in the pixel x is suppressed. ' ₂ (x), ..., v' _K (x)) is obtained. If the above calculation is repeated for all the pixels x, K images with reduced noise can be obtained.

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

  Each process described above may be performed in accordance with a flowchart as shown in FIGS. 4A and 4B, for example, and calculation / image display related to each process may be performed. In order to realize it on the actual X-ray CT apparatus 100, for example, as shown in FIG. 2, an image processing unit 110 including a variance value estimation unit 111, a weight calculation unit 112, and a pixel value calculation unit 113 is provided. This can be done.

Among these, the weight calculation unit 112 is configured to obtain the weight w1 (p (x, y)) directly from the pixel values v (x) and v (y) as described above. Therefore, the said calculating part 112 is an apparatus which calculates | requires a weight directly, without calculating | requiring the value of risk factor p1 (x, y) concretely. In addition, it is not a structure as mentioned above, the risk factor calculating part (fitness quantification part) which calculates | requires the value of risk factor p1 (x, y) concretely, and weight w1 (p (x, x, y) based on the output It is also possible to employ a two-step procedure called a weight calculation unit for obtaining y)). In any case, the weight calculator 112 calculates the weight w1 (p (x, y)) using the variance σ ² estimated by the variance value estimator 111 and v (x) and v (y).

  In addition, the pixel value calculation unit 113 uses the pixel values v (x) and v (y) and the weight w1 (p (x, y)) numerically calculated by the weight calculation unit 112 to generate a pixel value v ′ ( x) is calculated. That is, the calculation unit 113 actually performs processing for suppressing noise of the original image, that is, application of a coherent filter (hereinafter, this is expressed as “applying a coherent filter”).

  In the case of applying a coherent filter to a moving image composed of K still images in the dynamic coherent filter processing as described above, the processing in the image processing unit 110 is performed after once reconstructing all still images. These may be stored in the storage device 10M and subjected to a coherent filter later as post-processing. However, the present embodiment is not limited to such a form. In the flow of continuous projection data collection, continuous reconstruction, and continuous display, a process for applying a coherent filter may be performed in real time (hereinafter referred to as “real time coherent filter process”).

In a preferred embodiment of real-time coherent filtering, each time a new image is taken and reconstructed, the following processing is performed. Among the first obtained image (image number 1) to the latest image (image number M), common identical points on K still images having image numbers M, M-1,..., M-K + 1. The pixel values (scalar values) of the pixel x at (same coordinates) are arranged and the K-dimensional vector value v (x) = (v _M (x), v _M−1 (x),..., V _{M−K + 1} (x) ). In this way, a coherent filter can be applied in exactly the same manner as the “dynamic coherent filter processing” described above. However, the pixel value calculation unit 113 does not actually calculate all the elements of the pixel value v ′ (x), but only the scalar value v _M ′ (x) corresponding to the latest image (image number M). calculate. As a result, since the calculation speed is improved, the latest image in which noise is suppressed can be displayed in real time.

As another preferred embodiment of this “real-time coherent filtering”, when the first K images are obtained, v ₁ ′ (x) _,. ′ (X) is obtained, and thereafter, using K still images having K-dimensional vector values having image numbers M, M−1,..., M−K + 1, v (x) = (v _M (x ), V _M−1 ′ (x),..., V _{M−K + 1} ′ (x)), and the real-time coherent filter processing described above may be applied to this. It is convenient that the dimension K of the pixel value vector v (x) can be changed at any time by manual setting or automatic setting during the real-time coherent filter processing.

  In this way, the CBP study is performed using a CT image in which only noise is effectively suppressed without reducing the spatial and temporal resolution by the coherent filter, and the local blood flow dynamics in the tissue, that is, in the local tissue. By quantitatively analyzing the dynamics of the blood flow passing through the capillary blood vessels and obtaining indexes (CBP, CBV, MTT, Err) representing the local blood flow dynamics, improvement in accuracy and reliability can be expected.

As described above, the CBP study process is performed on an image in which a reduction in resolution is suppressed and noise is removed.
(CBP study)
As described above, in the CBP study, indexes of CBP, CBV, MTT, and Err that quantitatively represent the dynamics of “blood flow through capillaries” in brain tissue are obtained, and the spatial distribution of these indexes is expressed. Output the map.

CBP: Unit volume and blood flow per unit time in capillaries of brain tissue [ml / 100 ml / min]
CBV: Blood volume per unit volume in brain tissue [ml / 100ml]
MTT: Average blood passage time in capillaries [sec]
Err: Indicator of residual deviation of actual measurement value from analysis model Note that a low residual Err means that there is a high possibility of being controlled by the reference cerebral artery, and conversely a high residual Err. This means that it is less likely to be dominated by the reference cerebral artery.

  In the CBP study, a contrast agent having no cerebral vascular permeability, such as an iodine contrast agent, is used as a tracer. The iodine contrast 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 to the cerebral vein through the capillaries in the brain tissue. At this time, a contrast agent having no cerebral vascular permeability, for example, an iodinated contrast agent, passes through the capillary in normal brain tissue without leaking out of the blood vessel.

  The state of passage of the contrast agent is continuously photographed by dynamic CT, and the time density curve Ca (t) of the pixel on the cerebral artery and the time density curve Ci of the pixel on the brain tissue including the capillary are obtained from the continuous image. (T) The time density curve Csss (t) of the pixel on the cerebral vein is measured.

  In the CBP study, the blood concentration of the contrast agent is ideal between the time curve Ca (t) of the blood concentration of the cerebral blood vessels close to the brain tissue and the time curve Ci (t) of the blood concentration of the capillary blood vessels. Are used as analysis models. In other words, when a contrast medium is injected from a blood vessel just before entering the brain tissue, the time density curve in the brain tissue unit volume (1 pixel) including the capillaries 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 Function method).

  Using the cerebral artery blood time concentration curve Ca (t) as the input function and the brain tissue time concentration curve Ci (t) as the output function, the characteristics of the process of passing through the capillaries are obtained as a transfer function represented by a rectangular function. Can do.

(Specific steps)
5 and 6 show a typical procedure of the CBP study according to the present embodiment. First, bolus injection (contrast medium is administered all at once) is performed on a blood vessel such as an elbow vein, and dynamic CT (images of the same part are repeatedly taken) is performed immediately after or immediately before. As a most typical procedure, when bolus injection is performed on the elbow vein, imaging is repeated for approximately 20 to 40 seconds, for example, at intervals of 0.5 to 2 seconds. The CT value of each j-th pixel (x, y) in N CT images obtained by dynamic CT is represented by v (x, y, j). This is nothing but a sample of the time density curve (smooth curve) f (t, x, y) at this pixel (x, y).

  First, as preprocessing, in step S1, pixels that are clearly identified as tissues other than brain tissue are excluded from the analysis target from each CT image. That is, a pixel indicating a value that does not fall within a range that can be considered as a CT value of brain tissue (for example, a CT value of 10 to 60 HU) is a pixel corresponding to air, bone, fat, etc. So these can be ignored. This analysis range is set to 10 to 60 HU as a default, but can be arbitrarily set via the input unit 109.

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

And this value is b (x, y). For pixel values v (x, y, j) of each image of j = K + 1, K + 2,.
q (x, y, j) = v (x, y, j) -b (x, y)
About j <K
q (x, y, j) = 0
And it is sufficient. In order to simplify the processing, the same K may be adopted for any pixel. Q (x, y, j) obtained in this way is considered to be nothing but a sample of a smooth continuous curve time density curve q (t, x, y) sampled at t = t1, t2,... TN. it can. Quantitative analysis of cerebral blood flow is performed using this q (t, x, y).

  In quantitative analysis, first, the right brain area and the left brain area can be separated on a CT image. As described above, in the CBP study, the blood flow dynamic state of the capillary is obtained as the transfer function MTF of the brain tissue time concentration curve Ci (t) with respect to the cerebral artery time concentration curve Ca (t). If the brain tissue to be analyzed is not under the control of the cerebral artery of the reference curve Ca (t), it is useless to calculate. At least the left brain and the right brain are separately analyzed using the time concentration curves Ca (t) of the separate cerebral arteries. That is, the time concentration curve Ca (t) of the left cerebral artery is used only for analyzing the brain tissue of the same left brain. Similarly, the time concentration curve Ca (t) of the cerebral artery of the right brain is used only for the analysis of the same right brain brain tissue, which has the effect of reducing useless calculations.

  In order to divide the brain into the left brain area and the right brain area, as shown in FIG. 7, a dividing line is displayed on the screen as a figure superimposed on the CT image (S3). The dividing line may be initially displayed at the center of the image. The operator refers to the image, moves the dividing line, and moves the plurality of constituent points constituting the dividing line to arbitrarily bend, thereby dividing the left and right areas.

  Thus, by dividing the brain into a left brain area and a right brain area and limiting the analysis range to each area, the number of analysis processing steps can be reduced. That is, the time concentration curve Ca (t) of the left brain cerebral artery (left ACA, left MCA, left PCA) is used only for analysis of the left brain area (transfer function optimization processing), and similarly, the right brain cerebral artery (right ACA) , Right MCA, right PCA) time concentration curve Ca (t) is used only for analysis of the same brain tissue of the 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 regions, and the analysis processing may be limited to a narrower area.

  For this region division, a geometric projection, for example, Voronoi projection is employed. As is well known, the Voronoi projection is a technique that is often used in fields such as the optimal layout of facilities such as hospitals, stores, and fire stations, and is based on many points (corresponding to stores, etc., mother points) arranged on a plane. The plane is divided into a plurality of power areas according to the distance.

  As shown in FIG. 8, in this embodiment, the Voronoi projection is applied individually to the left brain area and the right brain area. The left brain area is divided into a power region of left ACA, a power region of left MCA, and a power region of left PCA with left ACA, left MCA, and left PCA as three generating points. A Voronoi point is set at the center of a circle passing through three generating points corresponding to the left ACA, left MCA, and left PCA. Centered 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 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 power zones. Similarly, the right brain area is divided into a right ACA power region, a right MCA power region, and a right PCA power region, with the right ACA, right MCA, and right PCA as three generating points.

  A 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 obtained for each pixel while being limited to the power region of the left ACA. Similarly, the transfer function MTF is obtained for each pixel by limiting the left MCA, the left PCA, the right ACA, the right MCA, and the right PCA to each power region.

  In this way, by dividing the left brain area and the right brain area into a plurality of power regions and limiting the analysis range to each power region, the number of analysis processing steps can be further reduced.

  Next, the cerebral artery ROI is set on the cerebral artery on the CT image. In order to improve and facilitate the accuracy of this setting, a support map is created by the ROI setting support unit 121. Separately or overlaid on the CT image (S4). Examples of the support map include an AT (appearance time) map, a PT (peak time) map, and a TT (transit time) map, as shown in FIGS. 9A to 9C. For each pixel, as shown in FIG. 10, time AT, time T0 from an arbitrary time before contrast (for example, data collection start time) T0 to a time when the contrast agent concentration reaches several percent (for example, 1 percent) of the peak peak The time (peak time) PT from the time until the contrast agent concentration reaches the peak, or TT that represents the moving time of the contrast agent, for example, by a half-value width, is calculated and generated and displayed as a map. By default, all types of the AT map, PT map, and TT map are generated and displayed, but the operator can select any one type or any two types.

  These values tend to appear at higher values in the cerebral artery compared to other tissues, so by displaying in color through a color look-up table that is set to display only pixels with values centered on those values. The location of the cerebral artery can be easily identified and the cerebral artery ROI can be accurately set (S5). Typically, in each of the left and right brain areas, three cerebral arteries ROI are set: an anterior cerebral artery (ACA), a middle cerebral artery (MCA), and a posterior cerebral artery (PCA).

  In the case of multi-slice imaging, for example, when four adjacent slices are to be analyzed, setting the cerebral artery ROI individually in each slice as shown in FIG. This is not necessary for analysis. Therefore, the cerebral artery ROI set in one arbitrary slice is also shared with other slices. Alternatively, a cerebral artery time concentration curve Ca (t) that can be commonly used in all slices may be created using a coherent regression method described later.

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

  Here, many cerebral arteries are very thin compared to the pixel size, and generally, they are not orthogonal to the CT imaging slice, so that every pixel on the image does not accurately represent the CT value of arterial blood. First, one pixel is composed of a mixture of cerebral arteries and other tissues, and because of its partial volume effect, it usually shows only a lower contrast effect. Further, image noise is large in these pixels including the artery as a partial volume. In particular, in an artery or the like at a site where cerebral infarction occurs, the effect of noise is enormous because the contrast effect is relatively small. Although the image noise is suppressed by the above-described coherent filter, the influence of the partial volume effect still remains.

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

  Further, according to this method, a time density curve of only one slice image corresponding to each artery can be obtained, and therefore, it can be used for analysis of an arbitrary part in all slices within the imaging range. Allows you to select the slice for which the time density curve of the cerebral artery is most clearly obtained for a particular artery, and apply the time density curve of that cerebral artery to all slices. Can be reduced.

(Coherent regression method)
In creating the time concentration curve, it is important to remove the effects of the partial volume effect and random noise. First, the “time density curve” is a curve representing a change with time of a CT value (density value) at a specific part in the dynamic CT image. In particular, in the above-mentioned medical image diagnostic apparatus, the temporal change of the contrast medium concentration in a specific tissue of a human body is used as a time concentration curve for the purpose of examining details such as blood flow dynamics and metabolic functions in the human tissue. Measuring is done. In astronomical observation or the like, a time concentration curve is used for the purpose of analyzing a change in luminous intensity of a specific celestial object. And more formally specified, i.e., the time-density curve, when the density value of a site at time _{t k} and _{d k,} the pair row _{{<t k, d k>} (k = 1,2, · .., K)}. Also, for many applications of time density curves, it is not always necessary to have an absolute value of d _k , but rather it is sufficient if only the increment (d _k −d ₁ ) relative to the first image 1 is obtained. 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 pair of columns {<t _k , A (d _k −d ₁ )> (k = 1, 2,..., K)} is the time density curve to be obtained.

In order to obtain such a time density curve, in principle, an attempt is made to measure the time density curve in each image k (k = 1, 2,..., K) constituting the dynamic CT image. Using the scalar value v _k (x) of the pixel x included in the region to be paired, {<t _k , v _k (x)>} or {<t _k , A (v _k (x) −v ₁ (x)) >>.

  However, in practice, there is a problem that the time density curve to be originally measured cannot be obtained accurately because random noise is included in the dynamic CT image taken by the medical image diagnostic apparatus or the like.

Furthermore, in practical use, a so-called “partial volume effect” occurs in these dynamic CT images. The partial volume effect means that an image of a minute object in a subject is represented by a small number of pixels on the image, and an image of an adjacent object in the subject is also included in these few pixels. This is a phenomenon in which the pixel values of these small number of pixels exhibit relatively small fluctuations (although they are proportional to the fluctuations of the density value to be measured). In other words, the pixel values of these few pixels contain only a few signals. Therefore, when the partial volume effect is generated, the pair of columns {<t _k , v _k (x)> (k = 1, 2,..., K)} is extremely high regardless of which pixel x is taken. Since the signal level is low and affected by the change of the concentration value in the tissue that is not originally measured, and there is random noise, the time concentration curve to be originally measured {<t _k , d _k There is a problem that it is not possible to accurately obtain>

  Therefore, in the past, temporal or spatial smoothing has been used to suppress random noise, but temporal resolution is lost when temporal averaging is performed, and when spatial averaging is performed, There has been a problem that a change with time in the concentration of a portion other than the portion to be measured is mixed into the measured value. In order to solve such a problem and obtain a more accurate time density curve, a coherent filter is employed.

First, the null hypothesis to be used in the coherent filter of this embodiment will be described. When it is assumed that the true time concentration curve at the site to be measured is {<t _k , d _k > (k = 1, 2,..., K)}, {<t _k , A (d _k −d ₁ )> (k = 1, 2,..., K)} (where A is an unknown coefficient), the part to be measured A set R of pixels substantially corresponding to is set. For any pixel x∈R that is an element of the set R, the condition Q: “If this pixel x well reflects the true time density curve, it is almost affected by changes in density over time in other parts. If not, the partial volume effect and random noise of the pixel value v (x) = (v ₁ (x), v ₂ (x),..., V _K (x)) as a vector value By considering the impact
v _k (x) = p (x) d _k + q (x) + γ _k (x) (11)
(K = 1, 2,..., K)
Can be assumed to hold. Here, p (x) and q (x) is different for each pixel x but are unknown coefficients that do not change depending on the image number k (i.e. corresponding to the imaging time t _k), a model of the partial volume effect It is. Γ _k (x) is a model of random noise, and has a value different for each pixel x and for each image number k, but its expected value is 0, and its statistical distribution is the pixel x. And image number k.

  According to the above assumption, if the proposition “both pixels x and y satisfy the above-mentioned condition Q” holds for any two pixels x and y that are elements of the set R, then It can be proved that the relation of the expression holds.

v _k (x) = a ₁ v _k (y) + a ₂ + ξ _k (k = 1, 2,..., K)
(12)
Here, a ₁ and a ₂ are unknown coefficients that differ for each pixel set x and y but do not change depending on the image number k (ie, the photographing time t _k ). Ξ _k is random noise, and the value is different for each pixel set x, y and for each image number k, but the expected value is zero.

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

By doing so, the equation (12) is derived. Here, a ₁ and a ₂ in the equation (16) are parameters representing the partial volume effect, and ξ _k in the equation (16) represents random noise.

From the above, the proposition that “pixels x and y both satisfy condition Q” is the null hypothesis H ₀ ′ “v _k (x) = a ₁ v _k (y) + a ₂ + ξ _k (k = 1,... , K) ”.

Next, the null hypothesis H ₀ ′ “v _k (x) = a ₁ v _k (y) + a ₂ + ξ _k (k = 1,..., K)” is substantially equivalent and actually. Describes how to convert the proposition to a form that can be tested. When this null hypothesis is described in mathematically exact expression, the null hypothesis H ₀ ′ “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 an average of 0 and variance (σ ² h (a ₁ )). Here, the coefficient h (a ₁ ) is
h (a ₁ ) = 1 + a ₁ ² (17)
It is. The equation (17) is immediately derived from the equation (16), which is the definition of a ₁ and ξ _k , and the general nature of the variance for random variables. Further, the value of the variance σ ² can be estimated easily and practically sufficiently accurately.

From the above, if the above constants a ₁ and a ₂ can be determined, it is possible to test the null hypothesis H ₀ ′. In practice, it is sufficient ^to obtain optimal estimates a ₁ ^to a ₂ ^to these constants.

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

Define The value of S (a) depends on the constant vector a = (a ₁ , a ₂ ), that is, the values of the above constants a ₁ and a ₂ . 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. In addition, as a specific calculation method of the linear least square method, various known methods can be used, and these known calculation methods are all very simple and require a very short calculation time. It is.

Thus, as a result of calculating the optimal estimated values a ₁ ^to a ₂ ^to the constants a ₁ and a ₂ , residuals defined by the following equations are calculated.
r _k ^to (x, y) = v _k (x) −a ₁ ^to v _k (y) −a ₂ ^to (19)
Can be calculated specifically. Therefore, using this residual r _k ^˜ , the above-described null hypothesis H ₀ ′ is substantially equivalent to the null hypothesis H ₀ ″ “r _k ^˜ (x, y) (k = 1,..., K ) mean 0, variance ^{_{(1+ (a 1 ~) 2}} ) σ 2 of normally distributed. "and can be translated. This is a concrete proposition that can actually perform the calculation of the test.

In addition, expression by vector

(However, the vectors a and ξ depend on the pixel set x, y.) And the following expression f (a ^~ , v (y)) = a ₁ ^~ v (y) + a ₂ ^~ ... ( 21)
In other words the null hypothesis _{H 0} 'using defined as a vector value function f, the null hypothesis _{H 0} "is" v (x) = f (a ~, v (y)) + ξ, ( provided that, xi] mean 0, variance _{^{(1+ (a 1 ~) 2}} ) follows the normal distribution of the sigma ^2.) ", and this is exactly the same format as the null hypothesis _{H 0} described above. That is, it is clear that this embodiment is a modification of the above-described coherent filter. Here, the ^{f (a ~, v (y} )) and is, i.e., the pixel value v (y) of the pixel y, optimally adjusted to the parameter a indicating the partial volume effect, a pixel x The pixel value v (x) is converted so as to have the highest matching degree.

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

The time density curve thus obtained is {<t _k , A (d _k −d ₁ )>} (which is a linear transformation of the true time density curve {<t _k , d _k >} at the pixel x. A is a measurement value approximating an unknown coefficient), and random noise is suppressed by the effect of the weighted average. In addition, as is apparent from the equation, a pixel value vector of another pixel y corrected for the influence of the partial volume effect is used. Furthermore, the present embodiment has a characteristic 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 the coherent filter. Therefore, according to the present embodiment, it is possible to obtain a time density curve in which the influence of the partial volume effect is suppressed and random noise is suppressed without impairing the time resolution. The method for obtaining the time density curve in this way is particularly referred to as “coherent regression method”.

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

In this application example, in many cases, arteries in human tissue are generally very thin, and thus an arterial image that appears on a tomographic image by CT produces a partial volume effect. Furthermore, it goes without saying that the image contains random noise. For this reason, in the conventional method, it is difficult to obtain a sufficiently accurate time concentration curve related to the artery, and if measurement is performed forcibly, it is a linear transformation of the true time concentration curve <t _k , D _k > related to the artery. <T _k , A (D _k −D ₁ )> (where D _k represents a pixel value (which is a scalar value) of a group of pixels corresponding to an artery image at time t _k , and k = 1. , 2,..., K) to some extent, only measured values <t _k , (v _k (x) −v ₁ (x))> were obtained. This measurement includes random noise. Also, the coefficient A remains unknown due to the effect of the partial volume effect.

Therefore, measured values <t _k , (v ′ _k (x) −v ′ ₁ (x))> (k = 1, 2,...) Sufficiently approximating <t _k , A (D _k −D ₁ )>. .., K) can be obtained. On the other hand, some veins that can be observed on the same tomographic image are considerably thick. Therefore, for those veins, a sufficiently good approximation value <t _k , (J _k− J ₁ )> (k = 1, 2,..., K) can be obtained. Here, J _k represents a pixel value at a time t _k of a group of pixels corresponding to a vein image.

By the way, in the time concentration curve relating to blood circulation, the proposition S: “If the contrast agent concentration in the blood at time t ₁ is 0, the time concentration curve relating to which blood vessel d <t _k , (d _k -d ₁ )> is also known to have the property that the area under the curve (AUC) matches. The area under the curve here means the integration of the time density curve <t _k , (d _k -d ₁ )> with respect to time t.

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

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 can be calculated using equation (22). (J should be substituted for d.) Also, regarding the artery, if the time concentration curve {<t _k , (D _k −D ₁ )>} is known, the area under the curve AUC (D) is (18) can be calculated in the same way, and according to the proposition S,
AUC (D) ≈AUC (J) (24)
Should be true. However, in practice, the time concentration curve <t _k , (D _k −D ₁ )> is unknown, so AUC (D) cannot be calculated.

On the other hand, the time concentration curve <t _k , (v ′ _k (x) −v ′ ₁ (x))> obtained by the method according to the present embodiment is expressed as <t _k , A (D _k −D ₁ )>. The latter includes an unknown coefficient A. Therefore, the area under the curve AUC (v ′) that can be specifically calculated from {<t _k , (v ′ _k (x) −v ′ ₁ (x))>} using the expression (23) is AUC ( It must be exactly A times D). That is,
AUC (v ′) ≈A · AUC (D) (25)
It is. That is, from the equations (24) and (25),
A≈AUC (v ′) / AUC (J) (26)
This relationship holds. Since the right side of equation (26) can be specifically calculated using equation (23), the unknown value of coefficient A can be specifically determined. Therefore, if the value of the coefficient A is used to construct a time concentration curve <t _k , (v ′ _k (x) −v ′ ₁ (x)) / A>, this is equivalent to the arterial time concentration curve <t _{It is} nothing but an approximation of _k , (D _k −D ₁ )>. A method for constructing a time concentration curve in which the unknown constant A value is determined using the area under the curve in this way is referred to as an “AUC method”.

  From the above, in the clinical use of the time density curve in dynamic CT images obtained by dynamic CT imaging or the like, it is difficult or impossible to measure with the conventional method by combining the AUC method with the coherent regression method. Even for a fine arterial time-concentration curve, which is possible, a measurement value that eliminates the influence of the partial volume effect and random noise and does not include the unknown constant A is obtained.

Of course, the AUC method can also be applied to a time-concentration curve <t _k , (v ′ _k (x) −v ′ ₁ (x))> relating to an artery measured by a conventional method alone (randomly). Although the influence of noise and partial volume effect cannot be excluded, it is possible to construct a time-density curve in which the value of constant A, which has been unknown, is determined.

(Correction of cerebral artery time concentration curve using upper sagittal sinus time concentration curve (suppression of partial volume effect))
In order to suppress the influence of partial volume, the time concentration curve Ca (t) of the cerebral artery is corrected using the time concentration curve Csss (t) of the superior sagittal sinus instead of or in combination with this coherent regression. You may make it do.

  First, in step S7 of FIG. 5, as shown in FIG. 12, a large upper sagittal sinus ROI is set so as to surround the upper sagittal sinus on the CT image. Since the upper sagittal sinus is larger than the artery and the position is relatively fixed, the setting of the upper sagittal sinus ROI is easy. The large upper sagittal sinus ROI includes a plurality of pixels.

  Next, the upper sagittal sinus ROI is reduced so that all the pixels in the upper sagittal sinus ROI are included in the upper sagittal sinus over the entire area (S8). As the reduction processing, for example, first, threshold processing (binarization) is performed for each pixel in the upper sagittal sinus ROI, and a binary map (“0”, “1”) in the ROI is obtained. create. The threshold value is set to a value that separates the image of the superior sagittal sinus from the images of surrounding tissues and bones. “1” represents a pixel on the image of the upper sagittal sinus, and “0” represents a pixel on an image of surrounding tissue or bone. Each pixel (center pixel) of this binary map is replaced according to the value of its neighboring 4 or 8 pixels. Only when the central pixel is “1” and all the neighboring 4 or 8 pixels are “1”, the value of the central pixel is maintained at “1”. That is, when the central pixel is “0”, even if it is “1”, if any one of the neighboring 4 or 8 pixels indicates “0”, the value of the central pixel is changed. Replace with “0”. Accordingly, the upper sagittal sinus ROI is reduced by at least one pixel from the outline of the image of the upper sagittal sinus. As a result, the condition that all the pixels in the upper sagittal sinus ROI subjected to the reduction process are pixels on the upper sagittal sinus sinus image can be realized with high accuracy.

  Instead of this method, the upper sagittal sinus ROI may be corrected using the area AUC under the time density curve. In this case, the area under the curve AUC of the time density curve is calculated for each of the pixels in the large ROI as a search range. Due to the contrast effect, the area under the curve AUC of the pixel on the upper sagittal sinus image clearly shows a higher value than that of the surrounding pixels. Therefore, by executing threshold processing on the area AUC under the curve, only pixels on the upper sagittal sinus image can be selected from the ROI.

  In this way, the time density curve of each pixel is averaged for a plurality of highly accurate pixels that are on the upper sagittal sinus image picked up by the AND condition using either method or both methods in combination. Then, a time concentration curve Csss (t) of the superior sagittal sinus is created (S9).

  Here, since the iodine contrast agent does not pass through the blood brain barrier, in principle, the iodine concentration does not change between the cerebral artery and the cerebral vein, that is, the temporal concentration curve Csss of the upper sagittal sinus ( The area under the curve AUC of t) is substantially equivalent to the area under the curve AUC of the time concentration curve Ca (t) of the cerebral artery created in S6. Therefore, as shown in FIG. 13, the area under the curve AUCa of the time concentration curve Ca (t) of the cerebral artery created in S6 with respect to the area AUCsss of the time density curve Csss (t) of the upper sagittal sinus. Are corrected by multiplying each time value of the time concentration curve Ca (t) of the cerebral artery created in S6 by AUC (sss / AUCa) (S10).

  Next, using the cerebral artery time-concentration curve Ca (t) shown in FIG. 14A in which the noise and the partial volume effect are suppressed as described above, the state of the blood flow dynamics of the brain tissue (capillary blood vessel) is shown. Quantify. For this purpose, first, a time density curve Ci (t) shown in FIG. 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 separate cerebral artery time density curves Ca (t) in the left and right areas, and the brain tissue time density curve is used as an input function. Using Ci (t) as an output function, the characteristic of the process of the tracer passing through the capillary is obtained as the transfer function MTF. That is, for the time density curve Ci (t) of the brain tissue in the left area, the time density curve Ca (t) of the same cerebral artery in the left area is used, and the time density curve Ci (t) of the brain tissue in the right area. , The transfer function MTF is obtained using the time concentration curve Ca (t) of the cerebral artery in the same right area. Moreover, since the cerebral artery time concentration curve Ca (t) is created for each ACA, MCA, and PCA as described above, the calculation of the transfer function MTF is repeated for each Ca (t).

  Here, as shown in FIG. 15, the ideal relationship established between the time concentration curve Ca (t) of the cerebral artery and the time concentration curve Ci (t) of the capillary is used as an analysis model, and box-MTF. Apply the law.

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

Based on the transfer function box-MTF that minimizes the residual, as shown in FIG. 15, CBP, CBV, and MTT are calculated (S13), and the sum of squares of the residual 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 ′
It is corrected by. Here, Ht is the hematocrit value of the large blood vessel, and b * Ht is the hematocrit value of the peripheral blood vessel (generally, b is about 0.7).

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 a 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 approximating the transfer function, and is calculated, for example, as shown in the following equation.

Here, S is a constant, for example, S = N−p. p represents the degree of freedom, that is, the number of parameters included in the model f to be approximated. w (t _i ) is a weighting factor that determines the degree to which the residual at time t _i contributes to Err. For example, w (t _i ) may be a fixed value such as 1 without depending on i,
w (t _i ) = αe ^{−ti2 / β}
As described above, the weight w may be gradually reduced as | t _i | increases, so that the residual does not significantly contribute to Err with time.

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

  Next, as shown in FIGS. 18A to 18D, respective maps are generated from the CBP, CBV, MTT, and Err that have undergone output optimization (S16). The CBP, CBV, MTT, and Err indices are calculated individually for each anterior cerebral artery ACA, middle cerebral artery MCA, posterior cerebral artery PCA, and for each slice. Accordingly, the map is 4 × 3 = 12 even in one slice as shown in FIG. If the number of cerebral arteries to be set is increased, the number of maps increases by four times the increase. It is not realistic to evaluate many such maps comprehensively. Therefore, to reduce the number of maps, the maps are synthesized (S17).

  As shown in FIG. 20, the synthesis method includes a CBP map of the anterior cerebral artery ACA, a CBP map of the middle cerebral artery MCA, and a CBP map of the posterior cerebral artery PCA, the residual Err of the anterior cerebral artery ACA, Based on the residual Err of the cerebral artery MCA and the residual Err of the posterior cerebral artery PCA. For example, when the transfer function MTF is obtained from the time concentration curve Ca (t) of the anterior cerebral artery ACA and the time concentration curve Ci (t) of the brain tissue under its control, the residual Err is relatively small, and conversely When the transfer function MTF is obtained from the time concentration curve Ci (t) of the brain tissue that is not under control, the residual Err is relatively large. That is, the residual Err represents the dominance possibility of each cerebral artery.

  Therefore, for each pixel, a CBP value corresponding to the lowest residual Err is selected from the CBP value of the anterior cerebral artery ACA, the CBP value of the middle cerebral artery MCA, and the CBP value of the posterior cerebral artery PCA. Select as the value of that pixel. The map synthesized in this way is composed of CBP values of brain tissue that are likely to be under the control of the anterior cerebral artery ACA, middle cerebral artery MCA, and posterior cerebral artery PCA. The same applies to the map synthesis of other indexes CBV and MTT.

Here, the map composition will be described in detail below. The time density curve obtained from the pixel located at the position corresponding to the artery reflects the arterial blood contrast medium concentration, and the above-mentioned coherent regression method is applied to the time density curve of the correct arterial blood contrast medium concentration. Can be obtained. Such a time-concentration curve of the cerebral artery can be created for each artery, and there are differences depending on the blood circulation state. In particular, this difference may be significant in cases that have cerebrovascular disorders. For example, the time concentration curve of the cerebral artery obtained from K arteries is represented by A _k (t) (k = 1, 2,..., K).

By comparing the time-concentration curve of a tissue with the time-concentration curve of the cerebral artery of the artery that is feeding (dominating) the tissue, the microcirculation (structure and function of the capillary system in the tissue) ) Reflecting CBP or the like can be obtained. Since these indexes are calculated for each part (x, y, z), an image having the values as pixel values can be configured, and such an image is an index map. For example, when R types of indexes (typically four types of CBP, CBV, MTT, and Err as described above) are obtained, R maps can be configured. The R maps created in this way can be regarded as one map (vector value map) in which each pixel takes a vector value. That is,
V _k (x, y, z) = <P _{k, 1} (x, y, z), P _{k, 2} (x, y, z),..., P _{k, R} (x, y, z) >
It becomes.

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

Among the parts (x, y, z), those that are clearly not corresponding to the organ to be analyzed are excluded from the analysis from the beginning, and P _{k, r} (x, y) May be substituted with a special value indicating that it is not subject to analysis (steps S14 and S15). As such a value, it is convenient to use a negative value having a large absolute value. Alternatively, as yet another element to be added to the vector V _k (x, y, z)
P _{k, R + 1} (x, y, z) = (0 if (x, y, z) is not subject to analysis, 1 otherwise)
You can make a map. Such a map is called a “mask”.

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

  Of these, the cerebral artery time-concentration curve obtained from the artery in the right hemisphere is only for the analysis of the part (x, y, z) belonging to the right hemisphere, and the time of the cerebral artery obtained from the artery in the left hemisphere. The concentration curve should be used only for the analysis of the part (x, y, z) belonging to the left hemisphere. Therefore, it is desirable that the operator designates the boundary (midline) between the right hemisphere and the left hemisphere as a straight line, a curved line or a polygonal line, a plane, a curved surface, or the like, and creates a map for each hemisphere. However, there may still be as many cerebral artery time concentration curves as K = 3-10 for each hemisphere on one side.

Thus, when the number K of the cerebral artery time density curve A _k (k = 1, 2,..., K) is large, the resulting vector value map V _k (k = 1, 2,... ., K) is inconvenient to observe because of the large number of sheets. That is, if a normal grayscale image or color scale image is to be observed, one map is composed of R images, and there are K sets, so a total of K × R images must be compared. . Further, it is not always obvious which part is nourished by which artery, and which map V _k (k = 1, 2,..., K) is observed for each part using anatomical knowledge. You must decide what to do. In particular, in cases where cerebrovascular disorders such as cerebral infarction have occurred, the arteries that control the tissue do not necessarily match the anatomical knowledge, and abnormal control is often observed. Due to these problems, there is a problem that it is difficult to interpret the vector value map.

To solve this problem, map composition is performed. That is, the K vector value maps V _k (k = 1, 2,..., K) are aggregated into one vector value map V using the residual map. For example, if P _{k, R} (x, y, z) is a residual map,
V (x, y, z) = V _k (x, y, z) where k is the minimum of | P _{k, R} (x, y, z) | Let k be such that

In addition, a map for indicating which of k = 1, 2,...
P ₀ (x, y, z) = (k = 1,2,..., K, where | P _{k, R} (x, y, z) | is the smallest)
Can also be added.

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

According to this method, there is a possibility that a calculation result using _Ak is erroneously used at a site (x, y, z) that should be calculated using the time density curve A _j of the cerebral artery. However, in order to make such an error, as is apparent from the definition of V (x, y, z), | P _{k, R} (x, y, z) | <| P _{j, R} (x, y, z) ｜
This relationship only occurs when A _j and A _k are very similar. For this reason, V _k (x, y, z) and V _j (x, y, z) are considered to be similar to each other at the site (x, y, z). There is almost no possibility of errors in the interpretation of y, z).

When this method is actually applied, P ₀ (x, y, z) = k for each site in the tissue that seems to be almost uniform only when A _j and A _k are very similar. Or P ₀ (x, y, z) = j, and P _{k, r} (x, y, z) ≈P _{j, r} (x, y, z) (r = 1) , 2,..., R), and it is observed that the result is almost the same no matter which one is employed.

Conversely, a tissue that are dominated by a specific artery, in a case where time-density curve A _k of the cerebral artery corresponding to it is not similar to other curve, by using this method, of the tissue In the region (x, y, z), V _k (x, y, z) is selected almost certainly and automatically. Therefore, by observing the above P ₀ (x, y, z), it is possible to observe which tissue is controlled by which artery without anatomical knowledge.

  Returning now to FIG. As shown in FIGS. 21A to 21D, a plurality of pixels are included for the CBP map, the CBV map, the MTT map, and the Err map synthesized as described above or each cerebral artery alone. A region of interest ROI is set (S18), and an average value (CBP average value, CBV average value, MTT average value, Err average value) of pixel values (CBP value, CBV value, MTT value, Err value) in the ROI is set. The calculated value (S19) may be used as a diagnostic material. At the time of this averaging, an appropriate range is set for each of CBP, CBV, MTT, and Err in step S14, and a value within the range is maintained. A value outside the range is, for example, a minimum corresponding to black expression. Since the value is replaced with a value, if the average including the replaced value is averaged, the average value naturally includes an error. Therefore, in this averaging process, it is necessary to perform the averaging process by selecting only values within the appropriate range or excluding the replaced values.

  Further, in setting the region of interest ROI for averaging, if the region of interest ROI is set on any of the CBP map, CBV map, MTT map, and Err map, the ROI becomes another map. Are also used in common, simplifying ROI setting work, and enabling calculation of average values (CBP average value, CBV average value, MTT average value, Err average value) related to the same ROI Yes.

  Here, as described above, the minimum residual Err of the time concentration curve of a certain tissue with respect to the time concentration curve of a certain cerebral artery is the degree to which the cerebral artery dominates the tissue, that is, the cerebral artery has its It shows how much blood flow is supplied to the tissue. In other words, it represents how much the tissue is subordinate to the cerebral artery, that is, how much the brain tissue is supplied with blood flow from the cerebral artery. The cerebral artery corresponding to the small residual Err indicates that the pixel is likely to dominate the brain tissue, and the cerebral artery corresponding to the large residual Err is less likely to dominate the brain tissue of the pixel. Represents. Therefore, a map in which cerebral arteries having a high dominance for each pixel from the residual Err are distinguished by a label, that is, a region having a high dominance of the anterior cerebral artery ACA, a region having a high dominance of the middle cerebral artery MCA, It is possible to generate a dominance map that distinguishes the region of the posterior cerebral artery PCA that is highly likely to be dominated.

  As shown in FIG. 22, for each of the left and right areas of the brain, 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 compared for each pixel. This indicates that the cerebral artery (ACA, MCA or PCA) showing the smallest residual is most likely to dominate the brain tissue of the pixel. For each pixel, the cerebral artery that has the highest control possibility, that is, the smallest residual Err is specified. A label corresponding to the identified cerebral artery is given to each pixel.

  FIG. 23 shows an example of the generated dominance map. This dominance map is displayed by distinguishing labels by color and shading. Further, by filtering the index map with an arbitrary label, as shown in FIGS. 24 (a), 24 (b), and 24 (c), each cerebral artery (ACA, MCA, or PCA) is obtained from the index map. It is possible to extract a region having a high dominance.

  As described above, according to the present embodiment, by using a coherent filter or coherent regression, it is possible to suppress a reduction in spatial and temporal resolution and suppress noise, thereby improving the analysis accuracy of a CBP study.

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

Explanatory drawing of the principle of CBP study. The block diagram which shows the structure of the index calculating apparatus regarding the blood flow dynamics of the capillary in the brain tissue which concerns on embodiment of this invention. Explanatory drawing of the image process by the coherent filter of this embodiment. The flowchart which shows the flow of the noise suppression process by the coherent filter in this embodiment. The flowchart of the first half part of the whole index calculation process in this embodiment. The flowchart of the second half part of the whole index calculation process in this embodiment. The figure which shows an example of the dividing line of step S3 of FIG. The figure which shows the influence area | region of each cerebral artery on the Voronoi diagram used in order to reduce the process man-hour of step S12 of FIG. FIG. 6 is a diagram showing an AT map, a PT map, and a TT map in step S4 of FIG. FIG. 6 is a diagram showing AT, PT, and TT in step S4 of FIG. 5. FIG. 6 is a diagram showing a cerebral artery ROI that is common between slices in step S6 of FIG. 5. The figure which shows the upper sagittal sinus ROI set by step S7 of FIG. FIG. 6 is a supplementary diagram regarding correction of the cerebral artery time density curve in step S10 of FIG. 5; The figure which shows an example of the time density | concentration curve Ca (t) of the cerebral artery created in step S10, S11 of FIG. 5, and the time density | concentration curve Ci (t) of brain tissue. FIG. 7 is a diagram illustrating the principle of the box-MTF method in step S12 of FIG. Explanatory drawing of the box-MTF process of step S12 of FIG. The figure which shows an example of the output range setting screen of each index of step S14 of FIG. The figure which shows an example of each map of CBP, CBV, MTT, and Err created by step S16 of FIG. FIG. 7 is a diagram showing a list of CBP maps, CBV maps, MTT maps, and Err maps created for each cerebral artery in step S16 of FIG. The figure for demonstrating the map composition method of step S17 of FIG. The figure which shows the example of a display of the average value calculated by step S19 of FIG. The figure which shows the production | generation method of the map (domination map) which distinguished the cerebral artery (ACA, MCA, PCA) with high possibility of control for every pixel (local tissue) by the label in step S17 of FIG. The figure which shows the example of the domination map produced | generated by the production | generation method of FIG. The figure which shows the CBP map filtered using the dominance map.

Explanation of symbols

10 ... Gantry part,
20 ... Computer device,
30. Image processing device,
101 ... X-ray tube,
101a ... high voltage generator,
102 ... X-ray detector,
103 ... data collection unit,
104 ... Pre-processing unit,
105: Memory unit,
106: 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 ... Weight 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 concentration curve correction unit,
124 ... MTF processing unit,
125: Index calculation unit,
126 ... map creation part,
127: Map composition unit.

Claims (4)

Creating a first time concentration curve related to an artery in the specific region and a second time concentration curve related to a tissue in the specific region from a plurality of continuous images related to the specific region of the subject injected with the contrast agent;
By selecting one of the first time concentration curves that best fits to each of the second time concentration curves, one of the arteries most likely to be subject to local blood flow dynamics of each of the tissues is identified. ,
A method for creating a map for distinguishing a dependent region of an artery based on one of the arteries specified for each tissue. 2. The method of claim 1, wherein the map is displayed with the arterial dependent regions distinguished by color. 2. The method according to claim 1, wherein an index map relating to an index representing the local blood flow dynamics corresponding to the artery is synthesized based on the map. Time to create a first time concentration curve related to an artery in the specific region and a second time concentration curve related to a tissue in the specific region from a plurality of continuous images related to the specific region of the subject injected with the contrast agent A concentration curve generator,
Means for identifying one of the arteries most likely to be subordinate to each tissue by selecting one of the first time concentration curves that best fits to each of the second time concentration curves;
And a map creating unit that creates a map for distinguishing a dependent region of the artery based on one of the arteries specified for each tissue.

Images