WO2014185078A1 - Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置 - Google Patents

Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置 Download PDF

Info

Publication number
WO2014185078A1
WO2014185078A1 PCT/JP2014/002571 JP2014002571W WO2014185078A1 WO 2014185078 A1 WO2014185078 A1 WO 2014185078A1 JP 2014002571 W JP2014002571 W JP 2014002571W WO 2014185078 A1 WO2014185078 A1 WO 2014185078A1
Authority
WO
WIPO (PCT)
Prior art keywords
ray
substance
absorption coefficient
distribution
rays
Prior art date
Application number
PCT/JP2014/002571
Other languages
English (en)
French (fr)
Inventor
新一 前田
大悟 吉川
田中　匠
石井　信
Original Assignee
国立大学法人京都大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立大学法人京都大学 filed Critical 国立大学法人京都大学
Priority to EP14798390.2A priority Critical patent/EP2997900A4/en
Priority to US14/891,303 priority patent/US9662081B2/en
Priority to JP2015516922A priority patent/JPWO2014185078A1/ja
Publication of WO2014185078A1 publication Critical patent/WO2014185078A1/ja

Links

Images

Classifications

    • 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/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
    • 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
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • 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/42—Arrangements for detecting radiation specially adapted for radiation diagnosis
    • A61B6/4208—Arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a particular type of detector
    • A61B6/4241—Arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a particular type of detector using energy resolving detectors, e.g. photon counting
    • 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
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/416—Exact reconstruction

