X-ray CT image processing method, X-ray CT image processing program, and X-ray CT image device

Abstract

This invention provides an X-ray CT image processing method that allows more flexible expression by expressing X-ray absorption coefficients probabilistically, makes it possible to acquire reconstructed images that are comparable to those obtained by conventional methods but involve lower X-ray doses, and can reduce beam-hardening artifacts. A probability distribution for the observation of projected X-rays is set and statistical inference is performed. Said probability distribution is expressed in terms of the process of observing a multiple-X-ray sum resulting from multiple projected X-rays being incident upon a detector. Bayesian inference in which the expected value of the posterior distribution is used for statistical inference is performed on the basis of a prior distribution for X-ray absorption coefficients, said prior distribution having parameters for the material and the observation process in terms of which the multiple-X-ray sum is expressed.

Claims

1. An X-ray CT image processing method for visualizing the inside of a scanned object in a form of a spatial distribution of X-ray absorption coefficients based on projection X-ray data obtained by detection of at least one projected X-ray at one or more X-ray detectors under several X-ray projections on an observation object, comprising: setting a probability distribution of a projection X-ray at an observation process in a case of said observation process wherein projection X-rays detected by an X-ray detector include a sum of multiple projection X-rays incident on spatially different positions, the multiple projection X-rays being incident at temporally different timings and including multiple projection X-rays of different wavelengths; setting a prior probability distribution regarding an X-ray absorption coefficient of a material inside the observation object, using a parameter concerning the material; computing a joint posterior distribution of X-ray absorption coefficient and said material's parameter, the joint posterior distribution computed from said projection X-ray data, said prior probability distribution, and said observation process; and estimating an X-ray absorption coefficient based on the said joint posterior distribution.

2. An X-ray CT image processing method as set forth in claim 1, further comprising: approximating said posterior distribution of X-ray absorption coefficient with a test distribution using the said observation process and said prior distribution so that the distributional distance between said posterior distribution and the test distribution is small; and consequently computing an expected value of the test distribution regarding said X-ray absorption coefficient.

3. An X-ray CT image processing method as set forth in claim 1, wherein a sum of projected X-rays of different X-ray wavelengths is included in a sum of said multiple projection X-rays, and a prior probability distribution regarding the wavelength dependent X-ray absorption coefficients of each pixel is set using a parameter for expressing a distribution of an X-ray absorption coefficient for each said material and each X-ray wavelength.

4. An X-ray CT image processing method as set forth in claim 1, wherein a sum of the projection X-rays of different X-ray wavelengths is included in a sum of said multiple projection X-rays, and an X-ray absorption coefficient possessing a wavelength dependence of each pixel is set as a product of an X-ray absorption coefficient density not depending on the wavelength and an X-ray absorption coefficient ratio among X-ray wavelengths specified for each said material, not depending on a pixel, the method comprises using a parameter for expressing an X-ray absorption coefficient ratio among X-ray wavelengths specified for each said material, and the prior probability distribution of X-ray absorption coefficient possessing wavelength dependence of each pixel is set by a prior probability distribution conditional to said X-ray absorption coefficient density of each pixel under a condition of said material and a prior probability distribution for said material.

5. An X-ray CT image processing method as set forth in claim 3, wherein a prior probability distribution for said material is set by a prior probability distribution expressed by a parameter for representing the content percentage of each material and a parameter for showing an extent of spatial continuity of each material.

6. An X-ray CT image processing method as set forth in claim 4, wherein a conditional prior probability distribution of said X-ray absorption coefficient density specified for each said material is set by a conditional prior probability distribution of said X-ray absorption coefficient density specified for each subclass of each material and a mixed probability distribution resulting from a sum regarding all attainable subclasses of each material for a product of a conditional prior probability distribution of a subclass of each material, and a prior probability distribution of subclass for said each material is set by a parameter for showing appearance extent of a subclass for each material and a parameter for showing the spatial continuity extent of a subclass for each said material.

7. An X-ray CT equipment for imaging an inside of an observation object as a spatial distribution of an X-ray absorption coefficient based on projection X-ray data obtained by detecting a projected X-ray on an X-ray detector after projecting an X-ray on an observation object, the X-ray CT equipment comprising: a means for setting a probability distribution of projection X-rays at an observation process in a case where a projection X-ray detected by an X-ray detector is a sum of any or all of three types of multiple X-rays: (1) multiple projection X-rays incident on spatially different positions, (2) multiple projection X-rays incident at temporally different timings, and (3) multiple projection X-rays of different wavelengths; a means for setting a prior probability distribution regarding a parameter for expressing an X-ray absorption coefficient of a material constituting an inside of an observation object and a material; and a means for estimating a spatial distribution of said X-ray absorption coefficient by estimating a posterior probability distribution regarding said X-ray absorption coefficient and said material from said projection X-ray data, a probability distribution expressing a likelihood of said projection X-ray data and said prior probability distribution.

