Motion estimation and compensation in cone beam computed tomography (CBCT)

Abstract

In various embodiments, the present invention provides an improved system and method for motion estimation and motion compensation in cone beam computer tomography (CBCT). The present invention utilizes a method for alternating minimization between volume and motion sub-iterations and results in an improved system and method for motion estimation and compensation in CBCT.

Claims

1. A method for motion estimation and compensation using projection scan data, the method comprising: acquiring tomographic projection data of a moving object of interest; sorting the projection data according to a phase to generate phase-sorted projection data; loading the phase-sorted projection data; performing an iterative optimization procedure comprising; changing one or more regularization parameters during the iterative optimization procedure, wherein the one or more regularization parameters affect a strength of one or more regularizers of a cost functional and wherein the one or more regularization parameters includes an optical flow constraint (OFC) strength and a volume regularization strength; performing a volume sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional; performing a motion sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional; and repeating the iterative optimization procedure until a stopping criterion has been met.

2. The method of claim 1, wherein the strength of the OFC increases during the optimization procedure and the strength of the volume regularization decreases during the iterative optimization procedure.

3. The method of claim 1, wherein the one or more of the regularization parameters are changed non-monotonically during the iterative optimization procedure.

4. The method of claim 1, wherein the OFC strength is changed non-monotonically during the iterative optimization procedure.

5. The method of claim 1, wherein the motion sub-iteration is based on phase motion vector fields (MVF) to describe motion.

6. The method of claim 5, wherein phase MVF act between neighboring phase images and each MVF is specified on a regular grid to ensure the uniform application of the optical flow constraint.

7. The method of claim 1, wherein the motion sub-iteration is based on a global in time motion function to describe motion.

8. The method claim 1, wherein the projection data is selected from cone beam computed tomography (CBCT) projection data, dental CT projection data, clinical Computed Tomography (“CT”) data, industrial CT projection data, flat panel CT (e.g., C-arms) projection data, x-ray diffraction projection data, x-ray based imaging projection data and projection data not based on x-ray imaging.

9. A method for motion estimation and compensation using projection scan data, the method comprising: acquiring tomographic projection data of a moving object of interest; sorting the projection data according to a phase to generate phase-sorted projection data; loading the phase-sorted projection data; performing an iterative optimization procedure comprising; performing a volume sub-iteration of the phase-sorted projection data using one or more regularization parameters and one or more regularizers of a cost functional, wherein the volume sub-iteration is based on phase motion vector fields (MVF) to describe motion, wherein the one or more regularization parameters affect a strength of the one or more regularizers of the cost functional and wherein the one or more regularization parameters includes an optical flow constraint (OFC) strength and a volume regularization strength; performing a motion sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional, wherein the motion sub-iteration is based on a global in time deformation function to describe motion; and repeating the iterative optimization procedure until a stopping criterion has been met.

10. The method of claim 9, wherein one or more of the regularization parameters are changed during the iterative optimization procedure.

11. The method of claim 10, wherein the strength of the OFC increases during the iterative optimization procedure and the strength of the volume regularization decreases during the iterative optimization procedure.

12. The method of claim 10, wherein one or more of the regularization parameters changes non-monotonically during the iterative optimization procedure.

13. The method of claim 12, wherein the strength of the OFC changes non-monotonically during the iterative optimization procedure.

14. The method of claim 9, wherein the optical flow constraint (OFC) is used to regularize volume sub-iterations and wherein phase MVF act between neighboring phase images and each MVF is specified on a regular grid to ensure a uniform application of the OFC.

15. A method for motion estimation and compensation using projection scan data, the method comprising: acquiring tomographic projection data of a moving object of interest; sorting the projection data according to a phase to generate phase-sorted projection data; loading the phase-sorted projection data; performing an iterative optimization procedure comprising; performing a volume sub-iteration of the phase-sorted projection data using one or more regularization parameters and one or more regularizers of a cost functional, wherein the one or more regularization parameters affect a strength of the one or more regularizers of the cost functional and wherein the one or more regularization parameters includes an optical flow constraint (OFC) strength and a volume regularization strength; performing a motion sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional, wherein the motion sub-iteration is based on a global in time deformation function to describe motion, wherein the deformation function encodes the direction and phase of the motion of voxels in the object of interest and wherein the phase of the motion is used for motion regularization; and repeating the iterative optimization procedure until a stopping criterion has been met.

16. The method of claim 15, wherein the strength of the OFC changes non-monotonically during the iterative optimization procedure.

17. The method of claim 15, wherein the strength of the OFC increases during the iterative optimization procedure and the strength of the volume regularization decreases during the iterative optimization procedure.

18. The method of claim 15, wherein the volume sub-iterations use motion vector fields (MVF) between phase images to describe motion.

19. One or more non-transitory computer-readable media having computer-executable instructions for performing a method of running a software program on a computing device for motion estimation and compensation of a CT image, the computing device operating under an operating system, the method including issuing instructions from the software program comprising: acquiring tomographic projection data of a moving object of interest; sorting the projection data according to a phase to generate phase-sorted projection data; loading the phase-sorted projection data; performing an iterative optimization procedure comprising; changing one or more regularization parameters during the iterative optimization procedure, wherein the one or more regularization parameters affect a strength of one or more regularizers of a cost functional and wherein the one or more regularization parameters includes an optical flow constraint (OFC) strength and a volume regularization strength; performing a volume sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional; performing a motion sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional; and repeating the iterative optimization procedure until a stopping criterion has been met.

20. The media of claim 19, wherein the strength of the OFC increases during the optimization procedure and the strength of the volume regularization decreases during the iterative optimization procedure.

21. The media of claim 19, wherein the one or more of the regularization parameters are changed non-monotonically during the iterative optimization procedure.

22. The media of claim 19, wherein the motion sub-iteration is based on motion vector fields to describe motion.

23. The media of claim 19, wherein the motion sub-iteration is based on a global in time motion function to describe motion.

24. A system for motion estimation and compensation using projection data of an imaged object of interest, the system comprising: a memory for storing projection data of an object of interest; a data processor for estimating and compensating for motion in the projection data, wherein the data processor is adapted for performing the following operations: sorting the projection data according to a phase to generate phase-sorted projection data; loading the phase-sorted projection data; performing an iterative optimization procedure comprising; changing one or more regularization parameters during the iterative optimization procedure, wherein the one or more regularization parameters affect a strength of one or more regularizers of a cost functional and wherein the one or more regularization parameters includes an optical flow constraint (OFC) strength and a volume regularization strength; performing a volume sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional; performing a motion sub-iteration of the phase-sorted projection data using the one or more regularization parameters and the one or more regularizers of the cost functional; and repeating the iterative optimization procedure until a stopping criterion has been met.

25. The system of claim 24 wherein the motion sub-iteration is based on motion vector fields to describe motion or on a global in time motion function to describe motion.

Description

BRIEF DESCRIPTION OF THE FIGURES

(1) For a fuller understanding of the invention, reference should be made to the following detailed description, taken in connection with the accompanying drawings, in which:

(2) FIG. 1 is an illustration of an exemplary helical computer tomography (CT) scanner, in accordance with an embodiment of the present invention.

(3) FIG. 2 is a flow diagram illustrating the global iteration process for motion estimation, in accordance with an embodiment of the present invention.

(4) FIG. 3 is a flow diagram illustrating the volume sub-iteration process for motion estimation, in accordance with an embodiment of the present invention.

(5) FIG. 4 is a flow diagram illustrating the motion sub-iteration process for motion estimation, in accordance with an embodiment of the present invention.

(6) FIG. 5A is a graphical illustration of monotonically incrementing the parameters along the global iterations, in accordance with an embodiment of the present invention.

(7) FIG. 5B is a graphical illustration of non-monotonically incrementing the parameters along the global iterations, in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

(8) Motion of an object of interest commonly occurs during an imaging scan, such as the natural breathing motion of a patient during a computed tomography (CT) imaging of the patient's abdomen. In various embodiments, the present invention provides a system and method for estimating and compensating for this motion during the reconstruction of the image from the projection data acquired by the CT scanner.

(9) FIG. 1 illustrates an exemplary embodiment of a helical CT scanner system 100 according to the present invention. This embodiment is intended to be exemplary and the present invention will be described as relating to medical imaging. However, this is not intended to be limited and various embodiments of the present invention may also be used for other purposes, such as baggage scanning for security and material analysis. Additionally, while a circular CT scanner system is illustrated numerous other scanning devices for collecting projection data are within the scope of the present invention.

(10) As shown in FIG. 1, a circular CT scanner system 100 includes a gantry 102 which is rotatable around a rotational axis 105. The gantry 102 carries a radiation source 110 that forms a cone shaped radiation beam 115. The radiation beam 115 emitted from the radiation source 110 is focused on an object of interest 120 positioned in the center of the gantry 102. A radiation detector 125 is positioned on the gantry 102 opposite from the radiation source 110. The radiation detector 125 comprises a plurality of detector elements for measuring the intensity of the cone shaped radiation beam 115 after it passes through the object of interest 120.