Definitions

  • the present invention relates to an X-ray CT (Computed Tomography) image processing method, an X-ray CT image processing program, and an X-ray CT image apparatus in consideration of the tissue structure in the body.
  • X-ray CT Computed Tomography
  • X-ray CT (Computed Tomography) is a technique for imaging the inside of an observation object as a distribution of X-ray absorption coefficients from projection data obtained by irradiating the observation object with X-rays.
  • the X-ray CT has a configuration in which an X-ray detector is disposed opposite to an X-ray tube and a turntable is disposed therebetween. The rotation axis of the turntable is set to be orthogonal to the X-ray optical axis connecting the X-ray tube and the X-ray detector.
  • transmission X-ray data from a plurality of angles is collected with respect to a measurement object, and a CT tomographic image is reconstructed.
  • transmission X-ray data is collected by rotating 180 ° or more.
  • X-ray CT reconstruction is a technique for estimating the X-ray absorption coefficient distribution of an internal substance by obtaining the degree of attenuation of X-rays transmitted through a measurement object from projection images obtained at various projection angles.
  • the superabsorbent material strongly attenuates X-rays, so that almost no transmitted X-rays are detected and the observation noise is lost.
  • the FBP method is susceptible to noise, so a virtual image called a metal artifact is generated in the reconstructed image, and an accurate image is generated.
  • the FBP method is easily affected by noise.
  • the main noise includes shot noise according to the Poisson distribution, the signal-to-noise ratio can be improved by increasing the intensity of irradiated X-rays. This can reduce the influence of relative noise and improve the accuracy of image reconstruction, but on the other hand, the radiation exposure dose increases, and a probabilistic health risk such as carcinogenesis occurs.
  • FBP filtered back projection
  • PCLIS Principal Component Interference Deformation Method based on Linear Interpolation in the Sinogram
  • MLE maximum likelihood estimation
  • MAP Maximum A Posteriori
  • Bayesian estimation Bayesian estimation.
  • Maximum likelihood estimation (MLE) is a statistical estimation method that does not have prior knowledge about the reconstructed image to be estimated. It represents physical processes including uncertainty in observation by a probability model, and based on the probability model. The most likely reconstructed image is estimated. The physical process including uncertainty in observation includes fluctuation (shot noise) related to the number of transmitted X-ray photons to be observed.
  • the MAP estimation method and the Bayesian estimation method express prior knowledge about an estimation target in the form of a probability model and incorporate it into estimation. That is, in the MAP estimation method and the Bayesian estimation method, estimation is performed in consideration of both the probability model related to observation and the probability model related to the estimation target.
  • an estimation target is an X-ray absorption coefficient.
  • the maximum value of the posterior distribution of the X-ray absorption coefficient calculated from the probability model related to observation and the probability model related to the X-ray absorption coefficient of the object to be estimated is used as an estimated value, and the Bayesian estimation method is used. Then, let the expected value of the posterior distribution of the same X-ray absorption coefficient be an estimated value.
  • the Bayesian estimation method takes the expected value of the posterior distribution, the estimation is performed by averaging the population of the X-ray absorption coefficient that increases the posterior distribution. Do. The point at which the posterior distribution is maximum can change sensitively due to observation noise, but the expected value of the posterior distribution is not so sensitive to observation noise, so the Bayesian estimation method is more stable and accurate. It becomes.
  • a single absorption coefficient is estimated on the assumption that an X-ray to be irradiated is a line spectrum in image reconstruction.
  • actual X-rays are not single energy, but are composed of an energy spectrum that is continuously distributed to the lower side with the maximum energy determined by the highest tube voltage applied to the X-ray tube. For this reason, when X-rays are irradiated to a substance, X-rays with high energy are transmitted almost without being absorbed, but X-rays with low energy are greatly attenuated. A phenomenon called beam hardening that shifts to the angle occurs.
  • the phenomenon caused by the energy dependence (wavelength dependence) of the substance's X-ray absorption coefficient assumes that the X-ray source is a line spectrum, and the energy dependence (wavelength dependence) of the substance's X-ray absorption coefficient It cannot be expressed by a model of an observation process that is not taken into account, and incorrect estimation is performed. In this way, beam hardening artifacts are generated by estimation using a model that does not consider the energy dependency (wavelength dependency) of the X-ray absorption coefficient of a substance.
  • streak artifacts and dark band artifacts in which the X-ray absorption coefficient around the superabsorbent object is estimated to be smaller than actual may occur.
  • this is largely due to the fact that the effects of beam hardening are not taken into account. This will be described below.
  • X-rays are greatly attenuated when passing through a highly absorbing object such as metal having a very high X-ray absorption coefficient.
  • the property is expressed in the observation model, even if it is a model that does not consider the energy dependence (wavelength dependence) of the X-ray absorption coefficient of the substance, an object having a large X-ray absorption coefficient in the vicinity of the highly absorbing object It can be estimated that it will exist.
  • high-energy X-rays are easily transmitted even with a highly absorbing object. Therefore, in practice, some X-rays are observed even if many of the photons having the main energy of the X-ray source are absorbed while passing through the highly absorbing object.
  • the X-rays transmitted through all of the superabsorbent objects are high in a model that does not consider the energy dependency (wavelength dependency) of the X-ray absorption coefficient of the substance. It is expected that transmitted X-rays will be further attenuated when objects having absorption coefficients overlap. However, in practice, a certain amount of high-energy X-rays are observed even if a plurality of highly-absorbing objects are transmitted. For this reason, it is a model that does not take into account the energy dependence (wavelength dependence) of the X-ray absorption coefficient of a substance.
  • a method for reducing metal artifacts a method of separating projection data into a metal region and a non-metal region, reconstructing an image in each region, and superimposing the obtained images is already known. By doing so, it is possible to correct the contradiction of projection data in the metal region and reduce the occurrence of metal artifacts. Further, the inventors of the present invention provide an X-ray CT image processing method for performing statistical estimation with prior knowledge about an X-ray absorption coefficient to be a reconstructed image, wherein the prior knowledge is stored in each pixel region of the reconstructed image.
  • an X-ray CT image processing method capable of reducing metal artifacts by using what is expressed by a probability distribution characterized by see Patent Document 1.
  • these methods do not consider the energy dependence (wavelength dependence) of the X-ray absorption coefficient, the actual observation process and the observation process that does not consider the energy dependence of the X-ray absorption coefficient of the substance. It has not yet been possible to eliminate the beam hardening artifact caused by the deviation.
  • Non-Patent Documents 1 to 4 models taking into account the energy dependence (wavelength dependence) of the X-ray absorption coefficient are known (Non-Patent Documents 1 to 4).
  • a model that considers the energy dependence (wavelength dependence) of the X-ray absorption coefficient it is necessary to estimate the X-ray absorption coefficient for each X-ray energy. This is based on a more accurate physical process, but as an estimation problem, it causes a problem of poor setting. For example, let us consider the problem of discretizing X-ray energy for simplicity and estimating the X-ray absorption coefficient for two types of energy.
  • Non-Patent Documents 1 to 3 the energy dependency (wavelength dependency) of a substance is expressed by expressing the energy-dependent absorption coefficient in each pixel by a weighted sum of absorption coefficient functions of known substances. Yes.
  • the energy dependency (wavelength dependency) of the X-ray absorption coefficient varies depending on the material, but there is a common property among the materials such as a tendency to attenuate as the energy increases. If the function that outputs the X-ray absorption coefficient at that energy is called the absorption coefficient function when the energy of the X-ray is used as an input variable, even if different substances have similar absorption coefficient functions, the absorption of any substance
  • the coefficient function can also be well expressed by a weighted sum (basis function expression) of a small number, typically two absorption coefficient functions. For example, in Non-Patent Document 2, the substance class is determined based on an image obtained by FBP, and image reconstruction is performed based on that knowledge.
  • Non-Patent Document 4 the X-ray absorption coefficient for each pixel is estimated on the assumption that it can be expressed by a combination of X-ray absorption coefficients for each known substance type.
  • the method considering the energy dependency (wavelength dependency) of the X-ray absorption coefficient is a basis function expression using the wavelength dependency of the X-ray absorption coefficient of a known substance in order to suppress the defect setting property. Has been used to reduce the degree of freedom of parameters to be estimated. On the other hand, however, problems that impose too strong restrictions have been pointed out.
  • Non-Patent Documents 5 to 8 the image reconstruction method using DECT generally has a problem that the exposure dose is larger than that of the monochromatic X-ray CT, and the estimation accuracy is deteriorated if the exposure dose is reduced to the same level. Further, Non-Patent Document 5 is a method of correcting by image processing afterwards, and there is a possibility that a further virtual image may be generated by the correction.
  • the observation model that does not consider the energy dependence of the X-ray absorption coefficient of the substance assumes that the X-rays incident on the detector are not composed of a plurality of wavelengths but only a certain representative wavelength. This can be considered as a simplified process in which multiple X-ray spectra are incident on the detector.
  • the phenomenon that X-rays are incident on the detector multiple times occurs not only in the spectral space but also in time and space. That is, multiple X-rays are incident on the detector within a certain detection time, and multiple X-rays are incident on a detection surface having a certain area.
  • Non-Patent Document 9 the sum of the multiple X-rays assuming a plurality of light rays, not a single light ray (not ignoring the multiplicity of X-rays) It is desirable to capture.
  • the spectral distribution of X-rays obtained by projection is estimated using spectroscopy or a filter, the X-rays of each spectral intensity are measured, and the X-ray absorption coefficient of each spectrum is estimated based on the measurement.
  • Non-Patent Documents 1 and 2 There is a technique for reducing the degree of freedom of a parameter to be estimated by using a basis function expression using (see Non-Patent Documents 1 and 2).
  • Non-Patent Documents 1 and 2 in order to deterministically estimate the organization class in which the type of organization is represented by a discrete variable, when the assumed number of organization classes is not appropriate or the estimated organization class is incorrect In this case, the solutions are greatly different, and there is a problem that appropriate expression cannot be performed for an organization that is not expressed in a known organization class. In addition, since there are no restrictions on the tissue class to be estimated, there is a problem that the defect setting property is not sufficiently solved and the estimation of the tissue class tends to be unstable. In particular, in Non-Patent Document 1, the ratio of the X-ray absorption coefficient is variable between wavelengths, and different tissue classes are weighted. If the number of tissue classes is small, an X-ray absorption coefficient that cannot be expressed may appear, and if the number of tissue classes is large, defect setting properties are generated without meaning of constraint.
  • the X-ray absorption coefficient is expressed “probabilistically” through the substance class, which is a discrete variable representing the type of substance, as compared with the image reconstruction using the prior art, and the substance class is probabilistic.
  • the substance class is probabilistic.
  • An X-ray CT image processing method, an X-ray CT image processing program, and an X-ray CT image on which the program is mounted which can acquire a reconstructed image equivalent to the conventional one with a lower X-ray exposure and can reduce artifacts
  • An object is to provide a processing apparatus.
  • the substance class is not necessarily a living tissue, and is used to include a tissue class.
  • the X-ray incident on the detector needs to be regarded as a sum of multiple X-rays assuming a plurality of light rays, not a single light ray.
  • a bias is generated by ignoring the multiplicity of X-rays and representing it with a single ray.
  • the present invention sets the probability distribution regarding the observation in which the projected X-rays are incident on the detector multiple times and performs statistical estimation to thereby reduce the bias caused by representing with one ray. It is an object to provide a processing method, an X-ray CT image processing program, and an X-ray CT image processing apparatus equipped with the program.
  • the X-ray CT image processing method of the present invention sets a probability distribution related to observation of an X-ray projection image and performs statistical estimation.
  • the probability distribution described above represents the observation process of the multiple X-ray sum in which the projected X-rays are incident on the detector multiple times.
  • the multiple X-ray sum is such that the detection surface of the X-ray in the detector has a predetermined spread instead of one point spatially, the detection time of the X-ray has a predetermined width instead of a moment, and
  • the X-ray spectral distribution is not a line spectrum but has a predetermined spread, and a plurality of X-rays having different spatial positions, directions, and wavelengths are incident on the detector, and a weighted sum of these multiple X-rays. Is observed by the detector.
  • the statistical estimation is based on the expected value of the X-ray absorption coefficient and the posterior distribution of the substance. Is estimated.
  • the observation process expressed by the multiple X-ray sum has the following features (a) and (b).
  • (a) The X-ray absorption coefficient having wavelength dependency of each pixel is expressed by the product of the X-ray absorption coefficient density of each pixel and the X-ray absorption coefficient ratio between X-ray wavelengths determined for each substance. .
  • the X-ray absorption coefficient density is considered by assuming that there is only one X-ray wavelength region in the X-ray absorption coefficient having wavelength dependence.
  • (b) Represents the observation process of multiple X-ray sums in which projected X-rays are incident on the detector multiple times; Multiple X-ray sum is The X-ray detection surface of the detector has a predetermined spread instead of one point spatially; The detection time of X-rays in the detector has a predetermined width instead of a moment, The spectral distribution of the X-ray has a predetermined spread rather than a line spectrum, Thus, it is assumed that a plurality of X-rays having different spatial positions, directions, and wavelengths are incident on the detector in a multiple manner. Then, it is assumed that a projection is calculated for each of a plurality of assumed light beams, and the sum of (electrical) signals excited by these projections is observed by the detector.
  • the prior distribution relating to the X-ray absorption coefficient indicating what X-ray absorption coefficient the imaging target tends to take is determined via a variable called a substance class.
  • the substance class is a variable defined by each pixel in the imaging target that indicates which substance the imaging target is composed of.
  • Material classes probabilistically constrain the degree of freedom of wavelength-dependent X-ray absorption coefficients by adding probabilistic constraints on which material classes are likely to occur and how spatially smooth each material class is. To do.
  • the probability distribution of what kind of X-ray absorption coefficient density is likely to be taken when the substance class is determined is expressed as the sum of two gamma distributions, express to which of each gamma distribution it belongs Subclass variables of the substance class to be introduced, how easily the combination of the substance class and the subclass with respect to the combination of the substance class and the subclass, and how much the combination of the respective substance class and the subclass is spatially smooth.
  • the degree of freedom of the wavelength-dependent X-ray absorption coefficient is stochastically restricted by adding a probabilistic constraint.
  • the wavelength dependence of the X-ray absorption coefficient of a substance A probability distribution considering the observation process of multiple X-ray sums such as sex is set and X-ray CT image processing is performed.
  • a probability distribution considering the observation process of multiple X-ray sums such as sex is set and X-ray CT image processing is performed.
  • Bayesian estimation allows the probability of substance class belonging to each pixel and the X-ray absorption coefficient density of each pixel to be represented stochastically, making it more detailed and more flexible than selecting a substance class deterministically.
  • X-ray absorption coefficient can be expressed.
  • the X-ray CT image processing method of the present invention by treating the variable of the material class as a random variable, a material that is not expressed by a known material class can be expressed flexibly by its weighted sum. Further, in the X-ray CT image processing method of the present invention, it is possible to suppress defect setting by imposing constraints such as spatial smoothness as prior knowledge on the substance class. Further, in the X-ray CT image processing method of the present invention, by performing approximate Bayesian inference, convergence is ensured and estimation instability (overfitting) can be suppressed. That is, by using the X-ray CT image processing method of the present invention, it is possible to perform stable estimation that is more faithful to the physical process and highly accurate.
  • the substance class is used to include a tissue class.
  • a substance class is expressed by a discrete variable.
  • the existence ratio (density) of a substance class refers to what kind of substance class (normal cells such as muscle, soft cells such as fat, bone, metal, etc.) This is a non-negative scalar parameter.
  • the parameter expressing the X-ray absorption coefficient ratio fixed between wavelengths for each substance class is the stochastic X-ray absorption coefficient, assuming that the ratio of absorption coefficient between wavelengths is fixed, This is used for the probability distribution of knowledge to restrict the degree of freedom of the wavelength-dependent X-ray absorption coefficient.
  • the parameter expressing the degree of spatial continuity of the substance class is, for example, the degree of spatial continuity such as the substance class of air is easy to connect pixels (no other tissues enter the air). Since the restriction of high can be imposed, the degree of freedom of the wavelength-dependent X-ray absorption coefficient is restricted by using the degree of spatial continuity of the substance class in the probability distribution of prior knowledge.
  • a case is assumed in which a sum related to an X-ray wavelength is included in the multiple X-ray sum in the X-ray CT image processing method of the present invention.
  • the prior distribution with respect to the X-ray absorption coefficient having the wavelength dependency of each pixel is characterized by a parameter expressing the distribution of the X-ray absorption coefficient for each substance class and for each X-ray wavelength. It probabilistically limits whether an X-ray absorption coefficient is likely to occur.
  • the X-ray absorption coefficient having the wavelength dependency of each pixel is not dependent on the substance class.
  • the X-ray absorption coefficient having the wavelength dependency of each pixel is characterized by a parameter expressing an X-ray absorption coefficient ratio between X-ray wavelengths that is determined for each substance class.
  • the prior distribution of the X-ray absorption coefficient having the wavelength dependency of each pixel is expressed by a conditional prior distribution with respect to the X-ray absorption coefficient density of each pixel on the condition of the material class and a prior distribution with respect to the material class. That is, when the X-ray absorption coefficient is expressed by the product of the X-ray absorption coefficient density of each pixel and a constant dependent on the substance class-dependent wavelength, the prior distribution of the X-ray absorption coefficient is the prior distribution with respect to the X-ray absorption coefficient density of each pixel. And a prior distribution for the substance class.
  • the distribution of the X-ray absorption coefficient density determined for each substance class is expressed by a mixed gamma distribution.
  • the substance class is expressed by a Boltzmann distribution characterized by a parameter that expresses the degree of appearance and a parameter that expresses the degree that the substance class is spatially continuous. This is an example of the expression of the prior distribution of the X-ray absorption coefficient.
  • a pre-distribution of the X-ray absorption coefficient is determined by assuming a mixed gamma distribution for the X-ray absorption coefficient density subject to the substance class and assuming a Boltzmann distribution for the substance class variable.
  • the reason why it is expressed as (mixed) gamma distribution as the distribution of X-ray absorption coefficient density defined for each substance class is that the distribution shape can express a unimodal distribution and the value of the random variable is It focuses on the fact that it can always be guaranteed that the value is positive.
  • the reason why the projected X-rays are expressed in a Gaussian distribution as the signal distribution observed in the detector in the process of observing the multiple X-ray sums that are incident on the detector in multiple ways is photon scattering, quantization noise, etc.
  • Bayesian estimation is performed by estimating the posterior distribution for each of the substance class, the X-ray absorption coefficient density, and the subclass variable of the substance class.
  • the subclass variable belongs to which gamma distribution when the probability distribution of what X-ray absorption coefficient density is likely to be obtained when the substance class is determined is expressed by the sum of two gamma distributions.
  • Is a discrete variable that specifies The prior distribution of the X-ray absorption coefficient under conditions given a substance class is likely to take a specific value (for example, an average value), but is not always a symmetric probability distribution with respect to the average value.
  • the kurtosis of the probability distribution can be expected to vary in size, so that it can be expressed as a mixed probability distribution that is a weighted sum of two probability distributions so that such a complex prior distribution can be expressed. is there. Furthermore, since the X-ray absorption coefficient takes only a positive value, this mixing probability distribution is assumed to be expressed by a mixing probability distribution of a gamma distribution.
  • the subclass variable of the substance class is a binary variable of 0 or 1 introduced for convenience in order to specify which probability distribution belongs to the mixed probability distribution of two distributions.
  • the X-ray CT image processing apparatus of the present invention is an X-ray CT image processing apparatus that sets a probability distribution related to observation of an X-ray projection image and performs statistical estimation, and includes measurement conditions including an X-ray projection image and at least an X-ray intensity.
  • Input means a probability distribution related to observation, a prior distribution related to an X-ray absorption coefficient expressed through a variable indicating which substance exists in each pixel of a reconstructed image called a substance class, and a substance class
  • An estimation means for estimating the X-ray absorption coefficient and the substance class by Bayesian estimation based on the expected value of the posterior distribution regarding the X-ray absorption coefficient and the substance class under the prior distribution is provided.
  • the probability distribution regarding the observation is set as a probability distribution representing what kind of projection image is likely to be obtained through the observation process of the multiple X-ray sum in which the projected X-rays are incident on the detector multiple times.
  • the multiple X-ray sum is such that the X-ray detection surface of the detector has a predetermined spread instead of one point spatially, the X-ray detection time has a predetermined width instead of a moment, and X
  • the spectral distribution of the line is not a line spectrum but has a predetermined spread, and a plurality of X-rays having different spatial positions, directions, and wavelengths are incident on the detector in a multiplexed manner, and the weighted sum of these multiplexed X-rays is It is to be observed with a detector.
  • a substance that deviates from a previously assumed substance class by treating a variable of the substance class as a random variable using a parameter expressing the existence ratio (density) of the substance class.
  • approximate expression is possible by expressing it with a weighted sum of previously assumed substance classes.
  • a parameter that expresses the degree to which the substance class is spatially continuous it is possible to suppress defect setting by imposing a restriction on spatial smoothness as a probability distribution for the substance class. ⁇
  • an X-ray CT image processing program for performing statistical estimation with a probability distribution relating to an X-ray absorption coefficient, and executing the following steps 1) to 2) on a computer: It is a program to make it. 1) Input the measurement conditions including the X-ray projection image and at least the X-ray intensity 2) Based on the probability distribution relating to the observation of the X-ray projection image and the prior distribution relating to the X-ray absorption coefficient having the parameter relating to the substance. The step of estimating the X-ray absorption coefficient and the substance with the expected value of the posterior distribution regarding the X-ray absorption coefficient and the substance.
  • the multiple X-ray sum is that the X-ray detection surface of the detector has a predetermined spread instead of one point in space, and the X-ray detection time has a predetermined width instead of a moment.
  • the X-ray spectral distribution is not a line spectrum but has a predetermined spread, and a plurality of X-rays having different spatial positions, directions, and wavelengths are incident on the detector in a multiplexed manner. Weighted sum is detector It is to be observed at.
  • the present invention can provide an X-ray CT image processing apparatus equipped with the above X-ray CT image processing program.
  • the X-ray absorption coefficient is stochastically expressed as compared with the image reconstruction using the prior art, thereby enabling more flexible expression and at the same level as the conventional case with a lower X-ray exposure.
  • the reconstructed image can be obtained, and the beam hardening artifact can be reduced.
  • a bias is generated by representing a single light ray by setting a probability distribution related to the observation process of the multiple X-ray sum in which the projected X-rays are incident on the detector multiple times and performing statistical estimation. It is possible to reduce the effect.
  • a prior distribution relating to the X-ray absorption coefficient is set, and image reconstruction is performed using statistical estimation.
  • the beam hardening artifact is eliminated by estimating the energy-dependent (wavelength-dependent) X-ray absorption coefficient based on a statistical model.
  • a basis function expression is used to set the basis as an energy-dependent (wavelength-dependent) X-ray absorption coefficient of a known substance, and a discrete variable that expresses which substance class the weight belongs to, and an X-ray absorption Given the product of non-negative scalar continuous variables expressing the coefficient density, we impose a constraint on the solution by imposing a prior distribution of the X-ray absorption coefficient density and a prior distribution of the material class subject to the material class.
  • the human tissue in the case of a reconstructed image for human tissue, the human tissue is composed of limited substances such as fat, muscle, and bone, and the rough distribution of the X-ray absorption coefficient of each substance is known in advance. Therefore, the knowledge is expressed as prior distribution and used for estimation. For this estimation, it is necessary to estimate which substance class each pixel belongs to, but each substance has spatial continuity, and the proportion of each substance in the human body is roughly Therefore, by introducing prior knowledge such as spatial smoothness constraints and occupancy ratios to discrete variables representing substance classes, the failure setting is suppressed.
  • the observation process related to the observation of projected X-rays expresses the observation process of multiple X-ray sums in which the projected X-rays are incident on the detector multiple times. ) Represents that a plurality of X-rays having different wavelengths are incident on the detector.
  • the probability distribution related to the reconstructed image characterizes what value the X-ray absorption coefficient is likely to take in each pixel area of the reconstructed image, and each substance class and each subclass variable
  • a mixed gamma distribution that expresses what X-ray absorption coefficient density is likely to be taken for each X-ray wavelength, and a parameter that expresses how easily each substance class and each subclass variable combination appears This is expressed by a Boltzmann distribution consisting of parameters representing how much the combination of each substance class and each subclass variable is likely to be spatially smooth and continuous.
  • T is the number of projections
  • I is the number of detectors
  • y i (t) is the number of photons detected by the i-th detector. Since X-rays are exponentially attenuated when passing through a substance, they can be expressed as the following formula (1).
  • b i (t) represents the number of photons emitted from the X-ray source (the number of photons that can be observed when no object is placed).
  • L ij (t) corresponds 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 projections (see Figure 2).
  • x je is the X-ray absorption coefficient of the j-th pixel with respect to the X-ray of energy e.
  • Each element z jc is a binary variable that takes 0 or 1, and takes 1 when the j-th pixel belongs to the substance class c, and takes 0 otherwise.
  • the X-ray absorption coefficient x je can be expressed as the following mathematical formula (3).
  • z j ) can be expressed as the following formula (7) as a mixed distribution of two distributions so that asymmetry of distribution and various different kurtosis can be expressed.
  • the variable b j is a binary variable.
  • the X-ray absorption coefficient x j can be expressed as the following equation (8) assuming that it follows a gamma distribution.
  • u cd and v cd are parameters of a gamma distribution, and are average u cd / v cd and variance u cd / v 2 cd . These parameters are determined by empirical knowledge about the substance class.
  • variable z jc is a binary variable that takes 0 or 1 and satisfies the following mathematical formula (9).
  • the X-ray absorption coefficient X, the substance class Z, and the variable B are estimated as the posterior distribution of the observation data D.
  • D) of the X-ray absorption coefficient X, the substance class Z, and the variable B can be expressed by the following formula (12) by Bayes' theorem.
  • Bayesian estimation In Bayesian estimation, the expected value of the posterior distribution p (X, Z, B
  • the test distribution q (X, Z, B) can be arbitrarily selected as long as the above equation (12) can be minimized or approximated. it can.
  • D) is evaluated by a KL (Kullback-Leibler) distance, and the test distribution q (X , Z, B).
  • the KL distance can be expressed as the following formula (13).
  • ⁇ •> q (X, Z, B) is an operator representing performing an integral calculation with respect to the distribution q (X, Z, B).
  • ⁇ •> q (X, Z, B) represents taking an expected value with respect to q (X, Z, B).
  • test distribution q (X, Z, B) is given by the product of its peripheral distributions q (X) and q (Z, B), i.e. the variables X and ⁇ Z, B ⁇ Assume that there is independence between them. Further, it is assumed that each element of q (X) and q (Z, B) is independent between pixels, and further, each element q (x j ) of q (X) follows a gamma distribution. At this time, q (z j , b j ) inevitably has a multinomial distribution because the random variables z j and b j are discrete variables.
  • test distribution q (X, Z, B) can be expressed as the following mathematical formulas (14) and (15). However, since it is difficult to optimize the test distributions q (X) and q (Z, B) at the same time, the optimization is performed alternately as in the following formulas (16) and (17).
  • the SCG (Scaled Conjugate Gradient) method was used for optimization of the parameters ⁇ j and ⁇ j of q (x j ).
  • Test distribution q (X) and j th other than the pixel q (z j ', b j ') ( where, j'is different from j) is under a fixed, minimizing the KL distance q (z j , B j ) can be obtained analytically.
  • This analytical solution depends on q (z j ′, b j ′) (j′ ⁇ (j)) of neighboring pixels of j.
  • a subset of pixels having a checkered pattern is alternately optimized as shown in FIG. That is, the parameter q (Z, B) does not update all the pixels at the same time, but alternately optimizes a subset of the pixels with the checkered pattern shown in FIG. Figure 3 (a) q (z j , b j) corresponding to the pixel of the gray area After updating updates the q (z j, b j) corresponding to the pixel of the gray portion shown in FIG. 3 (b) .
  • the processing contents in the X-ray CT image processing method and the X-ray CT image processing program of the present invention have been described above.
  • a computer experiment is performed to examine the effectiveness of the present invention, and the screen reconstruction of the X-ray CT image processing method and the X-ray CT image processing program of the present invention is compared with the screen reconstruction according to the prior art. And evaluate.
  • the usefulness of the X-ray CT image processing method and X-ray CT image processing program of the present invention can be understood from the following evaluation experiment.
  • an appropriate true image is prepared so that an error can be measured between the reconstructed image and the true image, and an X-ray projection is simulated on the true image to artificially generate a projection image (sinogram).
  • X-ray projection in order to consider the energy dependency (wavelength dependency) of the X-ray absorption coefficient of a substance, the energy of incident X-rays is five types of 60 keV, 70 keV, 80 keV, 90 keV, and 100 keV (each photon number is , 6 ⁇ 10 4 , 2.5 ⁇ 10 4 , 1.5 ⁇ 10 4 , 7 ⁇ 10 3 , 5 ⁇ 10 3 ), and each was attenuated according to the X-ray absorption coefficient specific to the substance.
  • a noise noise generated according to Poisson noise is added to the attenuated X-ray photon.
  • the number of detectors is 95 (pieces)
  • the resolution of the image is 64 ⁇ 64 (pixels)
  • the projection is 1 ° to 180 ° at 1 ° intervals.
  • PSNR Peak Signal to Noise Ratio
  • the FBP method a method of performing a Bayesian estimation using a prior distribution that does not consider the X-ray energy distribution as a Gaussian distribution (hereinafter abbreviated as “conventional method”), and the X-ray of the present invention in which the X-ray energy is discretized into three.
  • Image reconstruction is performed by three types of CT image processing methods (hereinafter abbreviated as “method of the present invention”), and the comparison results are shown.
  • the prior distribution of metal and bone is as shown in FIG. FIG. 4 shows the prior distribution of metal and bone and the X-ray absorption coefficient at 60 keV of the true image.
  • the X-ray absorption coefficients of metal and bone are probability distributions expressed by gamma distributions with an average of 0.7 and 0.4, respectively.
  • the X-ray absorption coefficient of the true image is 5% smaller than the average value of the prior distribution at each X-ray energy.
  • Example 1 a computer experiment for image reconstruction for a phantom including four metals will be described with reference to FIG. 5A is a true image, FIG. 5B is a reconstructed image by the FBP method, FIG. 5C is a reconstructed image by the naive Bayes method, and FIG. 5D is a reconstructed image by the method of the present invention.
  • an image reconstruction experiment was performed on the true image shown in FIG. 5A by using the FBP method, the naive Bayes method, and the method of the present invention.
  • the total number of incident X-ray photons is 1.12 ⁇ 10 5
  • the X-ray energy is discretized into three of 60 keV, 80 keV, and 100 keV, and the number of photons belonging to each of them. Are 7.5 ⁇ 10 4 , 3 ⁇ 10 4 , and 10.7 ⁇ 10 4 .
  • the number of repetitions was set to 25 for comparison.
  • the estimation results are shown in FIGS. 5 (b) to (d).
  • PSNR in each of the FBP method, the naive Bayes method, and the method of the present invention is shown in Table 1 below.
  • Table 1 the PSNR, it can be understood that the method of the present invention is larger than the FBP method and the naive Bayes method, and the method of the present invention can perform estimation with higher accuracy.
  • Example 2 a computer experiment for image reconstruction for a tooth phantom including a metal implant will be described with reference to FIG. 6A is a true image, FIG. 6B is a reconstructed image by the FBP method, FIG. 6C is a reconstructed image by the naive Bayes method, and FIG. 6D is a reconstructed image by the method of the present invention.
  • an image reconstruction experiment was performed on the true image shown in FIG. 6A using the FBP method, the naive Bayes method, and the method of the present invention.
  • the total number of incident X-ray photons is 1.12 ⁇ 10 5
  • the X-ray energy is discretized into three of 60 keV, 80 keV, and 100 keV, and the number of photons belonging to each of them. Are 7.5 ⁇ 10 4 , 3 ⁇ 10 4 , and 10.7 ⁇ 10 4 .
  • the number of repetitions was set to 25 for comparison.
  • the estimation results are shown in FIGS. 6 (b) to (d).
  • the PSNR in each of the FBP method, the naive Bayes method, and the method of the present invention is shown in Table 2 below. Similar to the experimental results of Example 1, in PSNR, the method of the present invention is larger than the FBP method and the naive Bayes method, and it can be understood that the method of the present invention can perform estimation with higher accuracy.
  • the multiple X-ray sum is such that the X-ray detection surface of the detector has a spatially finite extent, and the X-ray detection time of the detector has a finite width. It is possible to express that a plurality of X-rays with different spatial positions, directions, and wavelengths are incident on the detector by having a wavelength region (spectral distribution) with a finite spread. To do.
  • a parallel beam will be described as an example with reference to FIG.
  • a projection by a parallel beam on a two-dimensional plane is assumed.
  • a direction in which X-rays are emitted is taken as a v-axis
  • a direction orthogonal to the direction is taken as a u-axis.
  • variables are defined as follows.
  • x (u, v, ⁇ ) be the X-ray absorption coefficient of the wavelength ⁇ of the substance at position (u, v).
  • y i ( ⁇ ) be the number of photons of wavelength ⁇ incident on the detection surface S i of the i-th detector per unit time.
  • d ( ⁇ ) be the number of photons of wavelength ⁇ emitted per unit area by a parallel beam X-ray source
  • ⁇ ( ⁇ ) be the average number of electrons excited by one photon of wavelength ⁇ .
  • the expected value of the current Y i observed by the i-th detector per unit time is expressed as the following formula (19). Furthermore, when a current for a certain period of time is stored and observed, integration over time is performed.
  • the line integral is obtained when there is a spatial spread of the detection surface or a spread of the X-ray spectral distribution. It can be seen that artifacts occur when using an estimation method that assumes. Similarly, it is not possible to reduce to line integration when considering rotation in the time that the detector is receiving light. Since calculation is difficult in the above-mentioned multiple integral form, it is considered to approximate the integral with a sum. This indicates that multiple X-ray sums will be considered. First, paying attention to only a specific wavelength ⁇ in the above equation (19), the integral relating to u and v in the following equation (20) is replaced with a sum.
  • Discretization is performed in order to approximately express the integral related to u and v on the computer (see FIG. 8).
  • the integral relating to v in the above equation (20) can be approximated as the following equation (21).
  • x j is the X-ray absorption coefficient of the j-th pixel
  • I (u) is a set of indices of the pixels at the u position on the u-axis
  • ⁇ L represents the length of one side of the pixel.
  • 1 / N indicates the width of the strip when the integral is cut into strips and approximated. This corresponds to considering the sum of projections expressed by line integration for each index k. From the above, it can be seen that multiple X-rays need to be considered when considering the spatial extent of the detector.
  • ⁇ N 1 / N
  • r (k) represents the position of the ray k on the detection plane
  • r (k) ⁇ S i represents a set of rays that reach the detector i.
  • I k represents a set of indices of pixels at which the ray k intersects
  • l kj represents a length at which the ray k intersects the pixel j.
  • the formula (25) can be considered as a case where l kj of the formula (26) is 1.
  • the above formula (26) is established by calculating l kj in consideration of the direction of light rays in addition to the parallel beam, fan beam, and cone beam. What can be understood from the above equation (26) is that when a projection observation model that considers the spread of the detection surface is considered, it is necessary to consider the projection for each ray passing through the detection surface.
  • the matrix L is (the number of detectors) ⁇ (the number of pixels of the reconstructed image), but when the detection surface expansion is considered, the matrix L is (The number of rays passing through the detector) ⁇ (the number of pixels of the reconstructed image). If the number of rays is large, the integral can be accurately approximated, but the calculation amount increases.
  • the photons y i, t incident on the detector at each time t are expressed as the sum of a plurality of light beams as in the above equation (28). Integration over time is equivalent to considering a sum over a new ray bundle to include additional sums. Therefore, when the set S i of ray indexes is set to S i, t so as to depend on time, and the set I k of pixel indexes through which each ray k passes is redefined, the measurement is performed at the detection time (exposure time) T.
  • the number of photons can be expressed by the following formula (30).
  • the present invention is useful not only for medical use but also for industrial use as an X-ray CT apparatus for reducing beam hardening artifacts. It is also useful as a belt conveyor type X-ray CT apparatus or dental X-ray CT apparatus in which movement and rotation by the belt conveyor cannot be ignored.

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Optics & Photonics (AREA)
  • Public Health (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Veterinary Medicine (AREA)
  • High Energy & Nuclear Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Surgery (AREA)
  • Biophysics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Pulmonology (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

　Ｘ線吸収係数を確率的に表現することで、より柔軟な表現を可能とし、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、ビームハードニングアーティファクトを低減可能なＸ線ＣＴ画像処理方法を提供する。投影Ｘ線の観測に関する確率分布を設定し統計推定を行う。確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現される。多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有すること、検出器におけるＸ線の検出時間が一瞬ではなく所定の幅を有すること、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有することによって、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射するとして表現される。そして、多重Ｘ線和を表現する観測過程と物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、統計推定が事後分布の期待値によるベイズ推定を行う。

Description

Ｘ線ＣＴ画像処理方法，Ｘ線ＣＴ画像処理プログラム及びＸ線ＣＴ画像装置
　本発明は、体内組織構成を考慮したＸ線ＣＴ（Computed Tomography）画像処理方法、Ｘ線ＣＴ画像処理プログラム及びＸ線ＣＴ画像装置に関する。
　Ｘ線ＣＴ（Computed Tomography）は、観測対象にＸ線を照射することで得られる投影データから観測対象内部をＸ線吸収係数の分布として画像化する技術である。
　一般に、Ｘ線ＣＴは、Ｘ線管にＸ線検出器を対向配置するとともに、これらの間にターンテーブルを配置した構成を取る。このターンテーブルの回転軸は、Ｘ線管とＸ線検出器を結ぶＸ線光軸に対して直交する向きとされる。Ｘ線ＣＴ装置においては、測定対象物に対して複数の角度からの透過Ｘ線データを収集してＣＴ断層画像の再構成を行う。一般的には、１８０°以上、回転させて透過Ｘ線データを収集する。
　従来から、より低いＸ線被曝量で従来と同程度のＣＴ画像を再構成する技術、メタルアーティファクトの低減技術など、Ｘ線ＣＴ再構成の高精度化を可能にするアルゴリズムが研究されている。上記のメタルアーティファクトは、測定対象物にＸ線を強く吸収する金属などの密度の高い高吸収物質があったときに、正確な再構成ができず再構成画像に人為的なノイズがのるというものである。Ｘ線ＣＴ再構成は、様々な投影角度において得られる投影像から、測定対象物を透過したＸ線の減衰の度合いを求め、内部の物質のＸ線吸収係数の分布を推定する手法である。しかし、高吸収物質はＸ線を強く減衰させるため、ほとんど透過Ｘ線が検出されず観測ノイズに紛れてしまう。特に既存手法であるフィルター補正逆投影法（ＦＢＰ：Filtered Back Projection）を用いた場合、ＦＢＰ法ではノイズの影響を受けやすいことから、再構成画像にはメタルアーティファクトと呼ばれる虚像が生じ、正確な画像の再構成ができないという問題がある。このようにＦＢＰ法はノイズの影響を受けやすいのであるが、主なノイズにはポアソン分布に従うショットノイズが含まれるため照射Ｘ線の強度を高くすることで信号雑音比を向上することができる。これによって相対的なノイズの影響を小さくし、画像再構成の精度を向上させることができるが、その一方で放射線被曝量が増大させてしまい、発癌などの確率的健康リスクを生じさせてしまう。
　このように、Ｘ線被曝量の低減とＳＮ比はトレードオフの関係があり、より低いＸ線被曝量で従来と同程度のＳＮ比の再構成画像を得る技術として、統計推定法が提案されている。統計推定法では、観測過程が確率的なノイズによる不確実性を含むことや、撮像対象の物体に関する事前知識を利用することで、低被曝環境下における推定の不良設定性を軽減することができる。 結果として、従来と同じＸ線強度や解像度でも再構成画像に含まれるアーティファクトを低減できる。
　現在、画像再構成法の主流であるフィルター補正逆投影法（ＦＢＰ：Filtered Back Projection）や、ＦＢＰを改良したＰＣＬＩＳ（Projection Completion Method based on a Linear Interpolation in the Sinogram）は、統計推定を行わないものである。
  統計推定法は、最尤推定法（Maximum Likelihood Estimation: ＭＬＥ），ＭＡＰ（Maximum A Posteriori）推定法，ベイズ推定法の３手法に大別される。
　最尤推定法（ＭＬＥ)は、推定対象である再構成画像に関する事前知識をおかない統計推定法で、観測における不確実性を含む物理過程を確率モデルで表現して、その確率モデルを基に最も尤もらしい再構成画像を推定するものである。観測における不確実性を含む物理過程には、観察される透過Ｘ線のフォトン数に関する揺らぎ(ショットノイズ)などが含まれる。
　ＭＡＰ推定法，ベイズ推定法は、推定対象に関する事前知識を確率モデルの形式で表現し、推定に取り入れるものである。すなわち、ＭＡＰ推定法やベイズ推定法は、観測に関する確率モデルと推定対象に関する確率モデルの両方を勘案して、推定が行われる。 Ｘ線CTにおいては、推定対象はＸ線吸収係数である。具体的には、ＭＡＰ推定法では、観測に関する確率モデルと推定対象となる物体のＸ線吸収係数に関する確率モデルから計算されるＸ線吸収係数の事後分布の最大値を推定値とし、ベイズ推定法では同じＸ線吸収係数の事後分布の期待値を推定値とする。ＭＡＰ推定法では、事後分布が最大となる一点のみを取り上げて推定を行うが、ベイズ推定法では事後分布の期待値をとるため、事後分布が大きくなるＸ線吸収係数の集団に対する平均によって推定を行う。 事後分布が最大となる点は、観測ノイズによって敏感に変化しうるが、事後分布の期待値は観測ノイズに対してそれほど敏感に影響されないため、ベイズ推定法のほうがより安定し、精度の高い推定となる。
　ここで、通常のＸ線ＣＴアルゴリズムの多くは画像再構成にあたって、照射するＸ線が線スペクトルであると仮定し単一の吸収係数を推定している。しかしながら、実際のＸ線は単一エネルギーではなく、Ｘ線管に加えられた最高管電圧で決まるエネルギーを最大として低い方へ連続的に分布するエネルギースペクトルで構成されている。そのため、物質にＸ線を照射するとエネルギーの高いＸ線はほとんど吸収されずに透過するが、エネルギーの低いＸ線は減弱が大きいため、物質を透過した後にはＸ線のスペクトル分布は高エネルギー側にシフトするビームハードニングと呼ばれる現象が生じる。
　物質のＸ線吸収係数のエネルギー依存性（波長依存性）の影響によって生じる現象は、Ｘ線源が線スペクトルであると仮定し、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しない観測過程のモデルでは表現できず、正しくない推定が行われてしまう。このように物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルによる推定でビームハードニングアーティファクトが生じる。
　特に、金属などの極端にＸ線吸収係数が高い高吸収物体があるとストリークアーティファクトと呼ばれる筋状のアーティファクトや高吸収物体周辺のＸ線吸収係数が実際より小さく見積もられるダークバンドアーティファクトが生じることが知られているが、これはビームハードニングの影響を考慮していないことが大きく影響していると考える。
　これについて、以下に説明する。
　金属など非常に高いＸ線吸収係数をもつ高吸収物体を透過すると、Ｘ線は大きく減弱する。その性質が観測モデルに表現されていれば、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルであっても、高吸収物体近辺に大きなＸ線吸収係数をもつ物体が存在するであろうことを推定できる。しかしながら、高吸収物体でも高エネルギーのＸ線は透過しやすい。そのため、実際にはＸ線源の主要なエネルギーをもつフォトンの多くが高吸収物体を透過中に吸収されていたとしても多少のＸ線が観測される。高吸収物体が離れた位置に複数ある場合、その複数の高吸収物体のいずれも透過するＸ線は、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルでは高いＸ線吸収係数をもつ物体が重なると透過Ｘ線がさらに大きく減衰することを予想する。 しかし、実際にはある一定量の高エネルギーのＸ線は高吸収物体を複数透過したとしても観測されてしまう。このため、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルで、単純に高吸収物体が存在すると仮定すると実際の観測と矛盾が生じてしまう。このようにして、物質のＸ線吸収係数のエネルギー依存性を考慮しない観測モデルと実際の観測とができるだけ矛盾しないようにＸ線吸収係数の推定を行うと、高吸収物体周辺のＸ線吸収係数が実際より小さく見積もられてしまうダークバンドアーティファクトが生じる。
　メタルアーティファクトを軽減する方法として、投影データを金属領域、非金属領域に分離し、それぞれの領域において画像再構成し、得た画像を重ね合わせるという方法が既に知られている。このようにすることで、金属領域における投影データの矛盾を正し、メタルアーティファクトの発生を軽減することができる。また、本願発明者らは、再構成画像となるＸ線吸収係数に関する事前知識をおいて統計推定を行うＸ線ＣＴ画像処理方法であって、事前知識を、再構成画像の各画素の領域において定義されるパラメータで、撮像対象となる人体の各組織の存在割合を表現するパラメータと、各組織のＸ線吸収係数を表現するパラメータと、各組織の空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現するものを用いることで、メタルアーティファクトを低減可能なＸ線ＣＴ画像処理方法を提案している（特許文献１を参照）。
　しかし、これらの方法は、Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮するわけではないため、実際の観測過程と物質のＸ線吸収係数のエネルギー依存性を考慮しない観測過程とのズレのために生じるビームハードニングアーティファクトを除去するには至っていない。
　そして、ビームハードニングアーティファクトを軽減する手法として、Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮したモデルが知られている（非特許文献１～４）。Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮するモデルでは、Ｘ線のエネルギー毎のＸ線吸収係数を推定する必要がある。これは、より正確な物理過程に基づくが、推定問題としては不良設定性の問題を生じさせる。たとえば、単純のためＸ線のエネルギーを離散化し、二種類のエネルギーにおけるＸ線吸収係数を推定する問題を考えたとする。このとき、観測される投影像は従来と同じであっても、推定しなければならないＸ線吸収係数は従来に比べて２倍に増えてしまう。このように解となるＸ線吸収係数の自由度が高まるため、実際の観測とそれほど矛盾を生じない解がいくつも生じてしまう。この場合、ちょっとした観測ノイズに推定結果が敏感に影響され安定しない。 そこで、解となるＸ線吸収係数の自由度を抑える制約が必要となる。 先ず、非特許文献１～３では、各画素におけるエネルギー依存の吸収係数を既知の物質の吸収係数関数の重み付き和で表現することで、物質のエネルギー依存性（波長依存性）を表現している。Ｘ線吸収係数のエネルギー依存性（波長依存性）は、物質毎に異なるが高エネルギーになるに従って減衰する傾向など各物質間で共通する性質が存在する。Ｘ線のエネルギーを入力変数とし、そのエネルギーでのＸ線吸収係数を出力する関数を吸収係数関数と呼ぶことにすると、異なる物質であっても類似した吸収係数関数をもつため、どの物質の吸収係数関数も少数の、典型的には二つの物質の吸収係数関数の重み付き和（基底関数表現）でよく表現できる。例えば、非特許文献２では、物質クラスの決定をＦＢＰで得られた画像をもとに行い、その知識をもとに画像再構成を行っている。非特許文献４では、画素毎のＸ線吸収係数は、既知の物質種毎のＸ線吸収係数の組み合わせで表現できるという仮定で推定を行っている。
　このように、Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮する方法は、その不良設定性を抑制するため、既知の物質のＸ線吸収係数の波長依存性を利用した基底関数表現を用いることで、推定するパラメータの自由度を減らすことが行われてきた。しかしながら、その一方で、強すぎる制約が課せられる問題が指摘されている。
　一方、ハードウェア的には、ＤＥＣＴ（Dual Energy Computed Tomography）の登場により、異なるＸ線管電圧で撮影し、Ｘ線吸収係数の違いから物質のエネルギー依存性（波長依存性）を推定できるようになった。このＤＥＣＴを用いた様々な画像再構成法が知られている（非特許文献５～８）。しかし、ＤＥＣＴを用いた画像再構成法は、一般的に単色Ｘ線ＣＴより被曝量が大きく、同程度の被曝量に抑えようとすると推定精度が悪くなってしまうなどの問題がある。また、非特許文献５は、事後的に画像処理によって補正する方法であり、その補正によって更なる虚像を生じさせる可能性もあり限界がある。
　物質のＸ線吸収係数のエネルギー依存性を考慮しない観測モデルは、検出器に入射するＸ線が複数の波長からなるのではなく、ある代表的な波長のみからなることを想定することになる。 これは、検出器に複数のＸ線スペクトルが多重に入射する過程を、単純化したものと考えることができる。このようにＸ線が多重に検出器に入射する現象は、スペクトル空間だけではなく、時間的、空間的にも生じる。すなわち、検出器はある一定の検出時間内に多重のＸ線が入射し、ある一定の面積の検出面内に多重のＸ線が入射する。これら投影したＸ線が多重に検出器に入射する観測過程において、検出器に入射するＸ線を一本の光線で代表させた場合、Ｘ線吸収係数のエネルギー依存性を考慮しない時にビームハードニングアーティファクトが生じるのと同様にＸ線吸収係数の推定にバイアスが生じることが知られている（非特許文献９）。つまり、投影したＸ線が多重に検出器に入射する観測過程においては、一本の光線ではなく（Ｘ線の多重性を無視しないで）、複数の光線を想定した多重のＸ線の和と捉えることが望ましい。
特開２０１１－１５６３０２号公報
Joseph A. O ’Sullivan, and Jasenka Benac, Alternating Minimization Algorithms for Transmission Tomography, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 26, NO. 3 (2007), 283-297 Idris A. Elbakri and Jeffrey A. Fessler, Senior Member, IEEE, Statistical Image Reconstruction for Polyenergetic X-Ray Computed Tomography, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 21, NO. 2 (2002), 89-99 Joseph A. O’Sullivan, Bruce R. Whiting, Donald L. Snyder, and Orville A. Earl, Image Reconstruction from Data Acquired With an X-Ray Computerized Tomographic System Having Energy-Integrating Detectors Chye Hwang Yan, Robert T. Whalen, Gary S. Beaupre, Shin Y. Yen, and Sandy Napel, Reconstruction Algorithm for Polychromatic CT Imaging: Application to Beam Hardening Correction, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 19, NO. 1, (2000) A J Coleman and M Sinclair, A beam-hardening correction using dual-energy computed tomography, Phys. Med. Biol., 1985, Vol 30, No. 11, (1985), 1251-1256 Jeffrey A. Fessler, Idris Elbakri, Predrag Sukovic, Neal H. Clinthorne, Maximum-likelihood dualenergy tomographic image reconstruction, Medical Imaging 2002: Image Processing, Milan Sonka, J. Michael Fitzpatrick,38 Editors, oceedings of SPIE Vol. 4684 (2002), 38-49 Joonki Noh, Jeffrey A. Fessler Paul E. Kinahan, Statistical Sinogram Restoration in Dual-Energy CT for PET Attenuation Correction, IEEE Trans Med Imaging, (2009) 1688-1702 山崎陽一, 戸田尚宏, 低被爆環境下におけるDual-Energy Ｘ線CT雑音低減アルゴリズム, IEICE Technical Report MBE2008-113, (2009), 101-106 Ruhrnschopf, E. P., G. Schwierz, et al. (1981). Nonlinearity and Inhomogenity Effects in Computerized Tomography Due to the Exponential Attenuation of Radiation. Mathematical Aspects of Computerized Tomography. G. T. Herman and F. Natterer, Springer Berlin Heidelberg. 8：252-269
　上述の如く、波長が短くエネルギーの強い硬Ｘ線は、波長が長くエネルギーの弱い軟Ｘ線が物体を透過する力に比べて、物体を透過しやすいため、物体に入射したＸ線のスペクトル分布は物体を通過後、軟Ｘ線が減衰し硬Ｘ線側にパワースペクトルのピークが偏ったものとなり、ビームハードニングと呼ばれる現象を生じている。
　理想的には、投影によって得られるＸ線のスペクトル分布を分光法やフィルターを用いて推定し、各スペクトル強度のＸ線を測定し、その測定をもとに各スペクトルのＸ線吸収係数の推定を行えば問題は解決するのであるが、実際に分光すると各スペクトルで得られる信号のＳＮ比が下がるため推定精度は悪化するという問題がある。 そのため従来から、情報処理アルゴリズムでビームハードニングアーティファクトの除去を行うことが要望されている。
　様々なビームハードニングアーティファクトの除去を行う情報処理アルゴリズムが提案されている中で、物質のＸ線吸収係数の波長依存性を考慮する方法があり、特に、組織のＸ線吸収係数の波長依存性を利用した基底関数表現を用いることで、推定するパラメータの自由度を減らす技術がある（非特許文献１，２を参照）。
　しかし、非特許文献１，２では、組織の種類を離散変数で表現した組織クラスを決定論的に推定するために、想定する組織クラス数が適切でない場合や、推定された組織クラスが誤っていた場合、解が大きく異なってしまい、既知の組織クラスで表現されない組織に対して適切な表現ができないといった問題がある。
　また、推定する組織クラスに関して制約がおかれなかったため、不良設定性の解消が十分でなく、組織クラスの推定が不安定になりがちであるといった問題がある。
　特に、非特許文献１では、波長間でＸ線吸収係数の比が可変であり、異なる組織クラスを重み付けする。組織クラス数が少ないと表現できないＸ線吸収係数が現れる可能性があり、組織クラス数が多いと拘束の意味をなさずに不良設定性が生じてしまう。
　そこで、本発明は、従来技術を用いた画像再構成と比べて、物質の種類を表す離散変数である物質クラスを介してＸ線吸収係数を「確率的」に表現し、物質クラスに確率的な拘束を設けることで不良設定性を抑制する一方で、各画素がどの物質クラスに所属するかを決定論的ではなく確率的に推定することでより柔軟なＸ線吸収係数の表現を可能とし、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、アーティファクトを低減可能なＸ線ＣＴ画像処理方法、Ｘ線ＣＴ画像処理プログラムおよび該プログラムが搭載されたＸ線ＣＴ画像処理装置を提供することを目的とする。　ここで、本発明において物質クラスは、必ずしも生体の組織とは限らないので組織クラスを包含する意味で用いている。
　また、検出器に入射するＸ線は、一本の光線ではなく、複数の光線を想定した多重のＸ線の和と捉えることが必要である。上述の如く、このＸ線の多重性を無視して、一本の光線で代表させることによりバイアスが生じる。
　そこで、本発明は、投影したＸ線が多重に検出器に入射する観測に関する確率分布を設定し統計推定を行うことで、一本の光線で代表させることによるバイアスを低減可能なＸ線ＣＴ画像処理方法、Ｘ線ＣＴ画像処理プログラムおよび該プログラムが搭載されたＸ線ＣＴ画像処理装置を提供することを目的とする。
　上記状況に鑑みて、本発明のＸ線ＣＴ画像処理方法は、Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行う。
　上記の確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を表現する。ここで、多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである。
　そして、多重Ｘ線和を表現する観測過程と物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、統計推定がＸ線吸収係数及び物質に関する事後分布の期待値によって線吸収係数及び物質を推定する。
　ここで、多重Ｘ線和で表現される観測過程は、以下の（ａ），（ｂ）の特徴を備える。
(ａ) 各画素の波長依存性を有するＸ線吸収係数は、各画素のＸ線吸収係数密度と、物質毎に定められているＸ線波長間のＸ線吸収係数比の積で表現される。Ｘ線吸収係数の波長依存性を考えない場合のＸ線吸収係数の推定は、波長依存性を有するＸ線吸収係数において想定するＸ線波長域が１つしかないと考えてＸ線吸収係数密度のみを推定することと等価となる。そのため、波長依存性を有するＸ線吸収係数の各画素のＸ線吸収係数密度と、物質毎に定められているＸ線波長間のＸ線吸収係数比の積による表現は、Ｘ線吸収係数の波長依存性を考えない場合のＸ線吸収係数の表現の拡張となっている。
(ｂ)投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を表現するものであり、
　多重Ｘ線和は、
　検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有すること、
　検出器におけるＸ線の検出時間が一瞬ではなく所定の幅を有すること、
　Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有すること、
によって、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射することを想定する。そして、想定した複数の光線のそれぞれで投影を計算し、それらの投影によって励起された（電気）信号の和が検出器で観測されるとする。
　また、撮像対象がどのようなＸ線吸収係数をとる傾向があるかを表すＸ線吸収係数に関する事前分布は、物質クラスという変数を介して決められる。 物質クラスは、撮像対象の内部がどの物質から構成されるかを表す撮像対象内部の各画素で定義される変数である。物質クラスは、どの物質クラスがどの程度生じやすいか、またそれぞれの物質クラスの空間的滑らかさの度合いについて確率的な制約を加えることにより波長依存のＸ線吸収係数の自由度を確率的に制約する。 
　また、物質クラスが定まった時に、どのようなＸ線吸収係数密度をとりやすいかという確率分布を２つのガンマ分布の和で表現される場合には、各ガンマ分布のいずれに所属するかを表現する物質クラスのサブクラス変数を導入し、物質クラスとサブクラスの組み合わせに対してどの物質クラスとサブクラスの組み合わせがどの程度生じやすいか、またそれぞれの物質クラスとサブクラスの組み合わせがどの程度、空間的に滑らかさであるか、の確率的な制約を加えることにより波長依存のＸ線吸収係数の自由度を確率的に制約する。
　本発明のＸ線ＣＴ画像処理方法によれば、ビームハードニングアーティファクトをはじめとする多重Ｘ線和の観測過程を無視することで生じるアーティファクトを除去するために、物質のＸ線吸収係数の波長依存性など多重Ｘ線和の観測過程を考慮した確率分布を設定しＸ線ＣＴ画像処理を行う。ここで、物質のＸ線吸収係数の波長依存性を考慮する際には、波長ごとのＸ線吸収係数を全て未知変数として推定すると、問題の不良設定性のために信頼のおける解が得られないことから、いくつかの既知の物質が特定の波長依存性をもったＸ線吸収係数を有することを表現する事前分布を利用することにし、再構成画像の各画素が既知の物質クラスのうちのいずれの物質クラスから構成されるかを推定するという、より易しい推定問題に変更することで、より信頼のおける解を得る。すなわち、波長ごとに正の実数変数であるＸ線吸収係数を推定する代わりに、波長に依存しない正の実数変数であるＸ線吸収係数密度と離散変数である物質クラスを推定する。この各画素の物質クラスは直接的に観測できない隠れ状態変数、すなわち物質クラスの変数を隠れた確率変数として扱う。
　この隠れ状態変数の物質クラスに対して、空間的滑らかさの確率的な制約を課すことで、波長依存のＸ線吸収係数の自由度を確率的に制約することとした。その上で隠れた変数をもつ確率分布の推論として、ベイズ推論を近似的に行うことで安定した精度の高い推定（収束が保証され、観測ノイズに対してあまり敏感に反応しない安定した推定）を行うのである。また、本発明におけるベイズ推定によれば、Ｘ線吸収係数と再構成画像の各画素の領域が属する物質クラスの両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。ベイズ推定により、各画素における物質クラスの所属確率と各画素のＸ線吸収係数密度を確率的に表現できるので、決定論的に物質クラスを選択する場合に比べてより柔軟性をもったきめ細やかなＸ線吸収係数の表現ができる。
　本発明のＸ線ＣＴ画像処理方法では、物質クラスの変数を確率変数として扱うことで、既知の物質クラスで表現されない物質に対しても、その重み付き和による柔軟な表現が可能になる。また、本発明のＸ線ＣＴ画像処理方法では、物質クラスに対して事前知識として空間的滑らかさなどの制約を課すことで不良設定性を抑制することができる。 また、本発明のＸ線ＣＴ画像処理方法では、近似ベイズ推論を行うことで、収束性が保証され、かつ、推定の不安定性(オーバーフィッティング)を抑えることができる。
　すなわち、本発明のＸ線ＣＴ画像処理方法を用いることにより、物理過程により忠実でかつ精度の高い安定した推定が行えるのである。
　本発明において、物質クラスは、組織クラスを包含する意味で用いている。物質クラスは、離散変数で表現される。物質クラスの存在割合（密度）とは、人体などの撮像対象の物質クラス分布に関する事前知識について、どのような物質クラス（筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど）がどの程度の割合で分布するのを示すもので、非負のスカラーのパラメータである。
　物質クラス毎に波長間で固定されたＸ線吸収係数比を表現するパラメータとは、確率的に表現されたＸ線吸収係数において、波長間で吸収係数の比は固定と仮定し、それを事前知識の確率分布に用いて、波長依存のＸ線吸収係数の自由度を制約するものである。
　また、物質クラスが空間的に連続する度合いを表現するパラメータは、例えば、空気の物質クラスは画素同士の連結がしやすい（空気中に他の組織は入らない）など空間的に連続する度合いが高いという制約を課すことができるので、物質クラスが空間的に連続する度合いを事前知識の確率分布に用いて、波長依存のＸ線吸収係数の自由度を制約するものである。
　上記のＸ線ＣＴ画像処理方法によれば、人体などの撮像対象の物質クラス分布に関する事前知識について、どのような物質クラス（筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど）がどの程度の割合で分布するか、また、それぞれの物質クラスがどの程度のＸ線吸収係数を有するか、また、それぞれの物質クラスがどの程度、空間的に連続して分布しやすいかの分布情報を確率分布の形で表現して、事後確率の期待値によるベイズ推定を用いて、画像再構成や物質推定を行う。
　なお、物質クラスに関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と物質クラスの推定精度のさらなる向上を図ることが可能である。
　本発明のＸ線ＣＴ画像処理方法における多重Ｘ線和にＸ線波長に関する和を含める場合を想定する。
　この場合、各画素の波長依存性を有するＸ線吸収係数に対する事前分布を、物質クラス毎およびＸ線波長毎のＸ線吸収係数の分布を表現するパラメータを介して特徴づけることによって、どのようなＸ線吸収係数が生じやすいかを確率的に限定するものである。
　他の観点から、本発明のＸ線ＣＴ画像処理方法における多重Ｘ線和にＸ線波長に関する和を含める場合、各画素の波長依存性を有するＸ線吸収係数は、物質クラスに依存しない各画素のＸ線吸収係数密度と、物質クラス毎に定められているＸ線波長間のＸ線吸収係数比の積で表現される。
　そして、各画素の波長依存性を有するＸ線吸収係数は、物質クラス毎に定められている、Ｘ線波長間のＸ線吸収係数比を表現するパラメータで特徴付けられる。
　そして、各画素の波長依存性を有するＸ線吸収係数の事前分布が物質クラスを条件とした各画素のＸ線吸収係数密度に対する条件付き事前分布と物質クラスに対する事前分布によって表現される。
　すなわち、Ｘ線吸収係数を各画素のＸ線吸収係数密度と物質クラス依存波長依存の定数の積で表現する場合、Ｘ線吸収係数の事前分布は、各画素のＸ線吸収係数密度に対する事前分布と物質クラスに対する事前分布から一意に決定されることになる。
　ここで、本発明のＸ線ＣＴ画像処理方法において、上記の物質クラス毎に定められているＸ線吸収係数密度の分布が、混合ガンマ分布で表現される。また、物質クラスが、出現する度合いを表現するパラメータと物質クラスが空間的に連続する度合いを表現するパラメータとで特徴付けられるボルツマン分布で表現される。
　これは、Ｘ線吸収係数の事前分布の表現についての一例である。
　物質クラスを条件としたＸ線吸収係数密度に対して混合ガンマ分布を仮定し、物質クラス変数に対してボルツマン分布を仮定することによって、Ｘ線吸収係数の事前分布が決定されるようにする。
　物質クラス毎に定められているＸ線吸収係数密度の分布として（混合）ガンマ分布で表現されるとした理由は、その分布形状が単峰な分布を表現でき、かつその確率変数のとる値が必ず正の値となることを保証できることに着目したものである。また、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程において検出器において観測される信号の分布としてガウス分布で表現されるとした理由は、光子散乱や量子化雑音などの雑音により、理想値からずれた値を示すことと、投影したＸ線が多重である場合にも計算が用意であることを考慮して、これらの雑音をガウス分布でよく表現されると仮定したものである。
　また、物質クラスが、出現する度合いを表現するパラメータと前記物質クラスが空間的に連続する度合いを表現するパラメータとで特徴付けられるボルツマン分布で表現されるとした理由は、同じ組織が空間的に集まりやすい(組織別に集まりやすさの調整が可能)ことと、標準とされる各組織の割合を表現しやすいことに着目したものである。
　また、本発明のＸ線ＣＴ画像処理方法において、ベイズ推定は、物質クラス、Ｘ線吸収係数密度、物質クラスのサブクラス変数、それぞれに関する事後分布の推定が行われる。
　ここで、サブクラス変数は、物質クラスが定まった時に、どのようなＸ線吸収係数密度をとりやすいかという確率分布を２つのガンマ分布の和で表現される場合に、いずれのガンマ分布に所属するかを指定する離散変数である。
　物質クラスが与えられた条件下でのＸ線吸収係数の事前分布は、特定の値(たとえば平均値)を取りやすいと考えられるものの、その平均値に対して対称な確率分布となるとは限らなかったり、確率分布の尖度が大小様々であったりすることが予想できることから、このような複雑な事前分布を表現できるように２つの確率分布の重み付き和である混合確率分布で表現したものである。さらにＸ線吸収係数は正の値しか取らないため、この混合確率分布をガンマ分布の混合確率分布で表現されると仮定したものである。物質クラスのサブクラス変数は、２つ分布の混合確率分布においてどちらの確率分布に所属するかを指定するために便宜的に導入される０又は１の二値変数である。
　本発明のＸ線ＣＴ画像処理装置は、Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理装置であって、Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力する入力手段と、観測に関する確率分布と、物質クラスという再構成画像の各画素にどの物質が存在するかを表す変数を介して表現されるＸ線吸収係数に関する事前分布と、物質クラスに関する事前分布のもとで、Ｘ線吸収係数及び物質クラスに関する事後分布の期待値によるベイズ推定でＸ線吸収係数及び物質クラスを推定する推定手段を備える。ここで、観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を通じてどのような投影像が得られやすいかを表す確率分布として設定される。また、多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである。
　上記のＸ線ＣＴ画像処理装置によれば、物質クラスの存在割合（密度）を表現するパラメータを用いて、物質クラスの変数を確率変数として扱うことで、予め想定していた物質クラスから外れる物質に対しても、予め想定した物質クラスの重み付き和によって表現することで、近似表現が可能になる。また、物質クラスが空間的に連続する度合いを表現するパラメータを用いることにより、物質クラスに対して確率分布として空間的滑らかさの制約を課すことで不良設定性を抑制することができる。 ベイズ推定を近似的に行うことで、収束性が保証され、かつ、推定精度が向上することが期待できる。
　また、本発明の他の観点によれば、Ｘ線吸収係数に関する確率分布をおいて統計推定を行うＸ線ＣＴ画像処理プログラムであって、コンピュータに、以下の１）～２）のステップを実行させるためのプログラムである。
１）Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力するステップ
２）上記のＸ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、Ｘ線吸収係数及び物質に関する事後分布の期待値でＸ線吸収係数及び物質を推定するステップ
　ここで、上記の確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現され、多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである。
　上記のＸ線ＣＴ画像処理プログラムによれば、物質の存在割合（密度）を表現するパラメータを用いて、物質の変数を確率変数として扱うことで、既知の物質で表現されない物質に対しても、その重み付き和による柔軟な表現が可能になる。また、物質が空間的に連続する度合いを表現するパラメータを用いることにより、物質に対して確率分布として空間的滑らかさの制約を課すことで不良設定性を抑制することができる。 ベイズ推定を近似的に行うことで、収束性が保証され、かつ、推定の不安定性を抑えることができる。
　また、本発明は、上記のＸ線ＣＴ画像処理プログラムを搭載したＸ線ＣＴ画像処理装置を提供することができる。
　本発明によれば、従来技術を用いた画像再構成と比べて、Ｘ線吸収係数を確率的に表現することで、より柔軟な表現を可能とし、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、ビームハードニングアーティファクトを低減できるといった効果を有する。
　また、本発明によれば、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程に関する確率分布を設定し統計推定を行うことで、一本の光線で代表させることによるバイアスを低減できるといった効果を有する。
本発明のＸ線ＣＴ画像処理方法の処理フロー Ｘ線ＣＴの説明図 更新の説明図 金属と骨のＸ線吸収係数の確率分布図 実施例１の計算機実験での画像再構成を示す図 実施例２の計算機実験での画像再構成を示す図 平行ビームにおける投影の幾何の説明図 平行ビームにおける投影と画素化の説明図 回転時の投影の幾何の説明図
　以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。
　本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、Ｘ線吸収係数に関する事前分布を設定し統計推定を用いて画像の再構成を行う。特に、Ｘ線スペクトルのエネルギー分布（波長分布）を考慮し、エネルギー依存（波長依存）のＸ線吸収係数を統計モデルに基づいて推定することにより、ビームハードニングアーティファクトを解消する。
　推定する変数増加による不良設定性に対処するため、適当な事前知識を用いて解に何らかの制約を課すことで解決を試みる必要がある。本発明では基底関数表現を用いて、その基底を既知の物質のエネルギー依存（波長依存）のＸ線吸収係数とし、その重みはどの物質クラスに所属するかを表現する離散変数と、Ｘ線吸収係数密度を表現する非負のスカラーの連続変数の積で与えられるとし、物質クラスを条件としたＸ線吸収係数密度の事前分布と物質クラスの事前分布を課して解に制約を課す。
　また、本発明では、人体組織を対象とする再構成画像の場合、人体組織は脂肪，筋肉，骨といった限られた物質から構成され、それぞれの物質のＸ線吸収係数の大まかな分布はあらかじめ分かっているので、その知識を事前分布として表現して推定に利用する。 この推定のためには、各画素がどの物質クラスに所属するかの推定を行う必要があるが、それぞれの物質は空間的な連続性を有しており、各物質が人体に占める割合も大まかに分かっているという状況を想定できることから、物質クラスを表す離散変数に空間的な滑らかさ制約や占有割合などの事前知識を導入することにより不良設定性の抑制を行う。
　投影Ｘ線の観測に関する観測過程については、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を表現するものであり、Ｘ線が有限の広がりを持つ波長域（スペクトル分布）を有することによって、波長が異なる複数のＸ線が検出器に入射することを表現する。
　再構成画像に関する確率分布については、再構成画像の各画素の領域において、Ｘ線吸収係数がどのような値をとりやすいかを表現するかを特徴付けるものであって、各物質クラスと各サブクラス変数の組み合わせがＸ線波長毎にどのようなＸ線吸収係数密度を取りやすいかを表現する混合ガンマ分布と、各物質クラスと各サブクラス変数の組み合わせがどの程度、出現しやすいかを表現するパラメータと、各物質クラスと各サブクラス変数の組み合わせがどの程度、空間的に滑らかに連続しやすいかを表現するパラメータからなるボルツマン分布で表現する。
　そして、Ｘ線吸収係数と物質クラスに関する事後分布の期待値推定による統計推定を行う。画像再構成において、この事後分布を求めることが必要であるが、各画素で定義される物質クラスという高次元の隠れ離散状態変数に関する和計算が含まれるため、解析的に計算を実行することは困難である。よって、本発明では、近似手法を適用することでこの計算困難さを克服する（図１のフローチャートを参照）。
（問題の定式化）
　Ｘ線ＣＴ画像のように、様々な方向から投影されたＴ個の投影データを、Ｄ＝｛Ｙ（１），・・・，Ｙ（Ｔ）｝と表すことにする。各々のデータＹ（ｔ）は、ｔ番目の投影によって検出器で検出されるデータの集合であり、Ｙ（ｔ）＝｛ｙ１ （１），・・・，ｙＩ （ｔ）｝となる。但し、Ｔは投影数、Ｉは検出器の数、ｙｉ （ｔ）はｉ番目の検出器で検出された光子の数を表す。
　そして、Ｘ線は物質を透過する際に指数的に減衰することから、下記数式（１）のように表すことができる。
Figure JPOXMLDOC01-appb-M000001
　ここで、ｘｊは観測対象を撮像対象のＸ線吸収係数をラスタスキャンして得られるＪ次元ベクトルｘ＝｛ｘ１，・・・，ｘＪ｝のｊ番目の画素のＸ線吸収係数である。ｂｉ （ｔ）は、Ｘ線源から放出される光子数（物体が何も置かれていないときに観測され得る光子数）を表す。また、ｌｉｊ （ｔ）は角度θ（ｔ）から投影された際のｉ番目の検出器で検出される投影線とｊ番目の画素が交差する距離に相当するものであり、ｌｉｊ （ｔ）xｊがｔ番目の投影のｉ番目の検出器に入射するＸ線に対するｊ番目の画素領域の実効的なＸ線吸収係数を表す（図２を参照）。
　ここで、Ｘ線のエネルギースペクトルをＥ通りに離散化して表現することを考える。入射Ｘ線の内、エネルギーｅ（ｅ＝１，・・・，Ｅ）に属するＸ線の光子数をＮｅ （ｔ）とすると、物質に透過した後に観測される光子数の平均値は、下記数式（２）のように表すことができる。
　ただし、ｘｊｅはエネルギーｅのＸ線に対するｊ番目の画素のＸ線吸収係数である。
Figure JPOXMLDOC01-appb-M000002
（Ｘ線吸収係数のエネルギー依存性について）
　ｊ番目の画素におけるエネルギーｅでのＸ線吸収係数ｘｊｅを、エネルギー毎、画素毎に自由に決定すると不良設定性が生じる。実際には、いずれの物質も高エネルギーになるにつれてＸ線吸収係数が小さくなるという物質間で共通した性質を有する。つまり、異なるエネルギー間のＸ線吸収係数の比に特定の傾向が存在する。そこで、物質種を区別する物質クラスを導入して、Ｘ線吸収係数を制約した形で表現する。
　観測対象の物質がＣ種類の物質クラスに分類されると仮定する。ここで、ｊ番目の画素の物質クラスを変数ｚｊを用いて表現する。ただし、ｚｊは、ｚｊ＝｛ｚｊ１，・・・, ｚｊＣ｝と表される。各要素ｚｊｃは０もしくは１をとる二値変数であり、ｊ番目の画素が物質クラスｃに属している場合に１をとり、他の場合は０をとる。この物質クラスに依存して一意に決まるエネルギー間のＸ線吸収係数の比を表現する変数として、変化率ｒｃｅ（＞０）を定義する。ただし、変化率ｒｃｅは、ｅ＝１の時ｒｃｅ１＝１と正規化されているものとする。ｅ＝１におけるｊ番目の画素のＸ線吸収係数をｘｊ（＞０）で表現する。以上の表現を用いて、Ｘ線吸収係数ｘｊｅは下記数式（３）のように表すことができる。
Figure JPOXMLDOC01-appb-M000003
　上記数式（３）の表現を用いた場合、各画素で推定する必要がある変数は、スカラーの連続変数となるｘｊとＣ種類のいずれかの値をとる離散変数ｚｊの２つだけとなる。
（観測モデルについて）
　各検出器で観測される光子数の理想値は、上述の数式（１）に示した通りである。しかし、実際には光子散乱や量子化雑音などの雑音によって理想値からずれた値を示すことがわかっている。これらの雑音が、ガウス分布でよく表現されると仮定すると、各検出器における観測データは下記数式（４）のように表すことができる。但し、Ｘ＝｛ｘ１，・・・，ｘＪ｝、Ｚ＝｛ｚ１，・・・，ｚＪ｝である。
　また、観測が、検出器毎、投影毎に独立であると仮定すると、投影データの集合Ｄ（観測データ）は下記数式（５）のように表すことができる。但し、Ｔは投影数、Ｉは検出器の数である。ここで、ガウス分布の標準偏差を表すσは、ショット雑音の性質を考慮して、下記数式のように観測ｙｉ （ｔ）に依存する関数として良い。
Figure JPOXMLDOC01-appb-M000004
Figure JPOXMLDOC01-appb-M000005
（事前分布について）
　事前分布については、観測対象の物質種ごとのＸ線吸収係数に関する事前知識を表現する。画素ｊでのＸ線吸収係数ｘｊは画素ｊにおける物質クラスｚｊｃのみに依存して決定されると仮定すると、事前分布ｐ（Ｘ｜Ｚ）は下記数式（６）のように表すことができる。
Figure JPOXMLDOC01-appb-M000006
　物質クラスｚｊｃが与えられたもとでは、Ｘ線吸収係数ｘｊは特定の値を取りやすいことが予想できる。ここでは事前分布ｐ（ｘｊ｜ｚｊ）は、分布の非対称性や様々な異なる尖度を表現できるよう２つの分布の混合分布として下記数式（７）のように表すことができる。但し、変数ｂｊは二値変数である。また、ベクトルＢをＢ＝｛ｂ１，・・・，ｂＪ｝で定義する。
Figure JPOXMLDOC01-appb-M000007
　画素毎の物質クラスｚｊｃと変数ｂｊｄが与えられたとき、Ｘ線吸収係数ｘｊはガンマ分布に従うと仮定すると下記数式（８）のように表すことができる。但し、ｕｃｄ，ｖｃｄはガンマ分布のパラメータであり、平均ｕｃｄ／ｖｃｄ、分散ｕｃｄ／ｖ２ ｃｄである。これらのパラメータは物質クラスに関する経験的な知識により定める。
Figure JPOXMLDOC01-appb-M000008
　また、事前分布において変数ｚｊｃは下記数式（９）を満たす、０もしくは１をとる二値変数とする。ガンマ分布を用いることでＸ線吸収係数が負にならない制約を表現できる。
Figure JPOXMLDOC01-appb-M000009
　ここで、変数Ｚ，Ｂの事前分布に関しては変数間の独立性を仮定せず、下記数式（１０）で表されるように、ボルツマン分布に従うと仮定する。但し、Ａは正規化項であり、エネルギーＨ（Ｚ，Ｂ）は下記数式（１１）のように定義する。下記数式（１１）では、２次元平面の近傍として縦横に隣接する４画素を考えており、Ｇｃｄ ｓｅｌｆ，Ｇｃｄ ｉｎｔｅｒはそれぞれ非負の定数である。
Figure JPOXMLDOC01-appb-M000010
Figure JPOXMLDOC01-appb-M000011
（事後分布について）
　本発明のＸ線ＣＴ画像の再構成では、Ｘ線吸収係数Ｘ、物質クラスＺ，変数Ｂは観測データＤの事後分布として推定される。Ｘ線吸収係数Ｘ、物質クラスＺ、変数Ｂの事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）はベイズの定理により下記数式（１２）のように表すことができる。
Figure JPOXMLDOC01-appb-M000012
（ベイズ推定）
　ベイズ推定では、上記数式（１２）で表される事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）の期待値が必要になるが、この事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）を解析的に計算するのは困難であることに鑑みて、試験分布ｑ（Ｘ，Ｚ，Ｂ）により事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）を近似する。試験分布ｑ（Ｘ，Ｚ，Ｂ）は、上記数式（１２）の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。そして、試験分布ｑ（Ｘ，Ｚ，Ｂ）と事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）の誤差をＫＬ（Kullback－Leibler）距離で評価し、ＫＬ距離を最小化する試験分布ｑ（Ｘ，Ｚ，Ｂ）を算定する。ここで、ＫＬ距離は下記数式（１３）のように表すことができる。
　但し、＜・＞ｑ（Ｘ，Ｚ，Ｂ）は分布ｑ（Ｘ，Ｚ，Ｂ）に関して積分計算を行うことを表す演算子である。ＫＬ距離は常に非負であり、ｑ＝ｐのときにのみ０（ゼロ）になる。＜・＞ｑ（Ｘ，Ｚ，Ｂ）はｑ（Ｘ，Ｚ，Ｂ）に関して期待値をとることを表している。
Figure JPOXMLDOC01-appb-M000013
　この最適化を容易にするため、試験分布ｑ（Ｘ，Ｚ，Ｂ）がその周辺分布ｑ（Ｘ）とｑ（Ｚ，Ｂ）の積で与えられる、すなわち変数Ｘと｛Ｚ，Ｂ｝の間に独立性があることを仮定する。また、ｑ（Ｘ），ｑ（Ｚ，Ｂ）の各要素について画素間で独立であると仮定し、さらにｑ（Ｘ）の各要素ｑ（ｘｊ）はガンマ分布に従うとする。
　このときｑ（ｚｊ、ｂｊ）は、確率変数ｚｊ, ｂｊが離散変数であるため必然的に多項分布となる。このことから、試験分布ｑ（Ｘ，Ｚ，Ｂ）は下記数式（１４），（１５）のように表すことができる。但し、試験分布ｑ（Ｘ），ｑ（Ｚ，Ｂ）について同時に最適化するのは困難であるため、下記数式（１６），（１７）のように交互に最適化を行うものとする。
Figure JPOXMLDOC01-appb-M000014
Figure JPOXMLDOC01-appb-M000015
Figure JPOXMLDOC01-appb-M000016
Figure JPOXMLDOC01-appb-M000017
　ここで、ｑ（ｘｊ）のパラメータαｊ，βｊの最適化にはＳＣＧ（Scaled Conjugate Gradient）法を用いた。試験分布ｑ（Ｘ）とｊ番目の画素以外の       ｑ（ｚｊ´，ｂｊ´) (但し、ｊ´はｊと異なる) が固定のもとで、ＫＬ距離を最小化するｑ（ｚｊ，ｂｊ)は解析的に求めることができる。この解析解はｊの近傍画素のｑ（ｚｊ´，ｂｊ´）（ｊ´∈η（ｊ）)に依存する。そこで、近傍画素の分布を変更せずにｑ（Ｚ，Ｂ）を更新するために、図３に示すように市松模様となる画素の部分集合を交互に最適化している。すなわち、パラメータｑ（Ｚ，Ｂ）は全画素を同時に更新するのではなく、図３に示す市松模様のパターンで画素の部分集合を交互に最適化している。図３（ａ）のグレー部分の画素に該当するｑ（ｚｊ，ｂｊ）を更新したのち、図３（ｂ）のグレー部分の画素に該当するｑ（ｚｊ，ｂｊ）を更新する。
　以上、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおける処理内容について説明した。以下の実施例では、本発明の有効性を調べるために計算機実験を行い、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムの画面再構成を、従来技術による画面再構成と比較して評価する。
　下記の評価実験から、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムの有用性が理解できるであろう。
　実施例では、再構成画像と真の画像と誤差を測れるように、適当な真の画像を用意し、それに対してＸ線投影をシミュレートして人工的に投影像（サイノグラム）を生成した。Ｘ線投影は、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮するために入射するＸ線のエネルギーを６０ｋｅＶ，７０ｋｅＶ，８０ｋｅＶ，９０ｋｅＶ，１００ｋｅＶの５通り（それぞれのフォトン数は、６×１０４，２．５×１０４，１．５×１０４，７×１０３，５×１０３）に分割し、それぞれで物質固有のＸ線吸収係数に従って減衰させた。検出器においてはその減衰したＸ線の光子にポアソンノイズに従って生成したノイズ（雑音）が加えられたものが観測されるとした。
　以下の実施例１，２において、検出器の個数は９５（個）、画像の解像度は６４×６４（画素）、投影は１°から１８０°の間で１°間隔としている。再構成画像の推定の良さを定量的に評価するために、下記数式（１８）で表されるＰＳＮＲ（Peak Signal to Noise Ratio) を用いた。
Figure JPOXMLDOC01-appb-M000018
　上記数式（１８）において、ＭＡＸは画素値の最大値であり、ＭＳＥは真の画像と再構成画像との平均２乗誤差である。ＰＳＮＲが高いほどＭＳＥが小さい。すなわち、より真の画像に近いといえる。以下、ＦＢＰ法、Ｘ線エネルギー分布を考慮しない事前分布をガウス分布としてベイズ推定により行う方法（以下、「従来方法」と略する）、Ｘ線エネルギーを３つに離散化した本発明のＸ線ＣＴ画像処理方法（以下、「本発明方法」と略する）の３通りで画像再構成を行い、それぞれ比較した結果を示す。
　なお、金属と骨の事前分布は図４に示す通りである。図４は、金属と骨の事前分布と真の画像の６０ｋｅＶにおけるＸ線吸収係数を示している。金属と骨のＸ線吸収係数は、それぞれ０．７，０．４を平均とするガンマ分布で表現される確率分布である。真の画像のＸ線吸収係数は、それぞれのＸ線エネルギーにおいて、事前分布の平均値より５％小さい値であるとしている。
　実施例１では、金属を４つ含むファントムに対する画像再構成についての計算機実験について図５を参照して説明する。
　図５（ａ）は真の画像、（ｂ）はＦＢＰ法による再構成画像、（ｃ）は単純ベイズ法による再構成画像、（ｄ）は本発明方法による再構成画像である。
　実施例１では、図５（ａ）に示す真の画像に対して、ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれの方法を用いて画像再構成実験を行った。入射Ｘ線の光子数の合計を１．１２×１０５とし、本発明のＸ線ＣＴ画像処理方法において、Ｘ線のエネルギーを６０ｋｅＶ，８０ｋｅＶ，１００ｋｅＶの３つに離散化し、それぞれに属する光子数を７．５×１０４，３×１０４，１０．７×１０４であるとした。繰り返し回数は比較のため共に２５回と設定した。推定結果を、図５（ｂ）～（ｄ）に示す。
　また、ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれにおけるＰＳＮＲを下記表１に示す。ＰＳＮＲでは、ＦＢＰ法、単純ベイズ法と比べて本発明方法が大きく、本発明方法がより精度の高い推定を行えることが理解できる。
Figure JPOXMLDOC01-appb-T000001
　実施例２では、金属インプラントを含む歯のファントムに対する画像再構成についての計算機実験について図６を参照して説明する。
　図６（ａ）は真の画像、（ｂ）はＦＢＰ法による再構成画像、（ｃ）は単純ベイズ法による再構成画像、（ｄ）は本発明方法による再構成画像である。
　実施例２では、図６（ａ）に示す真の画像に対して、ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれの方法を用いて画像再構成実験を行った。入射Ｘ線の光子数の合計を１．１２×１０５とし、本発明のＸ線ＣＴ画像処理方法において、Ｘ線のエネルギーを６０ｋｅＶ，８０ｋｅＶ，１００ｋｅＶの３つに離散化し、それぞれに属する光子数を７．５×１０４，３×１０４，１０．７×１０４であるとした。繰り返し回数は比較のため共に２５回と設定した。推定結果を、図６（ｂ）～（ｄ）に示す。
　ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれにおけるＰＳＮＲを下記表２に示す。実施例１の実験結果と同様に、ＰＳＮＲでは、ＦＢＰ法、単純ベイズ法と比べて本発明方法が大きく、本発明方法がより精度の高い推定を行えることが理解できる。
Figure JPOXMLDOC01-appb-T000002
　以上のように、Ｘ線スペクトルのエネルギー分布を考慮し、エネルギー依存性（波長依存性）のあるＸ線吸収係数を統計モデルに基づいて推定することにより、ビームハードニングの影響による様々なアーティファクトを軽減できることを示した。実施例１，２の計算機実験としてはＸ線エネルギーを３つに分割した場合のみを行ったが、このように実際のエネルギーの分割数より少ない数で離散表現しても、アーティファクトが大きく軽減できることがわかる。さらにエネルギーの分割数を増やすことでより良い推定を行うことが可能であろう。
　次に、投影Ｘ線の観測に関する観測過程について、多重Ｘ線和が、検出器におけるＸ線の検出面が空間的に有限の広がりを有すること、検出器におけるＸ線の検出時間が有限の幅を有すること、Ｘ線が有限の広がりを持つ波長域（スペクトル分布）を有すること、によって、空間的な位置と方向、波長が異なる複数のＸ線が検出器に入射することを表現できることを説明する。
　まず、平行ビームを例にとり、図７を参照して説明する。図７のように２次元平面で平行ビームによる投影を想定する。ここで、Ｘ線の放出される方向をｖ軸、それに直交する方向をｕ軸ととる。この座標系で以下のように変数を定義する。
　位置（ｕ，ｖ）における物質の波長λのＸ線の吸収係数をｘ（ｕ，ｖ，λ）とする。但し、物質は、ｖ∈Ｌに含まれるとする。
　ｉ番目の検出器の検出面Ｓｉに、単位時間あたり入射する波長λのフォトン数をｙｉ（λ）とする。平行ビームのＸ線源が、単位時間、単位面積あたり放出する波長λのフォトン数をｄ（λ）、波長λのフォトン１個が励起する電子の平均的な数をρ（λ）とする。
　このとき単位時間にi番目の検出器で観測される電流Ｙｉの期待値は、下記数式（１９）のように表される。さらに、ある一定時間の電流を貯めてそれを観測する場合は、時間に関する積分を行うことになる。
Figure JPOXMLDOC01-appb-M000019
　上記数式（１９）は逆投影法の仮定する線積分（上記のｖに関する積分）に帰着できないことから、検出面の空間的な広がりや、Ｘ線のスペクトル分布の拡がりが存在する場合は線積分を仮定する推定手法を用いるとアーティファクトが生じることがわかる。同様に、検出器が受光している時間における回転を考慮する場合も線積分に帰着できない。上記の多重積分の形では計算が難しいので積分を和で近似することを考える。これが多重Ｘ線和を考えることになることを示す。まず、上記数式（１９）において特定の波長λにのみ注目し、下記数式（２０）におけるｕ, ｖに関する積分を和に置き換える。
Figure JPOXMLDOC01-appb-M000020
　計算機上で上記のｕ, ｖに関する積分を近似的に表現するために離散化を行う（図８を参照）。　空間をｕ軸、ｖ軸にそって離散化すると、上記数式（２０）におけるｖに関する積分は下記数式（２１）のように近似できる。
Figure JPOXMLDOC01-appb-M000021
　ここで、ｘｊはｊ番目の画素のＸ線吸収係数であり、Ｉ（ｕ）はｕ軸上のｕの位置にある画素のインデックスの集合、ΔＬは画素の１辺の長さを表す。画素の１辺の長さをΔＬ＝１とすると、上記数式（２１）は下記数式（２２）として表現できる。
Figure JPOXMLDOC01-appb-M000022
　したがって、画素の１辺の長さが１のときは、ｙｉは下記数式（２３）となり、さらに同じ検出器に入る光線に関する積分（ｕに関する積分）を和で近似すると、下記数式（２４）のように表現できる。
Figure JPOXMLDOC01-appb-M000023
Figure JPOXMLDOC01-appb-M000024
　上記数式２４において、１／Ｎは積分を短冊で切って近似するときの短冊の幅を示している。
　これはインデックスｋごとに線積分で表現される投影の和を考えていることに相当する。
　以上から、検出器の空間的拡がりを考慮する場合においては多重のＸ線を考慮する必要があることがわかる。
　上記では、Ｘ線が画素に対してまっすぐ入射する場合を想定している。しかしながら、一般のＸ線ＣＴ装置などのように、Ｘ線源と検出器が対向しながら回転するような装置（図９を参照）では、画素に対して斜めに入射する状況を想定する必要がある。
　この場合、投影についての下記数式（２４）を、下記数式（２５）のように一般化することで表現できる。
Figure JPOXMLDOC01-appb-M000025
Figure JPOXMLDOC01-appb-M000026
　ここで、ΔＮ＝１／Ｎであり、ｒ（ｋ）は光線ｋの検出平面上の位置を表し、｛ｋ｜ｒ（ｋ）∈Ｓｉが検出器ｉに到達する光線の集合を表す。Ｉｋは光線ｋが交差する画素のインデックスの集合を表し、ｌｋｊは光線ｋが画素ｊと交差する長さを表している。上記数式（２５）は上記数式（２６）のｌｋｊを１とした場合と考えることができる。
　上記数式（２６）は、平行ビームのほか、ファンビーム、コーンビームでも光線の向きを考えて、ｌｋｊを計算することで同様の式が成立する。
　上記数式（２６）からわかることは、検出面の拡がりを考えた投影の観測モデルを考えた時には、検出面内を通る光線ごとの投影を考える必要がある点である。
　上記数式（２６）の中のｌｋｊｘｊの和を、検出器ｉや光線ｋに関してまとめて表記するならば、Ｌｘのような行列演算で表現できる。検出面の拡がりを考えなかった場合には、行列Ｌは（検出器の数）×（再構成画像の画素数）であったものが、検出面の拡がりを考えた場合には、行列Ｌは（検出器を通る光線の数）×（再構成画像の画素数）となる。光線の数は、多ければ積分を精度よく近似できる一方で、計算量が増す問題がある。
　次に、投影中の回転（時間に関する積分）を想定する。まず時刻ｔにおいて観測されたフォトン数をｙｉ，ｔとすると、下記数式（２７）で表現できる。
Figure JPOXMLDOC01-appb-M000027
　検出器の検出時間（露光時間）の間に計測されたフォトン数は、ｙｉ，ｔを時刻ｔに関して積分したものとなるので、下記数式（２８）で表現できる。
Figure JPOXMLDOC01-appb-M000028
　上記数式（２８）について、時刻に関する積分を和で置き換えると、下記数式（２９）のように表現できる。
Figure JPOXMLDOC01-appb-M000029
　前述のように、検出面の拡がりを考慮すると上記数式（２８）のように各時刻ｔで検出器に入射するフォトンｙｉ，ｔが複数の光線束の和で表現される。時間に関する積分は、さらなる和を含めるため新たな光線束に関する和を考えることに等価である。そこで、光線インデックスの集合Ｓｉを時刻に依存するようＳｉ，ｔとし、それぞれの光線ｋが通る画素インデックスの集合Ｉｋを改めて定義しなおすと、検出時間（露光時間）Ｔで計測されたフォトン数は、下記数式（３０）で表現できる。
Figure JPOXMLDOC01-appb-M000030
　上記数式（３０）を波長に関する積分を和で置き換えると、下記数式（３１）で表現できる。
Figure JPOXMLDOC01-appb-M000031
　以上をまとめると、検出面の拡がりと投影中の回転,Ｘ線の波長依存性を考慮した場合、下記数式（３２）で表現できる。ここで、ｂｍ＝ｄｍρｍである。
Figure JPOXMLDOC01-appb-M000032
　上記数式（３２）のように、ｍやｋに関する和が指数関数の外側でとられることを多重Ｘ線和と呼ぶ。これは、指数関数の外側に和（もしくは積分）があることを想定しない線積分では表現できない。このような多重Ｘ線和のある観測過程一般に対して、本明細書で開示したベイズ推定に基づくアルゴリズムを適用することができる。
　本発明は、ビームハードニングアーティファクトを低減するＸ線ＣＴ装置として、医療用のみならず産業用のＸ線ＣＴ装置に有用である。また、ベルトコンベアによる移動や回転が無視できないベルトコンベア型Ｘ線ＣＴ装置や歯科用Ｘ線ＣＴ装置としても有用である。
 