8. An X-ray CT equipment comprising: an X-ray detection mechanism including at least one X-ray detector; a computer which executes an image processing program that performs a method to assist imaging an inside of an observation object as a spatial distribution of X-ray absorption coefficients based on projection X-ray data obtained by detecting projected X-rays at the X-ray detection mechanism by projecting an X-ray on an observation object in an observation process, the method comprising (a) setting a probability distribution of a projection X-ray wherein a projection X-ray detected by the X-ray detection mechanism is a sum of multiple X-rays consisting of at least one projection X-ray chosen from multiple projection X-rays to be detected by the X-ray detection mechanism incident on a spatially different positions to a detection plane of the X-ray detection mechanism, multiple projection X-rays incident at temporally different timings and multiple projection X-rays of different wavelengths, (b) setting a prior probability distribution regarding a parameter for expressing an X-ray absorption coefficient of a material constituting an inside of an observation object and said material, and (c) estimating a spatial distribution of said X-ray absorption coefficients by estimating a posterior probability distribution regarding said X-ray absorption coefficient and said material from said projection X-ray data from said projection X-ray data and a probability distribution of a projection X-ray at said observation process and said prior probability distribution.

9. An X-ray CT equipment as set forth in claim 8, wherein said X-ray absorption coefficient estimation comprises calculating a test distribution regarding said X-ray absorption coefficient and said material from a probability distribution of a projection X-ray at said observation process and said prior probability distribution for making said posterior distribution and inter distribution distance thereof small; and the method comprises estimating an expected value of test distribution regarding said X-ray absorption coefficient.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 A process flow of the X-lay CT image processing method of the present invention

(2) FIG. 2 An explanatory drawing of the X-ray CT

(3) FIG. 3 An explanatory drawing of the update

(4) FIG. 4 A probability distribution drawing of the metal and bone X-ray absorption coefficients

(5) FIG. 5 A figure showing the image reconstruction at the computer experiment of embodiment 1

(6) FIG. 6 A figure showing the image reconstruction at the computer experiment of embodiment 2

(7) FIG. 7 An explanatory drawing of the projection geometry for the parallel beam

(8) FIG. 8 An explanatory drawing of the projection and the encoding for the parallel beam

(9) FIG. 9 An explanatory drawing of the projection geometry under rotation

BEST MODE FOR CARRYING OUT THE INVENTION

(10) Embodiments of the present invention will be described in detail below with reference to the drawings. The present invention is not limited to the following embodiment and examples of shown in the figure, and the present invention can be variously changed in design.

(11) In the X-ray CT image processing method and the X-ray CT image processing program according to the present invention, the image reconstruction is performed using a statistical inference along with the setting of the prior distribution regarding the X-ray absorption coefficient. Especially, the beam hardening artifacts are disbanded by the estimation of the energy dependent (wavelength dependent) X-ray absorption coefficient based on a statistics model considering the energy distribution (the wavelength distribution).

(12) To address the ill-posedness induced by the increase of variables to be estimated, it is necessary to solve the problem by imposing some kind of constraint on the solution by using some proper prior knowledge. In the present invention, by using a basis function expression, the basis thereof being assumed to be an X-ray absorption coefficient dependent on energy (dependent on wavelength) of a known material and being also assumed to be given by a product of a discrete variable for expressing the material class the weight belongs to and a continuous variable of a non-negative scalar for expressing an X-ray absorption coefficient density, a constraint is imposed on the solution by imposing a prior distribution of X-ray absorption coefficients under a condition of a material class and the prior distribution of the material class.

(13) Also, according to the present invention, in the case of a reconstruction image targeting human tissue, the human tissue consists of limited materials such as fat, muscle, bone, and because a rough distribution of the X-ray absorption coefficient of each material is known in advance, such knowledge as mentioned thus far is expressed as a prior distribution to be utilized for estimation. For this estimation, it is necessary to estimate the material class each pixel belongs to, however, based on the fact that each material possesses spatial continuity and a situation in which the percentage for each material to occupy the human body can be assumed, ill-posedness can be restrained by introducing the prior knowledge such as a spatial smoothness constraint and an occupancy ratio to the discrete variable for expressing the material class.

(14) The observation process of a projected X-ray expresses the observation process of multiple-X-ray sum of the projected X-ray incident on a detector in a multiple mode, and also expresses that multiple X-rays of different wavelengths are incident on a detector, the X-rays having wavelength regions (spectral distributions) of a finite width.

(15) The probability distribution regarding the reconstruction images characterizes the way for expressing the tendency of the X-ray absorption coefficient value to be settled at, in the region of each pixel of the reconstruction images, and it is expressed by a Boltzmann distribution consisting of a mixed gamma distribution expressing the tendency of the X-ray absorption coefficient density for the combination of each material class and each subclass variable to take, a parameter expressing the extent for the combination of each material class and each subclass variable to emerge and a parameter expressing the extent for the combination of each material class and each subclass variable to be spatially continuous.

(16) And the statistical inference by the expectation estimation of the posterior distribution regarding the X-ray absorption coefficient and the material class is performed. At image reconstruction, it is necessary to find this posterior distribution. However, it is difficult to execute the calculation analytically because the sum calculation regarding the high order hidden discrete state variable which is a material class defined at each pixel is included. Therefore, in this invention, this calculation difficulty is overcome by applying an approximation method. (Refer to a flowchart in FIG. 1)

(17) (Formalization of the Problem)

(18) As in the X-ray CT image, T projection data projected from various directions are to be expressed by D={Y.sup.(1), ***, Y.sup.(T)}. Each data Y.sup.(t) is a set of data detected by the detector through the t-th projection and becomes Y.sup.(t)={y.sub.1.sup.(1), ***, y.sub.1.sup.(t)}. Note that T is a projection number, I a number of the detector and y.sub.i.sup.(t) is the photon number detected at the i-th detector. And the equation (1) below is formed, because the X-ray is attenuated exponentially when transmitted through a material.

