Method and device for iterative reconstruction

Abstract

A method for stepwise reconstructing an image of an object, wherein the method comprises the steps: accessing measurement data associated with the image; and using a simultaneous algorithm to reconstruct the image, wherein using the simultaneous algorithm to reconstruct the image comprises the steps: determining a voxel-dependent scaling factor; and applying the voxel-dependent scaling factor to a gradient of an objective function to reconstruct the image, wherein determining a voxel-dependent scaling factor includes the step of calculating the scaling factor directly based on an acquisition geometry.

Description

BACKGROUND OF THE INVENTIONThe invention described herein relates generally to imaging systems and, more particularly, to a method and apparatus for reconstructing an image using iterative techniques.Conventionally, images have been reconstructed from computed tomography (CT) data using immediate reconstruction algorithms such as filtered back projection (FBP) or convolutional back projection (CBP). Recently, model-based iterative reconstruction (MBIR) algorithms for reconstructing CT images have been commercially available. An advantage of MBIR algorithms is based on the fact that MBIR algorithms can more accurately simulate the measured values obtained with the aid of CT systems. This is particularly true for spiral CT systems with multilayer detectors because spiral CT systems generate projection measurements that obliquely traverse the reconstructed two-dimensional (2D) image planes, and because the acquired data inherently includes noise. By more accurately simulating these projections, MBIR algorithms can produce higher quality (e.g., resolution) reconstructions, lower interference level, and fewer artifacts. As a result, MBIR algorithms can be used as a tool to substantially reduce the dose in the case of CT scans while maintaining diagnostic image quality.However, the computational time and the computational resources required to perform a reconstruction are a relatively great problem in MBIR. MBIR algorithms usually reconstruct an image by first constructing a target function involving an accurate system model, a statistical noise model, and a conventional model. With the present objective function, the image is subsequently reconstructed by calculating an approximation value that minimizes the objective function, which is usually performed by means of an iterative optimization algorithm. Examples of some such iterative optimization algorithms include: iterative coordinated gradient descent (ICD), expected maximization (EM) variations, conjugate gradient (CG) methods, and ordered subsets (OS). However, due to the complexity of the MBIR objective function and the associated iterative solution, some iterative algorithms may require a relatively high number of iterations to achieve the final approximate value. As a result, known iterative algorithms that solve the MBIR objective function may take a lot of time to reconstruct an image.For example, Ho Lee et al., 2012, Phys. Med. Biol. 57, 2287 a method for improved reconstruction of an image in cone beam computed tomography is described. Constraints derived from previous images are used here. The article Jang Cho et al., Proc. SPIE 8313, Medical Imaging 2012: Physics of Medical Imaging, an ordered subset method for image reconstruction can be obtained, which achieves an acceleration of the method with the aid of a double surrogate function. Zhou Yu et al. in IEEE Transactions on Image Processing, Vol. 20, No. 1, Jan. 2011, introduces a fast, model-based method for X-ray CT image reconstruction with the aid of an IC optimization, wherein the optimization is spatially nonhomogeneous.BRIEF DESCRIPTION OF THE INVENTIONThe present invention relates to a method having the features of claim 1 and to an imaging system having the features of claim 13. In one embodiment, a method for gradually reconstructing an image of an object is provided. The method includes the steps of: accessing measurement data associated with the image and using a simultaneous algorithm to reconstruct the image. Using the simultaneous algorithm to reconstruct the image includes: determining a scaling factor that is voxel dependent; and applying the voxel dependent scaling factor to a gradient of a target function to reconstruct the image.In yet another embodiment, an imaging system includes a detector array and a computer coupled to the detector array. The computer is configured to access measurement data associated with an image of an object and use a simultaneous algorithm to reconstruct the image. The computer is configured to use the simultaneous algorithm for reconstruction of the image by determining a scaling factor that is voxel dependent and applying the voxel dependent scaling factor to a gradient of a target function to reconstruct the image.In yet another embodiment, a method for gradually reconstructing an image of an object is provided. The method includes the steps of: accessing measurement data associated with the image and using a simultaneous algorithm to reconstruct the image. The simultaneous algorithm includes the step of calculating a gradient of a regularization function. Using the simultaneous algorithm to reconstruct the image includes the steps of: storing a plurality of gradient samples in a look-up table; accessing the gradient samples in the look-up table; and calculating an approximated actual gradient value using the gradient samples.In yet another embodiment, a method for gradually reconstructing an image of an object is provided. The method includes the steps of: accessing measurement data associated with the image and using a simultaneous algorithm to reconstruct the image. The simultaneous algorithm includes the step of calculating a gradient of a regularization function. Using the simultaneous algorithm to reconstruct the image includes: calculating an actual gradient at an iteration of the reconstruction of the image, adding a compensation term to the actual gradient, and using the actual gradient and the compensation term as gradients for a subsequent iteration of the reconstruction of the image.BRIEF DESCRIPTION OF THE DRAWINGSFIG. 1 shows a flow diagram of an exemplary method for reconstructing an image of an object according to various embodiments. FIG. 2 pictorially illustrates an exemplary coefficient table generated using the method illustrated in FIG. 1. FIG. 3 graphically illustrates an example convergence of an example of a separable parabolic substitution function (OS SPS) algorithm using an ordered subset in conjunction with the method of FIG. 1. FIG. 4 is a pictorial representation of a plurality of images generated using the exemplary OS SPS algorithm of FIG. 3. FIG. 5 illustrates, on the basis of a further flow chart, an exemplary method for reconstructing an image of an object according to various exemplary embodiments. FIG. 6 shows a pictorial representation of an exemplary helical geometry and an exemplary reconstructed image in the helical geometry. FIG. 7 graphically illustrates exemplary mean and standard deviations of an OS SPS algorithm associated with the method of FIG. 5. FIG. 8 illustrates, on the basis of a further flow chart, an exemplary method for reconstructing an image of an object according to various exemplary embodiments. FIG. 9 illustrates, on the basis of a further flow chart, an exemplary method for reconstructing an image of an object according to various exemplary embodiments. FIG. 10 illustrates multiple images generated using an example OS SPS algorithm in conjunction with the method of FIG. 9. FIG. 11 graphically illustrates convergence rates associated with the exemplary OS SPS algorithm associated with the method of FIG. 9. FIG. 12 shows further graphs of convergence rates associated with the exemplary OS SPS algorithm associated with the method of FIG. 9. FIG. 13 illustrates multiple images generated using the example OS SPS algorithm associated with the method of FIG. 9. FIG. 14 shows an illustrative view of an exemplary multi-modal imaging system formed in accordance with various embodiments. Figure 15 is a block diagram of the system illustrated in Figure 14.DETAILED DESCRIPTION OF THE INVENTIONThe foregoing summary and the following detailed description of various embodiments will become more fully understood after reading in conjunction with the accompanying drawings. As long as the figures illustrate diagrams of the functional blocks of the various embodiments, the functional blocks do not necessarily identify the division of the hardware circuit. Thus, for example, one or more functional blocks (e.g., processors or memories) may be implemented in a single hardware component (e.g., in a general purpose signal processor or in a RAM memory / hard disk block or the like) or in multiple hardware elements. Similarly, the programs may be based on stand-alone programs, may be used as subroutines in an operating system, may be functions in an installed software package, and the like. It is understood that the various embodiments are not limited to the arrangements and functionalities shown in the figures.Various embodiments for the stepwise reconstruction of an image using a simultaneous algorithm are described herein. The image is composed of a plurality of picture elements. For purposes described herein, the term "pixel" is intended to refer to an element of an object in the image space matrix. For example, a pixel may include an image pixel or pixel that may correspond to a single voxel or volume element in a two-dimensional (2D) or three-dimensional (3D) reconstruction. The image is reconstructed by optimizing a target function. As used herein, a target function or cost function generally includes a model of the imaging system used to acquire the imaging data, a model of noise associated with the imaging system, and a model of the type of image to be reconstructed. In the following, the terms "optimize" and "minimize" or "target function" and "cost function" are interchangeable.Various embodiments provide a new approach to calculating an "adjustable" voxel-dependent scaling factor. The voxel dependent scaling factor is applied to a gradient of the simultaneous algorithm to reconstruct the image. Various embodiments provide a new approach to relaxing (or scaling) updates depending on the acquisition (e.g., scanner / reconstruction) geometry on the image. The updates for the image may be relaxed (or scaled) by either appropriately scaling the gradient, or by modifying the step size applied to the gradient to calculate the image update. The step size is usually a fixed value multiplier applied to the gradient. The calculation of the optimal step size may produce an optimal convergence. Various exemplary embodiments provide new proposals for shortening the time duration for the calculation of the gradient of a regularization function without adversely affecting the monotonic behavior of an underlying monotonic algorithm. At least one technical effect of a variety of exemplary embodiments is based on a reduction in transient processes and an improvement in the image quality given a computing time. At least one technical effect of various embodiments is based on a more rapid convergence towards solving the original cost function by means of a smaller number of iterations than with conventional simultaneous algorithms. For example, for useful image quality, different embodiments may converge to similar values of a convergence measure in five to ten iterations as compared to 50 or 100 iterations required by conventional algorithms published in the art.There is great interest in developing model-based iterative reconstruction (MBIR) algorithms for computed tomography (CT) imaging because MBIR algorithms have been found to improve image quality by reducing noise and improving resolution. MBIR algorithms usually perform their task by first forming a cost function that uses a data statistics model, a system model that describes the system geometry, and a conventional model that describes the desired characteristics of the reconstructed image. Since the cost function is usually too complex to allow immediate reconstruction of the image from the one pass measurement, the image is approximated by minimizing the cost function in multiple steps. Thus, there is a need for iterative optimization algorithms to solve such model-based cost functions and are referred to as MBIR algorithms.Let μ ∈ RN be the image and y ∈ RM be the measured projection data. The maximum posterior (MAP) approximation problem can be formulated as the sum of a logarithmized likelihood term of data L and a past term R. The logarithmized likelihood term is formulated based on statistics of the data and system geometry. The former term ensures that the image approximation values belong to some probability density function (pdf), and the β is the former scale, or former strength, with which this constraint is to be applied.Depending on the model of the data statistics, there may be several choices of logarithmized likelihood term L. One such choice of L is a quadratic form resulting in where W is an M×M diagonal matrix calculated as the approximated inverse of the variance of the measured data and A ∈ RM×N is the system geometry model. One end of the spectrum of iterative optimization algorithms is the iterative coordinated gradient descent (ICD). ICD updates a picture element (referred to as a "volume element") at a time while all other volume elements are maintained unchanged, creating a series of one-dimensional (1D) problems, with the gradient of the cost function updated after each voxel update. Calculating the gradient and Hessian Web log update of an ID cost function is relatively unproblematic, and ICD has a relatively fast convergence rate given a usable set of initial conditions. A non-homogeneous ICD (NH-ICD), a variant of an ICD, updates certain volume elements more frequently than others, resulting in more rapid convergence. Voxels are selected based on a table and the table is calculated from historical data of all updates. A non-homogeneous ICD is due to the observation that those volume elements with extensive updates are furthest from the PWLS solution (point-weighted least squares) and thus require more frequent updates. However, ICD (or modifications thereof) is difficult to parallelize relatively efficiently and relatively thoroughly on highly parallelizable computing platforms.At the other end of the spectrum of iterative optimization algorithms are simultaneous algorithms. As used herein, a "simultaneous algorithm" is an algorithm that simultaneously updates two or more volume elements. Furthermore, the term "optimization" in connection with an iterative algorithm for the reconstruction of an image is synonymous with "updating". Simultaneous algorithms for MBIR can solve an n-dimensional problem in each step. Simultaneous algorithms are available for parallelization on relatively highly parallel computing platforms. An iteration of such algorithms can be faster than ICD in terms of calculation, but simultaneous optimization algorithms usually require a substantially higher number of iterations to converge than ICD. Therefore, the problem is based on finding an optimization algorithm which both converges relatively quickly and is accessible to relatively efficient parallelization.Most iterative algorithms for finding a minimator of the penalty likelihood function Φ use the gradient ≅Φ of the likelihood function Φ. Many such algorithms can be expressed in the following form: where α n >0 equals relaxation factor (or step size), where the N×N matrix D(μ) is a voxel-dependent scaling matrix. The scaling matrix D(μ) may also be referred to herein as a "scaling function" and / or a "scaling factor.". D may be diagonal or non-diagonal depending on the optimization algorithm. In a conventional steepest descent (gradient) algorithm, D is a unit matrix and the appropriate choice of an guarantees convergence and an is calculated by means of a 1D optimization called row search. The relatively slow convergence of the steepest descent gradient method, however, makes this method unattractive for practical purposes. A conjugate gradient (CG) method is superior to steepest descent gradient methods in finding a better search direction. Conditioning the CG algorithm beforehand leads to an acceleration of the convergence rate as a function of the iteration. The ideal preconditioner is the inverse of the Hessian Web Protocol of the cost function, in which case it becomes a Newton's method and can be solved in a single step in the case of quadratic problems. However, in the case of tomographic problems, the Hessian web protocol is spatially invariant and relatively large, and is difficult to invert. Detailed diagonal and circular approximation methods have been proposed for the Hessian Web protocol of the cost function. Implicit preconditioners may also be derived.The problem can also be solved by the principles of a transfer of the optimization and the generation of a replacement function φ(..; μ n) which can be more easily maximized to calculate μ n+1 by using the current iteration μ n. The substitute function is formulated to have the following properties:The optimization transfer algorithm, for example, the separable parabolic equivalent function (SPS) algorithm, alternating minimization algorithms, and expectation maximization algorithms, use an equivalent function.Block-iterative methods have been proposed (also referred to as "ordered subset" methods). Such algorithms can substantially speed up the convergence of many simultaneous algorithms with suitable block dependent scaling. Algorithms using ordered subsets use range decomposition of the forward projector to speed up the underlying algorithm. Each subset is based on a subset of the data, and acceleration is achieved by approximating the entire gradient using the subset of the data. This results in an acceleration of about the number of subsets. However, algorithms accelerated by ordered subsets are not monotonic and enter an undesired limit cycle. The size of the limit cycle depends on the number of subsets, and larger subsets result in a larger limit cycle. A larger limit cycle is not desired and may result in a considerable deviation from the actual desired solution upon completion of the algorithm.The simultaneous and "block-iterative" methods are usually highly parallelizable, but do not have sufficient flexibility in affecting the updates based on historical data, as proposed in NH-ICD. Such algorithms usually have a relatively fast low frequency convergence and a rather slow high frequency convergence. It has been found that NH-ICD significantly speeds up the convergence of the underlying ICD algorithm by updating those voxels that are further from the convergence more frequently. In the method 100 illustrated in FIG. 1, a new approach to calculating an "adjustable" voxel-dependent scaling factor D(μ) (equation 4) for the updates is discussed, which results in more rapid convergence than the conventional gradient-based optimization algorithms. In known reconstruction methods, the term "rescaling" D has been used to speed up the convergence of "blockitive" algorithms, however, the resscale factor is a fixed factor determined by the geometry, does not change with the iteration index, and does not adapt to the data. In the methods described and / or illustrated herein, the proposed scaling factor may be adjusted to the data, may change with the iteration number, and may be calculated based on factors such as historical updates, acquisition geometry, data values, pixels, and / or the like.In geometries such as spiral CT, the volume elements are not uniformly scanned by the scan trajectory. Usually, a three-dimensional (3D) region of interest (ROI) is reconstructed, which is specified by the user. The number of views is determined by the ROI, and voxels outside the ROI may be scanned by a fewer number of projection views than the number of volume elements inside the ROI. During the convergence process, these voxels may have considerable settling artifacts and degrade the image quality of intermediate approximation results. In the method 200 illustrated in FIG. 5, a new approach for relaxing (or scaling) the updates μ n+1- μ n( see equation 4) as a function of the acquisition geometry is described and illustrated. The updates can be relaxed (or scaled) either by suitably scaling the gradient ≅φ, or by modifying the step size αn. The method 200 may reduce transients and improve image quality. In the literature for "blockitive" algorithms, various configurations for the relaxation factor αnhave been proposed to ensure convergence. However, the scaling factor in that context is an iteration dependent sequence of scalars. The scaling factors described herein are spatially varying values that may depend on factors such as acquisition geometry, historical updates, data values, pixels, and / or the like.Referring to equations 1 and 4, the iterative optimization algorithms for the MBIR cost function require that the gradient of the regularizer R(μ) be calculated at each iteration. The regularization function is usually based on performance, logarithm and trigonometric functions. Consequently, the time required to calculate the gradient of the regularizer may be considerable. In the methods 300 and 400 shown in FIGS. 8 and 9, respectively, new proposals are made to shorten the time for the calculation of the regularizer gradient without compromising the monotone behavior of an underlying monotone algorithm or the convergence rate of optimization of the algorithm.In the method 100 illustrated in FIG. 1, various approaches are discussed that may be used to efficiently generate an "adjustable" voxel-dependent table that is used to scale the diagonal matrix D(μ) in Equation 4, resulting in further acceleration of iterative optimization algorithms for the MBIR cost function. Some suggestions are discussed to improve the adaptability of the table and also to more efficiently create the table. An example of rescaling the diagonal matrix D(μ) using the voxel dependent table is shown using an SPS (OS SPS) algorithm using ordered subsets.In the method 200 illustrated in FIG. 5, techniques for appropriately relaxing (or scaling) the updates (μ n+1- μ n) are described that achieve decreased transients and improved image quality. An example of calculating a suitable factor for relaxation (or scaling) from the acquisition geometry is described for the OS SPS algorithm.FIG. 1 shows a flow chart of the exemplary method 100 for the stepwise reconstruction of an image of an object. Measurement data associated with the image is acquired using an exemplary imaging system. In step 102, the method 100 includes the step of accessing the measurement data. In step 104, a simultaneous algorithm is used to reconstruct the image. The use of the simultaneous algorithm in step 104 for reconstruction of the image includes in step 106 determining a scaling factor D(μ) that is voxel-dependent. In the embodiment of the method 100, the scaling factor is determined in step 106 by generating a coefficient table in step 108. Thibault et al. found that more frequent updating of voxels farther from the solution may speed up convergence. Therefore, it is necessary to identify such volume elements, and the "coefficient table" should be based on such volume elements. In NH-ICD, the coefficient table may be directly integrated to choose the order of updating voxels. For simultaneous algorithms, the diagonal scaling matrix D(μ) may be modified to integrate the coefficient table. Care should be taken to ensure that the monotone behavior of the algorithms is not changed.An obvious choice for such a coefficient table is the difference of at least two iterations or at least two sub-iterations. The at least two iterations or sub-iterations may or may not follow one another. In one embodiment, the coefficient table is the difference of two consecutive iterations:The difference in equation 6 should emphasize the voxels that change most and are thus likely to be farther from the solution. This difference varies with the iteration. Therefore, the coefficient table may be a function of the iteration and may also be adapted to the data. The coefficient table u nj may be transformed based on the histogram of the difference in Equation 6. The transformation may be discrete or continuous. The ideal table is locally continuous and the final table can therefore be low-pass filtered on the basis of an intermediate table. However, the coefficient table may not be available before completion of a complete iteration, so an initial coefficient table may be generated from the initial image. The initial image is usually reconstructed using the filtered backprojection (FBP) algorithm, but the choice is not limited to FBP. The initial table may be generated by applying an edge and / or intensity determination detector to the initial image, and may also depend on the scan geometry. Given that the initial image has an adequate low frequency component, the edge and intensity determination detector may be an appropriate choice for improving high frequency convergence.The coefficient table provides a useful basis for adaptively selecting the scale for the image volume. The decision to create the coefficient table may be augmented using several factors, such as, but not limited to, prior knowledge of the initial image, data about the acquisition geometry, image elements such as features of the current image (e.g., but not limited to, edges, homogeneous regions, and / or the like), and / or the like. For example, FBP algorithms exhibit artifacts due to the break-off of the scanned object outside the scan field of view, and this data can be used to additionally scale the updates in the broken-off area. Another example is due to the fact that the simultaneous algorithms have a relatively slow high frequency convergence.Therefore, the coefficient table may be augmented with a table obtained from a high pass filtered image of μn, unat each iteration index.An example of a coefficient table u 4j: u nj is shown in FIG. 2, wherein the coefficient table is designed discretely by means of the histogram of the table. The largest 5% volume elements are assigned a value of 16, and the next 10%, 20% and 40% volume elements are assigned values of 8, 4 and 2, respectively. The remaining voxels are assigned a value of 1, followed by low pass filtering. Adaptively scaled SPS updates perform larger updates on the voxels that are bright in that figure, whereas ordinary SPS equivalently updates all voxels.Referring again to FIG. 1, a coefficient table is used after the generation of the coefficient table to influence the diagonal scaling matrix D(μ). Accordingly, in step 110, the method 100 includes using the coefficient table to calculate the voxel-dependent scaling factor D(μ). After the voxel-dependent scaling factor D(μ) is determined, the method 100 includes, in step 112, applying the voxel-dependent scaling factor D(μ) to the gradient ≅Φ of the objective function to reconstruct the image.The coefficient table should influence the voxel dependent scale factor D(μ) without compromising the monotonic behavior of the original algorithm. It is relatively unproblematic when the original algorithm, for example the "block-repetitive" algorithm, is not monotone. A practical aspect that should be considered unconditionally is the efficiency of using the coefficient table for the scaling of the diagonal scaling matrix in the execution. Some possible proposals are described below.The suggestion of using the coefficient table to scale the diagonal scaling matrix D(μ) can be used for most simultaneous algorithms, for example, but not limited to, CG algorithms with implicit or apparent preconditioners. An exemplary OS SPS algorithm is described below, illustrating an embodiment of method 100. The OS SPS algorithm is merely an example of rescaling a diagonal matrix D(μ) using a voxel-dependent coefficient table.The constrained least squares weighted cost function is expressed by: where y is equal to the projection data with M elements, w is equal to the statistical weighting matrix, A is equal to the system matrix, μ is equal to the image with N elements, C is equal to the finite differentiating matrix, and ψ is equal to any Markov Network (MRF) predecessor that contains directional coefficients. The separable replacement algorithm is motivated on the basis of the convex decomposition lemma proposed by De Pierce, which permits a comparatively unproblematic inversion of the Hessian Web protocol of the replacement function and a very simple incorporation of the restriction to the non-negative numerical space. The replacement function for the data adaptation term and for the regularization term can be derived by means of the following transformations:By means of equation 8, the substitute function can be established as:The use of and the maximum curve for the replacement function of the regularizer results in a form of the SPS algorithm which can be expressed by:In Equation 10, the diagonal scaling matrix is the inverse of the Hessian Web Protocol of the replacement cost function (Equation 9). The Hessian web protocol is a diagonal matrix and thus can be relatively easily inverted. The choice of this scaling matrix merely makes it dependent on both the statistics and the geometry of the data. Therefore, the scaling matrix is fixed and does not vary with the iteration index.The non-uniform table u j can be integrated into the scaling matrix by modification of π ij and κ rj. By employing updates for SPS with non-uniform updates (NU-SPS) can be expressed by:The inclusion of u j results in an adjustable d .,nj which is dependent on the iteration and can adapt to the data. Equation 11 uses a full forward and back projection to calculate d .,nj using the coefficient table u nj. FIG. 3(a) depicts the square-average (RMS) difference (in HU) from the converged image versus the iteration with patient data, and shows that adaptively scaling the diagonal matrix approximately triples the convergence rate. Figure 4 also confirms that the reconstructed image with the new scaling matrix is closer to the converged image than that without the new scaling. FIG. 4( a) illustrates an initial FBP image. Fig. 4(b) shows the converged image. FIGS. 4(c) and 4(d) illustrate reconstructed images after 20 iterations using (c) usual OS SPS (82 subsets) and (d) adaptively scaled OS SPS (82 subsets). Moreover, FIG. 3( b) shows that the convergence rates of high-frequency components are more than tripled. (The convergence of high frequency terms is measured by calculating the full width at half maximum (FWHM) of a high contrast wire.)A practical aspect that should be considered unconditionally is the efficiency of utilizing the coefficient table for the scaling of the diagonal scaling matrix. The suggestions are specific to the implementation and may vary from algorithm to algorithm. Some possible proposals for reducing the computation paths are: 1. the frequency of calculation of the coefficient table can be reduced and it can be calculated after a few iterations (e.g. after every nth iteration) in each case instead of every individual iteration. 2. the coefficient table may be down sampled to reduce the time for the forward and back projections. The coefficient table is locally steady and down sampling thus cannot significantly affect convergence. 3. the updates in equation 11 require a computation of the gradient ≅ q requiring forward and back projection. The forward and back projection for calculation of πijmay be performed simultaneously with the calculation of the gradient. This can reduce the time required to calculate the diagonal scaling matrix D(μ). For example, the gradient is calculated in OS algorithms by means of a partial forward and back projection. In this case, the new scaling matrix can be computed with the gradient, but the full computational path would require to traverse a portion, if not all, of the subsets.FIG. 5 shows a flow diagram of a further exemplary method 200 for stepwise reconstruction of an image of an object. In some acquisition geometries, e.g., multi-slice helical CT, the image volume is sampled non-uniformly. As shown in FIG. 6( b), the volume elements within the three-dimensional region of interest (3D ROI) are scanned by a greater number of projection views than the voxels outside the ROI. These subsampled voxels can slow down the convergence of algorithms considerably and cause transient artifacts if the updates to the current image (Equation 4) are not properly scaled.Referring again to FIG. 5, measurement data associated with the image is acquired using an example imaging system. In step 202, the method 200 includes the step of accessing the measurement data. In this case, the measurement data are sampled non-uniformly in a 3D acquisition geometry. In step 204, a simultaneous algorithm is used to reconstruct the image. The use of the simultaneous algorithm in step 204 for reconstruction of the image includes in step 206 determining a scaling factor D(μ) that is voxel-dependent. In the exemplary embodiment of method 200, the voxel-dependent scaling factor is generated on the basis of the acquisition geometry. Accordingly, in step 208, the method 200 includes calculating the voxel-dependent scaling factor directly based on the acquisition geometry.After the voxel-dependent scaling factor D(μ) is determined, the method 200 includes, in step 210, applying the voxel-dependent scaling factor D(μ) to the gradient ≅Φ (see equation 4) of the objective function to reconstruct the image. The gradient ≅Φmay be scaled, or the step width αnmay be modified in terms of stability. For example, in some embodiments, applying the voxel-dependent scaling factor to the gradient in step 210 includes scaling the gradient using the voxel-dependent scaling factor. The scaling can be achieved by means of AT WA (see equation 15). The scaling may be associated with the contribution of the subsets to each voxel, which results in appropriate scaling of the updates. In further embodiments, applying the voxel-dependent scaling factor to the gradient in step 210 includes calculating a relaxation factor for each layer of the image data. The scanning of the image volume is non-uniform, and slices within the ROI are scanned through a greater number of views than slices outside the ROI. The step size (or relaxation factor) is usually a scalar, but may be modified to be a vector for each layer. This may allow the slices outside the ROI to relax more than the slices inside the ROI. However, this can increase the size of the search space for α and increase the computation time.The suggestion of relaxing (or scaling) the updates for subsampled volume elements is applicable to all simultaneous algorithms. An exemplary OS SPS algorithm is described below, illustrating an embodiment of method 200. The OS SPS algorithm is merely an example of the structure of the scaling factor and of the use of the scaling factor for scaling the gradient ≅Φ.The update for the SPS algorithm can be expressed using equation 10: where dj is equal to the scaling function (see equation 4) determined from the Hessian Web protocol of the replacement cost function. dj contains the data adjustment term and the regularizer term expressed in equation 10. Now, in OS-SPS, the gradients are calculated using the subset of views and a complete iteration corresponds to the S number of all subsets being traversed. Each iteration uses multiple sub-iterations and each sub-iteration requires a calculation of the gradient of the regularizer and the data adjustment term using the subset of the data. The OS SPS updates can be expressed by: where s is equal to the subset index and n is equal to the iteration index. If the gradient computed using the subset of views is properly scaled, it may match the gradient computed using all the views. This accelerates the convergence by a factor S compared to the original SPS algorithm. However, this ideal scaling factor cannot be calculated in practice. The conventional practice of OS SPS is to use a constant scaling factor γ n,sj= S corresponding to the number of subsets. However, this choice leads to poor results if the ideal scaling factor γ n,sj is clearly non-uniform. This non-uniformity in scaling may result in the algorithm having significant transition artifacts that are not desired. One of the reasons for non-uniformity may be the geometry of the acquisition. The geometry of the acquisition can result in non-uniform sampling of the image volume. Voxel-dependent scaling of the gradient can be developed based on the geometry of the acquisition. One such non-uniform table choice is: with I{B} equal to 1 if B is true and equal to 0 if B is false. In Equation 15, B is true when a given volume element indexed j holds contributions from views associated with the subset indexed s. The non-uniform table and denominator may be calculated simultaneously.FIG. 6( b) shows that the correctly scaled OS algorithm provides improved image quality compared to the usual OS algorithm that does not converge very quickly outside the ROI. The instability that appeared outside the ROI and affects the image quality inside the ROI (as seen from the standard deviations of the noise in FIG. 6( b)) was reduced by the proposed relaxation (or scaling). Figure 6(a) illustrates a diagram of a spiral geometry. A region drawn in dashed lines denotes the detector, which acquires measurement data and records contributions of both volume elements within the ROI and volume elements outside the ROI. Figure 6(b) illustrates the effect of gradient scaling in the OS SPS spiral geometry algorithm. Each image is reconstructed after 30 iterations of the OS algorithm are run through with 328 subsets using the usual and proposed scaling approaches. For comparison, the standard deviation σ of a uniform area (in the box) is calculated. (Some iterations of a convergent algorithm are shown as a reference value.).FIG. 7 also demonstrates that the correctly scaled OS within the ROI is robust, whereas the usual OS algorithm is corrupted due to its instability outside the ROI within the ROI. FIG. 7 illustrates a mean and standard deviation in a uniform region of an end slice image of the ROI versus iteration, showing the instability of the common OS approach compared to the proposed OS approach (328 subsets).FIG. 8 shows a flow diagram of a further exemplary method 300 for gradually reconstructing an image of an object. FIG. 8 illustrates an embodiment of the calculation of the regularizer gradient. Iterative optimization algorithms for the MBIR cost function require a calculation of the gradient of the regularization function to update the current image approximation (see equation 4). The time for calculating the gradient of the regularizer may be considerable depending on the regularizer potential function Φ. Quite generally, the potential function includes powers, logarithms and trigonometric functions, which require a great amount of computing effort. Examples of such potential functions include, but are not limited to:To reduce the cost required to calculate the gradient of the regularization function, an approximate regularizer gradient value ≅ Φ R( μ) may be used that is based on gradient samples in a look-up table as described in method 300. In a variation, a replacement function as described in method 400 of FIG. 9 may be established for the regularization function, wherein the gradient of the replacement function may be calculated more efficiently than the gradient of the current regularization function.With reference to the method 300, measurement data associated with the image is acquired via an example imaging system. In step 302, the method 300 includes the step of accessing the measurement data. In step 304, a simultaneous algorithm is used to reconstruct the image. Using the simultaneous algorithm in step 304 to reconstruct the image includes storing multiple gradient samples in a look-up table in step 306. Optionally, the storing in step 306 includes sampling the gradient values at a fixed sampling frequency. One approach to storing the gradient samples in the look-up table will typically be based on sampling, i.e., tabulating, the values of ≅ Φ R( μ) in step 306. If the points μ are k= kΔ, where Δ is a positive integer indicating the sampling frequency, the following indexing function may be used for the table (where [.] denotes a rounding function):In step 308, the method 300 includes accessing the gradient samples in the look-up table. In step 310, an approximated actual gradient value is then calculated from the look-up table using the gradient samples. In step 310, the approximated actual gradient value may be calculated using any interpolation method, such as, but not limited to, a linear interpolation method, a bi-cubic interpolation method, and / or the like. In one embodiment, calculating the actual approximated gradient value in step 310 includes calculating an actual approximated gradient value for each iteration of the reconstruction.FIG. 9 shows a flow diagram of a further exemplary method 400 for gradually reconstructing an image of an object. FIG. 9 illustrates an embodiment of the calculation of a regularizer gradient, wherein, as briefly described above, a replacement function is created for the regularizer function. Measurement data associated with the image is acquired using an exemplary imaging system. The method 400 includes, in step 402, accessing the measurement data. In step 404, a simultaneous algorithm is used to reconstruct the image.Using the simultaneous algorithm to reconstruct the image in step 404 includes, in step 406, calculating an actual gradient at an iteration of the reconstruction of the image. The replacement function for the regularization function is then constructed in step 408 by adding a compensation term to the actual gradient. The compensation term may include relatively simple mathematical operations. For example, in some embodiments, the compensation term includes only at least one of addition, subtraction, division, and / or multiplication.The actual gradient and the compensation term are then used in step 410 as the gradient for a subsequent iteration of the reconstruction of the image. For each subsequent iteration of the reconstruction of the image, the actual gradient calculated in step 406 is added to a compensation term for use other than the gradient for the particular subsequent iteration. In one embodiment, the gradient used in step 410 for a particular iteration is used as the gradient for all sub-iterations of the particular iteration. That is, in some embodiments, the gradient used in step 410 for a particular iteration is not recomputed for the sub-iterations of the particular iteration, but rather is reused for each sub-iteration of the particular iteration.The suggestion of a substitute function is applicable to any simultaneous reconstruction algorithm. An exemplary OS SPS algorithm is described below, illustrating an embodiment of method 400. The OS SPS algorithm uses a replacement function for the data adjustment term and the regularizer term. Above the existing replacement function, another replacement function is introduced for the regularizer to minimize the computation time for the gradient calculation; thus, the OS SPS algorithm may be referred to as an "ordered subsets with dual replacement functions" algorithm.Assume a general PL objective function of the following form: where Φ L. equals the data adjustment term and Φ R equals the regularizer. It is assumed that the data adjustment term.phi.L has a quadratic substitute function of the following form: with a suitable diagonal matrix DL. A typical choice is D L= diag(d Q,nj). It is further assumed that the regularizer Φ R( μ) contains a quadratic substitute function of the following shape: with a suitable diagonal matrix D L. A typical choice is D R( μ)=diag(d R,nj). The following "double substitute" function is then defined:By design, this quadratic replacement function has the following properties: For a method using conventional regularized ordered subsets, the minimization step for each subset is given by: where m=1,..., M denotes the index of the subset. Each iteration indexed by an n ordered subset uses M sub-iterations. This iteration is undesirably slow because it calculates the regularization gradient ≅φ R for each subset. To reduce this effort, it is proposed to take advantage of the generality of the double-substitute function equation 23 by using the following novel update: with Dφ(μ( n,m-1), μ (n)) = [ D L( μ (n,m-1)) + D R( μ (n))]-1. This new iteration equation 27 uses the same regularizer gradient for all subsets. Compared to equation 24, the updates in equation 27 are very similar except for an additional term that compensates for the absence of the regularizer gradient update. However, this additional term has simpler operations (e.g., additions, subtractions, divisions, and / or multiplications) and can therefore be calculated more quickly. In the above description, the regularizer gradient was not updated until all subsets were updated. This involves calculating the gradient once at the beginning of each iteration and re-using this calculation for all subsets, rather than calculating the gradient M times for each iteration in the original OS method described in Equation 24. It is understood that the regularizer gradient can be updated as many times as needed, and the update frequency is denoted by Uf. Updating the regularizer gradient less frequently reduces computational costs at the expense of the convergence rate in early iterations.The above-described OS-SPS algorithm has been studied in view of the problem of three-dimensional cone beam CT image reconstruction with limited viewing angles.To evaluate the convergence rate of method 400, the quadratic mean difference between the image approximation at the nth iteration μ(n) was calculated with both the fully converged solution μ(∞) and the true image μ(true).For this type of "relaxed" OS, convergence is ensured because the final phase uses only 1 subset for which equation 24 is convergent. FIG. 10 illustrates the true phantom, FBP reconstruction, and fully converged (μ(∞)) images from cone beam CT data with 164 projection views.The projections were divided into 41 subsets, which corresponds to 4 views per subset. This is a rather aggressive selection compared to conventional selections which attempt to significantly speed up convergence.The regularizer gradient was updated at different frequencies to visualize the effect on convergence and compute time. FIG. 11 illustrates convergence rates at each iteration for different frequencies of regularizer update. In FIG. 11, OS-41-DS-n means an OS having 41 subsets, and Uf=n and "n=all" means that only one update is performed after all subset updates are performed. FIG. 11( a) illustrates μ(n) looking at μ(∞), and FIG. 11( b) illustrates μ(n) looking at μ(true). FIG. 11 shows that the exemplary OS SPS algorithm associated with method 400 yields differences in root mean square (RMS) similar to the conventional OS, even though the regularizer is updated relatively rarely per iteration.On the other hand, FIG. 12 shows that the computational effort required to achieve the same level of RMS differences has been reduced by the example OS SPS algorithm associated with method 400. FIG. 12 illustrates the convergence rate versus computation time for different frequencies of regularizer update. FIG. 12( a) illustrates μ(n) looking at μ(∞), and FIG. 12( b) illustrates μ(n) looking at μ(true). With Uf=13, which in this case yielded the best result, the exemplary OS SPS algorithm in connection with method 400 converges about 3 times faster than the conventional OS.FIG. 13 compares the convergence rate of OSDS with different update frequencies. In the left column, FIG. 13 illustrates images located at the same time (4000 seconds after initialization). In the right column, FIG. 13 shows absolute difference images with respect to μ(∞). Observation of the reconstructed image at the same time is clearly seen from FIG. 13 that the exemplary OS SPS algorithm associated with method 400 converges more quickly. Tradeoffs between convergence rates and computational effort are to be made, and in the present case, calculating the regularizer gradient for every 13 subset updates yielded the most efficient results on the processing hardware used for the trial. For other problems, the optimal frequency of updates may vary. However, it should be appreciated that the example OS SPS algorithm associated with method 400 converges faster than the conventional OS regardless of update frequency. As the problem and the number of subsets grow, the computational effort required to calculate the gradient of the regularizer may be significantly more significant. Significant advantages of the method 400 are therefore to be expected.The methods and algorithms described herein are used to reconstruct an image of an object in a stepwise manner. The methods and algorithms can be executed as a set of instructions stored on a computer and performed using, for example, a module 530, shown in FIG. 13, software, hardware, a combination thereof, and / or a computer readable tangible, non-transitory medium. In one embodiment, a tangible, non-transitory computer readable medium excludes signals.As used herein, the concept of using coefficient tables to improve convergence and image quality (e.g., as discussed in the paragraph) is not limited to model-based cost functions, but rather is applicable to all simultaneous algorithms that solve model-based as well as non-model-based cost functions.FIG. 14 shows an illustrative view of an exemplary imaging system 500 formed in accordance with various embodiments. FIG. 15 is a block diagram of a portion of the multi-mode imaging system 500 shown in FIG. 14. The imaging system may be embodied as, but is not limited to: 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, a C-arm interventional tomography imaging system, a dedicated CT system, e.g., for limb or chest scanning, and combinations thereof. In the exemplary embodiment, method 100 is described in connection with a CT imaging system. Moreover, the cost function in the embodiment is used to explain the different embodiments described herein.Although various embodiments are described in the context of an exemplary dual modality imaging system that includes a computed tomography (CT) imaging system, as well as a positron emission tomography (PET) imaging system or a single photon emission computed tomography (SPECT) system, it should be understood that other imaging systems capable of performing the functions described herein are also contemplated for use.The multi-mode imaging system 500 is illustrated and includes a CT imaging system 502 and a PET imaging system 504. The imaging system 500 allows multiple scans in different imaging modalities to enhance diagnostic capabilities over single modality systems. In an embodiment, the example multi-modal imaging system 500 is a CT / PET imaging system 500. Optionally, other imaging modalities are used in conjunction with imaging system 500 than CT and PET. For example, the imaging system 500 may include, but is not limited to, a stand-alone CT imaging system, a stand-alone 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, a C-arm interventional tomography system, a dedicated CT system, for example for scanning extremities or the breast, and combinations thereof.The CT imaging system 502 includes a gantry 510 having an x-ray source 512 that projects a beam of x-rays toward a detector array 514 disposed on the opposite side of the gantry 510. Detector array 514 includes a number of detector elements 516 arranged in rows and channels that together sense the projected x-rays traversing an object, e.g., patient 506. The imaging system 500 also includes a computer 520 that captures the projection data from the detector array 514 and processes the projection data to reconstruct an image of the patient 506.In operation, control commands and parameters input by the operator are used by the computer 520 to output control signals and data to repeatedly position an electrically driven table 522. More specifically, the electrically powered table 522 is used to move the patient 506 into and out of the gantry 510. Specifically, the table 522 moves at least a portion of the patient 506 through a gantry tunnel 524 that extends through the gantry 510.The imaging system 500 also includes a module 530 configured to perform a variety of methods and algorithms described herein. The module 530 may be implemented as a hardware component installed on the computer 520. Optionally, the module 530 may be performed as a set of instructions installed on the computer 520. The instruction sets may be based on, among other things, stand-alone programs, as subroutines incorporated in an operating system installed on the computer 520, and / or include functions in a software package installed on the computer 520. It is understood that the various embodiments are not limited to the arrangements and functionalities shown in the figures.As discussed above, detector 514 includes a number of detector elements 516. Each detector element 516 generates an electrical signal or output signal indicative of the strength of an incident x-ray beam, thereby allowing the attenuation of the beam experienced by it to be determined as it passes through the patient 506. During a scan to acquire the x-ray projection data, the gantry 510 and the components mounted thereon rotate about an axis of rotation 540. FIG. 15 shows only a single row of detector elements 516 (i.e., a detector row). However, the multi-slice detector matrix 514 has a plurality of parallel detector rows of detector elements 516, so that projection data corresponding to a plurality of slices can be acquired simultaneously during a scan.Rotation of the gantry 510 and operation of the x-ray source 512 are controlled by a controller 542. The controller 542 includes an x-ray controller 544 that outputs power and timing signals to the x-ray source 512 and a gantry drive controller 546 that controls / controls the rotational speed and position of the gantry 510. A data acquisition system (DAS) 548 in the controller 542 samples analog data from detector elements 516 and converts the data to digital signals for subsequent processing. For example, subsequent processing may include utilizing module 530 to perform the various methods described herein. An image reconstructor 550 receives the scanned and digitized x-ray data from the DAS 548 and performs high speed image reconstruction. The reconstructed images are input to the computer 520, which stores the image in a storage device 552.Optionally, the computer 520 may receive the scanned and digitized x-ray data from the DAS 548 and perform a variety of methods described herein using the module 530. Via a console 560 having a keyboard, the computer 520 also receives control commands and scan parameters from a user. An associated display unit 562 allows the operator to observe the reconstructed image and other data output by the computer.The operator input control commands and parameters are used by the computer 520 to output control signals and data to the DAS 548, the x-ray controller 544, and the gantry drive controller 546. In addition, the computer 520 controls a table drive controller 564 that controls / controls the electrically powered table 522 to position the patient 506 within the gantry 510. Specifically, as shown in FIG. 14, the table 522 moves at least a portion of the patient 506 through the gantry opening 524.Referring again to FIG. 15, in one embodiment, the computer 520 includes a device 570, such as a floppy disk drive, a CD-ROM drive, a DVD drive, a Magneto Optical Disc (MOD) drive, or any other digital device, such as a network connection device such as an Ethernet device, for reading commands and / or data from a computer readable medium 572, such as a floppy disk, a CD-ROM, a DVD, or other digital source, such as a network or the Internet, as well as from digital devices to be developed in the future. In yet another embodiment, 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 only 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.In the embodiment, the x-ray source 512 and detector array 514 are rotated together with the gantry 510 in the imaging plane and about the patient 506 to be imaged so that the angle at which an x-ray 574 intersects the patient 506 constantly changes. A set of x-ray attenuation measurements, i.e., projection data, from the detector array 514 at a gantry angle is referred to as a "view.". A "scan" of the patient 506 includes a set of views generated at different gantry angles or view angles during one revolution of the x-ray source 512 and the detector 514. In a CT scan, the projection data is processed to reconstruct an image corresponding to a two-dimensional slice image taken by the patient 506.Embodiments of a multi-mode imaging system are described in detail above. The illustrated components of the multi-modal imaging system are not limited to the specific embodiments described herein, but rather, components of each multi-modal imaging system may be used independently and separately from other components described herein. For example, the components of the multi-mode imaging system described above may also be used in conjunction with other imaging systems.It should be noted that the various embodiments may be implemented as hardware, software, or a combination thereof. The various embodiments and / or components, for example the modules or components and control devices contained therein, can also be realized as part of one or more computers or processors. The computer or processor can contain a computer, 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 data transfer bus. The computer or processor may further include a memory. The memory may include a RAM memory (RAM) and / or a read only memory (ROM). The computer or processor may further include a storage device, which may be a hard disk drive or a removable storage drive, such as a solid state drive (SSD) memory card, an optical drive, and / or the like. The storage device may also be any other similar means for loading computer programs or other instructions into the computer or processor.As used herein, the term "computer" may include any processor- or microprocessor-based system, such as systems including microcontrollers, reduced instruction set computers (RISC), application specific integrated circuits (ASICs), logic circuits GPUs, FPGAs, and other circuits or processors capable of performing the functions described herein. The examples mentioned above are merely exemplary and are therefore not intended to limit the definition and / or meaning of the term "computer" in any way. The computer or processor executes a set of instructions stored in one or more storage elements to process input data. The storage elements may also store data or other information as desired or as needed. The storage element may be in the form of a data source or a physical memory element in a processing machine.The set of instructions may include a variety of control commands that cause 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 implemented in a variety of ways, such as system software or application software, which may be a non-transitory computer readable medium. Moreover, the software may be present as a collection of independent programs, as a program module within a larger program, or as a portion of a program module. The software can also be based on modular programming in the form of object-oriented programming. Processing of input data by the processing engine may be responsive to user commands, or responsive to results of previous processing, or responsive to a request made by another processing engine.As used herein, elements or steps mentioned in the singular, preceded by the indefinite article, should be understood in the sense that the plural elements or steps are not excluded, unless such exclusion is expressly mentioned. Further, reference to "one embodiment" of the present invention is not to be interpreted as excluding the existence of additional embodiments that also embody the recited features. Moreover, unless expressly stated otherwise, embodiments that "comprise" or "contain" one or more elements having a particular characteristic may contain additional elements that do not have the characteristic in question.Moreover, as used herein, the term "reconstructing an image" is not intended to exclude embodiments of the present invention in which data representing an image is generated while a viewable image is not generated. Thus, as used herein, the term "image" broadly refers to both viewable images and data representing a viewable image. However, many embodiments generate at least one viewable image or are configured to do so.As used herein, the terms "software" and "firmware" are interchangeable and include any computer program stored in memory, such as RAM, ROM, EPROM, EEPROM, or non-transitory RAM (NVRAM) memory, for execution by a computer. The types of memories mentioned above are merely exemplary and are therefore not limiting with respect to the types of memories suitable for storing a computer program.It is to be understood that the above description is illustrative and is not intended to be limiting. For example, the above-described embodiments (and / or aspects thereof) may be used in combination with each other. Moreover, many modifications may be made to adapt a particular situation or material to the teachings of the invention without departing from the scope of the invention. While the dimensions and types of materials described herein are intended to define the parameters of the invention, they are by no means limiting and are merely illustrative. Many further embodiments will be apparent to those skilled in the art after reading the above-mentioned description. The scope of the invention should therefore be determined from the appended claims, along with the full scope of equivalent forms to which such claims are entitled. In the appended claims, the terms "include / include / have" and "wherein" are used as plain-text equivalents of the corresponding terms "comprise" and "wherein.". Moreover, the words "first", "second", "third" and the like in the following claims are merely for designation and are not intended to numerically define the objects.The present description uses examples to describe the various embodiments of the invention that include the best mode, and also to enable any person skilled in the art to put the various embodiments of the invention into practice, for example, to make and use any devices and systems and to perform any methods associated therewith. The patentable scope of the various embodiments of the invention is defined by the claims and may include other examples that occur to those skilled in the art. Such other examples are intended to fall within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insubstantial differences from the literal language of the claims.A method for iteratively reconstructing an image of an object is provided. The method includes the steps of: accessing measurement data associated with the image and using a simultaneous algorithm to reconstruct the image. Using the simultaneous algorithm to reconstruct the image includes the steps of: determining a scaling factor that is voxel dependent; and applying the voxel dependent scaling factor to a gradient of a target function to reconstruct the image.