(11) During a scan of the object of interest 120, the radiation source 110 and the radiation detector 125 are rotated with the gantry 102 in a direction indicated 145. The object of interest 120 is additionally positioned on a stationary table 130.

(12) Following the circular scanning of the object of interest 120, the radiation detector 125 provides the collected helical CT scan data to a data processor 135. The data processor 135 is adapted for reconstructing an image from the measurements provided by the radiation detector 125. The image generated by the data processor 135 may then be provided to a display 140 for subsequent viewing.

(13) The data processor 135 is additionally adapted to perform motion estimation and motion compensation to correct for motion in the scan data provided by the radiation detector 125. Motion estimation and compensation may be performed by the data processor 135 as detailed below.

(14) Throughout the detailed description of the invention, the following notations will be used:

(15) i.sub.q qth axis volume grid index

(16) j.sub.q axis velocity grid index

(17) k detector grid index

(18) m source position (projection data) index

(19) n.sub.g global iteration step index

(20) n.sub.s sub-iteration step index

(21) p phase bin index

(22) q spatial coordinate axis index

(23) K total number of detector pixels

(24) M total number of source positions (projection data)

(25) N total number of volume voxels

(26) Ñ total number of velocity voxels

(27) N.sub.g total number of global iteration steps

(28) N.sub.p total phase of bins

(29) N.sub.q total number of volume voxels along qth axis

(30) Ñ.sub.q total number of velocity voxels along qth axis

(31) N.sub.s total number of volume sub-iteration steps

(32) Ñ.sub.s total number of velocity sub-iteration steps

(33) In the method of the present invention, the reconstruction volume grid is given by:
{right arrow over (x)}.sub.1=(h.sub.1i.sub.1,h.sub.2i.sub.2,h.sub.3i.sub.3),i=(i.sub.1,i.sub.2,i.sub.3), for 0≤i.sub.q<N.sub.q,q=1,2,3.  (1.1)

(34) If the values of a function ƒ are given on the grid (1.1): ƒ.sub.i=ƒ({right arrow over (x)}.sub.i), then the interpolated function becomes:

(35) $\begin{array}{cc}f\left(\stackrel{.fwdarw.}{x}\right)\approx \underset{i}{.Math.}{f}_i{\phi }_h\left(\stackrel{.fwdarw.}{x}-{\stackrel{.fwdarw.}{x}}_i\right).& \left(1.2\right)\end{array}$

(36) The deformation (or, ‘motion’) of the object is described by the function {right arrow over (ψ)}(s, {right arrow over (x)}). Here s is time of phase point of a breathing cycle of a patient being imaged.
{right arrow over (y)}={right arrow over (ψ)}(s,{right arrow over (x)}),  (1.3)

(37) The Equation (1.3) means that the physical point, which at current time s is located at {right arrow over (x)}, was located at {right arrow over (y)} at reference time s.sub.ref. It is assumed, of course, that {right arrow over (ψ)}(s.sub.ref, {right arrow over (x)})={right arrow over (x)}. Hence, the dynamic attenuation coefficient is given by:
ƒ(s,{right arrow over (x)})=ƒ.sub.0({right arrow over (ψ)}(s,{right arrow over (x)})),  (1.4)

(38) where f.sub.0({right arrow over (x)}) is the spatial distribution of the attenuation coefficient at reference time.

(39) Additionally, the linear interpolation kernel is used for volume grid (φ.sub.h), and the cubic B-spline kernel for the motion grid ({tilde over (φ)}.sub.h). An exact interpolation of the tri-cubic B-spline is also used for volume-to-motion grid transformation.

(40) To describe the motion function, a motion grid is introduced:
{right arrow over (m)}.sub.j=({tilde over (h)}.sub.1j.sub.1,{tilde over (h)}.sub.2j.sub.2,{tilde over (h)}.sub.3j.sub.3),j=(j.sub.1,j.sub.2,j.sub.3), for 0≤j.sub.q<Ñ.sub.q,q=1,2,3.  (1.5)

(41) To enforce smoothness of the motion function, {tilde over (h)}.sub.q>>h.sub.q and Ñ.sub.q>>N.sub.q are selected and c.sub.j(s) is some function of time. To make sure the motion grid (1.5) covers the same region as the reconstruction grid (1.1), it is required that:
{tilde over (h)}.sub.qÑ.sub.q=h.sub.qN.sub.q=: H.sub.q.  (1.6)

(42) Then

(43) $\begin{array}{cc}\stackrel{.fwdarw.}{\phi }\left(s,\stackrel{.fwdarw.}{x}\right)=\underset{i\le j_q<{\tilde{N}}_q}{.Math.}{\stackrel{.fwdarw.}{c}}_j\left(s\right){\tilde{\phi }}_h\left(\stackrel{.fwdarw.}{x}-{\stackrel{.fwdarw.}{m}}_j\right).& \left(1.7\right)\end{array}$

(44) For the data points outside the actual data region, the clamped boundary condition is used, i.e., the data point outside the domain has the same value as the closest boundary value.

(45) In accordance with the present invention, the method and associated algorithm solves the motion estimation problem by finding functions ƒ(s, {right arrow over (x)})) and {right arrow over (φ)}(s, {right arrow over (x)}), so that the following equation holds:
ƒ(s,{right arrow over (x)})=ƒ.sub.0({right arrow over (ψ)}(s,{right arrow over (x)})),  (2.1)

(46) where f.sub.0({right arrow over (x)}) is the volume of the attenuation coefficient at reference time. Introduce a function {right arrow over (u)}(s, {right arrow over (x)}), which is the inverse of {right arrow over (φ)}(s, {right arrow over (x)}). More precisely:
{right arrow over (μ)}(s,{right arrow over (ψ)}(s,{right arrow over (x)}))={right arrow over (x)},{right arrow over (ψ)}(s,{right arrow over (μ)}(s,{right arrow over (x)}))={right arrow over (x)}.  (2.2)
Equation:
{right arrow over (y)}={right arrow over (μ)}(s,{right arrow over (x)})  (2.3)

(47) means that the physical point, which was located at {tilde over (x)} at reference time s.sub.ref, is located at {tilde over (y)} at current time s. In particular, for a fixed {right arrow over (x)}, {right arrow over (u)}(s, {right arrow over (x)}) is just the curve describing the motion of the particle due to breathing. Hence, the derivative:
{right arrow over (v)}(s,{right arrow over (x)})=∂{right arrow over (μ)}(s,{right arrow over (x)})/∂s  (2.4)

(48) is simply the velocity of the particle. Replacing {right arrow over (x)} with {right arrow over (u)}(s, {right arrow over (x)}) in (2.1) and use (2.2) to find:
ƒ(s,{right arrow over (μ)}(s,{right arrow over (x)}))=ƒ.sub.0({right arrow over (x)}),  (2.5)

(49) i.e., the left side of (2.5) is independent of s. Differentiate (2.5) with respect to s and use (2.4) gives:

(50) $\begin{array}{cc}\frac{\partial f\left(s,\stackrel{.fwdarw.}{\mu }\left(s,\stackrel{.fwdarw.}{x}\right)\right)}{\partial s}+\nabla f\left(s,\stackrel{.fwdarw.}{\mu }\left(s,\stackrel{.fwdarw.}{x}\right)\right).Math.\stackrel{.fwdarw.}{\upsilon }\left(s,\stackrel{.fwdarw.}{x}\right)=0.& \left(2.6\right)\end{array}$

(51) Denote {right arrow over (y)}={right arrow over (u)}(s, {right arrow over (x)}). Then, because of (2.2), {right arrow over (x)}={right arrow over (ψ)}(s, {right arrow over (y)}). Using in (2.6) gives:

(52) $\begin{array}{cc}\frac{\partial f\left(s,\stackrel{.fwdarw.}{y}\right)}{\partial s}+f\left(s,\stackrel{.fwdarw.}{y}\right).Math.\stackrel{.fwdarw.}{v}\left(s,\stackrel{.fwdarw.}{\psi }\left(s,\stackrel{.fwdarw.}{y}\right)\right)=0.& \left(2.7\right)\end{array}$

(53) Denote:
{right arrow over (v)}*(s,{right arrow over (y)}):={right arrow over (v)}(s,{right arrow over (ψ)}(s,{right arrow over (y)})).  (2.8)

(54) Then (2.7) simplifies as follows:

(55) $\begin{array}{cc}\frac{\partial f\left(s,\stackrel{.fwdarw.}{y}\right)}{\partial s}+\nabla f\left(s,\stackrel{.fwdarw.}{y}\right).Math.{\stackrel{.fwdarw.}{v}}_{*}\left(s,\stackrel{.fwdarw.}{y}\right)=0.& \left(2.7\right)\end{array}$

(56) Thus, (2.9) is a necessary condition for a function ƒ(s, {right arrow over (y)}) to satisfy (2.5).

(57) Following is a more detailed description of the proposed method of the invention.