(19) $\begin{array}{cc}\left[\mathrm{Equation}1\right]{\hat{y}}_i^{\left(t\right)}=b_i^{\left(t\right)}\exp \left\{-\underset{j=1}{\overset{J}{.Math.}}{l}_{ij}^{\left(t\right)}{x}_j\right\}& \left(1\right)\end{array}$

(20) Here, x.sub.j is an X-ray absorption coefficient of the j-th pixel in the J-dimensional vector x={X.sub.1, ***, X.sub.J} obtained by a raster scan of the X-ray absorption coefficient of the imaging object on an observation object. b.sub.i.sup.(t) expresses the photon number emitted from the X-ray source (the number of photons observable when no object is placed). Also, l.sub.ij.sup.(t) corresponds to the intersection distance between the projection line detected by the i-th detector when projected from the angle .sup.(t) and the i-th pixel, and l.sub.i j.sup.(t) x.sub.j expresses the effective X-ray absorption coefficient of the j-th picture region against the X-ray incident on the i-th detector for the t-th projection. (Refer to FIG. 2)

(21) Here, the expression for discretization of the X-ray energy spectrum into E is considered. If the photon number of the X-ray belonging to energy e (e=1, ***, E) within the incident X-ray is assumed to be N e.sup.(t), the average photon number observed after transmitting a material can be expressed by the equation (2) below.

(22) Note that x.sub.je is the X-ray absorption coefficient of the j-th pixel against the X-ray of energy e.

(23) $\begin{array}{cc}\left[\mathrm{Equation}2\right]{\hat{y}}_i^{\left(t\right)}=\underset{e=1}{\overset{E}{.Math.}}{N}_e^{\left(t\right)}\exp \left\{-\underset{j=1}{\overset{J}{.Math.}}{l}_{ij}^{\left(t\right)}{x}_{\mathrm{je}}\right\}& \left(2\right)\end{array}$

(24) (Energy Dependence of the X-Ray Absorption Coefficient)

(25) The ill-posedness is generated if the X-ray absorption coefficient x.sub.je of j-th pixel at an energy e is freely determined by each energy and each pixel. Practically, every material has a character, roughly common among materials, that the X-ray absorption coefficient becomes smaller with the increase of energy. Namely, the ratio of X-ray absorption coefficients of different energy shows a specific tendency.

(26) Therefore, the X-ray absorption coefficient is expressed under a constrained way by introducing a material class for classifying materials.

(27) Assume that the material for observation is classified into C kinds of material classes. Here, the material class of the j-th pixel is expressed by using a variable z.sub.j. Note that Z is expressed as Z.sub.j={Z.sub.j1,***, Z.sub.jC}. Each element z.sub.jC is a binary variable taking either 0 or 1, taking 1 when the j-th pixel belongs to the material class C and taking 0 for other cases. As a variable for expressing a ratio of X-ray absorption coefficient among energies uniquely determined depending on this material class, a rate of variability r.sub.cs (>0) is defined. Note that the rate of variability r.sub.cs is assumed to be normalized as r.sub.cs=1 when e=1. The X-ray absorption coefficient of the j-th pixel under the condition e=1, is expressed by X.sub.j(0>). By using the expression above, the absorption coefficient x.sub.js can be expressed as shown in the equation (3) below.

(28) $\begin{array}{cc}\left[\mathrm{Equation}3\right]& \\ x_{\mathrm{je}}=\underset{c=1}{\overset{C}{.Math.}}{r}_{\mathrm{ce}}{z}_{\mathrm{jc}}{x}_j& \left(3\right)\end{array}$

(29) When the equation (3) above is used, the variables to be estimated at each pixel become only two, the one is an x.sub.j scalar continuous variable and another is a discrete variable z taking one of the values of C kinds.

(30) (Regarding the Observation Model)

(31) The ideal value of the photons observed at each detector follows the equation (1) mentioned above. However practically, it is known that the observation value is off the ideal value due to noises such as photon scattering and quantization noise. If these noises are assumed to be well expressed by the Gaussian distribution, the observation data at each detector can be expressed by the equation (4) below. Note that X={x1,***, X.sub.j}, and Z={z, ***, z.sub.j}.

(32) Also, if the observation is assumed to be independent for each detector and for each projection, a set D of the projection data (observation data) can be expressed by the equation (5) below. Note that T is the number of projection and I is the number of detector. Here, that represents a standard deviation of the Gaussian distribution can be treated as a function to be dependent on the observation y.sub.i.sup.(t) as shown by the equation (4) below, considering the nature of the shot noise.

(33) $\begin{array}{cc}\left[\mathrm{Equation}4\right]& \\ p\left(y_i^{\left(t\right)}.Math.X,Z\right)=\frac{1}{\sqrt{2}}\exp \left\{-\frac{{\left({\hat{y}}_i^{\left(t\right)}-y_i^{\left(t\right)}\right)}^2}{2{}^2}\right\}& \left(4\right)\end{array}$
={square root over (y.sub.i.sup.(t)+v)}

(34) (Note v is a constant)

