WO2014119664A1 - モデルに基づく反復的再構成のための逆投影と順投影との少なくとも一方におけるシステムオプティクス - Google Patents

モデルに基づく反復的再構成のための逆投影と順投影との少なくとも一方におけるシステムオプティクス Download PDF

Info

Publication number
WO2014119664A1
WO2014119664A1 PCT/JP2014/052115 JP2014052115W WO2014119664A1 WO 2014119664 A1 WO2014119664 A1 WO 2014119664A1 JP 2014052115 W JP2014052115 W JP 2014052115W WO 2014119664 A1 WO2014119664 A1 WO 2014119664A1
Authority
WO
WIPO (PCT)
Prior art keywords
model
system optical
image reconstruction
voxel
optical model
Prior art date
Application number
PCT/JP2014/052115
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 CN201480006906.2A priority Critical patent/CN104968275B/zh
Publication of WO2014119664A1 publication Critical patent/WO2014119664A1/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
    • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/416—Exact reconstruction
    • 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 of the present invention relate to iterative reconstruction based on system models, and more particularly to the use of system optical models in both backprojection and forward projection with a predetermined iterative reconstruction algorithm.
  • IR algorithms are divided into two categories.
  • the first category includes system optical models and is referred to as model-based IR (MBIR).
  • the second category does not include the system optical model.
  • the MBIR algorithm theoretically further improves the performance of the IR algorithm. The reason is that errors and statistics are substantially corrected. One source of error arises from system optics.
  • the MBIR algorithm further includes a system optical model (SOM).
  • SOM system optical model
  • An object of the embodiment is to provide an image reconstruction method and an image reconstruction system capable of reconstructing an image having excellent image quality.
  • An image reconstruction method includes a collection step of scanning a subject with a X-ray by a scanner and collecting raw data relating to the subject, an estimation step of estimating an image based on the raw data, and a system A forward projection step for forward projection to generate reprojection data using one of an optical model and a non-system optical model, and updating an image based on the reprojection data using the system optical model And an iterative process of repeatedly performing a backprojecting process of backprojecting to the above.
  • FIG. 11 shows a one-dimensional profile along the X axis through a 0.2 mm bead after 4000 iterations under the four conditions of FIG.
  • a full system optical model that includes optics in both forward projection and backprojection (FPBP) in one embodiment according to the present invention.
  • this SOM-FPBP is used in combination with a predetermined iterative (IR) algorithm to reconstruct an image.
  • the image data is reconstructed based on a full system optical model in both forward and back projection using an iterative algorithm (IR).
  • IR iterative algorithm
  • the reconstruction based on the full system optical model in both forward projection and backprojection using an iterative algorithm is expressed by the abbreviation IR-SOM-FPBP as necessary.
  • system optical model is a comprehensive term including individual system optical models.
  • the individual system optical models are divided into four groups including 1) voxel models, 2) x-ray source models, 3) detector element models, and 4) rotational blur or view integration models.
  • the voxel model further includes specific methods using multiple sub-voxels (microrays), multiple blobs, and multiple splines.
  • the X-ray source model further includes a specific method using a plurality of sub-sources (micro-rays), a plurality of low-pass filters, and a plurality of separable footprints.
  • the detector element model further includes specific methods using sub-detectors (microrays) and separable footprints.
  • the rotation blur or view integration model further includes specific methods using subviews (microrays) and low pass filters. Models other than the specific methods described above are considered non-system optical models. In other words, non-system optical models include methods such as single ray or distance-driven methods.
  • the distance drive method includes certain characteristics or aspects of the system geometry, but is not considered part of the system optical model in this embodiment.
  • FIG. 1 shows a system optical model of IR-SOM-FPBP according to the present embodiment.
  • the system optical model of FIG. 1 includes a plurality of detector element models 110 corresponding to the plurality of detector elements included in the detector 100 and an X-ray source model 310 corresponding to the X-ray source 300, respectively. Further, between the detector element model 110 and the X-ray source model 310, a plurality of voxel models 210 respectively corresponding to a plurality of voxels constituting the predetermined volume 200 are conceptually arranged.
  • a plurality of micro rays are incorporated in the system optical model.
  • the system optical model according to the present embodiment is further specified as a system optical model based on micro rays.
  • system optics are used to define a three-dimensional X-ray beam.
  • the micro-ray based system optical model takes into account the geometry for a given detector, the geometry for a given voxel, and the geometry for a given x-ray source in order to accurately define a 3D x-ray beam.
  • the three-dimensional X-ray beam B includes a predetermined number of individual microrays such as R1, R2, R3, and R4. If the 3D X-ray beam B is filled with a sufficient number of micro-rays, it samples the object inside the 3D X-ray beam more accurately than the distance driven method according to the prior art, resulting in an image with a higher resolution. be able to.
  • a three-dimensional X-ray beam is used in both the forward and back projection steps of a predetermined iterative reconstruction method. The three-dimensional X-ray beam can be applied to any iterative reconstruction algorithm according to this embodiment.
  • the X-ray source model and the detector model are related, and in the back projection process, the X-ray source model and the voxel model are related.
  • the detector element is conceptually divided into a plurality of micro-detector points
  • the voxel is conceptually divided into a plurality of micro-voxel points
  • the X-ray source is a plurality of micro-X Conceptually divided into source points.
  • a micro ray connects a micro detector point, a micro voxel point, and a micro X-ray source point. The micro ray is used for forward projection and back projection.
  • system optics is applied in both forward projection and backprojection. If necessary, system optics may be applied only to backprojection. When system optics is applied only to backprojection, X-ray source deconvolution is performed, providing better resolution than prior art methods where system optics is only included in forward projection. Resolution is improved in IR using full system optical model in backprojection in addition to forward projection using OS-SART algorithm.
  • the X-ray source model 310 will be described with reference to FIG. 2A.
  • the X-ray source model 310 and the detector element model are incorporated into the forward projection.
  • the geometry of the X-ray source model 310 includes an anode angle ⁇ anode , an X-ray source width Sx in the Xs direction, and an X-ray source height Sz in the Zs direction.
  • the surface 310A of the X-ray focal point 300 is sampled by N ⁇ S micro X-ray source points located at the centers of the respective grids.
  • Each of the grids in this embodiment is equally defined according to a predetermined dimension along the Xs axis and the Ys axis.
  • the N ⁇ S micro X-ray source points are a collection of a predetermined number of discrete individual micro X-ray source points labeled ⁇ SP1, ⁇ SP2,..., ⁇ SPN. If necessary, N ⁇ s micro X-ray source points are distributed in an arbitrary manner. In other words, the sampling pattern is not limited to a rectangular grid, and includes an arbitrary pattern such as a Gaussian distribution or a random pattern. Each micro X-ray source point from ⁇ SP1 to ⁇ SPN has an associated spatial X-ray source response weight value SWu.
  • the detector model 110 according to the present embodiment will be described with reference to FIG. 2B.
  • the detector model 110 and the X-ray source element model are incorporated into the forward projection.
  • the detector geometry of the curved detector 100 includes an angular channel width ⁇ and a segment width Wseg.
  • the surface 110A of the detector 100 is sampled by N ⁇ D micro-detector points located at the center of each grid. Each of the grids is equally defined according to a predetermined dimension along the ⁇ axis and the Z D axis.
  • the N ⁇ D microdetector points are a collection of a predetermined number of discrete individual microdetector points labeled ⁇ DP1, ⁇ DP2,..., ⁇ DPN.
  • the N ⁇ D microdetector points are distributed in any manner.
  • the sampling pattern is not limited to a rectangular grid, and includes an arbitrary pattern such as a Gaussian distribution or a random pattern.
  • each microdetector point from ⁇ DP1 to ⁇ DPN has an associated spatial X-ray source response weight value DWu.
  • each microdetector element from ⁇ DP1 to ⁇ DPN includes a dead space (pitch) DS to accommodate the inactive space between the individual microdetector elements.
  • dead space DS is between the difference between the active segment width Wseg AP along the segment width Wseg and Z D axis, the active angular channel width [Delta] [gamma] AP along a corner channel width [Delta] [gamma] and ⁇ -axis Defined by the difference between That is, the internal dead space DS is active aperture of the detector elements is defined by the active segment width Wseg AP and the angular channel width [Delta] [gamma] AP.
  • the voxel model 210 will be described with reference to FIG.
  • the voxel model 210 and the X-ray source element model are incorporated into the back projection.
  • the voxel geometry of the voxel 200 includes a voxel width VxWIDTH in the Xv direction, a voxel length VyLENGTH in the Yv direction, and a voxel height VzHEIGHT in the Zv direction.
  • the voxel 200 is sampled by N ⁇ v micro voxel points located at the intersection. Each of the intersection points is equally defined along the Xv axis, the Yz axis, and the Zv axis.
  • the N ⁇ v micro voxel points are a collection of a predetermined number of discrete individual micro voxel points labeled ⁇ VP1, ⁇ VP2,..., ⁇ VPN.
  • N ⁇ v micro voxel points are distributed in an arbitrary manner.
  • the sampling pattern is not limited to a rectangular grid, and includes an arbitrary pattern such as a Gaussian distribution or a random pattern.
  • Each micro voxel point from ⁇ VP1 to ⁇ VPN has an associated spatial X-ray source response weight value VWui.
  • the voxel geometry in FIG. 3 is a cube, but may be a voxel geometry of a different shape.
  • the full system optical model includes an X-ray source 300, a voxel 200A, and a detector 100.
  • the voxel 200A is divided into micro voxels having micro voxel elements with a total number of N ⁇ v .
  • the voxels are not limited to the cube shown in FIG.
  • the shape of the voxel according to the present embodiment is a cube and is sampled by a grid, but the voxel is not limited to a specific shape. Voxels that are three-dimensional cubes sampled by a given grid are inefficient due to the presence of microray overlap.
  • a voxel 200A shown in FIG. 4 has a cylindrical shape and has a rotating two-dimensional micro voxel plane 200B.
  • the voxel 200A having a cylindrical shape has a radius Rv and a height Hv.
  • the voxel 200A has four micro voxel points from ⁇ VP1 to ⁇ VP4, and reduces the number of micro voxels.
  • the rotating two-dimensional micro voxel plane 200B rotates in accordance with the viewing angle and the voxel position, and extends from the X-ray source center SC to the voxel center VC.
  • the rotating 2D micro voxel plane 200B thus reduces voxel sampling from 3D to 2D.
  • each microray inside the three-dimensional X-ray beam is defined by a corresponding pair of a start point and an end point.
  • the micro X-ray source point is the starting point for both forward and back projection.
  • the micro detector point is the end point for forward projection and the micro voxel point is the end point for back projection.
  • the microrays may be distributed in parallel or in a conical shape.
  • the microrays are preferably distributed in a conical shape in order to correctly sample the voxels.
  • the micro rays may be distributed in a conical shape or in a parallel manner.
  • the total number of microrays NuRays from ⁇ R1 to ⁇ R4 is equal to the total number of micro X-ray source points NuSrc from ⁇ SP1 to ⁇ SP4 and the total number NuDet of microdetector points from ⁇ DP1 to ⁇ DP4. That is, there is a one-to-one connection between the micro X-ray source point and the micro detector point applied to the forward projection.
  • the total number of microrays NuRays from ⁇ R1 to ⁇ RN is the product of the total number of micro X-ray source points NuSrc from ⁇ SP1 to ⁇ SP4 and the total number NuDet of microdetector points from ⁇ DP1 to ⁇ DP4. equal. That is, each of the micro X-ray source points is connected to all of the micro detector points applied to back projection.
  • the micro ray is dispersed in all the micro detector points from ⁇ DP1 to ⁇ DP4.
  • the micro ray bundle from ⁇ Rs2 to ⁇ Rs4 four microrays are connected from each of the micro X-ray source points from ⁇ SP2 to ⁇ SP4 to all the microdetector points from ⁇ DP1 to ⁇ DP4. Note that the total number of NuSrc and NuDet may be different from each other.
  • each of the plurality of micro rays u included in the three-dimensional X-ray beam is individually forward-projected.
  • the final forward projection value for a detector element value is the average of all NuRays forward projections contained in the 3D x-ray beam.
  • the dispersion of the three-dimensional X-ray beam is parallel or conical as required.
  • the forward projection ray tracing of the micro ray u is defined by the following equation (1) using the response weight value.
  • FP ′ u is the final forward projection value for the detector element value
  • SW u is the associated spatial X-ray source response weight value
  • VW u is the associated space.
  • VW is the typical voxel response weight value
  • DW u is the associated spatial detector response weight value
  • L nu is the path length of ray u through voxel nu
  • V nu intersects ray u It is a voxel.
  • the final forward projection value for detector elements ch, seg is the average of the individual forward projections, as defined by equation (2) below.
  • FP ch and seg are the final forward projection values for the detector elements specified by the channel ch and the segment seg, and NuRays is the total number of micro rays.
  • any ray tracing algorithm such as Siddon for calculating a line integral value can be used.
  • the body axis cross section is divided into four quadrants.
  • the image volume is upsampled by two factors before forward projection.
  • either the X or Y plane is used as an “anchor” plane to reduce the 3D configuration to a 2D configuration.
  • this graph shows that the microray given abscissa angle ⁇ is between 45 degrees and 135 degrees. Since the given body axis crossing angle ⁇ is 45 ° ⁇ ⁇ 135 °, the closest one is fixed to the Y plane according to the following plane.
  • the raytracing forward projection of the micro ray u may be defined by the following equation (4).
  • any micro X-ray source point is x us , yus , z us
  • any micro detector point is x ud , y ud , z ud
  • the corresponding body axis transverse angle is ⁇
  • micro ray u The length of the crossing point of the body axis is L ⁇ , and the voxel is constant.
  • the forward projection on the micro ray passing through the image volume from the micro X-ray source point us to the micro detector point ud is defined by the following (4).
  • n represents the voxel where the micro ray intersects
  • N is the total number of voxels where the micro ray intersects.
  • the overall forward projection value for detector element m is defined by the following equation (5).
  • N uS is the total number of micro X-ray source points (NSx ⁇ NSz)
  • N uD is the total number of micro-detector point (NuCh ⁇ NuSeg).
  • NSx is the number of micro X-ray source elements along the X direction of the X-ray source model
  • NSz is the number of micro X-ray source elements along the Z direction of the X-ray source model. It is shown.
  • FIG. 2B shows that NuCh is the number of micro-detector elements along the ⁇ direction of the detector model
  • NuSeg is the number of micro-detector elements along the Z direction of the detector model. .
  • microrays related to the X-ray source and the voxel in back projection will be described.
  • the X-ray source model is the same as in forward projection, and the detector model is replaced by a voxel model.
  • each microray u in the three-dimensional X-ray beam is backprojected individually.
  • a voxel is divided into a predetermined number of micro voxels or micro voxel points.
  • Each of the micro voxel points from voxel [1] to voxel [n] in the voxel model 210 receives NuRays microrays averaged over the x-ray source 300 and the voxel 200.
  • the detector 100 has a projection or reprojection value PD [ch u , seg u ] at each of the micro detector points.
  • the updated voxel value is an average value of back projection with NuRays microrays included in the three-dimensional X-ray beam.
  • back projection a conical beam is used.
  • the ray-tracing forward projection of microray u is defined by equation (6) with response weight values.
  • ch u is the channel position of micro ray u
  • seg u is the segment position
  • PD [ch u , seg u ] is the projection value or reprojection value for ch u and seg u
  • SW u is the associated spatial X-ray source response weight value
  • V W u is the associated spatial voxel response weight value
  • DW u is the associated spatial detector response weight value.
  • the nearest value is used, and in another embodiment, bilinear interpolation is used.
  • the voxel is a micro voxel or micro voxel point with a cylinder having a radius of R Voxel and a height of H Voxel and a total of NuV positions.
  • a two-dimensional backprojection plane having The plane 200B rotates to have a viewing angle and a voxel position and is always perpendicular to the center of the x-ray source to the voxel center ray.
  • N uS ⁇ N uV microrays are backprojected.
  • the back projection values of N uS ⁇ N uV microrays are averaged over the x-ray source and voxel as defined by the following equation (7).
  • BP n is the voxel value of the back-projected voxel n
  • m represents a detector element crossed by the microrays included in the three-dimensional X-ray beam
  • M is a detection detected by crossing the three-dimensional X-ray beam
  • a n, m is a back projection operator.
  • step S100 raw data is collected using a scanner using a predetermined modality such as a computed tomography (CT) apparatus.
  • step S110 based on the collected raw data, a seed image is estimated according to a predetermined algorithm. That is, in step S110, a seed image is reconstructed initially based on the raw data.
  • step S120 forward projection value data (OS-SART: ordered-subset simultaneous algebraic reconstruction technique) is used from the estimated seed image using a predetermined sequential reconstruction algorithm (such as ordered-subset simultaneous algebraic reconstruction technique). Reprojection value data).
  • OS-SART ordered-subset simultaneous algebraic reconstruction technique
  • the estimated seed image is forward projected using either a system optical model or a non-system optical model to generate reprojection value data.
  • step S130 the reprojection value data in the full system optical model reconstruction is backprojected to update the seed image with the system optical model.
  • the above-described forward projection and backprojection in steps S120 and S130 are repeatedly applied until a predetermined condition is satisfied. This condition is determined in step S140.
  • step S140 if a condition such as a predetermined number of iterations is satisfied, the process ends. On the other hand, if the condition is not met in step S140, the process returns to step S120.
  • the fan beam simulation includes a full three-dimensional geometry (ie, the fan beam has a thickness due to the x-ray source and the detector and is not simply a zero-dimensional two-dimensional plane). .
  • simulation (1) beads having a diameter of 1.0 mm and an X-ray source having a diameter of 6.0 mm are used.
  • simulation (2) 0.2 mm diameter beads and 1.1 mm X-ray source are used.
  • simulation (1) by incorporating an x-ray source that is much larger than the bead, including system optics compensates for x-ray source blurring and can recover objects that are much smaller than the x-ray source Evaluate whether or not.
  • the resolution in the transverse direction of the body axis is evaluated by expressing a realistic focus size using beads smaller than the size of the detector.
  • the beads were reconstructed using 4 reconstruction algorithms. That is, standard filtered back projection (FBPJ) using the LAKS kernel, IR using the full system optics described above (IR-SOM-FPBP), pencil beam IR without system optics (IR-P), and IR (IR-SOM-FP) with system optics only for forward projection.
  • FBPJ standard filtered back projection
  • IR-SOM-FPBP full system optics described above
  • IR-P pencil beam IR without system optics
  • IR-SOM-FP IR-SOM-FP
  • Table 1 summarizes the parameters used for the two simulations and the standard fan beam configuration.
  • FIG. 9 an iteratively reconstructed image containing full system optics is shown for comparison under certain conditions according to this embodiment.
  • the reconstructed image is reconstructed under the condition of a 6.0 mm X-ray source and 1.0 mm beads.
  • the combination of FBPJ and IR is used to reconstruct and normalize the image.
  • Three reconstructed images under the conditions of FBPJ, IR-P, and IR-SOM-FP are rendered with distorted beads because the exaggerated X-ray source size is not fully sampled. ing. Only images reconstructed under IR-SOM-FP conditions are the most sensitive and most distorted. Under the conditions of IR-SOM-FPBP, the beads are accurately reconstituted by fully recovering the bead size.
  • IR-SOM-FPBP image and plot show some small side lobes, which may be due to a mismatch between forward and back projection.
  • Images reconstructed under the conditions of IR-SOM-FPBP are filtered backprojection method (FBPJ), pencil beam IR method without system model (IR-P), IR using system model only for forward projection It has substantially better resolution than images reconstructed under legal conditions.
  • FIG. 10 is a diagram showing a one-dimensional profile along the x-axis passing through 1.0 mm beads after 4000 iterations under the four conditions of FIG.
  • the profiles under the conditions of IR-SOM-FPBP are filtered backprojection (FBPJ), pencil beam IR without system model (IR-P), and forward projection only. It has substantially better resolution than an image reconstructed under the conditions of the IR method using a system model.
  • IR-SOM-FPBP has the highest resolution and FWHM is 0.48 mm. This is almost half of FBPJ (0.87 mm), which is considerably better than IR-P (0.63 mm) and IR-SOM-FPBP (0.59 mm). Regarding the exaggerated X-ray source simulation, IR-SOM-FPBP shows side lobes.
  • FIG. 12 is a graph illustrating a one-dimensional profile along the x-axis passing through a 0.2 mm bead after 4000 iterations under the four conditions described above according to the present invention.
  • profiles under the conditions of IR-SOM-FPBP are filtered backprojection (FBPJ), pencil beam IR without system model (IR-P), and forward projection only. The resolution is substantially better than the reconstructed image under the conditions of the IR method using the system model.
  • the profile under IR-SOM-FPBP conditions after 10,000 iterations is plotted using a dashed line.
  • an image having excellent image quality can be reconstructed.
  • Detector 110 ... Detector element model, 200 ... Voxel, 210 ... Voxel model, 300 ... X-ray source, 310 ... X-ray source model, R1, R2, R3, R4 ... Microray

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Abstract

