Controlled Cardiac Computed Tomography

Cardiac computed tomography (CT) has been a hot topic for years because of the clinical importance of cardiac diseases and the rapid evolution of CT systems. In this paper, we propose a novel strategy for controlled cardiac CT that may effectively reduce image artifacts due to cardiac and respiratory motions. Our approach is radically different from existing ones and is based on controlling the X-ray source rotation velocity and powering status in reference to the cardiac motion. We theoretically show that by such a control-based intervention the data acquisition process can be optimized for cardiac CT in the cases of periodic and quasiperiodic cardiac motions. Specifically, we formulate the corresponding coordination/control schemes for either exact or approximate matches between the ideal and actual source positions, and report representative simulation results that support our analytic findings.


INTRODUCTION
Cardiovascular disease is the number one killer in the Western world. It is responsible for 1 of every 2.6 deaths in the United States. Cardiovascular disease was listed as a primary or contributing cause on about 1 408 000 death certificates annually. One in five persons has some form of cardiovascular diseases in the United States, account for 64 400 000 Americans (http://www.americanheart.org/). In 2004, the estimated cost of cardiovascular diseases is $368.4 billion. Needless to say, early detection of diseases is vitally important.
Electron-beam CT (EBCT) is a unique mode of medical X-ray CT, which is for early screening of coronary artery diseases [1][2][3][4]. EBCT is more accurate than relying on some traditional indicators, for example, high cholesterol levels, high blood pressure, obesity, diabetes, and so on, for determining an individual's risk for heart disease [5][6][7][8][9][10]. Specifically, EBCT can quantify calcification in the arteries, which is a sign of atherosclerosis. An EBCT-based coronary calcium score is, according to the American Heart Association (http://www.americanheart.org/), of significant predictive value of heart attacks or need for bypass or angioplasty over the next year or two. In addition, the American Heart Association has noted that EBCT is also useful in evaluating bypass graft potency, intra and congenital cardiac lesions, and quantifying ventricular muscle mass, chamber volumes, and systolic and diastolic functions. Hence, it has been recognized as the gold standard for detecting early heart diseases. In addition to its remarkable applications for dynamic anatomical imaging of cardiac structures, EBCT is also a powerful tool for physiological imaging [11]. However, there are three major limitations with the current EBCT techniques. First, it is not in cone-beam geometry, and hence it can only acquire a limited number of transverse slices through the heart. It has become clear that the multislice/cone-beam scanning is much more advantageous for a good portion or even the whole volume of the heart to be contained in a cone-beam, and reconstructed from a rapidly acquired data set for both high temporal resolution and excellent anatomical consistency. Second, the X-ray spot in EBCT is not sufficiently intensive to produce the image quality that the mechanical rotation-based scanners can achieve, which compromises contrast resolution in reconstructed images. Third, the EBCT scanner is both expensive and monstrous. As such, EBCT is far less accessible and less cost-effective than the main stream medical CT scanner that utilizes a mechanically rotated X-ray source.
Given the limitations of the EBCT, multislice/cone-beam CT has been growing into a strong competitor of the EBCT. Research on single-slice cardiac CT was initiated in [12,13]. Multislice/cone-beam cardiac CT has been actively studied since then [14][15][16][17][18][19][20][21]. To obtain usable images of the cardiac or coronary features, a CT scan must be performed in coordination with an ECG signal or equivalent. There are two types of synchronization mechanisms: prospective and retrospective gating [22]. In the prospective mode, the heart is scanned at a pre-specified delay relative to a temporal landmark, such as after the peak of an R wave. Multiple prospectively triggered acquisitions can reveal different phases of the cardiac motion. In the retrospective mode, the heart is continuously scanned while the ECG is simultaneously recorded. Then, images are reconstructed from a retrospectively selected cardiac phase. Most of these cardiac imaging methods rely on retrospective gating to lock into consistent phases of the cardiac motion. An alternative to the ECG is referred to as the kymogram, which is based on an analysis on the trajectory of the center of the mass of the beating heart [18]. A knowledge-based cardiac CT algorithm was also developed [21]. Advantages of the cardiac phase-based approach clearly outperformed the classic half-scan methods. The manufacturers have now already switched from half-scan or partial-scan to multisector reconstruction algorithms for cardiac CT.
There are however two major problems with current gating-based cardiac CT methods [3,[23][24][25]. First, these existing methods are passive in their nature, and require that the cardiac motion and the X-ray source rotation must be at favorable relative frequencies; otherwise, the data sectors to be assembled for a complete data set would span a wide range of the cardiac status leading to a compromised image quality, or the data acquisition time would be too long to be practical (in a most unfavorable case, it could be impossible to acquire a complete data set) [17,21]. Second, even in the favorable cases, retrospectively reconstructed cardiac images still suffer from substantial motion blurring because in practice each projection sector covers a projection angular range of a substantial length. Within such an angular range, the heart will move appreciably, especially when it is not in a relative stationary phase. As a benchmark, we routinely achieve 0.3 mm spatial resolution in spiral CT of the temporal bone where the motion magnitude is much less than that of the heart [21,26]. On the other hand, the spatial resolution with cardiac CT is at best in the millimeter domain.
It is also important to comment that to understand etiology and pathogenesis of the human cardiovascular diseases as well as to develop prevention and treatment strategies, small animals have become common laboratory models. To use these animal models fully, the extension of cardiovascular imaging from patients to small animals, such as mice and rats, is imperative [27]. Because of the small size and high heart rate of a mouse, high spatial and temporal resolutions are required. Currently, an X-ray micro-CT scanner takes at least 20 seconds to acquire a full data set [27,28]. Hence, cardiac micro-CT represents a daunting challenge. The recent small animal micro-CT work performed at Duke University showed the feasibility of respiratory and cardiacgated micro-CT [29], but their system runs very slowly (not well suited for high throughput [29]), and requires that the mouse is placed vertically and rotated asynchronously, violating the natural physiological conditions.
All the above problems demonstrate themselves as motion-induced image artifacts. To this end, we propose a new approach, which is radically different from any existing technique. Our primary insight is that cardiac imaging can be optimally performed by integrating modern control and imaging techniques, that is, by adaptively controlling the imaging device according to the real-time analysis and prediction of an individualized respiratory and cardiac motion functions, little artifacts would be left. It is expected that the spatial resolution with the proposed controlled cardiac CT should be significantly better than that with current cardiac CT, and eventually made to be comparable to that with temporal bone CT.
Two questions come to the mind immediately. Is controlling source rotation velocity of a CT a realistic task and how to coordinate prediction and control? This paper discusses the prediction and coordination questions, and sets the theoretical foundation for controlled cardiac CT/micro-CT. We only make a comment on the control hardware. Control technology is a matured area and a number of advanced control methods, including robust, adaptive, nonlinear, fuzzy, and intelligent controls, have been proposed, theoretically justified and practically tested. One can easily find a precision control including hardware and software everywhere in every day life from airplanes, electronics to microsurgery. A CT with controllable source rotation velocity is completely feasible which will be demonstrated in the paper.

PROBLEM STATEMENTS
To describe coordination of prediction and control in a precise term, we define a few terminologies.
(i) Let v(t) be the volume or phase of the object, for example, a heart, to be scanned at time t. It is a function of time and could be periodic or nonperiodic. (ii) Let a l = (l−1)/ p·2π, l = 1, 2, . . . , p, be p evenly spaced angles. Evenly spaced is for simplicity only. Non-evenly spaced angles can be considered similarly. (iii) Let s(t) represent the angle, at time t, of the CT X-ray source. It is a function of time and its rotational velocity is assumed to be controllable.
In this paper, we assume that at any given time, only one projection can be taken. The results can be easily extended to multiple source/sensor cases.
First consider a simple case of a periodic v(t) with period P, that is, Define the volume levels of the object and the corresponding two sequences Erwei Bai et al. Here, it is recognized that for the same cardiac volume level, the structure of the heart could be different in the expanding and contracting phases.
Volume levels are the levels where the object is required to be imaged. Figure 1 shows the interaction of v(t i )'s and a l 's for q = 2 and p = 3. Now, for the prescribed angles a l 's and the levels v(t i )'s, solving the problem of CT imaging without motion interference is equivalent to solving the following problem.
In other words, the source angle is a l when the volume level is at v(t i ) for all i, l. In short, the source covers all the angles a l , l = 1, . . . , p, exactly at all the levels v(t i ), i = 1, . . . , q.
(b) Among all the solutions s(t)'s, find one that requires the minimum time to cover all the angles for all the levels.
In general, a constant source rotation velocity s(t) = αt, for some constant α, does not solve Problem 1.

Periodic v(t)
To solve Problem 1, we now consider controlled source rotation velocity s(t). Initially, intuition seems to suggest that a solution of Problem 1 would need a very fast and complicated profile of s(t). This is however not the case. To this end, we make two observations. (i) At each level i, p projections at p different angles a l , l = 1, . . . , p, have to be taken. (ii) Only times that projections can be taken exactly at the level i are Because of p angles at each level and the assumption that at any given time, only one projection can be made, the minimum time needed for the level i is t i + (p − 1)P. This leads to the following lemma.

Lemma 1. The minimum time needed to cover all the angles at all the levels for any solution s(t) is at least
This is a lower bound on the minimum time. Any solution s(t) of Problem 1(a) achieving the lower bound (6) solves Problem 1(b). Base on these observations, we have the following Theorem 1. Up to a possible permutation on a l , a necessary and sufficient condition for s(t) to be a solution of Problem 1(b) is This is an interpolation problem and many solutions, for example, polynomial or spline interpolations, exist. To compare the results with the gating method for a constant source rotation velocity s(t), we give a numerical example. First note the gating method is to specify the X-ray source rotation velocity and data acquisition timing based on the cardiac rate so that data acquisition timing is at the consistent heart volume levels while the X-ray source is distributed in an appropriate angular range. Now, let the heart volume be v(t) = r 1 + r 2 cos(2πt) ( 8 ) for some constants r 1 , r 2 . v(t) is periodic with period P = 1 (sec). Suppose the level that we are interested in is at t 1 = 5/8 or v t 1 = r 1 + r 2 cos 2π 5 8 + m .
Since v(t) is periodic, v(t 1 + m) = v(t 1 ) for any integer m. Now, without loss of generality, let the angles be a l = 2π(l − 1)/3, l = 1, 2, 3. To have exact angles a l = 2π(l − 1)/3 at the level v(t 1 ), we must have for some m's, It can be shown that there does not exist a constant source rotation velocity s(t) = αt, so the above three equations can be satisfied simultaneously. Approximate solutions do exist.
To quantify the error, we define the minimum error between the given angles and all possible X-ray source angles at the given heart volume level v(t 1 ) as Table 1, converted to degrees, shows e and the corresponding optimal α * and m * . For instance, if the error is required to be no larger than 0.631 • , the best α ∈ (0, ∞) is 0.316π that results in the minimum error 0.631 • at time t = 82+5/8 (sec).  Clearly, the more accurate the solution is required, the longer time the scan has to take. Now, let the X-ray source velocity be controllable as defined by with It is easily verified that this X-ray source s(t) covers three angles exactly at the required volume level v(5/8 + m) for any integer k > 0, that is, s(5/8 + 3k) = 0(+2π), s(5/8 + 4k) = 2 * pi/3(+2π), and s(5/8 + 5k) = 4 * π/3(+2π). k balances the time required and the maximum velocity and acceleration. A small k gives rise to a small maximum velocity and acceleration but at the same time a longer time to cover three angles. A large k results in a short time with a price of a large At maximum velocity and acceleration. Figure 2 shows the relationship between the time needed to cover three angles exactly and the corresponding maximum velocity and acceleration as a function of k. It is important to note that for a variable source rotation velocity s(t) to work, the maximum velocity and acceleration do not have to be large. Whether or not there is control of the source rotation velocity of s(t) is vital.
If the angle accuracy is within 0.05 • , the time required as shown in Table 1 is (433 + 5/8) (sec) for the best constant source rotation velocity s(t). By using a variable source rotation velocity s(t) as defined above, say k = 5, it takes (25+5/8) (sec) with the maximum velocity 0.42 (rad/sec) and the maximum acceleration 0.0064 (rad/sec 2 ). This kind of velocity and acceleration imposes no problem to hardware at all for most commercially available servo systems. For instance, suppose the X-ray source and sensor weigh 100 kg and the gantry diameter is 1m, an acceleration of 0.0064 (rad/sec 2 ) results in a torque 0.16 Nm. The Panasonic servo motor series A, E, and S and other commercially available products can easily meet this kind of requirements.
Clearly, a constant source rotation velocity s(t) is not a solution and this adds complexity to the control of s(t). Thus, it is interesting to note that the following modified Problem 2 allows an constant but not prefixed velocity solution.
where θ i is independent of l. The interpretation is that at any level i, projections have to be taken at p evenly spaced angles a i,l = a l + θ i = ((l − 1)/ p)2π + θ i , l = 1, . . . , p, with the same offset angle θ i . However, the offset angle θ i at level i may not be the same as the offset angle θ j at level j.
(b) Among all the solutions s(t)'s, find one that requires the minimum time.
The difference between Problems 1 and 2 is that Problem 1 requires that at any level v(t i ), samples are taken exactly at p evenly spaced angles a l 's while Problem 2 does not require that p evenly spaced angles at level i are the same as that at level j, as illustrated in Figure 3. This makes sense for practical applications where an accurate image of level v(t i ) Erwei Bai et al.

5
is reconstructed as long as p is large enough independent of whether the same angles are used for level j = i. Surprisingly, the seeming complicated Problem 2 is solved by a constant source rotation velocity s(t).

Theorem 2. (1) Let m ≥ 0 and k > 0 be any integers and consider
where P c = (p/(1 + mp))Pk is the period of s(t). Then, (2) The X-ray source s(t) defined above with k = 1 solves Problem 2. In other words, s(t) covers the angles a i,l 's exactly at precise volume levels v(t i ) with the minimum time.
because 2πm(l − 1) is an integer multiple of 2π. Thus, s(t) covers all the angles a i,l at all the levels i. Moreover, with k = 1, it covers all the angles at all the levels with the minimum time t q + (p − 1)P. This completes the proof.

Remark 1.
(1) The reasons why the X-ray source rotation s(t) of (15) can cover all the angles for all the levels are the following. (a) It is constant but selected by the user depending on the period of the heart volume v(t) and not refixed. (b) Different offset angles θ i 's are allowed for different levels. This is key difference between s(t) of (15) and methods currently available in the literature.
(2) Clearly, the minimum time is achieved with the minimum k. On the other hand, the small values of k imply high velocities of s(t). There is a tradeoff between the minimum time and the maximum velocity. The balance can be achieved by adjusting k. Note that the velocity of s(t) is 2π/(p/(1 + mp))Pk and moreover the acceleration of s(t) = (2π/P c )t is always zero for any m, k.
(3) If the imaging device has multiple sources/detectors, say n sources/detectors, the minimum time could be reduced by a factor of n.
(4) The optimal s(t) does not have to be fast. For instance, let m = 0 and p = 360, that is, there are 360 evenly spaced angles. Then, The CT rotates 360 k times slower than the motion of the object. In general, P c = p kP if m = 0. For instance, a heart rate is 60 cycles per minute and the CT could rotate as slow as 1 cycle in every 6 minutes when k = 1 and even slower when k > 1.
(5) Though s(t) does not have to be fast, it has to be accurate. For instance, let k = 1 and write a constant source rotation velocity s(t) as, for some d, At k = 1 and t = t i + (l − 1)P, d determines the increment. If s(t) is not exactly a constant where Δd(t) indicates variation of the velocity in s(t), then, at time t = t i + (l − 1)P, Obviously, the last term controls the accuracy. (6) The maximum difference between θ i and θ j for any k and m is Now, consider the example (8) again. Suppose three heart volume leaves v(t 1 ) = v(0), v(t 2 ) = v(1/2), and v(t 3 ) = v(5/8) are interested and at each level, three angles a l = (l − 1)2π/3, l = 1, 2, 3, are required. We can simply set s(t) = (2π/P c )t. Since P = 1 and p = 3, with m = 0 and k = 1, we have s(t) = (2π/3)t and that result in This s(t) covers three angles exactly at precise levels as shown in Figure 4. Also, note that the heart volume period is 1 (sec) and the X-ray source rotation period is 3 (sec), three times slower than the heart rate but still able to cover everything exactly with the minimum time. Theorems 1 and 2 show that artifacts can be completely eliminated if the motion is periodic and the source rotation velocity of s(t) is controllable. For a periodic motion, determination of the period is trivial. It can be estimated and predicted based on the analysis of the first a few cycles of ECG signals.

Quasiperiodic v(t)
To deal with a nonperiodic v(t), let us recall the properties of a periodic function. A function v(t) is periodic with period P if and only if it is completely determined by the first period In CT imaging, we assume that the dominant motion is cardiac and respiratory movement. Then, we observe that such motion can be nonperiodic, that is, the time that takes to finish one cycle may be different from the time to finish another cycle. Independent of time, the cycle however always starts at the minimum volume level extending to the maximum volume level, and then ends at the minimum volume level. Such a motion can be characterized by a class of functions that we call quasiperiodic. Let be a monotone increasing sequence. The interval P l+1 − P l is the lth "period" of v(t), that is, the time v(t) finishes a cycle from the minimum volume to the maximum volume and then back to the minimum volume level. In general, v(t) may be completely determined by its first "period" v(t) = v 1 (t), t ∈ [0, P 1 ). More precisely, let Erwei Bai et al. Then, for t ∈ [(l − 1)P, lP),
The intuition is that v(t) is periodic. However, in each "period" [P l−1 , P l ), v(t) is an expanded or compressed version of the first "period" v 1 (t), t ∈ [0, P 1 ). We give an example as shown in Figure 5: Because v(t) is no longer periodic, the definition of the volume level has to be modified. For the first period v(t) = v 1 (t), t ∈ [0, P 1 ), the levels are well defined: For the lth "period," we note P l−1 ≤ α l t i + P l−1 < P l and at t = α l t i + P l−1 , Thus, the ith volume level can be defined as (34) or in time domain where a i,l was defined before.
(b) Among all the solutions s(t)'s, find one that requires the minimum time.
Therefore, s(t) of (37) covers all the angles a i,l at all the levels v(α l t i + P l−1 ). The proof of the minimum time is similar as before.

Remark 2.
(1) Unlike the periodic case, any constant s(t) fails to cover all angles a i,l 's at all levels with the minimum time, even with slack variables θ i 's added.

8
International Journal of Biomedical Imaging (2) The idea of s(t) in (37) is again to interpolate a i,l at t = α l t i + P l−1 . Thus, the solution of s(t) is not unique and any other form of s(t) that interpolates a i,l at t = α l t i + P l−1 is also a solution.
(3) The source rotation velocity of s(t) in t ∈ [P l , P l−1 ) is In general if the "frequency" of v(t) is high, that is, the "period" P l+1 − P l is small. This implies that the velocity s (t) is small, provided that |β 1l /α l | is not large. In the case that the maximum velocity |s (t)| is of concern, interpolation can be done for both s(α l t i + P l−1 ) and the constraints on |s (α l t i + P l−1 )|. There exists a large volume of works in the literature in this direction, referred to as Hermite-or Fejertype interpolations.
(4) s(t) of Problem 3 may have discontinuities at P l 's. If these discontinuities are of concern, spline interpolations can be used. The result is a smooth piecewise polynomial at a price of increased computational complexity.
Control Problem 3 is solved in Theorem 3 for quasiperiodic v(t)'s provided the period P l − P l−1 is available. Estimation of the period P l − P l−1 is rather straightforward for quasiperiodic functions. v(t) is quasiperiodic, so the form of v(t) can be determined from its first a few cycles. Then, in theory, the period can be calculated by observing some points of v(t). To be precise, take v(t) of (31) as an example. Suppose the form of v(t) has been determined from its first cycle and we want to find the period P 2 − P 1 . Note v(t) = sin(c(t − 2π)), t ∈ (P 1 , P 2 ) for some unknown constant c. If c is determined, the period P 2 − P 1 = 2π/c. To estimate c, what we need is one observation v(t),t ∈ (P 1 , P 2 ), which leads to In implementation, because of noises, it is preferred to have several observations and the estimate c takes the form of the average of individual estimates.
To illustrate the idea, we consider an example: , a 2 , 120 • + (1080 • ) Figure 6: Controlled X-ray source for a quasiperiodic function. This is no longer periodic but quasiperiodic. Suppose three volume levels v(t) = 0.64, 1, 1.5 and three angles a l = (l − 1) * 2 * π/3, l = 1, 2, 3, are interested. Figure 6 shows the computer simulation results based on the polynomial interpolation and the estimates of the unknown periods. Clearly, three angles at precise three volume levels are imaged.
The advantage of the approach is that no other signal is needed. If other signals, for example, ECG signals, are available, the period can also be calculated from the gating of ECG signals.

SOURCE ROTATION CONTROL FOR NONEXACT MATCH
In practice, the minimum time to cover all the angles at all the levels can be important. In the previous discussions, the exact level v(t i ) and exact angle a l have to be matched to take images. This is referred to as the exact match and in this situation, the minimum time needed for any solution s(t) to cover all the angles at all the levels is α p t q +P p−1 . From a practical point of view, the exact match is not necessary as long as errors are small enough. We discuss the nonexact match in both angle and level. By nonexact matches, we mean that, at each level i, images do not have to be taken exactly at v(α l t i + P l−1 ) but in a small neighborhood of v(α l t i + P l−1 ). Similarly, angle a l does not have to be exact but within a small neighborhood.
To this end, let, at level i, the error bound δ > 0 and the corresponding times t i,l ≤ t i ≤t i,l be given that satisfy v(t) , τ i,l = α l t i,l + P l−1 , α lti,l + P l−1 .

9
The neighborhood of v(α l t i +P l−1 ) is then defined as Δ i,l with the corresponding time interval τ i,l . Clearly, v α l t i + P l−1 ∈ Δ i,l , max and τ i,l is the time interval in which the images for any angle a l , l = 1, . . . , p, can be taken at the level i. By the same token, the volume level i in a nonexact match case can be defined as or in time domain In the nonexact match case, elements of L i and T i and intervals not points. For simplicity, we assume τ i,l and τ j,k do not overlap for i = j and/or l = k. Now, the problem can be modified as to find a source angle profile s(t) with the minimum time satisfying s(t) − a l ≤ , l = 1, . . . , p, t ∈ τ i,k , for some k, where > 0 is the given angle error bound. If the velocity and acceleration of s(t) are of concern, the constraints |s (t)| ≤ M 1 , |s (t)| ≤ M 2 can be added in optimization, where M 1 and M 2 reflect the constraints of the physically achievable velocity and acceleration of the source angle s(t).

CONCLUDING REMARKS
Our primary insight is that cardiac imaging can be optimally performed by integrating modern control and imaging techniques, that is, by adaptively controlling the imaging device according to the real-time analysis and prediction of an individualized respiratory and cardiac motion functions. The work intends to synchronize the X-ray source rotation and cardiac CT for unprecedented performance. Clinically, our work may induce significant advancement of the field of cardiology.
Although the studies have been focused on the cardiac motion, our methodology can be in principles applied to compensate for either the respiratory motion effect or the combined effect due to both cardiac and respiratory motions. Also, with the controlled data acquisition process the X-ray source positions will be distributed quite randomly along the scanning circle. Therefore, new CT image reconstruction methods should be developed for optimal image quality in this context. Some Feldkamp-type/Katsevich-type algorithms can be developed for that purpose, similarly to what reported in the literature [30][31][32][33][34]. Further efforts are needed to formulate such integrated control strategies and associated image reconstruction methods, and evaluate their performance systematically.
It should be underlined that there is no difficulty in implementation of the proposed control schemes in terms of hardware limitations. The hardware requirements can be easily met by commercially available products, for instance, the Panasonic servo motor and driver series. The fact is that the control theory and techniques have been well developed over the past decades, but the problem or the opportunity is that control and imaging techniques have never been combined before for cardiac CT.
In conclusion, we have proposed an innovative idea and specific schemes to eliminate CT image artifacts due to cardiac and respiratory motions, and demonstrated the feasibility and importance of controlled cardiac CT. We feel that the idea reported here may open a new frontier of biomedical imaging. We will continue to work along this direction.