(35) $\begin{array}{cc}\left[\mathrm{Equation}5\right]& \\ \begin{array}{c}p\left(D.Math.X,Z\right)=\underset{t=1}{\overset{T}{.Math.}}\underset{i=1}{\overset{I}{.Math.}}p\left(y_i^{\left(t\right)}.Math.X,Z\right)\\ =\underset{t=1}{\overset{T}{.Math.}}\underset{i=1}{\overset{I}{.Math.}}\frac{1}{\sqrt{2}}\exp \left\{-\frac{{\left({\hat{y}}_i^{\left(t\right)}-y_i^{\left(t\right)}\right)}^2}{2{}^2}\right\}\end{array}& \left(5\right)\end{array}$

(36) (Regarding the Prior Distribution)

(37) With regard to the prior distribution, the prior knowledge regarding the X-ray absorption coefficient for each material class of the observation object is expressed. If the X-ray absorption coefficient x.sub.j at the pixel j is assumed to be determined depending only on the material class z.sub.jc at the pixel j, the prior distribution p (X/Z) can be expressed by the equation (6) below.

(38) $\begin{array}{cc}\left[\mathrm{Equation}6\right]& \\ p\left(X.Math.Z\right)=\underset{j=1}{\overset{J}{.Math.}}p\left(x_j.Math.z_j\right)& \left(6\right)\end{array}$

(39) When the material class Z.sub.jc is given, the expectation for the X-ray absorption coefficient x.sub.j tends to take some specific value. Here, the prior distribution p (s.sub.j/z.sub.j) can be expressed according to the equation (7) below as a mixture distribution of 2 distributions so that the asymmetry of the distribution and various different kurtosis can be expressed. Note that the variable b is a binary variable. Also the vector B is defined as B={b.sub.1, ***, b.sub.J}.

(40) $\begin{array}{cc}\left[\mathrm{Equation}7\right]& \\ p\left(x_j.Math.z_j\right)=\underset{b_j}{.Math.}p\left(x_j.Math.z_j,b_j\right)p\left(b_j\right)& \left(7\right)\end{array}$
b.sub.j=[b.sub.j1,b.sub.j2]{[0,1],[1,0]}

(41) When the material class z.sub.jc and the variable b.sub.jc for each pixel are given, the X-ray absorption coefficient x.sub.j can be expressed, under an assumption that it follows the gamma distribution, by the equation (8) below. Note that u.sub.cd, and v.sub.cd are parameters of gamma distribution, u.sub.cd/v.sub.cd being an average and u.sub.cd/v.sup.2.sub.cd being dispersion. These parameters are set by empirical knowledge regarding the material class.

(42) $\begin{array}{cc}\left[\mathrm{Equation}8\right]& \\ p\left(x_j.Math.z_j,b_j\right)=\underset{c=1}{\overset{C}{.Math.}}\underset{d=1}{\overset{2}{.Math.}}{\left[\frac{1}{\left(u_{cd}\right)}{v}_{cd}^{u_{cd}}{x}_j^{u_{cd}-1}{}^{-v_{cd}{x}_j}\right]}^{z_{\mathrm{jc}}{b}_{\mathrm{jd}}}& \left(8\right)\end{array}$

(43) Also, in the prior distribution, the variable Z.sub.jc satisfying the equation (9) below is a binary variable taking 0 or 1. Employment of the gamma distribution can express that the X-ray absorption coefficient does not become a negative value.
[Equation 9]
.sub.cz.sub.jc=1(9)

(44) Here, regarding the prior distribution of the variables Z and B, independency among variables is not assumed and follow the gamma distribution as is expressed by the following equation (10). Note that A is a normalization term and the energy H (Z, B) is defined by the following equation (11). In the formula (11) below, 4 pixels neighboring lengthwise and breadthwise are considered to be the vicinity of a 2 dimensional plane and each G.sub.cd.sup.self and G.sub.cd.sup.inter is a non-negative constant.

(45) $\begin{array}{cc}\left[\mathrm{Equation}10\right]& \\ p\left(Z,B\right)=\frac{1}{A}\exp \left\{-H\left(Z,B\right)\right\}& \left(10\right)\\ \left[\mathrm{Equation}10\right]& \\ p\left(Z,B\right)=\frac{1}{A}\exp \left\{-H\left(Z,B\right)\right\}& \left(10\right)\\ \left[\mathrm{Equation}11\right]& \\ H\left(Z,B\right)=-\underset{j=1}{\overset{J}{.Math.}}\underset{c=1}{\overset{C}{.Math.}}\underset{d=1}{\overset{2}{.Math.}}{G}_{cd}^{\mathrm{self}}{z}_{jc}{b}_{\mathrm{jd}}-\underset{j=1}{\overset{J}{.Math.}}\underset{s\left(j\right)}{.Math.}\underset{c=1}{\overset{C}{.Math.}}\underset{d=1}{\overset{2}{.Math.}}\underset{k=1}{\overset{2}{.Math.}}{G}_{cd}^{\mathrm{inter}}{z}_{\mathrm{jc}}{z}_{\mathrm{sc}}{b}_{\mathrm{jd}}{b}_{\mathrm{sk}}& \left(11\right)\end{array}$

(46) (j): A set of neighboring pixels for the j'th pixel

(47) (About Posterior Distribution)