Claims (8)

  1. 　Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理方法であって、
    　Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
    　前記Ｘ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、前記統計推定が前記Ｘ線吸収係数及び前記物質に関する事後分布の期待値によって前記Ｘ線吸収係数及び前記物質を推定するステップと、
    　を備え、
    　前記Ｘ線投影画像の観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現され、
    　前記多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである、
    　ことを特徴とするＸ線ＣＴ画像処理方法。
  2. 　前記多重Ｘ線和に、Ｘ線波長に関する和を含める場合、
    　各画素の波長依存性を有するＸ線吸収係数に関する事前分布は、
    　前記物質毎およびＸ線波長毎のＸ線吸収係数の分布を表現するパラメータを用いて表現される、
    　ことを特徴とする請求項１に記載のＸ線ＣＴ画像処理方法。
  3. 　前記多重Ｘ線和に、Ｘ線波長に関する和を含める場合、
    　各画素の波長依存性を有するＸ線吸収係数は、
    　波長に依存しない各画素のＸ線吸収係数密度と、画素に依存せず前記物質毎に定められるＸ線波長間のＸ線吸収係数比の積として表現され、
    　前記物質毎に定められる、Ｘ線波長間のＸ線吸収係数比を表現するパラメータを用いて表現され、
      各画素の波長依存性を有するＸ線吸収係数の事前分布が前記物質を条件とした各画素の前記Ｘ線吸収係数密度に対する条件付き事前分布と前記物質に対する事前分布によって表現される、
    　ことを特徴とする請求項１に記載のＸ線ＣＴ画像処理方法。
  4. 　前記物質に対する事前分布が、それぞれの物質の出現する度合いを示すパラメータとそれぞれの物質が空間的に連続する度合いを示すパラメータとで表現される事前分布で設定されることを特徴とする請求項３に記載のＸ線ＣＴ画像処理方法。
  5. 　前記物質毎に定められている前記Ｘ線吸収係数密度の条件付き事前分布が各物質のサブクラス毎に定められている前記Ｘ線吸収係数密度の条件付き事前分布と各物質のサブクラスの事前分布の積を各物質のとりえる全てのサブクラスに関して和をとった混合確率分布で表現され、前記各物質のサブクラスの事前分布が、各物質のサブクラスの出現する度合いを示すパラメータと前記各物質のサブクラスが空間的に連続する度合いを示すパラメータとで設定されることを特徴とする請求項３に記載のＸ線ＣＴ画像処理方法。
  6. 　Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理装置であって、
    　Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力する入力手段と、
    　前記Ｘ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、前記Ｘ線吸収係数及び前記物質に関する事後分布の期待値によって前記Ｘ線吸収係数及び前記物質を推定する推定手段と、
    　を備え、
    　前記Ｘ線投影画像の観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として設定され、
    　前記多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである、ことを特徴とするＸ線ＣＴ画像処理装置。
  7. 　Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理プログラムであって、
    　コンピュータに、
    　Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
    　前記Ｘ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、前記統計推定が前記Ｘ線吸収係数及び前記物質(及び物質のサブクラス)に関する事後分布の期待値によって前記Ｘ線吸収係数及び前記物質を推定するステップと、
    　を実行させ、
    　前記Ｘ線投影画像の観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現され、
    　前記多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである、ことを特徴とするＸ線ＣＴ画像処理プログラム。
  8. 　請求項７に記載のＸ線ＣＴ画像処理プログラムが搭載されたＸ線ＣＴ画像装置。
     
