WO2013094186A1 - 動き追従x線ct画像処理方法および動き追従x線ct画像処理装置 - Google Patents
動き追従x線ct画像処理方法および動き追従x線ct画像処理装置 Download PDFInfo
- Publication number
- WO2013094186A1 WO2013094186A1 PCT/JP2012/008088 JP2012008088W WO2013094186A1 WO 2013094186 A1 WO2013094186 A1 WO 2013094186A1 JP 2012008088 W JP2012008088 W JP 2012008088W WO 2013094186 A1 WO2013094186 A1 WO 2013094186A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- ray
- image
- probability
- estimation
- distribution
- Prior art date
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 34
- 238000012545 processing Methods 0.000 title claims abstract description 23
- 230000033001 locomotion Effects 0.000 claims abstract description 133
- 238000010521 absorption reaction Methods 0.000 claims abstract description 45
- 238000005259 measurement Methods 0.000 claims abstract description 31
- 238000009826 distribution Methods 0.000 claims description 195
- 238000000034 method Methods 0.000 claims description 42
- 230000002123 temporal effect Effects 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 14
- 230000001419 dependent effect Effects 0.000 abstract description 2
- 230000006698 induction Effects 0.000 abstract 1
- 210000001519 tissue Anatomy 0.000 description 97
- 210000002216 heart Anatomy 0.000 description 20
- 238000012360 testing method Methods 0.000 description 16
- 230000008520 organization Effects 0.000 description 15
- 239000013598 vector Substances 0.000 description 15
- 238000004364 calculation method Methods 0.000 description 12
- 230000006870 function Effects 0.000 description 9
- 238000005457 optimization Methods 0.000 description 9
- 230000000694 effects Effects 0.000 description 8
- 238000003384 imaging method Methods 0.000 description 8
- 238000007796 conventional method Methods 0.000 description 7
- 210000000988 bone and bone Anatomy 0.000 description 6
- 230000006872 improvement Effects 0.000 description 6
- 210000003205 muscle Anatomy 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 230000000747 cardiac effect Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 230000009467 reduction Effects 0.000 description 4
- 238000007476 Maximum Likelihood Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 208000019622 heart disease Diseases 0.000 description 3
- 239000000126 substance Substances 0.000 description 3
- 210000004204 blood vessel Anatomy 0.000 description 2
- 230000001066 destructive effect Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 238000007689 inspection Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 239000002184 metal Substances 0.000 description 2
- 229910052751 metal Inorganic materials 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000029058 respiratory gaseous exchange Effects 0.000 description 2
- 238000013179 statistical model Methods 0.000 description 2
- 206010002383 Angina Pectoris Diseases 0.000 description 1
- 210000001015 abdomen Anatomy 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 230000001364 causal effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000002939 conjugate gradient method Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000003925 fat Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 239000007943 implant Substances 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 210000004072 lung Anatomy 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 150000002739 metals Chemical class 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 208000010125 myocardial infarction Diseases 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/277—Analysis of motion involving stochastic approaches, e.g. using Kalman filters
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
- A61B6/503—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of the heart
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20076—Probabilistic image processing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30048—Heart; Cardiac
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/412—Dynamic
Definitions
- the present invention relates to an X-ray CT algorithm capable of following the movement of an object to be measured, a program having the algorithm, and an X-ray CT apparatus equipped with the program.
- X-ray CT is widely used in medical applications and industrial applications such as defect inspection as a powerful technique capable of knowing the internal structure in a non-invasive and non-destructive manner.
- This reconstruction of the internal structure (CT image) by X-ray CT has generally been performed on the assumption that the object whose internal structure is to be known is stationary.
- heart diseases such as myocardial infarction and angina have been the second leading cause of death in the past 10 years, and heart diseases are the leading cause of death in the United States. Since one of the most effective countermeasures is the early detection of abnormalities in heart disease, accurate CT image reconstruction especially in cardiac X-ray CT is required.
- a technique for obtaining a reconstructed image with the same S / N ratio at a lower X-ray exposure as in the prior art is a technique for improving the resolution so that a high-resolution image can be obtained even with the same X-ray intensity and the number of images to be taken.
- the technology to shorten the imaging time by reducing the number of images to be taken and the exposure dose to reduce the exposure dose by reducing the X-ray intensity even with the same resolution as before
- Technology to reduce noise contained in the reconstructed image even with the same X-ray intensity and resolution as before, and the contrast resolution of the reconstructed image with the same X-ray intensity and resolution as before There is a technique for improving the contrast resolution that can be expressed in a tone, that is, the accuracy of the X-ray absorption coefficient.
- a filtered back projection method (FBP: Filtered Back Projection) that is a mainstream of the image reconstruction method.
- the FBP irradiates the imaging target with X-rays from various directions, obtains a projection image (a set of X-ray absorption images and projection images is called a sinogram), and then converts the obtained projection image in the reverse direction.
- This is an image reconstruction method in which a blur is removed by filtering the sinogram on the frequency domain when reconstructing an image by projecting the image.
- a projection image of the same cardiac phase is acquired at the timing when an electrocardiogram can be taken, and based on an X-ray CT image reconstructed from the acquired projection image,
- the length of one cycle of the heart is allowed to change, but it is assumed that the same volume change is performed every time within one cycle, at any time. Identify the cardiac phase.
- it is actually difficult to adjust the phase accurately in a fast moving part and there is a problem that the heart does not move in the same way every time.
- the present invention provides a motion tracking X-ray CT image processing method and motion tracking X-ray CT image processing that can reduce motion artifacts by simultaneously estimating motion and CT images as compared to image reconstruction using the prior art.
- An object is to provide an apparatus.
- the motion follow-up X-ray CT image processing method of the present invention is an X-ray CT image processing method for performing statistical estimation with prior knowledge about the motion of an object to be measured and the X-ray absorption coefficient, Assuming that the object to be measured changes smoothly with time, the first probability model for motion (prior knowledge of the object to be measured at all times) and the second probability model for observation (of the projection image at all times) (Observation model) and performing statistical estimation depending on both the first probability model and the second probability model.
- the first probability model (prior knowledge of the object to be measured at all times) regarding the movement of the object to be measured and the second probability model for observation (all By defining a time projection image observation model) and making inferences based on those probability models, motion and CT images can be estimated simultaneously without requiring synchronization of an electrocardiogram or the like.
- the second probability model for observation (all By defining a time projection image observation model) and making inferences based on those probability models, motion and CT images can be estimated simultaneously without requiring synchronization of an electrocardiogram or the like.
- inference methods established in the field of statistics can be applied, and the convergence of the algorithm is guaranteed, and excellent optimization is easy to perform.
- the motion is estimated only from the obtained CT reconstruction image depending on the model related to the motion (the conventional method is a deterministic model), and only the observation model is obtained from the obtained motion.
- the present invention estimates both the CT reconstructed image and the motion at the same time in consideration of both the probabilistic model for motion and the stochastic model for observation. This makes it possible to estimate with higher accuracy.
- a probability model that performs statistical estimation is used, other prior knowledge can be easily incorporated, and the quality of image reconstruction can be improved.
- the first probability model (prior knowledge of the measurement target object at all times) is a parameter defined for each pixel of the reconstructed image, and expresses the degree to which the reconstructed image is temporally continuous. It is expressed by a first probability distribution characterized by parameters, and a second probability model (an observation model of a projected image at all times) assumes that the reconstructed image is already known. It is preferable that the projection image is expressed by a second probability distribution describing what kind of projection image is expected to be observed.
- the first probability model different reconstructed images are assumed at different times and are different from those represented by the first probability distribution characterized by the parameters representing the degree of temporal succession of the reconstructed images. It is expressed in the form of a probability distribution that a reconstructed image of time is obtained by being deformed smoothly in time.
- the first probability distribution is further characterized by adding a term expressing a spatial correlation between adjacent pixels of the reconstructed image.
- the probability distribution representing that the value of each pixel of the reconstructed image is likely to take a specific value depending on the tissue and the probability distribution representing that the tissue is similarly displaced continuously in time, Further characterizing the first probability distribution.
- the probability distribution when the reconstructed image at each time is assumed to be known, information on what projection image is expected to be observed from the reconstructed image is the probability distribution. Expressed in shape.
- the probability distribution and the tissue expressing that the value of each pixel of the reconstructed image is likely to take a specific value depending on the tissue are similarly continuous in time.
- the first probability distribution described above is further characterized by a probability distribution representing displacement.
- the tissue class is expressed by a hidden state variable given to each pixel of the reconstructed image, and the first probability distribution is further characterized. That is, the first probability distribution is further characterized by a random variable representing the tissue class to further characterize the reconstructed image.
- a random variable representing a tissue class is determined for each pixel of a reconstructed image, for example, what type of tissue (normal cells such as muscle, soft cells such as fat, bones, etc.) Express whether it has a linear absorption coefficient.
- tissue normal cells such as muscle, soft cells such as fat, bones, etc.
- the representation is a representation of the probability distribution that assumes that the reconstructed image at each time has a direct correlation in time.
- it is possible to perform indirect expression through the organization class by assuming that the organization class at each time is smoothly displaced in time.
- distribution information on how much of these tissues are distributed and how easily each of the tissues is spatially continuously distributed is expressed in the form of a probability distribution.
- Prior knowledge about tissue distribution may use certain fixed average parameters, but it is possible to recreate the image by changing it appropriately according to individual differences such as physique, past history, gender, and age, and the imaging site. It is possible to further improve the accuracy of configuration and tissue estimation.
- the first probability distribution is further characterized by a hidden state variable expressing a tissue class belonging to each pixel. Due to the hidden state variable that controls the distribution of the tissue class belonging to each pixel, the assumed value of the X-ray absorption coefficient of the tissue class assumed in advance is greatly deviated from the X-ray absorption coefficient of the actual tissue class or is not assumed Can be made less affected by the presence of an organizational class. In addition, by statistically estimating this hidden state variable, it is possible to control the degree of influence of how much prior knowledge that a tissue-dependent CT value exists affects the reconstructed image that is the estimated solution. Can do.
- the motion tracking X-ray CT image processing method of the present invention it is preferable to use a Boltzmann distribution as the prior distribution of the tissue class.
- the reason why the Boltzmann distribution is used as the prior distribution of each tissue is that the measured object changes smoothly with time, so that the tissue changes smoothly with time even between adjacent frames.
- the assumption is easy to express by the temporal correlation of the organization class between adjacent frames.
- there are characteristics that can be expressed at the same time as the spatial continuity of the organization which is a property to be satisfied by the organization, and which organization is likely to exist at what ratio.
- Parameters representing the strength of time correlation, spatial correlation, and standard existence ratio can be adjusted by manual adjustment by experts for each organizational class or by learning based on how well it is applied to actual projection data.
- a Gaussian distribution can be suitably used as a parameter distribution that expresses a typical X-ray absorption coefficient value taken by each tissue, and a parameter distribution that expresses the degree of spatial continuity of each tissue.
- a Boltzmann distribution can be preferably used.
- a mixed Gaussian distribution that is a weighted sum of the Gaussian distribution can be suitably used as the parameter distribution that expresses the X-ray absorption coefficient of each tissue.
- the statistical estimation in the motion tracking X-ray CT image processing method of the present invention described above is specifically an estimation based on a value that maximizes the posterior probability (MAP estimation) or a Bayes estimation based on an expected value of the posterior probability.
- the MAP estimation estimates the most likely combination in the posterior distribution related to both the tissue to which each pixel region of the reconstructed image belongs and the X-ray absorption coefficient.
- Estimating the combination with the tissue to which the region of each pixel of the reconstructed image belongs includes the tissue class dependency.
- the tissue class of the air can include knowledge depending on the tissue class such that pixels are easily connected to each other (no other tissue enters the air). .
- the integral calculation necessary for obtaining the posterior probability of the reconstructed image can be omitted, so that the optimization can be accelerated and the reconstructed image can be obtained at a higher speed.
- Bayesian estimation uses a test distribution that can approximate the posterior distribution well in estimating the posterior distribution for both the X-ray absorption coefficient and the tissue to which each pixel region of the reconstructed image belongs.
- Bayesian estimation as in the case of MAP estimation, organization class dependency can be included. Thereby, distribution expression with more flexibility becomes possible.
- Bayesian estimation it is possible to make an inference that takes into account the uncertainty of estimation, it is easy to learn parameters, and X-ray CT image processing can be performed more accurately.
- the motion tracking X-ray CT image processing method of the present invention includes the following steps 1) to 5).
- a measurement condition including a projection image and at least an X-ray intensity is input.
- a physical process related to a photon observed in the projection image is defined by a probability distribution.
- a probabilistic model (prior knowledge of the object to be measured at all times) is characterized by parameters that define the degree to which each tissue continues in time, assuming that the object to be measured changes smoothly with time.
- Step 5 Information about what projection image is expected to be observed from the image is defined in the form of a second probability distribution. Step 5) Both the first probability model and the second probability model described above Exist to, steps of image reconstruction and tissue class estimation is performed by the expected value estimation of the posterior probability (Bayesian estimation) or maximize estimated posterior probability (MAP estimation)
- the probability distribution of the physical process related to photons in 2) above, specifically, the probability model related to photons observed in the projected image is expressed by Poisson distribution or the like.
- a physical process closer to the actual observation process can be expressed, and the uncertainty of observation due to photon noise or the like can be expressed.
- the description of 3) to 5) is the same as the above description.
- the motion-following X-ray CT image processing apparatus of the present invention is an X-ray CT image processing apparatus that performs statistical estimation with prior knowledge about the motion of an object to be measured and an X-ray absorption coefficient, and the object to be measured is time
- the first probability model preliminary knowledge of the object to be measured at all times
- the second probability model observation model of the projection image at all times
- the first probability model (prior knowledge of the measurement target object at all times) is a parameter defined for each pixel of the reconstructed image, and expresses the degree to which the reconstructed image is temporally continuous. It is expressed by a first probability distribution characterized by parameters, and a second probability model (an observation model of a projected image at all times) assumes that the reconstructed image is already known. It is preferable that the projection image is expressed by a second probability distribution describing what kind of projection image is expected to be observed.
- the first probability model different reconstructed images are assumed at different times and are different from those represented by the first probability distribution characterized by the parameters representing the degree of temporal succession of the reconstructed images. It is expressed in the form of a probability distribution that a reconstructed image of time is obtained by being deformed smoothly in time.
- the first probability distribution is further characterized by adding a term expressing a spatial correlation between adjacent pixels of the reconstructed image.
- the probability distribution representing that the value of each pixel of the reconstructed image is likely to take a specific value depending on the tissue and the probability distribution representing that the tissue is similarly displaced continuously in time, Further characterizing the first probability distribution.
- the probability distribution when the reconstructed image at each time is assumed to be known, information on what projection image is expected to be observed from the reconstructed image is the probability distribution. Expressed in shape.
- the value of each pixel of the reconstructed image takes a specific value depending on the tissue, as in the motion tracking X-ray CT image processing method of the present invention described above. It is preferable that the first probability distribution is further characterized by a probability distribution expressing that it is easy and a probability distribution indicating that the tissue is similarly displaced continuously in time. More preferably, in the motion following X-ray CT image processing apparatus of the present invention, the tissue class is expressed by a hidden state variable given to each pixel of the reconstructed image, and the first probability distribution is further characterized. That is, the first probability distribution is further characterized by a random variable representing the tissue class to further characterize the reconstructed image.
- a random variable representing a tissue class is determined for each pixel of a reconstructed image, for example, what type of tissue (normal cells such as muscle, soft cells such as fat, bones, etc.) Express whether it has a linear absorption coefficient.
- tissue normal cells such as muscle, soft cells such as fat, bones, etc.
- the representation is a representation of the probability distribution that assumes that the reconstructed image at each time has a direct correlation in time.
- it is possible to perform indirect expression through the organization class by assuming that the organization class at each time is smoothly displaced in time.
- distribution information on how much of these tissues are distributed and how easily each of the tissues is spatially continuously distributed is expressed in the form of a probability distribution.
- Prior knowledge about tissue distribution may use certain fixed average parameters, but it is possible to recreate the image by changing it appropriately according to individual differences such as physique, past history, gender, and age, and the imaging site. It is possible to further improve the accuracy of configuration and tissue estimation.
- the motion following X-ray CT image processing apparatus of the present invention includes the following means 1) to 5). 1) Means for inputting a measurement condition including a projection image and at least an X-ray intensity 2) Means for defining a physical process related to photons observed in the projection image by a probability distribution 3) A first probability model relating to the movement of an object to be measured A first probability characterized by a parameter that defines the degree of time continuity, assuming that the subject to be measured changes smoothly over time Means for defining distribution 4)
- the second probability model observation model of the projection image at all times regarding observation assumes that the reconstructed image is already known, what projection image is obtained from the reconstructed image?
- the probability distribution of the physical process related to photons in 2) above, specifically, a probability model related to photons observed in the projected image is expressed by Poisson distribution or the like.
- a physical process closer to the actual observation process can be expressed, and the uncertainty of observation due to photon noise or the like can be expressed.
- the motion tracking X-ray CT image processing program of the present invention is an X-ray CT image processing program for performing statistical estimation with prior knowledge about the movement of an object to be measured and the X-ray absorption coefficient.
- Steps 1) to 5) are executed.
- a measurement condition including a projection image and at least an X-ray intensity is input 2
- a physical process relating to photons observed in the projection image is represented by a probability distribution 3)
- the present invention can provide an X-ray CT image processing apparatus equipped with the above X-ray CT image processing program.
- motion and CT images are estimated at the same time, and motion artifacts are reduced as compared with image reconstruction using a conventional technique. It has the effect of being able to.
- the estimation method is based on the assumption that the measurement target object is deformed smoothly in time, there is no need to synchronize with an electrocardiogram, etc.
- image reconstruction is performed using statistical estimation with prior knowledge regarding motion and X-ray absorption coefficients.
- a normal CT model and a CT model assuming motion will be described with reference to FIG.
- FIG. 1 (1) in the case of the conventional filter corrected back projection (FBP) CT model, it is assumed that the object is stationary even if it is an object to be measured such as a moving heart. Thought that the projection could be obtained from one object (CT image).
- FIG. 1 (2) it is assumed that the object to be measured is moving from the beginning, and the projection is obtained from these different objects.
- step S11 the projection image and measurement conditions such as the X-ray intensity are input (step S11).
- step S12 the prior knowledge regarding the X-ray absorption coefficient is expressed on the assumption that the measurement target object changes smoothly with time, and the degree to which the reconstructed image continues in time.
- step S13 Expressed by the first probability model (prior knowledge of the object to be measured at all times) and the second probability model (observation model of the projection image at all times) characterized by the parameters to be performed (step S13), and the posterior probability In the expected value estimation (Bayesian estimation), motion and tissue class are estimated and a plane image is reconstructed (step S14).
- the above-described first state is represented by a hidden state variable representing a tissue class belonging to each pixel.
- the probability distribution of 1 is further characterized (step S24). Since image reconstruction has defect setting properties, it is solved by imposing some restrictions on the reconstructed image (X-ray absorption coefficient) as a solution using appropriate prior knowledge. Of these constraints, the key to the present invention is the time constraint imposed on the reconstructed image between adjacent frames. Since it is assumed that the non-measuring object is deformed smoothly in time, the reconstructed image between adjacent frames as shown in FIG.
- the reconstructed image in the ( k ⁇ 1) th frame deformed according to the motion m j (k) should have a strong temporal correlation with the reconstructed image in the kth frame.
- the temporal constraints on the reconstructed images are both the direct time correlation between the reconstructed images and the indirect time correlation via the material / tissue that is assumed to exist potentially behind the reconstructed image described later. Or by one of them in the form of a probability distribution.
- Constraints through materials / organizations are as follows. As shown in FIG. 5, the human body is composed of limited substances / tissues such as fat, water, blood, muscles, bones, and metals (for example, implants), and each substance / tissue (hereinafter simply referred to as tissue).
- tissue a substance / tissue
- the restriction that it is easy to take a typical X-ray absorption coefficient can be characterized in the form of a probability distribution.
- the deformation of the tissue along the movement of the measurement target can be expressed in a form similar to the temporal correlation of the reconstructed image between adjacent frames, and this constraint is characterized by a probability distribution.
- each organization has spatial continuity, and knowledge that the proportion of each organization in the human body is roughly known is similarly characterized in the form of a probability distribution as a constraint on the organization.
- FIG. 6 specifically shows the causal relationship between each variable in FIG. 1 (2) according to the probability model of the present invention.
- the scan time is divided into frames having a sufficiently short time width, and the object is represented by the same CT image in each frame (in FIG. 6, this is represented by x (1) , x (2) ,. As such).
- x (1) , x (2) In an actual scan, hundreds to thousands of projections are obtained in a time as short as 0.3 seconds, so that a plurality of projections y are obtained from x of each frame.
- B is a hidden state variable that controls the shape of the prior distribution, and makes the prior distribution more robust.
- FIG. 7 is a diagram showing the effect of introducing the hidden state variable B. If there is no hidden state variable B that controls the shape of the prior distribution, for example, if the variance of the prior distribution is reduced to reduce the ambiguity, the effect of prior knowledge also becomes stronger depending on it. By introducing B, the prior distribution can be expressed as a mixed distribution, and the ambiguity of prior knowledge and the strength of the effect can be controlled separately.
- Image reconstruction can be performed by obtaining a value (MAP estimation) or expected value (Bayesian estimation) that maximizes the posterior distribution for the reconstructed image.
- MAP estimation a value
- Bayesian estimation that calculates the expected value of the posterior probability is more robust to disturbances, but the calculation of the posterior probability involves integration calculations for high-dimensional hidden state variables, and is difficult to perform analytically. Therefore, this difficulty of calculation is overcome by applying an approximation method that approximately obtains the posterior probability.
- this calculation is difficult by applying MAP estimation by maximizing simultaneous posterior probabilities including both hidden state variables and reconstructed images without performing integral calculation on hidden state variables. You can overcome this.
- each projection is divided into K frames, and a projection image is obtained from the same X-ray CT image within each frame.
- a set of projection data indexes belonging to the k-th frame is represented as S (k), and an X-ray CT image in the k-th frame is represented as X (k) . Since X-rays are exponentially attenuated when passing through a substance, the average y i (t) of X-ray photons expected to be observed can be expressed as Equation 1 below.
- v represents the number of photons emitted from the X-ray source (number of photons that can be observed when no object is placed).
- l ij (t) is proportional to the distance at which the projection line detected by the i th detector and the j th pixel intersect when projected from the angle ⁇ (t)
- l ij (t ) x j represents the effective X-ray absorption coefficient of the j th pixel area with respect to the X-rays incident on t th i th detector projection.
- a physical model related to photons observed in a projection image a prior distribution related to tissue as prior knowledge, that is, an imaging target Parameter distribution expressing the existence ratio of each tissue, parameter distribution expressing typical X-ray absorption coefficient that each tissue can easily take, parameter distribution expressing the degree of spatial continuity of each tissue, time of each tissue
- an imaging target Parameter distribution expressing the existence ratio of each tissue a parameter distribution expressing typical X-ray absorption coefficient that each tissue can easily take
- parameter distribution expressing the degree of spatial continuity of each tissue time of each tissue
- Equation 2 is expressed as a product of independent Poisson distribution as a whole.
- S (k) represents a set of projection indices included in the kth frame.
- Equation 3 (Prior distribution regarding movement) Assuming that the distribution of the tissue class under the motion ⁇ is expressed by the Boltzmann distribution, the prior distribution related to the tissue class is modeled as shown in the following Equation 3.
- a 1 is a normalization constant
- H 1 (Z) is an energy function representing what combination of tissue classes is likely.
- the energy function H 1 (Z) is expressed by the following mathematical formula 4.
- z j (k) [z j1 (k) ,..., Z jC (k) ] ⁇ ⁇ 0, 1 ⁇ C represents an organization class to which the j-th pixel belongs in the k-th frame.
- Z jc (k) 1 if the vector belongs to the organization class c, and 0 for the other components.
- the uppercase Z represents a set of lowercase variables z j (k) .
- Equation 4 ⁇ (j) is 4 neighborhoods of the jth pixel, and N (j) is the surrounding 9 ⁇ 9 (pixels) including the jth pixel.
- 1 item is the tissue ratio
- 2 item is the spatial smoothness (ease of connection) of the tissue in the same frame
- 3 item is the smooth temporal displacement / deformation of the tissue between adjacent frames.
- G (i, j, m (k) ( ⁇ )) is the weight of the three items H 1 (Z), and is expressed by Equation 5 below.
- each pixel in the ( k ⁇ 1) -th frame moves according to the current motion m (k) ( ⁇ ), and the distance to each pixel in the k-th frame is calculated, and the distance The closer the is, the greater the weight.
- Equation 6 (A prior distribution regarding the reconstructed image X) The prior distribution of the reconstructed image X is shown in Equation 6 below.
- the X-ray absorption coefficient is expressed as a Gaussian distribution.
- the average of the Gaussian distribution is set for each tissue class.
- the degree of approach to the value of the typical X-ray absorption coefficient in the tissue class is controlled by the hidden state variable b.
- b j1 (k) and b j2 (k) are variables in which one takes 1 and one takes 0.
- the dispersion parameters ⁇ c1 2 and ⁇ c2 2 are fixed values for each tissue class c, but this is a parameter in which one has a large positive value and the other has a small positive value.
- Equation 6 determines the degree of proximity of the X-ray absorption coefficient represented by the reconstructed image to a typical X-ray absorption coefficient value.
- Equation 6 are terms representing that the reconstructed image changes smoothly in space and that the image is smoothly deformed in time based on the motion vector m.
- G x spatial and G x teemporal Is a parameter that controls the strength of spatial correlation and the strength of temporal correlation.
- the hidden state variable B follows a Boltzmann distribution (the following equation 7)
- the energy function H 2 (B) representing what combination of the hidden state variables B is likely is expressed by the following equation 8.
- one item of Expression 8 represents which of binary values B is a set of binary hidden variables, and two items are spatial smoothness of distribution (same in adjacent pixels). (Ease of distribution) is controlled.
- r i represents the position vector of the i-th pixel.
- the X-ray absorption coefficient x, the tissue class z, and the hidden state variable b are expressed as the posterior distribution p (X, Z, B
- D, ⁇ ) can be decomposed by Bayes' theorem, and is expressed as the following Equation 10 in a form proportional to the product of each prior distribution and likelihood.
- Equation 10 in a form proportional to the product of each prior distribution and likelihood.
- test distribution q (X, Z, B) that approximates a posterior distribution p (X, Z, B
- Bayesian estimation itself is a method of calculating the posterior distribution and using the expected value as an estimated value. If the posterior distribution and the expected value can be calculated analytically, the test distribution is not necessary.
- the test distribution q (X, Z, B) is determined so as to minimize an index representing the closeness of the distribution called KL distance.
- KL distance here is defined by Equation 11 below.
- ⁇ •> q (X, Z, B) is an operator representing performing integral calculation with respect to the test distribution q (X, Z, B). That is, it means that an expected value is taken with respect to the test distribution q (X, Z, B).
- a search is made for a set of computable distributions that is closest to the posterior distribution.
- the test distribution q (X, Z, B) is arbitrarily selected as long as the KL distance (the above formula 11) can be minimized or approximated. be able to.
- the test distribution q (X, Z, B) has a variable called a factorization assumption with respect to the test distribution q (X, Z, B) as shown in Equation 12 below.
- a factorization assumption with respect to the test distribution q (X, Z, B) as shown in Equation 12 below.
- X, Z, and B are independent, and further independent for each frame and pixel.
- x j on the distribution q (x j) is a Gaussian distribution, modeled as shown in following equation 13.
- the optimal test distribution is obtained by using alternate optimization that alternately optimizes q (X), q (Z), and q (B).
- the image shown in FIG. 8 is an artificially generated deformable object (64 ⁇ 64 pixels), which simulates a moving object.
- the center part of each circular screen represents an object, which expands with time (a) ⁇ (b) ⁇ (c) ⁇ (d) ⁇ (e). It was.
- the deformation of the object is assumed to expand at different ratios in the vertical and horizontal directions.
- the result of reconstruction of the X-ray CT image is shown in FIG. FIG.
- PSNR Peak signal to noise ratio
- Equation 14 MSE represents the mean square error between the pixel values of the true image and the reconstructed image, and max i y i represents the maximum pixel value of the true image.
- the PSNR value is appended to the bottom of each figure in FIG. From the results of FIG. 9, it can be seen that the present invention provides an average improvement of about 3.8 (dB) in PSNR. An improvement of about 3.8 (dB) means a reduction of about 42% in terms of mean square error. It can be seen from the obtained image that the movement of the object whose center deforms can be well expressed.
- the image shown in FIG. 10 is an image of a chest model with heart motion, and simulates the chest with heart motion.
- the ellipse on the center corresponding to the heart is gradually increased, that is, (a) ⁇ (b) ⁇ (c) ⁇ (d) ⁇ (e) corresponding to the heart with time.
- the ellipse was assumed to expand.
- the deformation of the ellipse was assumed to expand at different ratios in the vertical and horizontal directions.
- the image size of the chest model was 256 ⁇ 256 (pixels), and the number of projections was 360 between 360 °.
- the result of reconstruction of the X-ray CT image is shown in FIG. FIG.
- FIG. 11 shows the result of image reconstruction by the filtered back projection method (FBP), which is most commonly used as a conventional method that does not assume motion, and the result of image reconstruction (proposed method) according to the present invention. Is a comparison. From the results of FIG. 11, it can be seen that the present invention provides an average improvement in PSNR of about 1.1 (dB). An improvement of about 1.1 (dB) means that the mean square error is reduced by about 22%. It can be seen from the obtained image that the heart motion can be well expressed.
- FBP filtered back projection method
- the motion tracking X-ray CT image processing method of the present invention uses a stochastic model that performs statistical estimation, so other prior knowledge can be easily incorporated, and the quality of image reconstruction is improved. it can.
- the estimation accuracy can be improved by performing estimation taking into account the representative X-ray absorption coefficient of each tissue.
- the estimation accuracy can be improved and robust estimation can be performed by including a hidden state variable for controlling the distribution of each tissue class. An experiment for examining the effect will be described below.
- FIG. 12 shows the X-ray absorption coefficient of the assumed tissue class and the value of the X-ray absorption coefficient of the object actually used in the experiment.
- the X-ray absorption coefficient of the fat and bone class is 10 as shown by the line a in the figure. %.
- the X-ray absorption coefficient was assumed to be high.
- the line indicated by symbol b in the figure is the value of the X-ray absorption coefficient taken by other outliers.
- X-ray CT was performed on an object having a value 5% larger than a typical X-ray absorption coefficient assumed for fat and muscle.
- FIG. 13 compares the results of X-ray CT reconstruction under the above-described condition setting.
- (a) is a true image
- (b) is a reconstructed image by FBP
- (c) is a reconstructed image according to the present invention, and no hidden state variable B is introduced
- (D) is a reconstructed image according to the present invention with the introduction of a hidden state variable B.
- (e) represents the estimated hidden state variable B.
- the PSNR was improved by 0.8 (dB) (the mean square error was improved by about 17%).
- the improvement is 17.6 (dB).
- the case where the X-ray absorption coefficient in each tissue deviates from the assumption is considered, only a small shift of 5% is assumed. For this reason, it is expected that a greater effect can be obtained when a larger shift occurs.
- Y (t) [y 1 (t) ,..., Y 1 (t) ] (where I is the number of detectors) and summarize all time projections.
- D ⁇ Y (1) ,..., Y (T) ⁇ .
- L (t ) in the equation (16) is an mxn matrix
- the (i, j) element of the irradiated X-rays is the ray j incident on the detector i and the pixel j of the reconstructed image. It only shows how to intersect.
- b in the equation (16) is the number of photons (number of irradiated photons) at the time of blank scanning.
- a Gaussian distribution can be adopted as the noise distribution.
- Variable B ⁇ B (1) ,..., B (K) ⁇ for controlling the distribution shape.
- the variable b j (k) that controls the distribution shape of the pixel j in the k-th frame is a two-dimensional vector of the following Expression 19. In order to express a more complex distribution, it is possible to increase the number of values that d can take to an arbitrary number larger than two.
- the parameter ⁇ is defined as a parameter representing motion, and the prior distribution is expressed by the following Equation 20.
- E (X, Z, B, ⁇ ) in Expression 20 is expressed by Expression 21 below.
- G x, cd spatial it is possible to control how smooth the pixel value x should be when the tissue class of adjacent pixels and the distribution shape control variable are the same. Further, by changing G x, c spatial , it is possible to control how smooth the pixel value x should be when the tissue classes of adjacent pixels are the same regardless of the distribution shape control variable. Similarly, G z, c spatial controls how easily adjacent pixels are in the same tissue class, and G b, cd spatial controls how much the same distribution shape control variable is taken when adjacent pixels are in the same tissue class. G cd self represents what combination of the tissue class variable and the distribution shape control variable is likely to occur.
- G x temporal is a variable that controls whether G z temporal
- Each variable x, z how much temporal easy transitions along the motion parameter theta.
- ⁇ (j) and N (j) represent a set of pixel indexes in the vicinity of the pixel j. Each of them determines how much the smoothness is affected, and it is preferable to make it small in order to save the calculation amount.
- the function G (i, j, m (k) , ( ⁇ )) is defined by the following formula 24.
- r i is a position vector of the pixel i
- m i (k) ( ⁇ ) is a motion vector representing a movement amount of the pixel i from the frame k to k + 1. That is, it expresses that the pixel i of the k-th frame is expected to move from the position r i to the position r i + m i (k) ( ⁇ ) in the (k + 1) -th frame.
- ⁇ is a certain positive constant.
- the function G (i, j, m (k) , ( ⁇ )) has a larger value for the function G (i, j, m (k) , ( ⁇ )) as the pixel is closer to the expected movement amount. It works as a weighting function.
- the motion vector m i (k) representing the movement amount is defined by the following Equation 25. If exp ( ⁇
- Equation 27 p (Z, B
- E (X, Z, B, ⁇ ) in the above equation 26 is defined by the following equation 30.
- Equation 31 u cd and v cd are parameters of the gamma distribution.
- Equation 33 The log likelihood logp (y
- Equation 37 the logarithmic likelihood can be expressed as the following Equation 37 using free energy.
- the left side in Expression 37 is an amount that does not depend on the test distribution q (x)
- the right side also does not depend on how the test distribution is taken.
- Equation 39 holds. From this equation, the right side is the minimum value ⁇ logp (when KL [q (x)
- Equation 43 the independence represented by the following Equation 43 was assumed. Considering that optimization is performed by alternating optimization under this assumption, the following Expression 44 is obtained. And the solution of Numerical formula 45 can be analytically calculated
- Equation 49 can be evaluated analytically. Further, in Equation 49, if the parameter (q (x j (k) ) of q (X) follows a Gaussian distribution, the parameter ( ⁇ cd , ⁇ cd ) assumes an average ⁇ and variance ⁇ 2 , and a gamma distribution. Can be optimized by a nonlinear optimization method such as a conjugate gradient method, etc. Note that the gamma distribution has a preferable property that can express the non-negativeness of x.
- the present invention is useful for medical and industrial X-ray CT apparatuses.
- Medical images, especially X-ray CT of the heart, is an important application, but even in industrial non-destructive inspection, when measuring minute objects, the minute movement of the stage becomes a problem or fluid is present. Since deformation of fluid becomes a problem, it is expected that the reduction of motion artifacts of the present invention is useful in various applications.
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Medical Informatics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Health & Medical Sciences (AREA)
- Radiology & Medical Imaging (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Animal Behavior & Ethology (AREA)
- High Energy & Nuclear Physics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Biophysics (AREA)
- Multimedia (AREA)
- Pulmonology (AREA)
- Quality & Reliability (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
Description
日本では、ここ10年ほど、心筋梗塞、狭心症などの心疾患が死因の2位を占めており、アメリカでは心疾患が死因の1位であるなど心疾患は国内外で対応が求められている疾患であり、その最も有効な対応策の一つが心疾患の異常の早期発見にあるとされるため、特に心臓のX線CTにおける正確なCT画像再構成が求められている。
統計推定を行う確率モデルに基づくことにより、統計の分野で確立された推論手法が適用でき、またアルゴリズムの収束性が保証されるなど、優れた最適化が行いやすい。
すなわち、従来の技術においては、得られたCT再構成画像のみから動きに関するモデル(従来法は決定論的モデル)に依存して動きを推定しており、また得られた動きから観測モデルのみに依存してCT再構成画像を推定しているのに対して、本発明では、動きに関する確率的モデルおよび観測に関する確率モデルの両方を考慮して、CT再構成画像と動きの両方を同時に推定することでより精度の高い推定が行える。
また、統計推定を行う確率モデルを用いているために、他の事前知識も容易に組み込むことができ、画像再構築の品質を向上できる。
ここで、好ましくは、再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第1の確率分布をさらに特徴付ける。
また、好ましくは、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付ける。
さらに好ましくは、本発明の動き追従X線CT画像処理方法において、再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスを表現し、上記の第1の確率分布をさらに特徴付ける。すなわち、再構成画像をさらに特徴付けるために組織クラスを表す確率変数によって、第1の確率分布をさらに特徴付ける。組織クラスを表す確率変数は、再構成画像の画素毎に、例えば、人体などの撮像対象に関し、どのような組織(筋肉などの通常細胞、脂肪などの軟細胞、骨など)がどのようなX線吸収係数を有するかを表現する。再構成画像は時間的に滑らかに変形することが第1の確率分布に表現されるが、その表現は各時刻の再構成画像が時間的に直接相関をもつことを仮定する確率分布による表現の他、各時刻の組織クラスが時間的に滑らかに変位すること仮定することで組織クラスを通じた間接的な表現を行うことができる。また、これらの組織がどの程度の割合で分布するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいか、の分布情報が確率分布の形で表現される。
なお、組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることができる。
画素毎に所属する組織クラスの分布を制御する隠れ状態変数により、事前に想定した組織クラスのX線吸収係数の想定値が実際の組織クラスのX線吸収係数から大きくズレたり、想定していなかった組織クラスが存在していた場合に影響を受け難くすることができる。また、この隠れ状態変数を統計推定するにより、組織依存のCT値が存在するという事前知識がどの程度、推定される解である再構成画像へ影響するか、という影響度の度合いを制御することができる。
各組織の事前分布としてボルツマン分布を用いることとした理由は、被計測対象物が時間と共に滑らかに変化することを仮定したことから、隣接するフレーム間でも組織は時間と共に滑らかに変化することを仮定するべきであるが、その仮定を隣接するフレーム間の組織クラスの時間相関によっての表現しやすいことに注目したものである。
その他、組織が満たすべき性質である、組織の空間的連続性や、どの組織がどの程度の割合で存在しやすいかも同時に表現できる特性をもつ。 時間相関や空間相関の強さ、標準的な存在割合を表すパラメータは、組織クラスごとに専門家による人手による調整や実際の投影データへの当てはまりの良さを基準とした学習による調整が可能である。
なお、上記の各組織のとる典型的なX線吸収係数の値を表現するパラメータ分布としてはガウス分布を好適に用いることができ、また、各組織の空間的に連続する度合いを表現するパラメータ分布としてはボルツマン分布を好適に用いることができる。組織クラスの分布を制御する隠れ状態変数を導入した際には、各組織のX線吸収係数を表現するパラメータ分布はガウス分布の重み付き和である混合ガウス分布を好適に用いることができる。
ここで、MAP推定は、再構成画像の各画素の領域が属する組織とX線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定する。再構成画像の各画素の領域が属する組織との組合せを推定することは、組織クラス依存性を含めることになる。この組織クラス依存性を含めることで、例えば、空気の組織クラスはピクセル同士の連結がしやすい(空気中に他の組織は入らない)などの組織クラスに依存した知識を含めることができるのである。これにより、より現実の状況に沿った分布表現が可能となり結果的に解の推定精度を高めることが期待される。MAP推定によれば、再構成画像の事後確率を求めるのに必要な積分計算が省略することができるため、最適化を速めることができ、より高速で再構成画像を求めることができる。
1)投影画像と少なくともX線強度を含む計測条件が入力されるステップ
2)投影画像に観測されるフォトンに関する物理過程が確率分布で定義されるステップ
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、各組織が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第1の確率分布の形で定義されるステップ
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについての情報が第2の確率分布の形で定義されるステップ
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定が行われるステップ
上記3)~5)についての説明は、上述の説明と同様である。
本発明の動き追従X線CT画像処理装置は、被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理装置であって、被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第2の確率モデル(全時刻の投影像の観測モデル)を定義するモデル定義手段と、上記の第1の確率モデルと第2の確率モデルの両方に依存した統計推定を行う統計推定手段を備える。
ここで、好ましくは、再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第1の確率分布をさらに特徴付ける。
また、好ましくは、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付ける。
さらに好ましくは、本発明の動き追従X線CT画像処理装置において、再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスを表現し、上記の第1の確率分布をさらに特徴付ける。すなわち、再構成画像をさらに特徴付けるために組織クラスを表す確率変数によって、第1の確率分布をさらに特徴付ける。組織クラスを表す確率変数は、再構成画像の画素毎に、例えば、人体などの撮像対象に関し、どのような組織(筋肉などの通常細胞、脂肪などの軟細胞、骨など)がどのようなX線吸収係数を有するかを表現する。再構成画像は時間的に滑らかに変形することが第1の確率分布に表現されるが、その表現は各時刻の再構成画像が時間的に直接相関をもつことを仮定する確率分布による表現の他、各時刻の組織クラスが時間的に滑らかに変位すること仮定することで組織クラスを通じた間接的な表現を行うことができる。また、これらの組織がどの程度の割合で分布するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいか、の分布情報が確率分布の形で表現される。
なお、組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることができる。
1)投影画像と少なくともX線強度を含む計測条件を入力する手段
2)投影画像に観測されるフォトンに関する物理過程を確率分布で定義する手段
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、その時間的に連続する度合いを定義するパラメータによって特徴付けられる第1の確率分布を定義する手段
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布を定義する手段
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定を行う手段
1)投影画像と少なくともX線強度を含む計測条件を入力するステップ
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップ
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)が、各時刻の再構成画像が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第1の確率分布で定義するステップ
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時にその再構成画像から、どのような投影像が観測されると期待されるかについて記述した第2の確率分布で定義するステップ
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップ
また、被計測対象物が時間的に滑らかに変形するという仮定に基づいた推定方法であるため、心電図などによる同期を取る必要がなく、心電図などの同期が取れない予測不可能な体動や内臓運動に対応が可能な他、全時刻の投影像を利用した推定が行えるため被曝量の低減が期待される。
図1を参照し、通常のCTモデルと動きを想定したCTモデルの説明を行う。図1(1)に示すように、従来のフィルター補正逆投影法(FBP)のCTモデルの場合、動きのある心臓などの被計測対象物であっても、対象は静止していると仮定して1つの対象(CT像)から投影が得られると考えていた。本発明の場合、図1(2)に示すように、被計測対象物は初めから動いているものと仮定し、投影はこれらの異なる対象から得られるものと考える。
被計測対象物は時間と共に滑らかに変化することを仮定するとは、CT画像の画素毎に動き(変形)パラメータを設け、画素毎に隣接フレーム間での動き(変形)を推定することである。具体的には、図2に示すように、j番目の画素の異なる時刻間の再構成画像(k-1番目のフレームとk番目のフレーム)の間での、動きをベクトルmj (k)で表す。
以下では、時間と共に滑らかに変化する被計測対象物に対して投影を行っていることを統計モデルによって記述し、この統計モデルに基づきベイズ推定を行うことで画像再構成できることを説明する。
本発明の動き追従X線CT画像処理方法のフローの一実施形態としては、図3のフローに示すように、投影画像とX線強度などの計測条件を入力し(ステップS11)、投影画像の観測フォトンに関して物理モデルで表現し(ステップS12)、X線吸収係数に関する事前知識を、被計測対象物が時間と共に滑らかに変化することを仮定し、再構成画像が時間的に連続する度合いを表現するパラメータで特徴付けた第1の確率モデル(全時刻の被計測対象物の事前知識)と第2の確率モデル(全時刻の投影像の観測モデル)とで表現し(ステップS13)、事後確率の期待値推定(ベイズ推定)で、動き,組織クラスを推定および面像再構成する(ステップS14)。
画像再構成は不良設定性を有しているため、適当な事前知識を用いて解である再構成画像(X線吸収係数)に何らかの制約を課すことで解決する。この制約のうち、本発明の肝となるのが隣接するフレーム間の再構成画像に与えられる時間的制約である。非計測対象物が時間的に滑らかに変形することを想定しているため、図2のとおり隣接するフレーム間の再構成画像も各画素jごとに動きmj (k)に応じた動きを生じる。したがって、k-1番目のフレームにおける再構成画像を動きmj (k)にしたがって変形させたものは、k番目のフレームにおける再構成画像と強い時間相関をもつべきである。
この再構成画像に対する時間的制約は、再構成画像間の直接の時間相関と、後述の再構成画像の背後に潜在的に存在すると仮定される物質・組織を介した間接的な時間相関の両方、もしくは一方によって確率分布の形で表現される。
先ず、X線CT画像のように、様々な方向から投影されたT個の投影データを、D={Y(1),・・・,Y(T)}と表すことにする。各々のデータY(t)は、t番目の投影によって検出器で検出されるデータの集合であり、Y(t)={y1 (1),・・・,yI (t)}となる。yi (t)は、i番目の検出器で検出された光子の数を表す。
短時間に十分、多くの投影が得られているときには、その短時間内では、X線CT画像が同一であると想定したほうが推定すべきX線CT画像を減らすことができ計算量を抑えることができる。そこで、各投影をK通りのフレームに分け、各フレーム内では同一のX線CT画像から投影像が得られると仮定する。k番目のフレームに所属する投影データのインデックスの集合をS(k)と表し、k番目のフレームにおけるX線CT画像をX(k)と表すことにする。X線は物質を透過する際に指数的に減衰することから、観測されると期待されるX線光子の平均yi (t)は、下記数式1のように表すことができる。
X線CTにおける観測データの主なノイズの要因は、検出される光子(フォトン)のショットノイズであると考えられる。そこで、光子(フォトン)の観測データは、投影毎、検出器毎に独立なポアソン分布に従って生成されるとしてモデル化する。
なお、ポアソン分布に特に限定されるものではなく、物理過程をより良く表現できる他の物理モデルも適用可能である。
上述の数式1は、ポアソン分布の平均を表しており、数式1の結果を用いて、下記数式2で表わす。ノイズは投影、画素(ピクセル)毎に独立であることから、数式2は全体として独立なポアソン分布の積で表されている。数式2において、S(k)は、k番目のフレームに含まれる投影のインデックスの集合を表している。
動きθが与えられたもとでの組織クラスの分布はボルツマン分布で表されると仮定して、組織クラスに関する事前分布を下記数式3に示すようにモデル化する。ここで、A1は正規化定数であり、H1(Z)はどのような組織クラスの組み合わせが尤もらしいか表すエネルギー関数である。エネルギー関数H1(Z)は、下記数式4で表される。
ここで、zj (k) = [zj1 (k), ・・・, zjC (k) ] ∈ {0,1}Cは、k番目のフレームにおいてj番目のピクセルが属する組織クラスを表すベクトルであり、組織クラスcに属する場合はzjc (k) =1、その他の成分は0となる。大文字のZは小文字の変数zj (k)の集合を表す。
再構成画像Xの事前分布を下記数6に示す。組織クラスZと隠れ状態変数Bが与えられたもとで、X線吸収係数はガウス分布として表現している。以下に説明するように、ガウス分布の平均は組織クラスごとに設定する。また、その組織クラスで典型的なX線吸収係数の値への接近の度合いは隠れ状態変数bによって制御される。各フレームk、各画素jにおいてbj1 (k), bj2 (k)はいずれか一方が1をとりいずれか一方が0をとる変数である。組織クラスcごとに与えられるその組織クラスcに典型的なX線吸収係数の値vc1, vc2は通常、同じ値vc1=vc2=vcが設定される。たとえば、平均vkの値は、各々の組織(空気,脂肪,筋肉,骨,金属)ごとにそれぞれ、v1=0.000,v2=0.018,v3=0.022,v4=0.040,v5=0.120と設定する。また、分散パラメータεc1 2, εc2 2は組織クラスcごとに固定した値であるが、これは一方が大きく、一方が小さい正の値をとるパラメータである。vc1, vc2,εc1 2, εc2 2はいずれも専門家による人手あるいは実際に得られた多くの投影像への当てはまりの良さを基準としたパラメータ学習によって与えられる。
数式6の1項目は典型的なX線吸収係数の値への再構成画像が表すX線吸収係数の近接度の度合いを決定する。隠れ状態変数bj1 (k), bj2 (k)はいずれか一方が1、いずれか一方が0をとる変数であるため、再構成画像の組織クラスcの典型的なX線吸収係数への値の近接度は組織クラスcのみによって決定されるのではなく各フレームk、各画素jごとにbj1 (k), bj2 (k)のとる値によっても制御されることになる。このような組織クラスcごとに固定された分散パラメータεc1 2, εc2 2の他に、与えられた投影像から推定される隠れ状態変数bj1 (k), bj2 (k)を用意し推定することで各フレーム、各画素ごとに適応的に近接度を制御可能とし、
これによって組織クラスごとにどのX線吸収係数をとりやすいかに関する不確かさとその影響の強さを別々に制御することができる。数式6の2項目と3項目はそれぞれ再構成画像が空間的に滑らかに変化すること、動きベクトルmに基づいて時間的に滑らかに変形することを表す項であり、Gx spatialとGx teemporalは空間相関の強さと時間相関の強さを制御するパラメータである。
CT画像の再構成では、X線吸収係数xと組織クラスzと隠れ状態変数bは、観測データDと被計測対象物の動きθの事後分布p(X,Z,B|D,θ)として推定される。事後分布p(X,Z,B|D,θ)は、ベイズの定理により分解でき、それぞれの事前分布と尤度の積に比例する形で、下記数式10のように表される。MAP推定においては事後分布p(X,Z,B|D,θ)の最大値を求めることにより、高次元の隠れ状態変数に関する和計算の必要がなくなる。
試験分布q(X,Z,B)は、KL距離と呼ばれる分布の近さを表す指標を最小化するように決定される。ここでのKL距離は下記数式11で定義される。
試験分布q(X,Z,B)は、KL距離(上記数式11)の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。そのため、計算を簡単にするために、試験分布q(X,Z,B)には、下記数式12に示すように、試験分布q(X,Z,B)に対して因子化仮定と呼ばれる変数xとzとbの間の独立性の仮定を行う。すなわち、X,Z,Bは独立で、さらにフレーム、ピクセル毎にも独立であるとする。 さらに、xjに関する分布q(xj)がガウス分布になるものと仮定して、下記数式13に示されるようにモデル化する。これらの仮定の下で、q(X),q(Z),q(B)を交互に最適化する交互最適化を用い最適試験分布を求めている。
q(Z)とq(B)を現在、得られている推定値に固定しても、q(X)は解析解が得られないので、上記数式13におけるq(X)のパラメータμ,σは数値最適化によって求める。これらの分布q(X), q(Z), q(B)は、事前に決められた所定の収束条件を満たすまで最適化が繰り返され、解が求められる。
先ず、図8に示す画像は、人工的に生成した変形する物体(64×64ピクセル)であり、動きのある物体を模擬するものである。図8に示すように、それぞれの円形の画面の中央部分が物体を表しており、これが時刻と共に(a)→(b)→(c)→(d)→(e)と膨張していくものとした。物体の変形は、縦と横にそれぞれ異なる比率で膨張していくものとした。
X線CT画像の再構成の結果を図9に示す。図9は、真の画像と、動きを想定しない従来法として最もよく用いられているフィルター補正逆投影法(FBP)による画像再構成の結果と、本発明による画像再構成(提案法)の結果とを比較したものである。
それぞれの手法による再構成画像と真の画像との差をPSNR(Peak signal to noise ratio)によって評価した。PSNRは下記数式14のように定義されるものである。
図9の結果から、本発明はPSNRで平均して約3.8(dB)の改善をもたらしていることがわかる。約3.8(dB)の改善は平均二乗誤差でいうと約42%の誤差を低減していることを意味する。得られる画像からも中心の変形する物体の動きをよく表現できていることがわかる。
X線CT画像の再構成の結果を図11に示す。図11は、真の画像と、動きを想定しない従来法として最もよく用いられているフィルター補正逆投影法(FBP)による画像再構成の結果と、本発明による画像再構成(提案法)の結果とを比較したものである。
図11の結果から、本発明はPSNRで平均して約1.1(dB)の改善をもたらしていることがわかる。約1.1(dB)の改善は平均二乗誤差でいうと約22%の誤差を低減していることを意味する。得られる画像からも心臓の動きをよく表現できていることがわかる。
実際の投影では、想定した代表値からX線吸収係数の値のズレが生じることを考えて、図中に符号aで示した線のように、脂肪と骨のクラスのX線吸収係数は10%だけX線吸収係数が高いものとした。また図中に符号bで示した線は、その他の外れ値のとるX線吸収係数の値である。各組織クラスのうち、脂肪と筋肉に関して想定した代表的なX線吸収係数から5%大きな値をとる物体に対してX線CTを行った。
隠れ状態変数Bを導入した場合と導入しない場合とで、PSNRで0.8(dB)の改善(平均二乗誤差で17%程度の改善)となった。従来法であるフィルター補正逆投影法(FBP)と比較すると17.6(dB)の改善となっている。また、各組織におけるX線吸収係数が想定からズレた場合を考えているが、5%と小さいシフトしか想定していない。このため、より大きなズレが生じる場合には、より大きな効果が得られると期待される。
(1)確率モデルについて
(1-a)尤度
再構成の対象となる物体のCT画像を1次元ベクトルX=[x1,・・・,xJ](ここで,Jは再構成画像の画素数)で表現する。 ここでは、時間に応じて変形もしくは移動する物体に対する投影を考えるため、時間をK個のフレームS(k), k∈{1,・・・,K}に分割し、k番目のフレームでのCT画像をX(k)=[x1 (k),・・・,xJ (k)]で表現する。全フレームのX(k)(k∈{1,・・・,K})を、纏めてX={X(1),・・・,X(k)}で表現する。
時刻tにおける物体に対する投影による像をY(t)=[y1 (t),・・・,y1 (t)] (ここで、Iは検出器数)とし、全時刻の投影を、纏めてD={Y(1),・・・,Y(T)}で表現する。時刻tがt∈S(k)を満たす場合、CT画像X(k)の投影を得たと想定するので、光電変換で生じるノイズが独立なポアソン分布に従うことを考えたとき、下記数式15,数式16で表現できる。
事前分布を下記数式17で表現する。
分布形状を制御する変数B={B(1),・・・,B(K)}である。さらに、k番目のフレームにおける組織クラス, 分布形状を制御する変数は、それぞれZ(k)=[z1 (k),・・・,zJ (k)],B(k)=[b1 (k),・・・,bJ (k)]と表現される。想定する組織クラスの数をCとしたとき、各画素の組織クラスは、C次元ベクトルzj (k)=[zj1 (k),・・・,zjC (k)]で表現される。zjC (k)は、そのk番目のフレームの画素jがとる組織クラスcに対してのみzjC (k)=1をとり、その他の要素はゼロをとる。
したがって、下記数式18の制約を満足する。また、k番目のフレームの画素jの分布形状を制御する変数bj (k)は、は、下記数式19の2次元ベクトルとしている。なお、より複雑な分布を表現するためには、dの取りうる値の数を2より大きな任意の数に増やすことで対応可能である。
数式20におけるE(X,Z,B,θ)は、下記数式21で表現される。
まず、下記数式22で表現されるパラメータは、zjc (k)=zsc´ (k)=bjd (k)=bsd´ (k)=1のときのみ、値Gx,cd spatialをとる。また、下記数式23で表現されるパラメータは、zjc (k)=zsc´ (k)=1のときのみ、値Gx,c spatialをとる。
(2)自由エネルギー最小化に基づく推定方法について
表記を簡単にするために観測変数yから隠れ変数xを推定する一般的な問題を考える。上記の確率モデルでは、観測変数をD,隠れ変数をX,Z,Bと考えれば対応がつく。また、自由エネルギーFは、下記数式32で定義される。数式32において、q(x)は試験分布と呼ばれる分布である。
ここでは、数式41を求めるために、数式42に示す2つのステップを繰り返す。
試験分布q(x)の最適値は、事後分布p(x|y,θ)で与えられるが、一般に事後分布の解析的な計算は困難である。そこで、q(x)の候補を制約することで最適化を可能にする。
そして、変分法によって、数式45の解を解析的に求めることができる。一方、数式46に関しては解析的に求まらないので、q(xj (k))として、下記数式47のガウス分布を仮定するか、或いは、下記数式48のガンマ分布を仮定する。
Claims (18)
- 被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理方法であって、
被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第2の確率モデル(全時刻の投影像の観測モデル)を定義するステップと、
上記の第1の確率モデルと第2の確率モデルの両方に依存した統計推定を行うステップと、
を備えることを特徴とする動き追従X線CT画像処理方法。 - 1)前記第1の確率モデルが、
再構成画像の画素毎に定義されるパラメータであって、
再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第1の確率分布で表現され、
2)前記第2の確率モデルが、
再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布で表現される、
ことを特徴とする請求項1に記載の動き追従X線CT画像処理方法。 - 再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第1の確率分布をさらに特徴付ける、ことを特徴とする請求項2に記載の動き追従X線CT画像処理方法。
- 再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付ける、ことを特徴とする請求項2に記載の動き追従X線CT画像処理方法。
- 再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスの事前分布を制御し、上記の第1の確率分布をさらに特徴付ける、ことを特徴とする請求項4に記載の動き追従X線CT画像処理方法。
- 前記時間的に連続する度合いを表現するパラメータの分布としてボルツマン分布を用いる、ことを特徴とする請求項2に記載の動き追従X線CT画像処理方法。
- 前記統計推定は、事後確率の最大化による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定で行う、ことを特徴とする請求項1~6のいずれかに記載の動き追従X線CT画像処理方法。
- 1)投影画像と少なくともX線強度を含む計測条件が入力されるステップと、
2)前記投影画像に観測されるフォトンに関する物理過程が確率分布で定義されるステップと、
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、その時間的に連続する度合いが定義されるパラメータによって特徴付けられる第1の確率分布で定義されるステップと、
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布で定義されるステップと、
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定が行われるステップと、
を備えたことを特徴とする動き追跡X線CT画像処理方法。 - 被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理装置であって、
被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第2の確率モデル(全時刻の投影像の観測モデル)を定義するモデル定義手段と、
上記の第1の確率モデルと第2の確率モデルの両方に依存した統計推定を行う統計推定手段と、
を備えることを特徴とする動き追従X線CT画像処理装置。 - 1)前記第1の確率モデルが、
再構成画像の画素毎に定義されるパラメータであって、
再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第1の確率分布で表現され、
2)前記第2の確率モデルが、
再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布で表現される、
ことを特徴とする請求項9に記載の動き追従X線CT画像処理装置。 - 再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第1の確率分布をさらに特徴付ける、ことを特徴とする請求項9に記載の動き追従X線CT画像処理装置。
- 再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付ける、ことを特徴とする請求項9に記載の動き追従X線CT画像処理装置。
- 再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスの事前分布を制御し、上記の第1の確率分布をさらに特徴付ける、ことを特徴とする請求項12に記載の動き追従X線CT画像処理装置。
- 前記時間的に連続する度合いを表現するパラメータの分布としてボルツマン分布を用いる、ことを特徴とする請求項9に記載の動き追従X線CT画像処理装置。
- 前記統計推定手段は、事後確率の最大化による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定で行う、ことを特徴とする請求項9~14のいずれかに記載の動き追従X線CT画像処理装置。
- 1)投影画像と少なくともX線強度を含む計測条件を入力する手段と、
2)前記投影画像に観測されるフォトンに関する物理過程を確率分布で定義する手段と、
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、その時間的に連続する度合いを定義するパラメータによって特徴付けられる第1の確率分布を定義する手段と、
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布を定義する手段と、
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定を行う手段と、
を備えたことを特徴とする動き追跡X線CT画像処理装置。 - 被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理プログラムであって、
コンピュータに、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)が、各時刻の再構成画像が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第1の確率分布で定義するステップと、
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布で定義するステップと、
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップと、
を実行させるための動き追従X線CT画像処理プログラム。 - 請求項17の動き追従X線CT画像処理プログラムを搭載したX線CT画像処理装置。
Priority Applications (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US14/366,037 US9129431B2 (en) | 2011-12-18 | 2012-12-18 | Motion-tracking X-ray CT image processing method and motion-tracking X-ray CT image processing device |
EP12859055.1A EP2792303A4 (en) | 2011-12-18 | 2012-12-18 | MOTION-TRACKING X-RAY CT IMAGE PROCESSING AND MOTION-TRACKING X-RAY CT IMAGE PROCESSING DEVICE |
JP2013550113A JP6044046B2 (ja) | 2011-12-18 | 2012-12-18 | 動き追従x線ct画像処理方法および動き追従x線ct画像処理装置 |
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2011-276545 | 2011-12-18 | ||
JP2011276545 | 2011-12-18 |
Publications (1)
Publication Number | Publication Date |
---|---|
WO2013094186A1 true WO2013094186A1 (ja) | 2013-06-27 |
Family
ID=48668104
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
PCT/JP2012/008088 WO2013094186A1 (ja) | 2011-12-18 | 2012-12-18 | 動き追従x線ct画像処理方法および動き追従x線ct画像処理装置 |
Country Status (4)
Country | Link |
---|---|
US (1) | US9129431B2 (ja) |
EP (1) | EP2792303A4 (ja) |
JP (1) | JP6044046B2 (ja) |
WO (1) | WO2013094186A1 (ja) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2019212016A1 (ja) * | 2018-05-01 | 2019-11-07 | 国立大学法人東北大学 | 画像処理装置,画像処理方法および画像処理プログラム |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2014185078A1 (ja) | 2013-05-15 | 2014-11-20 | 国立大学法人京都大学 | X線ct画像処理方法,x線ct画像処理プログラム及びx線ct画像装置 |
US9495770B2 (en) * | 2013-08-14 | 2016-11-15 | University Of Utah Research Foundation | Practical model based CT construction |
US10026015B2 (en) * | 2014-04-01 | 2018-07-17 | Case Western Reserve University | Imaging control to facilitate tracking objects and/or perform real-time intervention |
US10085703B2 (en) * | 2015-01-27 | 2018-10-02 | Septimiu Edmund Salcudean | Dynamic computed tomography imaging of elasticity |
US9848844B2 (en) * | 2015-10-09 | 2017-12-26 | Carestream Health, Inc. | Iterative reconstruction process |
EP3710109B1 (en) * | 2016-11-21 | 2023-08-23 | Leo Cancer Care, Inc. | Three-dimensional tracking of a target in a body |
CN108257196A (zh) * | 2018-01-11 | 2018-07-06 | 苏州润心医疗器械有限公司 | 一种基于心脏ct图像的血管拉直重建方法 |
US11204409B2 (en) * | 2018-10-11 | 2021-12-21 | University Of Virginia Patent Foundation | Systems and methods for motion-compensated reconstruction of magnetic resonance images |
DE102019201079B4 (de) * | 2019-01-29 | 2021-06-17 | Siemens Healthcare Gmbh | Verfahren und bildgebendes Gerät zum Erzeugen eines bewegungskompensierten Bildes, Computerprogramm und Speichermedium |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2005532137A (ja) * | 2002-07-10 | 2005-10-27 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | 画像の情報内容を向上する方法及びシステム |
US20070242796A1 (en) * | 2004-07-19 | 2007-10-18 | Deutsches Krebsforschungszentrum | Method for Producing X-Ray Computer Tomography Images From Limited Data of an Image Object |
JP2009500115A (ja) * | 2005-07-08 | 2009-01-08 | ウイスコンシン アラムナイ リサーチ フオンデーシヨン | Ct画像の逆投影再構成法 |
JP2011156302A (ja) * | 2010-02-03 | 2011-08-18 | Kyoto Univ | X線ct画像処理方法,x線ctプログラムおよび該プログラムが搭載されたx線ct装置 |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7099435B2 (en) * | 2003-11-15 | 2006-08-29 | Agilent Technologies, Inc | Highly constrained tomography for automated inspection of area arrays |
US7215732B2 (en) * | 2004-09-30 | 2007-05-08 | General Electric Company | Method and system for CT reconstruction with pre-correction |
US7881505B2 (en) * | 2006-09-29 | 2011-02-01 | Pittsburgh Pattern Recognition, Inc. | Video retrieval system for human face content |
WO2008144908A1 (en) * | 2007-05-29 | 2008-12-04 | Christopher Mott | Methods and systems for circadian physiology predictions |
WO2011002874A1 (en) | 2009-06-30 | 2011-01-06 | University Of Utah Research Foundation | Image reconstruction incorporating organ motion |
DE102010019016B4 (de) * | 2010-05-03 | 2017-03-02 | Siemens Healthcare Gmbh | Verfahren zur Rekonstruktion von Bilddaten eines bewegten Untersuchungsobjektes aus Messdaten nebst zugehöriger Gegenstände |
RU2013136488A (ru) * | 2011-01-05 | 2015-02-10 | Конинклейке Филипс Электроникс Н.В. | Способ и устройство детектирования и коррекции движения в данных позитронно-эмиссионной томографии в режиме списка с использованием синхронизированного сигнала |
-
2012
- 2012-12-18 JP JP2013550113A patent/JP6044046B2/ja not_active Expired - Fee Related
- 2012-12-18 US US14/366,037 patent/US9129431B2/en not_active Expired - Fee Related
- 2012-12-18 WO PCT/JP2012/008088 patent/WO2013094186A1/ja active Application Filing
- 2012-12-18 EP EP12859055.1A patent/EP2792303A4/en not_active Withdrawn
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2005532137A (ja) * | 2002-07-10 | 2005-10-27 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | 画像の情報内容を向上する方法及びシステム |
US20070242796A1 (en) * | 2004-07-19 | 2007-10-18 | Deutsches Krebsforschungszentrum | Method for Producing X-Ray Computer Tomography Images From Limited Data of an Image Object |
JP2009500115A (ja) * | 2005-07-08 | 2009-01-08 | ウイスコンシン アラムナイ リサーチ フオンデーシヨン | Ct画像の逆投影再構成法 |
JP2011156302A (ja) * | 2010-02-03 | 2011-08-18 | Kyoto Univ | X線ct画像処理方法,x線ctプログラムおよび該プログラムが搭載されたx線ct装置 |
Non-Patent Citations (13)
Title |
---|
AA. ISOLA ET AL.: "Fully automatic nonrigid registration-based local motion estimation for motion-corrected iterative cardiac CT reconstruction", MED PHYS., vol. 37, no. 3, 2010, pages 1093 - 1109 |
AA. ISOLA ET AL.: "Motion compensated iterative reconstruction of a region of interest in cardiac cone-beam CT", COMPUT. MED. IMAG. GRAP., vol. 34, no. 2, 2010, pages 149 - 159 |
B. OHNESORGE ET AL.: "Cardiac imaging by means of electrocardiographically gated multisection spiral CT: initial experience", RADIOLOGY, vol. 217, 2000, pages 564 - 571 |
BING FENG ET AL.: "Use of Three-Dimensional Gaussian Interpolation in the Projector/ Backprojector Pair of Iterative Reconstruction for Compensation of Known Rigid-Body Motion in SPECT", IEEE TRANSACTIONS ON MEDICAL IMAGING, vol. 25, no. 7, July 2006 (2006-07-01), pages 838 - 844, XP001545917 * |
C. C. MOREHOUSE ET AL.: "Gated cardiac computed tomography with a motion phantom", RADIOLOGY, vol. 134, 1980, pages 134 - 137 |
C. J. RITCHIE ET AL.: "Correction of computed tomography motion artifacts using pixel specific back-projection", IEEE TRANS. MED. IMAGING, vol. 15, 1996, pages 333 - 342 |
DANAI LAKSAMEETHANASAN ET AL.: "A Bayesian Reconstruction Method with Marginalized Uncertainty Model for Camera Motion in Microrotaion Imaging", IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, vol. 57, no. 7, July 2010 (2010-07-01), pages 1719 - 1728, XP011343285 * |
G. WANG ET AL.: "A knowledge-based cone-beam x-ray CT algorithm for dynamic volumetric cardiac imaging", MED PHYS., vol. 29, no. 8, 2002, pages 1807 - 1822 |
H. I. GOLDBERG ET AL.: "Evaluation of ultrafast CT scanning of the adult abdomen", INVEST. RADIOL., vol. 24, 1989, pages 537 - 543 |
K. TAGUCHI ET AL.: "Toward time resolved cardiac CT images with patient dose reduction: image-based motion estimation", NUCLEAR SCIENCE SYMPOSIUM CONFERENCE RECORD, vol. 4, 2006, pages 2029 - 2032 |
P. M. JOSEPH; J. WHITLEY: "Experimental simulation evaluation of ECG-gated heart scans with a small number of views", MED. PHYS., vol. 10, 1983, pages 444 - 449 |
See also references of EP2792303A4 |
WATARU FUKUDA ET AL.: "Bayesian X-ray Computed Tomography Using a Mixture Prior", TECHNICAL REPORT OF IEICE, NC2009-133, March 2010 (2010-03-01), pages 267 - 272, XP008173885 * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2019212016A1 (ja) * | 2018-05-01 | 2019-11-07 | 国立大学法人東北大学 | 画像処理装置,画像処理方法および画像処理プログラム |
JPWO2019212016A1 (ja) * | 2018-05-01 | 2021-05-13 | 国立大学法人東北大学 | 画像処理装置,画像処理方法および画像処理プログラム |
JP7273416B2 (ja) | 2018-05-01 | 2023-05-15 | 国立大学法人東北大学 | 画像処理装置,画像処理方法および画像処理プログラム |
Also Published As
Publication number | Publication date |
---|---|
EP2792303A1 (en) | 2014-10-22 |
US9129431B2 (en) | 2015-09-08 |
EP2792303A4 (en) | 2015-08-05 |
US20140334705A1 (en) | 2014-11-13 |
JP6044046B2 (ja) | 2016-12-14 |
JPWO2013094186A1 (ja) | 2015-04-27 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP6044046B2 (ja) | 動き追従x線ct画像処理方法および動き追従x線ct画像処理装置 | |
Würfl et al. | Deep learning computed tomography: Learning projection-domain weights from image domain in limited angle problems | |
US11126914B2 (en) | Image generation using machine learning | |
Arndt et al. | Deep learning CT image reconstruction in clinical practice | |
JP6855223B2 (ja) | 医用画像処理装置、x線コンピュータ断層撮像装置及び医用画像処理方法 | |
Liu et al. | Total variation-stokes strategy for sparse-view X-ray CT image reconstruction | |
Mohammadinejad et al. | CT noise-reduction methods for lower-dose scanning: strengths and weaknesses of iterative reconstruction algorithms and new techniques | |
EP3716214A1 (en) | Medical image processing apparatus and method for acquiring training images | |
Zhang et al. | Regularization strategies in statistical image reconstruction of low‐dose x‐ray CT: A review | |
Wang et al. | Deep learning-based image quality improvement for low-dose computed tomography simulation in radiation therapy | |
US20180018757A1 (en) | Transforming projection data in tomography by means of machine learning | |
JP2020168352A (ja) | 医用装置及びプログラム | |
EP2150918B1 (en) | Methods and systems for improving spatial and temporal resolution of computed images of moving objects | |
US20130051516A1 (en) | Noise suppression for low x-ray dose cone-beam image reconstruction | |
Isola et al. | Fully automatic nonrigid registration‐based local motion estimation for motion‐corrected iterative cardiac CT reconstruction | |
JP2020036877A (ja) | 反復的画像再構成フレームワーク | |
JP2020168353A (ja) | 医用装置及びプログラム | |
JP5590548B2 (ja) | X線ct画像処理方法,x線ctプログラムおよび該プログラムが搭載されたx線ct装置 | |
EP3047391B1 (en) | Method and system for statistical modeling of data using a quadratic likelihood functional | |
US10475215B2 (en) | CBCT image processing method | |
JP2021013725A (ja) | 医用装置 | |
CN111540025A (zh) | 预测用于图像处理的图像 | |
Gajera et al. | CT-scan denoising using a charbonnier loss generative adversarial network | |
Gao et al. | CoreDiff: Contextual error-modulated generalized diffusion model for low-dose CT denoising and generalization | |
Chen et al. | Low-dose CT image denoising model based on sparse representation by stationarily classified sub-dictionaries |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 12859055 Country of ref document: EP Kind code of ref document: A1 |
|
ENP | Entry into the national phase |
Ref document number: 2013550113 Country of ref document: JP Kind code of ref document: A |
|
NENP | Non-entry into the national phase |
Ref country code: DE |
|
WWE | Wipo information: entry into national phase |
Ref document number: 14366037 Country of ref document: US |
|
WWE | Wipo information: entry into national phase |
Ref document number: 2012859055 Country of ref document: EP |