WO2014021349A1 - Ｘ線コンピュータ断層撮影装置及び画像再構成方法 - Google Patents

Ｘ線コンピュータ断層撮影装置及び画像再構成方法 Download PDF

Info

Publication number
WO2014021349A1
WO2014021349A1 PCT/JP2013/070659 JP2013070659W WO2014021349A1 WO 2014021349 A1 WO2014021349 A1 WO 2014021349A1 JP 2013070659 W JP2013070659 W JP 2013070659W WO 2014021349 A1 WO2014021349 A1 WO 2014021349A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
noise
reconstruction
reconstruction filter
computed tomography
Prior art date
Application number
PCT/JP2013/070659
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 EP13825580.7A priority Critical patent/EP2881039A4/en
Priority to KR1020137031515A priority patent/KR101598265B1/ko
Priority to CN201380001535.4A priority patent/CN103732147B/zh
Publication of WO2014021349A1 publication Critical patent/WO2014021349A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • 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]
    • 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
    • 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
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/421—Filtered back projection [FBP]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative

Definitions

  • Embodiments described herein relate generally to an X-ray computed tomography apparatus and an image reconstruction method.
  • the filtered back projection (FBP) algorithm is simple and efficient and is used to reconstruct images in nuclear medicine, X-ray CT, and MRI.
  • the FBP algorithm is useful in X-ray CT and nuclear medicine image reconstruction because of its high computational efficiency.
  • FBP generates images that contain a lot of noise when applied to data collected with low-dose x-rays.
  • the FBP algorithm is gradually replaced by an iterative algorithm.
  • the FBP algorithm generates an image containing noise.
  • the FBP algorithm cannot incorporate a noise model to reduce the noise level.
  • the iterative algorithm can generate an image with less noise than the FBP algorithm by incorporating a projection noise model.
  • iterative algorithms require intensive calculations.
  • iterative algorithms can generate images that balance noise and resolution using maximum a posteriori probability (MAP).
  • MAP maximum a posteriori probability
  • FBP algorithm cannot take advantage of MAP or prior image information.
  • An object of the present invention is to provide an X-ray computed tomography apparatus and an image reconstruction method that use the FBP algorithm and can substantially reduce noise in an image without losing calculation efficiency.
  • the X-ray computed tomography apparatus includes at least one reconstruction based on a determination unit that determines a weight value for weighting projection data according to at least a predetermined noise model, and a parameter including the weight value
  • a shaping unit that shapes a filter
  • a generation unit that generates an image based on the reconstruction filter.
  • FIG. 1 is a diagram showing a configuration of an X-ray computed tomography apparatus according to the present embodiment.
  • FIG. 2 is a diagram illustrating a configuration of the reconstruction unit in FIG. 1.
  • FIG. 3 is a diagram showing a typical flow of processing by the reconstruction unit in FIG.
  • FIG. 4 is a diagram illustrating a noise weight value determination method for each view in the image reconstruction method according to the present embodiment.
  • FIG. 5 is a diagram illustrating a method for determining a noise weight value for each projection line in the image reconstruction method according to the present embodiment.
  • FIG. 6A is a diagram illustrating a typical flow of an image reconstruction method including an additional step for generating an image based on a pre-shaped reconstruction filter according to the present embodiment.
  • FIG. 1 is a diagram showing a configuration of an X-ray computed tomography apparatus according to the present embodiment.
  • FIG. 2 is a diagram illustrating a configuration of the reconstruction unit in FIG. 1.
  • FIG. 3 is a diagram showing a
  • FIG. 6B is a diagram illustrating a typical flow of another image reconstruction method including an additional step for generating an image based on a pre-shaped reconstruction filter according to the present embodiment.
  • FIG. 7 is a diagram showing a typical flow of an image reconstruction method using a noise weighted ramp filter in the FBP algorithm according to the present embodiment.
  • FIG. 8 is a diagram showing a typical flow of an image reconstruction method using a set of noise weighted ramp filters in the edge preserving FBP-MAP algorithm according to the present embodiment.
  • FIG. 9 is a diagram showing a typical flow of the edge map generation processing of FIG.
  • FIG. 10 is a diagram illustrating a typical flow of an image reconstruction method based on the FBP-MAP algorithm using a prior image according to the present embodiment.
  • FIG. 11 is a diagram showing a typical flow of an image reconstruction method using the noise weighted iterative emulation FBP algorithm according to the present embodiment.
  • FIG. 12 is a diagram illustrating a concept of the influence of the above-described emulated iteration in the noise weighted FBP algorithm according to the present embodiment.
  • FIG. 13A shows an image of a transverse cut in a cadaver torso region reconstructed based on a conventional FBP felt-kamp algorithm.
  • FIG. 13D shows an image of a transverse cut in the torso region of a cadaver reconstructed based on an rFBP-MAP reconstruction with an edge-preserving prior, which is substantially the weighted addition of FIGS. 13B and 13C.
  • FIG. FIG. 13E is a diagram showing an edge map used in edge preserving FBP-MAP reconstruction.
  • FIG. 13F shows an image of a transverse cut in the torso region of a cadaver reconstructed based on vFBP-MAP reconstruction.
  • This embodiment relates to an X-ray computed tomography apparatus and a reconstruction processing method for processing projection data using a noise weighted filter-corrected backprojection method.
  • FIG. 1 is a diagram showing a configuration of an X-ray computed tomography apparatus according to the present embodiment.
  • the X-ray computed tomography apparatus according to this embodiment is a multi-slice X-ray CT apparatus or a multi-slice X-ray CT scanner.
  • the X-ray computed tomography apparatus according to this embodiment is equipped with a gantry 100.
  • the gantry 100 includes an X-ray tube 101, an annular frame 102, and a multi-row type or two-dimensional array type X-ray detector 103.
  • the X-ray tube 101 and the X-ray detector 103 are installed in the diameter direction across the subject S on the frame 102 that rotates around the axis RA.
  • the rotating unit 107 rotates the frame 102 at a high speed such as 0.4 seconds per rotation while the subject S is moved in the illustrated page or outside the illustrated page along the axis RA.
  • the X-ray computed tomography apparatus includes a high voltage generator 109.
  • the high voltage generator 109 applies a tube voltage to the X-ray tube 101 so that the X-ray tube 101 generates X-rays.
  • the high voltage generator 109 is installed on the frame 102.
  • X-rays are emitted toward the subject S.
  • the X-ray detector 103 is located on the opposite side of the X-ray tube 101 across the subject S in order to detect X-rays transmitted through the subject S.
  • the X-ray computed tomography apparatus further includes another device for processing a signal detected from the X-ray detector 103.
  • the data acquisition circuit (DAS) 104 converts the signal output from the X-ray detector 103 into a voltage signal for each channel, amplifies the voltage signal, and converts the amplified voltage signal into a digital signal.
  • X-ray detector 103 and DAS 104 are configured to process a predetermined total number of projections (TPRP) per revolution.
  • the above-described data is transmitted to the preprocessing unit 106 stored in the console outside the gantry 100 via the non-contact data transmission unit 105.
  • the preprocessing unit 106 performs preprocessing such as sensitivity correction on the raw data.
  • the preprocessed data is called projection data.
  • the storage unit 112 stores projection data.
  • the storage unit 112 is connected to the system control unit 110 via the data / control bus together with the reconstruction unit 114, the display unit 116, the input unit 115, and the scan plan support apparatus 200.
  • the scan plan support apparatus 200 has a function for supporting an imaging technique in order to construct a scan plan.
  • the reconstruction unit 114 reconstructs an image from the projection data stored in the storage unit 112 based on a filtered back projection (FBP) technique using noise weighting.
  • the reconstructing unit 114 may reconstruct an image from projection data based on filtered weighted back projection and filtered correction back projection that uses a prior distribution input such as a reference image.
  • Any one of the above two embodiments of the reconstructor 114 can be applied to a filtered backprojection method with an additional feature that emulates a specific iteration result in a predetermined number of iterations according to a predetermined iterative reconstruction algorithm. Based on this, an image is reconstructed from the projection data.
  • the reconfiguration unit 114 is configured by a combination of software and hardware, and is not limited to a specific implementation form.
  • the reconstruction unit 114 according to the present embodiment can be implemented in other modalities such as nuclear medicine and magnetic resonance imaging (MRI).
  • FIG. 2 is a diagram illustrating a configuration of the reconstruction unit 114 according to the present embodiment.
  • the reconstruction unit 114 includes a noise weight determination unit 114A, a filter shaping unit 114B, and an image generation unit 11C.
  • the noise weight determination unit 114A determines a weight value for weighting projection data according to at least a predetermined noise model.
  • the filter shaping unit 114B shapes at least one reconstruction filter based on the parameter including the determined weight value.
  • the image generation unit 114C generates an image based on the shaped reconstruction filter. It should be noted that the implementation form of the reconfiguration unit 114 is merely an example.
  • the noise weight determination unit 114A, the filter shaping unit 114B, and the image generation unit 11C will be described in detail.
  • FIG. 3 is a diagram showing a typical flow of processing by the reconfiguration unit 114.
  • the noise weight determination unit 114A determines a weight value for weighting the projection data according to at least a predetermined noise model. The weight value is determined for each predetermined data unit for the projection data.
  • the data unit may be a view or a projection line.
  • the filter shaping unit 114B shapes at least one reconstruction filter based on the parameter including the weight value determined in step S100. For the shaping according to this embodiment, a plurality of reconstruction filters are optionally used, while these parameters are not limited to weight values.
  • the image generation unit 114C generates an image based on the reconstruction filter shaped in step S120.
  • Step S120 shapes the filter based on a predetermined parameter.
  • one exemplary parameter is the variance of noise components included in the projection data (hereinafter referred to as noise variance).
  • the parameter including the weight value according to the present embodiment is not limited to the noise parameter.
  • the reconstruction filter may be shaped based on a system matrix. More specifically, the system matrix specifies non-uniform sampling of projection data.
  • the conventional FBP algorithm assumes that the objects are sampled uniformly, it is possible to use a variable sampling strategy. For example, signals within a certain angular range containing important structures are sampled at a first angular interval. In contrast, signals within another angular range outside the region are sampled at a second angular interval that is greater than the first angular interval.
  • step S100 does not have to be performed or completed before step S120 when determining the weight value for each of the predetermined data portions of the projection data. That is, the weight value may be optionally determined each time the reconstruction filter is shaped.
  • step S140 the image generation unit 114C executes a predetermined set of sub-steps while applying the reconstruction filter shaped in step S120.
  • the image is reconstructed by applying a shaped filter to the projection data in the frequency domain after Fourier transform in a given FBP technique, as described in further detail.
  • the image is reconstructed by applying a filter customized as a convolution in a predetermined FBP technique to the projection data in the spatial domain.
  • FIG. 4 is a diagram illustrating a noise weight value determination method for each view in the image reconstruction method according to the present embodiment.
  • each view is generally assigned a weight value according to a predetermined function, such as maximum noise variance or average noise variance in a particular view.
  • FIG. 4 illustrates an X-ray source at three positions, namely positions A, B, and C, and an X-ray detector at three corresponding positions at positions A ′, B ′, and C ′. For example, when the X-ray source is at position A, the X-ray detector is at the corresponding position A ′.
  • the projection data collected by the X-ray detector at the position A ′ comprises a view, and the weight value w ( ⁇ 1 ) is determined according to the view angle ⁇ 1 in accordance with a predetermined weighting method for the view.
  • the determination of per-view weight values is performed and completed in one exemplary process prior to shaping the reconstruction filter.
  • the determination of this per-view weight value is determined each time the reconstruction filter is shaped in another exemplary process.
  • the noise variance included in the projection data can be modeled using a weighting method in units of views.
  • a typical approach for reconstructing an image while compensating for the effects of noise is to minimize an objective function weighted according to noise variance, defined by equation (1).
  • A is a projection matrix
  • X is an image array described as a column vector
  • P is a projection array described as a column vector
  • W is a noise weighting matrix that defines a weighting norm. is there.
  • an iterative Landweber algorithm defined in equation (2) is used.
  • Equation (2) AT is a back projection matrix
  • X (k) is an image estimated in the k-th iteration
  • ⁇ is a step size
  • equation (3) is simplified as equation (4).
  • Equation (4) corresponds to a noise weighted FBP algorithm, and the filter function of Equation (4) in the frequency domain is expressed as Equation (5).
  • ⁇ is frequency dispersion
  • w (view) is a weight value determined in view units. That is, w (view) is a noise-related weight function at a specific view angle.
  • ⁇ is a parameter that emulates the step size in the iterative algorithm
  • k is a parameter that emulates the number of iterations in the iterative algorithm.
  • the view-by-view weighted FBP algorithm denoted as vFBP, has a shift invariant point response function (PRF).
  • the reconstruction filter is shaped to emulate a specific iteration result in a predetermined number k iterations according to a predetermined iterative reconstruction algorithm.
  • the above features for emulating iterations in a reconstruction filter are not necessarily included in other embodiments for reconstructing an image from projection data using a per-view weighted FBP algorithm.
  • FIG. 5 is a diagram illustrating a method for determining a noise weight value for each projection line in the image reconstruction method according to the present embodiment.
  • a weight value is assigned to each projection line according to a predetermined function such as maximum noise variance or average noise variance in a specific view.
  • FIG. 5 illustrates an X-ray source at a predetermined position and an X-ray detector at a corresponding position. When the projection line is projected from the X-ray source, the projection line reaches the X-ray detector at the corresponding position.
  • the projection data consists of a plurality of projection lines within the view angle at ⁇ .
  • the weight value w ( ⁇ , t 1 ) is determined for a specific projection line t 1 according to the view angle ⁇ and according to a predetermined weighting method for the projection line. As described above, the determination of projection line weight values is performed and completed in one embodiment prior to shaping the reconstruction filter. The projection line weight value is determined as a reconstruction filter in other embodiments.
  • the weight value that is, the weighting factor w (ray) is a weighting factor determined by the projection line.
  • w (ray) is determined according to the noise variance of the associated projection line.
  • k is a parameter indicating a predetermined number of iterations for emulating a particular iteration result according to a predetermined iterative reconstruction algorithm. If the inverse Fourier transform of H ( ⁇ , ray) is h (t, ray), h (t, ray) is the kernel in the spatial domain of the filter.
  • p (r, ⁇ ) be the projection at view ⁇ and projection line r
  • q (r, ⁇ ) be the filtered projection.
  • q (r, ⁇ ) is defined by integration as shown in the following equation (7).
  • Kernel h (t, r) depends on r, so this integration is not convolution.
  • the computational cost is equivalent to convolution.
  • the backprojection step in this new rFBP algorithm is the same as the backprojection step of the conventional FBP algorithm.
  • the rFBP algorithm is therefore computationally efficient.
  • image domain filtering is performed in the two-dimensional Fourier domain using a transfer function, as described later in the implementation for the noise weighted FBP-MAP (rFBP-MAP) algorithm for each projection line.
  • a typical method for assigning the weighting factor is to make w (ray) the reciprocal of the noise variance of the projection line measurement. This approach is supplemented using a likelihood function (ie, joint probability density function) as the objective function for the optimization problem.
  • a likelihood function ie, joint probability density function
  • a weighting factor can be arbitrarily assigned under the restriction that (ray) is assigned.
  • Equation (7) The easiest and more efficient way to implement the rFBP algorithm according to this embodiment is to use equation (7) to filter the projection in the spatial domain.
  • equation (7) An alternative way to calculate q (r, ⁇ ) for each view ⁇ is to perform equation (7) in the frequency domain as follows:
  • Step 1 Find a one-dimensional Fourier transform of p (r, ⁇ ) with respect to r and obtain P ( ⁇ , ⁇ ).
  • Step 2 For each projection line “ray”, a weighting factor w (ray) is assigned, and a frequency domain transfer function is formed as represented by equation (6).
  • ⁇ is a frequency index and takes an integer.
  • Step 2. b: Find the product, Q ray ( ⁇ , ⁇ ) P ( ⁇ , ⁇ ) H ( ⁇ , ray).
  • Step 2. c Find the one-dimensional inverse Fourier transform of Q ray ( ⁇ , ⁇ ) with respect to ⁇ , and obtain q (ray, ⁇ ).
  • step 1 processes all r in each view column, while step 2. From a. c processes only one projection line at a time.
  • a single weight value i.e., a weighting factor w (view)
  • a separate weight value i.e., a weighting factor w (ray)
  • the noise weighting scheme is more accurate in the per-projection weighted FBP (rFBP) algorithm than the per-view weighted FBP (vFBP) algorithm.
  • the present embodiment is based on at least a parameter including a weight value before generating an image.
  • Shape one reconstruction filter is shaped to emulate a specific iteration result in a predetermined number of iterations denoted k, according to a predetermined iterative reconstruction algorithm.
  • the above features for emulating iterations in a reconstruction filter are included in other embodiments for reconstructing an image from projection data using a weighted per-view FBP algorithm or a weighted per-projection FBP algorithm. Not always.
  • noise reduction in image reconstruction is based on the concept that image reconstruction, particularly backprojection, is an addition process and that appropriate weighting of projection data substantially reduces backprojection variance.
  • FIG. 6A is a diagram illustrating a typical flow of an image reconstruction method including an additional step for generating an image based on a pre-shaped reconstruction filter.
  • a reconstruction filter is applied to the projection data before backprojection is performed.
  • the reconstruction filter Prior to step S140A, the reconstruction filter is shaped in advance based on parameters including noise weight values.
  • the reconstruction filter shaped in step S140A is applied to the projection data. If the projection data is subjected to Fourier transformation before step S140A, the filtered projection data undergoes inverse Fourier transformation before step S140B.
  • the projection data filtered in step S140B is backprojected. Back projection produces an image based on the filtered projection data.
  • FIG. 6B is a diagram showing a typical flow of another image reconstruction method including an additional step for generating an image based on a pre-shaped reconstruction filter.
  • a reconstruction filter is applied to the projection data after the backprojection is performed. If the projection data is subjected to Fourier transform before step S140B, the projection data undergoes inverse Fourier transform before step S140B. In step S140B, the projection data is backprojected to generate an image. The shaped reconstruction filter is applied to the backprojected data. Prior to step S140A, the reconstruction filter is pre-shaped based on parameters including noise weight values. In step S140A, the shaped reconstruction filter is applied to the projection data.
  • FIG. 7 is a diagram showing a typical flow of an image reconstruction method using a noise weighted ramp filter in the FBP algorithm.
  • noise weighting filtering is performed in the Fourier domain before the back projection.
  • the projection data PD is collected in a plurality of views.
  • the projection data PD is processed in units of views or projection lines.
  • the projection data PD is subjected to the Fourier transform FT so that the transformed data exists in the Fourier domain or the frequency domain.
  • the ramp filter is shaped according to a predetermined set of parameters including noise weight values. Thereby, the weighted ramp filter WF is formed.
  • the plurality of ramp filters are shaped to be a plurality of weighted ramp filters WF.
  • the transformed projection data is filtered according to a noise weighted ramp filter NWF in a predetermined FBP algorithm.
  • the data filtered in the frequency domain undergoes an inverse Fourier transform IFT so that the projection data is backprojected in the spatial domain at the backprojection BP.
  • another noise weighted ramp filter NWF2 may be applied to the image. In this case, the noise weighted ramp filter NWF2 applied to the image is different from the noise weighted ramp filter NWF applied to the converted data.
  • Backprojection BP normalizes the variable angle sampling described above using a weight function during image generation. For example, signals within a certain angular range containing important structures are sampled at a first angular interval. Signals in other angular ranges outside the region are sampled at a second angular interval that is greater than the first angular interval.
  • the backprojection BP is an integral over the sampling angle, and angular sampling density compensation is achieved by using a normalization factor that is essentially a Jacobian factor in the backprojection integral. Mathematically, the backprojection image is expressed as in equation (8).
  • equation (8) when the projection is first backprojected, p is the projection at the angle ⁇ .
  • p is the projection filtered at the angle ⁇ and t indicates the position of the detector bin.
  • Equation (9) The discrete expression of equation (8) is expressed as in equation (9).
  • n is the detector position index
  • m is the projection angle index
  • M is the total number of views from which projections were collected
  • INT is used to indicate nearest neighbor interpolation.
  • the backprojection depends on a density function d (0) that is the number of views per unit angle. For example, for 0 ⁇ ⁇ / 4, if the sampling interval is 1 ° and 3 ° for ⁇ / 4 ⁇ ⁇ , the density function is expressed by equation (10) and backprojection (9) The expression is changed to the expression (11).
  • the sampling density function is a function of the angle index in, instead of the actual angle ⁇ .
  • FIG. 8 is a diagram showing a typical flow of an image reconstruction method using a set of noise weighted ramp filters in a predetermined edge-preserving FBP-MAP algorithm according to the present embodiment.
  • noise weighting filtering is performed in the Fourier domain before the image is reconstructed by back projecting a set of projection data. Finally, a final image is generated based on the two images.
  • the projection data PD is collected in a plurality of views.
  • the projection data PD may be processed in view units or may be processed in projection line units.
  • the projection data PD undergoes a Fourier transform FT so that the transformed data is in the Fourier domain or the frequency domain.
  • the two ramp filters are individually shaped according to a predetermined set of parameters including noise weight values and shaped into weighted ramp filters WF- ⁇ 0 and WF- ⁇ .
  • the ramp filter can be further shaped by normalization.
  • the filter is a function of the value ⁇ .
  • the larger the value ⁇ the stronger the smoothing of the image. For example, when smoothing is not performed or when there is no blur, the ⁇ value is selected to be zero.
  • the converted projection data may be noise-weighted in a predetermined FBP algorithm.
  • the noise weighted projection data is filtered by normalized ramp filters WF- ⁇ 0 and WF- ⁇ , respectively. Thereafter, the projection data filtered in the frequency domain is subjected to an inverse Fourier transform IFT.
  • the projection data after the inverse Fourier transform is back-projected in the spatial domain so that two images X 0 and X ⁇ are generated.
  • the image X 0 back-projected without smoothing or back-projected with minimal smoothing is used to extract edge map information B in the prepared edge map EMAP.
  • the A final image FIMG is generated according to certain rules based on the minimally smoothed image X 0 , the filtered image X ⁇ , and the extracted edge information.
  • edges are preserved during denoising using Bayesian normalization terms in the objective function shown in equation (12) below.
  • the penalty function g (X) in equation (12) measures the step of the pixel value in the vicinity of the target pixel.
  • A is a projection matrix
  • X is an image array written as a column vector
  • P is a projection array written as a column vector
  • W is a diagonal composed of a weight function for each projection.
  • I is a matrix
  • b is a fidelity term that is the first term of equation (12).
  • I a relative weighting factor that adjusts the importance of the Bayesian term g (X) for.
  • Penalty g (X) is a quadratic function of the pixel value step.
  • the pixel value step is calculated as the difference between the pixel value of the pixel and the average value of neighboring pixels.
  • a large pixel value step corresponds to a large penalty.
  • the penalty function results in a smoothed image.
  • the penalty function g (X) should have different characteristics with respect to the edge and non-edge regions.
  • a Hoover function having a threshold is used as a penalty function having such ability.
  • the Hoover function is a quadratic expression and is smoothed. If the pixel value step is greater than the threshold and the pixel value step is classified as an edge, the Hoover function is linear and performs a smaller smoothing.
  • TV (total variation) measurement Another common prior distribution that encourages constant images per part is TV (total variation) measurement.
  • the gradient of TV measurement of g (X) takes only three values. That is, if g (X) is greater than its left and right neighborhoods, it is a positive value, and if g (X) is less than its left and right neighborhoods, it is a negative value and g (X) is between its left and right neighborhoods. In this case, it is 0.
  • a positive gradient value pushes down the value of g (X)
  • a negative gradient value pushes up the value of g (X)
  • a gradient value of 0 does not change anything Means not done.
  • the TV prior distribution encourages monotonic functions, suppresses amplitude, and preserves sharp edges.
  • the image reconstruction method applies an edge-preserving FBP-MAP algorithm to projection data.
  • ⁇ is set to ⁇ 0 for one noise weighted ramp filter, ie zero.
  • ⁇ 0 g (X) is zero and no smoothing is performed. Thereafter, the image X 0 is reconstructed without decorated any smoothing.
  • ⁇ is set to ⁇ for other noise weighted ramp filters, ie non-zero. A certain amount of smoothing is performed by ⁇ , and the image X ⁇ is reconstructed. Thus, two images X beta and X 0 is reconstructed.
  • the first image X ⁇ is reconstructed based on rFBP-MAP reconstruction using a non-zero ⁇ value
  • the second image X 0 is reconstructed based on rFBP reconstruction using a zero ⁇ value. Is done.
  • FIG. 9 is a diagram showing a typical flow of processing for generating the edge map B of FIG.
  • the image X 0 is reconstructed by the rFBP reconstruction method using the zero ⁇ value.
  • the pixel value step is evaluated by edge detection using essentially a Sobel kernel which is the norm of the image gradient.
  • the edge map generated by the Sobel kernel is then converted to a binary image using value 0 or value 1 according to a predetermined threshold in step S220.
  • An image value exceeding the threshold is set to a value 1, and an image value less than the threshold is set to a value 0.
  • the binary image is blurred by a predetermined low-pass filter, and an edge map B is generated.
  • the reason for blurring the binary image is to smooth the boundary between the edge region and the non-edge region included in the binary image.
  • Final image as shown in (13) below, is generated by a linear combination of the X b and X 0.
  • FIG. 10 is a diagram illustrating a typical flow of an image reconstruction method based on the FBP-MAP algorithm using a prior image. In general, a given prior image
  • the objective function (8) is generalized by the following equation (14).
  • the equivalent FBP-MAP algorithm has two Fourier domain filter functions expressed by the following equations (16) and (17).
  • the FBP-MAP algorithm has two inputs: projection data PD and a prior image PI.
  • the final image FI is an addition image based on the first image derived from the projection data PD and the second image derived from the preliminary image PI.
  • the FBP-MAP reconstruction generates a first reconstructed image based on a detector Fourier domain filter such as equation (16) and a subsequent backprojection.
  • the image Fourier domain processing IFD generates a second image using an image Fourier domain filter as shown in equation (17).
  • the image Fourier domain filter may be two-dimensional or three-dimensional.
  • FIG. 11 is a diagram showing a typical flow of an image reconstruction method using the noise weighted iterative emulation FBP algorithm according to the present embodiment.
  • a reconstruction filter mimics a specific iteration result in a predetermined number of iterations.
  • the reconstruction filter is shaped based on parameters including weight values according to a predetermined iterative reconstruction algorithm.
  • the noise weighted iterative emulation FBP algorithm 200D collects projection data at step 201.
  • Projection data is collected in the form of a Radon transform that represents the line integral of the object being imaged.
  • the projection data is collected by an arbitrary modality such as an X-ray computed tomography apparatus, a PET apparatus, or a SPECT apparatus as exemplified in FIG.
  • At least one reconstruction filter is shaped based on a combination of parameters including the number of iterations, noise weight values, system matrix, and the like. These parameters can be arbitrarily selected.
  • the parameter k is selected to emulate a specific iteration result in a predetermined number of iterations according to a predetermined iterative reconstruction algorithm.
  • a predetermined weight value W may be used together with an iterative parameter k for shaping a reconstruction filter applied to projection data.
  • Other filter parameters include the step size in the iterative algorithm and the relative weighting of Bayesian conditions so that different applications can be addressed.
  • step S202 After the reconstruction filter is shaped according to the combination of parameters in step S202, the shaped reconstruction filter is applied to the projection data in step 203.
  • the filtered data is backprojected in step 204.
  • An image is generated by back projection.
  • step S205 it is determined whether to repeat the emulation. When the emulation is repeated, the process proceeds to step S202 again for the emulation repetition. If the emulation is not repeated in step S205, the noise weighted iterative emulation FBP algorithm according to the present embodiment ends.
  • FIG. 12 is a diagram showing the concept of the influence of the above-described emulated iteration in the noise weighted FBP algorithm according to the present embodiment.
  • noise control is performed using pre-filtering or post-filtering.
  • noise is substantially reduced in image reconstruction using a noise weighting and iteratively emulated FBP algorithm.
  • Noise reduction in image reconstruction is based on the philosophy that image reconstruction (particularly backprojection) is an addition process, and that dispersion of backprojection can be reduced by performing weighting appropriate to projection.
  • two projections that is, measurement data is represented by two projection lines P1 and P2
  • a solution that is, an image I is an intersection of the two projection lines P1 and P2.
  • the horizontal axis X 1 and the vertical axis Y 1 represent two variables in the image I.
  • the two image pixels can take values at any location along a single straight line.
  • the values of the two pixels that are the intersection of the two straight lines are determined. If a third measurement is made and all measurements contain noise, the three lines may not intersect at one point. In general, if an image consists of many pixels and measurements with many contradictory noises are made, these measurements may not intersect at a single point, so there may be no solution to the imaging problem .
  • the FBP algorithm can provide a single noisy solution from a set of noisy projections.
  • the FBP algorithm is modified to emulate multiple iteration steps of the iterative algorithm to produce an image with less noise.
  • the iterative algorithm is usually terminated before a true solution is achieved. This is especially true when the true solution involves noise because the measurement data contains noise due to low photon counting.
  • the points NNW1 and NNW2 on the straight line show the intermediate results after the first predetermined number of iterations and the second predetermined number of iterations, respectively.
  • the second predetermined number of iterations is larger than the first predetermined number of iterations, and the intermediate result NNW2 is emulated after the second predetermined number of iterations.
  • the intermediate results NW1, NW2, and NW3 are located on a curve, which indicates that their emulated iterations are weighted by noise variance.
  • the intermediate results NW1, NW2, and NW3 are after the third predetermined number of iterations, the fourth predetermined number of iterations, and the fifth predetermined number of iterations, respectively.
  • the fifth number is greater than the third number of iterations and the fourth number of iterations, and the third number is the smallest.
  • the intermediate result NW3 is emulated after the fifth iteration.
  • This solution trajectory is determined by a set of weighting factors. One weighting factor is assigned to each measurement. A larger factor means that the associated measurement is more reliable.
  • a common approach is to assign a weighting factor that is proportional to the inverse of the noise variance of the measurement.
  • a pseudo solution with less noise in the vicinity of a solution with true noise is acceptable.
  • the pseudo-solution procedure is parameterized by a predetermined set of parameters including a noise variance weight value and an iteration index k.
  • a pseudo solution with a smaller k corresponds to a smoother image
  • a pseudo solution with a larger k corresponds to an image that is clearer but with more noise.
  • One embodiment utilizes noise weighted FBP per projection line, using a maximum posterior probability (MAP) algorithm specified by rFBP-MAP.
  • MAP maximum posterior probability
  • the derivation step for the rFBP-MAP algorithm is very similar to the derivation step for the per-view weighted FBP-MAP algorithm. Some main derivation steps are briefly described below.
  • the objective function is the same as the equation (12).
  • the Bayesian condition g (X) is quadratic
  • the gradient of g (X) is in the form of RX with respect to several matrices R, and the solution is obtained by the Landweber algorithm expressed by equation (18). It is done.
  • Equation (18) has an equivalent non-iterative formula.
  • the non-recursive formula is expressed by the formula (19).
  • a Fourier domain notation is derived from the determinant so that an FBP type algorithm is obtained.
  • the matrix A ⁇ 1 is processed such that one-dimensional ramp filtering is followed by backprojection, or backprojection is followed by two-dimensional ramp filtering.
  • determinant (21) is considered as a “filter after backprojection” algorithm
  • image domain filtering in (21) is performed in the two-dimensional Fourier domain using the transfer function shown in (22).
  • v x and v y are frequencies related to x and y, respectively.
  • I a two-dimensional frequency vector.
  • the projection operator A is a line integral in two-dimensional space (ie, Radon transform) and AT is an operator (ie, backprojection transform), then projection and backprojection
  • the combined operator A T A is a two-dimensional convolution with the two-dimensional kernel 1 / r of the original image. In this case, in xy Cartesian coordinates
  • the two-dimensional filtering of the image area is converted into a one-dimensional filtering of the projection area.
  • the transfer function in the Fourier domain of the two-dimensional image (equation (22)) is equal to the transfer function in the Fourier domain of the one-dimensional projection as shown in the following equation (23).
  • ⁇ is a frequency related to a variable along the detector.
  • L 1D is a central cross section of L 2D in the equation (22).
  • the weighting factor w depending on the noise in the equations (22) and (23) is not a constant.
  • w w (view) depending on the view angle.
  • w w (ray) according to the projection line.
  • a common technique for assigning weighting factors is to make w (ray) proportional to the inverse of noise variance in projection line measurement. This approach is supplemented by using a likelihood function (ie, joint probability density function) as the objective function for the optimization problem.
  • the fundamental principle is that measurements with less noise should be relied upon than measurements with greater noise.
  • the modified ramp filter (Equation (23)) becomes a conventional ramp filter, and the rFBP-MAP algorithm is reduced to the conventional FBP algorithm.
  • the inverse Fourier transform of is the spatial domain kernel of a one-dimensional modified ramp filter
  • Step 1 Find a one-dimensional Fourier transform of p (t, ⁇ ) with respect to t and obtain P ( ⁇ , ⁇ ).
  • Step 2 Assign a weighting factor w (t) to each projection line t to form a frequency domain transfer function (Equation (23)).
  • ⁇ is a discrete frequency index, 0, 1, 2,. . . Takes an integer value such as
  • Step 4 Find the one-dimensional inverse Fourier transform of Q 1 ( ⁇ , ⁇ ) with respect to ⁇ to obtain q (t, ⁇ ).
  • Step 1 processes all t in each view, while steps 2 through 4 process only one t at a time.
  • the above-described rFBP-MAP algorithm is not limited only to imaging geometry using parallel beams.
  • This rFBP-MAP algorithm can be extended to other imaging geometries such as fan beams, cone beams, plane integration, and attenuation measurements as long as the FBP algorithm exists.
  • the only modification to the FBP algorithm is projection data filtering. A weighting factor is assigned to each measurement, and the filter transfer function in the frequency domain is formed as equation (23) that is independent of the imaging geometry.
  • a cadaver torso was scanned using a low-dose X-ray computed tomography system. Images were reconstructed using the FBP algorithm as well as the rFBP-MAP algorithm.
  • the imaging geometry is a cone beam, and the X-ray source trajectory is a circle with a radius of 600 mm.
  • the X-ray detector had 64 rows, the row height was 0.5 mm, each row had 896 channels, and the fan angle was 49.2 °.
  • the tube voltage was 120 kV and the current was 60 mA. 1200 views were sampled uniformly over 360 °.
  • the data was first weighted using a cosine function and then a one-dimensional ramp filter was applied to each column of the cone beam projection. Finally, cone beam backprojection was used to generate 3D volume data.
  • the one-dimensional ramp filter was replaced with a ramp filter as represented by equation (23).
  • the iteration index k was selected as 20,000 and the step size ⁇ was set to 2.0.
  • the Bayesian penalty function g (X) was a common quadratic function of the difference between the central pixel value in the cross section and the average value of its four neighbors.
  • R is simply a Laplacian that is a second derivative using a kernel whose one-dimensional version is [ ⁇ 1 2 ⁇ 1].
  • w (t, ⁇ ) exp ( ⁇ p (t, ⁇ )).
  • the transmission measurement approximately follows a Poisson distribution, and the line integral p (t, ⁇ ) is a logarithm of the transmission measurement.
  • the edge map is acquired for each cross section using a 2D Sobel kernel.
  • the threshold for edge detection was set to 70% of the maximum pixel value.
  • a 3 ⁇ 3 smooth kernel is then applied to the binary image. By applying the smooth kernel, the edge included in the binary image is blurred, and the edge map B shown in Expression (13) is generated from the binary image.
  • the image volume data is reconstructed with a three-dimensional array of 512 ⁇ 512 ⁇ 512, and the non-center axis fragment is used for display. No other pre-filtering or post-filtering was applied to this reconstruction.
  • FIG. 13A shows an image of a transverse cut in a cadaver torso region reconstructed based on a conventional FBP felt-kamp algorithm. Due to the low dose, there is a significant streak artifact from left to right across the fuselage.
  • the noise weighting for each projection line in the rFBP algorithm effectively removed the streak artifact that appeared in FIG. 13A.
  • the streak artifact disappears due to noise weighting. Due to the prior distribution of the blur, the image looks excessively smooth, especially in the edge region.
  • FIG. 13D shows an image of a transverse cut in the torso region of a cadaver reconstructed based on an rFBP-MAP reconstruction with an edge-preserving prior, which is substantially the weighted addition of FIGS. 13B and 13C.
  • the edge map was a binary image.
  • a low pass filter was used to blur the edge map so that the transition from the edge region to the non-edge region is smooth.
  • the image of FIG. 13D has less noise than FIG. 13B, but retains the main edges of FIG. 13B in a non-smooth state.
  • FIG. 13E is a diagram showing an edge map B used in edge preservation FBP-MAP reconstruction.
  • FIG. 13F shows an image of a transverse cut in the torso region of a cadaver reconstructed based on vFBP-MAP reconstruction.
  • all parameters except the noise weighting function w were set to the same values as the reconstruction parameters in FIG. 13D.
  • the noise model is defined by Poisson statistics.
  • the noise model is defined by a composite Poisson statistic of photon intensity measurement that includes electronic noise described by a Gaussian distribution.
  • a composite Poisson distribution is a probability distribution of the sum of independent uniformly distributed random variables, and the number of variables follows the Poisson distribution.
  • the Poisson distribution is suitable for describing the number of incident photons.
  • the number N of photons reaching the detector during the storage time is N to Poisson ( ⁇ ), where ⁇ is the count rate.
  • Each detected photon has a random energy X i and its distribution is explained by the polychromatic spectrum of the X-ray tube and the detector efficiency.
  • the integral energy Y follows the composite Poisson statistics (mean value and variance) as shown in equations (25) and (26).
  • Equation (27) The concept of dispersion gain is indicated by W and is defined as shown in Equation (27).
  • the noise dispersion gain is proportional to the average effective energy of the spectrum of the incident X-rays to the X-ray detector.
  • coefficients that affect the average spectral energy there are the following coefficients that affect the average spectral energy:
  • Tube voltage (kVp). Usually, the tube voltage varies from 80 kVp to 140 kVp and determines the upper limit of the spectral energy.
  • Equation (29) g is a detector gain coefficient, and e is electronic noise with a variance V e .
  • V e is a detector gain coefficient
  • e is electronic noise with a variance V e .
  • the X-ray tube is switched off and data is collected to measure the dispersion.
  • the mean and variance are measured in the time direction independently for each detector pixel.
  • the statistics of the collected data are obtained by equation (30).
  • equation (31) is obtained.
  • DESCRIPTION OF SYMBOLS 100 ... Gantry, 101 ... X-ray tube, 102 ... Frame, 103 ... X-ray detector, 104 ... Data acquisition circuit, 105 ... Non-contact data transmission part, 106 ... Pre-processing part, 107 ... Rotation part, 108 ... Slip ring , 109 ... high voltage generation unit, 110 ... system control unit, 112 ... storage unit, 114 ... reconstruction unit, 114A ... noise weight determination unit, 114B ... filter shaping unit, 114C ... image generation unit, 115 ... input unit, 116 ... Display unit, 117 ... Current adjustment unit, 200 ... Scanning plan support device