PCT/JP2014/002571 2013-05-15 2014-05-15 Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置 WO2014185078A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
EP14798390.2A EP2997900A4 (en) 2013-05-15 2014-05-15 X-ray ct image processing method, x-ray ct image processing program, and x-ray ct image device
US14/891,303 US9662081B2 (en) 2013-05-15 2014-05-15 X-ray CT image processing method, X-ray CT image processing program, and X-ray CT image device
JP2015516922A JPWO2014185078A1 (ja) 2013-05-15 2014-05-15 Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2013-102646 2013-05-15
JP2013102646 2013-05-15

Publications (1)

Publication Number Publication Date
WO2014185078A1 true WO2014185078A1 (ja) 2014-11-20

Family

ID=51898066

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/002571 WO2014185078A1 (ja) 2013-05-15 2014-05-15 Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置

Country Status (4)

Country Link
US (1) US9662081B2 (ja)
EP (1) EP2997900A4 (ja)
JP (1) JPWO2014185078A1 (ja)
WO (1) WO2014185078A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019230741A1 (ja) * 2018-05-28 2019-12-05 国立研究開発法人理化学研究所 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム
JP2020192269A (ja) * 2019-05-30 2020-12-03 キヤノン株式会社 画像処理装置、画像処理方法及びプログラム
US11307153B2 (en) 2018-05-28 2022-04-19 Riken Method and device for acquiring tomographic image data by oversampling, and control program

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6321405B2 (ja) * 2014-02-27 2018-05-09 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像生成装置、放射線断層撮影装置及び画像生成方法並びにプログラム
US9964499B2 (en) * 2014-11-04 2018-05-08 Toshiba Medical Systems Corporation Method of, and apparatus for, material classification in multi-energy image data
US10080539B2 (en) * 2015-11-25 2018-09-25 General Electric Company Systems and methods for X-ray image pasting
JP6780948B2 (ja) * 2016-03-28 2020-11-04 キヤノンメディカルシステムズ株式会社 Ｘ線ｃｔ装置
US9770221B1 (en) * 2016-07-01 2017-09-26 Siemens Medical Solutions Usa, Inc. Imaging using multiple energy levels
US10559100B2 (en) 2017-05-22 2020-02-11 Prismatic Sensors Ab Method and devices for image reconstruction
WO2019019199A1 (en) 2017-07-28 2019-01-31 Shenzhen United Imaging Healthcare Co., Ltd. SYSTEM AND METHOD FOR IMAGE CONVERSION
CN107886068B (zh) * 2017-11-10 2021-12-07 哈尔滨工程大学 一种基于时空域线谱增强器的未知线谱目标检测方法
GB201909987D0 (en) 2019-07-11 2019-08-28 Univ Surrey Apparatus and methods for non-reciprocal broken ray tomography
CN112834535B (zh) * 2020-12-30 2022-05-31 山西大学 镀膜样品膜层形态三维可视化定量ct检测方法
WO2022220721A1 (en) * 2021-04-13 2022-10-20 Prismatic Sensors Ab Determining a confidence indication for deep-learning image reconstruction in computed tomography
CN114129912B (zh) * 2021-11-25 2023-01-17 北京理工大学 一种放疗剂量验证的方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009118887A (ja) * 2007-11-12 2009-06-04 Ge Medical Systems Global Technology Co Llc Ｘ線ｃｔ装置及び画像処理装置
JP2011156302A (ja) 2010-02-03 2011-08-18 Kyoto Univ Ｘ線ｃｔ画像処理方法，ｘ線ｃｔプログラムおよび該プログラムが搭載されたｘ線ｃｔ装置
JP2013005840A (ja) * 2011-06-22 2013-01-10 Univ Of Tokyo 画像再構成装置及びプログラム

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002067201A1 (en) * 2001-02-15 2002-08-29 The Regents Of The University Of Michigan Statistically reconstructing an x-ray computed tomography image with beam hardening corrector
US9113839B2 (en) * 2003-04-25 2015-08-25 Rapiscon Systems, Inc. X-ray inspection system and method
DE102005049586A1 (de) * 2005-10-17 2007-04-26 Siemens Ag Verfahren zur Erzeugung von CT-Darstellungen in der Röntgen-Computertomographie
CN101076218B (zh) * 2006-05-19 2011-05-11 清华大学 产生具有不同能量的x射线的设备、方法及材料识别系统
DE102006026945B4 (de) * 2006-06-09 2008-04-10 Siemens Ag Computertomographisches Bildaufnahmeverfahren, Verfahren zur Bestimmung der ortsabhängigen Konzentration einer Anzahl vorgegebener Stoffe in einem Untersuchungsobjekt und zugehöriger Computertomograph
DE102010038544A1 (de) * 2009-10-19 2011-04-21 Robert Bosch Gmbh Sensorvorrichtung für eine Verpackungsmaschine
FR2953603A1 (fr) * 2009-12-09 2011-06-10 Commissariat Energie Atomique Procede et dispositif de reconnaissance d'un materiau a l'aide de sa fonction de transmission
FR2961904B1 (fr) * 2010-06-29 2012-08-17 Commissariat Energie Atomique Procede d'identification de materiaux a partir de radiographies x multi energies
US8693626B2 (en) * 2011-06-17 2014-04-08 Uop Llc Solid material characterization with X-ray spectra in both transmission and fluoresence modes
WO2013094186A1 (ja) 2011-12-18 2013-06-27 国立大学法人京都大学 動き追従ｘ線ｃｔ画像処理方法および動き追従ｘ線ｃｔ画像処理装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009118887A (ja) * 2007-11-12 2009-06-04 Ge Medical Systems Global Technology Co Llc Ｘ線ｃｔ装置及び画像処理装置
JP2011156302A (ja) 2010-02-03 2011-08-18 Kyoto Univ Ｘ線ｃｔ画像処理方法，ｘ線ｃｔプログラムおよび該プログラムが搭載されたｘ線ｃｔ装置
JP2013005840A (ja) * 2011-06-22 2013-01-10 Univ Of Tokyo 画像再構成装置及びプログラム