Claims

A method for gradually reconstructing an image of an object, the method comprising the steps of: accessing measurement data associated with the image; and using a simultaneous algorithm to reconstruct the image, wherein using the simultaneous algorithm to reconstruct the image comprises the steps of: determining a scaling factor that is voxel dependent; and applying the voxel dependent scaling factor to a gradient of a target function to reconstruct the image, wherein determining a scaling factor that is voxel dependent comprises the step of calculating the scaling factor immediately based on an acquisition geometry.The method of claim 1, wherein determining a scaling factor that is voxel dependent includes the step of calculating the scaling factor directly based on an acquisition geometry, wherein the measurement data is non-uniformly sampled in a three-dimensional (3D) acquisition geometry.The method of claim 1, wherein determining a scaling factor that is voxel dependent includes the steps of generating a coefficient table and using the coefficient table to calculate the voxel dependent scaling factor.The method of claim 3, wherein generating a coefficient table includes the step of transforming an update obtained from a difference of either at least two iterations or at least two sub-iterations.The method of claim 3, wherein generating a coefficient table includes the step of generating a coefficient table that is locally steady and that is low pass filterable.The method of claim 3, wherein generating a coefficient table includes the step of obtaining the coefficient table from an initial image.The method of claim 3, wherein generating a coefficient table includes the step of expanding the coefficient table with at least past data about a starting image, about a feature of a current image, and / or about an acquisition geometry.The method of claim 1, wherein applying the voxel-dependent scaling factor to a gradient of a target function to reconstruct the image includes the step of scaling the gradient using the voxel-dependent scaling factor.The method of claim 8, wherein the scaling is at least one of obtained by ATWA and / or associated with the contribution of the subsets to each voxel.The method of claim 1, wherein applying the voxel dependent scaling factor to a gradient of a target function to reconstruct the image includes the step of calculating a relaxation factor for each slice of the measurement data, wherein the measurement data is non-uniformly sampled in a three-dimensional (3D) acquisition geometry, wherein slices within a region of interest (ROI) are sampled by a different number of views than slices outside the ROI.The method of claim 1, wherein applying the voxel-dependent scaling factor to a gradient of a target function to reconstruct the image includes updating the voxel-dependent scaling factor at a frequency that is decreased at each nth iteration rather than at each individual iteration.The method of claim 3, wherein the coefficient table is down sampled.An imaging system comprising: a detector array; and a computer coupled to the detector array, the computer configured to perform the steps of: accessing measurement data associated with an image of an object; and using a simultaneous algorithm to reconstruct the image, the computer configured to use the simultaneous algorithm to reconstruct the image by determining a scaling factor that is voxel dependent and applying the voxel dependent scaling factor to a gradient of a target function to reconstruct the image, wherein the computer is configured to determine the scaling factor that is voxel dependent by calculating the scaling factor directly based on an acquisition geometry.The imaging system of claim 13, wherein the computer is configured to calculate the scaling factor that is voxel dependent by generating a coefficient table and using the coefficient table to determine the voxel dependent scaling factor.The imaging system of claim 14, wherein the computer is configured to generate the coefficient table by transforming an update obtained from a difference of either at least two iterations or at least two sub-iterations, wherein the transformation is either discrete or continuous.The imaging system of claim 13, wherein the computer is configured to apply the voxel-dependent scaling factor to the gradient of the objective function by scaling the gradient using the voxel-dependent scaling factor.The imaging system of claim 13, wherein the computer is configured to apply the voxel-dependent scaling factor to the gradient of the target function by calculating a relaxation factor for each slice of the measurement data, wherein the measurement data is non-uniformly sampled in a three-dimensional (3D) acquisition geometry, wherein slices within a region of interest (ROI) are sampled by a different number of views than slices outside the ROI.

Images