(58) Let s.sub.m, 0≤m<M, be the collection of times corresponding to the measured views, and L.sub.mk, 0≤k<K, be the collection of lines corresponding to the mth radiation source position. In other words, L.sub.mk is the line through the mth source position and the kth detector pixel. It is assumed that the detector consists of K pixels. The measured data are denoted G.sub.mk. Finding f(s, {right arrow over (y)}) and {right arrow over (v)}*(s, {right arrow over (y)}) can be done by minimizing the following functional:

(59) $\begin{array}{cc}\Phi \left(f,{\stackrel{.fwdarw.}{v}}_{*}\right)=\frac{1}{2}\underset{\mathrm{mk}}{.Math.}{w_{\mathrm{mk}}\left[G_{\mathrm{mk}}-{\int }_{ℒ\mathrm{mk}}f\left(s_m,\stackrel{.fwdarw.}{x}\right)d\ell \right]}^2+\frac{\lambda }{2}\underset{i}{.Math.}{\int }_0^{\Delta S}{\left[\frac{\partial f\left(s,{\stackrel{.fwdarw.}{x}}_i\right)}{\partial s}+\nabla f\left(s,{\stackrel{.fwdarw.}{x}}_i\right).Math.{\stackrel{.fwdarw.}{v}}_{*}\left(s,{\stackrel{.fwdarw.}{x}}_i\right)\right]}^2\mathrm{ds}+\left(\mathrm{penalty}\mathrm{term}\mathrm{enforcing}\mathrm{smoothness}\mathrm{of}f\right)+\left(\mathrm{penalty}\mathrm{term}\mathrm{enforcing}\mathrm{smoothness}\mathrm{of}{\stackrel{.fwdarw.}{v}}_{*}\right).& \left(2.10\right)\end{array}$

(60) Here λ>0 needs to be a large number to enforce the OFC (2.9). In general, volume image is considered as a non-smooth function, while the velocity is considered to be smooth. To match the characteristics of each image, an edge-preserving regularizer for the volume penalty term (hyperbolic potential), and the Tikhonov regularizer for the velocity penalty term are used.

(61) In a particular embodiment, an alternating minimization is used to demonstrate the performance of the proposed method, but other numerical minimization techniques may be applied to find the optimal solution.

(62) In performing the alternating minimization and dynamic parameter selection, let there be N.sub.p discrete phase bins. For computational efficiency, G.sub.mk is sorted into a group of subsets G.sub.p so that the pth subset corresponds to the data within the pth phase bin, i.e., s.sub.p≤s.sub.m<s.sub.p+1. Let f and v denote phase and velocity volumes unrolled into 1D vectors, respectively. Then (2.10) can be rewritten in a matrix-vector form as follows:
Φ(f,v)=Φ.sub.L(f)+λΦ.sub.fv(f,v)+k.sub.fΦ.sub.R.sub.f(f)+k.sub.vΦ.sub.R.sub.v(v),  (3.1)

(63) where Φ.sub.L(f) is the data fidelity function, Φ.sub.fv(f,v) is the OFC penalty, Φ.sub.R.sub.f(f) is the volume regularizer, and Φ.sub.R.sub.v(v) is the velocity regularizer defined as:

(64) $\begin{array}{cc}\begin{array}{c}{\Phi }_L\left(f\right)\text{=:}\mathrm{ }\frac{1}{2}{\left(G-\mathrm{Af}\right)}^{T}W\left(G-\mathrm{Af}\right)\\ =\frac{1}{2}\underset{p}{.Math.}{\left(G_p-A_p{f}_p\right)}^T{W}_p\left(G_p-A_p{f}_p\right)\end{array}& \left(3.2\right)\\ \begin{array}{c}{\Phi }_{\mathrm{fv}}\left(f,v\right)\text{=:}\mathrm{ }\frac{1}{2}{.Math.\mathrm{Df}+\underset{q}{.Math.}\mathrm{diag}\left\{{\mathrm{Bv}}_q\right\}{C}_{q}f.Math.}_2^2,\\ =\mathrm{ }\frac{1}{2}\underset{p}{.Math.}{.Math.{\left[\mathrm{Df}\right]}_p+\underset{q}{.Math.}\mathrm{diag}\left\{{\mathrm{Bv}}_{p,q}\right\}{C}_q{f}_p.Math.}_2^2,\end{array}& \left(3.3\right)\\ \begin{array}{c}{\Phi }_{R_4}\left(f\right)\text{=:}\mathrm{ }\underset{q}{.Math.}\underset{i}{.Math.}{\varphi }_R\left({\left[C_{q}f\right]}_i,\delta \right),\\ =\mathrm{ }\underset{p}{.Math.}\underset{q}{.Math.}\underset{i}{.Math.}{\varphi }_R\left({\left[C_q{f}_q\right]}_i,\delta \right),\end{array}& \left(3.4\right)\\ \begin{array}{c}{\Phi }_{R_v}\left(v\right)\text{=:}\mathrm{ }\frac{1}{2}\underset{q}{.Math.}\underset{r}{.Math.}{.Math.C_r{v}_q.Math.}_2^2\\ =\mathrm{ }\frac{1}{2}\underset{p}{.Math.}\underset{q}{.Math.}\underset{r}{.Math.}{.Math.C_r{v}_{p,q}.Math.}_2^2,\end{array}& \left(3.5\right)\end{array}$

(65) where

(66) a. G is the entire set of measured data, and G.sub.p is the subset of measured data corresponding to the pth phase bin,

(67) b. f is the entire set of phase volumes in vector form, and f.sub.p is the pth phase volume in vector form,

(68) c. v is the entire set of phase velocity form, and v.sub.p,q is the qth axis component of velocity vector {right arrow over (v)}* at pth phase bin in vector form,

(69) d. A is the entire projection matrix that projects f onto the data domain corresponding to G, and A.sub.p is the projection matrix that projects f.sub.p onto the data domain corresponding to G.sub.p,

(70) e. W is the entire set of weighting matrices, and W.sub.p is the weighting matrix corresponding to G.sub.p,

(71) f. B is the upsampling matrix by cubic B-spline interpolation with the upsampling rate N.sub.q/Ñ.sub.q along the qth axis,

(72) g. D is the temporal finite difference matrix along the phase bins with periodic boundary conditions,

(73) h. C.sub.q is the spatial finite difference matrix along the qth axis with clamped boundaries. Its dimension is determined by the vector to which it applies,

(74) i. diag {t} is the diagonal matrix with the elements of vector t on the main diagonal,

(75) j. [t].sub.i is the ith element (scalar, vector, or matrix) along the first dimension of ector or matrix t.

(76) k. ϕ.sub.R(t, δ) is the element-wise hyperbolic potential defined as:
ϕ.sub.R(t,δ)=δ.sup.2(√{square root over (1+(t/δ).sup.2)}−1),  (3.6)

(77) and δ is the shape parameter defining the “edge”.

(78) Note that even though every linear operation is expressed in matrix form in (3.1)-(3.5), the actual computation does not store those matrices in memory, but instead calculates each element when needed. Also note that Φ.sub.L, Φ.sub.R.sub.f and Φ.sub.R.sub.v are phase-separable, i.e., computing the pth terms does not require to access data G.sub.l or parameters f.sub.l and v.sub.l for all l≠p, which is preferred for multi-device parallel computing. However, Φ.sub.fv is phase-coupled, and it requires accessing neighboring phases. This data exchange latency can be reduced by computing other phase-separable terms while copying data.

(79) The present invention employs a method of alternating minimization, with respect to volume and velocity. The total cost functional in (3.1) is nonlinear, nonconvex and ill-posed. Minimizing (3.1) with respect to volume and velocity simultaneously often falls into undesirable local minima. The method of alternating minimization, in accordance with the present invention, helps to solve the problem by splitting the problem into two simpler sub-problems and guiding the solution to the desired local minimum.

(80) With the alternating minimization scheme, each global iteration step consists of a set of volume sub-iterations and a set of velocity sub-iterations. In each sub-iteration, the cost functional is minimized with respect to f and with respect to v in an alternating fashion while keeping the other variable fixed. This way, each sub-problem becomes convex, which is easier to solve. Also, parameters of the minimization problem are tuned based on the properties of the currently computed solution or using some other consideration. FIG. 2, FIG. 3 and FIG. 4 are flowcharts illustrating a global iteration, volume sub-iteration and motion sub-iteration, respectively.

(81) The method steps 200 of the present invention for global iteration are illustrated with reference to FIG. 2. Assuming the projection data of the object has been previously collected by a scanner, the method begins by sorting the projection data according to the phase 205. At a next step 210, the phase-sorted projection data is loaded into the data processor, and the global iteration counter is set to zero in the following step 215. The method proceeds by computing one or more regularization parameters at step 220. At step 225, a volume sub-iteration is performed followed by a motion sub-iteration at step 230. The method then checks at step 235 to see if the global iteration counter is equal to a total number of global iteration steps. If the total number of global iteration steps has not yet been reached, at step 245, the iteration counter is incremented and steps 225 and 230 are repeated. The process iteratively repeats until the total number of global iteration steps has been met and the motion has been estimated. The estimated motion can then be used to perform a separate motion compensated reconstruction of the image from the projection data. Alternatively, the set of phase volumes obtained at the end of the last volume sub-iteration can be used by the practitioners.