　優れた画質を有する画像を再構成すること。　収集工程Ｓ１００において、スキャナにより被検体をＸ線でスキャンして被検体に関する生データを収集する。推定工程Ｓ１１０において、生データに基づいて画像を推定する。順投影工程Ｓ１２０において、システム光学モデルと非システム光学モデルとのうちの一方を用いて再投影データを発生するために順投影する。逆投影Ｓ１３０において、システム光学モデルを用いて再投影データに基づいて画像を更新するために逆投影する。順投影工程Ｓ１３０と逆投影工程Ｓ１４０とは反復される。

Description

モデルに基づく反復的再構成のための逆投影と順投影との少なくとも一方におけるシステムオプティクス
　本発明の実施形態は、システムモデルに基づく反復的再構成に関し、より詳細には、所定の反復的再構成アルゴリズムでの逆投影と順投影との両方におけるシステム光学モデルの使用に関する。
　反復的再構成（ＩＲ：iterative reconstruction）アルゴリズムの標準的なフィルタ補正逆投影法（ＦＢＰ：filtered BackProjection）よりも優れている２つの長所は、解像度の改善とノイズ性能の向上とである。このため、ＩＲアルゴリズムによると、標準的なＦＢＰを用いる場合に従来要求されたものよりも低い患者被爆線量の使用が可能になる。
　ＩＲアルゴリズムは、２つのカテゴリに分けられる。第１のカテゴリはシステム光学モデルを含んでおり、モデルに基づくＩＲ（ＭＢＩＲ：model-based IR）と称される。これに対し、第２のカテゴリは、システム光学モデルを含まないものである。ＭＢＩＲアルゴリズムは、理論的には、ＩＲアルゴリズムの性能を更に改善する。その理由は、誤差と統計とが実質的に訂正されるからである。誤差の原因の１つには、システムオプティクスから生じるものがある。ＭＢＩＲアルゴリズムは、更に、システム光学モデル（ＳＯＭ: system optics model）を含む。このように、優れた解像度とノイズ性能とを有するＣＴ画像を、低被曝線量で収集されたデータから再構成するためには、ＭＢＩＲにおけるＳＯＭの利用を改善することが望まれる。
　実施形態の目的は、優れた画質を有する画像を再構成可能な画像再構成方法及び画像再構成システムを提供することにある。
　本実施形態に係る画像再構成方法は、スキャナにより被検体をＸ線でスキャンして前記被検体に関する生データを収集する収集工程と、前記生データに基づいて画像を推定する推定工程と、システム光学モデルと非システム光学モデルとのうちの一方を用いて再投影データを発生するために順投影する順投影工程と、前記システム光学モデルを用いて前記再投影データに基づいて画像を更新するために逆投影する逆投影工程とを反復的に実行する反復工程と、を具備する。
　優れた画質を有する画像を再構成すること。
本実施形態に係る順投影と逆投影との両方でオプティクスを含むフルシステム光学モデルを示す図 本実施形態に係るＸ線源とそれに対応するＸ線源モデルとを示す図 本実施形態に係るＸ線検出器とそれに対応するＸ線検出器モデルとを示す図 本実施形態に係るボクセルとそれに対応するボクセルモデルとを示す図 本実施形態に係る順投影と逆投影との両方を用いるフルシステム光学モデルに関するボクセルとそれに対応するボクセルモデルとを示す他の図 本実施形態に係る順投影と逆投影との両方を用いるフルシステム光学モデルにおける平行に分布するマイクロレイを示す図 本実施形態に係るシステム光学モデルにおける円錐状に分布するマイクロレイを示す図 本実施形態に係る３次元最隣接レイトレーシングアルゴリズムを用いた順投影を示す図 本実施形態に係る逆投影におけるＸ線源とボクセルとに関するマイクロレイを示す図 本実施形態に係る再構成プロセスの流れを示す図 本実施形態に係るフルシステムオプティクスを用いて反復的に再構成された画像を示す図 図９の４つの条件の下での４０００回の反復の後の１．０ｍｍのビーズを通過するＸ軸に沿った１次元のプロファイルを示す図 本実施形態に係るフルシステムオプティクスを用いて図９とは異なる条件で反復的に再構成された画像を示す図 図１１の４つの条件の下での４０００回の反復の後の０．２ｍｍのビーズを通過するＸ軸に沿った１次元のプロファイルを示す図
　図１を参照すると、本発明によるある実施形態において順投影と逆投影との両方（ＦＰＢＰ： forward projection and back projection）の両方でオプティクスを含むフルシステム光学モデル（ＳＯＭ）が示されている。図では示されてはいないが、このＳＯＭ－ＦＰＢＰは、画像を再構成するための所定の反復的（ＩＲ）アルゴリズムと組み合わせて用いられる。本実施形態において画像データは、反復的アルゴリズム（ＩＲ）を用いて、順投影と逆投影との両方でのフルシステム光学モデルに基づいて再構成される。以下、必要に応じて、反復的アルゴリズムを用いた順投影と逆投影との両方でのフルシステム光学モデルに基づく再構成をＩＲ－ＳＯＭ－ＦＰＢＰという略語で表現する。
　本実施形態においてシステム光学モデルという用語は、個別のシステム光学モデルを含む包括的な用語である。一般に、個別のシステム光学モデルは、１）ボクセルモデル、２）Ｘ線源モデル、３）検出器要素モデル、及び４）回転ブラー（rotational blur）またはビュー統合モデルを含む４つのグループに分類される。ボクセルモデルは、更に、複数のサブボクセル（マイクロレイ）、複数のブロブ（blob）、複数のスプライン（spline）を用いる特定の方法を含む。Ｘ線源モデルは、更に、複数のサブ線源（マイクロレイ）、複数のローパスフィルタ、及び分離可能な複数のフットプリントを用いる特定の方法を含む。検出器要素モデルは、更に、サブ検出器（マイクロレイ）と分離可能なフットプリントとを用いる特定の方法を含む。回転ブラーまたはビュー統合モデルは、更に、サブビュー（マイクロレイ）とローパスフィルタとを用いる特定の方法を含む。上述した特定の方法以外のモデルは、非システム光学モデルと考える。換言すると、非システム光学モデルは、単一レイまたは距離駆動法（distance driven method）などの方法を含む。距離駆動法は、システムジオメトリのある特性または様相を含むが、本実施形態においては、システム光学モデルの一部であるとは考えない。
　図１は、本実施形態に係るＩＲ－ＳＯＭ－ＦＰＢＰのシステム光学モデルを示している。図１のシステム光学モデルは、検出器１００に含まれる複数の検出器要素にそれぞれ対応する複数の検出器要素モデル１１０とＸ線源３００に対応するＸ線源モデル３１０とを含んでいる。また、検出器要素モデル１１０とＸ線源モデル３１０との間には、所定のボリューム２００を構成する複数のボクセルにそれぞれ対応する複数のボクセルモデル２１０が概念的に配置されている。
　図１に示すように、複数のマイクロレイがシステム光学モデルに組み込まれている。本実施形態に係るシステム光学モデルは、更に、マイクロレイに基づくシステム光学モデルとして特定される。マイクロレイに基づくシステム光学モデルでは、システムオプティクスが３次元Ｘ線ビームを定義することに用いられる。更に、マイクロレイに基づくシステム光学モデルでは、３次元Ｘ線ビームを正確に定義するために、所定の検出器に関するジオメトリと所定のボクセルに関するジオメトリと所定のＸ線源に関するジオメトリとを考慮する。Ｘ線源モデル３１０と検出器要素モデル１１０とボクセルモデル２１０とをサンプリングすることにより３次元Ｘ線ビームを満たす複数のマイクロレイが設定される。以下の説明のため、３次元Ｘ線ビームＢは、Ｒ１、Ｒ２、Ｒ３、及びＲ４などの所定の数の個別のマイクロレイを含む。３次元Ｘ線ビームＢが十分な数のマイクロレイで満たされている場合、従来技術による距離駆動法よりも正確に３次元Ｘ線ビーム内部のオブジェクトをサンプリングし、より高い解像度を有する画像を生じることができる。３次元Ｘ線ビームは、所定の反復的再構成法の順投影工程と逆投影工程との両方で用いられる。３次元Ｘ線ビームは、本実施形態に係る任意の反復的再構成アルゴリズムに適用可能である。
　順投影工程ではＸ線源モデルと検出器モデルとが関係し、逆投影工程ではＸ線源モデルとボクセルモデルとが関係する。マイクロレイに基づくシステム光学モデルを用いる場合、検出器要素が複数のマイクロ検出器点に概念的に分割され、ボクセルが複数のマイクロボクセル点に概念的に分割され、Ｘ線源が複数のマイクロＸ線源点に概念的に分割される。マイクロレイはマイクロ検出器点とマイクロボクセル点とマイクロＸ線源点とを結ぶ。マイクロレイは順投影や逆投影に供される。
　本実施形態においてシステムオプティクスは、順投影と逆投影との両方において適用される。必要に応じてシステムオプティクスは、逆投影のみに適用されても良い。システムオプティクスが逆投影のみに適用される場合、Ｘ線源のデコンボリューションが行われ、システムオプティクスが順投影だけに含まれる従来技術による方法よりも優れた解像度を提供する。解像度は、ＯＳ－ＳＡＲＴアルゴリズムを用いる順投影に加えて、逆投影におけるフルシステム光学モデルを用いるＩＲにおいて改善される。
　次に図２Ａを参照して本実施形態に係るＸ線源モデル３１０を説明する。Ｘ線源モデル３１０と検出器要素モデルとは順投影に組み込まれる。Ｘ線源モデル３１０のジオメトリには、アノード角θanode、Ｘｓ方向のＸ線源幅Ｓｘ、及びＺｓ方向のＸ線源高さＳｚが含まれる。Ｘ線焦点３００の表面３１０Ａは、本実施形態では、それぞれのグリッドの中心に位置するＮμS個のマイクロＸ線源点によってサンプリングされる。この実施形態でのグリッドの各々は、Ｘｓ軸およびＹｓ軸に沿った所定の寸法に従って等しく画定されている。ＮμS個のマイクロＸ線源点は、μＳＰ１、μＳＰ２、・・・、μＳＰＮというラベルが付された所定の数の離散的な個別のマイクロＸ線源点の集合物である。必要に応じて、Ｎμs個のマイクロＸ線源点は、任意の態様で分布している。換言すると、サンプリングパターンは、矩形のグリッドには限定されることはなく、ガウス分布やランダムパターンなど任意のパターンを含む。μＳＰ１からμＳＰＮまでのマイクロＸ線源点の各々は、関連する空間的なＸ線源応答重み値ＳＷｕを有する。
　次に図２Ｂを参照して本実施形態に係る検出器モデル１１０を説明する。検出器モデル１１０とＸ線源要素モデルとは順投影に組み込まれる。曲線状の検出器１００の検出器ジオメトリには、角チャネル幅Δγと、セグメント幅Ｗｓｅｇとが含まれる。検出器１００の表面１１０Ａは、それぞれのグリッドの中心に位置するＮμD個のマイクロ検出器点によってサンプリングされる。このグリッドの各々は、γ軸とＺD軸に沿った所定の寸法に従って等しく画定されている。ＮμD個のマイクロ検出器点は、μＤＰ１、μＤＰ２、・・・、μＤＰＮというラベルが付された所定の数の離散的な個別のマイクロ検出器点の集合物である。ＮμD個のマイクロ検出器点は、任意の態様で分布される。換言すると、サンプリングパターンは、矩形のグリッドに限定されることはなく、ガウス分布やランダムパターンなどの任意のパターンを含む。いずれにしても、μＤＰ１からμＤＰＮまでのマイクロ検出器点の各々は、関連する空間的なＸ線源応答重み値ＤＷｕを有する。
　図２Ｂに示すように、μＤＰ１からμＤＰＮまでのマイクロ検出器要素の各々は、個別のマイクロ検出器要素の間のアクティブでない空間に対応するためのデッドスペース（ピッチ）ＤＳを含む。例えば、デッドスペースＤＳは、セグメント幅ＷｓｅｇとＺD軸に沿ったアクティブなセグメント幅ＷｓｅｇＡＰとの間の差と、角チャネル幅Δγとγ軸に沿ったアクティブな角チャネル幅ΔγＡＰとの間の差とによって画定される。すなわち、デッドスペースＤＳの内部は、アクティブなセグメント幅ＷｓｅｇＡＰと角チャネル幅ΔγＡＰとにより画定される検出器要素のアクティブな開口部である。
　次に図３を参照して本実施形態に係るボクセルモデル２１０を説明する。ボクセルモデル２１０とＸ線源要素モデルとは逆投影に組み込まれる。ボクセル２００のボクセルジオメトリには、Ｘｖ方向のボクセル幅ＶｘＷＩＤＴＨ、Ｙｖ方向のボクセル長さＶｙＬＥＮＧＴＨ、及びＺｖ方向のボクセル高さＶｚＨＥＩＧＨＴが含まれる。ボクセル２００は、交点に位置するＮμv個のマイクロボクセル点によってサンプリングされる。交点の各々は、Ｘｖ軸、Ｙｚ軸、及びＺｖ軸に沿って等しく画定される。Ｎμv個のマイクロボクセル点は、μＶＰ１、μＶＰ２、・・・、μＶＰＮというラベルが付された所定の数の離散的な個別のマイクロボクセル点の集合物である。本実施形態では、Ｎμv個のマイクロボクセル点が任意の態様で分布される。換言すると、サンプリングパターンは、矩形のグリッドに限定されることはなく、ガウス分布やランダムパターンなどの任意のパターンを含む。μＶＰ１からμＶＰＮまでのマイクロボクセル点の各々は、関連する空間的なＸ線源応答重み値ＶＷｕｉを有する。図３のボクセルジオメトリは立方体であるが、異なる形状のボクセルジオメトリでも良い。
　次に図４を参照しながら、順投影と逆投影との両方を備えたフルシステム光学モデルに関するボクセルモデル２１０を説明する。このフルシステム光学モデルには、Ｘ線源３００、ボクセル２００Ａ、及び検出器１００が含まれる。ボクセル２００Ａは、総数がＮμvであるマイクロボクセル要素を有するマイクロボクセルに分割される。
　更に図４を参照すると、ボクセルは、図３に示す立方体に限定されない。本実施形態に係るボクセルの形状は、立方体でありグリッドによってサンプリングされるが、ボクセルが特定の形状には限定されない。所定のグリッドによってサンプリングされる３次元の立方体であるボクセルは、マイクロレイの重複が存在するため、非効率的である。図４に示されるボクセル２００Ａは、円柱形状を有しており、回転する２次元のマイクロボクセル平面２００Ｂを有している。円柱形状を有するボクセル２００Ａは、半径Ｒｖと高さＨｖとを有する。ボクセル２００Ａは、μＶＰ１からμＶＰ４までの４つのマイクロボクセル点を有し、マイクロボクセルの数を減少させる。回転する２次元のマイクロボクセル平面２００Ｂは、視野角とボクセル位置とに応じて回転し、Ｘ線源の中心ＳＣからボクセルの中心ＶＣまで延長するＸ線源ボクセル中心レイ
Figure JPOXMLDOC01-appb-M000001
に対して常に垂直である。回転する２次元マイクロボクセル平面２００Ｂは、このように、ボクセルのサンプリングを３次元から２次元に減少させる。
　本実施形態に係る順投影と逆投影との両方を備えたフルシステム光学モデルにおいては、３次元Ｘ線ビームの内部のマイクロレイの各々は、始点と終点との対応する対によって画定される。マイクロＸ線源点は、順投影と逆投影との両方に対する始点である。マイクロ検出器点は順投影に対する終点であり、マイクロボクセル点は逆投影に対する終点である。更に、マイクロレイは、平行に分布しても良いし、円錐状に分布しても良い。逆投影においてマイクロレイは、ボクセルを正しくサンプリングするため、円錐状に分布するのが好ましい。順投影においてマイクロレイは、円錐状に分布しても良いし平行に分布しても良い。
　次に図５Ａを参照しながら本実施形態に係る順投影と逆投影との両方を備えたフルシステム光学モデルにおける平行に分布するマイクロレイを説明する。図５Ａに示すように、μＲ１からμＲ４までのマイクロレイの総数ＮｕＲａｙｓは、μＳＰ１からμＳＰ４までのマイクロＸ線源点の総数ＮｕＳｒｃと、μＤＰ１からμＤＰ４までのマイクロ検出器点の総数ＮｕＤｅｔとに等しい。すなわち、順投影に適用されるマイクロＸ線源点とマイクロ検出器点との間には、１対１の接続が存在する。
　次に図５Ｂを参照しながら本実施形態に係るシステム光学モデルにおける円錐状に分布するマイクロレイを説明する。図５Ｂに示すように、μＲ１からμＲＮまでのマイクロレイの総数ＮｕＲａｙｓは、μＳＰ１からμＳＰ４までのマイクロＸ線源点の総数ＮｕＳｒｃと、μＤＰ１からμＤＰ４までのマイクロ検出器点の総数ＮｕＤｅｔとの積に等しい。すなわち、マイクロＸ線源点の各々は、逆投影に適用されるマイクロ検出器点の全てに接続されている。例えば、マイクロＸ線源点μＳＰ１からは、レイ束μＲｓ１において見られるように、μＤＰ１からμＤＰ４までのマイクロ検出器点の全てに、マイクロレイが分散される。同様に、μＲｓ２からμＲｓ４までのレイ束の各々は、μＳＰ２からμＳＰ４までのマイクロＸ線源点の各々から、μＤＰ１からμＤＰ４までのマイクロ検出器点の全てに４つのマイクロレイが接続する。なお、ＮｕＳｒｃとＮｕＤｅｔとの総数は、相互に異なっても良い。
　順投影では、３次元Ｘ線ビームに含まれる複数のマイクロレイｕの各々は個別に順投影される。ある検出器要素の値に対する最終的な順投影値は、３次元Ｘ線ビームに含まれる全てのＮｕＲａｙｓ個の順投影の平均である。３次元Ｘ線ビームの分散は、必要に応じて、平行または円錐状である。マイクロレイｕの順投影レイトレーシングは、応答重み値を用いて、次の方程式（１）により規定される。
Figure JPOXMLDOC01-appb-M000002
　なお、方程式（１）において、ＦＰ’uは検出器要素の値に対する最終的な順投影値であり、ＳＷuは関連する空間的なＸ線源応答重み値であり、ＶＷuは関連する空間的なボクセル応答重み値であり、ＤＷuは関連する空間的な検出器応答重み値であり、Ｌnuはボクセルｎｕを通過するレイｕの経路の長さであり、Ｖnuはレイｕと交差するボクセルである。検出器要素ｃｈ，ｓｅｇに対する最終的な順投影値は、次の方程式（２）によって定義されるように、個別的な順投影の平均である。
Figure JPOXMLDOC01-appb-M000003
　なお、方程式（２）において、ＦＰch，segはチャネルｃｈとセグメントｓｅｇとにより特定される検出器要素に対する最終的な順投影値であり、ＮuRaysはマイクロレイの総数である。本実施形態では、線積分値を計算するＳｉｄｄｏｎなど、任意のレイトレーシングアルゴリズムが利用可能である。
　次に図６を参照しながら、本実施形態に係る３次元最隣接レイトレーシングアルゴリズムを用いた順投影を説明する。この３次元最近接レイトレーシングアルゴリズムによると、体軸横断面が４つの象限に分割される。画像ボリュームは、順投影の前に２つの因数によりアップサンプリングされる。マイクロレイの体軸横断角度θに応じて、Ｘ平面またはＹ平面のいずれかが、３次元構成を２次元構成に縮小するための「アンカー（anchor）」平面として用いられる。例えば、このグラフには、マイクロレイの与えられている体軸横断角度θは４５度から１３５度までの間であることが示されている。与えられている体軸横断角度θが、４５°＜θ≦１３５°であるから、最近接のものは、次の平面に従って、Ｙ平面に固定される。
　－４５＜ｑ＜４５　ｘｍａｘ，・・・，ｘｍｉｎのアンカー平面を用いる　
　　４５＜ｑ＜１３５　ｙｍｉｎ，・・・，ｙｍａｘのアンカー平面を用いる　
　１３５＜ｑ＜２２５　ｘｍｉｎ，・・・，ｘｍａｘのアンカー平面を用いる　
　２２５＜ｑ＜３１５　ｙｍｉｎ，・・・，ｙｍａｘのアンカー平面を用いる
　ここで、ｘｍｉｎ、ｘｍａｘ、ｙｍｉｎ、ｙｍａｘは、画像体積を用いて、マイクロレイの入口点と出口点とから決定される。結果的に、上述の方程式（１）は、交点の長さＬθを用いて次の方程式（３）として定義される。なお、Ｌθは、特定のレイに対する定数である。そして、マイクロレイｕに対する線積分は、次のようになる。
Figure JPOXMLDOC01-appb-M000004
　マイクロレイｕのレイトレーシング順投影は、下記の方程式（４）により規定されても良い。なお、任意のマイクロＸ線源点をｘus、ｙus、ｚusとし、任意のマイクロ検出器点をｘud、ｙud、ｚudとし、対応する体軸横断角度をφとし、マイクロレイｕの体軸横断交点の長さをＬφとし、ボクセルを一定とする。この場合、マイクロＸ線源点ｕｓからマイクロ検出器点ｕｄまでの画像ボリュームを通過するマイクロレイに対する順投影は、以下の（４）により規定される。
Figure JPOXMLDOC01-appb-M000005
　ここで、ｎはマイクロレイが交差するボクセルを表し、Ｎはマイクロレイが交差するボクセルの総数である。検出器要素ｍに対する全体の順投影値は、次の方程式（５）により規定される。
Figure JPOXMLDOC01-appb-M000006
　ここで、ＮuSはマイクロＸ線源点（ＮＳｘ・ＮＳｚ）の総数であり、ＮuDはマイクロ検出器点（ＮｕＣｈ・ＮｕＳｅｇ）の総数である。図２Ａには、ＮＳｘがＸ線源モデルのＸ方向に沿ったマイクロＸ線源要素の数であり、ＮＳｚがＸ線源モデルのＺ方向に沿ったマイクロＸ線源要素の数であることが示されている。図２Ｂには、ＮｕＣｈが検出器モデルのγ方向に沿ったマイクロ検出器要素の数であり、ＮｕＳｅｇが検出器モデルのＺ方向に沿ったマイクロ検出器要素の数であることが示されている。
　次に図７を参照しながら本実施形態に係る逆投影におけるＸ線源とボクセルとに関するマイクロレイを説明する。逆投影については、Ｘ線源モデルは順投影の場合と同じであり、検出器モデルはボクセルモデルによって置き換えられる。逆投影では、３次元Ｘ線ビームの中の各々のマイクロレイｕは個別的に逆投影される。Ｘ線源モデル及び検出器モデルと同じように、ボクセルは、所定の数のマイクロボクセルまたはマイクロボクセル点に分割される。ボクセルモデル２１０におけるｖｏｘｅｌ［１］からｖｏｘｅｌ［ｎ］までのマイクロボクセル点の各々は、Ｘ線源３００とボクセル２００とに亘って平均化されたＮuRays個のマイクロレイを受け取る。検出器１００は、マイクロ検出器点の各々において投影値又は再投影値ＰＤ［ｃｈu，ｓｅｇu］を有している。
　更新されたボクセル値は、３次元Ｘ線ビームに含まれるＮuRays個のマイクロレイでの逆投影の平均値である。逆投影では、円錐状のビームが用いられる。この点で、マイクロレイｕのレイトレーシング順投影は、応答重み値を有する方程式（６）により規定される。
Figure JPOXMLDOC01-appb-M000007
　ここで、ｃｈuはマイクロレイｕのチャネル位置であり、ｓｅｇuはセグメント位置であり、ＰＤ［ｃｈu，ｓｅｇu］はｃｈu及びｓｅｇuに対する投影値又は再投影値である。ＳＷuは関連する空間的Ｘ線源応答重み値であり、ＶＷuは関連する空間的ボクセル応答重み値であり、ＤＷuは関連する空間的検出器応答重み値である。ある実施形態では最近接の値が用いられ、別の実施形態では双直線補間が用いられる。
　効率的な逆投影のためには、図４に示すように、ボクセルは、半径がＲVoxelで高さがＨVoxelである円柱と、全部でＮｕＶ個の位置を備えたマイクロボクセルまたはマイクロボクセル点を有する２次元逆投影平面とによってモデル化される。平面２００Ｂは、ある視野角とあるボクセル位置とを有するように回転し、ボクセルの中心レイへのＸ線源の中心に対して常に垂直である。Ｘ線源の中心からボクセルの中心までの単一のレイを通過する逆投影ではなく、ＮuS・ＮuV本のマイクロレイが逆投影される。ＮuS・ＮuV本のマイクロレイの逆投影値は、次の方程式（７）により規定されるように、Ｘ線源及びボクセルに亘り平均化される。
Figure JPOXMLDOC01-appb-M000008
　ここで、ＢＰｎは逆投影されたボクセルｎのボクセル値であり、ｍは３次元Ｘ線ビームに含まれるマイクロレイが交差した検出器要素を表し、Ｍは３次元Ｘ線ビームが交差した検出器要素の総数であり、ａn,mは逆投影演算子である。
　ＩＲ－ＳＯＭを用いた再構成においてビーズの直径が回復可能か否かを確認するため、シミュレートされた投影値又は再投影値が、誇張された６．０ｍｍのＸ線源と１．０ｍｍのビーズを用いて発生された。また、実際的な条件下で解像度の改善を調べるために、上記の条件に加え、１．１ｍｍのＸ線源と０．２ｍｍのビーズとで投影値又は再投影値をシミュレーションにより発生した。
　次に、図８を参照しながら、本実施形態に係る再構成プロセスの流れを説明する。ステップＳ１００では、コンピュータ断層撮影（ＣＴ）装置などの所定のモダリティを用いるスキャナを用いて生データが収集される。ステップＳ１１０では、収集された生データに基づいて、所定のアルゴリズムに従ってシード画像を推定する。すなわち、ステップＳ１１０においては、生データに基づいて初期的にシード画像が再構成される。ステップＳ１２０では、推定されたシード画像から、オーダードサブセット同時代数的再構成法（ＯＳ－ＳＡＲＴ：ordered-subset simultaneous algebraic reconstruction technique）などの所定の逐次再構成アルゴリズムを用いて順投影値のデータ（再投影値データ）を発生する。具体的には、推定されたシード画像を、システム光学モデルと非システム光学モデルとの何れかを用いて順投影して再投影値データを発生する。ステップＳ１３０では、フルシステム光学モデル再構成において再投影値データは、システム光学モデルを用いてシード画像を更新するために逆投影される。ステップＳ１２０およびＳ１３０における上述した順投影および逆投影は、所定の条件が満たされるまで反復的に適用される。この条件は、ステップＳ１４０で判断される。ステップＳ１４０において、所定の反復回数などの条件が満たされていると、プロセスは終了する。他方で、ステップＳ１４０において条件が満たされていない場合には、プロセスは、ステップＳ１２０に戻る。
　システム光学モデルを用いる反復的再構成（ＩＲ－ＳＯＭ）は、ＯＳ－ＳＡＲＴ反復的再構成アルゴリズムに組み込まれており、緩和パラメータλ＝０．５である。ある例では、ファンビームだけが実装され、体軸横断方向のＸＹ解像度だけが調べられた。ファンビームのシミュレーションではフル３次元のジオメトリが含まれている（すなわち、ファンビームはＸ線源と検出器とに起因する厚さを有し、単純に厚さがゼロの２次元平面ではない）。
　２つのシミュレーションの結果を以下で示す。シミュレーション（１）では、直径が１．０ｍｍのビーズと６．０ｍｍのＸ線源とが用いられる。シミュレーション（２）では、直径が０．２ｍｍのビーズと１．１ｍｍのＸ線源とが用いられる。シミュレーション（１）では、ビーズよりもサイズがはるかに大きなＸ線源を組み込むことにより、システムオプティクスを含むことでＸ線源のブラーイングが補償されＸ線源よりもはるかに小さなオブジェクトを回復できるか否かを評価する。シミュレーション（２）では、検出器のサイズよりも小さなビーズを用いて現実的な焦点サイズを表すことにより、体軸横断方向の解像度を評価する。
　４つの再構成アルゴリズムを用いてビーズが再構成された。すなわち、ＬＡＫＳカーネルを用いた標準的なフィルタ補正逆投影（ＦＢＰＪ）、上述したフルシステムオプティクスを用いたＩＲ（ＩＲ－ＳＯＭ－ＦＰＢＰ）、システムオプティクスを用いないペンシルビームＩＲ（ＩＲ－Ｐ）、および順投影にだけシステムオプティクスを備えているＩＲ（ＩＲ－ＳＯＭ－ＦＰ）である。
　表１には、２つのシミュレーションと標準的なファンビーム構成の場合とに用いられたパラメータがまとめられている。
Figure JPOXMLDOC01-appb-T000009
　次に図９を参照すると、フルシステムオプティクスを含む反復的に再構成された画像が、本実施形態に係る所定の条件の下での比較のために示されている。再構成された画像は、６．０ｍｍのＸ線源と１．０ｍｍのビーズとを条件に再構成されている。ＦＢＰＪとＩＲとの組み合わせは、画像を再構成し正規化するために用いられている。ＦＢＰＪとＩＲ－ＰとＩＲ－ＳＯＭ－ＦＰとの条件の下での３つの再構成された画像は、誇張されたＸ線源のサイズを十分にサンプリングしておらず、ビーズが歪んで描出されている。ＩＲ－ＳＯＭ－ＦＰの条件の下で再構成された画像だけが、最も感度が高く、最も歪んでいる。ＩＲ－ＳＯＭ－ＦＰＢＰの条件の下では、ビーズのサイズを完全に回復することによって、ビーズを正確に再構成している。このＩＲ－ＳＯＭ－ＦＰＢＰ画像およびプロットは、いくつかの小さなサイドローブを示しているが、これは、順投影と逆投影との不一致に起因する可能性がある。ＩＲ－ＳＯＭ－ＦＰＢＰの条件の下で再構成された画像は、フィルタ補正逆投影法（ＦＢＰＪ）、システムモデルを用いないペンシルビームＩＲ法（ＩＲ－Ｐ）、順投影だけにシステムモデルを用いるＩＲ法の条件の下で再構成された画像よりも実質的に優れた解像度を有する。
　図１０は、図９の４つの条件の下での４０００回の反復の後の１．０ｍｍのビーズを通過するｘ軸に沿った１次元のプロファイルを示す図である。図９に関して上述したように、ＩＲ－ＳＯＭ－ＦＰＢＰの条件の下でのプロファイルは、フィルタ補正逆投影法（ＦＢＰＪ）、システムモデルを用いないペンシルビームＩＲ法（ＩＲ－Ｐ）、及び順投影だけにシステムモデルを用いるＩＲ法の条件の下で再構成された画像よりも、実質的に優れた解像度を有する。
　次に図１１を参照しながら、本実施形態に係るフルシステムオプティクスを用いて図９とは異なる条件で反復的に再構成された画像を説明する。再構成された画像は、１．１ｍｍのＸ線源と０．２ｍｍのビーズとを条件に再構成されている。ＦＢＰＪとＩＲとの組み合わせによって、画像が再構成され、正規化されている。ＩＲ－ＳＯＭ－ＦＰＢＰが最高の解像度を有しており、ＦＷＨＭは０．４８ｍｍである。これは、ＦＢＰＪ（０．８７ｍｍ）のほとんど半分であり、ＩＲ－Ｐ（０．６３ｍｍ）やＩＲ－ＳＯＭ－ＦＰＢＰ（０．５９ｍｍ）よりも相当に優れている。誇張されたＸ線源のシミュレーションに関しては、ＩＲ－ＳＯＭ－ＦＰＢＰがサイドローブを示している。更に、ＩＲ－ＰとＩＲ－ＳＯＭ－ＦＰによる再構成ではビーズの真ん中に「ディンプル」が示されているが、ＩＲ－ＳＯＭ－ＦＰＢＰの場合はそのようなことはない。ＩＲ－ＳＯＭ－ＦＰＢＰのために反復の回数を１０，０００回まで増加させたことによる効果もまた調べられたが、その結果、ＦＷＨＭの解像度が０．４１ｍｍまで更に改善し、サイドローブのレベルがほぼ一定に留まっていることが判明している。
　図１２は、本発明による上述した４つの条件の下での４０００回反復した後の０．２ｍｍのビーズを通過するｘ軸に沿った１次元のプロファイルを図解するグラフである。図１１に関して上述したように、ＩＲ－ＳＯＭ－ＦＰＢＰの条件の下でのプロファイルは、フィルタ補正逆投影法（ＦＢＰＪ）、システムモデルを用いないペンシルビームＩＲ法（ＩＲ－Ｐ）、及び順投影だけにシステムモデルを用いるＩＲ法の条件の下での再構成された画像よりも、実質的に優れた解像度を有する。１０，０００回の反復の後でのＩＲ－ＳＯＭ－ＦＰＢＰの条件の下でのプロファイルは、破線を用いてプロットされている。
　かくして、本実施形態によれば、優れた画質を有する画像を再構成することが可能となる。
　しかし、本発明の多数の特性および効果に関して、本発明の構造および機能の詳細と共に以上の説明で論じられてきたが、ここでの開示は単に説明目的のものに過ぎないことを理解すべきである。また、詳細において、特に、部分の形状、サイズおよび構成や、ソフトウェア、ハードウェアまたはそれら両者の組み合わせによる実装に関して変更を行うことが可能であるが、そのような変更は、後述の特許請求の範囲に記載のある用語の広く一般的な意味によって指示される完全な範囲までを含む本発明の原理の範囲内にあることも理解すべきである。
　１００…検出器、１１０…検出器要素モデル、２００…ボクセル、２１０…ボクセルモデル、３００…Ｘ線源、３１０…Ｘ線源モデル、Ｒ１，Ｒ２，Ｒ３，Ｒ４…マイクロレイ