(48) The reconstruction of the X-ray CT image according to the present invention, an X-ray absorption coefficient X, a material class Z and a variable B are estimated as a posterior distribution. The X-ray absorption coefficient X, the posterior distribution p (X, Z, B/D) can be expressed by the equation (12) below by the Bayesian theorem.
[Equation 12]
p(X,Z,B|D)p(D|X,Z)p(X|Z,B)p(Z,B)(12)

(49) (Bayesian Inference)

(50) Although the expected value of the posterior distribution p (X, Z, B/D) expressed by the equation (12) above becomes necessary according to the Bayesian inference, the posterior p (X, Z, B/D) is approximated by the test distribution q (X, Z, B) in view of the difficulty of estimating this posterior distribution p (X, Z, B/D) analytically. The test distribution q (X, Z, B) can be selected arbitrarily as long as it can minimize or approximately minimize the equation (12) mentioned above. And the tolerance of the test distribution q (X, Z, B) and the posterior distribution p (X, Z, B/D) is evaluated by the KL (Kullback-Leibler) distance to compute the test distribution q (X, Z, B) that minimize the KL distance. Here, the KL distance can be expressed by the equation (13) below.

(51) Note that <*>.sub.q(X, Z, B) is an operator for performing integral computation regarding the distribution q (X, Z, B). The KL distance is always non-negative and becomes 0 (zero) only when q=p. <*>.sub.q(X, Z, B) expresses that it takes the expected value regarding q (X, Z, B).

(52) 0$\begin{array}{cc}\left[\mathrm{Equation}13\right]& \\ D_{\mathrm{KL}}={.Math.\log \frac{p\left(X,Z,B.Math.D\right)}{q\left(X,Z,B\right)}.Math.}_{q\left(X,Z,B\right)}& \left(13\right)\end{array}$

(53) To make this optimization easy, it is assumed that the test distribution q (X, Z, B) is given by the product of its marginal distribution q (X) and q (Z, B), namely there is dependence between the variable X and {Z, B}. Also, each element of q (X) and q (Z, B) are assumed to be independent among pixels, and further each element q (x.sub.j) is assumed to obey the gamma distribution.

(54) Under the conditions described above, q (Z.sub.j, b.sub.j) consequently becomes the multinomial distribution because the probability variables Z.sub.j and b.sub.j are discrete variables. From above, the test distribution q (X, Z, B) can be expressed by equation that are the equation (14) and the equation (15) below. Note that the alternate optimization as shown by the equation (16) below is performed because it is difficult to optimize test distributions a(x) and q (Z, B) simultaneously.

(55) $\begin{array}{cc}\left[\mathrm{Equation}14\right]& \\ \begin{array}{c}q\left(X,Z,B\right)=q\left(X\right)q\left(Z,B\right)\\ =\underset{j=1}{\overset{J}{.Math.}}q\left(x_j\right)q\left(z_j,b_j\right)\end{array}& \left(14\right)\\ \left[\mathrm{Equation}15\right]& \\ q\left(X\right)=\underset{j}{.Math.}\left[\frac{1}{\left({}_j\right)}{}_j^{{}_j}{x}_j^{{}_j-1}{}^{-{}_j}\right]& \left(15\right)\\ \left[\mathrm{Equation}16\right]& \\ q\left(X\right)=\arg \underset{q\left(X\right)}{\min }{D}_{\mathrm{KL}}& \left(16\right)\\ \left[\mathrm{Equation}17\right]& \\ q\left(z_j,b_j\right)=\arg \underset{q\left(z_j,b_j\right)}{\min }{D}_{\mathrm{KL}}& \left(17\right)\end{array}$

(56) Here, the SCG (Scaled Conjugate Gradient) method was employed for the optimization of the parameters q(x.sub.i) .sub.j, .sub.j in the q(x.sub.j). Under the fixation of the test distribution q (X) and q (z.sub.j, b.sub.j) (note that j differs from j) other than the j-th pixel, q (z.sub.j, b.sub.j) that minimizes the KL distance can be analytically found. This analytical solution depends on the neighboring pixel q(z.sub.j, b.sub.j) (j.sub.(j)) of j. Thus, in order to renovate the q (Z, B) without changing the distribution of neighboring pixels, the subsets of pixels to become checkered as shown in FIG. 3 are alternately optimized. Namely, the parameter q (Z, B) does not renovate all the pixels simultaneously but optimize the subset of pixels in the checkered patterns shown in the FIG. 3. After renovating q (zj, bj) corresponding to the gray portion in FIG. 3, q (z.sub.j, b.sub.j) corresponding to pixels in the gray portion is renovated.

(57) Thus far, the X-ray CT image processing method and the X-ray CT image processing program according to the present invention have been explained. In the embodiments below, a computer experiment to study the effectiveness of the present invention is conducted and the screen reconstruction of the X-ray image processing method and the X-ray image processing program are evaluated by comparing these with the screen reconstruction by the conventional technology.

(58) By the evaluation experiment described below, the usefulness of the X-ray CT image processing method and the X-ray CT image processing program according to the present invention will be understood.