(82) FIG. 3 further details the volume sub-iteration process 225 of the global iteration process 200 for motion estimation. As shown in FIG. 3, the volume sub-iteration process 225 begins at step 305 by resetting Nesterov's parameters. The method continues at step 310 by setting a volume sub-iteration counter to zero. At step 315, a volume preconditioner is computed. At step 320 a volume gradient of data fidelity is computed and at step 325 a volume gradient of regularizer is computed. The method continues at step 330 by computing a volume gradient of the OFC at then performing a Nesterov update at step 335. At step 340, if the volume sub-iteration counter is not equal to a total number of volume sub-iteration steps, the method continues at step 345 by incrementing the volume sub-iteration counter and iteratively repeating the steps of computing a volume preconditioner 315, computing a volume gradient of data fidelity 320, computing a volume gradient of regularizer 325, computing a volume gradient of the volume-velocity coupling 330 and performing a Nesterov update 335 until the volume sub-iteration counter is determined to be equal to the total number of volume sub-iteration steps at step 340. At step 350, the result of the volume sub-iterations is returned after the total number of volume sub-iteration steps has been met.

(83) FIG. 4 further details the motion sub-iteration process 230 of the global iteration process 200 for motion estimation. As shown in FIG. 4, the motion sub-iteration process 230 begins at step 405 by resetting Nesterov's parameters. The method continues at step 410 by setting a motion sub-iteration counter to zero. At step 415 a velocity preconditioner is computed. A velocity gradient of the OFC is computed at step 420 and a velocity gradient of regularizer is computed at step 425. To remove the false constant motion, the voxel-wise phase mean velocity is subtracted in every motion sub-iteration and the velocity at the boundaries of the ROI (region of interest) is set to zero at step 430. A Nesterov update is then performed at step 435 and if a difference between consecutive sub-iteration values of the functional being minimized is not lower than a specified threshold value, at step 445 the method continues by incrementing the motion sub-iteration counter and iteratively repeating the steps of computing a motion preconditioner 415, computing a velocity gradient of the OFC 420, computing a velocity gradient of regularizer 425, subtracting the mean 430 and performing a Nesterov update 435. The method continues until the difference between consecutive sub-iteration values of the functional being minimized is determined to be lower than the specified threshold value at step 440. At step 450, the result of the motion sub-iteration is returned after the total number of specified threshold has been met.

(84) Starting volume sub-iterations with a large λ when the velocity is initialized as zero often results in an underestimation of the motion. In the extreme case, when λ.fwdarw.∞, phase volume images are averaged at the end of the initial set of volume sub-iterations. Once motion information is lost in the phase-volume images due to the averaging effect, it is hard to recover the motion in the following velocity sub-iterations. To appropriately utilize the OFC penalty, a changing λ is required. On the other hand, the data sampling rate of the phase subset data G.sub.p is very low compared to the full angle CT, thus strict phase gating methods often create strong sparse data streak artifacts that adversely affect velocity estimation.

(85) To address this problem, the global iteration process illustrated in FIG. 1 is started with a weak OFC coefficient λ to avoid motion underestimation and strong volume regularizer ϰ.sub.f to reduce sparse view-sampling streaks. Then, λ is gradually ramped-up and ϰ.sub.f is simultaneously reduced to restore details. For example, the following values for λ and ϰ.sub.f may be used at each global iteration step n.sub.g using parameters λ.sub.max, ϰ.sub.max and ϰ.sub.min:

(86) $\begin{array}{cc}\lambda \left(n_g\right)=\min \left({\lambda }_{\max }\frac{n_g}{N_g-1},{\lambda }_{\max }\right),& \left(3.7\right)\\ {\kappa }_f\left(n_g\right)=\max \left({\kappa }_{\max }{2}^{-n_g},{\kappa }_{\min }\right),& \left(3.8\right)\end{array}$

(87) where N.sub.g is the total number of global iteration steps. The estimated velocity is merely an approximation, which may or may not contain errors. To avoid overly-relying on the estimated velocity, λ is set to be no greater than λ.sub.max to limit the maximum volume-velocity constraint strength. Also, ϰ.sub.f is set to be no less than ϰ.sub.min to maintain the denoising property of the regularizer. ϰ.sub.min needs to be chosen depending on the intensity of noise in the measured data. In another embodiment, the strength of the OFC may change non-monotonically during global iterations, e.g. in an oscillatory way.

(88) In the present invention, a preconditioner and Nesterov's momentum acceleration method are applied to each sub-iteration step to increase the convergence speed. No acceleration of global iterations is applied, due to the risk of overshooting the instability, considering the global functional is not convex.

(89) In the volume-iteration of the present invention, an approximated Newton's method is used with separable quadratic surrogate preconditioner. First, the gradient of the cost functional with respect to the volume parameters is found as:
f.sub.p.sup.(n.sup.g.sup.,n.sup.s.sup.+1)=f.sub.p.sup.(n.sup.g.sup.,n.sup.s.sup.)−Ĥ.sub.f.sub.p.sup.−1∇.sub.f.sub.pΦ(f.sup.(n.sup.g.sup.,n.sup.s.sup.),v.sup.(n.sup.g.sup.−1)),  (3.9)

(90) where n.sub.s is the sub-iteration step, v.sup.(n.sup.g.sup.−1) is the result of the velocity sub-iterations in (n.sub.g−1)th global iteration step, ∇.sub.t is the gradient operator with respect to t, and Ĥ.sub.t is the preconditioner for the update of t. Volume gradient of each component of Φ is as follows:

(91) $\begin{array}{cc}{\nabla }_{f_p}{\Phi }_L\left(f\right)=A_p^T{W}_p\left(A_p{f}_p-G_p\right),& \left(3.10\right)\\ \nabla f_p{\Phi }_{\mathrm{fv}}\left(f,v\right)=\left({\left[D^{T}g\right]}_p+\underset{q}{.Math.}{C}_q^T\left(\mathrm{diag}\left\{{\mathrm{Bv}}_{p,q}\right\}{g}_p\right\}\right)& \left(3.11\right)\\ \nabla f_p{\Phi }_{R_f}\left(f\right)=\underset{p}{.Math.}\underset{q}{.Math.}{C}_q^T{\varphi }_R\left(C_q{f}_p,\delta \right),& \left(3.12\right)\end{array}$

(92) where g in (3.11) is defined as:

(93) 0$\begin{array}{cc}g\text{=:}\mathrm{Df}+\underset{q}{.Math.}\mathrm{diag}\left\{{\mathrm{Bv}}_q\right\}{C}_{q}f,& \left(3.13\right)\end{array}$

(94) and ϕ.sub.R(t, δ) is the first derivative of ϕ.sub.R with respect to t.

(95) To design the preconditioner Ĥ.sub.f, the true Hessian must first be found. Let H.sub.f denote the true Hessian of Φ(f,v) with respect to f:

(96) $\begin{array}{cc}{H}_f\text{=:}{\nabla }_f^2\Phi \left(f,v\right)={\nabla }_f^2{\Phi }_L\left(f\right)+\lambda {\nabla }_f^2{\Phi }_{\mathrm{fv}}\left(f,v\right)+{\kappa }_f{\nabla }_f^2{\Phi }_{R_f}\left(f\right),& \left(3.14\right)\\ H_f^L\text{=:}{\nabla }_f^2{\Phi }_L\left(f\right)=A^T\mathrm{WA},& \left(3.15\right)\\ H_f^{\mathrm{fv}}\text{=:}{\nabla }_f^2{\Phi }_{\mathrm{fv}}\left(f,v\right)=D^{T}D+D^T\underset{q}{.Math.}\mathrm{diag}\left\{{\mathrm{Bv}}_q\right\}{C}_q+\underset{q}{.Math.}{C}_q^T\mathrm{diag}\left\{{\mathrm{Bv}}_q\right\}D+\underset{q}{.Math.}{C}_q^T\mathrm{diag}{\left\{{\mathrm{Bv}}_q\right\}}^2{C}_q,& \left(3.16\right)\\ H_f^{R_f}\text{=:}{\nabla }_f^2{\Phi }_{R_f}\left(f\right)=\underset{q}{.Math.}{C}_q^T\mathrm{diag}\left\{{\stackrel{\mathrm{.Math.}}{\varphi }}_R\left(C_{q}f,\delta \right)\right\}{C}_q,& \left(3.17\right)\end{array}$

(97) H.sub.f is not phase-separable, and the size of H.sub.f is very large. Inverting H.sub.f is time-consuming, so it is desirable to have a separable (diagonal) preconditioner for each phase bin volume Ĥ.sub.f.sub.p such that Ĥ.sub.f.sub.pH.sub.f.sub.p to avoid large matrix inversion. The same applies to other Hessians as well:
Ĥ.sub.f.sub.p=: Ĥ.sub.f.sub.p.sup.L+λĤ.sub.f.sub.p.sup.fv+k.sub.fĤ.sub.f.sub.p.sup.R.sup.f.  (3.18)

(98) Beginning with the following estimates:

(99) $\begin{array}{cc}{H}_{f_p}^L⪯\mathrm{diag}\left\{A_p^T{W}_p{A}_p{1}_N\right\},& \left(3.19\right)\\ H_{f_p}^{f_v}⪯{\rho \left(D\right)}^2{I}_N+2\rho \left(D\right)\rho \left(C_q\right)\mathrm{diag}\left\{\mathrm{Bvec}\left\{\underset{q}{.Math.}\underset{i\in {𝒫}_p}{\max }.Math.v_{i,q,j}.Math.\right\}\right\}+{\rho \left(C_q\right)}^2\mathrm{diag}\left\{\mathrm{Bvec}\left\{\underset{q}{.Math.}\underset{i\in {𝒫}_p}{\max }{\left(v_{i,q,j}\right)}^2\right\}\right\}=4I_N+8\mathrm{diag}\left\{\mathrm{Bvec}\left\{\underset{q}{.Math.}\underset{i\in {𝒫}_p}{\max }.Math.v_{i,q,j}.Math.\right\}\right\}+4\mathrm{diag}\left\{\mathrm{Bvec}\left\{\underset{q}{.Math.}\underset{i\in {𝒫}_p}{\max }{\left(v_{i,\mathrm{qJ},j}\right)}^2\right\}\right\},& \left(3.20\right)\\ H_{f_p}^{R_f}⪯3{\rho \left(C_q\right)}^2{I}_N=12I_N,& \left(3.21\right)\end{array}$

(100) where ρ(T) is the spectral radius of a matrix T, vec(t.sub.1) is the unrolled 1D vector of data t in 3D coordinate I, .sub.p is the set of adjacent phase bins .sub.p={l|p−1≤l≤p+1}, 1.sub.N is the all-ones vector of size N=Π.sub.qN.sub.q, and I.sub.N is the identity matrix of size N×N. Storing the preconditioner (3.19) requires a large memory space. Assuming the angular sampling frequency of the phase-gated data G.sub.p is approximately the same for every phase bin, (3.19) can be approximated using the averaged preconditioner:

(101) $\begin{array}{cc}{\hat{H}}_{f_p}^L\text{=:}\frac{2M_p}{M}\mathrm{diag}\left\{\underset{l}{.Math.}{A}_l^T{W}_l{A}_l{1}_N\right\},& \left(3.22\right)\end{array}$

(102) where M.sub.p is the total number of projection data in the pth phase bin. The factor 2 in the numerator is the safety factor in case the angular sampling frequency of each phase bin projection data is not similar to each other.

(103) Also, storing the preconditioner (3.20) requires a considerable amount of memory. To reduce memory storage, (3.20) can be further simplified using a single maximum value:

(104) $\begin{array}{cc}{\hat{H}}_{f_p}^{\mathrm{fv}}\text{=:}\left(4+24\underset{l,q,j}{\max }.Math.v_{l,q,j}.Math.+12\underset{l,q,j}{\max }{\left(v_{l,q,j}\right)}^2\right)I_N.& \left(3.23\right)\end{array}$

(105) Finally, for the volume regularization term is used:
Ĥ.sub.f.sub.p.sup.R.sup.f=: 12I.sub.N.  (3.24)

(106) The volume gradient descent update in (3.9) with the final preconditioner (3.18) is only linearly converging. To speed-up the convergence, Nesterov's momentum acceleration method is applied. If each phase-gated data has significantly biased angular sampling, (3.22) might lead to instability. To avoid the instability, Nesterov acceleration is reset whenever the cost functional value increases during the iteration.

(107) As previously described, in addition to the volume sub-iteration, the present invention also employs a velocity sub-iteration. Similar to the volume sub-iteration step, an approximated Newton's method is used with a separable quadratic surrogate preconditioner. The velocity gradient is obtained as follows:
v.sub.p,q.sup.(n.sup.g.sup.,n.sup.s.sup.+1)=v.sub.p,q.sup.(n.sup.g.sup.,ñ.sup.s.sup.)−Ĥ.sub.v.sub.p,q.sup.−1∇.sub.v.sub.p,qΦ(f.sup.(n.sup.g.sup.−1),v.sup.(n.sup.g.sup.,n.sup.s.sup.)),  (3.25)
∇.sub.v.sub.p,qΦ.sub.fv(f,v)=B.sup.T(diag{C.sub.qf.sub.p}g.sub.p),  (3.26)
∇.sub.v.sub.p,qΦ.sub.R.sub.v(v)=C.sub.q.sup.TC.sub.qv.sub.p,q,  (3.27)

(108) where f.sup.(n.sup.g.sup.−1) is the result of the volume sub-iterations in (n.sub.g−1)th global iteration step. Similar to (3.23) and (3.24), the separable preconditioner for the volume update can be obtained as follows:

(109) $\begin{array}{cc}{\hat{H}}_{v_{p,q}}\text{=:}\lambda {\hat{H}}_{v_{p,q}}^{\mathrm{fv}}+{\kappa }_v{\hat{H}}_{v_{p,q}}^{R_v},& \left(3.28\right)\\ {\hat{H}}_{v_{p,q}}^{\mathrm{fv}}\text{=:}{B}^T\mathrm{diag}\left\{\underset{r}{.Math.}{.Math.C_r{f}_p.Math.}_e\right\}{.Math.C_q{f}_p.Math.}_e,& \left(3.29\right)\\ {\hat{H}}_{v_{p,q}}^{\mathrm{Rv}}\text{=:}12I_{\stackrel{\_}{N}},& \left(3.30\right)\end{array}$

(110) where |t|.sub.e is the element-wise absolute value operator. Nesterov's momentum acceleration is also applied to the velocity gradient update to speed-up the convergence.

(111) While volume sub-iterations do not need to be fully converged at each global iteration step, as details can be restored at later global iterations, stopping velocity sub-iterations while the cost functional is not fully converged typically leads to a loss of small, detailed motions. To avoid this condition, each set of velocity sub-iterations is run until the difference between the values of the cost functional at consecutive sub-iterations is lower than a specified threshold. In contrast, each set of volume sub-iterations is run for a fixed number of iterations.

(112) After a set of motion sub-iterations converges for the first time, motion sub-iterations of the subsequent global iteration steps converge much faster. On the other hand, volume sub-iterations slow down as the global iteration number increases. This is due to the competition between the data fidelity and OFC terms. The latter becomes stronger as the global iterations progress. To accommodate the slower rate of convergence of the volume sub-iterations, the total number of volume sub-iterations may be increased per set by 100 every global iteration.

(113) Motion estimation in the region where there is no object is vulnerable to instability. Although the regularizer prevents the unstable, behavior, the estimated velocity often shows false constant motion. To remove the false constant motion, the voxel-wise phase mean velocity is subtracted in every motion sub-iteration. This mean subtraction process does not improve image quality significantly as most of the false constant motion happens in empty regions, but it helps the convergence speed and stability of the motion iterations as the maximum velocity is reduced. Also, the velocity at the boundaries of the ROI (region of interest) is set to zero, because nothing is moving there. Going further, if the location of a stationary object, such as a patient couch, is known, velocity can be set to zero inside those regions.

(114) In an additional embodiment, motion estimation can be based upon a parameterized cyclic deformation model. In this embodiment, the OFC approach is reformulated using the deformation function instead of the MVFs. Using (2.5), one way to incorporate the OFC (which is an alternative to Equation (2.9)) is via:

(115) $\begin{array}{cc}{\Phi }_{f\mu }\left(f,\stackrel{.fwdarw.}{\mu }\right)=\underset{i}{.Math.}\underset{p}{.Math.}{.Math.f\left(s_p,\stackrel{.fwdarw.}{\mu }\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right)\right)-\frac{1}{M}\underset{p^{\prime }}{.Math.}f\left(s_{p^{\prime },}\stackrel{.fwdarw.}{\mu }\left(s_{p^{\prime }},{\stackrel{.fwdarw.}{x}}_i\right)\right).Math.}^m,& \left(4.1\right)\end{array}$
for some m>0. The choices m=1, 2 are popular ones. Another option is:

(116) $\begin{array}{cc}{\Phi }_{f\mu }\left(f,\stackrel{.fwdarw.}{\mu }\right)=\underset{i}{.Math.}\underset{p}{.Math.}{.Math.f\left(s_{p+1},\stackrel{.fwdarw.}{\mu }\left(s_{p+1},{\stackrel{.fwdarw.}{x}}_i\right)\right)-f\left(s_{p,}\stackrel{.fwdarw.}{\mu }\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right)\right).Math.}^m.& \left(4.2\right)\end{array}$
Eq. (2.9) describes the OFC in the differential form, and Eqs. (4.1), (4.2) describe the OFC in integral form.

(117) The functional Φ.sub.fμ should replace Φ.sub.fv in (2.10). Similarly, the regularization of the velocity term should be replaced by the regularization of the motion function {right arrow over (μ)}. The resulting functional can also be minimized using the alternating minimization approach. It is still linear with respect to f(s, {right arrow over (x)}). Even though the functional is non-linear with respect to {right arrow over (μ)}(s, {right arrow over (x)}), iterations with respect to {right arrow over (μ)}(s, {right arrow over (x)}) are fairly easy and not time consuming, since they do not require access to raw data.

(118) Using the deformation-based approach, several advantages may be realized: 1) It can be enforced that each particle moves periodically along a line segment (or, along an elliptical path). 2) Such a description will have fewer parameters for compute, which will make the overall ME/MC (motion estimation/motion compensation) problem more stable. 3) Since the proposed description has fewer parameters, some of these parameters can be allowed to slowly change during the scan to account for irregular breathing. 4) Initially, for improved stability, it can be assumed that the phase of each voxel is the same as the phase observed by the tracking device on the patient's chest. At later iterations, this requirement can be relaxed. 5) Rigid body constraints can be enforced inside high-attenuating regions, such as bones. This can be achieved for example by thresholding the reconstruction volume (starting after sufficiently many global iterations) and enforcing that neighboring voxels within highly attenuating regions have the same motion function and by adjusting the strength parameter of motion regularization separately at each data point of motion grid based on the corresponding volume to softly encourage the rigid motion at highly attenuating regions. Generally, it is not advisable to add more constraints to the optimization problem, because it may dilute the effect of the data fidelity term. Here, such a problem does not arise, because volume iterations will not see these additional constraints. They will be formulated only in terms of the bone mask (in approach (a)) or in terms of the solution f.sup.(n.sup.g.sup.) at some fixed intermediate iteration n.sub.g (in approach (b)). 6) Use the data fidelity weight w.sub.mk to account for the reliability of the estimated phase of each measure view. 7) Generally, it is easier to incorporate any other physical constraint. 8) Having easy access to the phase information (as part of the deformation function), one can reduce the progressive wave that causes false circular motion by using appropriate regularization, and leave only the stationary wave that is more physical. Similarly to item (5), the requirement that nearby points have approximately the same phase is invisible to volume iterations. 9) An advantage of (4.1) and (4.2) over (2.10) is that the former do not use the derivatives of f, hence the requirements on the smoothness of f are relaxed. This may lead to improved spatial and temporal resolution of the method. 10) The proposed method automatically guarantees the cyclical motion requirement in the reconstruction. 11) The proposed method does not require computing the inverse of the motion {right arrow over (ψ)}={right arrow over (μ)}.sup.−1 during iterations.

(119) In an exemplary embodiment, by adopting the second option (4.2) with m=2 multiplied by the conventional constant ½, we get the OFC constraint

(120) $\begin{array}{cc}{\Phi }_{f\mu }\left(f,\stackrel{.fwdarw.}{\mu }\right)=\underset{i}{.Math.}\underset{p}{.Math.}{.Math.f\left(s_{p+1},\stackrel{.fwdarw.}{\mu }\left(s_{p+1},{\stackrel{.fwdarw.}{x}}_i\right)\right)-f\left(s_{p,}\stackrel{.fwdarw.}{\mu }\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right)\right).Math.}^2.& \left(4.3\right)\end{array}$

(121) The motion model is also chosen as follows, so that each voxel moves along a linear path:

(122) $\begin{array}{cc}\stackrel{.fwdarw.}{\mu }\left(s,\stackrel{.fwdarw.}{x}\right)=\stackrel{.fwdarw.}{x}+\left[\cos \left(2\pi \frac{s-\stackrel{\_}{s}\left(\stackrel{.fwdarw.}{x}\right)}{\Delta S}\right)-\cos \left(2\pi \frac{\stackrel{\_}{s}\left(\stackrel{.fwdarw.}{x}\right)}{\Delta S}\right)\right]\stackrel{.fwdarw.}{u}\left(\stackrel{.fwdarw.}{x}\right),& \left(4.4\right)\end{array}$

(123) where amplitude {right arrow over (u)}({right arrow over (x)}) and phase {right arrow over (s)}({right arrow over (x)}) are given by

(124) 0$\begin{array}{cc}\stackrel{.fwdarw.}{u}\left(\stackrel{.fwdarw.}{x}\right)=\underset{j}{.Math.}{\stackrel{.fwdarw.}{u}}_j\hat{\phi }\left(\stackrel{.fwdarw.}{x}-{\stackrel{.fwdarw.}{m}}_j\right),& \left(4.5\right)\\ \stackrel{\_}{s}\left(\stackrel{.fwdarw.}{x}\right)=\underset{k}{.Math.}{\stackrel{\_}{s}}_k\stackrel{\bigvee }{\phi }\left(\stackrel{.fwdarw.}{x}-{\stackrel{.fwdarw.}{n}}_k\right).& \left(4.6\right)\end{array}$

(125) Note that the grids {right arrow over (x)}.sub.i, {right arrow over (m)}.sub.j, and {right arrow over (n)}.sub.k do not have to be the same, as well as the interpolation functions φ(.Math.), {circumflex over (φ)}(.Math.) and {hacek over (φ)}(.Math.) do not have to be the same. In general, {right arrow over (u)} is smooth, and s is even smoother. Thus, a sparse grid can be used with a smooth interpolation function for {right arrow over (u)}, and even more sparse grid and a smoother interpolation function for s.

(126) After replacing (3.3) with (4.3) and adding regularization terms for the motion parameters u and s, the total functional is as follows:
Φ(f,u,s)=Φ.sub.L(f)+λΦ.sub.f.sub.μ(f,u,s)+k.sub.fΦ.sub.R.sub.f(f,δ)+k.sub.uΦ.sub.R.sub.u(u)+k.sub.sΦ.sub.R.sub.s(s),  (4.7)

(127) For Φ.sub.R.sub.u and Φ.sub.R.sub.s, the Tikhonov regularization is used, analogously to Φ.sub.R.sub.v in the previous section.

(128) As in the previous section, an approximated Newton's method with separable surrogate function is used. Here, only the gradient descent part (computation of the gradient and Hessian) of the new motion model is described. All the other terms, such as data fidelity and volume regularization terms, are the same as in the previously described method.

(129) Hereinafter, the following shorthand notations are used for brevity:
{right arrow over (μ)}.sub.p,i={right arrow over (μ)}(s.sub.p,{right arrow over (x)}.sub.i),
μ.sub.p,i,q=μ.sub.q(s.sub.p,{right arrow over (x)}.sub.i),
u.sub.(i),q=u.sub.q({right arrow over (x)}.sub.i),
s.sub.(i)=s({right arrow over (x)}.sub.i),
.sub.p,1=f(s.sub.p,{right arrow over (μ)}.sub.p,1).

(130) In the motion sub-iteration step, the gradient of Φ.sub.f.sub.μ with respect to the motion vector parameter u is obtained as follows:

(131) $\begin{array}{cc}\frac{\partial {\Phi }_{f_{\mu }}}{\partial u_{j,q}}=\underset{i\in \hat{I}\left(j\right)}{.Math.}\underset{p}{.Math.}\left[\left({ℱ}_{p+1,i}-{ℱ}_{p,1}\right)\left(\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{\partial {\mu }_{p+1,i,q}}{\partial u_{j,q}}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{\partial {\mu }_{p,i,q}}{\partial u_{j,q}}\right)\right],& \left(4.8\right)\end{array}$
where {hacek over (I)}(j)={i|{circumflex over (φ)}({right arrow over (x)}.sub.i−{right arrow over (m)}.sub.j)≠0}, and

(132) $\begin{array}{cc}\frac{\partial {\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_k}=\frac{2\pi }{\Delta S}\left(\sin \left(2\pi \frac{s_p-{\stackrel{\_}{s}}_{\left(i\right)}}{\Delta S}\right)+\sin \left(2\pi \frac{{\stackrel{\_}{s}}_{\left(i\right)}}{\Delta S}\right)\right)u_{\left(i\right),q}.Math.\stackrel{\bigvee }{\phi }\left({\stackrel{.fwdarw.}{x}}_i-{\stackrel{.fwdarw.}{n}}_k\right).& \left(4.9\right)\end{array}$

(133) With bilinear interpolation, the analytic form of

(134) $\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}$
cannot be computed, as φ is not differentiable. Instead, we compute it numerically by applying the central difference first, then find the value at μ.sub.p+1,i,q. Similarly, the gradient of Φ.sub.f.sub.μ with respect to the motion phase parameters is obtained as follows:

(135) $\begin{array}{cc}\frac{\partial {\Phi }_{f_{\mu }}}{\partial {\stackrel{\_}{s}}_k}=\underset{i\in \stackrel{\bigvee }{I}\left(k\right)}{.Math.}\underset{p}{.Math.}[\left({ℱ}_{p+1,i}-{ℱ}_{p,i}\right).Math.\underset{q}{.Math.}\left(\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{\partial {\mu }_{p+1,i,q}}{\partial {\stackrel{\_}{s}}_k}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{\partial {\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_k}\right)],& \left(4.10\right)\end{array}$
where {hacek over (I)}(k)={i|{circumflex over (φ)}({right arrow over (x)}.sub.i−{right arrow over (m)}.sub.k)≠0}, and

(136) $\begin{array}{cc}\frac{\partial {\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_k}=\frac{2\pi }{\Delta S}\left(\sin \left(2\pi \frac{s_p-{\stackrel{\_}{s}}_{\left(i\right)}}{\Delta S}\right)+\sin \left(2\pi \frac{{\stackrel{\_}{s}}_{\left(i\right)}}{\Delta S}\right)\right)u_{\left(i\right),q}.Math.\stackrel{\bigvee }{\phi }\left({\stackrel{.fwdarw.}{x}}_i-{\stackrel{.fwdarw.}{n}}_k\right).& \left(4.11\right)\end{array}$

(137) The Hessian of Φ.sub.f.sub.μ is obtained as follows:

(138) $\begin{array}{cc}\frac{{\partial }^2{\Phi }_{f_{\mu }}}{\partial u_{j,q}\partial u_{j^{\prime },q^{\prime }}}=\underset{i\in \hat{I}\left(j\right).Math.\hat{I}\left(j^{\prime }\right)}{.Math.}\underset{p}{.Math.}\left[\left(\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{\partial {\mu }_{p+1,i,q}}{\partial u_{j,q}}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{\partial {\mu }_{p,i,q}}{\partial u_{j,q}}\right).Math.\left(\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q^{\prime }}}\frac{\partial {\mu }_{p+1,i,q^{\prime }}}{\partial u_{j^{\prime },q^{\prime }}}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q^{\prime }}}\frac{\partial {\mu }_{p,i,q^{\prime }}}{\partial u_{j^{\prime },q^{\prime }}}\right)+\left({ℱ}_{p+1,i}-{ℱ}_{p,i}\right).Math.\left(\frac{{\partial }^2{ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}\partial {\mu }_{p+1,i,q^{\prime }}}\frac{\partial {\mu }_{p+1,i,q}}{\partial u_{j,q}}\frac{\partial {\mu }_{p+1,i,q^{\prime }}}{\partial u_{j^{\prime },q^{\prime }}}-\frac{{\partial }^2{ℱ}_{p,i}}{\partial {\mu }_{p,i,q}\partial {\mu }_{p,i,q^{\prime }}}\frac{\partial {\mu }_{p,i,q}}{\partial u_{j,q}}\frac{\partial {\mu }_{p,i,q^{\prime }}}{\partial u_{j^{\prime },q^{\prime }}}\right)\right].& \left(4.12\right)\\ \frac{{\partial }^2{\Phi }_{f_{\mu }}}{\partial {\stackrel{\_}{s}}_k\partial {\stackrel{\_}{s}}_{k^{\prime }}}=\underset{i\in \stackrel{\bigvee }{I}\left(k\right).Math.\stackrel{\bigvee }{I}\left(k^{\prime }\right)}{.Math.}\underset{p}{.Math.}\left[\left(\underset{q}{.Math.}\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{\partial {\mu }_{p+1,i,q}}{\partial {\stackrel{\_}{s}}_k}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{\partial {\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_k}\right).Math.\left(\underset{q}{.Math.}\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{\partial {\mu }_{p+1,i,q}}{\partial {\stackrel{\_}{s}}_{k^{\prime }}}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{\partial {\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_{k^{\prime }}}\right)+\left({ℱ}_{p+1,i}-{ℱ}_{p,i}\right).Math.\left(\underset{q^{\prime }}{.Math.}\left(\frac{{\partial }^2{ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}\partial {\mu }_{p+1,i,q^{\prime }}}\frac{\partial {\mu }_{p+1,i,q}}{\partial {\stackrel{\_}{s}}_k}\frac{\partial {\mu }_{p+1,i,q^{\prime }}}{\partial {\stackrel{\_}{s}}_{k^{\prime }}}-\frac{{\partial }^2{ℱ}_{p,i}}{\partial {\mu }_{p,i,q}\partial {\mu }_{p,i,q^{\prime }}}\frac{\partial {\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_k}\frac{\partial {\mu }_{p,i,q^{\prime }}}{\partial {\stackrel{\_}{s}}_{k^{\prime }}}\right)+\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{{\partial }^2{\mu }_{p+1,i,q}}{\partial {\stackrel{\_}{s}}_k\partial {\stackrel{\_}{s}}_{k^{\prime }}}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{{\partial }^2{\mu }_{p,i,q}}{\partial {\stackrel{\_}{s}}_k\partial {\stackrel{\_}{s}}_{k^{\prime }}}\right)\right].& \left(4.13\right)\\ \frac{{\partial }^2{\Phi }_{f_{\mu }}}{\partial u_{j,q}\partial {\stackrel{\_}{s}}_k}=\underset{i\in \hat{I}\left(j\right).Math.\stackrel{\bigvee }{I}\left(k\right)}{.Math.}\underset{p}{.Math.}\left[\left(\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{\partial {\mu }_{p+1,i,q}}{\partial u_{j,q}}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{\partial {\mu }_{p,i,q}}{\partial u_{j,q}}\right).Math.\underset{q^{\prime }}{.Math.}\left(\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q^{\prime }}}\frac{\partial {\mu }_{p+1,i,q^{\prime }}}{\partial {\stackrel{\_}{s}}_k}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q^{\prime }}}\frac{\partial {\mu }_{p,i,q^{\prime }}}{\partial {\stackrel{\_}{s}}_k}\right)+\left({ℱ}_{p+1,i}-{ℱ}_{p,i}\right).Math.\left(\underset{q^{\prime }}{.Math.}\left(\frac{{\partial }^2{ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}\partial {\mu }_{p+1,i,q^{\prime }}}\frac{\partial {\mu }_{p+1,i,q}}{\partial u_{j,q}}\frac{\partial {\mu }_{p+1,i,q^{\prime }}}{\partial {\stackrel{\_}{s}}_k}-\frac{{\partial }^2{ℱ}_{p,i}}{\partial {\mu }_{p,i,q}\partial {\mu }_{p,i,q^{\prime }}}\frac{\partial {\mu }_{p,i,q}}{\partial u_{j,q}}\frac{\partial {\mu }_{p,i,q^{\prime }}}{\partial {\stackrel{\_}{s}}_k}\right)+\frac{\partial {ℱ}_{p+1,i}}{\partial {\mu }_{p+1,i,q}}\frac{{\partial }^2{\mu }_{p+1,i,q}}{\partial u_{j,q}\partial {\stackrel{\_}{s}}_k}-\frac{\partial {ℱ}_{p,i}}{\partial {\mu }_{p,i,q}}\frac{{\partial }^2{\mu }_{p,i,q}}{\partial u_{j,q}\partial {\stackrel{\_}{s}}_k}\right)\right].& \left(4.14\right)\end{array}$

(139) Elements of the Hessian matrix (4.12)-(4.14) are not pixel-independent. One can create a diagonal upper bound matrix (surrogate function) by summing the absolute values of all the coupled terms.

(140) In the volume sub-iteration of the new method, the gradient of Φ.sub.f.sub.μ with respect to the volume parameter f is obtained as follows:

(141) $\begin{array}{cc}\frac{\partial {\Phi }_{f_{\mu }}}{\partial f_{p,i}}=\underset{t\in T\left(p,i\right)}{.Math.}\left(-{ℱ}_{p+1,t}+2{ℱ}_{p,t}-{ℱ}_{p-1,t}\right)\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}},\mathrm{where}T\left(p,i\right)=\left\{t|\phi \left({\stackrel{.fwdarw.}{\mu }}_{p,t}-{\stackrel{.fwdarw.}{x}}_i\right)\ne 0\right\},\mathrm{and}& \left(4.15\right)\\ \frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}}=\phi \left({\stackrel{.fwdarw.}{\mu }}_{p,t}-{\stackrel{.fwdarw.}{x}}_i\right).& \left(4.16\right)\end{array}$

(142) The Hessian of Φ.sub.f.sub.μ with respect to the volume parameter f is obtained as follows:

(143) $\begin{array}{cc}\frac{{\partial }^2{\Phi }_{f_{\mu }}}{\partial f_{p,i}\partial f_{p^{\prime },i^{\prime }}}=\{\begin{array}{cc}\underset{t\in T\left(p,i\right).Math.T\left(p,i^{\prime }\right)}{.Math.}2\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}}\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i^{\prime }}}& \mathrm{if}{p}^{\prime }=p\\ \underset{t\in T\left(p,i\right).Math.T\left(p^{\prime },i^{\prime }\right)}{.Math.}-\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}}\frac{\partial {ℱ}_{p^{\prime },t}}{\partial f_{p^{\prime },i^{\prime }}}& \mathrm{if}{p}^{\prime }=p-1\mathrm{or}p+1\\ 0& \mathrm{otherwise},\end{array}& \left(4.17\right)\end{array}$

(144) To find the separable surrogate function, the following interpolation properties can be used:

(145) $\begin{array}{cc}0⩽\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}}⩽1\forall t,& \left(4.18\right)\\ \underset{i\in I\left(p,t\right)}{.Math.}\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}}=1\forall t.& \left(4.19\right)\end{array}$

(146) A diagonal upper bound of the Hessian can be found using the interpolation properties (4.18) and (4.19) as follows:

(147) 0$\begin{array}{cc}\begin{array}{c}{\left[\frac{{\partial }^2{\Phi }_{f_{\mu }}}{\partial f_{p,i}\partial f_{p^{\prime },i^{\prime }}}\right]}_{\left(p\times i,p^{\prime }\times i^{\prime }\right)}\preccurlyeq \mathrm{diag}\left\{\underset{p}{.Math.}\underset{p^{\prime }}{.Math.}\underset{i^{\prime }}{.Math.}\left|\frac{{\partial }^2{\Phi }_{f_{\mu }}}{\partial f_{p,i}\partial f_{p^{\prime },i^{\prime }}}\right|\right\}\\ =\mathrm{diag}\left\{4\underset{t}{.Math.}\frac{\partial {ℱ}_{p,t}}{\partial f_{p,i}}\right\}.\end{array}& \left(4.20\right)\end{array}$

(148) Note that not all interpolators satisfy the non-negativity property (4.18). One example of such interpolator is the high order polynomial interpolator. The Hessian with interpolators that do not satisfy (4.18) may not be diagonalized using (4.20).

(149) For massively parallel computing, there are two options to compute (4.15) and (4.20): a. Assign threads along t and loop over all i ∈I(p,t), where I(p,t)={i|φ({right arrow over (μ)}.sub.p,t−{right arrow over (x)}.sub.i)≠0}, b. Assign threads along i and loop over all t ∈T(p,i).

(150) Finding I(p,t) for each t is easier than finding T(p,i) for each i, because {right arrow over (μ)}.sub.p,t is fixed for each thread. On the other hand, approach (b) requires searching along t to find all {right arrow over (μ)}.sub.p,t such that t ∈T(p,i). Further more, the size of the support of φ is known prior to computation and in general small, so each t-th thread only needs to access a few f.sub.p,is.

(151) However, approach (a) can cause so called, “race condition” as different threads need to read and update the same i location simultaneously. Fortunately, this race condition can be avoided by using atomic operations in CUDA, with some sacrifice of the parallelism. Even better, the race condition does not happen very often when the bilinear interpolation is used as the size of the support is fairly small. Therefore, approach (a) has been chosen to compute (4.15) with bilinear interpolation.

(152) While in the previously described embodiment, the strength of the OFC λ was monotonically increased along the global iterations, it was found that this caused the motion to be consistently underestimated. As such, in this embodiment, instead of a monotonic increase, the OFC strength is changed to oscillate throughout the global iterations with the period of two iterations, i.e. λ=1 for odd iterations, and λ=150 for even iterations.

(153) Odd iterations add more details to the volumes and estimate the motion based on the added details, and even iterations impose the motion on the volumes more strictly. This way, global iterations can keep updating the motion even when the motion is underestimated in earlier iterations.

(154) Additionally, in this embodiment, instead of decreasing the volume regularization strength parameter ϰ.sub.f, the volume regularization shape parameter δ.sub.f is decreased, thereby providing two effects, reduction of regularization strength and transition from Tikhonov type regularization to TV-like regularization as the iterations continue. FIG. SA illustrates the assigned values of λ, ϰ.sub.f and δ.sub.f for each global iteration step in the previously described embodiment, wherein the parameters are increased monotonically. FIG. 5B illustrates the assigned values of λ, ϰ.sub.f and δ.sub.f for each global iteration step in the currently described embodiment, wherein the parameters are increased non-monotonically.

(155) While the results for both discrete and continuous motion data are promising, the estimated motion in the continuous case creates slight errors at the boundary. The reason is that even if the reconstructed motion matches the average motion within each phase bin, some projections exceed the deformation range of the estimated motion. As a result, the motion is slightly underestimated and the projections outside the estimate motion create slight streak artifacts in the air account the body being imaged.

(156) The results show that when there is inconsistency of motion within projection data in each phase bin, more studies may be required to achieve the best result. The results also show that artifacts due to replacing continuous motion with discrete motion in simulated experiments are likely negligible compared with artifacts due to significant imperfections in real data, such as breathing motion is not periodic, imperfect binning of data, etc.

(157) Additionally, while the exemplary method and associated algorithm work will with simulated data when there is no residual motion and no noise, applying the same algorithm to clinical data may lead to undesirable speckle artifacts. While the OFC encourages incompressibility, if the volume is indeed incompressible, the data density on the volume grid should be unchanged in every phase. However, the estimated volume usually shows some compressibility (whether it is real or not), and the data density at the expanding region becomes sparse, causing the OFC to be applied unevenly. Consequently, the OFC fails to regularize the expanding or compressing regions properly, resulting in speckle artifacts when the data are noisy and not perfectly matching.

(158) To address these speckle artifacts, MVFs between pairs of consecutive phase volumes on a regular grid are estimated from the motion function {right arrow over (μ)}(s.sub.p, {right arrow over (x)}.sub.i) as follows:

(159) $\begin{array}{cc}{\stackrel{.fwdarw.}{V}}_p=\underset{\stackrel{\_}{\upsilon }}{\arg \min }\frac{1}{2}\underset{i}{.Math.}|\stackrel{.fwdarw.}{\upsilon }\left(\stackrel{.fwdarw.}{\mu }\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right)\right)-\left(\stackrel{.fwdarw.}{\mu }\left(s_{p+1},{\stackrel{.fwdarw.}{x}}_i\right)-\stackrel{.fwdarw.}{\mu }\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right)\right){|}^2& \left(5.1\right)\end{array}$

(160) Then the estimated regular grid MVFs is used to impose the OFC in an even fashion (the magnitudes of MVFs are small, so with one grid being regular, its image at the next phase volume is close to a regular grid as well):

(161) $\begin{array}{cc}{\Phi }_{f_{\mu }}\left(f,\stackrel{.fwdarw.}{\mu }\right)=\frac{1}{2}\underset{i}{.Math.}\underset{p}{.Math.}|f\left(s_{p+1},{\stackrel{.fwdarw.}{x}}_i+{\stackrel{.fwdarw.}{V}}_p\left({\stackrel{.fwdarw.}{x}}_i\right)\right)-f\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right){|}^2.& \left(5.2\right)\end{array}$

(162) In the original OFC, pairs of consecutive phase volumes are compared using a grid, that was regular in the reference phase. This grid is far from regular after deformation:

(163) $\begin{array}{cc}{\Phi }_{f_{\mu }}\left(f,\stackrel{.fwdarw.}{\mu }\right)=\frac{1}{2}\underset{i}{.Math.}\underset{p}{.Math.}|f\left(s_{p+1},\stackrel{.fwdarw.}{\mu }\left(s_{p+1},{\stackrel{.fwdarw.}{x}}_i\right)\right)-f\left(s_p,\stackrel{.fwdarw.}{\mu }\left(s_p,{\stackrel{.fwdarw.}{x}}_i\right)\right){|}^2.& \left(5.3\right)\end{array}$

(164) Accordingly, is this embodiment, the OFC is modified to use MVFs on a regular grid to remove the speckle artifacts.

(165) The present invention may be embodied on various computing platforms that perform actions responsive to software-based instructions. The following provides an antecedent basis for the information technology that may be utilized to enable the invention.

(166) The computer readable medium described in the claims below may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain or store a program for use by or in connection with an instruction execution system, apparatus, or device.

(167) A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system. The present invention may be embodied on various computing platforms that perform actions responsive to software-based instructions. The following provides an antecedent basis for the information technology that may be utilized to enable the invention.

(168) The computer readable medium described in the claims below may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain or store a program for use by or in connection with an instruction execution system, apparatus, or device.

(169) Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wire-line, optical fiber cable, radio frequency, etc., or any suitable combination of the foregoing. Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, C#, C++ or the like and conventional procedural programming languages, such as the “C” programming language or similar programming languages.

(170) Aspects of the present invention are described below with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

(171) These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.

(172) The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

(173) It will be seen that the advantages set forth above, and those made apparent from the foregoing description, are efficiently attained and since certain changes may be made in the above construction without departing from the scope of the invention, it is intended that all matters contained in the foregoing description or shown in the accompanying drawings shall be interpreted as illustrative and not in a limiting sense.

(174) It is also to be understood that the following claims are intended to cover all of the generic and specific features of the invention herein described, and all statements of the scope of the invention which, as a matter of language, might be said to fall therebetween.