Claims (22)

  1. 　スキャナにより被検体をＸ線でスキャンして前記被検体に関する生データを収集する収集工程と、
    　前記生データに基づいて画像を推定する推定工程と、
    　システム光学モデルと非システム光学モデルとのうちの一方を用いて再投影データを発生するために順投影する順投影工程と、前記システム光学モデルを用いて前記再投影データに基づいて画像を更新するために逆投影する逆投影工程とを反復的に実行する反復工程と、
    　を具備する画像再構成方法。
  2. 　前記システム光学モデルは、ボクセルモデル、Ｘ線焦点モデル、検出器要素モデル、及び回転ブラーモデルの少なくとも一つを含む、請求項１記載の画像再構成方法。
  3. 　前記システム光学モデルは、ボクセルモデル、Ｘ線焦点モデル、検出器要素モデル、及び回転ブラーモデルから構成される、請求項１記載の画像再構成方法。
  4. 　前記順投影工程と前記逆投影工程とにおいて前記システム光学モデルが用いられる、請求項１記載の画像再構成方法。
  5. 　前記順投影工程において前記非システム光学モデルが用いられ、前記逆投影工程において前記システム光学モデルが用いられる、請求項１記載の画像再構成方法。
  6. 　前記システム光学モデルは複数のマイクロレイを用いる、請求項１記載の画像再構成方法。
  7. 　前記複数のマイクロレイは平行に投影される、請求項６記載の画像再構成方法。
  8. 　前記複数のマイクロレイは円錐状に投影される、請求項６記載の画像再構成方法。
  9. 　前記システム光学モデルはＸ線焦点モデルを含み、
    　前記Ｘ線焦点モデルは空間的に離散的に設定された複数の第１のサンプル点を含み、
    　前記複数の第１のサンプル点には重み値が個別に設定される、
    　請求項６記載の画像再構成方法。
  10. 　前記システム光学モデルは検出器モデルを含み、
    　前記検出器モデルは空間的に離散的に設定された複数の第２のサンプル点を含み、
    　前記複数の第２のサンプル点には重み値が個別に設定される、
    　請求項９記載の画像再構成方法。
  11. 　前記システム光学モデルはボクセルモデルを含み、
    　前記ボクセルモデルは空間的に離散的に設定された複数の第３のサンプル点を含み、
    　前記複数の第３のサンプル点には重み値が個別に設定される、
    　請求項１０記載の画像再構成方法。
  12. 　被検体をＸ線でスキャンして前記被検体に関する生データを収集するスキャナと、
    　前記生データに基づいて画像データを再構成する再構成部であって、前記生データに基づいて画像データを推定し、システム光学モデルと非システム光学モデルとのうちの一方を用いて再投影データを発生するための順投影処理と前記システム光学モデルを用いて前記再投影データに基づいて画像データを更新するための逆投影処理とを反復的に実行する再構成部と、
    　を具備する画像再構成システム。
  13. 　前記システム光学モデルは、ボクセルモデル、Ｘ線焦点モデル、検出器要素モデル、及び回転ブラーモデルの少なくとも一つを含む、請求項１２記載の画像再構成システム。
  14. 　前記システム光学モデルは、ボクセルモデル、Ｘ線焦点モデル、検出器要素モデル、及び回転ブラーモデルから構成される、請求項１２記載の画像再構成システム。
  15. 　前記再構成部は、前記順投影処理と前記逆投影処理とにおいて前記システム光学モデルを用いる、請求項１２記載の画像再構成システム。
  16. 　前記再構成部は、前記順投影処理において前記非システム光学モデルを用い、前記逆投影処理において前記システム光学モデルを用いる、請求項１２記載の画像再構成システム。
  17. 　前記システム光学モデルは複数のマイクロレイを用いる、請求項１２記載の画像再構成システム。
  18. 　前記複数のマイクロレイは平行に投影される、請求項１７記載の画像再構成システム。
  19. 　前記複数のマイクロレイは円錐状に投影される、請求項１７記載の画像再構成システム。
  20. 　前記システム光学モデルはＸ線焦点モデルを含み、
    　前記Ｘ線焦点モデルは空間的に離散的に設定された複数の第１のサンプル点を含み、
    　前記複数の第１のサンプル点には重み値が個別に設定される、
    　請求項１７記載の画像再構成システム。
  21. 　前記システム光学モデルは検出器モデルを含み、
    　前記検出器モデルは空間的に離散的に設定された複数の第２のサンプル点を含み、
    　前記複数の第２のサンプル点には重み値が個別に設定される、
    　請求項２０記載の画像再構成システム。
  22. 　前記システム光学モデルはボクセルモデルを含み、
    　前記ボクセルモデルは空間的に離散的に設定された複数の第３のサンプル点を含み、
    　前記複数の第３のサンプル点には重み値が個別に設定される、
    　請求項２１記載の画像再構成システム。