Landscapes

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

Abstract

　ＦＢＰアルゴリズムを使用し、計算効率を失わずに画像内のノイズを実質的に低減する。ノイズ重み決定部１１４Ａは、少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する。フィルタ整形部１１４Ｂは、重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する。画像発生部１１４Ｃは、再構成フィルタに基づいて画像を発生する。

Description

Ｘ線コンピュータ断層撮影装置及び画像再構成方法
　本発明の実施形態は、Ｘ線コンピュータ断層撮影装置及び画像再構成方法に関する。
　フィルタ補正逆投影（ＦＢＰ）アルゴリズムは、簡単かつ効率的であり、核医学、Ｘ線ＣＴ、及びＭＲＩにおいて画像を再構成するために使用される。ＦＢＰアルゴリズムは、計算効率が高いため、Ｘ線ＣＴ及び核医学の画像再構成において重宝されている。ＦＢＰは、低線量のＸ線で収集されたデータに適用されると多くのノイズを含む画像を発生してしまう。
　一般的な傾向として、ＦＢＰアルゴリズムは、反復アルゴリズムに徐々に置換されている。ＦＢＰアルゴリズムは、ノイズを含む画像を発生してしまう。さらに、ＦＢＰアルゴリズムは、ノイズレベルを低減するためにノイズモデルを組み込むことが可能でないことが知られている。この点で、反復アルゴリズムは、投影ノイズモデルを組み込むことにより、ＦＢＰアルゴリズムよりもノイズがより少ない画像を発生することができる。しかしながら、反復アルゴリズムは、集中的な計算を必要とする。
　計算要件にかかわらず、反復アルゴリズムは、最大事後確率（ＭＡＰ）を使用してノイズと分解能とのバランスがとれた画像を発生することができる。対照的に、ＦＢＰアルゴリズムは、ＭＡＰまたは事前の画像情報を活用することができないことが知られている。
　目的は、ＦＢＰアルゴリズムを使用し、計算効率を失わずに画像内のノイズを実質的に低減可能なＸ線コンピュータ断層撮影装置及び画像再構成方法を提供することにある。
　本実施形態に係るＸ線コンピュータ断層撮影装置は、少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、前記重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する整形部と、前記再構成フィルタに基づいて画像を発生する発生部と、を具備する。
　ＦＢＰアルゴリズムを使用し、計算効率を失わずに画像内のノイズを実質的に低減すること。
図１は、本実施形態に係るＸ線コンピュータ断層撮影装置の構成を示す図である。 図２は、図１の再構成部の構成を示す図である。 図３は、図２の再構成部による処理の典型的な流れを示す図である。 図４は、本実施形態に係る画像再構成方法におけるビュー毎のノイズ重み値の決定方式を示す図である。 図５は、本実施形態に係る画像再構成方法における投影線毎のノイズ重み値の決定方式を示す図である。 図６Ａは、本実施形態に係る、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む画像再構成方法の典型的な流れを示す図である。 図６Ｂは、本実施形態に係る、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む他の画像再構成方法の典型的な流れを示す図である。 図７は、本実施形態に係る、ＦＢＰアルゴリズムでノイズ重み付けランプフィルタを使用した画像再構成方法の典型的な流れを示す図である。 図８は、本実施形態に係る、エッジ保存ＦＢＰ－ＭＡＰアルゴリズムでノイズ重み付けランプフィルタの組を使用した画像再構成方法の典型的な流れを示す図である。 図９は、図８のエッジマップの発生処理の典型的な流れを示す図である。 図１０は、本実施形態に係る、事前画像を使用したＦＢＰ－ＭＡＰアルゴリズムに基づく画像再構成方法の典型的な流れを示す図である。 図１１は、本実施形態に係るノイズ重み付け反復エミュレーションＦＢＰアルゴリズムを用いて画像再構成方法の典型的な流れを示す図である。 図１２は、本実施形態に係るノイズ重み付けＦＢＰアルゴリズムにおける上述のエミュレートされた反復の影響の概念を示す図である。 図１３Ａは、従来のＦＢＰフェルトカンプアルゴリズムに基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。 図１３Ｂは、β＝０を用いたｒＦＢＰ－ＭＡＰ再構成であるｒＦＢＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。 図１３Ｃは、二次平滑化事前分布を使用し、β＝０．０００５を用いたｒＦＢＰ－ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。 図１３Ｄは、実質的に図１３Ｂ及び図１３Ｃの重み付け加算である、エッジ保存事前分布を用いたｒＦＢＰ－ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。 図１３Ｅは、エッジ保存ＦＢＰ－ＭＡＰ再構成において使用されるエッジマップを示す図である。 図１３Ｆは、ｖＦＢＰ－ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。
　本実施形態は、ノイズ重み付けフィルタ補正逆投影（filtered back projection）法を使用して投影データを処理するためのＸ線コンピュータ断層撮影装置及び再構成処理方法に関する。
　以下、図面を参照しながら、本実施形態に係わるＸ線コンピュータ断層撮影装置及び画像再構成方法について説明する。
　図１は、本実施形態に係るＸ線コンピュータ断層撮影装置の構成を示す図である。図１に示すように、本実施形態に係るＸ線コンピュータ断層撮影装置は、マルチスライスＸ線ＣＴ装置またはマルチスライスＸ線ＣＴスキャナである。本実施形態に係るＸ線コンピュータ断層撮影装置は、ガントリ１００を装備している。ガントリ１００は、Ｘ線管１０１、円環状フレーム１０２、及び多列タイプあるいは二次元配列タイプのＸ線検出器１０３を含む。Ｘ線管１０１とＸ線検出器１０３とは、軸ＲＡの周りを回転するフレーム１０２上に被検体Ｓを横断して直径方向に設置される。回転部１０７は、被検体Ｓが軸ＲＡに沿って例示されたページ内または例示されたページ外に移動される間に、一回転当たり０．４秒等の高速でフレーム１０２を回転させる。
　本実施形態に係るＸ線コンピュータ断層撮影装置は、高電圧発生装置１０９を備える。高電圧発生装置１０９は、Ｘ線管１０１がＸ線を発生するように、管電圧をＸ線管１０１に印加する。例えば、高電圧発生器１０９はフレーム１０２上に設置される。Ｘ線は、被検体Ｓに向けて放射される。Ｘ線検出器１０３は、被検体Ｓを透過したＸ線を検出するために、被検体Ｓを横断してＸ線管１０１の逆側に位置する。
　さらに図１に示すように、本実施形態に係るＸ線コンピュータ断層撮影装置は、Ｘ線検出器１０３から検出された信号を処理するための他のデバイスをさらに含む。データ収集回路（ＤＡＳ）１０４は、Ｘ線検出器１０３からの信号出力をチャンネル毎に電圧信号に変換して、当該電圧信号を増幅し、増幅後の電圧信号をデジタル信号に変換する。Ｘ線検出器１０３とＤＡＳ１０４とは、１回転ごとに所定の総投影数（ＴＰＲＰ）を処理するように構成される。
　上述のデータは、非接触データ伝送部１０５を介してガントリ１００の外部のコンソール内に収納された前処理部１０６に伝送される。前処理部１０６は、生データに対して感度補正等の前処理を施す。前処理が施されたデータは、投影データと呼ばれている。記憶部１１２は、投影データを記憶する。記憶部１１２は、再構成部１１４、表示部１１６、入力部１１５、及びスキャン計画サポート装置２００と共に、データ／制御バスを介してシステム制御部１１０に接続される。スキャン計画サポート装置２００は、スキャン計画を構築するために、撮影技法をサポートするための機能を有している。
　本実施形態によれば、再構成部１１４は、ノイズ重み付けを用いたフィルタ補正逆投影（ＦＢＰ）技法に基づいて、記憶部１１２に記憶された投影データから画像を再構成する。あるいは、再構成部１１４は、ノイズ重み付けと基準画像等の事前分布入力（prior in）とも用いたフィルタ補正逆投影に基づいて、投影データから画像を再構成しても良い。再構成部１１４の上記の２つの実施形態のいずれか１つは、所定の反復再構成アルゴリズムに従って所定数の反復において特定の反復結果をエミュレートする追加の特徴を用いたフィルタ補正逆投影法に基づいて、投影データから画像を再構成する。
　再構成部１１４は、ソフトウェアとハードウェアとの組合せにより構成され、特定の実装形態に限定されない。本実施形態に係る再構成部１１４は、核医学及び磁気共鳴撮像（ＭＲＩ）等の他のモダリティにも実装可能である。
　図２は、本実施形態に係る再構成部１１４の構成を示す図である。図２に示すように再構成１１４は、ノイズ重み決定部１１４Ａ、フィルタ整形部１１４Ｂ、及び画像発生部１１Ｃを有する。ノイズ重み決定部１１４Ａは、少なくとも所定のノイズモデルに従って投影データに重み付けを施すための重み値を決定する。フィルタ整形部１１４Ｂは、決定された重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する。画像発生部１１４Ｃは、整形された再構成フィルタに基づいて画像を発生する。なお、この再構成部１１４の実装形態は単なる例示であることを言及しておく。以下、ノイズ重み決定部１１４Ａ、フィルタ整形部１１４Ｂ、及び画像発生部１１Ｃについて詳細に説明する。
　図３は、再構成部１１４による処理の典型的な流れを示す図である。ステップＳ１００においてノイズ重み決定部１１４Ａは、少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定する。重み値は、投影データについての所定のデータ単位毎に決定される。データ単位としては、ビューであっても良いし、投影線であっても良い。ステップＳ１２０においてフィルタ整形部１１４Ｂは、ステップＳ１００で決定された重み値を含むパラメータに基づいて、少なくとも１つの再構成フィルタを整形する。本実施形態に係る整形のために、複数の再構成フィルタがオプションで使用されると同時に、これらのパラメータは重み値に限定されない。ステップＳ１２０において画像発生部１１４Ｃは、ステップＳ１２０において整形された再構成フィルタに基づいて画像を発生する。
　ステップＳ１２０は、所定のパラメータに基づいてフィルタを整形する。上記のように、１つの例示的なパラメータは、投影データに含まれるノイズ成分の分散（以下、ノイズ分散と呼ぶ）である。しかしながら、本実施形態に係る重み値を含むパラメータは、ノイズパラメータに限定されない。例えば、再構成フィルタは、システムマトリクスに基づいて整形されても良い。より詳細には、システムマトリクスは、投影データの不均一サンプリングを指定する。従来のＦＢＰアルゴリズムはオブジェクトが均一にサンプリングされると仮定するが、可変のサンプリング戦略を使用することも可能である。例えば、重要な構造を含むある角度範囲内の信号は、第１の角度間隔でサンプリングされる。対照的に、当該領域外の別の角度範囲内の信号は、第１の角度間隔よりも大きい第２の角度間隔でサンプリングされる。
　なおステップＳ１００は、投影データの所定のデータ部分の各々に関する重み値を決定する際に、ステップＳ１２０の前に実行または完了されなくても良い。すなわち、再構成フィルタが整形される度にオプションで重み値が決定されても良い。
　さらに図３に示すように、ステップＳ１４０において画像発生部１１４Ｃは、ステップＳ１２０において整形された再構成フィルタを適用しながら、サブステップの所定のセットを実行する。１つの例示的なプロセスにおいては、さらに詳細に説明されるように、所定のＦＢＰ技法においてフーリエ変換の後に、整形されたフィルタを周波数領域内で投影データに適用することによって画像が再構成される。別の例示的なプロセスにおいては、所定のＦＢＰ技法において畳み込みとしてカスタマイズされたフィルタを空間領域内で投影データに適用することにより画像が再構成される。
　図４は、本実施形態に係る画像再構成方法におけるビュー毎のノイズ重み値の決定方式を示す図である。ビュー毎の重み付け方式において各ビューは、概して、特定のビューにおける最大ノイズ分散または平均ノイズ分散など、所定の関数に従って、重み値が割り当てられる。図４は、３つの位置、すなわち、位置Ａ、Ｂ、及びＣにおけるＸ線源、ならびに位置Ａ’、Ｂ’、及びＣ’における３つの対応する位置におけるＸ線検出器を例示する。例えば、Ｘ線源が位置Ａにあるとき、Ｘ線検出器は、対応する位置Ａ’にある。位置Ａ’においてＸ線検出器により収集された投影データはビューからなり、重み値ｗ（θ1）は、当該ビューに関する所定の重み付け方式に従ってビュー角度θ1に応じて決定される。上記のように、ビュー単位の重み値の決定は、１つの例示的なプロセスでは、再構成フィルタの整形に先立って実行および完了される。このビュー単位の重み値の決定は、他の例示的なプロセスでは、再構成フィルタが整形される度に決定される。
　本実施形態において投影データに含まれるノイズ分散は、ビュー単位の重み付け方式を使用してモデル形成可能である。一般に、ノイズの影響を補償する間に画像を再構成するための典型的な手法は、（１）式により規定される、ノイズ分散に応じて重み付けされた目的関数を最小化することである。
Figure JPOXMLDOC01-appb-M000003
　（１）式中、Ａは投影行列であり、Ｘは列ベクトルとして記述される画像配列であり、Ｐは列ベクトルとして記述される投影配列であり、Ｗは重み付けノルムを定義するノイズ重み付け行列である。この目的関数を最小化するために、（２）式に規定される反復ランドウェーバー（Landweber）アルゴリズムが使用される。
Figure JPOXMLDOC01-appb-M000004
　（２）式中、ＡTは逆投影行列であり、Ｘ(k)は第ｋ番目の反復において推定された画像であり、αはステップサイズであり、α＞０である。この再帰的表現は、（３）式のように非再帰的表現に書き直すことができる。
Figure JPOXMLDOC01-appb-M000005
　（ＡTＷＡ）-1＝Ａ-1Ｗ-1（ＡT）-1が存在すると仮定すると、（３）式は、（４）式のように簡素化される。
Figure JPOXMLDOC01-appb-M000006
　（４）式は、ノイズ重み付けＦＢＰアルゴリズムに相当し、周波数領域内における（４）式のフィルタ関数は、（５）式のように表される。
Figure JPOXMLDOC01-appb-M000007
　（５）式中、ωは周波数分散であり、ｗ（ｖｉｅｗ）はビュー単位で決定される重み値である。すなわち、ｗ（ｖｉｅｗ）は、特定のビュー角度におけるノイズ関連の重み関数である。αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、ｋは反復アルゴリズムにおいて反復数をエミュレートするパラメータである。ｖＦＢＰとして示されるビュー毎の重み付けＦＢＰアルゴリズムは、シフト不変の点応答関数（ＰＲＦ）を有する。
　上記の実施形態において再構成フィルタは、所定の反復再構成アルゴリズムに従って、所定数ｋ回の反復において特定の反復結果をエミュレートするように整形される。再構成フィルタにおいて反復をエミュレートするための上記の特徴は、ビュー毎の重み付けＦＢＰアルゴリズムを使用して投影データから画像を再構成するための他の実施形態に含まれるとは限らない。
　図５は、本実施形態に係る画像再構成方法における投影線毎のノイズ重み値の決定方式を示す図である。投影線毎の重み付け方式では、特定のビュー内の最大ノイズ分散または平均ノイズ分散など、所定の関数に従って、各投影線に重み値が割り当てられる。図５は、所定の位置におけるＸ線源、ならびに、対応する位置におけるＸ線検出器を例示する。投影線がＸ線源から投影されるとき、投影線は対応する位置におけるＸ線検出器に達する。投影データは、θにおけるビュー角度内の複数の投影線からなる。例えば、重み値ｗ（θ，ｔ1）は、ビュー角度θに応じ、かつ当該投影線に関する所定の重み付け方式に従って、特定の投影線ｔ1について決定される。上記の通り、投影線の重み値の決定は、１実施形態において、再構成フィルタを整形するのに先立って実行および完了される。投影線重み値は、他の実施形態において、再構成フィルタとして決定される。
　シフト不変ＰＲＦに関する要件が除去されると、投影線毎の重み付け方式が構築され、その周波数領域における窓化されたランプフィルタが上記の（５）式を変形することにより、（６）式のように記述される。
Figure JPOXMLDOC01-appb-M000008
　（６）式中、重み値、すなわち、重み係数ｗ（ｒａｙ）は、投影線で決定される重み計数である。ｗ（ｒａｙ）は、関連する投影線のノイズ分散に応じて決定される。ｋは、所定の反復再構成アルゴリズムに従って、特定の反復結果をエミュレートするための所定数の反復を示すパラメータである。Ｈ（ω，ｒａｙ）の逆フーリエ変換がｈ（ｔ，ｒａｙ）であるとすると、ｈ（ｔ，ｒａｙ）はフィルタの空間領域におけるカーネルである。ｐ（ｒ，θ）をビューθ及び投影線ｒにおける投影であるとし、ｑ（ｒ，θ）をフィルタリング後の投影であるとする。この場合、ｑ（ｒ，θ）は、以下の（７）式のような積分で規定される。
Figure JPOXMLDOC01-appb-M000009
　カーネルｈ（ｔ，ｒ）はｒに依存するため、この積分は畳み込み（コンボルーション）ではない。フィルタリングが（７）式のように空間領域で実施される場合、計算コストは畳み込みと同等である。この新しいｒＦＢＰアルゴリズムにおける逆投影ステップは、従来のＦＢＰアルゴリズムの逆投影ステップと同じである。従ってｒＦＢＰアルゴリズムは、計算的に効率的である。他方で、画像領域フィルタリングは、投影線毎のノイズ重み付けＦＢＰ－ＭＡＰ（ｒＦＢＰ－ＭＡＰ）アルゴリズムに関する実装形態で後に説明されるように、伝達関数を用いて２次元フーリエ領域内で実施される。
　重み係数を割り当てるための典型的な手法は、ｗ（ｒａｙ）を投影線計測のノイズ分散の逆数にすることである。この手法は、最適化問題に関する目的関数として、尤度関数（すなわち、結合確率密度関数）を使用して補われる。ノイズを含む計測に比してノイズを含まない計測を信頼すべきであるという根本原理を採用する場合、ノイズがより少ない計測により大きなｗ（ｒａｙ）が割り当てられ、ノイズがよい多い計測により小さなｗ（ｒａｙ）が割り当てられという制限下において任意に重み係数を割り当てることができる。
　本実施形態に係るｒＦＢＰアルゴリズムを実施するための最も容易かつより効率的な様式は、空間領域において投影をフィルタリングするために、（７）式を使用することである。現在、本発明者らは積分カーネルｈ（ｔ，ｒ）に関する解析式を有していない。各ビューθに関してｑ（ｒ，θ）を計算するための代替様式は、周波数領域において、以下のように方程式（７）を実施することである。
　ステップ１：ｒに関するｐ（ｒ，θ）の１次元フーリエ変換を見出して、Ｐ（ω，θ）を取得する。
　ステップ２．ａ：各投影線「ｒａｙ」に関して、重み係数ｗ（ｒａｙ）を割り当て、方程式（６）で表されるように周波数領域伝達関数を形成する。実装形態において、ωは周波数指数であり、整数をとる。
　ステップ２．ｂ：積、Ｑray（ω，θ）＝Ｐ（ω，θ）Ｈ（ω，ｒａｙ）を見出す。
　ステップ２．ｃ：ωに関するＱray（ω，θ）の１次元逆フーリエ変換を見出して、ｑ（ｒａｙ，θ）を取得する。
　上記の実施形態においてステップ１は、各ビューの列内の全てのｒを処理し、一方、ステップ２．ａから２．ｃは一度に１つの投影線のみを処理する。ビュー毎の重み付けＦＢＰ（ｖＦＢＰ）アルゴリズムでは、単一の重み値、すなわち、重み係数ｗ（ｖｉｅｗ）がビュー内の全ての投影線に割り当てられる。対照的に、投影線毎の重み付けＦＢＰ（ｒＦＢＰ）アルゴリズムでは、別個の重み値、すなわち、重み係数ｗ（ｒａｙ）がビュー内の投影線の各々に割り当てられる。結果として、ノイズ重み付け方式は、ビュー毎の重み付けＦＢＰ（ｖＦＢＰ）アルゴリズムよりも、投影線毎の重み付けＦＢＰ（ｒＦＢＰ）アルゴリズムにおいてより正確である。
　ビュー毎の重み付けＦＢＰ（ｖＦＢＰ）アルゴリズム及び投影線毎の重み付けＦＢＰ（ｒＦＢＰ）アルゴリズムに関して上で説明されたように、本実施形態は、画像を発生する前に、重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する。さらに、上記の実施形態では、再構成フィルタは所定の反復再構成アルゴリズムに従って、ｋで記された所定数の反復において特定の反復結果をエミュレートするように整形される。再構成フィルタにおいて反復をエミュレートするための上記の特徴は、ビュー毎の重み付けＦＢＰアルゴリズムまたは投影線ごとの重み付けＦＢＰアルゴリズムを使用して投影データから画像を再構成するための他の実施形態に含まれるとは限らない。
　その後、整形された再構成フィルタに基づいて、画像を発生するための工程が実行される。一般に、画像再構成におけるノイズ低減は、画像再構成、特に逆投影が加算処理であり、投影データの適切な重み付けが実質的に逆投影の分散を削減するという概念に基づいている。
　図６Ａは、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む画像再構成方法の典型的な流れを示す図である。ＦＢＰアルゴリズムが投影データに関するタスクを実行すると仮定すると、逆投影が実行される前に、再構成フィルタが投影データに適用される。ステップＳ１４０Ａに先立って再構成フィルタは、ノイズ重み値を含むパラメータに基づいて予め整形されている。ステップＳ１４０Ａにおいて整形された再構成フィルタが投影データに適用される。ステップＳ１４０Ａの前に投影データがフーリエ変換に供される場合、フィルタリングされた投影データは、ステップＳ１４０Ｂの前に逆フーリエ変換を受ける。ステップＳ１４０Ｂにおいてフィルタリングされた投影データは、逆投影される。逆投影により、フィルタリングされた投影データに基づいて画像が発生される。
　図６Ｂは、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む他の画像再構成方法の典型的な流れを示す図である。ＦＢＰアルゴリズムが投影データに関するタスクを実行すると仮定すると、逆投影が実行された後に、再構成フィルタが投影データに適用される。ステップＳ１４０Ｂの前に投影データがフーリエ変換に供される場合、投影データは、ステップＳ１４０Ｂの前に逆フーリエ変換を受ける。ステップＳ１４０Ｂにおいて投影データは、逆投影され、画像が発生される。整形された再構成フィルタは、逆投影されたデータに適用される。ステップＳ１４０Ａに先立って、再構成フィルタはノイズ重み値を含むパラメータに基づいて予め整形されている。ステップＳ１４０Ａにおいて、整形された再構成フィルタが投影データに適用される。
　図７は、ＦＢＰアルゴリズムでノイズ重み付けランプフィルタを使用した画像再構成方法の典型的な流れを示す図である。図７の処理においては、逆投影の前段においてノイズ重み付けフィルタリングがフーリエ領域で行わると仮定する。図７に示すように、投影データＰＤは、複数のビューにおいて収集される。投影データＰＤは、ビュー又は投影線単位で処理される。投影データＰＤは、変換されたデータがフーリエ領域または周波数領域に存在するように、フーリエ変換ＦＴに施される。ランプフィルタは、ノイズ重み値を含む所定パラメータのセットに従って整形される。これにより、重み付けランプフィルタＷＦが形成される。いくつかの実施形態では、複数のランプフィルタは、複数の重み付けランプフィルタＷＦになるように整形される。これらパラメータは、ノイズ重みパラメータに限定されず、他のパラメータを含んでも良い。この例示的な実施形態では、変換後の投影データは、所定のＦＢＰアルゴリズムにおいてノイズ重み付けランプフィルタＮＷＦに従ってフィルタリングされる。周波数領域内でフィルタリングされたデータは、投影データが逆投影ＢＰにおいて空間領域内で逆投影されるように、逆フーリエ変換ＩＦＴを受ける。あるいは、周波数領域内でノイズ重み付けランプフィルタＮＷＦを変換後のデータに適用する代わりに、別のノイズ重み付けランプフィルタＮＷＦ２が画像に適用されても良い。この場合、画像に適用されるノイズ重み付けランプフィルタＮＷＦ２は、変換後のデータに適用されるノイズ重み付けランプフィルタＮＷＦとは異なる。
　逆投影ＢＰは、画像発生の間、重み関数を使用して上述の可変角度サンプリングを正規化する。例えば、重要な構造を含む一定の角度範囲内の信号は、第１の角度間隔でサンプリングされる。当該領域外の他の角度範囲内の信号は、第１の角度間隔よりもより大きい第２の角度間隔でサンプリングされる。逆投影ＢＰは、サンプリング角度全体に亘る積分であり、角度サンプリング密度補償は、逆投影積分において本質的にヤコビアン係数である正規化係数を使用することによって達成される。数学的に、逆投影画像は、（８）式のように表現される。
Figure JPOXMLDOC01-appb-M000010
　（８）式中、初めに投影が逆投影される場合、ｐは角度θにおける投影である。あるいは、ＦＢＰアルゴリズムが使用される場合、ｐは、角度θにおいてフィルタリングされた投影であり、ｔは、検出器ビンの位置を示す。方程式（８）の離散的な表現は、（９）式のように表される。
Figure JPOXMLDOC01-appb-M000011
　（９）式中、ｎは検出器の位置指数であり、ｍは投影角度指数であり、Ｍは投影が収集されたビューの総数であり、「ＩＮＴ」は最近傍補間を示すために使用される。実際には、「ＩＮＴ」関数は使用せず、２つの近傍検出器ビンの間の線形補間が使用される。
　均一でない角度サンプリングに関して、（９）式で表されるような簡易的な逆投影が施される場合、当該逆投影は、単位角度当たりのビュー数である密度関数ｄ（０）に依存する。例えば、０＜θ＜π／４に関して、サンプリング間隔が１°であり、π／４＜θ＜πに関して３°である場合、密度関数は、（１０）式により表現され、逆投影（９）式は、（１１）式のように変更される。
Figure JPOXMLDOC01-appb-M000012
Figure JPOXMLDOC01-appb-M000013
　この場合、サンプリング密度関数は、実際の角度θの代わりに、角度指数ｉｎの関数である。
　図８は、本実施形態に係る、所定のエッジ保存ＦＢＰ－ＭＡＰアルゴリズムでノイズ重み付けランプフィルタの組を使用した画像再構成方法の典型的な流れを示す図である。図８の画像再構成方法は、投影データの組を逆投影して画像を再構成する前段に、フーリエ領域においてノイズ重み付けフィルタリングが実施される。最終的に、２つの画像に基づいて最終画像が発生される。
　図８に示すように、投影データＰＤは、複数のビューにおいて収集される。投影データＰＤは、ビュー単位で処理されても良いし、投影線単位で処理されても良い。投影データＰＤは、変換後のデータがフーリエ領域または周波数領域にあるようにフーリエ変換ＦＴを受ける。二つのランプフィルタは、ノイズ重み値を含む所定のパラメータのセットに従って個別に整形され、重み付けランプフィルタＷＦ－β0及びＷＦ－βに整形される。ランプフィルタは、正規化によってさらに整形可能である。このようにフィルタは、値βの関数である。値βが大きいほど、画像の平滑化がより強く行われる。例えば、平滑化が行われない場合あるいはボケの無い場合、β値はゼロなるように選択される。変換後の投影データは、所定のＦＢＰアルゴリズムにおいてノイズ重み付けされると良い。ノイズ重みづけされた投影データは、正規化されたランプフィルタＷＦ－β0及びＷＦ－βによってそれぞれフィルタリングされる。その後、周波数領域においてフィルタリングされた投影データは、逆フーリエ変換ＩＦＴに施される。逆フーリエ変換後の投影データは、２つの画像Ｘ0及びＸβが発生されるように空間領域内でそれぞれ逆投影される。平滑化を用いないで逆投影された、あるいは、最小限の平滑化を用いて逆投影された画像Ｘ0は、予め用意されたエッジマップＥＭＡＰ内のエッジマップ情報Ｂを抽出するために使用される。最小限に平滑化された画像Ｘ0、フィルタリングされた画像Ｘβ、及び抽出されたエッジ情報に基づいて一定の規則に従って最終画像ＦＩＭＧが発生される。
　図８に例示される実施形態においてエッジは、下記の（１２）式に示す目的関数においてベイジアン正規化項を使用してノイズ除去の間に保存される。（１２）式のペナルティ関数ｇ（Ｘ）は、注目画素近傍における画素値の段差を計測する。
Figure JPOXMLDOC01-appb-M000014
　（１２）式中、Ａは投影行列であり、Ｘは列ベクトルとして書き込まれた画像配列であり、Ｐは列ベクトルとして書き込まれた投影配列であり、Ｗは各投影に関する重み関数からなる対角行列であり、ｂは、（１２）式の第１項である忠実度項
Figure JPOXMLDOC01-appb-M000015
に対するベイジアン項ｇ（Ｘ）の重要性を調整する相対的な重み係数である。
　ペナルティｇ（Ｘ）は、画素値段差の二次関数である。当該画素値段差は、当該画素の画素値とその近傍の画素の平均値との間の差として計算される。大きな画素値段差は大きなペナルティに対応する。当該ペナルティ関数は平滑化画像をもたらす。
　エッジを過度に平滑化せずに画像ノイズを抑えるために、ペナルティ関数ｇ（Ｘ）は、エッジ及び非エッジ領域に関して、異なる特性を有すべきである。そのような能力を有するペナルティ関数として、閾値を有するフーバ関数が用いられる。画素値段差が閾値より小さく、画素値段差がノイズとして分類される場合、フーバ関数は二次式であり、平滑化を施す。画素値段差が閾値より大きく、画素値段差がエッジとして分類される場合、フーバ関数は線形であり、より小さい平滑化を施す。
　部分毎の一定画像を奨励する別の一般的な事前分布は、ＴＶ（全変動）計測である。離散的な一次元関数ｇ（Ｘ）の場合、ｇ（Ｘ）のＴＶ計測の勾配は３つの値だけをとる。すなわち、ｇ（Ｘ）がその左右の近傍より大きい場合、正の値、ｇ（Ｘ）がその左右の近傍よりも小さい場合、負の値、ｇ（Ｘ）がその左右の近傍の間である場合、０である。勾配タイプの最適化アルゴリズムでは、正の勾配値は、ｇ（Ｘ）の値を下方に押さえ込み、負の勾配値はｇ（Ｘ）の値を上方に押し上げ、０の勾配値は何の変更も行われていないことを意味する。従ってＴＶ事前分布は、単調関数を奨励し、振幅を抑え、鮮明なエッジを保存する。
　図８に例示されるように、実施形態に係る画像再構成方法は、投影データにエッジ保存ＦＢＰ－ＭＡＰアルゴリズムを適用する。βは、１つのノイズ重み付けランプフィルタについてのβ0、すなわち、ゼロに設定される。（１２）式を参照すると、β0ｇ（Ｘ）はゼロであり、何の平滑化も施されない。その後、何の平滑化も施されずに画像Ｘ0が再構成される。他方、他のパラメータセットの場合、βは、他のノイズ重み付けランプフィルタに関するβ、すなわり、非ゼロに設定される。βにより一定量の平滑化が施され、画像Ｘβが再構成される。これにより、２つの画像Ｘβ及びＸ0が再構成される。すなわち、第１の画像Ｘβは、非ゼロβ値を用いたｒＦＢＰ－ＭＡＰ再構成に基づいて再構成され、第２の画像Ｘ0はゼロβ値を用いたｒＦＢＰ再構成に基づいて再構成される。
　図９は、図８のエッジマップＢの発生処理の典型的な流れを示す図である。画像Ｘ0はゼロβ値を用いたｒＦＢＰ再構成法により再構成される。ステップＳ２００において、本質的に、画像勾配のノルムであるソーベルカーネルを使用して、エッジ検出によって画素値段差が評価される。ソーベルカーネルによって生成されたエッジマップは、次いで、ステップＳ２２０において、所定の閾値に従って、値０または値１を用いて二値画像に変換される。閾値を超える画像値は、値１に設定され、閾値未満の画像値は値０に設定される。その後、ステップＳ２４０において、所定のローパスフィルタによって二値画像をボケさせ、エッジマップＢを発生させる。二値画像をボケさせる理由は、二値画像に含まれるエッジ領域と非エッジ領域との境界を滑らかにするためである。最終画像は、下記の（１３）式に示すような、ＸbとＸ0との一次結合により発生される。
Figure JPOXMLDOC01-appb-M000016
　（１４）式中、「．*」は、ＭＡＴＬＡＢ（登録商標）で定義される点毎の乗算を表す。
　図１０は、事前画像を使用したＦＢＰ－ＭＡＰアルゴリズムに基づく画像再構成方法の典型的な流れを示す図である。一般に、所与の事前画像
Figure JPOXMLDOC01-appb-M000017
を加えることによって、ＦＢＰ－ＭＡＰアルゴリズムに基づいて画像を再構成するために、目的関数（８）は、下記の（１４）式により一般化される。
Figure JPOXMLDOC01-appb-M000018
　反復的な中間解は、（１５）式のように表される。
Figure JPOXMLDOC01-appb-M000019
　等価ＦＢＰ－ＭＡＰアルゴリズムは、下記の（１６）式と（１７）式とに表される２つのフーリエ領域フィルタ関数を有する。
Figure JPOXMLDOC01-appb-M000020
Figure JPOXMLDOC01-appb-M000021
　図１０に示すように、ＦＢＰ－ＭＡＰアルゴリズムは、投影データＰＤと事前画像ＰＩとの２つの入力を有する。最終画像ＦＩは、投影データＰＤに由来する第１の画像と、事前画像ＰＩに由来する第２の画像とに基づく加算画像である。ＦＢＰ－ＭＡＰ再構成は、（１６）式のような検出器フーリエ領域フィルタとその後に行われる逆投影とに基づいて第１の再構成画像を発生する。画像フーリエ領域処理ＩＦＤは、（１７）式のように画像フーリエ領域フィルタを使用して第２の画像を発生する。画像フーリエ領域フィルタは、２次元であっても良いし３次元であっても良い。
　図１１は、本実施形態に係るノイズ重み付け反復エミュレーションＦＢＰアルゴリズムを用いて画像再構成方法の典型的な流れを示す図である。一般に、再構成フィルタは、所定数の反復において特定の反復結果を模倣（エミュレート）する。再構成フィルタは、所定の反復再構成アルゴリズムに従って、重み値を含むパラメータに基づいて整形される。ノイズ重み付け反復エミュレーションＦＢＰアルゴリズム２００Ｄは、ステップ２０１において投影データが収集される。投影データは、撮影されているオブジェクトの線積分を表すラドン変換の形式で収集される。投影データは、図１に例示されるようなＸ線コンピュータ断層撮影装置やＰＥＴ装置、ＳＰＥＣＴ装置等の任意のモダリティにより収集される。
　ステップ２０２において反復数やノイズ重み値、及びシステムマトリクス等を含むパラメータの組合せに基づいて、少なくとも１つの再構成フィルタが整形される。これらのパラメータは、任意に選択可能である。例えば、パラメータｋは、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするために選択される。また、投影データに適用される再構成フィルタを整形するための反復パラメータｋと共に、予め決定された重み値Ｗが使用されても良い。その他のフィルタパラメータは、様々なアプリケーションに対処できるように、反復アルゴリズム内のステップサイズと、ベイジアン条件の相対的重み付けとを含む。
　ステップＳ２０２においてパラメータの組合せに従って再構成フィルタが整形された後、ステップ２０３において、整形された再構成フィルタを投影データに適用する。フィルタが施されたデータは、ステップ２０４において逆投影される。逆投影により画像が発生する。ステップＳ２０５においてエミュレーションを反復するか否かが決定される。エミュレーションを反復する場合、エミュレーションの反復のために再びステップＳ２０２に進む。ステップＳ２０５においてエミュレーションを反復しない場合、本実施形態に係るノイズ重み付け反復エミュレーションＦＢＰアルゴリズムは終了する。
　図１２は、本実施形態に係るノイズ重み付けＦＢＰアルゴリズムにおける上述のエミュレートされた反復の影響の概念を示す図である。従来のＣＴ撮影においてノイズ制御は、事前フィルタリングまたは事後フィルタリングを使用して実行されている。対照的に本実施形態においてノイズは、ノイズ重み付け及び反復エミュレートされたＦＢＰアルゴリズムを使用して画像再構成において実質的に低減される。画像再構成におけるノイズ低減は、画像再構成（特に、逆投影）が加算処理であり、投影に適切な重み付けを実行することにより逆投影の分散を削減できるという理念に基づいている。
　図１２に示すように、２つの投影、すなわち、計測データが２本の投影線Ｐ１及びＰ２によって表され、解、すなわち画像Ｉは２本の投影線Ｐ１及びＰ２の交点である。横軸Ｘ1及び縦軸Ｙ1は、画像Ｉにおける２つの変数を表す。１つの計測が画像から行われる場合、２つの画像画素は１本の直線に沿ったいずれかの場所の値をとることができる。２つの独立した計測が行われる場合、それら２本の直線の交点である、それら２つの画素の値が決定される。第３の計測が行われ、すべての計測がノイズを含む場合、それらの３本の直線は、一点において交差しない可能性がある。一般に、画像が多くの画素からなり、多くの矛盾するノイズを伴う計測が行われた場合、それらの計測は単一点で交わらない可能性があるため、撮像問題の解が存在しない可能性がある。
　数学では、ノイズ計測が矛盾する場合、それらの計測は「投影演算子の範囲内にない」と見なされる。他方で、計測が矛盾しない場合、それらの計測は投影演算子の範囲内にマッピングされる。実際には、投影演算子の範囲内のデータ成分だけが画像領域内に逆投影されるため、このマッピングはＦＢＰアルゴリズムにおいて逆投影によって保証される。すなわち、この投影の矛盾はＦＢＰ再構成画像に寄与しない。一般性を失わずに、直線によって示されるように、多くの画素を有する画像の場合ですら、計測は矛盾しないものの、依然として雑音を伴い、投影演算子の範囲内にある状況を仮定することができる。反復アルゴリズムとは異なり、ＦＢＰアルゴリズムは、ノイズを伴う投影のセットから、ノイズを伴う単一の解を提供することが可能である。
　図１２に示すように、ノイズがより少ない画像を発生するため、反復アルゴリズムの複数の反復ステップをエミュレートするようにＦＢＰアルゴリズムは修正される。解軌道はｋ＝０の初期値からｋ＝∞の最終的な解に向けて示される。ＦＢＰアルゴリズムのその後の反復は、結果として生じる、より少ないノイズを伴って画像を移動させる。反復アルゴリズムは、通常、真の解に達成する前に終了する。これは、計測データが低い光子計数によるノイズを含むため、真の解決策がノイズを伴うものである場合に特に当てはまる。
　ｋ＝０からの１つの解軌道は、真の解決策、すなわち、画像Ｉに向けた直線であり、この直線は、エミュレートされた反復がノイズ分散によって重み付けされないことを示す。直線上の点ＮＮＷ１及びＮＮＷ２は、それぞれ、第１の所定数の反復および第２の所定数の反復の後の中間結果を示す。この例では、第２の所定数の反復は、第１の所定数の反復よりも大きく、中間結果ＮＮＷ２は、第２の所定数の反復の後にエミュレートされる。他方で、中間結果ＮＷ１、ＮＷ２、及びＮＷ３は曲線上に位置し、この曲線は、それらのエミュレートされた反復がノイズ分散によって重み付けされていることを示す。中間結果ＮＷ１、ＮＷ２、及びＮＷ３は、それぞれ、第３の所定数の反復、第４の所定数の反復、及び第５の所定数の反復の後である。この例では、第５の数は、第３の数の反復及び第４の数の反復より大きく、第３の数は最小である。中間結果ＮＷ３は、第５の反復の後にエミュレートされる。この解軌道は、重み係数のセットによって決定される。１つの重み係数がそれぞれの計測に割り当てられる。より大きな係数は、関連する計測がより信頼できるものであることを意味する。一般的な手法は、計測のノイズ分散の逆数に比例する重み係数を割り当てることである。
　実際には、真のノイズを伴う解の近傍にあるノイズがより少ない疑似解が受け入れられる。疑似解の手順は、ノイズ分散重み値と反復指数ｋとを含むパラメータの所定のセットによってパラメータ化される。一般に、より小さいｋを伴う疑似解は、より平滑な画像に対応し、より大きなｋを伴う疑似解は、より鮮明であるが、より多くのノイズを伴う画像に対応する。
　以下、本実施形態に係る上記の特徴の特定の組合せに関する追加の実施例が記述される。一実施例は、ｒＦＢＰ－ＭＡＰによって指定される、最大事後確率（ＭＡＰ）アルゴリズムを用いた、投影線毎のノイズ重み付けＦＢＰが利用される。ｒＦＢＰ－ＭＡＰアルゴリズムに関する導出ステップは、ビューごとの重み付けＦＢＰ－ＭＡＰアルゴリズムに関する導出ステップと非常に類似している。幾つかの主な導出ステップについて以下に簡単に説明する。目的関数は、（１２）式と同じである。ベイジアン条件ｇ（Ｘ）が二次式である場合、ｇ（Ｘ）の勾配は、幾つかの行列Ｒに関してＲＸの形であり、解は、（１８）式により表されるランドウェーバーアルゴリズムによって得られる。
Figure JPOXMLDOC01-appb-M000022
　（１８）式に示される反復式は、等価の非反復式を有する。初期状態がゼロであるとき、非再帰的数式は、（１９）式により表される。
Figure JPOXMLDOC01-appb-M000023
　（ＡTＷＡ＋βＲ）-1が存在する場合、（１９）式中に示される総和は、（２０）式により表される閉形式を有する。
Figure JPOXMLDOC01-appb-M000024
　（ＡTＷＡ）-1が存在する場合、（ＡTＷＡ＋βＲ）-1＝［Ｉ＋β（ＡTＷＡ）-1Ｒ］-1（ＡTＷＡ）-1である。十分なサンプリングが逆ラドン行列Ａ-1の存在をもたらすことが仮定される場合、（ＡTＷＡ）-1＝Ａ-1Ｗ-1（ＡT）-1及び（ＡTＷＡ）-1ＡTＷ＝Ａ-1となる。従って（２０）式は、（２１）式のように簡略化可能である。
Figure JPOXMLDOC01-appb-M000025
　不十分なビュー数等の理由でＡ-1が存在しない場合、ＦＢＰアルゴリズムは、投影演算子Ａの疑似逆数であり、（ＡTＷＡ）+ＡTＷ＝Ａ+であることが確認でき、式中、「＋」は疑似逆数を示し、Ａ+はランプフィルタリングの後に逆投影ＡTが続くことである。従って（２１）式は、依然として、疑似逆数「＋」によって置換することによって逆数「－１」を可とする。以下、この逆数及び疑似逆数を識別しないこととする。「Ａ-1Ｐ」の表記は、従来のＦＢＰ再構成を表すために使用されることになる。
　（２１）式において、［Ｉ－（Ｉ－αＡTＷＡ－αβＲ）k］は、反復アルゴリズムの第ｋ番目の反復におけるウィンドウ効果を表す。ｋが無限大になる傾向があるため［Ｉ－（Ｉ－αＡTＷＡ－αβＲ）k］は、識別行列になる。［Ｉ＋β（ＡTＷＡ）-1Ｒ］-1は、ベイジアン正規化の効果を表す。β＝０のとき、ベイジアン修正は効果的ではない。（２１）式は、反復アルゴリズムがノイズ重み付け及びベイジアン正規化をどのように処理するかの見識と理解を提供し、（２１）式は、同じタスクを実行するためのＦＢＰタイプのアルゴリズムの構築をもたらし得る。
　最後に、ＦＢＰタイプのアルゴリズムが取得されるように、フーリエ領域表記が行列式から導出される。行列Ａ-1は、１次元ランプフィルタリングの後に逆投影が続くか、又は逆投影の後に２次元ランプフィルタリングが続くように処理される。行列式（２１）が「逆投影の後にフィルタリングする」アルゴリズムと見なされるとき、（２１）式における画像領域フィルタリングは、（２２）式に示す伝達関数を用いて２次元フーリエ領域において実行される。
Figure JPOXMLDOC01-appb-M000026
　（２２）式中、ｖx及びｖyは、それぞれ、ｘ及びｙに関する周波数である。
Figure JPOXMLDOC01-appb-M000027
は２次元周波数ベクトルである。
Figure JPOXMLDOC01-appb-M000028
　よく知られているように、投影演算子Ａが２次元空間内の線積分（すなわち、ラドン変換）であり、ＡTが演算子（すなわち、逆投影変換）である場合、投影と逆投影との組み合わされた演算子ＡTＡは、元の画像の２次元カーネル１／ｒとの２次元畳み込みである。この場合、ｘ－ｙデカルト座標において
Figure JPOXMLDOC01-appb-M000029
である。１／ｒの２次元フーリエ変換は、
Figure JPOXMLDOC01-appb-M000030
である。（２２）式においては、Ｒのフーリエ領域等価（Fourier domain equivalence）はＬ2Dによって示される。
　中央断面定理を使用することによって、画像領域の２次元フィルタリングは、投影領域の１次元フィルタリングに変換される。２次元画像のフーリエ領域における伝達関数（（２２）式）は、下記の（２３）式のように１次元投影のフーリエ領域における伝達関数に等しい。
Figure JPOXMLDOC01-appb-M000031
　（２３）式中、ωは検出器に沿った変数に関する周波数である。ω＝０であるとき、
Figure JPOXMLDOC01-appb-M000032
を定義する。（２３）式において、Ｌ1Dは（２２）式のＬ2Dの中央断面である。ｒＦＢＰ－ＭＡＰルゴリズムでは、
Figure JPOXMLDOC01-appb-M000033
は、１次元フーリエ領域で表現される修正ランプフィルタである。反復効果は、
Figure JPOXMLDOC01-appb-M000034
によって特徴付けられる。ベイジアン効果は、
Figure JPOXMLDOC01-appb-M000035
によって特徴付けられる。最終的な（ｋ＝∞）解は、β＝０の場合、ノイズ重み付けｗに依存しないが、β≠０の場合、ノイズ重み付けｗに依存する。ｒＦＢＰ－ＭＡＰアルゴリズムに関する逆投影について修正は必要とされない。
　（２２）式及び（２３）式のノイズに依存する重み係数ｗは、定数ではない。ビュー単にのノイズ重み付けの場合、ｗはビュー角度に応じ、ｗ＝ｗ（ｖｉｅｗ）である。投影線単位のノイズ重み付けの場合、ｗは投影線に応じ、ｗ＝ｗ（ｒａｙ）である。重み係数を割り当てるための一般的な手法は、ｗ（ｒａｙ）を投影線計測のノイズ分散の逆数に比例させることである。この手法は、最適化問題に関する目的関数として、尤度関数（すなわち、結合確率密度関数）を使用することによって補われる。根本原理は、ノイズのより大きい計測よりもノイズのより少ない計測を信頼すべきであるということである。ｋが無限大になる傾向があり、βがゼロになる傾向があるとき、修正ランプフィルタ（（２３）式）は従来のランプフィルタになり、ｒＦＢＰ－ＭＡＰアルゴリズムは、従来のＦＢＰアルゴリズムに低減される。
　ｒＦＢＰ－ＭＡＰアルゴリズムを実施する際、
Figure JPOXMLDOC01-appb-M000036
の逆フーリエ変換は、１次元修正ランプフィルタの空間領域カーネルである
Figure JPOXMLDOC01-appb-M000037
になると仮定される。ｒＦＢＰ－ＭＡＰアルゴリズムでは、ｗは投影線の関数であるため、ｗ＝ｗ（ｔ，θ）である。ｐ（ｔ，θ）をビューθ及び検出器上の位置ｔにおける投影にし、ｑ（ｔ，θ）をフィルタリングされた投影にする。この場合、ｑ（ｔ，θ）は、カーネル
Figure JPOXMLDOC01-appb-M000038
がθ及びｔに依存するため、畳み込みではない、以下の（２４）式に示す積分により規定される。
Figure JPOXMLDOC01-appb-M000039
　フィルタリング処理が（２４）式のように空間領域で実施される場合、計算コストは畳み込みと同じである。従って、ｒＦＢＰ－ＭＡＰアルゴリズムは計算的に効率的である。ｒＦＢＰ－ＭＡＰアルゴリズムを実施するための最も容易で最も効率的な方法は、空間領域で投影をフィルタリングするために（２４）式を使用することである。
　現在、積分カーネル
Figure JPOXMLDOC01-appb-M000040
に関する解析式は存在しない。ｑ（ｔ，θ）を計算する代替的な方法は、以下のように、周波数領域内の各ビュー角度θにおいて（２４）式を実施することである。
　ステップ１：ｔに関するｐ（ｔ，θ）の１次元フーリエ変換を見出して、Ｐ（ω，θ）を取得する。
　ステップ２：各投影線ｔについて重み係数ｗ（ｔ）を割り当て、周波数領域伝達関数（（２３）式）を形成する。ωは離散的な周波数指数であり、０、１、２、．．．などの整数値をとる。
　ステップ３：積
Figure JPOXMLDOC01-appb-M000041
を見出す。
　ステップ４：ωに関するＱ1（ω，θ）の１次元逆フーリエ変換を見出して、ｑ（ｔ，θ）を取得する。
　ステップ１は、各ビュー内の全てのｔを処理し、一方、ステップ２から４は、一度に１つのｔだけを処理する。
　上述のｒＦＢＰ－ＭＡＰアルゴリズムは、平行ビームを用いる撮像ジオメトリのみに限定されない。このｒＦＢＰ－ＭＡＰアルゴリズムは、ＦＢＰアルゴリズムが存在する限り、ファンビーム、コーンビーム、平面積分、及び減衰計測等の他の撮像ジオメトリに拡張可能である。ＦＢＰアルゴリズムに対する唯一の修正は、投影データフィルタリングである。各計測に重み係数が割り当てられ、周波数領域におけるフィルタ伝達関数は、撮像ジオメトリとは無関係な（２３）式として形成される。
　ｒＦＢＰ－ＭＡＰアルゴリズムの効果を例示するために、低線量使用のＸ線コンピュータ断層撮影装置を使用して、死体の胴体が走査された。画像は、ＦＢＰアルゴリズム、ならびにｒＦＢＰ－ＭＡＰアルゴリズムを用いて再構成された。撮像ジオメトリはコーンビームであり、Ｘ線源軌道は半径６００ｍｍの円である。Ｘ線検出器は６４列を有し、列の高さは０．５ｍｍであり、各列は８９６チャネルを有し、ファン角度は４９．２°であった。管電圧は１２０ｋＶであり、電流は６０ｍＡであった。３６０°に亘って１２００個のビューが均一にサンプリングされた。
　ＦＢＰアルゴリズムとしてフェルトカンプアルゴリズムを使用し、データは、まず、余弦関数を用いて重み付けされ、次いで、１次元ランプフィルタがコーンビーム投影の各列に適用された。最後に、３次元ボリュームデータを発生するためにコーンビーム逆投影が使用された。ｒＦＢＰ－ＭＡＰアルゴリズムを実施するために、１次元ランプフィルタは、（２３）式により表されるようにランプフィルタと置換された。反復指数ｋは、２０，０００と選択され、ステップサイズαは２．０に設定された。エッジ保存ＭＡＰ再構成の場合、２つのβ値、すなわち、β＝０．０００５及びβ＝０が使用された。ベイジアンペナルティ関数ｇ（Ｘ）は、横断断面内の中央画素値とその４近傍の平均値との間の差の共通二次関数であった。すなわち、Ｒは、その１次元バージョンが［－１　２　－１］のカーネルを用いた二次導関数である単なるラプラシアンである。ノイズ重み付け関数は、ｗ（ｔ，θ）＝ｅｘｐ（－ｐ（ｔ，θ））によって定義された。本実施形態においては、伝送計測が近似的にポワソン分布に従い、線積分ｐ（ｔ，θ）が伝送計測の対数であると仮定された。
　エッジマップは、２次元ソーベルカーネルを使用して断面毎に取得される。エッジ検出のための閾値は、最大画素値の７０％に設定された。次いで、３×３平滑カーネルを二値画像に適用する。平滑カーネルの適用により二値画像に含まれるエッジがボケ、二値画像から（１３）式に示すエッジマップＢが発生される。
　画像ボリュームデータは、５１２×５１２×５１２の３次元配列で再構成され、表示のために非中央軸断片が使用される。その他の事前フィルタリング又は事後フィルタリングはこの再構成に適用されなかった。
　図１３Ａは、従来のＦＢＰフェルトカンプアルゴリズムに基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。低線量のため、胴体を横切って左から右への著しいストリークアーチファクトが発生している。図１３Ｂは、β＝０を用いたｒＦＢＰ－ＭＡＰ再構成であるｒＦＢＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。ｒＦＢＰアルゴリズムにおける投影線単位のノイズ重み付けは、図１３Ａに出現したストリークアーチファクトを効果的に除去した。図１３Ｃは、二次平滑化事前分布を使用し、β＝０．０００５を用いたｒＦＢＰ－ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。図１３Ｃに示すように、ノイズ重み付けにより、ストリークアーチファクトは消えている。ボケの事前分布により、画像は特にエッジ領域で過度に平滑に見える。図１３Ｄは、実質的に図１３Ｂ及び図１３Ｃの重み付け加算である、エッジ保存事前分布を用いたｒＦＢＰ－ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。エッジ検出ソーベルカーネルを用いて、重み付け加算画像ＢがＸ0から抽出された。エッジマップは二値画像であった。エッジ領域から非エッジ領域への遷移が平滑であるように、エッジマップをボケさせるためにローパスフィルタが使用された。図１３Ｄの画像は、図１３Ｂよりもノイズがより少ないが、図１３Ｂの主なエッジを平滑でない状態で保持している。図１３Ｅは、エッジ保存ＦＢＰ－ＭＡＰ再構成において使用されるエッジマップＢを示す図である。
　計算効率を改善するために、多くのアプリケーションに関して、ビュー単位のノイズ重み付けｖＦＢＰ－ＭＡＰアルゴリズムを使用することが可能である。ノイズ重み付けは、各ビューにおける最大ノイズ分散によって決定される。図１３Ｆは、ｖＦＢＰ－ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。ｖＦＢＰ－ＭＡＰ再構成は、ノイズ重み付け関数ｗを除いて、全てのパラメータが図１３Ｄにおける再構成のパラメータと同一値に設定された。
　本実施形態に係る画像再構成においてノイズモデルは、ポワソン統計により規定される。あるいは、ノイズモデルは、ガウス分布により記述される電子ノイズを含む光子強度計測の複合ポワソン統計により規定される。確率論において複合ポワソン分布は、独立した一様に分布したランダム変数の和の確率分布であり、変数の数は、ポワソン分布に従う。ポワソン分布は、入射光子の数の記述に適している。蓄電時間の間に検出器に達する光子数Ｎは、Ｎ～ポワソン（λ）であり、この場合、λは計数率である。検出された各光子は、ランダムエネルギーＸiを有し、その分布はＸ線管の多色スペクトルと検出器効率とによって説明される。エネルギー積分検出器の場合、積分エネルギーは、Ｙ＝ΣNＸiによって与えられる。積分エネルギーＹは、（２５）式及び（２６）式のように複合ポワソン統計（平均値および分散）に従う。
Figure JPOXMLDOC01-appb-M000042
Figure JPOXMLDOC01-appb-M000043
　分散利得の概念は、Ｗによって示され、（２７）式のように規定される。
Figure JPOXMLDOC01-appb-M000044
　ポワソン分散に従うランダム変数の場合、Ｗ＝１である点に留意されたい。複合ポワソン変数の場合、Ｗは（２８）式に従うとする。
Figure JPOXMLDOC01-appb-M000045
　従ってノイズ分散利得は、Ｘ線検出器への入射Ｘ線のスペクトルの平均有効エネルギーに比例する。一般に、平均スペクトルエネルギーに影響を及ぼす以下の係数が存在する。
　１）管電圧（ｋＶｐ）。通常、管電圧は８０ｋＶｐから１４０ｋＶｐまで様々であり、スペクトルエネルギーの上限を決定する。
　２）ボウタイフィルタ及びスペクトルフィルタ、ヒール効果。これらは、ＣＴデータ内に常に存在する。ボウタイフィルタ及びヒール効果によりスペクトルは、各画素において異なる。従って各検出器画素に関して、独立して計算を行う必要がある。
　３）オブジェクトのサイズ及び構造。オブジェクトがビームを硬化させること、すなわち、Ｅ［Ｘ］を増大することはよく知られている。ビームハードニング効果は、経路長と減衰とに依存する。
　複合ポワソンモデルにおいてノイズ分散利得に関する平均スペクトルエネルギーの効果が利用される。収集データは、（２９）式のように表される。
Figure JPOXMLDOC01-appb-M000046
　（２９）式中、ｇは、検出器利得係数であり、ｅは分散Ｖeを伴う電子ノイズである。Ｖeを計測するために、Ｘ線管のスイッチをオフにしてデータを収集して、分散を計測する。平均及び分散は、各検出器画素に関して独立して時間方向に計測される。収集データの統計は、（３０）式により得られる。
Figure JPOXMLDOC01-appb-M000047
　すなわち、（３１）式が得られる。
Figure JPOXMLDOC01-appb-M000048
　オブジェクトに関する投影データを収集した後で、平均Ｅ［Ｚ］及び分散Ｖａｒ［Ｚ］を計測して、ｇＷを計算することができる。利得ｇは投影データから容易に取得できない点に留意されたい。未知のスケールｇの効果を補償するために、空気のスキャンデータ（ＷOBJ＝Ｗ／ＷAIR）により投影データを較正する。ノイズ分散に関する複合ポワソンモデルは、（３２）式により書き直される。
Figure JPOXMLDOC01-appb-M000049
　（３２）式中、Ｗは（２８）式により与えられる。
　本発明の構造および機能の詳細と共に、本発明の多数の特徴および利点が前述の説明に記載されているものの、本開示は、単なる例示であり、特に、部品の形状、サイズ、および構成、ならびにソフトウェア、ハードウェア、またはそれら両方の組合せの形の実装形態に関して詳細な変更を行うことは可能であるが、それらの変更は、添付の請求項が表現される条件の一般的な意味によって示される完全範囲まで本発明の原理内である点を理解されたい。
　本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
　１００…ガントリ、１０１…Ｘ線管、１０２…フレーム、１０３…Ｘ線検出器、１０４…データ収集回路、１０５…非接触データ伝送部、１０６…前処理部、１０７…回転部、１０８…スリップリング、１０９…高電圧発生部、１１０…システム制御部、１１２…記憶部、１１４…再構成部、１１４Ａ…ノイズ重み決定部、１１４Ｂ…フィルタ整形部、１１４Ｃ…画像発生部、１１５…入力部、１１６…表示部、１１７…電流調整部、２００…スキャン計画サポート装置