(59) In the embodiment, in order to be able to measure the error between the reconstruction image and the genuine image, a proper genuine image was prepared, and against that proper genuine image, the X-ray projection simulation was performed to artificially generate a projection image (sinogram). For the X-ray projection, the energy of the incident X-ray is divided into 5 values, namely 60 keV, 70 keV, 80 keV, 90 keV and 100 keV (the photon numbers for each are 610.sup.4, 2.510.sup.4, 1.510.sup.3, 710.sup.3 and 510.sup.3), in order to consider the energy dependence (wavelength dependence) of the material X-ray absorption coefficient and each divided X-ray was attenuated according to the X-ray absorption coefficient specific for each material. It was assumed that the photon of the attenuated X-ray with addition of a noise (noise) generated according to the Poisson noise was observed at the detector.

(60) In the embodiments 1 and 2 below, the followings were assumed, namely, the number of the detector was 95 (pieces), the resolution of the image was 6464 (pixels) and the projection interval was 1 between 1 and 180. In order to quantitatively evaluate the inference quality of the reconstruction image, PSNR (Peak Signal to Noise Ratio) expressed by the equation (18) below was used.

(61) $\begin{array}{cc}\left[\mathrm{Equation}18\right]& \\ \mathrm{PSNR}=10{\log }_{10}\left(\frac{{\mathrm{MAX}}^2}{\mathrm{MSE}}\right)\left(\mathrm{dB}\right)& \left(18\right)\end{array}$

(62) In the equation (18) above, MAX is the maximum value of the pixel value and MSE is an average squared error between a true image and a reconstruction image. MSE is smaller with higher PSNR. Namely, the image is nearer to the true image. What is described below is a comparison result of image reconstruction by the following three methods that are, the FBP method, the method to reconstitute the image by the Bayesian inference assuming the prior distribution not considering X-ray energy distribution to follow the Gaussian distribution (the Conventional method; hereafter) and the X-ray CT image processing method by discretizing the X-ray energy into 3 energies, according to the present invention (The present invention method hereafter), followed by comparison results of each other.

(63) Note that the prior distribution of metal and bone is shown in FIG. 4. FIG. 4 shows the prior distribution of metal and bone and the X-ray absorption coefficient of the true image at 60 keV. The X-ray absorption coefficients of metal and bone are probability distribution expressed by the gamma distribution with 0.7 and 0.4 as average values. The X-ray absorption coefficient of the genuine image is set at a value 5% smaller than the average of the prior distribution.

Embodiment 1

(64) For embodiment 1, a computer experiment regarding the image reconstruction against a phantom including 4 metals is explained referencing FIG. 5.

(65) FIG. 5 (a) is a genuine image, (b) is a reconstruction image by the FBP method, (c) is a reconstruction image by the simple Bayesian method and (d) is the reconstruction image according to the present invention.

(66) In embodiment 1, image reconstruction experiments were conducted on the genuine image as shown in FIG. 5 (a), employing each method, namely the FBP method, the simple Bayesian method and the method in this invention. The sum of the photon number for the incident X-ray was assumed to be 1.1210.sup.5, and in the X-ray CT image processing method of the present invention, the X-ray energy was discretized to three intensities, that were 60 keV, 80 keV and 100 keV and the photon numbers belonging to each were set at 7.510.sup.4, 310.sup.4 and 10 0.710.sup.4, respectively. The repetition number was set both at 25, for comparison. The inference result is shown in FIG. 5 (b) to (d).

(67) Also, PSNR by each FBP method, the simple Bayesian method and the method of the present invention is shown in Table 1 below. It is understood from the table that the PSNR for the FBP method and the simple Bayesian method are larger than the PSNR for the method of the present invention and the present invention method can conduct higher precision inference.

(68) TABLE-US-00001 TABLE 1 Simple Present Bayesian invention FBP method method method PSNR (dB) 16.7489 22.3419 27.3353

Embodiment 2

(69) For embodiment 2, a computer experiment regarding the image reconstruction against a phantom of teeth including a metal implant is explained referencing FIG. 6.

(70) FIG. 6 (a) is an genuine image, (b) is a reconstruction image by the FBP method, (c) is a reconstruction image by the simple Bayesian method and (d) is the reconstruction image according to the present invention.

(71) In embodiment 2, image reconstruction experiments were conducted on the genuine image as shown in FIG. 5 (a), employing each method, namely the FBP method, the simple Bayesian method and the method in this invention. The sum of the photon number for the incident X-ray was assumed to be 1.1210.sup.5, and in the X-ray CT image processing method of the present invention, the X-ray energy was discretized to three intensities, that were 60 keV, 80 keV and 100 keV and the photon numbers belonging to each were set at 7.510.sup.4, 310.sup.4 and 10 0.710.sup.4, respectively. The repetition number was set both at 25, for comparison. The inference result is shown in FIG. 5 (b) to (d).

(72) Also, PSNR at each FBP method such as, simple Bayesian method and the method of the present invention is shown in Table 1 below. It is understood from the table that the PSNR for FBP method and simple Bayesian method is larger than the PSNR for the method of the present invention and the method of present invention can conduct higher precision inference.

(73) TABLE-US-00002 TABLE 2 Simple Present Bayesian invention FBP method method method PSNR (dB) 16.6327 21.5578 25.1628

(74) As was explained above, it was shown that various artifacts originating to the influence of the beam hardening can be reduced by considering the energy distribution of the X-ray spectrum and estimating X-ray absorption coefficient with energy dependence (wavelength dependence) based on the statistic model. As the computer experiments according to embodiments 1 and 2, only those cases where the X-ray energy was divided into three energies were selected. However, it can be seen that the artifacts were extensively reduced even when the number of the discrete representation is smaller than the partition number of the accrual energy. Furthermore, it would be possible to conduct a better estimation by further increasing the partition number of the energy.

Embodiment 3

(75) Next, regarding the observation process of the projected X-ray observation, it will be explained that the multiple-X-ray sum can express the phenomena in which multiple X-rays of different spatial positions, directions and wavelengths are incident on the detector under the conditions that the detection surface of the detector possesses a spatially finite area, the detection time of the X-ray at the detector has a finite width and the X-ray possesses a wavelength region (a spectral distribution) having a finite expanse.

(76) Firstly, an explanation is given by referencing FIG. 7 with a parallel beam as an example. A projection by a parallel beam in a 2 dimensional plane as shown in FIG. 7 is assumed. Here, the direction of the X-ray radiation is defined as the y-axis and the direction orthogonal thereto is defined as the u-axis. Variables are defined as below in this coordinate system.

(77) The X-ray absorption coefficient of a material at a position (u, v) for a wavelength is to be expressed by x (u, v, ). Note that the material is assumed to be included in vL.

(78) Let the photon number incident on the detection surface S.sub.i of the i-th detector be y.sub.i (). Let the photon number of a wavelength radiated per unit time and unit area by the X-ray source of a parallel beam be d () and let an average number of electrons excited by one photon of a wavelength be ().

(79) Under this condition, the expected value of the current Y.sub.i observed by the i-th detector per unit time is expressed by the equation (19) below. Furthermore, in a case where the current accumulated for a specific time duration is observed, the time integral of the current is to be performed.
[Equation 19]
Y.sub.i=()y.sub.i()d
y.sub.i()=.sub.uS.sub.id()exp(.sub.vLx(u,v,)dv)du(19)

(80) It is understood that artifacts are generated if the estimation method that assumes a line integral is employed in a case where a spatial expanse of the detection surface and an expanse of X-ray spectral distribution exist because the formula (19) above cannot arrive at a line integral (an integral regarding v described above) assumed by a reverse projection method according to the formula (19) above. Similarly, in a case where rotation is considered at a time when the detector receives a light, the formula (19) prevents the arrival at the line integral. Approximation of integral by sum is considered because the calculation by the multiple integral mentioned above is difficult. It is shown that the above leads to a consideration on a multiple-X-ray sum. Firstly, the integral regarding u, v in the formula (20) below is replaced to a sum, noting only on a specific wavelength in the formula (19) above.
[Equation 20]
y.sub.i=.sub.uS.sub.idexp(.sub.vLx(u,v)dv)du(20)
[Equation 20]
y.sub.i=.sub.uS.sub.idexp(.sub.vLx(u,v)dv)du(20)

(81) In order to approximately express the integral regarding u and v above on computer, a discretization is performed. (Refer to FIG. 8) When the space is discretized along u axis and v axis, the integral regarding v in the equation (20) above can be approximated following the equation (21) below.
[Equation 21]
.sub.vLx(u,v)dv.sub.jI(u)x.sub.jL(21)

(82) Here, X.sub.j is the X-ray absorption coefficient of the i-th pixel, I(u) expresses the set of pixel index placed at the position u on the u axis and L expresses the length of one side of the pixel. When the length of one side of the pixel is expressed by L=1, the equation (21) above is expressed by the equation (22) below.
[Equation 22]
.sub.vLx(u,v)dv.sub.jI(u)x.sub.j(22)

(83) Accordingly, when the length of one side of the pixel is 1, y.sub.i is expressed by the equation (23) below and further, when an integral regarding a ray (an integral regarding u) incident on the same detector is approximated by the sum, y.sub.i is expressed by the equation (24) below.

(84) $\begin{array}{cc}\left[\mathrm{Equation}23\right]& \\ y_i={}_{n{S}_i}\exp \left(-{}_{vL}x\left(u,v\right)v\right)u{}_{u{S}_i}\exp \left(-\underset{jI\left(u\right)}{.Math.}{x}_j\right)u& \left(23\right)\\ \left[\mathrm{Equation}24\right]& \\ \underset{\left\{k.Math.u_0+\frac{k}{N}{S}_i\right\}}{.Math.}\frac{1}{N}d\exp \left(-\underset{jI\left(u_0+\frac{k}{N}\right)}{.Math.}{x}_j\right)& \left(24\right)\end{array}$ 1/N is the rectangle width for the approximation of integral by thin rectangles.)

(85) In the equation (24) above, 1/N represents the width of a rectangle when the integral is approximated by cutting out rectangles.

(86) This corresponds to thinking of the sum of projection expressed by a linear integral about the index k.

(87) As has been explained thus far, it is known that the multiple X-ray needs to be considered when the spatial expanse of the detector is taken into account.

(88) The above assumes a case where the X-ray is incident on the pixel directly. However, for an equipment in which the X-ray source and the detector rotate facing each other (refer to FIG. 9) such as in an ordinary X-ray CT equipment, it is necessary to anticipate a situation where the X-ray is radiated at the pixel obliquely.

(89) In this case, the equation (24) above about the projection can be expressed by generalizing the equation (24) unto the following equation (25) below.