PCT/JP2014/052115 2013-01-31 2014-01-30 モデルに基づく反復的再構成のための逆投影と順投影との少なくとも一方におけるシステムオプティクス WO2014119664A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201480006906.2A CN104968275B (zh) 2013-01-31 2014-01-30 基于模型的用于迭代重建的反投影和正投影的至少一方中的系统光学

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US13/756,082 2013-01-31
US13/756,082 US9153048B2 (en) 2013-01-31 2013-01-31 System optics in at least in one of backprojection and forward projection for model-based iterative reconstruction

Publications (1)

Publication Number Publication Date
WO2014119664A1 true WO2014119664A1 (ja) 2014-08-07

Family

ID=51223001

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/052115 WO2014119664A1 (ja) 2013-01-31 2014-01-30 モデルに基づく反復的再構成のための逆投影と順投影との少なくとも一方におけるシステムオプティクス

Country Status (4)

Country Link
US (1) US9153048B2 (ja)
JP (1) JP6351986B2 (ja)
CN (1) CN104968275B (ja)
WO (1) WO2014119664A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6001198B2 (ja) * 2014-10-01 2016-10-05 日本碍子株式会社 層状複水酸化物を用いた電池

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10169909B2 (en) * 2014-08-07 2019-01-01 Pixar Generating a volumetric projection for an object
US9709513B2 (en) * 2014-09-30 2017-07-18 Hexagon Metrology, Inc. System and method for measuring an object using X-ray projections
US10417795B2 (en) * 2015-04-08 2019-09-17 Canon Medical Systems Corporation Iterative reconstruction with system optics modeling using filters
JP2017099755A (ja) * 2015-12-03 2017-06-08 キヤノン株式会社 情報処理装置、画像再構成方法、およびプログラム
WO2017166187A1 (en) * 2016-03-31 2017-10-05 Shanghai United Imaging Healthcare Co., Ltd. System and method for image reconstruction
US20170309060A1 (en) * 2016-04-21 2017-10-26 Honeywell International Inc. Cockpit display for degraded visual environment (dve) using millimeter wave radar (mmwr)
US11335038B2 (en) * 2019-11-04 2022-05-17 Uih America, Inc. System and method for computed tomographic imaging
CN115330900B (zh) * 2022-10-13 2022-12-27 中加健康工程研究院(合肥)有限公司 一种用于pet图像重构的系统矩阵快速迭代计算方法
CN118022200A (zh) * 2022-11-11 2024-05-14 中硼(厦门)医疗器械有限公司 治疗计划系统、重叠自动检查方法及治疗计划的制定方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050286749A1 (en) * 2004-06-23 2005-12-29 General Electric Company System and method for image reconstruction

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5414623A (en) * 1992-05-08 1995-05-09 Iowa State University Research Foundation Optoelectronic system for implementation of iterative computer tomography algorithms
US7227982B2 (en) 2002-04-15 2007-06-05 General Electric Company Three-dimensional reprojection and backprojection methods and algorithms for implementation thereof
US7356113B2 (en) * 2003-02-12 2008-04-08 Brandeis University Tomosynthesis imaging system and method
DE102005051620A1 (de) * 2005-10-27 2007-05-03 Siemens Ag Verfahren zur Rekonstruktion einer tomographischen Darstellung eines Objektes
US8135186B2 (en) * 2008-01-25 2012-03-13 Purdue Research Foundation Method and system for image reconstruction
US8116426B2 (en) 2008-11-11 2012-02-14 Kabushiki Kaisha Toshiba Computed tomography device and method using circular-pixel position-adaptive interpolation
WO2011007270A1 (en) * 2009-07-14 2011-01-20 Koninklijke Philips Electronics, N.V. Image reconstruction including shift-variant blur compensation
US8761478B2 (en) * 2009-12-15 2014-06-24 General Electric Company System and method for tomographic data acquisition and image reconstruction
US8731266B2 (en) * 2009-12-17 2014-05-20 General Electric Company Method and system for correcting artifacts in image reconstruction
DE102010006774A1 (de) * 2010-02-04 2011-08-04 Siemens Aktiengesellschaft, 80333 CT-Messung mit Mehrfachröntgenquellen
US8971599B2 (en) * 2010-12-20 2015-03-03 General Electric Company Tomographic iterative reconstruction
US9076255B2 (en) * 2011-05-31 2015-07-07 General Electric Company Method and system for reconstruction of tomographic images
US8416914B2 (en) * 2011-07-08 2013-04-09 General Electric Company System and method of iterative image reconstruction for computed tomography
US8923583B2 (en) * 2012-06-22 2014-12-30 General Electric Company Methods and systems for performing model-based iterative reconstruction

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050286749A1 (en) * 2004-06-23 2005-12-29 General Electric Company System and method for image reconstruction

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
DE MAN: "Distance-driven projection and backprojection", NUCLEAR SCIENCE SYMPOSIUM CONFERENCE RECORD, vol. 3, 2002, pages 1477 - 1480 *
HSIEH J.: "Recent advances in CT image reconstruction", CURRENT RADIOLOGY REPORTS, vol. 1, no. 1, 16 January 2013 (2013-01-16), pages 39 - 51 *
ILMAR A. HEIN ET AL.: "System optics in both backprojection and forward projection for model -based iterative reconstruction", PROC. OF SPIE, vol. 8313, pages 83133L - 1 - 83133L-8 *
LONG, Y.: "3D forward and back-projection for X-ray CT using separable footprints", IEEE TRANSACTIONS ON MEDICAL IMAGING, vol. 29, no. 11, 2010, pages 1839 - 1850 *
WANG, GE: "Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART", JOURNAL OF X-RAY SCIENCE AND TECHNOLOGY, vol. 12, no. 3, 2004, pages 169 - 177 *
WUNDERLICH ADAM: "Exact and efficient computation of noise covariance for fan-beam FBP reconstructions that use rebinning to parallel-beam geometry", SPIE MEDICAL IMAGING, vol. 8313, 2012, pages 831323 - 1 - 831323-9 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6001198B2 (ja) * 2014-10-01 2016-10-05 日本碍子株式会社 層状複水酸化物を用いた電池
JPWO2016051934A1 (ja) * 2014-10-01 2017-04-27 日本碍子株式会社 層状複水酸化物を用いた電池