Non-Patent Citations (10)

* Cited by examiner, † Cited by third party
Title
A J COLEMAN; M SINCLAIR: "A beam-hardening correction using dual-energy computed tomography", PHYS. MED. BIOL., 1 985, vol. 30, no. 11, 1985, pages 1251 - 1256, XP020022878, DOI: doi:10.1088/0031-9155/30/11/007
CHYE HWANG YAN; ROBERT T. WHALEN; GARY S. BE AUPRE; SHIN Y. YEN; SANDY NAPEL: "Reconstruction Algorithm for Poly chromatic CT Imaging: Application to Beam Hardening Correction", IEEE T RANSACTIONS ON MEDICAL IMAGING, vol. 19, no. 1, 2000
IDRIS A. ELBAKRI; JEFFREY A. FESSLER: "Sen ior Member, IEEE, Statistical Image Reconstruction for Polyenergetic X -Ray Computed Tomography", IEEE TRANSACTIONS ON MEDICAL IMAGING, vol. 2 1, no. 2, 2002, pages 89 - 99
JEFFREY A. FESSLER; IDRIS ELBAKRI; PREDRAG S UKOVIC; NEAL H. CLINTHORNE: "oceedings of SPIE", vol. 4684, 2002, article "Maximum-likelihood dualenergy tomographic image reconstruction, Medical Imaging 2002: Image Processing", pages: 38 - 49
JOONKI NOH; JEFFREY A. FESSLER; PAUL E. KINAH AN: "Statistical Sinogram Restoration in Dual-Energy CT for PET Attenua tion Correction", IEEE TRANS MED IMAGING, 2009, pages 1688 - 1702
JOSEPH A. 0 'SULLIVAN; JASENKA BENAC: "Al ternating Minimization Algorithms for Transmission Tomography", IEEE T RANSACTIONS ON MEDICAL IMAGING, vol. 26, no. 3, 2007, pages 283 - 297
JOSEPH A. O'SULLIVAN; BRUCE R. WHITING; DONA LD L. SNYDER; ORVILLE A. EARL, IMAGE RECONSTRUCTION FROM DATA ACQU IRED WITH AN X-RAY COMPUTERIZED TOMOGRAPHIC SYSTEM HAVING ENERGY-INTEG RATING DETECTORS
RUHRNSCHOPF, E. P.; G. SCHWIERZ ET AL.: "Mathematical Aspects o f Computerized Tomography", vol. 8, 1981, SPRINGER, article "Nonlinearity and Inhomogenity Effects in Computerized Tomography D ue to the Exponential Attenuation of Radiation", pages: 252 - 269
See also references of EP2997900A4 *
Y.YAMAZAKI; N.TODA: "Dual-Energy X-ray CT noi se reduction algorithm under low exposure environment", IEICE TECHNICAL REPORT MBE2008-113, 2009, pages 101 - 106

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019230741A1 (ja) * 2018-05-28 2019-12-05 国立研究開発法人理化学研究所 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム
JPWO2019230741A1 (ja) * 2018-05-28 2021-07-29 国立研究開発法人理化学研究所 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム
US11307153B2 (en) 2018-05-28 2022-04-19 Riken Method and device for acquiring tomographic image data by oversampling, and control program
US11375964B2 (en) 2018-05-28 2022-07-05 Riken Acquisition method, acquisition device, and control program for tomographic image data by means of angular offset
JP7273272B2 (ja) 2018-05-28 2023-05-15 国立研究開発法人理化学研究所 角度オフセットによる断層画像データの取得方法、取得装置、および制御プログラム
JP2020192269A (ja) * 2019-05-30 2020-12-03 キヤノン株式会社 画像処理装置、画像処理方法及びプログラム
WO2020241030A1 (ja) * 2019-05-30 2020-12-03 キヤノン株式会社 画像処理装置、画像処理方法及びプログラム
JP7352382B2 (ja) 2019-05-30 2023-09-28 キヤノン株式会社 画像処理装置、画像処理方法及びプログラム

Also Published As

Publication number Publication date
EP2997900A4 (en) 2017-01-04
US9662081B2 (en) 2017-05-30
US20160120493A1 (en) 2016-05-05
EP2997900A1 (en) 2016-03-23
JPWO2014185078A1 (ja) 2017-02-23

Similar Documents

Publication Publication Date Title
WO2014185078A1 (ja) Ｘ線ｃｔ画像処理方法，ｘ線ｃｔ画像処理プログラム及びｘ線ｃｔ画像装置
JP7073134B2 (ja) 医用画像生成装置
US11769277B2 (en) Deep learning based scatter correction
Jin et al. A model-based image reconstruction algorithm with simultaneous beam hardening correction for X-ray CT
US9959640B2 (en) Iterative image reconstruction with a sharpness driven regularization parameter
US20200311490A1 (en) Apparatus and method for sinogram restoration in computed tomography (ct) using adaptive filtering with deep learning (dl)
EP2732431B1 (en) Image processing for spectral ct
US20180330233A1 (en) Machine learning based scatter correction
Lauzier et al. Characterization of statistical prior image constrained compressed sensing (PICCS): II. Application to dose reduction
CN103619259B (zh) 图像重构装置及图像重构方法
US20130202079A1 (en) System and Method for Controlling Radiation Dose for Radiological Applications
JP2014509037A (ja) 画像データのモデルに基づく処理
Wu et al. Using uncertainty in deep learning reconstruction for cone-beam CT of the brain
EP3745959B1 (en) Non-spectral computed tomography (ct) scanner configured to generate spectral volumetric image data
Zhao et al. Robust beam hardening artifacts reduction for computed tomography using spectrum modeling
Gomi et al. Evaluation of digital tomosynthesis reconstruction algorithms used to reduce metal artifacts for arthroplasty: A phantom study
Chen et al. Line integral alternating minimization algorithm for dual-energy x-ray CT image reconstruction
US10363002B2 (en) X-ray CT image reconstruction device, X-ray CT image reconstruction method, and X-ray CT apparatus
Mou et al. A beam hardening correction method based on HL consistency
Hsieh et al. An algorithm to estimate the object support in truncated images
Gomi et al. Development of a denoising convolutional neural network-based algorithm for metal artifact reduction in digital tomosynthesis for arthroplasty: A phantom study
Min et al. A weighted rebinned backprojection‐filtration algorithm from partially beam‐blocked data for a single‐scan cone‐beam CT with hybrid type scatter correction
Jacobson et al. Abbreviated on-treatment CBCT using roughness penalized mono-energization of kV-MV data and a multi-layer MV imager
Hyun et al. Deep learning for dental cone-beam computed tomography
Zhao et al. Multi-materials beam hardening artifacts correction for computed tomography (CT) based on X-ray spectrum estimation

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2015516922

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14891303

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2014798390

Country of ref document: EP