(90) $\begin{array}{cc}\left[\mathrm{Equation}25\right]& \\ y_i\underset{\left\{k.Math.u_0+\frac{k}{N}{S}_i\right\}}{.Math.}\frac{1}{N}d\exp \left(-\underset{jI\left(u_0+\frac{k}{N}\right)}{.Math.}{x}_j\right)& \left(25\right)\\ \left[\mathrm{Equation}26\right]& \\ y_i\underset{\left\{k.Math.r\left(k\right)S_i\right\}}{.Math.}d\exp \left(-\underset{j{I}_k}{.Math.}{l}_{\mathrm{kj}}{x}_j\right)N& \left(26\right)\end{array}$

(91) Here, N=1/N, r(k) expresses the position of a beam k on the detection surface and {k|r(k)S.sub.i} expresses the aggregation of beams reaching the detector j. I.sub.k expresses the index aggregation of the pixel the beam k intersects and l.sub.kj expresses the intersection length of the beam k and the pixel j. The equation (25) above can be considered to be the case that l.sub.kj in the equation (26) above is assumed to be 1.

(92) Regarding the equation (26) above, a similar equation holds by computing l.sub.kj considering the beam direction even for a fan beam or a corn beam besides a parallel beam.

(93) What the above formula (26) suggests is the necessity to consider the projection by each beam that transmits the detection surface when an observation model of the projection considering the detection surface expansion is taken into account.

(94) When the sum of l.sub.kjx.sub.j in the equation (26) above is written altogether regarding the detector i and the beam k, it can be expressed by a matrix operation such as Lx. The case in which the matrix L is (the number of the detector)(the pixel number of the reconstruction image) when the detection surface expansion is not considered changes to another case in which the matrix L becomes (the number of beams that transmits the detector) x (the pixel number of the reconstruction image) when the detection surface expansion is considered. While the greater the beam number is, the higher the integral precision becomes, a problem arises that the calculation volume increases.

(95) Next, a rotation during the projection (an integral about time) is assumed. First of all, if the photon number observed at a time t is assumed to be y.sub.i,t, it can be expressed by the equation (27) below.

(96) $\begin{array}{cc}\left[\mathrm{Equation}27\right]& \\ y_{i,t}=\underset{\left\{k.Math.r\left(k\right)S_i\right\}}{.Math.}d\exp \left(-\underset{j{I}_k}{.Math.}{l}_{\mathrm{kj}}^{\left(t\right)}{x}_j\right)N& \left(27\right)\end{array}$

(97) The photon number detected during the detection time (the exposure time) of the detector becomes the integral of y.sub.i,t about time and accordingly can be expressed by the equation (28) below.
[Equation 28]
.sub.tTy.sub.i,tdt(28)

(98) Regarding the equation (28) above, when the integral about time is replaced by the sum, the expression by the equation (29) below becomes available.
[Equation 29]
.sub.{k|t(k)T}y.sub.i,t(k)t(29)

(99) The photon y.sub.i,t incident on the detector at each time t is expressed by the sum of bundles of multiple rays as is shown in the equation (28) above when the detection surface expansion is considered. The integral about time is equivalent to thinking of a sum regarding a new bundle of rays because a further sum is included. Therefore, the photon number measured during the detection time T (the exposure time) can be expressed by the equation (30) below when the set S of ray indices is expressed by S.sub.j,t which is time dependent and the set I.sub.k of pixel indices each ray k transmits through is redefined.

(100) $\begin{array}{cc}\left[\mathrm{Equation}30\right]& \\ \underset{\left\{k.Math.r\left(k\right)S_{i,T}\right\}}{.Math.}d\exp \left(-\underset{j{I}_k}{.Math.}{l}_{\mathrm{kj}}{x}_j\right)N& \left(30\right)\end{array}$

(101) The expression by the equation (31) below holds when the integral about the wavelength in the equation (30) above is replaced by a sum.
[Equation 31]
Y.sub.i=.sub.my.sub.i,t(.sub.m)(31)

(102) When the explanation given thus far is summed up, it is known that the number of photons observed can be expressed by the equation (32) below, when the surface expansion, the rotation during the projection and the X-ray wavelength dependence are taken into account. Here, b.sub.m=d.sub.m.sub.m.

(103) $\begin{array}{cc}\left[\mathrm{Equation}32\right]& \\ {\stackrel{\_}{Y}}_i\underset{m}{.Math.}\underset{\left\{k.Math.r\left(k\right)S_{i,t}\right\}}{.Math.}{b}_m\exp \left(-\underset{j{I}_k}{.Math.}{l}_{\mathrm{kj}}{x}_{\mathrm{jm}}\right)N& \left(32\right)\end{array}$

(104) Summing about m and k outside the exponential as expressed in the above formula (32) is called the multiple-X-ray sum. This cannot be expressed by a liner integral that does not assume the existence of a sum (or integral) outside the exponential. The algorithm based on the Bayesian inference disclosed in the specification of the present invention can be applied to observation process, in general, including such multiple X-ray summation.

INDUSTRIAL APPLICABILITY

(105) The present invention is useful not only for medical purpose but for industrial X-ray CT equipment that decreases beam hardening artifacts. Also, it is useful for belt conveyer type X-ray CT equipment that cannot ignore the shift and the rotation by a belt conveyer and the dental X-ray CT equipment.