Claims (16)

  1. 　少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
    　前記重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する整形部と、
    　前記再構成フィルタに基づいて画像を発生する発生部と、
    　を具備するＸ線コンピュータ断層撮影装置。
  2. 　前記決定部は、ビュー毎に前記投影データに基づいて前記重み値を決定する、請求項１記載のＸ線コンピュータ断層撮影装置。
  3. 　前記再構成フィルタは、下記数式のＨ（ｗ，ｖｉｅｗ）により規定され、
    　ここで、ωは周波数分散であり、ｗ（ｖｉｅｗ）はビュー単位で決定される重み値であり、αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、ｋは反復アルゴリズムにおいて反復数をエミュレートするパラメータである、
    　請求項２記載のＸ線コンピュータ断層撮影装置。
  4. 　前記決定部は、投影線毎に前記投影データに基づいて前記重み値を決定する、請求項１記載のＸ線コンピュータ断層撮影装置。
  5. 　前記再構成フィルタは、下記数式のＨ（ｗ，ｒａｙ）により規定され、
    Figure JPOXMLDOC01-appb-M000002
     　ここで、ωは周波数分散であり、ｗ（ｒａｙ）は投影線単位で決定される重み値であり、αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、ｋは反復アルゴリズムにおいて反復数をエミュレートするパラメータである、
    　請求項４記載のＸ線コンピュータ断層撮影装置。
  6. 　前記発生部は、逆投影が実行される前に、前記再構成フィルタを前記投影データに適用する、請求項１記載のＸ線コンピュータ断層撮影装置。
  7. 　前記発生部は、逆投影が実行された後に、前記再構成フィルタを前記投影データに適用する、請求項１記載のＸ線コンピュータ断層撮影装置。
  8. 　前記ノイズモデルは、ポワソン統計を含む、請求項１記載のＸ線コンピュータ断層撮影装置。
  9. 　前記ノイズモデルは、ガウス分布に従う電子ノイズを含む光子強度計測の複合ポワソン統計によって記述される、請求項１記載のＸ線コンピュータ断層撮影装置。
  10. 　前記再構成フィルタを整形するために正規化する正規化部をさらに備える、請求項１記載のＸ線コンピュータ断層撮影装置。
  11. 　前記画像が基準画像に基づいてさらに発生される、請求項１０記載のＸ線コンピュータ断層撮影装置。
  12. 　前記整形部は、前記再構成フィルタとして第１の再構成フィルタと第２の再構成フィルタとを成形し、
    　前記第１の再構成フィルタは、エッジを保存するための第１の画像を発生するように整形され、第２の再構成フィルタは、ノイズを低減してエッジマップを発生するための第２の画像を発生するように整形され、
    　前記画像は、前記第１の画像と前記第２の画像と前記エッジマップとに基づいて発生される、
    　請求項１０記載のＸ線コンピュータ断層撮影装置。
  13. 　前記整形部は、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするように前記再構成フィルタを整形する、請求項１記載のＸ線コンピュータ断層撮影装置。
  14. 　前記成形部は、前記再構成フィルタをシステムマトリクスに基づいて整形する、請求項１記載のＸ線コンピュータ断層撮影装置。
  15. 　前記システムマトリクスが前記投影データの不均一サンプリングを指定する、請求項１４記載のＸ線コンピュータ断層撮影装置。
  16. 　少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
    　前記重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形し、
    　前記再構成フィルタに基づいて画像を発生する、
    　ことを具備する画像再構成方法。