Also Published As

Publication number Publication date
US20140212018A1 (en) 2014-07-31
US9153048B2 (en) 2015-10-06
CN104968275A (zh) 2015-10-07
JP2014147772A (ja) 2014-08-21
JP6351986B2 (ja) 2018-07-04
CN104968275B (zh) 2018-01-02

Similar Documents

Publication Publication Date Title
JP6351986B2 (ja) モデルに基づく反復的再構成のための逆投影と順投影との少なくとも一方におけるシステムオプティクス
Thibault et al. A three‐dimensional statistical approach to improved image quality for multislice helical CT
Xu et al. A practical cone-beam CT scatter correction method with optimized Monte Carlo simulations for image-guided radiation therapy
Hsieh et al. Recent advances in CT image reconstruction
US8897528B2 (en) System and method for iterative image reconstruction
JP5348855B2 (ja) 対象の画像再構成方法およびその方法を実施するための装置
US20070297656A1 (en) System and method for iterative image reconstruction
EP3673457B1 (en) A method of generating an enhanced tomographic image of an object
US10722178B2 (en) Method and apparatus for motion correction in CT imaging
EP3447731A1 (en) A method of generating an enhanced tomographic image of an object
US11935160B2 (en) Method of generating an enhanced tomographic image of an object
US6987829B2 (en) Non-iterative algebraic reconstruction technique for tomosynthesis
US8335358B2 (en) Method and system for reconstructing a medical image of an object
WO2010011676A2 (en) Incorporation of mathematical constraints in methods for dose reduction and image enhancement in tomography
Lee et al. Interior tomography using 1D generalized total variation. Part II: Multiscale implementation
Sunnegårdh et al. Regularized iterative weighted filtered backprojection for helical cone‐beam CT
KR102348139B1 (ko) 이중 해상도의 관심 영역 내외 투영 데이터를 이용한 체내 단층 촬영 방법 및 시스템
Liu et al. Singular value decomposition-based 2D image reconstruction for computed tomography
Varslot et al. Fast high-resolution micro-CT with exact reconstruction methods
Sunnegårdh Iterative filtered backprojection methods for helical cone-beam CT
US7440602B2 (en) Methods, apparatus, and software to compensate for failed or degraded components
Vlasov et al. An a priori information based algorithm for artifact preventive reconstruction in few-view computed tomography
TWI613998B (zh) 斷層合成影像邊緣假影抑制方法
Fu et al. Modeling and estimation of detector response and focal spot profile for high-resolution iterative CT reconstruction
Zeng Detector blurring and detector sensitivity compensation for a spinning slat collimator

Legal Events

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

Ref document number: 14746257

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14746257

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP