US8885975B2 - Method and apparatus for iterative reconstruction - Google Patents

Method and apparatus for iterative reconstruction Download PDF

Info

Publication number
US8885975B2
US8885975B2 US13/531,082 US201213531082A US8885975B2 US 8885975 B2 US8885975 B2 US 8885975B2 US 201213531082 A US201213531082 A US 201213531082A US 8885975 B2 US8885975 B2 US 8885975B2
Authority
US
United States
Prior art keywords
image
objective function
auxiliary variable
sub
sequence
Prior art date
Legal status (The legal status 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 status listed.)
Active, expires
Application number
US13/531,082
Other languages
English (en)
Other versions
US20130343672A1 (en
Inventor
Zhou Yu
Bruno Kristiaan Bernard De Man
Jean-Baptiste Thibault
Debashish Pal
Lin Fu
Charles A. Bouman
Ken Sauer
Sathish Ramani
Jeffrey A. Fessler
Somesh Srivastava
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
General Electric Co
University of Michigan
University of Notre Dame
Purdue Research Foundation
Original Assignee
General Electric Co
University of Michigan
University of Notre Dame
Purdue Research Foundation
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
Priority to US13/531,082 priority Critical patent/US8885975B2/en
Application filed by General Electric Co, University of Michigan, University of Notre Dame , Purdue Research Foundation filed Critical General Electric Co
Assigned to GENERAL ELECTRIC COMPANY reassignment GENERAL ELECTRIC COMPANY ASSIGNMENT OF ASSIGNORS INTEREST (SEE DOCUMENT FOR DETAILS). Assignors: DE MAN, BRUNO KRISTIAAN BERNARD, FU, LIN, PAL, DEBASHISH, SRIVASTAVA, SOMESH, THIBAULT, JEAN-BAPTISTE, YU, ZHOU
Assigned to THE REGENTS OF THE UNIVERSITY OF MICHIGAN reassignment THE REGENTS OF THE UNIVERSITY OF MICHIGAN ASSIGNMENT OF ASSIGNORS INTEREST (SEE DOCUMENT FOR DETAILS). Assignors: FESSLER, JEFFREY A., RAMANI, SATHISH
Assigned to NATIONAL INSTITUTES OF HEALTH (NIH), U.S. DEPT. OF HEALTH AND HUMAN SERVICES (DHHS), U.S. GOVERNMENT reassignment NATIONAL INSTITUTES OF HEALTH (NIH), U.S. DEPT. OF HEALTH AND HUMAN SERVICES (DHHS), U.S. GOVERNMENT CONFIRMATORY LICENSE (SEE DOCUMENT FOR DETAILS). Assignors: UNIVERSITY OF MICHIGAN
Assigned to UNIVERSITY OF NOTRE DAME DU LAC reassignment UNIVERSITY OF NOTRE DAME DU LAC ASSIGNMENT OF ASSIGNORS INTEREST (SEE DOCUMENT FOR DETAILS). Assignors: SAUER, KEN
Assigned to PURDUE RESEARCH FOUNDATION reassignment PURDUE RESEARCH FOUNDATION ASSIGNMENT OF ASSIGNORS INTEREST (SEE DOCUMENT FOR DETAILS). Assignors: BOUMAN, CHARLES A.
Priority to JP2013126219A priority patent/JP6280700B2/ja
Priority to NL2011005A priority patent/NL2011005B1/en
Priority to DE102013106467.1A priority patent/DE102013106467A1/de
Priority to CN201310249820.XA priority patent/CN103514629B/zh
Publication of US20130343672A1 publication Critical patent/US20130343672A1/en
Publication of US8885975B2 publication Critical patent/US8885975B2/en
Application granted granted Critical
Active legal-status Critical Current
Adjusted expiration legal-status Critical

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
    • 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/424—Iterative

Definitions

  • the subject matter disclosed herein relates generally to imaging systems, and more particularly, to a method and apparatus for reconstructing an image using iterative techniques.
  • MBIR model based iterative reconstruction
  • MBIR algorithms typically reconstruct an image by first forming an objective function that incorporates an accurate system model, statistical noise model, and prior model. With the objective function in hand, the image is then reconstructed by computing an estimate that minimizes the objective function, which is typically performed using an iterative optimization algorithm. Examples of some of such iterative optimization algorithms include iterative coordinate descent (ICD), variations of expectation maximization (EM), conjugate gradients (CG), and ordered subsets (OS).
  • ICD iterative coordinate descent
  • EM expectation maximization
  • CG conjugate gradients
  • OS ordered subsets
  • a method for reconstructing an image of an object that includes image elements.
  • the method includes accessing measurement data associated with the image elements, introducing an auxiliary variable to transform an original problem of reconstructing the image to a constrained optimization problem, and solving the constrained optimization problem using a method of multipliers to create a sequence of sub-problems and solve the sequence of sub-problems.
  • Solving the sequence of sub-problems includes reconstructing the image by optimizing a first objective function.
  • the first objective function is optimized by iteratively solving a nested sequence of approximate optimization problems.
  • An inner loop iteratively optimizes a second objective function approximating the first objective function.
  • An outer loop utilizes the solution of the second objective function to optimize the first objective function.
  • a non-transitory computer readable medium is provided.
  • the non-transitory computer readable medium is programmed to instruct a computer to access measurement data associated with image elements, and create a sequence of sub-problems from an original problem of reconstructing an image of an object.
  • the sub-problems include a reconstruction of the image, an auxiliary variable, and a transform of a Lagrangian multiplier.
  • the computer is also instructed to update the reconstruction of the image and the auxiliary variable.
  • the computer is also instructed to update the transform of the Lagrangian multiplier using the updated reconstruction of the image and the updated auxiliary variable.
  • the computer is also instructed to iteratively update the reconstruction of the image using the updated auxiliary variable and the updated transform of the Lagrangian multiplier.
  • an imaging system in another embodiment, includes a detector array and a computer coupled to the detector array.
  • the computer is configured to access measurement data associated with image elements, introduce an auxiliary variable to transform an original problem of reconstructing an image of an object to a constrained optimization problem, and solve the constrained optimization problem using a method of multipliers to create a sequence of sub-problems and solve the sequence of sub-problems.
  • Solving the sequence of sub-problems includes reconstructing the image by optimizing a first objective function.
  • the first objective function is optimized by iteratively solving a nested sequence of approximate optimization problems.
  • An inner loop iteratively optimizes a second objective function approximating the first objective function.
  • An outer loop utilizes the solution of the second objective function to optimize the first objective function.
  • FIG. 1 is a flowchart of an exemplary method for reconstructing an image of an object in accordance with various embodiments.
  • FIG. 2 is another flowchart illustrating an exemplary method of performing a solving step of the method shown in FIG. 1 in accordance with various embodiments.
  • FIG. 3 is another flowchart illustrating an exemplary method of performing an iterative updating step of the method shown in FIG. 2 .
  • FIG. 4 is another flowchart illustrating an exemplary method of performing an optimizing step of the method shown in FIG. 3 .
  • FIG. 5 is a visual representation of exemplary difference images.
  • FIG. 6 is a pictorial view of an exemplary multi-modality imaging system formed in accordance with various embodiments.
  • FIG. 7 is a block schematic diagram of the system illustrated in FIG. 6 .
  • the functional blocks are not necessarily indicative of the division between hardware circuitry.
  • one or more of the functional blocks e.g., processors or memories
  • the programs may be stand alone programs, may be incorporated as subroutines in an operating system, may be functions in an installed software package, and the like. It should be understood that the various embodiments are not limited to the arrangements and instrumentality shown in the drawings.
  • an image element shall refer to an element of an object within the image space array.
  • an image element may include an image pixel or picture element that can correspond to a single voxel or volume element in a three-dimensional (3D) reconstruction.
  • the image is reconstructed by optimizing an objective function.
  • an objective function, or cost function generally includes a model of the imaging system used to acquire the imaging data, of the noise associated with the imaging system, and of the properties of the image to be reconstructed.
  • the terms ‘optimize’, ‘minimize’, and ‘maximize’ are interchangeable.
  • the terms ‘objective function’ and ‘cost function’ are interchangeable.
  • Various embodiments provide a relatively fast convergent parallel model based iterative reconstruction (MBIR) algorithm that uses a combination of variable splitting technique and an implicit preconditioner method to speed up the convergence of simultaneous algorithms for MBIR.
  • the variable splitting technique creates a sequence of sub-problems, where in each sub-problem the Hessian matrix may be approximately shift-invariant.
  • Preconditioners may be used to accelerate the convergence of the sub-problem.
  • the implicit preconditioner method is used to create an approximation to the sub-problem in each of a plurality of outer loop iterations.
  • the implicit preconditioner method works by iteratively solving the equation that relates image and measurement data without explicitly calculating the inverse of the Hessian.
  • the approximated sub-problem is then solved by inner loop iterations that have relatively low computational cost.
  • the inner loop calculation improves the convergence speed of the outer loop iteration.
  • the optimization algorithm for the MBIR objective function is composed of nested iteration loops, wherein the inner loop is formulated to iteratively optimize an approximate objective function defined to approximate the cost function iteratively optimized by the outer loop; and the outer loop iterations optimize the original MBIR cost function using the solution generated using the inner loop.
  • At least one technical effect of various embodiments is that the outer and inner loop calculations converge to the solution of the original cost function in fewer iterations than conventional simultaneous algorithms.
  • various embodiments described herein may converge to acceptable image quality in five to ten iterations as compared to 50 or 100 iterations required by conventional algorithms published in the art.
  • FIG. 1 is a flowchart of an exemplary method 100 for reconstructing an image of an object.
  • the method includes acquiring projection, or measurement, data using an exemplary imaging system.
  • an objective function is formulated.
  • the objective function is formulated at 104 based on the measurement data and the prior knowledge of the image representing the scanned object. For example, let x denote the image and y denote the measurement data, wherein both x and y represent vectors.
  • ⁇ (x) is the objective function, in which, J(x, y) is the data mismatch term that penalizes the inconsistency between the image and the measurement; and ⁇ (x) is the regularization function that penalizes differences between adjacent elements in the image.
  • A is the forward system matrix
  • W is a diagonal weighting matrix.
  • the i th diagonal entry of the matrix W, denoted by w i is inversely proportional to an estimate of the variance in the measurement y i .
  • other forms of data mismatch terms may be used with the method 100 , such as, but not limited to, Poisson log likelihood function and/or the like.
  • the optimization problem could have additional equality or inequality constraints on the variables. For example, non-negativity constraints are commonly applied to image pixels.
  • the method 100 includes transforming the original problem (i.e., Equation 1) to a sequence of sub-problems.
  • the sub-problems are created in two steps.
  • a set of auxiliary variables is introduced to add additional terms in the objective functions and additional constraints to the optimization problem.
  • the introduction at 106 of the set of auxiliary variables transforms the original problem to a constrained optimization problem.
  • the constrained optimization problem is then solved, at 108 , using a method of multipliers.
  • the method of multipliers turns the constrained optimization problem into a sequence of sub-problems (e.g., the steps 110 - 112 described below), where each sub-problem is relatively easier to solve than the original problem.
  • Equation 3 the Hessian matrix of the data mismatch term is given by A t WA. Due to the relatively large dynamic range in W, it may be difficult to design a preconditioner P that that approximates (A t WA) ⁇ 1 relatively well. On the other hand, the matrix A t A for tomographic reconstruction problems can be well approximated by shift-invariant image space operators, which are relatively easy to invert. Therefore, the method 100 is used to solve the original problem through a sequence of sub-problems, while in the sub-problem it is only needed to invert the A t A matrix.
  • the auxiliary variable u corresponds to the forward projection of the image.
  • the constrained optimization problem is solved using a method of multipliers.
  • the method of multipliers turns the constrained optimization problem into a sequence of sub-problems, and then solves the sequence of sub-problems iteratively.
  • solving the constrained optimization problem using a method of multipliers includes solving the original problem by solving a the sequence of sub-problems.
  • the constrained optimization problem can be solved using various methods.
  • solving at 108 the constrained optimization problem includes solving the constrained optimization problem using an augmented Lagrangian (AL) method.
  • AL Lagrangian
  • Other methods may be used in addition or alternative to the AL method.
  • FIG. 2 is another flowchart illustrating an exemplary method of solving at 108 the constrained optimization problem using a method of multipliers.
  • the sequence of sub-problems created at 108 may include, for example, the reconstruction of the image x at 110 a , the auxiliary variable u (which corresponds to the forward projection of x) at 110 b , and a transform of the Lagrangian multiplier ⁇ at 112 .
  • the reconstruction of the image x may be initialized with filtered back projection (FBP) reconstruction.
  • the auxiliary variable u may be initialized by forward projecting the image x, while ⁇ may be initialized as zeros.
  • the method of solving at 108 includes updating the sub-problems of the variables x, u, and ⁇ alternatively.
  • the updates of the sinogram variables ⁇ and u have closed form solutions, while the update of x is solved iteratively.
  • the method of solving at 108 includes updating the reconstruction of the image x and the auxiliary variable u using the initial value of ⁇ .
  • the reconstruction of the image x and the auxiliary variable u are updated at 110 alternatively.
  • the reconstruction of the image x is updated at 110 a before the auxiliary variable u is updated at 110 b .
  • the reconstruction of the image x and the auxiliary variable u may be updated at 110 in any order.
  • the auxiliary variable u may be updated before the reconstruction of the image x is updated. Examples of updating the reconstruction of the image x at 110 a and updating the auxiliary variable u at 110 b include, but are not limited to, using the Equations 13 and 12 (described below), respectively.
  • the transform of the Lagrangian multiplier ⁇ is then updated, at 112 , using the updated reconstruction of the image x and the updated auxiliary variable u.
  • An example of updating the transform of the Lagrangian multiplier ⁇ at 112 includes, but is not limited to, using Equation 9 (described below).
  • the method of solving at 108 includes iteratively updating the reconstruction of the image x using the auxiliary variable u and the updated transform of the Lagrangian multiplier ⁇ .
  • the reconstruction of the image x may be iteratively updated at 114 until the objective function ⁇ (x) is reduced to a predetermined level and/or minimized.
  • Each iteration of the AL method is composed of two steps.
  • the update of x e.g., the step 110 a of FIG. 2
  • the update of u e.g., the step 110 b of FIG. 2
  • the joint optimization problem is solved by updating x and u alternatively using:
  • x ( j + 1 ) arg ⁇ ⁇ min x ⁇ L ⁇ ( x , u ( j ) , ⁇ ( j ) ) Equation ⁇ ⁇ 10
  • u ( j + 1 ) arg ⁇ ⁇ min u ⁇ L ⁇ ( x ( j + 1 ) , u , ⁇ ( j ) ) Equation ⁇ ⁇ 11
  • Equation 10 Equation 10
  • x ( j + 1 ) arg ⁇ ⁇ min x ⁇ ⁇ 2 ⁇ ⁇ Ax - u ( j + 1 ) + ⁇ ( j + 1 ) ⁇ ⁇ 2 + ⁇ ⁇ ( x ) Equation ⁇ ⁇ 13
  • FIG. 3 is another flowchart illustrating an exemplary method of iteratively updating at 114 the reconstruction of the image x using the updated forward projection of the image u and the updated transform of the Lagrangian multiplier ⁇ .
  • an implicit preconditioner method is used to create an approximation to the sub-problem in each of a plurality of outer loop iterations.
  • the method 114 includes an outer loop that iteratively minimizes the original cost function.
  • an approximate cost function is first formulated. The approximate cost function is then optimized using inner loop iterations.
  • the method 114 includes nested iteration loops, wherein the inner loop is formulated to iteratively optimize an approximate objective function defined to approximate the cost function iteratively optimized by the outer loop.
  • the result of the inner loop is used as preconditioned gradient direction to update the image.
  • the outer loop iterations optimize the original cost function using the solution generated using the inner loop.
  • the method 114 includes various steps 116 - 120 that define the outer loop. Specifically, at 116 , the method 114 includes formulating an approximate objective function. At 118 , the approximate objective function is optimized using inner loop iterations shown in FIG. 4 .
  • the step 118 may be referred to herein as the “inner loop”.
  • the inner loop 118 of the algorithm is utilized to increase convergence speed of the algorithm, or to decrease computational complexity for fast implementation on computer hardware.
  • the inner loop 118 includes steps 124 - 132 (described below and illustrated in FIG. 4 ) that are used to solve the inner loop problem, which is then utilized to solve the problem in the outer loop.
  • the solution to the approximate objective function minimized by the inner loop 118 is used by the outer loop during each iteration.
  • an update is computed using the solution of the inner loop 118 .
  • the image x is updated using the update computed at 120 .
  • the steps 116 - 122 are performed iteratively until the cost function is reduced to a predetermined level and/or minimized.
  • FIG. 4 is another flowchart illustrating an exemplary method of optimizing at 118 that approximate objective function of the outer loop.
  • the inner loop is initialized.
  • a local approximate function is formulated at 126 .
  • the local approximate function may be solved directly or iteratively. Solving at 128 the local approximate function may include computing an initial search direction. Moreover solving the local approximate function at 128 may further include computing an improved search direction. An optimal step size may be computed along the new improved search direction.
  • the method 118 further includes computing an update using the solution to the local approximate function solved at 128 .
  • the update is computed at 130 using the improved search direction and/or the optimal step size.
  • the updates from steps 124 - 130 are applied to the image x.
  • the output from the inner loop calculations are then used by the outer loop. Specifically, and referring again to FIG. 3 , after the input has been received from the inner loop calculations, an optimized approximate objective function is computed at 118 by solving the inner loop problem. A new search direction may be computed for the outer loop.
  • the update is computed using the solution of the inner loop 118 .
  • the image x is updated using the update computed at 120 .
  • the steps 116 - 122 may be performed iteratively until the cost function is reduced to a predetermined level and/or minimized.
  • the problem of solving the update of x in general does not have a closed form solution, and therefore needs to be solved iteratively.
  • the method does not require solving the optimization problems fully, therefore only a relatively small number of iterations (e.g., less than 20, equal to or less than 10, and/or the like) needs to be run.
  • the sub-problem of updating x has similar structure to that of the original problem (i.e., Equation 1).
  • the sub-problem of updating x is composed of a data mismatch term and a regularization term.
  • the data mismatch term is replaced by u (j+1) ⁇ (j+1) , and the statistical weighting matrix W is replaced by ⁇ .
  • ⁇ may be chosen to make A t ⁇ A an approximately shift-invariant operator.
  • ⁇ may also be chosen to improve the conditioning of the A t ⁇ A matrix.
  • ⁇ is chosen to be identity matrix to simplify the sub-problem of updating x. Therefore the Hessian of the data mismatch term is A t ⁇ A instead of A t WA.
  • the sub-problem of updating x may be easier to solve because: (1) A t ⁇ A does not depend on the data and therefore it is possible to precompute and store an approximation of A t ⁇ A or (A t ⁇ A) ⁇ 1 for each scan protocol; and (2) A t ⁇ A can be approximately modeled as shift invariant and therefore it may be easier to design shift invariant preconditioners for the sub-problem of updating x.
  • Equation 13 In the sub-problem of updating x using Equation 13, if we choose ⁇ to be the identity matrix, the eigenvalues of A t A associating with lower frequencies are relatively large, while the eigenvalues associating with higher frequencies are relatively small. Therefore, conventional algorithms such as conjugate gradients (CG) or gradient descent (GD) algorithms will suffer from slower convergence of higher frequencies.
  • An explicitly defined preconditioner matrix P that approximates (A t ⁇ A) ⁇ 1 can be applied.
  • a relatively good preconditioner design may satisfy two goals: (1) P must be relatively easy to apply; and (2) P needs to be a relatively close approximation to (A t ⁇ A) ⁇ 1 .
  • such a preconditioner can be designed by approximating A t A as a shift invariant operator. Therefore, P ⁇ (A t ⁇ A) ⁇ 1 can be computed using only a relatively few fast Fourier transforms (FFTs) and image sealing. Alternatively, one can select P to be an approximation to the inverse of the Hessian of the entire inner objective function in Equation 13.
  • FFTs fast Fourier transforms
  • One drawback of explicit preconditioner design is that it may be difficult to include the Hessian of the regularization term into the design, for example because the Hessian of the regularization term is not shift-invariant in general and is a function of x when ⁇ (x) is non-quadratic.
  • One way to address such a problem is to split the regularization term and data mismatch term by introducing another set of auxiliary variables.
  • an auxiliary variable v that corresponds to the image x or to a function of the image x is introduced to further split the regularization term in the sub-problem of updating x.
  • the method 114 solves the sub-problem of updating x using an implicitly defined preconditioner.
  • An implicit preconditioner differs from an explicit preconditioner in both design and application.
  • a matrix M that approximates the Hessian matrix H is found.
  • P is designed to directly to approximate H ⁇ 1 .
  • approximating H instead of H ⁇ 1 has at least two benefits.
  • approximating H instead of H ⁇ 1 may provide the flexibility to use a more sophisticated and/or more accurate model of the Hessian.
  • the Hessian matrix changes every iteration.
  • H the Hessian matrix
  • the implicit and explicit preconditioner also differ in the way the implicit and explicit preconditioner are applied.
  • the preconditioned gradient direction d (n) is computed by solving an approximate cost function using iterative loops. Therefore, the algorithm using the implicit preconditioner is composed of nested iterative loops.
  • the outer loop iterations can be a preconditioned gradient-based algorithm, in which the preconditioned gradient direction is supplied by the inner loop iterations. In each inner loop, we minimize an approximate cost function.
  • x (n) is the image value at n th outer loop iteration
  • ⁇ (n) A t ⁇ A(Ax (n) ⁇ u (j) +n (j) )
  • c is a constant.
  • a t A is approximated with a shift-invarient operator K.
  • K can be constructed by sampling the A t A operator at the iso center. Notice that the proposed operator is a pure image space operator; and it is relatively easy to compute because it only requires filtering the image using FFTs.
  • ⁇ is not identity matrix
  • the Fessler and Rogers approximation to A t ⁇ A operator is used according to: A t ⁇ A ⁇ DA t AD Equation 16 where D is a diagonal matrix, with the i th diagonal element given by
  • d i ⁇ j ⁇ a ij 2 ⁇ ⁇ j ⁇ j ⁇ a ij 2 .
  • a t A is not assumed to be shift-invariant. Instead, an approximation to A t A is found as a product of sparse matrices.
  • M 1 One choice for M 1 is to precompute and store the A t A matrix.
  • a t A is not sparse, storing A t A may be unpractical.
  • a method has been proposed by Cao, Bouman and Webb to compute a sparse matrix representation for a non-sparse operator.
  • the sparsified matrices are relatively small enough to be stored and are relatively easy to apply.
  • a method has been proposed by Wei to compute relatively fast space-varying convolutions that could potentially be used to approximate the A t A operator.
  • the calculation of A itself can be simplified by using a simpler forward model, for example denoted by ⁇ .
  • simpler forward models include, but are not limited to, pixel driven forward and back projections, ray driven forward and back projections, and/or the like.
  • the computation time of iterative reconstruction depends on the accuracy of the system model being used. More specifically, an accurate system model that incorporates the exact 3D scanning geometry, the effect of finite focal spot size, detector blurring, or other physical effects, increases computation time. Accordingly, various methods of simplification and approximations may be utilized to reduce computation time.
  • Siddon's method may be used to efficiently compute line integrals through a 3D pixel array with relatively compact computer code.
  • other simplified interpolation schemes may be used to make forward and back projection operators more computationally efficient.
  • ⁇ can be computed using a lower-resolution image grid of detectors.
  • ⁇ may be constructed by ignoring certain physical effects (positron range, crystal penetration, etc.), using a numerically more sparse system matrix, or configured based on texture mapping hardware on the GPU.
  • CT reconstruction is that the outer loop may compute ⁇ using an accurate model, such as the Distance Driven model, while the inner loop is computed using Pixel Driven or Ray Driven models with coarse sampling in focal spot and detector cells.
  • the outer loop may compute ⁇ using exact 3D cone beam geometry, while the inner loop ⁇ is computed using a parallel or semi-parallel beam geometry using rebinned or re-sliced sinogram.
  • ⁇ may also be computed using data types with lower precision and faster computation than ⁇ .
  • ⁇ may also be computed using a Fast Fourier Transform (FFT) based technique.
  • FFT Fast Fourier Transform
  • “fast” refers to faster computation than for the full original model A.
  • “fast” can refer to a computational complexity on the order of N 2 log(N) instead of N 3 in 2D reconstruction or N 3 log(N) instead of N 4 in 3D reconstruction.
  • ⁇ may also be computed with lower resolution, but higher-order image basis functions from A. For example, instead of using rectangular voxels, ⁇ may be computed using a piece-wise linear or higher-order parametric representation of voxels (for example, a z-slope voxel model or blobs).
  • FIG. 5 is a visual representation of exemplary difference images.
  • FIG. 5 illustrates preliminary results of one example of the method 100 .
  • the sub-problems were created according to Equation 10.
  • the sub-problem of updating the image x was solved using the implicit preconditioner method.
  • the implicit preconditioner method twenty iterations of a preconditioned conjugate gradient algorithm were used to solve the inner loop problem. Wires and resolution bars were used to test the spatial resolution of the reconstruction, and uniform region to measure the noise.
  • FIG. 5 compares the results of the example of the method 100 with a fully converged reference image using an ICD algorithm.
  • the difference images of FIG. 5 are shown in a window between ⁇ 100 to 100 Hounsfield units.
  • FIG. 5 illustrates the difference image between the FBP image and the reference, in which the differences near the wire and the uniform region can be seen, which illustrates the difference in spatial resolution and noise, respectively.
  • FIG. 5( b ) shows the difference between the image of three iterations of the example of the method 100 and the reference, in which the differences are significantly reduced. After ten iterations, as shown in FIG. 5( c ), the difference diminishes to the point that the difference is barely noticeable. The results qualitatively show that the method 100 converges to the solution within ten iterations.
  • the methods and algorithms described herein are used to reconstruct an image of an object.
  • the methods and algorithms may be embodied as a set of instructions that are stored on a computer and implemented using, for example, a module 530 , shown in FIG. 6 , software, hardware, a combination thereof, and/or a tangible non-transitory computer readable medium.
  • a tangible non-transitory computer readable medium excludes signals.
  • FIG. 6 is a pictorial view of an exemplary imaging system 500 that is formed in accordance with various embodiments.
  • FIG. 7 is a block schematic diagram of a portion of the multi-modality imaging system 500 shown in FIG. 6 .
  • the imaging system may be embodied as a computed tomography (CT) imaging system, a positron emission tomography (PET) imaging system, a magnetic resonance imaging (MRI) system, an ultrasound imaging system, an x-ray imaging system, a single photon emission computed tomography (SPECT) imaging system, an interventional C-Arm tomography imaging system, a CT system for a dedicated purpose such as extremity or breast scanning, and combinations thereof, among others.
  • CT computed tomography
  • PET positron emission tomography
  • MRI magnetic resonance imaging
  • SPECT single photon emission computed tomography
  • interventional C-Arm tomography imaging system a CT system for a dedicated purpose such as extremity or breast scanning, and combinations thereof, among others.
  • CT computed tomography
  • PET positron emission tomography
  • the multi-modality imaging system 500 is illustrated, and includes a CT imaging system 502 and a PET imaging system 504 .
  • the imaging system 500 allows for multiple scans in different modalities to facilitate an increased diagnostic capability over single modality systems.
  • the exemplary multi-modality imaging system 500 is a CT/PET imaging system 500 .
  • modalities other than CT and PET are employed with the imaging system 500 .
  • the imaging system 500 may be a standalone CT imaging system, a standalone PET imaging system, a magnetic resonance imaging (MRI) system, an ultrasound imaging system, an x-ray imaging system, and/or a single photon emission computed tomography (SPECT) imaging system, interventional C-Arm tomography, CT systems for a dedicated purpose such as extremity or breast scanning, and combinations thereof, among others.
  • MRI magnetic resonance imaging
  • SPECT single photon emission computed tomography
  • interventional C-Arm tomography CT systems for a dedicated purpose such as extremity or breast scanning, and combinations thereof, among others.
  • the CT imaging system 502 includes a gantry 510 that has an x-ray source 512 that projects a beam of x-rays toward a detector array 514 on the opposite side of the gantry 510 .
  • the detector array 514 includes a plurality of detector elements 516 that are arranged in rows and channels that together sense the projected x-rays that pass through an object, such as the subject 506 .
  • the imaging system 500 also includes a computer 520 that receives the projection data from the detector array 514 and processes the projection data to reconstruct an image of the subject 506 . In operation, operator supplied commands and parameters are used by the computer 520 to provide control signals and information to reposition a motorized table 522 .
  • the motorized table 522 is utilized to move the subject 506 into and out of the gantry 510 .
  • the table 522 moves at least a portion of the subject 506 through a gantry opening 524 that extends through the gantry 510 .
  • the imaging system 500 also includes a module 530 that is configured to implement various methods and algorithms described herein.
  • the module 530 may be implemented as a piece of hardware that is installed in the computer 520 .
  • the module 530 may be implemented as a set of instructions that are installed on the computer 520 .
  • the set of instructions may be stand alone programs, may be incorporated as subroutines in an operating system installed on the computer 520 , may be functions in an installed software package on the computer 520 , and the like. It should be understood that the various embodiments are not limited to the arrangements and instrumentality shown in the drawings.
  • the detector 514 includes a plurality of detector elements 516 .
  • Each detector element 516 produces an electrical signal, or output, that represents the intensity of an impinging x-ray beam and hence allows estimation of the attenuation of the beam as it passes through the subject 506 .
  • the gantry 510 and the components mounted thereon rotate about a center of rotation 540 .
  • FIG. 7 shows only a single row of detector elements 516 (i.e., a detector row).
  • the multislice detector array 514 includes a plurality of parallel detector rows of detector elements 516 such that projection data corresponding to a plurality of slices can be acquired simultaneously during a scan.
  • the control mechanism 542 includes an x-ray controller 544 that provides power and timing signals to the x-ray source 512 and a gantry motor controller 546 that controls the rotational speed and position of the gantry 510 .
  • a data acquisition system (DAS) 548 in the control mechanism 542 samples analog data from detector elements 516 and converts the data to digital signals for subsequent processing. For example, the subsequent processing may include utilizing the module 530 to implement the various methods described herein.
  • An image reconstructor 550 receives the sampled and digitized x-ray data from the DAS 548 and performs high-speed image reconstruction.
  • the reconstructed images are input to the computer 520 that stores the image in a storage device 552 .
  • the computer 520 may receive the sampled and digitized x-ray data from the DAS 548 and perform various methods described herein using the module 530 .
  • the computer 520 also receives commands and scanning parameters from an operator via a console 560 that has a keyboard.
  • An associated visual display unit 562 allows the operator to observe the reconstructed image and other data from computer.
  • the operator supplied commands and parameters are used by the computer 520 to provide control signals and information to the DAS 548 , the x-ray controller 544 and the gantry motor controller 546 .
  • the computer 520 operates a table motor controller 564 that controls the motorized table 522 to position the subject 506 in the gantry 510 .
  • the table 522 moves at least a portion of the subject 506 through the gantry opening 524 as shown in FIG. 6 .
  • the computer 520 includes a device 570 , for example, a floppy disk drive, CD-ROM drive, DVD drive, magnetic optical disk (MOD) device, or any other digital device including a network connecting device such as an Ethernet device for reading instructions and/or data from a computer-readable medium 572 , such as a floppy disk, a CD-ROM, a DVD or an other digital source such as a network or the Internet, as well as yet to be developed digital means.
  • the computer 520 executes instructions stored in firmware (not shown).
  • the computer 520 is programmed to perform functions described herein, and as used herein, the term computer is not limited to just those integrated circuits referred to in the art as computers, but broadly refers to computers, processors, microcontrollers, microcomputers, programmable logic controllers, application specific integrated circuits, and other programmable circuits, and these terms are used interchangeably herein.
  • the x-ray source 512 and the detector array 514 are rotated with the gantry 510 within the imaging plane and around the subject 506 to be imaged such that the angle at which an x-ray beam 574 intersects the subject 506 constantly changes.
  • a group of x-ray attenuation measurements, i.e., projection data, from the detector array 514 at one gantry angle is referred to as a “view”.
  • a “scan” of the subject 506 comprises a set of views made at different gantry angles, or view angles, during one revolution of the x-ray source 512 and the detector 514 .
  • the projection data is processed to reconstruct an image that corresponds to a two dimensional slice taken through the subject 506 .
  • multi-modality imaging system Exemplary embodiments of a multi-modality imaging system are described above in detail.
  • the multi-modality imaging system components illustrated are not limited to the specific embodiments described herein, but rather, components of each multi-modality imaging system may be utilized independently and separately from other components described herein.
  • the multi-modality imaging system components described above may also be used in combination with other imaging systems.
  • the various embodiments may be implemented in hardware, software or a combination thereof.
  • the various embodiments and/or components also may be implemented as part of one or more computers or processors.
  • the computer or processor may include a computing device, an input device, a display unit and an interface, for example, for accessing the Internet.
  • the computer or processor may include a microprocessor.
  • the microprocessor may be connected to a communication bus.
  • the computer or processor may also include a memory.
  • the memory may include Random Access Memory (RAM) and Read Only Memory (ROM).
  • the computer or processor further may include a storage device, which may be a hard disk drive or a removable storage drive such as a solid state drive, optical drive, and/or the like.
  • the storage device may also be other similar means for loading computer programs or other instructions into the computer or processor.
  • the term “computer” may include any processor-based or microprocessor-based system including systems using microcontrollers, reduced instruction set computers (RISC), application specific integrated circuits (ASICs), logic circuits, GPUs, FPGAs, and any other circuit or processor capable of executing the functions described herein.
  • RISC reduced instruction set computers
  • ASICs application specific integrated circuits
  • the above examples are exemplary only, and are thus not intended to limit in any way the definition and/or meaning of the term “computer”.
  • the computer or processor executes a set of instructions that are stored in one or more storage elements, in order to process input data.
  • the storage elements may also store data or other information as desired or needed.
  • the storage element may be in the form of an information source or a physical memory element within a processing machine.
  • the set of instructions may include various commands that instruct the computer or processor as a processing machine to perform specific operations such as the methods and processes of the various embodiments of the invention.
  • the set of instructions may be in the form of a software program.
  • the software may be in various forms such as system software or application software, which may be a non-transitory computer readable medium. Further, the software may be in the form of a collection of separate programs, a program module within a larger program or a portion of a program module.
  • the software also may include modular programming in the form of object-oriented programming.
  • the processing of input data by the processing machine may be in response to user commands, or in response to results of previous processing, or in response to a request made by another processing machine.
  • the phrase “reconstructing an image” is not intended to exclude embodiments of the present invention in which data representing an image is generated, but a viewable image is not. Therefore, as used herein the term “image” broadly refers to both viewable images and data representing a viewable image. However, many embodiments generate, or are configured to generate, at least one viewable image.
  • the terms “software” and “firmware” are interchangeable, and include any computer program stored in memory for execution by a computer, including RAM memory, ROM memory, EPROM memory, EEPROM memory, and non-volatile RAM (NVRAM) memory.
  • RAM memory random access memory
  • ROM memory read-only memory
  • EPROM memory erasable programmable read-only memory
  • EEPROM memory electrically erasable programmable read-only memory
  • NVRAM non-volatile RAM

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (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)
  • Image Analysis (AREA)
US13/531,082 2012-06-22 2012-06-22 Method and apparatus for iterative reconstruction Active 2032-12-15 US8885975B2 (en)

Priority Applications (5)

Application Number Priority Date Filing Date Title
US13/531,082 US8885975B2 (en) 2012-06-22 2012-06-22 Method and apparatus for iterative reconstruction
JP2013126219A JP6280700B2 (ja) 2012-06-22 2013-06-17 反復式再構成の方法、非一時的なコンピュータ可読の媒体及びイメージング・システム
NL2011005A NL2011005B1 (en) 2012-06-22 2013-06-19 Method and apparatus for iterative reconstruction.
DE102013106467.1A DE102013106467A1 (de) 2012-06-22 2013-06-20 Verfahren und Vorrichtung zur iterativen Rekonstruktion
CN201310249820.XA CN103514629B (zh) 2012-06-22 2013-06-21 用于迭代重建的方法和设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
US13/531,082 US8885975B2 (en) 2012-06-22 2012-06-22 Method and apparatus for iterative reconstruction

Publications (2)

Publication Number Publication Date
US20130343672A1 US20130343672A1 (en) 2013-12-26
US8885975B2 true US8885975B2 (en) 2014-11-11

Family

ID=49517577

Family Applications (1)

Application Number Title Priority Date Filing Date
US13/531,082 Active 2032-12-15 US8885975B2 (en) 2012-06-22 2012-06-22 Method and apparatus for iterative reconstruction

Country Status (5)

Country Link
US (1) US8885975B2 (zh)
JP (1) JP6280700B2 (zh)
CN (1) CN103514629B (zh)
DE (1) DE102013106467A1 (zh)
NL (1) NL2011005B1 (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140294147A1 (en) * 2013-03-15 2014-10-02 Varian Medical Systems, Inc. Systems and methods for multi-view imaging and tomography
US9824468B2 (en) 2015-09-29 2017-11-21 General Electric Company Dictionary learning based image reconstruction
JP2017209502A (ja) * 2016-05-27 2017-11-30 東芝メディカルシステムズ株式会社 Ｘ線コンピュータ断層撮影装置及び医用画像処理装置
US20190209015A1 (en) * 2017-12-18 2019-07-11 Arizona Board Of Regents On Behalf Of The University Of Arizona Multi-field miniaturized micro-endoscope

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10062135B2 (en) 2013-07-31 2018-08-28 National Technology & Engineering Solutions Of Sandia, Llc Graphics processing unit management system for computed tomography
US9495770B2 (en) 2013-08-14 2016-11-15 University Of Utah Research Foundation Practical model based CT construction
US9524567B1 (en) * 2014-06-22 2016-12-20 InstaRecon Method and system for iterative computed tomography reconstruction
US9761019B2 (en) * 2015-01-08 2017-09-12 Toshiba Medical Systems Corporation Computed tomography using simultaneous image reconstruction with measurements having multiple distinct system matrices
CN105212956B (zh) * 2015-08-25 2018-03-16 浙江大学 一种基于ist的次晶体级pet系统时间修正方法
CN105869075A (zh) * 2016-04-19 2016-08-17 东南大学 一种冷热电联供型微型能源网经济优化调度方法
CN107784684B (zh) * 2016-08-24 2021-05-25 深圳先进技术研究院 一种锥束ct三维重建方法及系统
US10217269B1 (en) * 2017-11-29 2019-02-26 Siemens Healthcare Gmbh Compressive sensing of light transport matrix
CN108416723B (zh) * 2018-02-07 2022-02-18 南京理工大学 基于全变分正则化和变量分裂的无透镜成像快速重构方法
US10902650B2 (en) * 2018-04-19 2021-01-26 Fei Company X-ray beam-hardening correction in tomographic reconstruction using the Alvarez-Macovski attenuation model
WO2020009752A1 (en) * 2018-07-02 2020-01-09 Exxonmobil Upstream Research Company Full wavefield inversion with an image-gather-flatness constraint
CN109614574A (zh) * 2018-11-23 2019-04-12 成都景中教育软件有限公司 一种动态几何软件中迭代的实现方法
EP3660790A1 (en) 2018-11-28 2020-06-03 Koninklijke Philips N.V. System for reconstructing an image of an object
CN110455719A (zh) * 2019-08-16 2019-11-15 中国科学技术大学 三维光声成像系统及方法
US11335038B2 (en) * 2019-11-04 2022-05-17 Uih America, Inc. System and method for computed tomographic imaging
CN113538649B (zh) * 2021-07-14 2022-09-16 深圳信息职业技术学院 一种超分辨率的三维纹理重建方法、装置及其设备
CN113641956B (zh) * 2021-08-05 2023-05-30 中国科学院软件研究所 面向SW26010-Pro处理器的1、2级BLAS函数库的高性能实现方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8310233B2 (en) * 2009-02-18 2012-11-13 Mayo Foundation For Medical Education And Research Method for image reconstruction from undersampled medical imaging data

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6980681B1 (en) * 2000-04-24 2005-12-27 Ge Medical Systems Global Technology Company, Llc Methods and apparatus for helical reconstruction for multislice CT scan
US8897528B2 (en) * 2006-06-26 2014-11-25 General Electric Company System and method for iterative image reconstruction
US8571287B2 (en) * 2006-06-26 2013-10-29 General Electric Company System and method for iterative image reconstruction
EP2149284B1 (en) * 2007-05-31 2014-08-13 General Electric Company Methods and systems to facilitate correcting gain fluctuations in image
US8660330B2 (en) * 2008-06-27 2014-02-25 Wolfram Jarisch High efficiency computed tomography with optimized recursions
US8761478B2 (en) * 2009-12-15 2014-06-24 General Electric Company System and method for tomographic data acquisition and image reconstruction
US20110164031A1 (en) * 2010-01-06 2011-07-07 Kabushiki Kaisha Toshiba Novel implementation of total variation (tv) minimization iterative reconstruction algorithm suitable for parallel computation
CN101794440B (zh) * 2010-03-12 2012-04-18 东南大学 图像序列的加权自适应超分辨率重建方法
CN101923056B (zh) * 2010-07-20 2012-10-03 南昌航空大学 一种非完全数据下的光学层析重建计算方法
CN102163338B (zh) * 2011-04-08 2014-09-03 哈尔滨工业大学 一种压缩感知系统中的高效重建方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8310233B2 (en) * 2009-02-18 2012-11-13 Mayo Foundation For Medical Education And Research Method for image reconstruction from undersampled medical imaging data

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140294147A1 (en) * 2013-03-15 2014-10-02 Varian Medical Systems, Inc. Systems and methods for multi-view imaging and tomography
US9778391B2 (en) * 2013-03-15 2017-10-03 Varex Imaging Corporation Systems and methods for multi-view imaging and tomography
US9824468B2 (en) 2015-09-29 2017-11-21 General Electric Company Dictionary learning based image reconstruction
JP2017209502A (ja) * 2016-05-27 2017-11-30 東芝メディカルシステムズ株式会社 Ｘ線コンピュータ断層撮影装置及び医用画像処理装置
US20190209015A1 (en) * 2017-12-18 2019-07-11 Arizona Board Of Regents On Behalf Of The University Of Arizona Multi-field miniaturized micro-endoscope

Also Published As

Publication number Publication date
CN103514629A (zh) 2014-01-15
NL2011005B1 (en) 2015-12-15
US20130343672A1 (en) 2013-12-26
JP2014004358A (ja) 2014-01-16
DE102013106467A1 (de) 2014-05-22
NL2011005A (en) 2013-12-24
JP6280700B2 (ja) 2018-02-14
CN103514629B (zh) 2017-11-07

Similar Documents

Publication Publication Date Title
US8885975B2 (en) Method and apparatus for iterative reconstruction
US9020230B2 (en) Method and apparatus for iterative reconstruction
US8958660B2 (en) Method and apparatus for iterative reconstruction
US8135186B2 (en) Method and system for image reconstruction
US7885371B2 (en) Method and system for image reconstruction
US6907102B1 (en) Iterative reconstruction methods for multi-slice computed tomography
EP2410491B1 (en) System and method for reconstruction of X-ray images
US7706497B2 (en) Methods and apparatus for noise estimation for multi-resolution anisotropic diffusion filtering
US8897528B2 (en) System and method for iterative image reconstruction
US8885903B2 (en) Method and apparatus for statistical iterative reconstruction
US9600924B2 (en) Iterative reconstruction of image data in CT
US9524567B1 (en) Method and system for iterative computed tomography reconstruction
US20060215891A1 (en) Method and system for controlling image reconstruction
EP1478274A2 (en) Methods and apparatus for divergent fast beam tomography
US8687869B2 (en) System and method for acceleration of image reconstruction
US20160350944A1 (en) Systems and methods for parallel processing of imaging information
Park et al. A fully GPU-based ray-driven backprojector via a ray-culling scheme with voxel-level parallelization for cone-beam CT reconstruction
US9965875B2 (en) Virtual projection image method
US8379948B2 (en) Methods and systems for fast iterative reconstruction using separable system models
US10307114B1 (en) Iterative volume image reconstruction using synthetic projection images
US20190206094A1 (en) Fast iterative image reconstruction method for emission tomography
Amrehn et al. Portability of TV-regularized reconstruction parameters to varying data sets

Legal Events

Date Code Title Description
AS Assignment

Owner name: GENERAL ELECTRIC COMPANY, NEW YORK

Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNORS:YU, ZHOU;DE MAN, BRUNO KRISTIAAN BERNARD;THIBAULT, JEAN-BAPTISTE;AND OTHERS;REEL/FRAME:030236/0196

Effective date: 20130416

AS Assignment

Owner name: THE REGENTS OF THE UNIVERSITY OF MICHIGAN, MICHIGA

Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNORS:FESSLER, JEFFREY A.;RAMANI, SATHISH;REEL/FRAME:030265/0148

Effective date: 20120418

AS Assignment

Owner name: NATIONAL INSTITUTES OF HEALTH (NIH), U.S. DEPT. OF

Free format text: CONFIRMATORY LICENSE;ASSIGNOR:UNIVERSITY OF MICHIGAN;REEL/FRAME:030304/0268

Effective date: 20130426

AS Assignment

Owner name: UNIVERSITY OF NOTRE DAME DU LAC, INDIANA

Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNOR:SAUER, KEN;REEL/FRAME:030397/0471

Effective date: 20130501

AS Assignment

Owner name: PURDUE RESEARCH FOUNDATION, INDIANA

Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNOR:BOUMAN, CHARLES A.;REEL/FRAME:030479/0936

Effective date: 20130524

STCF Information on status: patent grant

Free format text: PATENTED CASE

MAFP Maintenance fee payment

Free format text: PAYMENT OF MAINTENANCE FEE, 4TH YEAR, LARGE ENTITY (ORIGINAL EVENT CODE: M1551)

Year of fee payment: 4

MAFP Maintenance fee payment

Free format text: PAYMENT OF MAINTENANCE FEE, 8TH YEAR, LARGE ENTITY (ORIGINAL EVENT CODE: M1552); ENTITY STATUS OF PATENT OWNER: LARGE ENTITY

Year of fee payment: 8