PCT/JP2013/070659 2012-07-30 2013-07-30 Ｘ線コンピュータ断層撮影装置及び画像再構成方法 WO2014021349A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
EP13825580.7A EP2881039A4 (en) 2012-07-30 2013-07-30 IMAGE RECORDING DEVICE FOR X-RAY COMPUTER TOMOGRAPHY AND IMAGE RECONSTRUCTION PROCESS
KR1020137031515A KR101598265B1 (ko) 2012-07-30 2013-07-30 X선 컴퓨터 단층 촬영 장치 및 화상 재구성 방법
CN201380001535.4A CN103732147B (zh) 2012-07-30 2013-07-30 X射线计算机断层摄影装置以及图像重建方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US13/561,674 US10304217B2 (en) 2012-07-30 2012-07-30 Method and system for generating image using filtered backprojection with noise weighting and or prior in
US13/561,674 2012-07-30

Publications (1)

Publication Number Publication Date
WO2014021349A1 true WO2014021349A1 (ja) 2014-02-06

Family

ID=49994937

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2013/070659 WO2014021349A1 (ja) 2012-07-30 2013-07-30 Ｘ線コンピュータ断層撮影装置及び画像再構成方法

Country Status (6)

Country Link
US (1) US10304217B2 (ja)
EP (1) EP2881039A4 (ja)
JP (1) JP6301081B2 (ja)
KR (1) KR101598265B1 (ja)
CN (1) CN103732147B (ja)
WO (1) WO2014021349A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2018500082A (ja) * 2014-12-16 2018-01-11 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 対応確率マップ主導の視覚化
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
WO2022158575A1 (ja) 2021-01-21 2022-07-28 国立研究開発法人理化学研究所 断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体

Families Citing this family (29)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9208588B2 (en) * 2013-09-25 2015-12-08 Wisconsin Alumni Research Foundation Fast statistical imaging reconstruction via denoised ordered-subset statistically-penalized algebraic reconstruction technique
WO2015053787A1 (en) * 2013-10-11 2015-04-16 Analogic Corporation Tomosynthesis imaging
US10157483B2 (en) * 2013-12-19 2018-12-18 Board Of Regents, The University Of Texas System Backprojection approach for photoacoustic image reconstruction
JP2016540440A (ja) * 2014-06-12 2016-12-22 エスゼット ディージェイアイ テクノロジー カンパニー リミテッドＳｚ Ｄｊｉ Ｔｅｃｈｎｏｌｏｇｙ Ｃｏ．，Ｌｔｄ ピクチャー処理方法、装置
CN104318536B (zh) * 2014-10-21 2018-03-20 沈阳东软医疗系统有限公司 Ct图像的校正方法及装置
EP3286736B1 (en) 2015-02-03 2018-09-05 Koninklijke Philips N.V. Image reconstruction system, method, and computer program
JP6665119B2 (ja) * 2015-02-12 2020-03-13 株式会社日立製作所 Ｘ線ｃｔ装置及び画像処理装置
KR101812659B1 (ko) 2015-09-09 2017-12-27 삼성전자주식회사 단층 촬영 장치 및 그에 따른 단층 영상 복원 방법
EP3254623B1 (en) * 2016-06-09 2020-09-16 Agfa Nv Method and system for correcting geometric misalignment during image reconstruction in chest tomosynthesis
US10204425B2 (en) * 2016-07-29 2019-02-12 Toshiba Medical Systems Corporation Fast iterative reconstruction with one backprojection and no forward projection
US10198812B2 (en) * 2016-09-30 2019-02-05 Toshiba Medical Systems Corporation Data fidelity weight design for iterative reconstruction
JP6414236B2 (ja) * 2017-01-12 2018-10-31 オムロン株式会社 画像処理装置及び画像処理方法
US10706595B2 (en) * 2017-03-13 2020-07-07 General Electric Company System and method for reconstructing an object via tomography
WO2019012686A1 (ja) * 2017-07-14 2019-01-17 株式会社日立製作所 粒子線治療装置およびｄｒｒ画像作成方法
WO2019033390A1 (en) * 2017-08-18 2019-02-21 Shenzhen United Imaging Healthcare Co., Ltd. SYSTEM AND METHOD FOR IMAGE RECONSTRUCTION
CN107705261B (zh) * 2017-10-09 2020-03-17 东软医疗系统股份有限公司 一种图像重建方法和装置
CN107895387B (zh) * 2017-10-11 2020-10-02 四川大学 Mri图像重建方法及装置
US11210804B2 (en) * 2018-02-23 2021-12-28 Sony Group Corporation Methods, devices and computer program products for global bundle adjustment of 3D images
US10636128B2 (en) * 2018-06-18 2020-04-28 Analogic Corporation Recursive filter for space-variant images
US11049294B2 (en) * 2018-10-02 2021-06-29 Canon Medical Systems Corporation Activity-dependent, spatially-varying regularization parameter design for regularized image reconstruction
KR102297972B1 (ko) * 2018-12-28 2021-09-03 연세대학교 산학협력단 전역 변형 저감화 기법을 이용한 저선량 콘빔 전산화 단층 촬영 시스템
WO2020191435A1 (en) * 2019-03-22 2020-10-01 Southern Innovation International Pty Ltd Radiation detection with non-parametric decompounding of pulse pile-up
CN109949411B (zh) * 2019-03-22 2022-12-27 电子科技大学 一种基于三维加权滤波反投影和统计迭代的图像重建方法
US11176428B2 (en) * 2019-04-01 2021-11-16 Canon Medical Systems Corporation Apparatus and method for sinogram restoration in computed tomography (CT) using adaptive filtering with deep learning (DL)
US11763498B2 (en) * 2019-11-12 2023-09-19 Shanghai United Imaging Healthcare Co., Ltd. Systems and methods for image reconstruction
BR112022010084A2 (pt) * 2019-11-27 2022-08-30 Ka Imaging Inc Método e aparelho para ajuste espectral em geração de imagem de raio-x digital
CN112365479B (zh) * 2020-11-13 2023-07-25 上海联影医疗科技股份有限公司 Pet参数图像处理方法、装置、计算机设备及存储介质
CN112560243B (zh) * 2020-12-07 2022-11-15 桂林电子科技大学 一种改进频域临界采样图滤波器组的设计方法
CN113112561B (zh) * 2021-04-16 2021-12-07 赛诺威盛科技(北京)股份有限公司 图像重建方法、装置和电子设备

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012070814A (ja) * 2010-09-28 2012-04-12 Ge Medical Systems Global Technology Co Llc 画像処理装置およびプログラム並びにｘ線ｃｔ装置

Family Cites Families (38)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4991093A (en) * 1985-08-21 1991-02-05 Exxon Research And Engineering Company Method for producing tomographic images using direct Fourier inversion
US6438195B1 (en) * 2001-01-26 2002-08-20 Ge Medical Systems Global Technology Company, Llc Methods and apparatus for compensating for view aliasing artifacts
JP2004523037A (ja) * 2001-01-29 2004-07-29 ドイチェス クレブスフォルシュンクスツェントルム スチフトゥング デス エッフェントリヒェン レヒツ 収集された投影から空間体積の画像再構成をする方法および装置
US6870963B2 (en) * 2001-06-15 2005-03-22 Qualcomm, Inc. Configurable pattern optimizer
US7187794B2 (en) * 2001-10-18 2007-03-06 Research Foundation Of State University Of New York Noise treatment of low-dose computed tomography projections and images
US7206440B2 (en) * 2002-02-14 2007-04-17 Carnegie Mellon University Image smoothing with decoupled regularization
US7831296B2 (en) * 2002-11-27 2010-11-09 Hologic, Inc. X-ray mammography with tomosynthesis
US20040248066A1 (en) * 2003-06-03 2004-12-09 Recigno David T. Dental appliance on-line ordering including display of end product image and mold three-dimensional scanning for digital transmission
US7209535B2 (en) * 2003-06-20 2007-04-24 Wisconsin Alumni Research Foundation Fourier space tomographic image reconstruction method
US7653229B2 (en) * 2003-12-23 2010-01-26 General Electric Company Methods and apparatus for reconstruction of volume data from projection data
CN100563570C (zh) * 2004-07-07 2009-12-02 皇家飞利浦电子股份有限公司 心脏锥面光束ct重建中条纹伪影的减少
JP4535795B2 (ja) 2004-07-12 2010-09-01 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像処理装置及びｘ線ｃｔシステム
EP1828985A1 (en) * 2004-11-24 2007-09-05 Wisconsin Alumni Research Foundation Fan-beam and cone-beam image reconstruction using filtered backprojection of differentiated projection data
JP5248010B2 (ja) * 2006-02-17 2013-07-31 株式会社東芝 データ補正装置、データ補正方法、磁気共鳴イメージング装置およびｘ線ｃｔ装置
US7778386B2 (en) * 2006-08-28 2010-08-17 General Electric Company Methods for analytic reconstruction for mult-source inverse geometry CT
US20080118128A1 (en) * 2006-11-21 2008-05-22 Thomas Louis Toth Methods and systems for enhanced accuracy image noise addition
US7734076B2 (en) * 2006-12-11 2010-06-08 General Electric Company Material decomposition image noise reduction
KR100825048B1 (ko) 2006-12-22 2008-04-24 한국전기연구원 단층 영상을 고속으로 재구성하는 방법
US8374412B2 (en) * 2007-06-07 2013-02-12 Kabushiki Kaisha Toshiba Data processing apparatus, medical diagnostic apparatus, data processing method and medical diagnostic method
EP2216241B1 (en) * 2007-08-31 2016-04-06 Southwest University A cam self-adapting automatic speed-varying hub
DE102007061935A1 (de) * 2007-12-21 2009-06-25 Siemens Aktiengesellschaft Verfahren zur Qualitätssteigerung von computertomographischen Aufnahmeserien durch Bildverarbeitung und CT-System mit Recheneinheit
US8135186B2 (en) * 2008-01-25 2012-03-13 Purdue Research Foundation Method and system for image reconstruction
JP5335313B2 (ja) * 2008-08-05 2013-11-06 キヤノン株式会社 Ｘ線画像撮影装置、ｘ線画像撮影システム、ｘ線撮影制御装置、制御方法及びプログラム
CN101877124A (zh) 2009-04-28 2010-11-03 北京捷科惠康科技有限公司 一种医学图像的滤波方法及系统
EP2427868B1 (en) * 2009-05-07 2013-06-19 Koninklijke Philips Electronics N.V. System and method for generating a tomographic reconstruction filter
US8555616B2 (en) * 2009-07-09 2013-10-15 GM Global Technology Operations LLC Identifying ammonia non-slip conditions in a selective catalytic reduction application
DE102009039987A1 (de) * 2009-09-03 2011-03-17 Siemens Aktiengesellschaft Iterativer CT-Bildfilter zur Rauschreduktion
US8798353B2 (en) * 2009-09-08 2014-08-05 General Electric Company Apparatus and method for two-view tomosynthesis imaging
RU2557466C2 (ru) * 2009-11-03 2015-07-20 Конинклейке Филипс Электроникс Н.В. Устройство компьютерной томографии
US8731266B2 (en) * 2009-12-17 2014-05-20 General Electric Company Method and system for correcting artifacts in image reconstruction
CH702961A2 (de) * 2010-04-08 2011-10-14 Universitaetsklinik Fuer Nuklearmedizin Verfahren zum Abgleich von mit unterschiedlichen Systemen, z.B. Tomografen, aufgenommenen Bilddaten.
JP5379913B2 (ja) * 2010-05-14 2013-12-25 三菱電機株式会社 駐車支援装置、駐車支援システム、および駐車支援カメラユニット
CN103229212A (zh) 2010-11-30 2013-07-31 皇家飞利浦电子股份有限公司 具有基于恒定方差的加权因子的迭代重建算法
US8855394B2 (en) * 2011-07-01 2014-10-07 Carestream Health, Inc. Methods and apparatus for texture based filter fusion for CBCT system and cone-beam image reconstruction
GB201112359D0 (en) * 2011-07-19 2011-08-31 Univ Antwerpen Filter for tomographic reconstruction
JP5864958B2 (ja) * 2011-08-31 2016-02-17 キヤノン株式会社 画像処理装置、画像処理方法、プログラムおよびコンピュータ記録媒体
US8908942B2 (en) * 2011-10-21 2014-12-09 University Of Utah Research Foundation Filtered backprojection image reconstruction with characteristics of an iterative map algorithm
US20130202079A1 (en) * 2012-02-07 2013-08-08 Lifeng Yu System and Method for Controlling Radiation Dose for Radiological Applications

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012070814A (ja) * 2010-09-28 2012-04-12 Ge Medical Systems Global Technology Co Llc 画像処理装置およびプログラム並びにｘ線ｃｔ装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
AKIFUMI FUJITA: "Evaluation of Ultra-low-dose Thoracic Multi-detector-row CT using Different Reconstruction Kernels and New Reconstruction Algorithms", JAPAN RADIOLOGICAL SOCIETY ZASSHI, vol. 63, no. 9, 25 November 2003 (2003-11-25), pages 588 - 590, XP008173314 *
See also references of EP2881039A4 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2018500082A (ja) * 2014-12-16 2018-01-11 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 対応確率マップ主導の視覚化
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
WO2022158575A1 (ja) 2021-01-21 2022-07-28 国立研究開発法人理化学研究所 断層撮像のためのシステム、方法、プログラム、及びプログラムを記録した記録媒体

Also Published As

Publication number Publication date
KR101598265B1 (ko) 2016-02-26
EP2881039A1 (en) 2015-06-10
KR20140042809A (ko) 2014-04-07
US20140029819A1 (en) 2014-01-30
JP6301081B2 (ja) 2018-03-28
CN103732147A (zh) 2014-04-16
EP2881039A4 (en) 2016-08-03
CN103732147B (zh) 2016-10-26
JP2014023936A (ja) 2014-02-06
US10304217B2 (en) 2019-05-28

Similar Documents

Publication Publication Date Title
JP6301081B2 (ja) Ｘ線コンピュータ断層撮影装置、画像再構成方法及び再構成フィルタの構造
JP7139119B2 (ja) 医用画像生成装置及び医用画像生成方法
US10896486B2 (en) Denoising method and system for preserving clinically significant structures in reconstructed images using adaptively weighted anisotropic diffusion filter
US8731266B2 (en) Method and system for correcting artifacts in image reconstruction
La Riviere et al. Reduction of noise-induced streak artifacts in X-ray computed tomography through spline-based penalized-likelihood sinogram smoothing
Zhu et al. Improved compressed sensing‐based algorithm for sparse‐view CT image reconstruction
US10213176B2 (en) Apparatus and method for hybrid pre-log and post-log iterative image reconstruction for computed tomography
US9159122B2 (en) Image domain de-noising
US10111638B2 (en) Apparatus and method for registration and reprojection-based material decomposition for spectrally resolved computed tomography
CN109840927B (zh) 一种基于各向异性全变分的有限角度ct重建算法
US20100284596A1 (en) System and methods for fast implementation of equally-sloped tomography
US10204425B2 (en) Fast iterative reconstruction with one backprojection and no forward projection
Zhang et al. Statistical image reconstruction for low-dose CT using nonlocal means-based regularization. Part II: An adaptive approach
JP2007244871A (ja) サイノグラムから画像を再構成する方法およびこの方法を実行するコンピュータプログラム命令を格納するコンピュータ読み出し可能媒体およびサイノグラムから画像を再構成する装置
Levkovilz et al. The design and implementation of COSEN, an iterative algorithm for fully 3-D listmode data
US10475215B2 (en) CBCT image processing method
Karimi et al. A sinogram denoising algorithm for low-dose computed tomography
Xu et al. Statistical iterative reconstruction to improve image quality for digital breast tomosynthesis
CN110651302A (zh) 用于图像重建的方法和设备
US10198812B2 (en) Data fidelity weight design for iterative reconstruction
Bubba et al. Shearlet-based regularized reconstruction in region-of-interest computed tomography
Haase et al. Impact of the non‐negativity constraint in model‐based iterative reconstruction from CT data
KR102438068B1 (ko) 상호 정보량 기반 비국부 전역 변형 저감화 기법을 이용한 영상 화질 개선 방법 및 장치
US10515467B2 (en) Image reconstruction system, method, and computer program
Rajith et al. Edge preserved de-noising method for medical x-ray images using wavelet packet transformation

Legal Events

Date Code Title Description
ENP Entry into the national phase

Ref document number: 20137031515

Country of ref document: KR

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 2013825580

Country of ref document: EP

121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 13825580

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

NENP Non-entry into the national phase

Ref country code: JP