Keywords: Artis-Zeego, calibration, helical phantom, C-arm movement models, 3D reconstruction, motion compensation

C-arm CBCT is a widely used X-ray angiography imaging system. From these images, 3D reconstruction of the coronary arteries can provide more information in the interventional surgery of the cardiovascular diseases. For this dissertation, we have improved the 3D reconstruction quality of coronary artery by two ways. One way is modeling the data acquisition geometry. This calibration work will correct the non-ideal acquisition geometry. The geometry deviation may cause the artifact of the 3D reconstruction. The work has been applied on the advanced Artis-Zeego C-arm system. The other way is developing a 3D reconstruction method based on 2D motion compensation.

Chapter 2 describes the calibration work. We propose a new and complete data calibrating algorithm of Artis-Zeego C-arm system. First, we proposed an extended analytic algorithm based on the classical helical phantom to estimate the geometry parameters. The phantom geometries is easy to obtain in the experiment environment. All the geometries can be calculated from three intermediate vectors. In order to estimate the three intermediate vectors, we propose to calculate the projection of the axes paralleled to 𝑥 axis and more projection of the axes paralleled to 𝑦 axis. The experiment results indicate that the proposed algorithm is more accurate and increases the efficiency than the original algorithm.

Second, we convert the projection geometries to the C-arm coordinate system. We propose to estimate the nominal C-arm system by minimizing the oscillation results from the displacement of the phantom. The converted geometries can be independent of the placement of the calibration phantom. This procedure can simplify the alignment procedure of the iso-center before calibration. Last, Zeego system has more flexibility than the traditional C-arm. It provides more work positions, such as Head, Left, Right, Table15, Table30. We choose three representative work positions to realize the experiment. They are Head, Left, Table30 work positions. We proposed the movement models of the projection geometries by considering all the possible influencing factors objectively and systematically. The movement models can integrate the refinement and the prediction of the acquisition geometry. These movement models can simplify the calibration procedure in the clinical operations. The experiment results indicate that the proposed movement models have an acceptable precision to estimate the acquisition parameters. Head position has the least deviation and Left position has minor deviation. The two positions can be modeled as a single mean rigid motion. Table30 position has the most severe deviation and should be modeled with a residual translation part.

Chapter 3 describes the 3D reconstruction by 2D motion compensation. We propose a complete 2D motion compensation method. The method corrects the residual motion of each projection to increase the number of available projection images. The whole procedure includes 5 steps. They are maximum intensity forward projection, segmentation of the acquired projection, registration, initial and motion compensated 3D reconstruction. We adopt the simplified distance driven (𝑃 𝑆𝐷𝐷 ) to generate the maximum intensity forward projection. This projector is faster than distance driven and ray driven. Registration is the key part of the whole method. We adopt the mutual information (MI) with rigidity penalty to be the cost function. This cost criterion can avoid the extraction of the centerline of the vessels and do the registration accurately and robustly. To realize the optimization of the cost criterion,

Aujourd'hui encore, les pathologies cardiovasculaires représentent toujours la première cause de mortalité dans les pays développés (1/5 en Europe, 1/6 aux Etats-Unis). Parmi elles, nous pouvons citer les plus importantes, les coronaropathies, affectant les artères coronaires. L'une des principales causes est l'athérosclérose coronarienne. L'accumulation de plaques d'athéromes, c'est-à-dire des dépôts anormaux principalement de lipides et de calcium sur les parois des artères coronaires, entraîne la réduction de la lumière vasculaire. L'évolution de ces plaques athéromateuses peuvent engendrer des complications, telles que l'angor (ou angine de poitrine), l'infarctus du myocarde, l'insuffisance cardiaque, voire une mort subite.

Afin d'établir un diagnostic, l'imagerie médicale a une place déterminante notamment pour préciser la nature, la localisation et le degré de gravité des lésions coronariennes ou/et myocardiques. On peut citer les modalités d'imagerie, telles que l'imagerie par ultrasons (IVUS), par résonance magnétique (IRM, URMF), par rayons X (angiographie RX et tomodensitomètre X multidétecteurs), ou l'imagerie nucléaire (PET, SPECT).

Cependant, l'angiographie coronarienne RX (ou coronarographie) reste la modalité de référence permettant de déterminer avec précision le degré et le nombre de sténoses coronariennes ainsi que le nombre de vaisseaux atteints, malgré les importants progrès techniques réalisés ces dernières années des autres modalités d'imagerie. En effet, elle dispose de résolutions spatiale et temporelle supérieures aux autres systèmes d'acquisition, d'une visualisation temps réel du remplissage des artères grâce à la dynamique d'acquisition et d'une vision de l'ensemble de l'arbre des coronaires.

S'effectuant par cathétérisme (insertion d'un cathéter), la coronarographie génère des images RX de projection 2D des artères coronaires (gauche ou droite) enrichies de produit de contraste radio-opaque (riche en iode) selon plusieurs incidences (i.e. différents angles de projection).

Les types de traitement envisagés sont soient médicamenteux (bêtabloquants, inhibiteurs calciques, etc.), soient minimalement invasifs (angioplastie transluminale, avec ou sans pose d'un stent) relevant de la cardiologie interventionnelle et/ou, soit un acte de chirurgie cardiaque lourde (pontage coronarien). Seuls ces deux derniers traitements permettent la restauration du flux sanguin dans les coronaires.

C'est le traitement par angioplastie coronaire transluminale percutanée (PTCA) qui est le moins invasif, le moins couteux, le plus rapide, avec la récupération du patient la plus précoce et qui est donc le plus souvent préconisé. L'angioplastie consiste à insérer par voie endovasculaire un catheter muni d'un ballonnet à son extrémité jusqu'à la sténose (la zone à traiter) sous contrôle radiologique (selon une incidence donnée dans un plan 2D). Le gonflage du ballonnet va alors écraser la plaque d'athérome contre la paroi vasculaire. Le gonflage et dégonflage de ce ballonnet peuvent être répétés plusieurs fois jusqu'au retrait du catheter. Dans la très grande majorité des cas, l'angioplastie est complétée par la mise en place et le déploiement d'une endo-prothèse métallique ou stent (maintenant la lumière artérielle à son diamètre normal), ce qui améliorera le résultat initial et diminue le risque de re-sténose (reformation de la sténose). Plus d'un million d'angioplasties coronaires percutanées sont réalisées dans le monde chaque année pour traiter notamment l'angine de poitrine et l'infarctus du myocarde. Un patient à forte suspicion de lésions coronaires est souvent directement orienté vers une salle de cardiologie interventionnelle pour y passer une coronarographie, puis dans la continuité immédiate une PTCA. Le cardiologue interventionnel ne dispose donc pas de données préopératoires pour réaliser cette angioplastie (permettant la revascularisation coronaire). La représentation 3D des artères coronaires est, de ce fait, construite à partir des observations issues seulement de la coronarographie (acquises selon une incidence donnée dans un plan 2D).

Néanmoins, bien que des progrès aient été réalisés à la fois sur les instruments (guides, ballons coaxiaux et prothèses) et sur la technique (antérogrades, rétrogrades), le résultat de l'intervention reste très dépendant de l'expérience de l'opérateur ainsi que du tableau clinique. En effet, la nature 2D des observations angiographiques RX, la présence de calcification au niveau de la sténose, les conformations géométriques parfois irrégulières des vaisseaux (fortes tortuosités, longueur et localisation de la lésion) et le mouvement des structures (lié au battement cardiaque) sont à l'origine des difficultés rencontrées dans la réalisation du geste interventionnel. A cause de toutes ces difficultés, la précision de ce geste thérapeutique minimalement invasif est de ce fait impactée avec la réapparition à moyen terme de la lésion au niveau de l'endoprothèse.

Systèmes d'angiographie RX

Systèmes d'angiographie coronarienne RX conventionnelle ou monoplane

Le système d'angiographie mono-plan ou standard est composé de:

un arceau en forme de C (C-arm) formé d'un tube RX (émetteur) et d'un détecteur numérique plan d'image. Cet arceau est conçu de manière à couvrir une large gamme d'incidences autour du patient. un socle en L sur lequel repose le C-arm, fixé au sol, mais, pouvant tourner selon un axe vertical. une table motorisée de radiologie (sur laquelle repose le patient) pouvant se translater sur 3 degrés de liberté (voire 4).

Ce type d'angiographe permet l'acquisition de projections 2D des artères selon des incidences choisies par l'opérateur. Cette série d'images acquises suivant plusieurs angles de vue permettra au cardiologue de reconstruire mentalement le réseau coronarien. Néanmoins, le système mono-plan, se caractérise par une mobilité très réduite du couple émetteur/détecteur RX avec un seul mode statique d'acquisition.

Systèmes d'angiographie RX Bi-plan

Le système biplan, se caractérise par deux couples émetteur/détecteur RX, permettant une double acquisition RX simultanée. Il permet d'acquérir simultanément deux séquences d'images d'angiographie sous des incidences orthogonales (frontale et latérale). Ces deux séries d'images vont permettre au cardiologue interventionnel une représentation mentale 3D de la structure coronaire.

Système rotationnel Siemens Artis Zeego DynaCT

La radiologie interventionnelle associe l'imagerie radiologique avec une intervention chirurgicale. Assisté par un robot, le méolecim va déplacer la chaîne image autour du patient our visualizer l'organe d'interêt des image de projection RX.

Parmi les systèmes d'angiographie interventionnelle de dernière génération, le C-arm rotationnel Siemens Artis Zeego robotisé (voir Fig. 1) offre une très bonne souplesse de positionnement de par sa conception multi-axe (avec ses 6 degrés de liberté) et est particulièrement adapté aux salles de cardiologie interventionnelle et aux procédures exigeantes en termes de couverture et d'imagerie 3D avancée. Ce système fait partie intégrante de la plate-forme TheraImage du CHU de Rennes à Pontchaillou.

Limites et difficultés liées à l'angiographie coronarienne RX

Bien que l'angiographe RX (et quel que soit son type) permette de réaliser aussi bien le diagnostic que l'angioplastie, il fournit des observations angiographiques RX acquises selon une incidence donnée dans un plan 2D. Une interprétation mentale de ces observations en temps réel s'avère délicate mais nécessaire pour se représenter ces structures coronariennes 3D complexes. Une représentation imprécise peut entrainer des difficultés lors des phases diagnostique (localisation et caractérisation de la sténose), de planning (choix de ballon et stent, etc.) et d'intervention (montée du cathéter, déploiement du stent, etc.).

De par le principe d'acquisition RX par projection, les structures coronaires peuvent se superposées entre elles ou à d'autres structures environnantes (vertèbres, cage thoracique, cathéter d'injection, poumons, etc.) et alors confondus (effet de superposition) ou déformés (effet de raccourcissement des longueurs et variation du facteur d'agrandissement). L'aspect morphologique des artères coronaires peut donc être faussé.

Une autre limitation concerne le caractère invasif de ces examens. Le nombre d'acquisitions d'images 2D+T sous différentes incidences, et donc le nombre d'injections de produit de contraste afin de visualiser les artères coronaires, peut alors se révéler élevés pour un patient si son anatomie coronaire s'avère complexe. La dose d'irradiation du patient augmentera en fonction du nombre d'incidences envisagées et dans une autre mesure du personnel hospitalier réalisant cet examen plusieurs fois par jour et durant de nombreuses années. Le risque de complication lié à ce produit n'est pas à négliger (réaction allergique ou encore risque d'insuffisance rénale).

L'objectif de ce travail de thèse de doctorat concerne l'amélioration de la reconstruction 3D des coronaires afin d'améliorer le diagnostic, ainsi que la sécurité et la précision des interventions minimalement invasives. Dans un premier temps, une contribution majeure vise à améliorer l'étape de calibration du système d'imagerie rotationnelle R-X (Fig. 2).

Dans un second temps, une contribution majeure vise à proposer une nouvelle méthode de reconstruction des coronaires par compensation de mouvement (Fig. 3).

Chapitre 2 Calibrage du système d'angiographie rotationnel

Le chapitre 2 décrit l'étape de calibrage du système d'angiographie rotationnel. En effet, puisque la rotation du C-arm ne suit pas en routine clinique une courbe circulaire idéale autour du patient, des artéfacts significatifs, une détérioration de la résolution spatiale de la reconstruction et des informations anatomiques erronées surviennent. La géométrie d'acquisition est généralement représentée par la matrice de projection, qui peut être divisée en paramètres intrinsèques et extrinsèques. Ce sont les paramètres intrinsèques qui sont affectés par les déformations mécaniques du C-arm durant sa rotation. A cause de l'effet de couplage entre ces deux matrices intrinsèques et extrinsèques, il est difficile de les calibrer indépendamment de manière précise. Nous avons donc proposé un nouvel algorithme de calibrage basé objet inspiré de la méthode de Xu et al [23]. dédié au système robotisé Artis Zeego (Fig. 4.).

Premièrement, cet algorithme est basé sur l'estimation des paramètres géométriques du système à partir d'acquisitions d'un fantôme tubulaire classique de calibration (contenant des marqueurs radio-opaques, ici des billes d'acier de deux diamètres différents, réparties selon un ordre particulier et disposées sur une trajectoire hélicoïdale).

Les résultats expérimentaux montrent que l'algorithme proposé est plus précis et plus rapide que l'algorithme original proposé par Xu et al [23].

Deuxièmement, nous transposons les géométries de projection dans le système de coordonnées du C-arm. Puisque le système d'angiographie Siemens Artis Zeego est plus flexible que les systèmes C-arm traditionnels, il propose en effet davantage de positions d'acquisition, telles que Head, Left, Right, Table15, Table30. Parmi celles-là, nous en choisissons 3 représentatives pour évaluer notre méthode de calibrage. Nous avons proposé des modèles de mouvement des géométries de projection en considérant objectivement et systématiquement tous les facteurs possibles. Ces modèles de mouvement permettent de simplifier la procédure de calibration en routine clinique. The résultats expérimentaux indiquent que les modèles de mouvement proposés ont une précision acceptable afin estimer les paramètres d'acquisition. La position Head a le moins de déviation et la position Left a une déviation faible. Ces deux positions peuvent être modélisées par une simple transformation rigide. La position Table30 a la déviation la plus importante et peut être modélisé par une translation résiduelle. Nous avons évalué notre algorithme de calibrage seulement sur le fantôme physique tubulaire hélicoïdal.

Chapitre 3

Deux approches de reconstruction à partir d'angiographie coronaire rotationnelle sont généralement considérées : (i) La première envisage la reconstruction avec les projections qui sont associées à une même phase du cycle cardiaque (ou à une même position du coeur dans l'espace); le nombre de projections est de l'ordre de 4 à 6 (ce nombre peut augmenter avec le nombre de phases considérées en repos c.a.d jugées sans mouvement). Ce sous échantillonnage angulaire conduit à résoudre un problème inverse mal posé difficile à résoudre du fait de la complexité des images.

(ii) La deuxième catégorie de méthodes passe par une estimation préalable du mouvement des artères coronaires sur l'ensemble du cycle cardiaque, de manière à l'incorporer dans la procédure de reconstruction tomographique et exploiter toutes les projections (i.e. 80 -150 projections) dans la reconstruction de l'arbre coronaire, à un instant donné. Cela revient à réaliser une reconstruction compensée en mouvement.

Le chapitre 3 décrit la reconstruction 3D proposée par compensation de mouvement 2D. Elle corrige le déplacement résiduel à chaque projection afin d'augmenter le nombre d'images de projection disponible. Les étapes de reconstruction sont les suivantes (Fig. 5.): l'opérateur de projection, la segmentation des projections acquises, le recalage iconique, la reconstruction initiale et la reconstruction avec compensation de mouvement.

Nous avons opté pour l'opérateur de projection Distance Driven simplifié proposé par Oukili [29]. Cet opérateur de projection est plus rapide que celui de distance driven et de ray driven. Le recalage d'image est l'étape primordiale de cette méthode de reconstruction avec compensation de mouvement. Nous avons adopté comme mesure de similarité l'information Mutuelle (MI) avec un terme de rigidité. Cette fonction de coût permet d'éviter une étape de segmentation vasculaire et d'extraction des lignes centrales afin de réaliser un recalage précis et robuste. Pour optimiser cette fonction de coût, nous avons utilisé la méthode Advanced Adaptive Stochastic Gradient Descent (ASGD) [43]. Son coût calculatoire est plus faible et sa robustesse est meilleure qu'une descente de gradient classique (GD). La reconstruction synchronisée sur l'ECG initiale [27] utilise seulement les images de projection correspondant à la même phase cardiaque. Tandis que la reconstruction avec compensation de mouvement utilise des images projection des différentes phases cardiaques recalées. Nous avons adopté une reconstruction statistique itérative basée on MAP et l'a priori 𝐿 0 [7]. La méthode proposée a été évaluée sur des données simulées de référence CAVAREV [46] qualitativement puis quantitativement à travers différentes métriques, ainsi que sur des données réelles cliniques issues d'un angiographe rotationnel GE du CHU Pontchaillou de Rennes. Les résultats expérimentaux indiquent que la méthode proposée améliore visuellement la qualité de la reconstruction. Le contraste et les détails des reconstructions sont améliorés par compensation de mouvement. En augmentant la taille de la fenêtre temporelle associées aux phases cardiaques prises en compte dans la méthode de reconstruction avec compensation de mouvement, les artefacts d'image. L'amélioration de la qualité des images augmentent la visibilité des structures vasculaires ce qui permet de l'envisager en routine clinique pour le diagostic ou les actes minimalement invasifs des artères coronaires.

Chapitre 4

Le dernier chapitre reprend les principales contributions de cette thèse de doctorat et propose des perspectives de recherche concernant des futurs modèles de mouvement du système Artis Zeego. The heart is a muscle, about the same size as an adult human fist. Blood is pumped from the heart to the lungs, where it collects oxygen. This oxygen-rich blood is then pumped back to the heart and then to organs throughout the body through arteries. The blood then returns to the heat through the veins and is pumped to the lungs again. This is called circulation [1].

Coronary arteries are the heart's network of blood vessels. They exist on the surface of the heart, and they supply the heart muscle with oxygen. If the coronary arteries are narrow, the supply of oxygen-rich blood to the heart may become too low, especially during physical activity [1].

At first, this reduction in blood flow may not produce any symptoms, but as fatty deposits, or plaques, build up in the coronary arteries, signs and symptoms may emerge.

Coronary Heart Disease (CHD) is believed to start with injury or damage to the inner layer of a coronary artery. This damage causes fatty plaque deposits to build up at the site of the injury. These deposits consist of cholesterol and other cellular waste products. The accumulation is called atherosclerosis. See Fig. 1.1, this is the example of coronary atherosclerosis. It means the hardening of an artery specifically due to an atheromatous plaque. Atherosclerosis is asymptomatic for decades because the arteries enlarge at all plaque locations, thus there is no effect on blood flow. Even most plaque ruptures do not produce symptoms until enough narrowing or closure of an artery, due to clots, occurs. Signs and symptoms only occur after severe narrowing or closure impedes blood flow to different organs enough to induce symptoms. Most of the time, patients realize that they have the disease only when they experience other cardiovascular disorders such as stroke or heart attack. These symptoms, however, still vary depending on which artery or organ is affected [2].

Marked narrowing in the coronary arteries, which are responsible for bringing oxygenated blood to the heart, can produce symptoms such as the chest pain of angina and shortness of breath, sweating, nausea, dizziness of light-headedness, breathlessness or palpitations. Abnormal heart rhythms called arrhythmias are another consequence of ischemia.

Coronary heart disease is today the leading cause of death worldwide. It represents the first one among people over sixty and the second one among the youngest, with a number of deaths that yearly amounts to 3.8 million for men and 3.4 million for women. Coronary heart disease results from atherosclerosis, a vessel wall disease that leads over time to the accumulation of atheromatous plaques within the arterial walls. They involve morphological modifications of the lumen with the formation of stenosis that often results in thrombotic occlusion and acute ischemic syndromes (angina pectoris, myocardial infarction). The first symptom of atherosclerotic cardiovascular disease is a heart attack or sudden cardiac death (death within one hour of onset of the symptom). The symptoms of coronary heart disease include angina, shortness of breath (dyspnea) and heart attack [2].

CHD cannot be cured, but with today's technology, it can be managed effectively. Treatment involves lifestyle changes, and possibly some medical procedures and medications. Lifestyle recommendations include quitting smoking, eating a healthy diet, and exercising regularly. Surgery can open or replace blocked arteries, if the blood vessels have become very narrow, or if symptoms are not responding to medications. Laser surgery involves making several tiny holes in the heart muscle, which encourages the formation of new blood vessels. Coronary bypass surgery uses a blood vessel from another part of the body to create a graft that can bypass the blocked artery. The graft may come from the leg or an inner chest-wall artery. Angioplasty and stent placement use a catheter inserted into the narrowed part of the artery. A deflated balloon is passed through the catheter to the affected area. When the balloon is inflated, it compresses the fatty deposits against the artery walls. A stent, or mesh tube, maybe left in the artery to help keep it open [1].

Imaging Examination Technology

Imaging examination is a very important diagnostic technique, which plays a quite significant role in the early screening, therapy, assessment before and after operation of the coronary heart disease [3]. In the following section, we will introduce some important imaging examination technologies of CHD.

a. Non-invasive examination a.1 X-ray plain radiography X-ray plain radiography can show the pulmonary circulation of the patients with suspicious CHD. It also can examine that if there are left ventricle enlargement and myocardial infarction. X-ray plain radiography can reveal some complication of myocardial infarction, such as ventricular aneurysm. X-ray plain radiography has the advantages of simple operation and low cost. However, it can only be the auxiliary means of the clinical examination of CHD [3].

a.2 Cardiac CT imaging

The cardiac CT includes Electron Beam Computed Tomography (EBCT) and Multi-Slice Computed Tomography (MSCT). EBCT has the advantages of high dynamic resolution and high scan speed. However, EBCT has low spatial resolution, poor overall image quality. The equipment of EBCT is also too expensive. EBCT is gradually replaced by MSCT [4]. MSCT can clearly display the vessel lumen and wall of the coronary artery. Meanwhile, MSCT can be used to calculate the size of the calcified plaque. MSCT has high sensitivity and specificity in the evaluation of the coronary stenosis and can greatly reduce the motion artifacts and false positive. However, there are also some limitations in the clinical application. For example, the heart movement and breathing will lead to motion artifacts. For the patients whose heart rate is too fast, the examination can't be finished in one-breath-hold. Some blocker drugs should be injected to decrease the heart rate [5]. Currently, more advanced spiral CT has been developing. For example, 256-slices CT has the scanning time less than 1 cardiac cycle, which reduce the motion artifacts greatly. Now, MSCT has been recognized as the best non-invasive imaging examination technology of CHD [6].

In The signals will be used to detect the size of the vascular cavity, the structure of vessel wall, the character of atherosclerotic plaque. IVUS has obvious advantages in judging coronary atherosclerotic plaque stability, assessing coronary stenosis. IVUS is recognized as the gold standard for evaluating the vulnerability of coronary atherosclerotic plaque [7]. However, IVUS is invasive and not suitable to early screening. In addition, IVUS is limited by the diameter of the catheter. Therefore, there are many difficulties in coronary artery with small diameter and severe stenosis [8]. Coronary artery angiography (CAG) is the most widely used diagnosis technology of CHD. It can visualize the coronary arterial tree fully and reveal the location and abnormal degree of the lesions [9]. Currently, coronary angiography is the gold standard for the diagnosis of coronary atherosclerosis and the degree of vascular stenosis. CAG is indispensable for identifying the disease, determining the prognosis and developing the next treatment plan. CAG also can provide the basis and guidance for other coronary artery examination and interventional therapy. However, CAG can only show the shape of the lumen of the coronary artery and has limited information for the internal structure of the vessel wall and plaque characterization [9][10]. Currently, Quantitative coronary angiography (QCA) can make up the weakness of the visual evaluation. The 3D reconstruction technology of the coronary artery is helpful to improve the diagnosis accuracy. Coronary angiography is the topic of our dissertation.

Percutaneous Coronary Intervention

Percutaneous Coronary Intervention (PCI) is the therapeutic method that dredges the narrow and occlusive lumen by cardiac catheterization to improve the myocardial perfusion. To perform the interventional therapy, the physicians need to use a puncture needle with 1-2mm. The needle enters the vascular system by the superficial artery, then, the catheter is sent to the opening part of the coronary artery under the guidance of the coronary angiography imaging system. The vascular condition was revealed after the injection of the coronary agent by catheter. The physicians can use catheters to treat the lesions of the patients [11].

PCI mainly include Percutaneous Coronary Angioplasty (PTCA) and Coronary Artery Stenting (CASI) and so on [12]. See Fig. 1.2, for PTCA, the catheter is delivered to the ostium of the coronary artery that needs to be expanded by femoral artery or the radial artery. Then, the balloon with the corresponding size is transported along the guide wire to the narrow segment of the coronary artery. The balloon will be expanded with appropriate pressure and time according to the characteristics of the lesions to achieve the purpose of stenosis elimination. Recent years, PTCA has been constantly improved, which has been the basic method of interventional therapy of coronary heart disease.

Coronary Artery Stenting inserts the stent into the narrow segment of the coronary artery to support the vascular wall and to keep the blood flowing. The stents are usually a network with gaps and made of stainless steel or alloy material. CASI can reduce the elastic recoil after PTCA and close the coronary dissection. After stenting, the new endothelial cells can cover the appearance of the stent. The stent was eventually wrapped in the vascular wall. The vessels will keep the open state. CASI is a great progress in the cardiac interventional therapy. It can protect the safety of the patients with complex coronary artery anatomy and applicable to the patients with acute coronary artery occlusion [13]. PTCA and CASI have a significant effect in saving the lives of patients with acute myocardial infarction, improving the treating and curing rate, increasing the life quality of the patients with myocardial infarction. Coronary angiography imaging system plays a significant role in the aid and guidance of the interventional therapy.

Traditional Angiography System

For the traditional X-ray angiography system, an injection of contrast agent during acquisition permits to easily visualize the coronary tree, delineate its branching pattern, and outline the inner diameter of a coronary artery. Multiple images are acquired by changing the angles of the X-ray system to visualize the coronary tree from different perspectives [14]. However, this traditional 2D angiography produces a 2D silhouette that could misrepresent the true anatomic characteristics of the 3D vessel lumen.

There are two main limitations of this angiography [15]: vessel overlap and vessel foreshortening. Vessel overlap occurs with superimposition of two images and can be recognized by experienced angiographers. Foreshortening of a vessel segment could be unrecognized and result in missed lesions, underestimation of stenosis severity, and inaccurate lesion length measurements. Except the two limitations, vascular bending and non-optimal view are also huge challenges for the vascular pathological study of the traditional X-ray angiography system [16][17][18][19][20].

During coronary angiography, the selection of the optimal view is dependent on each patient's anatomy and the ability of the physician to identify the optimal imaging projection to better visualize a specific segment of the coronary artery [21][22].

See Fig. 1.3, to improve the limitation of traditional angiography, the biplane imaging system have been adopted to minimize artifacts through the acquisition of two orthogonal views of the vessel segment of interest [23][24]. However, traditional angiography has remained a subjective, operator-dependent imaging modality. To improve the diagnostic accuracy of angiography, multiple projections are necessary according to view angles. Angiographic techniques using automated rotation with 3D reconstruction have been invented [25][26]. Rotational angiography can provide more angiographic information than traditional angiography and the acquired projection images can be utilized for 3D volumetric reconstruction because of the acquisition technique. In general, images are acquired over a large angle using high-speed iso-centric rotation of a C-arm. Rotational angiography has several important imaging advantages over traditional angiography. The rotational angiography can provide a 3D image. It also provides a less operator-dependent method of acquiring images. It acquires markedly more information than traditional angiography that can be used for patient care decisions. Rotational angiography has been demonstrated to require less radiation and contrast than traditional angiography.

3D Rotational Angiography System

Over the last decade, the combination of well-defined rotational techniques and powerful computer technology has resulted in clinically useful 3D rotational angiography (3D-RA) systems [15]. 3D angiography improves the accuracy of traditional techniques. The contrast and radiation exposure will decrease. The interventions will be performed with precision guidance using this new imaging technique. In clinical, the physician with 3D X-ray angiography will be able to identify inaccurate 2D vessel lengths, quantify the tortuosity of a vessel, and identify the optimal angle to better estimate a vessel lesion [27][28]. 3D rotational angiography facilitates the performance of the procedure and provides critical information that modifies the treatment technique or eliminates unnecessary therapy. During interventional procedures, 3D angiography can be used to identify and assist the placement of an exactly sized coronary stent. 3D angiography has the potential to deeply impact patient care during the invasive and interventional procedures. This technique has the significant impact on the diagnosis and treatment of patients with all forms of vascular disease.

Another hand, 3D rotational angiography can be performed with acceptable patient radiation dose, comparable to cardiac CT. The modification to the current exposure modulation protocols can significantly reduce the ED (effective dose) of 3DRA in many patients while maintaining the high quality for clinical usefulness of imaging [29].

3D C-arm imaging system is a machine for the rotational angiography. It has been continuously improved over the years. It uses the 2D X-ray projection data acquired with detector to generate the 3D reconstruction. Although C-arm CT data acquisition is increasingly automated for ease of use, the following steps are usually involved [30]. See Fig. 1.4, First, the patient needs to be optimally positioned such that the region of interest is visible in all X-ray views acquired during a rotation around the patient; With the patient properly placed, the C-arm is initially driven into the end position, then a safety run is performed during which the C-arm is slowly moved into the actual C-arm CT start position. This safety run is required to rule out collision during the actual scan. After the C-arm has reached its start position, a short fluoroscopic X-ray pulse is applied to initialize the automatic exposure control. At this (a) Arits zee ceiling-mounted system (b) Artis zee floor-mounted system point, the system is ready to begin a 3D run during which the C-arm rotates from the start position to its end position. Then raw data acquisition is performed by activating a dead-man switch. During the acquisition, the contrast agent should be injected continuously.

3D reconstruction is based on the acquisition of rotational projection sequences with or without electrocardiogram (ECG) triggering. The reconstruction procedure is used to transform the projection into volume information. This volume data represents the X-ray absorption of the object of interest in each voxel position. The data consists of the CA-enhanced vascular structures as well as other objects with X-ray absorption capability surrounding the vessels. In combination with 2D fluoroscopic or radiographic imaging, the information provided by 3D C-arm imaging can be very valuable for therapy planning, guidance, and outcome assessment for complicated interventions [31].

Traditional 3D rotational C-arm system

A traditional 3D rotational C-arm CT system often comprises a stand and a C-arm gantry to which the flat panel detector (or image intensifier), X-ray tube, and collimator are attached. In theory, the C-arm keeps the X-ray tube, collimator, and detector exactly aligned under varying view angles. See Fig. 1.5, Artis zee family (Siemens AG, Healthcare Sector, Forchheim, Germany) are the state-of-the-art rotational C-arm systems. Fig. 1.5 (a) is the ceiling-mounted C-arm system. Fig. 1.5 (b) is the floor-mounted system. This kind of 3D C-arm CT system has two rotation axes and involves a mechanically fixed center of rotation commonly referred to as iso-center. The C-type gantry can be rotated and translated to increase patient coverage and access. In the following of this section, we introduce the main components of this C-arm system.

a. X-ray source

The X-ray tube, X-ray generator, and X-ray control system are crucial components of any C-arm imaging system. They determine tube voltage, tube current, and irradiation time, respectively. These exposure parameters are essential for X-ray imaging, since contrast-detail perceptibility and dose depend on them. b. Flat-panel detector Until the 1990s, C-arm systems for real-time angiography imaging used to rely on X-ray image intensifiers (XRIIs). This technology has a number of inherent disadvantages that limit its utility. For example, the convex input screen of XRIIs results in a non-homogeneous image quality across the output image. In addition, scatter processes of light and electrons within the image intensifier limit the contrast resolution. Another limitation is the patient access as well as the flexibility during angulation due to the large size of XRIIs. Nowadays, the advanced rotational C-arm system, for example, Artis-zee family is equipped with flat-panel detector (FD). Flat detector either avoids or at least reduces some of the major disadvantages of image intensifiers. The most important technical advantages of flat detectors are: homogeneous image quality across the entire image area resulting in distortion-free images; position-independent spatial resolution; tightly enclosed square or rectangular active imaging areas offering improved patient access [30].

c. C-arm gantry

The C-arm configuration owes its name to its C-shape. The detector and the image acquisition system are mounted on the top extreme of the C-arm gantry, and X-ray tube and collimator are mounted on the bottom extreme.

In clinical acquisition, we need to change the point of view to observe the vessel. The gantry system is able to provide two degrees of freedom. See Fig. 1.6, α is the rotation angle or primary angle (PA), when α > 0, we have a Right Anterior Oblique (RAO) view, and for α < 0, we obtain a Left Anterior Oblique (LAO) view. β is the angulation angle or secondary angle (SA), when β > 0, we have the Cranial (CRA) view, and when β < 0, we have the Caudal (CAU) view. When we set α = 0, β = 0, we call it as Anterior-Posterior (AP). There is a mechanically fixed center of rotation commonly referred to as the isocenter.

The distance between X-ray source and the detector is defined as focal length. The focal length can be adjusted, the X-ray source remains at a constant position while the detector moves up and down.

d. Scan protocol

The clinical C-arm CT scan protocols means a careful choice of imaging parameters such as system dose, scan time and pixel binning at the detector and so on. Pixel binning means that the outputs of neighboring detector pixels are combined into one reading.

Artis-Zeego C-arm

Artis-Zeego is one of the new Artis Zee family of interventional imaging systems. It is the first multi-axes system that can be positioned as the way you want. And it can be controlled with far greater ease and precision than a traditional floor or ceiling-mounted system. With its multi-axes movement, this system provides greater positional flexibility and broader coverage (more angulations, head-to-toe coverage, multiple work positions)-including large volume cross-sectional images up to 47 cm in diameter. The robotic technology integrated into the system makes it possible to position the C-arm exactly according to the view required, anywhere in a sphere around the patient. The movement of system can be coordinated with the operating table, which means that the physician is allowed to operate at an optimal position [32][33]. See Fig. 1.7, Artis-Zeego system is equipped with a flat panel detector, the X-ray source, the arc-support, the collimator, the robotic arm, the robotic shoulder and so on. In Fig. 1.7 (1)-( 12), we can see the movement of each axis. The C-arm keeps the detector exactly aligned with collimator and X-ray source, detector and collimator rotate in synchronization, X-ray source is built into the C-arm below the collimator [34].

The height of the table and the SDD (source to detector) can be adjusted, all the geometry can be read from DICOM acquisition data. See Fig. 1.8, A is DTO (table to object distance), the imaging region must be placed on the isocenter, the source, the isocenter and the middle of detector are collinear. G is SDD (source to detector distance), B is SAD (source to axis or isocenter distance), C is source to patient distance, D, E, F are the relative positions between isocenter and table system. Artis-Zeego offers more positioning flexibility than conventional systems. It provides more work positions [34]. In Fig. 1.9, we plot all the work positions: Head Side, Left Side, Right Side, Left Side with Table rotated by 30 °, Left Side with Table rotated by 15 °. We should choose one of them first when we do a 3D acquisition.

Its isocenter height can be adjusted to the operator, reducing the physician fatigue associated with long procedures. And Artis-Zeego can be parked compactly to give the clinical team easier access to the patient, making it an ideal solution for hybrid rooms, see Fig. 1.9(g) , by freeing up space in the angio lab and OR, Zeego enables your clinical teams to work faster and more efficiently [32][33].

This system provides an enhanced 3D imaging (syngo DynaCT) for greater precision. syngo DynaCT is a unique form of rotational angiography that uses special reconstruction algorithms to generate CT-like images in less than one minute directly. Artis zeego with syngo DynaCT enables large volume coverage that rivals a traditional CT scanner, ideal for obese patients; portrait volume imaging, ideal for liver, spine and carotids; with syngo DynaCT 360, it delivers large-organ, soft-tissue images in just 6 seconds. It also provides 3D images at the lowest possible dose levels, for neuro applications, it can achieve an effective dose of 0.3 mSv [START_REF] Fessler | A paraboloidal surrogates algorithm for convergent penalized-likelihood emission image reconstruction [C[END_REF].

Objective

For this dissertation, we try to improve the quality of the 3D reconstruction by two ways. First, we improve the calibration work that will reduce the 3D artifacts due to the non-ideal geometry. Second, we propose a new 3D reconstruction method that based on 2D motion compensation to improve the 3D imaging quality. We introduce the two research topics in the following.

Geometric Calibration

C-arm gantry rotation is subject to geometrical distortion, such as the motion of the X-ray source and detector differing significantly from a simple circular orbit. For example, the gravity-induced mechanical flex would cause such non-idealities. However, some CBCT reconstruction algorithms such as FDK filtered back-projection [36] assume a circular source-detector trajectory, it requires an accurate calibration method to account for geometric non-idealities (e.g. displacement of the x-ray source, rotation/tilt of the detector, changes in source-to-detector distance), such non-idealities in the source-detector orbit can result in a misregistration, a loss of detail, and image artifacts. So in order to reconstruct the volume data from projection images with high spatial resolution, the geometric calibration of the imaging system is required before the clinical data acquisition. The calibration should be performed periodically, approximately every 6 months or whenever hardware changed. In general, geometric calibration relates the 3D coordinates (𝑥, 𝑦, 𝑧) of voxels in the reconstructed image to the 2D coordinates (𝑢, 𝑣) of pixels in the projection image. Geometric calibration usually consists of two stages [37]: (i) characterization of pose through the range of source-detector orbit; (ii) correction of geometric non-idealities in the process of 3D reconstruction. Methods of pose characterization include the use of tracking systems and the image-based method, for the first one, a tracking system is used to monitor mechanical motion [38]; for the image-based calibration, we often operate on the projection acquired from a calibration phantom [39]. A classical calibration method is based on a helical phantom of BBs from which projection matrices are computed. A projection matrix is often used to describe the linear relationship between 3D voxel coordinates and 2D pixel coordinates, or a set of geometric parameters describing the imaging system such as source and detector positions, detector rotation angles and etc. Fig. 1.10. illustrates the angiography imaging system model of cone-beam imaging. Plane G is the focal plane, plane D is the image plane, 𝒔 is the X-ray source, distance F is the focal length. The optical axis is the line passing through 𝒔 perpendicular to G, the optical axis intersects D at point 𝑠 ′ , called the principal point. This is a perspective projection. We have three coordinate systems, the camera system (𝑋, 𝑌, 𝑍) built on the X-ray source, the world system(𝑋 𝑤 , 𝑌 𝑤 , 𝑍 𝑤 ) built arbitrarily and the detector pixel coordinate system(𝑢, 𝑣). Now, we transfer 3D world coordinates to detector pixel coordinates [40]. First, we perform a transformation from the 3D world coordinate system (𝑋 𝑤 , 𝑌 𝑤 , 𝑍 𝑤 ) to 3D camera coordinate system(𝑋, 𝑌, 𝑍):

[ 𝑋 𝑌 𝑍 ] = , 𝑅 𝑇 -[ 𝑋 𝑤 𝑌 𝑤 𝑍 𝑤 1 ] (1.1)
Where matrix R is a 3×3 rotation matrix and T is a translation vector. R and T constitute the extrinsic parameters. Second, a perspective projection is applied to the angiography context:

[ 𝑈 𝑉 𝑊 ] = 𝐼 [ 𝑋 𝑌 𝑍 ] (1.2) 𝑢 = 𝑈 𝑊 , 𝑣 = 𝑉 𝑊 when 𝑊 ≠ 0
Where matrix I is the matrix of the intrinsic parameters. An expression for intrinsic matrix I is:

𝐼 = [ 𝐹𝑘 0 𝑢 𝑠 0 𝐹𝑘 𝑣 𝑠 0 0 1 ] (1.3)
Where k is the inverse pixel size, F is the focal length and (𝑢 𝑠 , 𝑣 𝑠 ) are coordinates of the principal point on the image plane coordinates. Equations (1.1) and (1.2) can be combined as follows: The projection matrix P is defined as intrinsic and extrinsic matrices:

𝐸 = , 𝑅 𝑇 - 𝑃 = 𝐼𝐸 (1.5)

3D Reconstruction Method

3D reconstruction of coronary arteries from X-ray projections requires the acquisition of angiographic data along a circular short-scan trajectory (i.e., 180 °-220 °) during continuous CA injection. See Fig. 1.11,in theory, the vessel to reconstruct should be stationary and motionless during the acquisition procedure, though this assumption is not achievable in cardiac imaging. The breathing and patient movement can be eliminated by the setup and instruction. The cardiac motion of the heart must be resolved with complex strategies such as ECG-gating or motion compensation reconstruction techniques [31]. ECG-gated reconstruction uses only a subset of images which correspond to a particular cardiac phase and the 3D result of the coronary is for that specific time. Usually, for the end-systolic and late-diastolic cardiac phases, cardiac movement has the least motion. It is a good choice for the gating to provide a high reconstruction quality. A gating window is often used to do the gating reconstruction. Multiple images are taken from each heartbeat (i.e., 3-5 projection per gating window). This is an ill-posed problem. Iterative reconstruction methods can overcome this problem by incorporating a priori assumptions in the reconstruction process [41].

When the chosen window is too wide, motion artifacts will occur. Now, we can introduce motion compensated method. These methods try to estimate the cardiac motion and incorporate it into the tomographic reconstruction. There are 2D and 3D motion compensation approaches for 3D coronary reconstruction. 3D motion compensation uses the 3D centerline information to determine the motion-vector filed before reconstruction. 2D motion compensation tries to correct the residual motion of each projection images [42][43] and apply more projections in the reconstruction.

Proposed Method

Here, we describe the proposed methods of the two improvements briefly. For the calibration, we propose a new calibration algorithm of the data acquisition geometries of Zeego system. For the 3D reconstruction, we propose a new 3D reconstruction method by 2D motion compensation to improve the 3D quality.

For the calibration work, in order to estimate the acquisition geometry of Artis-Zeego C-arm system, we propose a new data calibrating method. First, we propose an extended analytic algorithm based on the classical helical phantom to estimate the acquisition geometry parameters. The algorithm adopt the invariant characteristics of the perspective projection to calculate three intermediate vectors, all the geometry parameters can be calculated from these three intermediate vectors. Then, we estimate the nominal C-arm system by minimizing the oscillation from the displacement of the phantom. All the geometries are transferred to this system that is independent of the phantom placement. Last, we consider separately the intrinsic and extrinsic parameters and propose different movement models by the estimated posterior information. Since this C-arm system is more complex and flexible than the traditional C-arm, we analyze three representative work positions: Head side, Left side, Left Side with Table rotated by 30 ° position (Table30 °briefly in the following) and propose different movement models.

For the 3D reconstruction work, we try to improve the reconstruction method by a new 2D motion compensation method. This method corrects the residual motion of each projection and makes more projection images with different cardiac phases available for the reconstruction. The whole procedure includes forward projection, preprocessing of the acquired projection, registration, 3D reconstruction. Each step will affect the final reconstruction. We try to improve each step in our algorithm. For the forward projection, we adopt the simplified distance driven projector 𝑃 𝑆𝐷𝐷 to generate a forward projection of the initial reconstruction. For the key problem of registration, we adopt the 2D registration algorithm based on the classical mutual information (MI) combined with a novel rigidity penalty. This cost criterion avoids the extraction of the centerline of the vessels and does the registration accurately and robustly. An advanced optimization method called adaptive stochastic gradient descent (ASGD) is applied to decrease the computation time of the registration part. For the 3D reconstruction algorithm, we combine the iterative reconstruction algorithm with the 2D motion compensation.

Outline

The main work of the proposed 3D reconstruction improvements are showed in Fig. 1.12. The main content of each chapter is described as follows:

Chapter 1: Introduction. This chapter introduces the main background of the dissertation. First, we introduce the clinical background of angiography CBCT system and the latest Artis-Zeego C-arm system. Second, we introduce the main problems to be solved, including calibration and 3D reconstruction. Third, we describe the main methods to resolve the problems. Last, we describe the outline of our dissertation.

Chapter 2: Calibration. See Fig. 1.12 (a), this chapter describes the work about calibration. For the first part, we introduce the background of calibration work. For the second part, we introduce some classical calibration algorithms and then we propose our extended calibration algorithm by a classical helical phantom. We show the evaluation results in this part. For the third part, we introduce some classical movement models of C-arm system, then we propose our movement models of this advanced Zeego system. For Zeego system, we transfer the estimated geometry to the nominal C-arm coordinates, then, our different predictive movement models are proposed on the above estimated geometry. We evaluate the proposed models in this part. For the last part, we discuss and conclude our work of this chapter.

Chapter 3: 3D Reconstruction by motion compensation. See Fig. 1.12 (b), this chapter describes the work about 3D reconstruction by motion compensation. First, we introduce the background of the work of motion compensated reconstruction. Second, we introduce some classical methods. Third, we propose a complete reconstruction method by 2D motion compensation. The main parts of this reconstruction procedure include forward projection, registration, 3D reconstruction, processing of the acquired projection. We describe each step in details, especially the registration part. Fourth, we evaluate this reconstruction algorithm. Last, we discuss and conclude our work of this chapter. Chapter 4: Summary and Prospect. We make a summary and propose some future prospects. First, we make a conclusion of the main work of Chapter 2 and Chapter 3. Second, we propose some plans of the future study to improve the calibration work and the reconstruction by 2D motion compensation.

Chapter 2 Geometric Calibration

Introduction

This chapter describes the 3D reference object-based calibration methods for determining the geometric alignment parameters for cone-beam scanners. Recently, cone-beam computed tomography (CBCT) with flat-panel detector has become a powerful technique in the interventional therapy. The imaging system rotates around the patient to obtain sequence of projection images. The projection images with the same cardiac phase are used to generate a 3D reconstruction of the coronary arteries.

When we do the 3D reconstruction, we often have the following model [1]:

𝐴𝒇 + 𝒃 = 𝒀 (2.1)
𝒇 is the 3D vector of the unknown image intensity, 𝒀 is the 2D vector of the projection image, 𝒃 means the acquisition noise, 𝐴 is the system matrix, the calculation of the system matrix will affect the reconstruction result directly [2]. Many approaches were developed to generate the system matrix, including voxel driven [3], separable footprint (SF-TT) [4], ray driven [5], and distance driven [6] and so on. The simplified distance driven [7] projector 𝑃 𝑠𝑑𝑑 was proved to be faster and comparable accurate than the classical ray driven (RD) and Separable Footprint (SF-TT), in this chapter, we use this projector in the following reconstruction. Fig. 2.1 represents the CBCT acquisition, 𝒇 is the 3D object which located between the X-ray source s and the flat-detector, the distance from source to detector (SDD) is called focal length. We define the C-arm coordinate system (𝑜 1 , 𝑥, 𝑦, 𝑧), the origin 𝑜 1 usually locates on the iso-center of the C-arm system, α is the rotation angle, x axis points to the initial position (α = 0) of s, z axis is along the rotation axis of the X-ray system, y axis is defined by right-handed criteria, the unit vector of each coordinate system is(𝑖, 𝑗, 𝑘), with‖𝑖‖ = ∆𝑥, ‖𝑗‖ = ∆𝑦,‖𝑘‖ = ∆𝑧, (∆𝑥, ∆𝑦, ∆𝑧) represent the voxel size. The distance from source s to origin 𝑜 1 is SAD. We define the detector coordinate system as (𝑜 2 , 𝑢, 𝑣). 𝑜 2 = (𝑢 𝑠 , 𝑣 𝑠 ) is the principal point of the system, 𝑢 is parallel to y axis, 𝑣 is antiparallel to z axis, (∆𝑢, ∆𝑣) represent the pixel size. The source system locates on the X-ray source s, 𝑥 𝑠 is anti-parallel to z axis, 𝑦 𝑠 is anti-parallel to y axis, 𝑧 𝑠 is anti-parallel to x axis.

There are three main steps for generating the system matrix by projector operator 𝑃 𝑠𝑑𝑑 based on simplified distance driven [7]. See Fig. 2.2, First, for each voxel n, we compute the conic projections of the eight cubic vertices. Second, we draw the smallest rectangles enclosing these projection points and we assume that the probability that a photon emitted by this voxel x will be detected in the rectangular surface (𝐴 𝑇 ) is 1. Last, we calculate the overlap surface between each pixel and the rectangular area (𝐴 𝑃 ), and the coefficients of the projection matrix is inferred as an area ratio such as: 𝑎,𝑢 𝑘 , 𝑣 𝑙 ; 𝒏-= 𝐴 𝑃 𝐴 𝑇 ⁄ . For the first step, we should use the acquisition geometry to calculate the cone beam perspective projection. We calculate the extrinsic matrix as follows: ]

𝑇 𝑐-𝑎𝑟𝑚 𝑠𝑜𝑢𝑟𝑐𝑒 = ,0,0, 𝑆𝐴𝐷-′ 𝑅 𝑐-𝑎𝑟𝑚 𝑠𝑜𝑢𝑟𝑐𝑒 = (𝑅 𝑠𝑜𝑢𝑟𝑐𝑒 𝑐-𝑎𝑟𝑚 ) ′ 𝐸 = ,𝑅 𝑐-𝑎𝑟𝑚 𝑠𝑜𝑢𝑟𝑐𝑒 , 𝑇 𝑐-𝑎𝑟𝑚 𝑠𝑜𝑢𝑟𝑐𝑒 - (2.2)
The superscript and subscript represent two coordinate systems respectively. The intrinsic matrix is:

𝐼 = [ 𝑆𝐷𝐷 0 0 0 𝑆𝐷𝐷 0 0 0 1 ] (2.
3)

The projection matrix is:

𝑃 = 𝐼𝐸 (2.4)
If we have one point 𝑀 = (𝑥, 𝑦, 𝑧, 1) 𝑇 on the 3D object, we can compute its projection point 𝑚 = (𝑢, 𝑣, 1) by 𝑃:

𝑠𝑚 = 𝑃𝑀 (2.5)
We deduce the formulation of 𝑢, 𝑣 by (2.2)-(2.5), we have (2.6):

𝑢 = 𝑆𝐷𝐷 • (-𝑧) (-𝑐𝑜𝑠(𝛼) • 𝑥 -𝑠𝑖𝑛(𝛼) • 𝑦 + 𝑆𝐴𝐷) ⁄ 𝑣 = 𝑆𝐷𝐷 • (𝑠𝑖𝑛(𝛼) • 𝑥 -𝑐𝑜𝑠(𝛼) • 𝑦) (-𝑐𝑜𝑠(𝛼) • 𝑥 -𝑠𝑖𝑛(𝛼) • 𝑦 + 𝑆𝐴𝐷) ⁄ (2.6
) 𝑢, 𝑣 are with the physical unit.

However, in practical applications, due to the errors associated with machine manufacturing system integration and gravity, the real path of the imaging system will deviate from the ideal circle. So these deviations will cause the inaccurate correspondence between 3D object and 2D projection, equation (2.6) will not be correct. This will result in some significant artifacts and the degeneration of spatial resolution of the reconstruction or the erroneous anatomical information.

A calibration procedure will take into account the geometric nonidealities and calculate the geometric parameters. These parameters should be considered in the reconstruction algorithm.

In computer vision, the acquisition geometry is generally represented by a projection matrix, which can be divided into intrinsic and extrinsic matrices [8]. The intrinsic parameters of C-arm can be affected by its mechanical deformation during the C-arm movement. Due to the coupling effect between intrinsic and extrinsic matrices, it is difficult to calibrate the two matrices independently with high accuracy. Multi-images [9] estimated the two matrices more accurate. However, the operation procedure by the algorithm was too tedious. Gorges [10] proposed another method to estimate the intrinsic parameters. However, the accuracy was not good enough. Another way to consider this problem is to calculate the defined geometry parameters of the cone-beam imaging instead of the whole projection matrix [11][12][13][14][15]. We follow this second way in our approach.

A dedicated phantom is often used in the calibration work. Strobel [8] and Rougée [16] used the classical projection matrix method on the classical phantom. The classical phantom carried beads regularly arranged along a helical trajectory. Zhang [17] proposed a calibration algorithm on a plate phantom. Cho [13] proposed an analytic algorithm based a cylinder phantom, which has 24 steel ball bearings in two plane parallel circles. Mennessier [15] proposed a direct analytic algorithm on a cubic phantom. Nine geometry parameters for each view can be determined by three intermediate parameters. Among these phantom, the classical helical phantom, plate phantom, cubic phantom are easy to obtain in the experiment environment. The algorithm based on the plate phantom can't estimate the intrinsic parameters at each angle separately. The cubic phantom has too few markers to prevent from the effect of acquisition noise. There are a lot of other methods with their designed phantoms [11,12,14].

As we know, it is not easy to place the phantom center to the iso-center of the acquisition system exactly. Ford [18] developed an algorithm that can estimate geometry related to iso-center and not require precise positioning of the phantom at the iso-center. This method is based on the algorithm of Cho [13]. However, this method will lose the accuracy when the title angles are large. We propose the algorithm to estimate the nominal C-arm coordinate system in our imaging system.

Another hand, some researchers built the movement models [9], [19][20][21][22] of the geometry parameters. These models improve the accuracy by the model refinement. The models can estimate and predict the geometry at each acquisition position and angle. These models may be based on some assumptions or the posteriori information of the motion by the calibration. Dumay [19] first proposed an ideal C-arm movement model, which was not accurate in the real condition. Gorges et al [9] evaluated the geometry by the multi-images algorithm. They analyzed the characteristics of the geometry parameters and take the estimated geometry as the posteriori information of their models. Li et al. [22] adopted some posteriori information proposed by Gorges et al [9] and proposed a new C-arm movement model. However, this method did not consider the residual deviation of the gantry and was difficult to calculate. There are many other movement models [20,21]. However, all of these are based on the traditional C-arm system. The simple assumptions of the above models are not suitable for Zeego system due to its complex movement.

The proposed data calibrating method is illustrated in Fig. 2.3. In part 2, we propose to estimate the acquisition geometries of a set of orientation and paths of the C-arm system. The proposed extended algorithm is more reliable and easier to perform than the original algorithm. In part 3, we transfer the geometry to the nominal C-arm system. Then, we proposed and calibrate the movement models of the projection geometries. In part 4, we have the evaluation result of the proposed method.

Geometric Calibration

Classical Method

We introduce some classical calibration algorithms for cone-beam imaging systems in this part. They are the classical projection matrix method, iterative parametric method, multi-images method and some advanced algorithms.

Classical projection matrix method

As described in equation (2.5), the projection matrix connects each 3D point 𝑀 𝑖 = ( 𝑋 𝑖 , 𝑌 𝑖 , 𝑍 𝑖 , 1) 𝑇 and it 2D projection 𝑚 𝑖 = ( 𝑢 𝑖 , 𝑣 𝑖 , 1) 𝑇 , 𝑖 = 1, … 𝑁. We have the projection matrix as:

𝑃 = [
𝑈 = [ 𝑀 1 𝑇 0 0 𝑀 1 𝑇 -𝑢 1 * 𝑀 1 𝑇 -𝑣 1 * 𝑀 1 𝑇 … 𝑀 𝑁 𝑇 … 0 … -𝑢 𝑁 * 𝑀 𝑁 𝑇 ]
Where U is a 2N×12 matrix, 𝑝 is a vector containing the 12 elements of 𝑃. Then, we set 𝑝 34 = 1 to decrease to 11 degrees of freedom. We need at least 6 non-planar points to resolve these equations. The solution to (2.8) is well known as the eigenvector of 𝑈 𝑇 𝑈 associated with the smallest eigenvalue. All the parameters can be extracted from 𝑃 as follows [8]:

𝒑 𝒊 = (𝑝 𝑖1 , 𝑝 𝑖2 , 𝑝 𝑖3 ) 𝑇 , 𝑖 = 1,2,3 𝑡 = 1/‖𝑝 3 ‖ 𝑢 𝑠 = 𝑡 2 × (𝒑 𝟏 • 𝒑 𝟑 ) 𝑣 𝑠 = 𝑡 2 × (𝒑 𝟐 • 𝒑 𝟑 ) 𝐹 = 𝑡 2 × ‖𝒑 𝟏 × 𝒑 𝟑 ‖ 𝒓 𝟏 = t/𝐹 × (𝒑 𝟏 -𝑢 𝑠 𝒑 𝟑 ) 𝒓 𝟐 = t/𝐹 × (𝒑 𝟐 -𝑣 𝑠 𝒑 𝟑 ) 𝒓 𝟑 = 𝑡 × 𝒑 𝟑 𝒕 𝒙 = 𝑡 × (𝑝 14 -𝑢 𝑠 )/𝐹 𝒕 𝒚 = 𝑡 × (𝑝 24 -𝑣 𝑠 )/𝐹 𝒕 𝒛 = 𝑡 𝐼 = [ 𝐹/𝑘𝑥 0 𝑢 𝑠 0 𝐹/𝑘𝑦 𝑣 𝑠 0 0 1 ] 𝑅 = ,𝒓 𝟏 , 𝒓 𝟐 , 𝒓 𝟑 -𝑇 𝑇 = [𝒕 𝒙 , 𝒕 𝒚 , 𝒕 𝒛 ] 𝑇 𝐸 = ,𝑅, 𝑇- (2.9)
𝑘𝑥, 𝑘𝑦 are the inverse of the pixel size. 𝐹 is the focal length. 𝑢 𝑠 , 𝑣 𝑠 are the principal point. 𝐼 is the intrinsic matrix. 𝐸 is the extrinsix matrix composed by the rotation matrix 𝑅 and translation vector 𝑇.

Iterative parametric method

An improvement of the classical projection matrix algorithm is the classical iterative parametric method [16]. The initial parameters are estimated by the projection matrix method. The refined parameters are obtained through the minimization of the projection error:

𝑃 = 𝐼𝐸 𝜀 𝑟 (𝑃) = 1 𝑛 ∑ ‖𝑃𝑀 𝑖 -𝑚 𝑖 ‖ 2 𝑛 𝑖=0 (2.10)
For the research results, there is little difference between the classical projection matrix method and this iterative method when the calibration markers are enough. However, this iterative method turned out to be superior when the markers are few [16]. The two classical methods are based on the classical helical phantom.

Multi-images method

Although the estimated projection matrix is known with a sub-pixel projection error, its decomposition into intrinsic and extrinsic parameters is known to be unstable [9]. Multi-images algorithm can provide more reliable intrinsic and extrinsic parameters. To reduce the statistical noise, calibration can be repeated with varying extrinsic parameters and fixed intrinsic parameters. Thereby, the inter-dependence between both sets of parameters is reduced, so that reliable intrinsic parameters are estimated.

For a given C-arm orientation, N images of the helix phantom are taken, with the phantom being moved in both rotation and translation between each image acquisition. This step is called multi-image calibration. The common intrinsic parameters I and extrinsic parameters ( 𝐸 𝑗 , 𝑗 = 1, … 𝑁 ) are then estimated simultaneously by minimizing the mean residual projection error 𝑅 𝑚 as follows, the initial parameters can be calculated by the classical projection matrix method:

𝑅 𝑚 = 1 𝑁 ⁄ ∑ 𝜀 𝑟 ( 𝑃 𝑗 ) 𝑁 𝑗=0 (2.11)
with

𝜀 𝑟 ( 𝑃 𝑗 ) = 1 𝑛 ∑ ‖ 𝑃 𝑗 𝑀 𝑖 -𝑚 𝑖 ‖ 2 𝑛 𝑖=0 𝑃 𝑗 = 𝐼 𝐸 𝑗 , 𝑗 = 1, … 𝑁
From (2.11), we obtain the robust and accurate intrinsic parameter 𝐼 and we choose one 𝐸 𝑗 to be the current extrinsic parameter.

Other methods

We introduce other two representative calibration algorithms in this section. These methods depend on their specific designed phantom.

a. Zhang's planar calibration algorithm

This algorithm [17] is one classical method in computer vision area. This technique was applied to the angiography acquisition system [21,22]. The method requires the acquisition of a planar pattern at a set of few orientations and the motion needs not to be known. Let 𝑚 ̃= ,𝑢, 𝑣, 1-𝑇 be an augmented 2D point, 𝑀 ̃= ,𝑋, 𝑌, 𝑍, 1-𝑇 is the 3D point, then the pinhole model can be represented as:

𝑠𝑚 ̃= 𝐼,𝑅 𝑡-𝑀 ̃ (2.12)
As the method adopts the planar phantom, the points with 𝑍 = 0 of the world coordinate can be denoted 𝑀 ̃= ,𝑋, 𝑌, 1-𝑇 . We rewrite (2.12) as:

𝑠𝑚 ̃= 𝐻𝑀 ̃ (2.
𝑉𝑑 = 0 (2.16)
Once 𝑏 is estimated, the intrinsic parameters that are considered to be constant and extrinsic parameters for each view can be obtained, see [17] for all details.

b. Mennessier's analytic algorithm

This algorithm proposed a direct analytic algorithm on a cubic phantom [15]. Each view is treated independently, and no restrictions are made on the position of the cone vertex, or on the position or orientation of the detector. The simple object consists of six points, two per orthogonal line. The algorithm can provide the unique analytic solution for our application in angiography system. See Fig. 2.4 in the following, they have the same imaging system.

First note that a location 𝑃 𝑟 is given by the (2.17) in the world coordinate system:

𝑷 𝒓 = 𝒔 + 𝐹𝒏 ̂+ (𝑢 ̃-𝑢 𝑠 )𝒖 ̂+ (𝑣 ̃-𝑣 𝑠 )𝒗 ̂ (2.17)
Every 3D point r that projects to the same detector location 𝑷 𝒓 belongs to the line passing through the source point 𝒔 and this detector location:

*𝒓|𝒓 = 𝒔 + 𝑡(𝑷 𝒓 -𝒔), 𝑡 ∈ 𝑅+ = *𝒓|𝒓 = 𝒔 + 𝑡((𝑢 ̃-𝑢 𝑠 )𝒖 ̂+ (𝑣 ̃-𝑣 𝑠 )𝒗 ̂+ 𝐹𝒏 ̂+ (2.18)
Then, the cone-beam projection of an arbitrary point r onto the detector is given by: 19) is converted into a set of simpler equations as follows, involving three intermediate unknown vectors 𝒂, 𝒃, 𝒄. The nine geometric parameters will be obtained directly from these three vectors.

𝑢 ̃(𝒓) = 𝑢 𝑠 + 𝐹 (𝒓 -𝒔)𝒖 (𝒓 -𝒔)𝒏 ̂ 𝑣 ̃(𝒓) = 𝑣 𝑠 + 𝐹 (𝒓 -𝒔)𝒗 (𝒓 -𝒔)𝒏 ̂ (2.19) (2.
First, the main simplification consists of shifting the origin of the detector pixel coordinates by:

𝑢 = 𝑢 ̃-𝑢 𝑐 𝑣 = 𝑣 ̃-𝑣 𝑐 (2.

𝑢 = 𝒂 • 𝒓

𝒄 • 𝒓 + 1 𝑣 = 𝒃 • 𝒓 𝒄 • 𝒓 + 1
(2.21) Consider a cubic calibration object with 6 points so that the raw detector images (𝑢 𝑐 , 𝑣 𝑐 ), (𝑢 ̃1, 𝑣 ̃1), … (𝑢 ̃6, 𝑣 ̃6) of the points 𝑟 0 , 𝑟 1 , … 𝑟 6 can be individually identified, where for some known positive constant k. 𝑟 1 = (𝑘, 0,0), 𝑟 2 = (-𝑘, 0,0) 𝑟 3 = (0, 𝑘, 0), 𝑟 4 = (0, -𝑘, 0) 𝑟 5 = (0,0, 𝑘), 𝑟 6 = (0,0, -𝑘)

(2.22)
After applying equation (2.22) in (2.21), we obtain six pairs of equations from which to determine 𝒂, 𝒃, 𝒄. See [15] for more details.

2.2.2Proposed Calibration Method

System geometry of cone-beam scanner

Mennessier et al. [15] proposed an analytical calibration algorithm from the cubic phantom. Xu et al. [23] applied the algorithm to the classical helical phantom. We extend this algorithm further to make it more robust and easier to implement. We build the world coordinate system on the phantom. 𝑧 axis is along the longitudinal direction of the phantom. We make a line from the 𝑛 0 th marker perpendicular to 𝑧 axis, the intersection is the origin 𝑶 𝒑 . 𝑥 axis is along this line, 𝑦 axis is the direction that subjects to right-handed coordinate system. The detector is a flat-panel detector without distortion and with known numbers of rows and columns of square pixels and known pixel size. See Fig. 2.4, the scanner geometry can be described by 9 parameters. Three scalars 𝑢 𝑠 , 𝑣 𝑠 , 𝐹 and four vectors 𝒖 ̂, 𝒗 ̂, 𝒏, ̂𝒔. They are all expressed in terms of the phantom system. Three unit vectors (𝒖 ̂, 𝒗 ̂, 𝒏 ̂) represent the detector orientation and they can be specified using three Euler angles. The vectors 𝒖 ̂, 𝒗 ̂ are parallel to the row and the column of the detector respectively. 𝒏 ̂ is the direction perpendicular to the detector plane. 𝒔 is the position of the X-ray source. A 3D point in the world system is denoted by 𝒓. The cone-beam projection of point 𝒓 in pixel units is (𝑢 ̃, 𝑣 ̃). Piercing point (𝑢 𝑐 , 𝑣 𝑐 ) is defined as the projection of the origin 𝑶 𝒑 . Principle point (𝑢 𝑠 , 𝑣 𝑠 ) is the orthogonal projection of the source position 𝒔 onto the detector. 𝐹 is the focal length, meaning the perpendicular distance from the source to the detector. For the calibration, the phantom is placed within the field of view. The longitudinal axis of the phantom is roughly parallel to the rotation axis of C-arm. The center of the phantom coincides roughly with the iso-center of the acquisition system. The phantom is imaged during the dynamic acquisition according to a detector orientation. The geometry parameters will be estimated at each acquisition angle independently and no requirement of the scanning motion and others.

Calibration phantom

The calibration phantom consists of 108 steel ball bearings (BBs) embedded in the plastic wall of a cylindrical phantom, see Fig. 2.4, the 108 spheres are made up of noncorrosive steel with two different radii arranged along a helical-shaped path, the large spheres have a diameter 3.2 mm and the small spheres have a diameter of 1.6 mm. The radius of the cylinder is 70 mm, there are 40 BBs spaced over the angle of 2π, the distance between two BBs is 1.3 mm, so the pitch of the helical curve is 52 mm.

As shown in Fig. 2.5, the distribution of BBs is designed with an 8 bit encoding according to a projection angle and phantom position. Indeed, we can label the alternating large and small spheres as a binary string. An arbitrary subsequence with the length more than 8 can provide enough information to identify the location of each sphere within the whole binary string. In this chapter, we set the secondary angle 𝛽 of the C-arm system to zero. From Fig. 2.5, we see that the projection of the spheres located from bottom to top. The large spheres have of course a bigger projection area than the small spheres. Most projection points are visible and some minor overlaps appear at the position with high curvature and noncircular. We use the encoding method to find the unique correspondence between each BB and their projection [24].

We skip the overlapped projection points and extract the centroid of each non-overlapped projection area. In this chapter, we adopt the ellipse fitting algorithm of Li et al. [25] to extract the centroid of the projected markers. We encode each extracted centroid by its projection area. The bigger area corresponds to 1 and the smaller one corresponds to 0. These subsequences can help us find the unique correspondence between the 3D and 2D. We label the 2D projected points by 𝑛 = 1, . .108.

Geometry estimation

We choose one work position of Zeego system. The center of the phantom is placed approximately near to the iso-center of the acquisition system. The longitudinal direction of the phantom is aligned roughly to the rotation axis.

As we described above, for each angle, we define the system geometry by 9 parameters, detector orientation (𝒖 ̂, 𝒗 ̂, 𝒏 ̂) which can be represented by Euler angles 𝜃, 𝛾, 𝜂, source to detector distance 𝐹, principal point (𝑢 𝑠 , 𝑣 𝑠 ), source position 𝒔. These parameters can be estimated by three intermediate vectors 𝒂, 𝒃, 𝒄. The projection coordinates (𝑢 ̃, 𝑣 ̃) and the vectors 𝒂, 𝒃, 𝒄 have the following relationship [15].

𝑢 = 𝒂 • 𝒓 𝒄 • 𝒓 + 1 𝑣 = 𝒃 • 𝒓 𝒄 • 𝒓 + 1
(2.23) where 𝑢, 𝑣, 𝒂, 𝒃, 𝒄 are defined as:

𝑢 = 𝑢 ̃-𝑢 𝑐 𝑣 = 𝑣 ̃-𝑣 𝑐 (2.24) 𝒂 = (𝑢 𝑠 -𝑢 𝑐 )𝒏 ̂+ 𝑓𝒖 -𝒏 ̂• 𝒔 𝒃 = (𝑣 𝑠 -𝑣 𝑐 )𝒏 ̂+ 𝑓𝒗 -𝒏 ̂• 𝒔 𝒄 = 𝒏 -𝒏 ̂• 𝒔 (2.
25) all the geometry parameters can be obtained by known intermediate vectors 𝒂, 𝒃, 𝒄:

𝜀 𝑓 = -𝑠 𝑝 𝑠𝑔𝑛(𝒂 × 𝒃 • 𝒄) 𝒖 ̂= 𝜀 𝑓 (𝒃 × 𝒄) ‖𝒃 × 𝒄‖ 𝒗 ̂= -𝜀 𝑓 (𝒂 × 𝒄) ‖𝒂 × 𝒄‖ 𝒏 ̂= -𝜀 𝑓 (𝒄) ‖𝒄‖ 𝐹 = 𝜀 𝑓 ‖𝒂 × 𝒄‖ ‖𝒄‖ 2 = 𝜀 𝑓 ‖𝒃 × 𝒄‖ ‖𝒄‖ 2 𝑢 𝑠 = (𝒂 • 𝒄) ‖𝒄‖ 2 + 𝑢 𝑐 𝑣 𝑠 = (𝒃 • 𝒄) ‖𝒄‖ 2 + 𝑣 𝑐 𝒔 = -𝑠 𝑝 ‖𝒄‖ ( 𝒂 • 𝒄 ‖𝒂 × 𝒄‖ 𝒖 ̂+ 𝒃 • 𝒄 ‖𝒃 × 𝒄‖ 𝒗 ̂-𝜀 𝑓 𝒏 ̂)
(2.26) 𝑠 𝑝 = -1 for the cone-beam imaging system. 𝒂, 𝒃, 𝒄 were estimated by the projection of the points located on the coordinate axes of the cubic phantom [15]. Our proposed algorithm applies these formulations to the helical phantom. Compared with the algorithm of Xu et al. [23], the proposed algorithm can estimate 𝒂, 𝒃, 𝒄 more robustly and eliminate the degeneration correction procedure of Xu's method.

a. Estimation of 𝒂, 𝒃, 𝒄

To estimate 𝒂, 𝒃, 𝒄 robustly, we try to estimate (𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥 ), (𝑎 𝑦 , 𝑏 𝑦 , 𝑐 𝑦 ), (𝑎 𝑧 , 𝑏 𝑧 , 𝑐 𝑧 ) separately. We use the same way as Xu [23] to calculate the projection of 𝑧 axis and 𝑢 𝑐 , 𝑣 𝑐 . We have the coordinate of 𝑛 𝑡 marker on the phantom as according to the helical arrangement formed a circular helix of radius 𝑅 and pitch 2 𝑏 0 : The label of each marker is:

𝒓 𝒏 = (𝑅 cos 𝑛 , 𝑅 sin 𝑛 , 𝑏 0 * 𝑛 ) with 𝑛 = (𝑛 -𝑛 0 ) * 2 𝑁 ⁄ 𝑏 0 = 𝑁 2 ⁄ *
𝑛 𝐴 = 𝑛, 𝑛 𝐵 = 𝑛 + 𝑁 𝑛 𝐶 = 𝑛 + 3𝑁 2 ⁄ 𝑛 𝐷 = 𝑛 + 𝑁 2 ⁄ (2.29) (𝑛 𝐴 , 𝑛 𝐵 , 𝑛 𝐶 , 𝑛 𝐷 ) ∈ (1, 𝑁𝑃), so we should seek 𝑛 ∈ (1, 𝑁𝑃 -3𝑁 2 ⁄ ).
There is a geometric basis of perspective projection. The projection of the intersection point of two lines is exactly the intersection point of the two projection lines. As seen in Fig. 2.6, the projection of 𝑟 𝐴 , 𝑟 𝐵 , 𝑟 𝐶 , 𝑟 𝐷 are 𝑟 𝐴′ , 𝑟 𝐵′ , 𝑟 𝐶′ , 𝑟 𝐷′ respectively, and the intersection point of 𝑟 𝐴 𝑟 𝐶 and 𝑟 𝐵 𝑟 𝐷 is:

𝒓 𝒏 𝒁 = (0,0, 𝑏 0 * ( 𝑛 + 3 2 ⁄ )) (2.30)
The intersection point of 𝑟 𝐴′ 𝑟 𝐶′ and 𝑟 𝐵′ 𝑟 𝐷′ is 𝑟 𝑛 𝑍′ , these points are the projection of 𝒓 𝒏 𝒁 and compose the projection of the 𝑧 axis. 𝑟 𝑛 𝑍′ can be calculated by the intersection of the two projection lines 𝑟 𝐴′ 𝑟 𝐶′ and 𝑟 𝐵′ 𝑟 𝐷′ . Given a set of 𝒓 𝒏 𝒁 = (0,0, 𝑍 𝑛 𝑍 ) and 𝑟 𝑛 𝑍′ = (𝑢 𝑛 𝑍′ , 𝑣 𝑛 𝑍′ ), equation (2.21) could be written as:

𝐴 𝑧 𝑝 𝑧 = 𝑏 (2.31)
Where

𝐴 𝑧 = * 𝑍 𝑛 𝑍 0 0 𝑍 𝑛 𝑍 -𝑢 𝑛 𝑍′ 𝑍 𝑛 𝑍 -𝑣 𝑛 𝑍′ 𝑍 𝑛 𝑍 … … … … … … + 𝑏 = ,𝑢 𝑛 𝑍′ 𝑣 𝑛 𝑍′ … …-′
We adopt the SVD method [26] to resolve (2.31). The steps are: (i) Find the SVD, 𝐴 𝑧 = 𝑈𝐷𝑉 𝑇 (ii)Set 𝑏 ′ = 𝑈 𝑇 𝑏 (iii)Find the vector 𝑦 defined by 𝑦 𝑖 = 𝑏 𝑖 ′ /𝑑 𝑖 , 𝑖 =1,2,3, where 𝑑 𝑖 is the 𝑖 𝑡 diagonal entry of 𝐷. (iv)The solution is 𝑝 𝑧 = 𝑉𝑦. a.2 Estimation of 𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥 a.2.1 Projection of the axes paralleled to 𝑥 axis See Fig. 2.7, we propose to calculate the projection of two axes paralleled to 𝑥 axis, differently from calculating the projection of 𝑥 axis proposed by Xu [23].This will not only improve the estimation accuracy of 𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥 , but also resolve the degeneration problem in Xu [23].

We try to find another group of parallelograms, each one is composed by 4 markers, see Fig. 2.7, we have the following expressions: 𝒓 𝑨 = (𝑅 cos 𝐴 , 𝑅 sin 𝐴 , 𝑏 0 * 𝐴 ) 𝒓 = (𝑅 cos 𝐵 , 𝑅 sin 𝐵 , 𝑏 0 * 𝐵 ) 𝒓 𝑪 = (𝑅 cos 𝐶 , 𝑅 sin 𝐶 , 𝑏 0 * 𝐶 ) 𝒓 𝑫 = (𝑅 cos 𝐷 , 𝑅 sin 𝐷 , 𝑏 0 * 𝐷 )

𝐴 = 𝑖 𝐵 = 𝑖 + 2 𝑘 𝐶 = -𝑖 + 2 𝑘 𝐷 = -𝑖 (2.32)
The label of each marker is: where The label of each marker is:

𝑛 𝐴 = 𝑖 𝑛 𝐵 = 𝑖 + 𝑁𝑘 𝑛 𝐶 = 2𝑛 0 -𝑖 + 𝑁𝑘 𝑛 𝐷 = 2𝑛 0 -𝑖 (
𝐴 𝑥 = [ 𝑋 𝑖 1 𝑋 0 0 𝑋 𝑖 1 𝑋 -𝑢 𝑖 1 𝑋′ 𝑋 𝑖 1 𝑋 -𝑣 𝑖 1 𝑋′ 𝑋 𝑖 1 𝑋 … 𝑋 𝑖 2 𝑋 0 … … 0 𝑋 𝑖 2 𝑋 … … -𝑢 𝑖 2 𝑋′ 𝑋 𝑖 2 𝑋 -𝑣 𝑖 2 𝑋′ 𝑋 𝑖 2 𝑋 … ] 𝑏 = [ 𝑏 𝑢𝑖 1 𝑏 𝑣𝑖 1 … 𝑏 𝑢𝑖 2 𝑏 𝑣𝑖 2 … ] ′ 𝑏 𝑢𝑖 1 =
𝑛 𝐴 = 𝑗 𝑛 𝐵 = 𝑗 + 𝑁𝑘 𝑛 𝐶 = 2𝑛 0 -𝑗 + 𝑁(2𝑛 + 1)𝑘/2 𝑛 𝐷 = 2𝑛 0 -𝑗 + 𝑁(2𝑛 -1)𝑘/2
(2.37)

𝑘 = 1, 𝑛 = 0,1, (𝑛 𝐴 , 𝑛 𝐵 , 𝑛 𝐶 , 𝑛 𝐷 ) ∈ (1, 𝑁𝑃). 𝑘 = 1, 𝑛 = 0, 𝑗 1 ∈ ,max (1,2𝑛 0 + 𝑁/2 -𝑁𝑃, 𝑛 0 -𝑁/4 + 1), min (𝑁𝑃 -𝑁, 2𝑛 0 - 𝑁/2 -1)-, 𝑗 1 ≠ 𝑛 0 + 𝑁/4. 𝑘 = 1, 𝑛 = 1, 𝑗 2 ∈ ,max (1,2𝑛 0 + 3𝑁/2 -𝑁𝑃, 𝑛 0 + 𝑁/4 + 1), min (𝑁𝑃 -𝑁, 2𝑛 0 + 𝑁/2 -1)-, 𝑗 2 ≠ 𝑛 0 + 3𝑁/4. 𝑘 = -1, 𝑛 = 0, 𝑗 3 ∈ ,max (1 + 𝑁, 2𝑛 0 + 𝑁/2 -𝑁𝑃), min (𝑁𝑃, 2𝑛 0 -𝑁/2 -1, 𝑛 0 + 𝑁/4 -1)-, 𝑗 3 ≠ 𝑛 0 -𝑁/4 𝑘 = -1, 𝑛 = 1, 𝑗 4 ∈ ,max (1 + 𝑁, 2𝑛 0 -𝑁/2 -𝑁𝑃), min (𝑁𝑃, 2𝑛 0 -3𝑁/2 - 1, 𝑛 0 -𝑁/4 -1)-, 𝑗 4 ≠ 𝑛 0 -3𝑁/4
The intersection point of 𝑟 𝐴 𝑟 𝐶 and 𝑟 𝐵 𝑟 𝐷 is: These projection of four axes provide more data. We use the same way as described in a.2 to calculate 𝑎 𝑦 , 𝑏 𝑦 , 𝑐 𝑦 .

𝒓 𝒋 𝟏 𝒀 = (0,

b. Degeneration elimination

The algorithm of Xu [23] may have two kinds of degeneration. First, the projection of 𝑥 axis may degenerate to one point. Second, the projection of the parallelogram may degenerate to a line. The degeneration correction is troublesome, decreasing the efficiency.

For the first kind of degeneration, the projection of x axis may degenerate to one point, then, 𝑎 𝑥 = 𝑏 𝑥 = 0, 𝑐 𝑥 can't be determined. Now, as seen in a.2, this problem is tackled by the projection of the axes paralleled to x axis.

For the second degeneration, it may reduce the available intersection points that used to determine the calibration parameters. This kind of degeneration only exists when we calculate the projection of the axes paralleled to x axis and y axis. Now, as seen in a.2 and a.3, we have more projection points to estimate 𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥 and 𝑎 𝑦 , 𝑏 𝑦 , 𝑐 𝑦 . Therefore, we propose a simple threshold method to deal with it. If 𝑑(𝑟 𝐴′ 𝑟 𝐵′ , 𝑟 𝐶′ 𝑟 𝐷′ ) < 𝑑 𝑇 , the second degeneration occurs, the correspondent quadrilaterals would not be adopted. [19] first proposed an ideal C-arm movement model. In Fig. 2.9, 𝑂 𝑠 is the X-ray source, α, 𝛽 are rotation axes, the world coordinate system 𝑋 𝑤 -𝑌 𝑤 -𝑍 𝑤 locates on the iso-center o, the camera coordinate system 𝑋 𝑠 -𝑌 𝑠 -𝑍 𝑠 . They have several ideal hypotheses: 1, 𝛼 axis and 𝛽 axis are orthogonal and intersect at the iso-center o. 2, the central X-ray beam passes the iso-center o and is orthogonal with the α axis and 𝛽 axis. 3, the projections of the 𝛼 axis and 𝛽 axis in the imaging plane are parallel with the horizontal u and vertical axis v of the image respectively. 4, the principle point (𝑢 𝑠 , 𝑣 𝑠 ) is at the center of the image plane and fixed. The focal length SDD and the distance from source to the rotation axis (SAD) are fixed.

This model represents the ideal acquisition geometry. We have the formulation of the intrinsic matrix 𝐼 and extrinsic matrix 𝐸 [19,27]: In practice, however, none of the ideal hypotheses for simplification purposes can be guaranteed due to the mechanical flaws. Therefore, Dumay's model contains systematic errors that cannot be ignored.

𝐼 = [ 𝑆𝐷𝐷 0 𝑢 𝑠 0 𝑆𝐷𝐷 𝑣 𝑠 0 0 1 ] 𝑅 𝛼 = [ 0 𝑠𝑖𝑛(𝛼) -cos(𝛼) 0 -cos(𝛼) -𝑠𝑖𝑛(𝛼) -1 0 0 ] 𝑅 𝛽 = [ 1 0 0 0 cos (β) -𝑠𝑖𝑛(β) 0 𝑠𝑖𝑛(β) cos (β) ]

Gorges's model

Gorges et al. [28] adopted the multi-image calibration algorithm [9] to assess the physical behavior of the C-arm for medical augmented reality application. From the knowledge of the main characteristics of the C-arm, realistic models of the acquisition geometry are proposed. As seen in Fig. 2.10, the main adopted characteristics of the C-arm system are: 1, the focal length SDD=F is constant, the principle point𝑂 𝑖 = (𝑢 𝑠 , 𝑣 𝑠 ), 𝑢 𝑠 changes with respect to the 𝛼 axis rotation, 𝑣 𝑠 is constant. The SDD and the angles 𝛼, 𝛽 measured by sensors are reliable. 2, the two rotation axes are stable and independent. 3, the motion around 𝛼 axis has a residual translation that proportional to this rotation axis vector.

The expression of the intrinsic I and extrinsic matrices E are: The refinement of the parameters can be obtained by minimizing the projection error:

𝐼 = [ 𝑆𝐷𝐷 0 𝑢 𝑠 0 𝑆𝐷𝐷 𝑣 𝑠 0 0 1 ] 𝑢 𝑠 = ∑ 𝑢 𝑖 3 𝑖=0 𝛼 𝑖 𝐸 = 𝐸 0 𝐸 𝛽 𝐸 𝛼 𝐸 𝛽 = [𝑅 𝛽 |(𝐼 -𝑅 𝛽 )𝑃 𝛽 ] λ = ∑ 𝜆 𝑖 3 𝑖=0 𝛼 𝑖
𝐺 = ∑ ‖𝑥 𝑖 𝑗 -𝑃𝑋 𝑖 𝑤 ‖ 𝑖,𝑗 (2.44)
i is the index of control points, j is the index of orientations, 𝑋 𝑖 𝑤 and 𝑥 𝑖 𝑗 are the world coordinates of the control points and its 2D projection position. The multi-images algorithm [9] of Gorges's model is too tedious to implement.

Li's model

Li et al. [22] adopted the posteriori information from Gorges [9] and proposed a new C-arm movement model objectively and systematically. See Fig. 2.11, 𝑂 𝑠 is the X-ray source, 𝛽 axis passes o and is perpendicular to the plane of circle o, there are world coordinate system 𝑋 𝑤 -𝑌 𝑤 -𝑍 𝑤 , camera coordinate system 𝑋 𝑠 -𝑌 𝑠 -𝑍 𝑠 , image coordinate system 𝑋 𝑖 -𝑌 𝑖 -𝑍 𝑖 , movement coordinate system 𝑋 𝑎 -𝑌 𝑎 -𝑍 𝑎 .

For this C-arm movement model, there were several new hypotheses: 1, 𝛼 axis and 𝛽 axis may neither intersect with each other nor be perpendicular. 2, the central X-ray beam may not pass o and naturally not be perpendicular to the 𝛼 axis and 𝛽 axis. 3, there is no direct connection between the image axes u v and the projection of the 𝛼 axis and 𝛽 axis. 4, the principal point 𝑂 𝑖 = (𝑢 𝑠 , 𝑣 𝑠 ), it changes with respect to the C-arm orientation. The focal length SDD is considered to be constant. (2.45)

𝑎 𝑗, 𝑏 𝑗 are the coefficients of the polynomial expansion of the principal point. See Fig. 2.12, 𝐸 0 is the transformation from world system to the reference coordinate system. 𝐸 𝛽 and 𝐸 𝛼 represent the rotation transformation about 𝛼 and 𝛽 axes, 𝐸 𝑎𝑠 is the transformation from movement system to camera system, 𝑟 is the radius of the circle. All of these parameters can be calibrated by this C-arm movement model 𝑀, see [22] for more details.

The refinement of the estimation can be obtained by minimizing the projection error:

𝐺 = ∑ ‖𝑥 𝑖 𝑗 -𝑃𝑋 𝑖 𝑤 ‖ 𝑖,𝑗
With 𝑃 = 𝐼𝐸 The movement model is:

𝑀 = {𝑟, 𝑟 𝛼 , 𝑃 𝛼 , 𝐸 𝑎𝑠 , 𝐸 0 , 𝐹, 𝑎 𝑗𝑗 , 𝑏 𝑗𝑗 } (2.46)
i is the index of control points, j is the index of orientations, 𝑋 𝑖 𝑤 and 𝑥 𝑖 𝑗 are the world coordinates of the control points and its 2D projection position.

There are many other proposed movement models. We compared 4 representative movement models in Table .2.1. We found that Gorges's model is easy to perform and the posterior information is more reliable. The disadvantage of Gorges's model is that the model only considers the traditional C-arm with a deviation along the rotation axis. Next, we proposed our movement model that will be suitable for the Zeego system.

Proposed Movement Model

2.3.2.1Recalculation to C-arm coordinate system

In part 2.2, the acquisition geometry is under the phantom coordinate system. The extrinsic matrix can be transformed to the nominal iso-center system (C-arm system), which is independent of the phantom placement. See Fig. 2.13, we see the geometry of Zeego system. This six-joint robot has a broader coverage and multi work positions. We plot 3 representative work positions: Head side, Left side, Table 30 in Fig. 2.13 (a), (b) and (c). In Fig. 2.13, we plot all the relative coordinate systems. The origin of the phantom coordinate system is 𝑶 𝒑 , the coordinate axes are defined as described in Dumay's model [19] Cañero's model [21,29] Gorges's model [9,28] Li's model [ 2.2.2.1. The table system is defined by the Zeego system, the origin 𝑶 𝒕 locates on the edge of the table, 𝒙 𝒕 points to the right side of the patient along the edge of the table, 𝒚 𝒕 is perpendicular to the table plane, 𝒛 𝒕 subjects to the right-handed system. The C-arm system with the origin is the nominal iso. 𝒙 𝒄 is antiparallel to 𝒚 𝒕 , 𝒚 𝒄 is parallel to 𝒙 𝒕 , 𝒛 𝒄 is parallel to 𝒛 𝒕 . The C-arm system will be coordinated with the operating table automatically. X-ray source 𝑶 𝒔 (same as 𝒔 defined above), iso, principle point (𝑢 𝑐 𝑐 , 𝑣 𝑐 𝑐 ) are collinear in the condition of ideal geometry. The source system locates at the X-ray source 𝑶 𝒔 , 𝒙 𝒔 is parallel to 𝒚 𝒄 , 𝒚 𝒔 is antiparallel to 𝒛 𝒄 , 𝒛 𝒔 is antiparallel to 𝒙 𝒄 . The detector coordinate system is defined as shown in Fig. 2.13 (d). See Fig. 2.13 (d), the ideal trajectory plane is a circle with the origin iso, the radius is the source to axis SAD, 𝒛 𝒄 is the normal vector of this plane. For the Head position, the table is placed in the horizontal direction, the plane of the C shaped gantry is perpendicular to the table plane. For the Left position, the plane of the C shaped gantry rotates, the table is fixed. For the Table30 position, the table plane is rotated by 30 °, and the C gantry plane also rotates so that the two planes are perpendicular.

We define 𝑹 𝒑 𝒄 , 𝑻 𝒑 𝒄 as the rotation and translation from C-arm system to phantom system. We have the following approximation.

The apparent displacement of the phantom as the gantry rotates corresponds to: 𝑑𝑥 = 𝑡 1 cos 𝛼 + 𝑡 2 sin 𝛼 𝑑𝑦 = -𝑡 1 sin 𝛼 + 𝑡 2 cos 𝛼 𝑑𝑧 = 𝑡 3 (2.47) 𝑡 1 , 𝑡 2 , 𝑡 3 are the components of the translation 𝑻 𝒑 𝒄 .

The magnification is defined as:

𝑆𝐴𝐷 𝑐 = 𝑆𝐷𝐷 -𝐷𝐴𝐷 𝑐 𝑀 𝑐 = 𝑆𝐷𝐷 * (𝑆𝐴𝐷 𝑐 ) -1 (2.48)
The superscript means the parameter under this system. DAD means detector to gantry axis distance, SAD means source to gantry axis distance, SDD source to detector distance.

The piercing point (𝑢 𝑐 𝑐 , 𝑣 𝑐 𝑐 ) is defined as:

𝑢 𝑐 𝑐 = 𝑢 𝑐 𝑝 + 𝑀 𝑐 • 𝑑𝑧 𝑣 𝑐 𝑐 = 𝑣 𝑐 𝑝 + 𝑀 𝑐 • 𝑑𝑦 (2.49)
The source position under C-arm system is defined as:

𝑷 𝑺 𝒄 = 𝑹 𝒑 𝒄 𝑷 𝑺 𝒑 + 𝑻 𝒑 𝒄 (2.50)
We minimize the apparent sinusoidal oscillations in detector offset and the apparent offset and tilt of the source trajectory plane to estimate the C-arm system. The objective function is: 𝑬 𝒄 : extrinsic matrix in C-arm system. T means the transposition.

𝑔 = 𝜎 𝑣 𝑐 2 +

Modeling C-arm movement

We propose different movement models based on the estimated geometry. As in Fig. 2.13, we have three representative work positions of the 3D acquisition: Head, Left and Table 30. The orientation is classically described by two angles: 𝛼 right/left anterior orientation (RAO/LAO), 𝛽 cranio-caudal (CC). The two angles can be measured in real-time by sensors. For this chapter, we set 𝛽 = 0.

The projection matrix at angle 𝛼 can be modeled as:

𝑷 𝜶 =
𝑢 𝑠 (𝛼) = ∑ 𝑝𝑢 𝑖 𝛼 𝑁-𝑖 𝑁 𝑖=1 𝑣 𝑠 (𝛼) = ∑ 𝑝𝑣 𝑖 𝛼 𝑁-𝑖 𝑁 𝑖=1
(2.58) 𝑝𝑢, 𝑝𝑣 are the vector coefficients of the polynomial expansion, N is the order.

𝑫 𝜶 can be modeled as a mean rotation motion for the Head and Left work positions:

𝑫 𝜶 = ,𝑹 𝜶 |(𝑰 -𝑹 𝜶 )𝑷 𝜶 - (2.
𝜃 0 , 0 , 𝛾 0 are three Euler angles of 𝑅 0 . 𝑇 0 is the translation vector of 𝐸 0 . 𝛼, 𝐹 are the ideal parameters read from sensor data.

Calibrating the models

The initial parameters can be estimated by the above parts. The accurate model parameters can be calibrated by minimizing the projection error𝐸 𝑓 .

For the models, we have the follows:

𝑘=0,1,2 = 𝑎𝑟𝑔 𝑚𝑖𝑛(𝐸 𝑓 )
With

𝐸 𝑓 = ∑ 𝜀(𝑃 𝑗 ) 𝑗𝜖𝑎𝑛𝑔𝑙𝑒 𝑛 ε(𝑃 𝑗 ) = ∑ ‖𝑃 𝑗 𝑋 𝑖 𝑐 -𝑞 𝑖 ‖ 𝑚 𝑖=0

𝑚

(2.64) 𝑃 𝑗 is the projection matrix build from the movement models. 𝑋 𝑖 𝑐 means the 3D coordinates of each marker under C-arm system. 𝑞 𝑖 is the 2D projection coordinates of each marker. 𝑚 is the number of markers, 𝑛 is the number of the calibration angles. 𝑘=0,1,2 mean the three models. We use Levenberg-Marquardt algorithm [31] to realize the optimization.

Implementation of the proposed movement models

We perform the acquisition of the phantom from some projection angles by a sampling interval 𝛼 𝑇 . First, we calculate the geometry by the proposed calibration algorithm. Second, we recalculate the geometry to the C-arm system. Third, we calibrate the proposed movement models. Finally, we generate all the needed geometries of the 3D reconstruction by the movement models.

Experimental Result

The experiments were performed on three representative imaging positions of Artis-Zeego: Head, Left, Table30. We set the focal length SDD = 𝐹 = 1199 mm, SAD = 785 mm, pixel size in each direction 𝑘 𝑥 = 𝑘 𝑦 = 0.616 mm/pixel. The projection image has the resolution of 480×616 pixels. The protocol for the 3D reconstruction is 8sDR. We divide the whole dataset into two parts, the calibration dataset and the test dataset. The sampling interval 𝛼 𝑇 =10 °.

The calibration dataset is defined as:

Geometry estimation

Geometry estimation of the proposed algorithm

We apply the proposed helical algorithm and Xu's algorithm [23] to the calibration dataset. We compare the two algorithms in two aspects: the precision and the efficiency. 𝑑 𝑇 = 2.

We adopt the mean projection error to evaluate the accuracy of the estimated geometries. The mean projection error (mpe) of each angle is calculated as:

𝑚𝑝𝑒 = ∑ ‖𝑃 𝑗 𝑋 𝑖 𝑝 -𝑞 𝑖 ‖ 𝑚 𝑖=0 𝑚 (2.

65) 𝑋 𝑖

𝑝 means the 3D coordinates of each marker under phantom system. 𝑞 𝑖 is the 2D projection coordinates of the center of each marker. 𝑚 is the number of markers. 𝑃 𝑗 is estimated projection matrix by the two algorithms. We plot the mpe over all the 2.2, we can see that the proposed algorithm has a lower projection error than Xu's algorithm. Meanwhile, the proposed algorithm decreases the computation time as it eliminates the degeneration correction procedure.

Intrinsic parameter analysis

We calculate the mean and standard deviation of each intrinsic parameter 𝑢 𝑠 , 𝑣 𝑠 , 𝐹 in Table 2.3. We use the polynomial function to fit the parameters according to different order. The fitting errors are also showed in Table 2.3.

Observing the results in table 2.3, we found that the deviation of 𝐹 is quite small related to its length, we can consider it to be constant. 𝑢 𝑠 , 𝑣 𝑠 are variable, considering the fitting error and the complexity of the order, we choose the third order fitting polynomial.

Estimation of C-arm coordinate system

Estimation of C-arm system

We set 𝐿 𝑡 =0.6 mm. We plot the source trajectories both under the phantom coordinate system and C-arm coordinate system in Fig. 2.15. Left column is the source trajectory under phantom system. Right column is the source trajectory under C-arm system, the blue represent the estimated source position, the red one is the ideal circle trajectory.

We can see that the source trajectory plane under C-arm system has a tilt angle to the ideal circle trajectory plane for Head and Left position. For Table position, the source trajectory under C-arm system has an oscillation, not located in the same plane.

Extrinsic parameter analysis

We analyze the extrinsic parameters related to C-arm system. We have the equation to describe it:

𝑬 = 𝑬 𝟎 𝑫 𝜶 (2.66)
𝑬 𝟎 is the extrinsic matrix at the reference position (𝛼 = 𝛽 = 0). 𝑫 𝜶 means the movement around 𝛼 axis. 𝑬, 𝑬 𝟎 can be estimated by the proposed method. 𝑫 𝜶 can be obtained by 𝑬 𝟎 -1 𝑬.

We suppose the motion 𝑫 𝜶 can be modeled as a mean rigid motion, defined as:

𝑫 𝜶 = ,𝑹 𝜶+𝜶 𝒆𝒓𝒓 |(𝑰 -𝑹 𝜶+𝜶 𝒆𝒓𝒓 )𝑷 𝜶 - (2.67)
𝑹 𝜶 is the rotation matrix, which can be represented by the rotation vector 𝒓 𝜶 and sensor angle 𝛼 by Rodrigue formulation. The translation 𝑻 𝜶 of 𝑫 𝜶 can be represented as 𝑹 𝜶 and 𝑷 𝜶 . 𝑷 𝜶 is a fixed point on the rotation vector 𝒓 𝜶 . 𝑷 𝜶 is a fixed point on the rotation vector 𝒓 𝜶 . We calculate the mean and standard deviation of the parameters 𝛼 𝑒𝑟𝑟 , 𝒓 𝜶 , 𝜃 𝛼 , 𝑷 𝜶 in table 2.4. 𝛼 𝑒𝑟𝑟 is the residual angle between the current calculated angle and the sensor angle. 𝒓 𝜶 is the rotation vector of current calibration angles. 𝜃 𝛼 is the deviation angle of the current rotation vector 𝒓 𝜶 to 𝒓 ̅ 𝜶 , 𝒓 ̅ 𝜶 is the mean rotation vector of all calibration angles. 𝑷 𝜶 is the fixed point on the current rotation vector 𝒓 𝜶 .

See table 2.4, 𝛼 𝑒𝑟𝑟 is small enough, so that the rotation angle is receivable by the sensor data. 𝒓 𝜶 is stable and 𝜃 𝛼 has small mean and standard deviation. The two values prove that the assumption of rotation matrix 𝑹 𝜶 is correct. 𝑷 𝜶 has a big deviation, we should consider more about the translation part.

Test of different possible movement models

We adopt the calibration dataset to calculate different movement models, and apply the calibrated movement models to the test dataset. Here, we consider intrinsic and extrinsic parameters respectively. We have 4 models: 𝑚 𝑖 , 𝑚 0 , 𝑚 1 , 𝑚 2 . 𝑚 𝑖 is the ideal geometry of the acquisition trajectory. 𝑚 0 has the ideal intrinsic parameters and the rigid extrinsic parameters. 𝑚 1 has the variable principle point and the rigid extrinsic parameters. 𝑚 2 has the variable principle point and the rigid extrinsic parameters added with residual translation.

We have the following formulations: We calculate the projection error by these models in table 2.5. For the Head position, 𝑚 1 has smaller error than others, 𝑚 𝑖 and 𝑚 0 has acceptable small error. This may indicate that the acquisition trajectory is stable and more close to a circle. For the Left side position, 𝑚 1 is the best model. 𝑚 𝑖 has much more error than Head position. This illustrates that the deviation results in the variation of the intrinsic parameters. For the last Table30 work position, this is the most complex motion. 𝑚 𝑖 ,𝑚 0 can't provide the acceptable error, 𝑚 2 can compensate this deviation. This indicates that there are some uncertain residual deviations. Next, we apply 𝑚 1 to Head and Left side position, 𝑚 2 to the Table30 position.

𝑚 𝑖 𝑖 =*

Geometry deviation of Zeego system

If the Zeego system has no geometry deviation, the obtained 3D object by DynaCT and the ideal geometry and the acquired projection are consistent [32]. Here, we calculate the 2D projection error observed for the real motion of Zeego system.

We generate a simulated 2D projection by the system matrix. The system matrix is generated by the simplified distance-driven [7] projector with the ideal geometry. Second, we extract the 2D centroid coordinates of the projected spheres from the projection images. They are the generated simulated projection and the clinical acquired projection. Last, we calculate the mean error between the two extracted 2D coordinates of each angle.

See Fig. 2.16, we plot the 2D projection error in horizontal (u) and vertical (v) directions. We found that the above consistence assumption is invalid. The acquisition system has obvious geometry deviations related to the ideal geometry. The deviation of Head position is minor. The machine motion along u direction has an obvious smooth deviation according to the rotation angles for the Left position. The geometry deviation has an oscillation for the Table position, especially along v direction.

We calculate the 2D projection error for the calibration and test angles. We show the results of some angles in table 2.6. For Head position, the mean error is 1.44 pixels (equals to 0.89 mm), and the biggest error is 1.94 pixels (equals to 1.19 mm). This means that Head position has the least deviation. For Left position, the mean error is 1.98 pixels (equals to 1.22 mm), and the biggest error is about 3.31 pixels (equals to 2.04 mm).This means that the Left position may have more system deviation. For Table30 position, the mean error is 4.85 pixels (equals to 2.99 mm), and the biggest error is about 10.02 pixels (equals to 6.17 mm). This means that the deviation for this position is the most severe.

Evaluation of the proposed movement models

3D reconstruction

We reconstruct the spheres by the estimated geometries and do the comparisons. We adopt the clinical acquired projections to realize the reconstruction. The acquisition angles include both calibration and test angles (39 angles). The projection matrix method [8] is a classical and widely used method. Here, the comparisons include three geometries, the ideal geometry, the estimated geometry by projection matrix method, the estimated geometry generated by the proposed movement model method.

First, we remove the background by some morphological operators and extract the projection of the spheres. Second, the system matrices for the 3D reconstruction are generated by the ideal and the two estimated geometries. Last, we use the iterative reconstruction method [1] to realize the reconstruction. This method adopts the MAP with 𝐿 0 prior for the ill-posed problem.

The ray-casting algorithm [33] is the classical volume rendering method. We adopt it to visualize the 3D reconstruction of the spheres. All the rendering parameters are identical for each visualization result. See Fig. 2.17 From Fig. 2.17 (a), we found that there are no obvious differences among the three geometries in (a). The Head position is close to an ideal circle. From Fig. 2.17 (b), we see that for the Left position, the 3D object without calibration (ideal geometry) has an obvious artifacts and degenerated spatial resolution. The 3D object with the estimated geometries eliminates the artifacts and improves the spatial resolution. The proposed method has comparable effect as the projection matrix method. From Fig. 2.17 (c), we found that for the Table position, it has the most severe geometry deviation. The 3D object without calibration has much artifacts and low spatial resolution. The estimated geometries by the proposed correct the geometry deviation well.

Efficiency of the algorithms

We calculate the computation time to generate all the projection matrices along the acquisition trajectory by the step of 2 °. The ideal projection matrices and the estimated projection matrices by the classical method are calibrated at each acquisition angles. The proposed movement models are calibrated by the calibration dataset, then, all the projection matrices are generated respectively by these models.

In Table 2.7, we see that the proposed movement model method has less computation time than the classical projection matrix method.

As a conclusion, we can prove that the proposed movement model method can estimate and predict the geometries with acceptable accuracy. Meanwhile it can increase the efficiency of the calibration process.

Discussion

We proposed a complete data acquisition modeling method of Artis-Zeego C-arm CBCT system. First, we proposed an extended calibration algorithm on a classical helical phantom. The phantom is easy to obtain in the experiment environment. For the helical algorithm, we propose to calculate the projection of the axes paralleled to 𝑥 axis and more projections of the axes paralleled to 𝑦 axis. These will not only improve the estimation accuracy of 𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥 and 𝑎 𝑦 , 𝑏 𝑦 , 𝑐 𝑦 , but also resolve the degeneration correction problem. The experiment proves that the proposed algorithm is more accurate and increases the efficiency than the original algorithm [23].

In this work, we set the secondary angle 𝛽 to be zero. However, the angle 𝛽 is often to be a non-zero angle in the clinical acquisition. The projection of the helical curve may be intersected under this condition. It is not easy to label these markers. In order to resolve this problem, we can use a rough prediction model by extending the proposed models with some ideal assumptions. The nearest projection to the rough predictive position is treated as the correspondent projection. This work will be done in the future.

Second, we convert the projection geometries to the C-arm coordinate system. The geometry parameters are independent of the placement of the calibration phantom. Our estimation method adopted the estimated parameters by the helical algorithm, minimizing the oscillation results from the displacement of the phantom. This procedure can simplify the alignment procedure of the iso-center before the calibration.

Last, as Zeego system has more flexibility than the traditional C-arm, it provides more work positions including Head, Left, Right, Table15,Table30. We choose three representative work position Head, Left, Table30 to perform the experiments. The complex mechanical motion of each work position was designed by the machine. The deviation of the machine may result from the gravity and other mechanical problems. We proposed the movement models by analyzing all the posterior information and considering all the possible factors objectively and systematically. The movement models can integrate the refinement and the prediction of the acquisition geometry. They will simplify the calibration procedure in the clinical operation. The experimental results can prove the validity of these models.

Many data acquisition models had been proposed for the traditional C-arm system, however, few models were proposed on the advanced Zeego system. This work fills this vacancy well.

Due to the lack of the clinical patient acquisition, we only apply the proposed algorithm to the phantom reconstruction. The application and evaluation of the clinical coronary artery reconstruction will be done in the future.

Conclusion

For this chapter, we proposed an extended calibration algorithm based on the classical helical phantom. Then, we consider the C-arm system transformation and model the movement based on the estimated posterior information. In this chapter, we consider three representative work positions: Head side, Left side and Table30. The experiment results indicate that the proposed helical calibration algorithm is more accurate and robust than the original algorithm [23]. The estimation of the C-arm system is reliable. The movement models have an acceptable precision to estimate the acquisition parameters. The models also increase the efficiency of the clinical calibration work. Head position has the least deviation and Left position has minor deviation. They can be modeled as a single mean rigid motion. Table30 position has

Chapter 3 3D Reconstruction by Motion Compensation

Introduction

C-arm based Percutaneous coronary intervention (PCI) is a common procedure for the diagnosis and treatment of coronary heart disease (CHD), in which X-ray coronary angiography is one of the most commonly utilized method to assess CHD and is still considered as the gold standard in clinical procedures [1]. To better assist PCI procedures, many efforts have been devoted to perform 3D tomographic reconstruction of high-contrast coronary arteries using angiograms rotationally acquired on a C-arm. However, because of the vessel motion induced by respiratory and heart beating, accurate 3D reconstruction of coronary vessels still remain a challenge [1]. Poor reconstruction affects the accuracy of quantitative coronary analysis (QCA) measurement, such as the length and diameter of the vessel branch or lesion. Inaccurate measurement may lead to incorrect selection of treatment apparatus and would reduce the effectiveness of the treatment consequently.

ECG-gated and motion compensation are two primary ways to handle vessel motion for 3D reconstruction of coronary arteries.

In ECG-gated reconstruction for a specific selected cardiac phase, a subset of angiograms is selected for reconstruction according to recorded ECG signal. The optimal phase should be the one with least motion phase, such as end-systolic and late-diastolic cardiac phases. In lack of the recorded ECG signal, the optimal phase should be determined. Rasche [2] and Blondel [3] adopt the image content to fine the optimal phase. Rasche reconstruct at each phase in one cardiac cycle and compare the intensity distribution. The reconstruction at the optimal phase should have a narrower intensity distribution because of the lack of blurring. However, this method needs too much computation effort. Another hand, Blondel used a line-integral method to estimate the optimal phase. For each raw projection image, each horizontal line is integrated to create one vertical line. They place all the vertical line of all the projection in the cardiac cycle to compose a picture. The composed picture indicates the variation of the intensity and can be used to estimate the optimal phase. This method is easy to perform and needs less computation effort.

Nearest-neighbor (NN) [4] and finite-window gating [5] are two commonly used strategies for selecting sinograms. The nearest neighbor gating [4] was the first proposed approach that used one image at the desired phase for each cycle. Next, the researchers used the finite gating window [5] around the desired phase. The shape of this window can be a rectangular or cosine-squared or a half-bell. If we use the standard filtered back-projection algorithm, the center images of each gating window contribute most and the other projections at the edges of the window contribute according to the weight from the window definition. Multiple images around the desired phase will result in better-gated reconstructions. The commonly used reconstruction method, such as Feldkamp's (FDK) method [6], is quite sensitive to the angular coverage of the angiorams. Small gating window causes under-sampling problem and the reconstruction will have severe streak artifacts. Increasing gating window can reduce the streak artifacts but will introduce more motion artifacts. Recent years, prior knowledge regularization iterative reconstruction has been proposed to handle this problem [7][8][9][10][11][12][13][14][15][16]. By imposing proper prior regularization can suppress the artifacts significantly even using quite few angiograms. These methods have shown promising results.

Motion compensation is a commonly used technique for 3D coronary artery reconstruction. Motion compensated reconstruction has two main categories: 3D motion compensated reconstruction and 2D motion compensated reconstruction.3D motion compensation works in the 3D image space, while 2D motion compensation works in the 2D projection space.

Blondel [17] and Tang [18] proposed the 3D motion compensation based on the 3D motion vector field. These approaches use the 3D centerline information of the cardiac vessel throughout multiple cardiac phases to determine the motion-vector. In theory, all the projection images can be used for reconstruction after motion correction. This strategy can significantly increase the overall signal-to-noise ratio while reducing artifacts caused by gating. However, the accurate estimation of the motion-vector field throughout different cardiac phases remains challenging and 3D motion estimation is also a strong ill-posed problem with high computation cost, Zeng [19] and Hansis [20] proposed the motion periodicity or regularization assumptions to tackle this problem.

In principle, fully 3D motion compensation can compensate for arbitrary motion, and more powerful than 2D compensation. However, the object of interest is spatially sparse and the motion inside the gating window is limited and smooth. Therefore, the 2D motion compensation also can perform well. There are three main approaches of 2D motion compensation reconstruction: landmark-based, model-based, and projection-based. Schafer [21] and Perrenot [22] proposed the landmark-based registration for the 3D reconstruction of stents. Two radiopaque marker balls on each side of the stent are used for tracking the stent motion in the projections during the cardiac cycle. The projections are then transformed accordingly to compensate for the motion. Lebois [23] proposed a model-based learning method to realize 2D motion compensation. They use an extensive projection dataset to learn the model. The residual motion of the actual projection image is corrected according to the model. Hansis [24] proposed the projection-based 2D motion compensation based on vessel centerline segmentation. The forward-projection of the 3D centerline is compared with the actual detected 2D centerlines to estimate the transformation. Then, the residual motion can be corrected by warping the actual projection to the reference-motion state. However, the segmentation of vasculature centerlines in the acquired projection data is also a difficult problem. Recently, Schwemmer proposed another projection-based 2D motion compensation reconstruction method [25]. The method used the normalized cross-correlation (NCC) [26] as the similarity metric of the registration and the modified FDK reconstruction algorithm [27] to realize the compensated reconstruction.

In this chapter, we propose a new projection-based 2D motion compensation reconstruction method. We correct the residual motion by a 2D-2D registration based on mutual information, avoiding the complex vessel centerline segmentation.

In section 2, we introduce some classical methods. In section 3, we describe each step of our proposed method. In section 4, we show the experiment results of our method. In section 5 and 6, we make some discussion and conclusion.

Classical reconstruction method

3D Motion Compensated Reconstruction

In this section, we introduce the 3D motion compensated reconstruction by Blondel et al. [17], which is based on the 4D motion vector field. They proposed a dynamic tomographic algorithm that compensated for the coronary artery motion along the cardiac contraction, by taking advantage of 4D motion field. The algorithm used all available images, independently from the cardiac phase at which they were acquired. We introduce this 3D motion compensation reconstruction briefly in the following.

a. 4D parametric motion field

A 4D motion field was pre-computed before the 3D motion compensation reconstruction. First, they define a reference time in the cardiac cycle time, then, the reconstruction of the coronary artery centerlines are performed at a particular cardiac cycle time t. The 3D reconstruction process typically results in a 3D centerline model. Then, they compute the 4D motion of the coronary arteries from the 3D centerline model.

The motion field was parameterized by a 4D tensor product of B-splines, which was adapted to cardiac motion fitting. The position of 3D point 𝑋 = (𝑥, 𝑦, 𝑧) at time t in 4D B-solid motion is given by the relationship: The motion model is then fitted to a given specific dataset using a large scale optimization process. The optimal motion maximizes an energy function combining an external energy, evaluating the superimposition of projected deformed 3D centerline points with vessels in the angiograms and an internal regularizing energy.

Φ(𝑝, 𝑋, 𝑡) = 𝑋 + ∑ 𝐵 𝑖 (𝑥) .∑ 𝐵 𝑗 (𝑦) .∑ 𝐵 𝑘 (𝑧)

The criterion we optimize is:

Ψ(𝑝) = ∑ ∑ 𝑟 𝑛 𝑋∈𝜒 𝑛∈𝑁 (𝑚 𝑛 (Φ(𝑝, 𝑋, 𝑡 𝑛 ))) + 𝛼Reg 𝑅 3 ×𝑅 Φ(𝑝, 𝑋, 𝑡 𝑛 ) (3.2)
𝜒 is the set of points describing the 3D centerline model, 𝑁 is the set of images from which the motion is estimated, 𝑚 𝑛 is the projection matrix associated with frame n, 𝑟 𝑛 is the multi-scale vessel detector, Reg 𝑅 3 ×𝑅 Φ(𝑝, 𝑋, 𝑡 𝑛 ) is the regularity measure on Φ in space and time.

b. Motion compensated 3D reconstruction

Taking advantage of the 4D parametric motion field, they designed a dynamic reconstruction algorithm. X is the position of a physical 3D point, t is the acquisition time, 0 is considered as the reference time, at which we want to reconstruct the linear attenuation map. 𝜇 𝑅 3 × 𝑅 → 𝑅 is the linear attenuation of any physical 3D point, varying along time. The motion application Φ 𝑡 𝑅 3 × 𝑅 → 𝑅 3 represents the position at a given time of the 3D point that was in a given position at the reference time.

Supposing the linear attenuation remains constant along the acquisition time, with neglecting the contrast agent propagation and diffusion affects. It is equivalent to:

𝜇(Φ 𝑡 (𝑋), 𝑡) = 𝜇(Φ 0 (𝑋), 0) = 𝜇(𝑋, 0) (3.3)
The discrete problem to be solved can be stated as:

𝑅 Φ • 𝜇 = 𝑑 (3.4)
𝑅 Φ is the matrix associated with the projection operator in motion. 𝑑 is the discrete subtracted sinogram data. The coefficient 𝑅 𝑗 𝑖,𝑘 Φ in matrix 𝑅 Φ be the contribution of the voxel 𝑘 to the pixel value 𝑗 𝑖 , belonging to frame 𝑖. Denoting 𝑆 𝑗 𝑖 the solid angle with vertices the corners of pixel 𝑗 𝑖 and the X-ray source position S(𝑡 𝑖 ). For the dynamic case, they take the motion field into account by replacing the voxel cube 𝐶 𝑘 by its image under the 3D motion field Φ 𝑡 𝑖 :

𝑅 𝑗 𝑖,𝑘 Φ = 𝑣𝑜𝑙 .𝑆 𝑗 𝑖 ∩ Φ 𝑡 𝑖 (𝐶 𝑘 )/ (3.5)
Let 𝑐 𝑘 be the center of the voxel cube 𝐶 𝑘 , 𝑚 𝑖 is the projection matrix associated with frame i, 𝑅 𝑗 𝑖,𝑘 Φ has a practical scheme: The motion only impacts the projection operator matrix 𝑅 Φ , many reconstruction techniques can be appropriate for this compensated 3D reconstruction.

2D Motion Compensated Reconstruction

We introduce the projection-based 2D motion compensation reconstruction method from rotational X-ray angiography by Hansis [24] shortly. We show the workflow of the projection-based 2D motion compensated reconstruction method in Fig. 3.1. First, we generate the initial ECG gated reconstruction of the coronary arteries for the desired cardiac phase. Then, we calculate the maximum intensity forward projection from the initial ECG gated reconstruction at each acquisition angle. Next, we correct the residual motion of each current projection image by the registration process. The registration process is performed between each forward projection and the corresponding pre-processed current projection. Last, we generate the compensated reconstruction by the motion corrected projection images. This procedure can be repeated for many cycles to obtain a refinement result.

a. Initial reconstruction and forward projection

An initial ECG-gated reference reconstruction is calculated [27] for the desired cardiac phase, using a suitable gating window. The maximum intensity forward projection 𝑃 𝑖 of the initial reconstruction is calculated according to the projection geometry of 𝑃 𝑖 at each acquisition angle.

b. Pre-processing of the forward and current projection

The vessel centerlines are detected in both the forward projection and the current projection by the proposed centerline filtering algorithm. The algorithm includes the multi-scale vessel filter, directional maxima detection, and hysteresis thresholding. The resulting centerline map is created after these three steps and it has a value at each centerline pixel 𝑝 ⃗ and zero everywhere else.

c. Registration

The current detected projection centerline 𝐶 𝐹 is used as model points 𝑀 = *𝑚 ⃗⃗⃗ 𝑖 + 𝑖=1 𝑁 𝑚 and the detected forward projection centerline 𝐶 𝑃 as data points 𝐷 = {𝑑 ⃗ 𝑖 } 𝑖=1 𝑁 𝑑 . 𝑁 𝑚 and 𝑁 𝑑 are the numbers of model points and data points respectively.

The Levenberg-Marquardt iterative closest point algorithm (LMICP) [28] was adopted to match the corresponding centerlines from 𝐶 𝑃 to 𝐶 𝐹 . The transformation parameters 𝑎 ⃗ of a transformation 𝑇 𝑎 ⃗⃗ can be found such that an error function 𝐸(𝑎 ⃗) is minimal. The error function is defined as:

𝐸(𝑎 ⃗) = ∑ 𝐸 𝑖 (𝑎 ⃗) 𝑁 𝑑 𝑖=1 = ∑ (𝑅 ̃𝑀 .𝑇 𝑎 ⃗⃗ (𝑑 ⃗ 𝑖 )/) 2 𝑁 𝑑 𝑖=1 (3.7) 𝑅 𝑀 (𝑝 ⃗) is the Euclidean distance of a point 𝑝 ⃗ to its closest point in 𝑀. 𝑅 ̃𝑀(𝑝 ⃗) is 𝑅 𝑀 (𝑝 ⃗) clipped to a threshold 𝑅 𝑚𝑎𝑥 : 𝑅 ̃𝑀(𝑝 ⃗) = { 𝑅 𝑀 (𝑝 ⃗), 𝑅 𝑀 (𝑝 ⃗) ≤ R 𝑚𝑎𝑥 𝑅 𝑚𝑎𝑥 , 𝑅 𝑀 (𝑝 ⃗) > R 𝑚𝑎𝑥 (3.8)
The well-known Levenberg-Marquard algorithm [28] was adopted to realize the minimization of (3.7). The LMICP algorithm requires the calculation of a Jacobian matrix 𝐽 ∈ 𝑅 𝑁 𝑑 ×𝑁 𝑎 . 𝑇 𝑎 ⃗⃗ is the bilinear transformation. The entries of 𝐽 are given by:

𝐽 𝑖𝑗 = 𝜕𝐸 𝑖 𝜕𝑎 𝑗 = [∇ ⃗ ⃗⃗ 𝑎 𝐸 𝑖 (𝑎 ⃗)] 𝑗 ⁄ (3.9)
𝑎 𝑗 are the elements of 𝑎 ⃗ and ∇ ⃗ ⃗⃗ 𝑎 is the gradient with respect to 𝑎 ⃗. ∇ ⃗ ⃗⃗ 𝑎 𝐸 𝑖 (𝑎 ⃗) can be write:

∇ ⃗ ⃗⃗ 𝑎 𝐸 𝑖 (𝑎 ⃗) = 2𝑅 ̃𝑀 .𝑇 𝑎 ⃗⃗ (𝑑 ⃗ 𝑖 )/ (∇ ⃗ ⃗⃗ 𝑎 𝑇 𝑎 ⃗⃗ (𝑑 ⃗ 𝑖 )) • (∇ ⃗ ⃗⃗ 𝑅 ̃𝑀 .𝑇 𝑎 ⃗⃗ (𝑑 ⃗ 𝑖 )/) (3.10)

d. Compensated reconstruction

The current projection data must be transformed according to the optimum transformation parameters 𝑎 ⃗ determined by the registration process. The resulting transformation 𝑇 𝑎 ⃗⃗ ′ are applied to generate the motion corrected projection data 𝑃 ̃. The value at point 𝑝 ⃗ of the corrected projection can be found as 𝑃 ̃(𝑝 ⃗) = 𝑃(𝑇 𝑎 ⃗⃗ ′ (𝑝 ⃗)).

Bilinear interpolation is taken in the current projection data for sampling point 𝑇 𝑎 ⃗⃗ ′ (𝑝 ⃗).

The corrected projection data will be the input of the compensated reconstruction. The motion compensation can be applied several times to achieve better results. After one motion compensation cycle, the improved reconstruction will facilitate the centerline determination in forward and current projections (a more accurate mating).

Proposed reconstruction method

In this chapter, we propose a new projection-based 2D motion compensation reconstruction method. The workflow can be seen in Fig. 3.2. First, we use a simple and automatic preprocessing method to extract the coronary arteries. Second, the intial reconstruction and motion compensated reconstruction are realized by the Maximum a Posteriori (MAP) iterative reconstruction algorithm [7], which is more powerful to deal with the ill-posed problem. Third, we adopt the simplified distance driven projector 𝑃 𝑆𝐷𝐷 [29] to generate the forward projection. Last, we adopt the mutual information (MI) [30] and rigidity penalty (RP) [31] as the cost function of the registration. This avoids the complex and difficult vessel centerline extraction process. We use the advanced adaptive stochastic gradient descent (ASGD) algorithm [32] to realize the optimization. The algorithm makes the registration process more fast and accurate. Therefore, the proposed reconstruction method is fully automatic, fast, and easy to implement. We describe each step in this section.

Preprocessing of the Current Projection

For the original acquired projection images, we perform a set of preprocessing steps to segment the vessel structure. We denote the original projection as 𝑝 𝑜𝑟𝑖 . We show the operations as follows: c. top-hat filtering morphological operator [33] with the neighborhood 10×10 to enhance the vessel structure: 𝑝 3 = 𝑡𝑜𝑝𝑎𝑡(𝑝 2 ). d. a threshold operator and then a dilation morphology operators [33] on 𝑝 3 to build the vessel mask 𝑝 𝑚𝑎𝑠𝑘 . e. We adopt 𝑝 𝑚𝑎𝑠𝑘 to estimate the background 𝑝 𝑏𝑔 by the joint interpolation algorithm [34]. f. the segmented projection 𝑝 𝑏𝑔𝑟 = 𝑙𝑛𝑝 𝑏𝑔 -𝑙𝑛𝑝 𝑜𝑟𝑖

Initial Reconstruction

We perform the advanced MAP iterative reconstruction method based on 𝐿 0 prior [7] to realize the initial and compensated reconstruction, see 3.3.2.2. This algorithm has a better performance to deal with the ill-posed inverse problem.

Reconstructed phase

Fig. 3.3 is the illustration of the ECG signal and the gating window. We suppose the patient has a regular cardiac motion. Let 𝑡 ∈ ,0,1be the reconstructed phase, which is the relative distance percentage between two consecutive R-peaks in the ECG signal. The parameterω ∈ ,0,1controls the width of the gating window. T is the cardiac cycle. The projection phases 𝜏 𝑖 ∈ (𝑡 -ω/2, 𝑡 + ω/2) will be used in the each cardiac cycle.

We perform the line-integral method [3] to visualize the cardiac motion and corresponding phases in the case of lacking the ECG record: a. For each original projection 𝑝 𝑜𝑟𝑖 , we segment the vessel as a binary image 𝑝 𝑜𝑟𝑖 𝑏 . We sum up pixels along each horizontal line of each binary image to generate one vertical line.

b. All the vertical line components from the segmented projections were placed together. The resultant image can illustrate the shifting of contrast during the cardiac motion. c. Some morphology operations can be applied to make the resultant image clear. d. The curve is almost periodic and the rest phase with less motion can be determined. The interval between two R-peaks is the period of the cardiac cycle T.

3D Reconstruction algorithm

We describe the reconstruction algorithm briefly. The problem can be formulated as [7]:

A f + b = Y (3.11)
A is the projection operator matrix referred as system matrix. f is the vector of the 3D unknown image volume, b is the acquisition noise, Y is the vector of the 2D projection data.

According to Bayesian theory, the estimation of f can be obtained by the objective function:

𝑓 ̃= 𝑎𝑟𝑔 𝑚𝑖𝑛 𝑓>0 *-𝐿(𝑌|𝑓) + λ𝑅(𝑓)+ (3.12)
𝐿(𝑌|𝑓) is the log-likelihood function reflecting the statistical features of the projection data. Generally, the statistical feature of the observed projection data can be modeled as the Gasussian distribution. Therefore, the projection data is supposed to be stationary during the acquisition time. 𝑅(𝑓) is the image prior information. λ is a positive weight parameter. 𝑓 has a positive hypothesis. I is the number of the measured projection images, the ideal projection data 𝑔:

𝑔 𝑖 = ,𝐴𝑓-𝑖 𝑔 = ,𝑔 1 , 𝑔 2 , … 𝑔 𝐼 - 𝑌 = ,𝑦 1 , 𝑦 2 , … 𝑦 𝐼 - (3.13)
𝜎 𝑖 2 is the variance, 𝐿(𝑌|𝑓) can be deduced by the probability density function of 𝑔 𝑖 as follows:

𝑝(𝑔 𝑖 ) = .1 √2 𝜎 𝑖 2 ⁄ / 𝑒𝑥𝑝(-(𝑦 𝑖 -𝑔 𝑖 ) 2 2𝜎 𝑖 2 ⁄ ) 𝑙 𝑖 (𝑔 𝑖 ) = -𝑙𝑛 𝑝(𝑔 𝑖 ) = (𝑦 𝑖 -𝑔 𝑖 ) 2 2𝜎 𝑖 2 ⁄ + 𝑙𝑛√2 𝜎 𝑖 2 -𝐿(𝑌|𝑓) = ∑ 𝑙 𝑖 (𝑔 𝑖 ) = 𝐼 𝑖=1 ∑ ((𝑦 𝑖 -𝑔 𝑖 ) 2 2𝜎 𝑖 2 ⁄ + 𝑙𝑛√2 𝜎 𝑖 2 ) 𝐼 𝑖=1 (3.14)
We often use the Gibbs prior to model 𝑅(𝑓):

𝑓 = [𝑓 1 , 𝑓 2 , … 𝑓 𝐽 ] 𝑅(𝑓) = ∑ ∑ 𝜔 𝑗𝑘 𝑘∈𝑁 𝑗 𝐽 𝑗=1 (𝑓 𝑗 -𝑓 𝑘 ) (3.15)
𝐽 is the length of the vector 𝑓. 𝑁 𝑗 is the neighborhood of voxel j, 𝜔 𝑗𝑘 is a positive weight parameter, is the symmetric potential function which penalize the pairwise differences between voxels.

We have the definition of 𝐿 0 prior [7]:

𝐿 0 (𝑡) = ‖𝑡‖ 0 (3.16)
‖𝑡‖ 0 means the number of the non-zero elements. However, the optimization of 𝐿 0 prior is a NP hard problem. In this algorithm, we use the following non-convex function to approach to 𝐿 0 prior [7]:

𝐿 0 (𝑡, 𝜌) = |𝑡| (|𝑡| + 𝜌) ⁄ (3.17)
We decrease the value of 𝜌 gradually according to the iteration. The Separable Paraboloidal Surrogate (SPS) algorithm [START_REF] Fessler | A paraboloidal surrogates algorithm for convergent penalized-likelihood emission image reconstruction [C[END_REF] was adopted to realize the optimization. We update the intensity of each voxel simultaneously as [7]:

𝑉(𝐴 𝑇 𝐴) = (∑ 𝑎 𝑖1 2 𝐼 𝑖=1 … ∑ 𝑎 𝑖𝐽 2 𝐼 𝑖=1 ) 𝑇 . 𝜕𝑅(𝑓 𝑗 ) 𝜕𝑓 𝑗 = kN j   𝜔 𝑘𝑗 ̇(𝑓 𝑗 -𝑓 𝑘 ) 𝑓 𝑘+1 = 𝑓 𝑘 + 𝐴 𝑇 (Y-𝐴𝑓 𝑘 ) 𝑉(𝐴 𝑇 𝐴) + 𝛽 𝜕𝑅(𝑓) 𝜕𝑓 𝑓=𝑓 𝑘 ( 𝑓 𝑘+1 ) 𝑗 = 𝑚𝑎𝑥((𝑓 𝑘+1 ) 𝑗 , 0) (3.18)
( 𝑓 𝑘+1 ) 𝑗 represents the j th voxel of the (k+1) th iteration.

When we reconstruct the phase 𝑡 with the gating window width ω for the initial reconstruction, Y is composed by the projection data 𝑝 𝑏𝑔𝑟 with the phases 𝜏 𝑖 ∈ (𝑡 -ω/2, 𝑡 + ω/2) in each cardiac cycle. We denote the reconstruction at phase 𝑡 as 𝑓 𝑡 .

Simplification of the system matrix

We adopt the simplified distance driven projector 𝑃 𝑆𝐷𝐷 [29] to generate the system matrix 𝐴 according to the acquisition geometries. This projector was proved to be faster and comparable accurate than the classical Ray-driven (RD) [36] and Separable Footprint (SF-TT) [37]. Xie [7] proposed to generate a 3D binary mask M of the image volume 𝑓 to reduce the storing space of 𝐴.

As we described in chapter 2, 𝑎,𝑢 𝑘 , 𝑣 𝑙 ; 𝒏is the coefficient of the system matrix 𝐴, 𝒏 is the voxel of the image volume 𝑓, (𝑢 𝑘 , 𝑣 𝑙 ) is the pixel position in the detector. We combine 𝑃 𝑆𝐷𝐷 with the 3D mask M to reduce the computation cost. We have:

{

𝑎,𝑢 𝑘 , 𝑣 𝑙 ; 𝒏-𝑏𝑦 𝑃 𝑆𝐷𝐷 𝑖𝑓 𝑀(𝒏) = 1 𝑎,𝑢 𝑘 , 𝑣 𝑙 ; 𝒏-= 0 𝑖𝑓 𝑀(𝒏) = 0 (3.19)

Maximum Intensity Forward Projection

a. Thresholding Generally, vascular structure has the high-contrast in the 3D image. Therefore, we can use a simple thresholding operation to remove the background structures in the initial 3D cardiac vasculature reconstruction. 𝑇 𝑟 ∈ ,0,1is the thresholding parameter, L is the largest voxel value of the initial reconstruction.The initial reconstruction after thresholding is denoted as 𝑓 𝑡 𝑇 . We have the following operation:

𝑓 𝑡 𝑇 (𝒏) = { 𝑓 𝑡 (𝒏) 𝑖𝑓 𝑓 𝑡 (𝒏) ≥ 𝑇 𝑟 * 𝐿 0 𝑖𝑓 𝑓 𝑡 (𝒏) < 𝑇 𝑟 * 𝐿 (3.20)

b. Forward projection

We also use the projector 𝑃 𝑆𝐷𝐷 [29] to generate the forward projection 𝑝 𝑓𝑤𝑝 of 𝑓 𝑡 𝑇 at each acquisition angle i. We have the expression as:

𝑝 𝑓𝑤𝑝𝑖 = 𝐴 𝑖 * 𝑓 𝑡 𝑇 (3.21)
𝐴 𝑖 corresponds to angle i and is the part of system matrix 𝐴.

Image Registration

Image registration is the significant step in the proposed 2D motion compensated reconstruction method. The registration of a fixed image (reference image) 𝑝 𝑓𝑤𝑝 Ω 𝑓 ∈ ℝ 2 → ℝ to the moving image (floating image) 𝑝 𝑏𝑔𝑟 Ω 𝑚 ∈ ℝ 2 → ℝ can be formulated as an optimization problem:

𝜇̂= 𝑎𝑟𝑔𝑚𝑖𝑛 𝜇 𝐶(𝑇 𝜇 ; 𝑝 𝑓𝑤𝑝 , 𝑝 𝑏𝑔𝑟 , Ω 𝑓 ) (3.22)
𝑇 𝜇 is a parameterized coordinate transformation. Ω 𝑓 → Ω 𝑚 often represents the overlapped domain. The transformation 𝑇 𝜇 is modeled by the parameters 𝜇. C is the cost function and the solution 𝜇̂ is the parameter vector that minimizes this cost function.

We plot the registration procedure in Fig. 3.4. The cost function C introduces the steps of sampler, transformation and interpolation. Registration means to resolve the optimization problem, multi-resolution scheme is often adopted. We describe the concrete algorithm of each step in the following.

Computing of cost function

Mutual information is a basic concept from information theory. It can measure the statistical dependence on information redundancy between two random variables or the amount of information that one variable contains about the other. The iconic registration based on mutual information (MI) was first proposed by Maes et al. [30] The registration criterion of mutual information is that if the MI of the two images in the overlapped domain is maximal, the two images are geometrically aligned. This iconic matching criterion is very powerful for our motion compensated reconstruction application [38]. We can apply it to the projection image without prior segmentation, feature extraction or other complex preprocessing.

a. Definition of mutual information

There are two random variables 𝐴 and B, 𝜌 𝐴 (𝑎), 𝜌 𝐵 (𝑏) are their marginal probability distributions, 𝜌 𝐴𝐵 (𝑎, 𝑏) is their joint probability distribution. If A and B are statistically independent, we have 𝜌 𝐴𝐵 (𝑎, 𝑏) = 𝜌 𝐴 (𝑎) • 𝜌 𝐵 (𝑏). If they are related by a one-to-one mapping (maximally dependent), 𝑇 is the transformation, we have 𝜌 𝐴 (𝑎) = 𝜌 𝐵 (𝑇(𝑎)) = 𝜌 𝐴𝐵 (𝑎, 𝑇(𝑎)). The mutual information 𝐼(𝐴, 𝐵) measures the degree of dependence of A and B by measuring the distance between the joint distribution 𝜌 𝐴𝐵 (𝑎, 𝑏) and the independence distribution 𝜌 𝐴 (𝑎) • 𝜌 𝐵 (𝑏) [30].

The definition of mutual information (MI) by means of the Kullback-Leibler measure is:

𝐼(𝐴, 𝐵) = , ab  𝜌 𝐴𝐵 (𝑎, 𝑏)𝑙𝑜𝑔 𝜌 𝐴𝐵 (𝑎,𝑏) 𝜌 𝐴 (𝑎)•𝜌 𝐵 (𝑏) (3.23)
Another definition of the mutual information is related to the entropy of images:

𝐼(𝐴, 𝐵) = 𝐻(𝐴) + 𝐻(𝐵) -𝐻(𝐴, 𝐵) = 𝐻(𝐴) -𝐻(𝐴|𝐵) = 𝐻(𝐵) -𝐻(𝐵|𝐴) (3.24)
𝐴, 𝐵 represent the image intensity. 𝐻(𝐴) and 𝐻(𝐵) being the Shannon entropy of A and B respectively, 𝐻(𝐴, 𝐵) is their joint entropy, 𝐻(𝐴|𝐵) and 𝐻(𝐵|𝐴) the conditional entropy of A given B and of B given A respectively:

𝐻(𝐴) = - a  𝜌 𝐴 (𝑎)𝑙𝑜𝑔𝜌 𝐴 (𝑎) 𝐻(𝐴, 𝐵) = - , ab  𝜌 𝐴𝐵 (𝑎, 𝑏)𝑙𝑜𝑔𝜌 𝐴𝐵 (𝑎, 𝑏) 𝐻(𝐴|𝐵) = - , ab  𝜌 𝐴 (𝑎)𝑙𝑜𝑔𝜌 𝐴|𝐵 (𝑎|𝑏) (3.25) b. Criterion
Let 𝐼 𝐹 (𝑥) denote the image intensity in the fixed image (reference image) 𝑝 𝑓𝑤𝑝 at pixel position 𝑥, 𝐼 𝑀 (𝑇 𝜇 (𝑥)) is the intensity at the transformed position 𝑇 𝜇 (𝑥) in the moving image (floating image) 𝑝 𝑏𝑔𝑟 . 𝑇 𝜇 (𝑥) is the parametric transformation. The image intensity values of two images to be registered are random variables. The registration criterion states that if the MI is maximal, the two images are geometrically aligned and the current 𝜇̂ is the optimal transformation parameters.

In the general application, the two images are only partially overlapped, the overlapped domain will change when 𝜇 is varied. Estimations of the marginal 𝜌 𝐹,𝜇 (𝐼 𝐹 ), 𝜌 𝑀,𝜇 (𝐼 𝑀 ) and joint distribution 𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) can be obtained by normalization of the respective histogram 𝜇 (•). We have the following formulations:

𝜌 𝐹,𝜇 (𝐼 𝐹 ) = 𝜇 (𝐼 𝐹 ) ∑ 𝜇 (𝐼 𝐹 ) 𝐼 𝐹 ⁄ 𝜌 𝑀,𝜇 (𝐼 𝑀 ) = 𝜇 (𝐼 𝑀 ) ∑ 𝜇 (𝐼 𝑀 ) 𝐼 𝑀 ⁄ 𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) = 𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) ∑ 𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) 𝐼 𝐹 ,𝐼 𝑀 ⁄ (3.26)
𝜇 (𝐼 𝐹 ), 𝜇 (𝐼 𝑀 ) are the image intensity histogram of the fixed image 𝐼 𝐹 and the moving image 𝐼 𝑀 . 𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) is the joint image intensity histogram. All the histograms can be computed by binning the image intensity values 𝐼 𝐹 (𝑥), 𝐼 𝑀 (𝑇 𝜇 (𝑥)), (𝐼 𝐹 (𝑥), 𝐼 𝑀 (𝑇 𝜇 (𝑥))) for the current overlapping domain 𝑥 ∈ Ω 𝜇  Ω 𝑓 . In order to do this efficiently, the moving and fixed image intensities are binned linearly rescaled to the range ,0, 𝑛 𝑀 -1and ,0, 𝑛 𝐹 -1respectively, 𝑛 𝑀 × 𝑛 𝐹 being the total number of bins in the joint histogram. Typically, we use 𝑛 𝑀 × 𝑛 𝐹 =256. However, the histogram method may not accurate to update the joint probability 𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) [39].

We adopt the B-spline Parzen window, which was first proposed by Viola [39] and Thévenaz [40] as a classical method to calculate the joint probabilities. We have the following formulation of the Parzen window joint histogram: 𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) = 𝜂 𝐹 and 𝜂 𝑀 represent the fixed and moving B-spline Parzen window kernels [39] used to distribute an intensity over the neighboring bins. The scaling constants 𝜎 𝐹 and 𝜎 𝑀 equal to the intensity bin widths. We can have the marginal discrete probabilities from the joint discrete probability [39]:

𝜌 𝐹,𝜇 (𝐼 𝐹 ) = M I  𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) 𝜌 𝑀,𝜇 (𝐼 𝑀 ) = I F  𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 ) (3.28)
The similarity function 𝐼(𝐼 𝐹 , 𝐼 𝑀 ) is then expressed by the expression:

𝐼(𝐼 𝐹 , 𝐼 𝑀 ) = 𝐻(𝐼 𝐹 ) + 𝐻(𝐼 𝑀 ) -𝐻(𝐼 𝐹 , 𝐼 𝑀 ) = -∑ 𝜌 𝐹,
The optimal registration parameter 𝜇̂ is found from:

𝜇̂= 𝑎𝑟𝑔 𝑚𝑖𝑛 𝜇 -𝐼(𝑇 𝜇 ; 𝐼 𝐹 , 𝐼 𝑀 , Ω 𝜇 )

The cost function to be minimized is equals to the mutual information similarity metric. However, the registration problem is ill-posed for the non-rigid transformation, a regularization term (penalty term) Ψ is often introduced to constrain 𝑇 𝜇 . Indeed, the cost function may have multiple local minima. The local minimum that selected as the solution 𝜇̂ depends on the optimization algorithm and the initial alignment of the images. This regularization term Ψ(𝑇 𝜇 ) is added to the cost function to penalize the undesirable deformation and consequently to reduce the number of local minima.

The registration formulation should be rewritten as follows:

𝜇̂= 𝑎𝑟𝑔 𝑚𝑖𝑛 𝜇 -𝐼(𝑇 𝜇 ; 𝐼 𝐹 , 𝐼 𝑀 , Ω 𝜇 ) + 𝛾Ψ(𝑇 𝜇 )

(3.31)
𝛾 is a constant that weights the similarity against regularity.

c. Regularization

We adopt the rigidity penalty proposed by Staring [31] in our non-rigid registration. The projection images that are to be registered often contain both structures that deform and ones that remain rigid. This regularization adds a local rigidity penalty term to penalize the deformation of rigid objects. Three conditions were combined in the penalty term to model the local rigid transformation.

For a displacement field 𝑇 𝜇 (𝑥) to be rigid, it must hold that: In order to distinguish between rigid and non-rigid tissue, the total penalty term is weighted by a rigidity coefficient 𝑐(𝑥) ∈ ,0,1of the tissue type at position 𝑥. The weights 𝑐 𝐴𝐶 , 𝑐 𝑂𝐶 and 𝑐 𝑃𝐶 determine the relative strength of each of the three terms. The rigidity coefficient 𝑐(𝑥) is set to 0 for pixels 𝑝 in completely non-rigid tissue, and 1 for rigid tissue, for other tissue types a value of 𝑐(𝑥) is chosen between 0 and 1. The rigidity coefficient image only has to be defined on the moving image, as we calculate the cost function from fixed image to moving image. For our application, unconstrained non-rigid registration (local deformation) severely compressed the vasculature, which can be avoided with the use of Ψ(𝑇 𝜇 ).

𝑇

Multi-resolution strategy

Some researches show that the subsampling of the high-resolution images can be applied without deteriorating the registration robustness [41][42]. Important speed-ups can be realized by using a multi-resolution optimization strategy, starting with a coarsely sampled image for efficiency and increasing the resolution during the optimization proceeds for accuracy. We adopt the Gaussian pyramid [41][42], which means that the images are applied smoothing by Gaussian filtering and down-sampling.

Approximation by subsampling

When we calculate the cost function based on mutual information, it is not necessary to loop over all the pixels 𝑥 ∈ Ω 𝜇 of the fixed image, a subset may suffice. In order to reduce the computation time, many subsampling stochastic gradient descent optimization method adopt the random subsampling, which selects a subset of pixels from the fixed image randomly for each iteration [43]. Random subsampling will reduce the computational cost substantially, without compromising registration accuracy.

Geometric Transformation

We apply firstly the affine transformation to model the rigid deformation field, secondly, the B-spline transformation to model the local deformation field on non-rigid deformation.

a. Affine transformation

The affine transformation is defined as:

𝑇 𝜇 (𝑥) = 𝐸(𝑥 -𝑐) + 𝑡 + 𝑐 (3.37)
𝑐 is the center of the image. 𝐸 is a matrix that includes the rotation, scaling, shearing.

𝑡 is the translation. The parameters 𝜇 are formed by 𝐸, 𝑡.

b. B-spline transformation

We use a cubic B-spline to model the local deformation [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF]:

𝑇 𝜇 (𝑥) = 𝑥 + xN x k   𝑙 𝑘 𝛽

B-spline Interpolation

In general, when we calculate the mutual information, 𝑇 𝜇 (𝑥) will not locate at a grid point of moving image M. 𝐼 𝑀 (𝑇 𝜇 (𝑥)) must be obtained by the interpolation. We use the 3 order B-spline to do the interpolation [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF], it is also be used to generate the final corrected projection after registration.

Optimization methods

To resolve the registration problem in equations (3.30) and (3.31), we commonly use an iterative optimization to obtain the optimal transformation parameter vector 𝜇.

We have the following formulation:

𝜇 𝑘+1 = 𝜇 𝑘 + 𝑠 𝑘 𝑑 𝑘 , 𝑘 = 0,1,2, … (3.39)
𝑑 𝑘 is the searching direction at iteration 𝑘. 𝑠 𝑘 is the scalar gain factor that controls the step size, this can be determined by a line search or a predefined function of 𝑘.

Classically, Powell's multidimensional direction set method was used combined with the Brent's one-dimensional line search algorithm [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF]. However, the specific cost function of mutual information can lead to an analytic expression of its gradient, according to this, researchers proposed many effective optimization [43]. We adopt the most advanced adaptive stochastic gradient descent (ASGD) algorithm [32] to do the optimization. ASGD algorithm is the development based on the gradient descent (GD) algorithm and the Robbins-Monro (RM). The algorithm needs less parameters and tends to be faster and more robust than GD and RM. We describe the three algorithms as follows: a. Gradient descent (GD) [43]:

𝜇 𝑘+1 = 𝜇 𝑘 -𝑠 𝑘 𝑔(𝜇 𝑘 ) (3.40)
Gradient descent takes steps in the direction of the negative gradient of the cost function, which is the derivative of the cost function evaluated at the current position 𝜇 𝑘 .

The gradient 𝑔(𝜇 𝑘 ) has a general formulation: The search direction 𝑑 𝑘 = -𝑔(𝜇 𝑘 ). More details of 𝑔(𝜇 𝑘 ) can be seen in [43].

𝐶(𝜇) = Φ (
𝑠 𝑘 can be determined by a line search algorithm [43].

b. Robbins-Monro (RM) [43]:

The function 𝑔(𝜇 𝑘 ) is replaced by an approximation 𝑔 ̃𝑘 resulting in the following: The search direction 𝑑 𝑘 = -𝑔 ̃𝑘 is the approximation of the gradient 𝑔(𝜇 𝑘 ), this approximation 𝑔 ̃𝑘 is done by randomly subsampling a subset 𝐵 𝑘 of the fixed image every iteration k. RM is much faster than GD because of the random sampling, while without compromising on accuracy [43]. However, 𝑠, 𝐻, α depend on the specific problem, the set of these parameters complicate the application of this algorithm. c. Adaptive stochastic gradient descent (ASGD) [32]:

𝜇 𝑘+1 = 𝜇 𝑘 -𝑠 𝑘 𝑔 ̃𝑘 𝑠 𝑘 = 𝑠 (𝑘 + 𝐻) 𝛼 ⁄ 𝑠 > 0, 𝐻 ≥ 1, 0≤ α ≤ 1 𝑔 ̃𝑘 =
𝜇 𝑘+1 = 𝜇 𝑘 -𝑠 𝑘 𝑔 ̃𝑘 𝑠 𝑘 = 𝑠(𝑚 𝑘 ) = 𝑠/(𝑚 𝑘 + 𝐴 𝑎 ) 𝑚 𝑘+1 = ,𝑚 𝑘 + (-𝑔 ̃𝑘 𝑇 𝑔 ̃𝑘-1 )-+ ,ξ-+ = 𝑚𝑎𝑥(ξ, 0) (ξ) = 𝑀𝐼𝑁 + 𝑀𝐴𝑋 -𝑀𝐼𝑁 1 -( 𝑀𝐴𝑋 𝑀𝐼𝑁 ⁄ )𝑒 -ξ/ϵ 𝑀𝐴𝑋 > 0, 𝑀𝐼𝑁 < 0, ϵ > 0 (3.43)
For , we define a general sigmoid shape with (0) = 0. The step size 𝑠 𝑘 is adaptive and evaluated at the iteration time 𝑚 𝑘 . 𝜇 0 , 𝑚 0 , 𝑚 1 are user-defined initial conditions. The new introduced parameters are 𝑠, 𝐴 𝑎 , 𝑀𝐴𝑋 , 𝑀𝐼𝑁 , ϵ . 𝐴 𝑎 is user defined, empirically, a reasonable setting is to use 𝐴 𝑎 =20. The ASGD algorithm proposed in [32] δ appears to be unaffected by the choice of the similarity measure. Empirically, we can set δ to be the average pixel size of the images.

Step 3: Generate 𝑁 𝜇 instances of 𝜇 𝑘 according to 𝜇 𝑘 ~𝑁(𝜇 0 , 𝜎 4 2 𝐼 𝑖 ) , 𝐼 𝑖 is the identity matrix.

Compute for each 𝜇 𝑘 the exact cost function derivative 𝑔, the approximated derivate 𝑔 ̃𝑘, and the approximation error 𝜀 𝑘 = 𝑔 -𝑔 ̃𝑘. Note that, to compute 𝑔 ̃𝑘, a new set of pixels 𝐵 𝑘 must be selected for each 𝜇 𝑘 .

Step 4: Compute 𝜎 1 , 𝜎 3

𝜎 1 2 = 1 𝑁 𝜇 1 N n    ‖𝑔(𝜇 𝑛 )‖ 2 /𝑡𝑟(𝐂) 𝜎 3 2 = 1 𝑁 𝜇 1 N n    ‖𝜀(𝜇 𝑛 )‖ 2 /𝑡𝑟(𝐂) (3.46)
Step

κ is a defined constant.

3. Start the optimization defined by (3.43). Convergence is assumed after K iterations:

𝜇̂= 𝜇 𝐾 (3.51)
The ASGD algorithm is more robust than the RM algorithm [32], because of its adaptive step size prediction.

Motion Compensated Reconstruction

After the registration of 𝑝 𝑓𝑤𝑝 and 𝑝 𝑏𝑔𝑟 , the corrected projection data 𝑝 𝑐𝑝𝑟 must be generated according to the obtained optimum transformation 𝑇 𝜇 . The intensity value of the corrected projection at point 𝑝 can be given by 𝑝 𝑐𝑝𝑟 (𝑥) = 𝑝 𝑏𝑔𝑟 (𝑇 𝜇 (𝑥)), 𝑝 𝑏𝑔𝑟 is the preprocessed current projection. We also adopt the 3 order B-spline to do this interpolation. We use the same reconstruction algorithm as initial reconstruction. When we reconstruct at the phase 𝑡 with the gating window width ω 𝑐 for the motion compensated reconstruction, Y is composed by the corrected projection data 𝑝 𝑐𝑝𝑟 with the phases 𝜏 𝑖 ∈ (𝑡 -ω 𝑐 /2, 𝑡 + ω 𝑐 /2) in each cardiac cycle.

Multiple Iterations

The motion compensated reconstruction procedure can be repeated for several cycles to achieve high refinement. After one motion compensation cycle, the quality of the 3D reconstruction has been improved, which generate a better maximum intensity forward and the more accurate transformation parameters.

Experiments and Results

In this section, we will present the results of each step described in 3.3. They are the cardiac phase selection, segmentation, image registration and the motion compensated reconstruction.

Simulation Data

CAVAREV (CArdiac VAsculature Reconstruction EValuation) [46] is a public and open platform for the evaluation of cardiac vasculature reconstruction algorithms from C-arm CT (rotational angiography). The platform proposes the evaluation measures which are independent from the reconstruction algorithm.

It has two simulated dynamic projection datasets based on the 4D XCAT phantom with contrasted coronary arteries which was derived from real patient CT/MR data. In the first dataset, the vasculature undergoes a continuous periodic motion. The second dataset contains aperiodic heart motion by including additional breathing motion.

We assume a breath hold acquisition, and we use the cardiac motion-only dataset for our experiment. The acquisition scenario and geometry calibration was obtained from a standard protocol of a clinical angiographic C-arm system (Artis Zee system). The acquisition time is T=5.3s with N=133 simulated projection images created from an equiangular 200 ° acquisition of a software phantom. The heart rate of the periodic cardiac motion is 80bpm. It means there are totally 7 cardiac cycles, about 19 projection images of one cycle. Phase zero corresponds to an R-peak in the ECG. The projection images have a height of 𝑆 𝑦 = 960 pixels and width of 𝑆 𝑥 = 960 pixels with an isotropic resolution of 0.32 mm. The isocenter-source distance (SAD) and the source-detector distances (SDD) of the imaging system are about 80cm and 120cm, respectively. The projection matrices 𝑃 𝑖 ∈ 𝑅 3×4 , 𝑖 ∈ *1, … 𝑁+ were obtained from an offline calibration of the clinical C-arm system. The simulation projections were generated by the projection matrices 𝑃 𝑖 .We adopted the calibrated geometries 𝑃 𝑖 in the reconstruction. All the reconstructions were calculated in a volume of 196×196×196 voxels, isotropic voxel size is 0.5 mm.

Cardiac motion phase selection

We adopt the method described in 3.3.2.1 to visualize the cardiac motion and corresponding phase, as shown in Fig. 3.5. See Fig. 3.5, we find that there are about 7 motion cycles, which can reflect the cardiac cycles indirectly. We select the phase with the least motion to be the rest phase, 𝑡 𝑟 = 0.4 (end-systolic). The phase with stronger motion is 𝑡 𝑠 = 0.7 . The two are the evaluation phases. The initial reconstruction and the motion compensated reconstruction can be performed with different gating window settings. A relatively small gating window 𝜔 for the initial reconstruction can avoid motion blur. We adopt a larger gating window ω 𝑐 for the compensated reconstruction to increase the image contrast and reduce artifacts.

For the rest phase 𝑡 𝑟 , the gating window of 𝜔 = 0.15 was adopted as the input of the initial reconstruction. We set the gating window widths ω 𝑐 = 0.15,0.3,0.45,0.6 for the compensated reconstruction. For the phase with stronger motion 𝑡 𝑠 , we adopt a narrow gating window of 𝜔 = 0.05 to realize its initial reconstruction. We set ω 𝑐 = 0.05,0.1,0.2,0.4 to realize the compensated reconstruction.

Segmentation and Image Registration

a. Segmentation of the original projection As described in 3.3.1, our approach does not require any complex segmentation of the coronaries and is fully automatic.

b. Image registration

In the motion compensated reconstruction, the parameters of the registration are set as: 𝑇 𝑟 =0.05, the intensity bin widths 𝜎 𝐹 = 𝜎 𝑀 =16, random subsampling of 2048 pixels, the order of B-spline transformation and interpolation is 3. Optimization parameters are 𝐴 𝑎 = 20, δ = 0.6, 𝑚 0 = 𝑚 1 = 0, 𝑁 𝜇 =10, κ = 0.1, K=500. The We designate the method proposed by Schwemmer [25] as method S. It is the latest research of 2D motion compensated reconstruction. The ECG-gated 𝐿 0 prior iterative reconstruction method [7] and method S [25] are taken as two reference methods. We use two iterative cycles in the motion compensated reconstruction.

a. Qualitative evaluation of the motion compensated reconstruction

The qualitative evaluation is carried out visually. We plot the reconstruction results at the two evaluation phases and also use the ray-casting algorithm [START_REF] Levoy | Display of surfaces from volume data[END_REF] to visualize the 3D reconstruction. The volume rendering parameters are identical for each reconstruction result.

We have 3 reconstruction results of each ω 𝑐 : the result by ECG-gated iterative reconstruction [7], the reconstruction by method S [25], the reconstruction by the proposed method. Fig. 3.8 shows the reconstruction at the rest phase 𝑡 𝑟 . Fig. 3.9 shows the reconstruction at phase 𝑡 𝑠 . We observe the reconstructions in Fig. 3.8 and Fig. 3.9. We find that the ECG-gated reconstruction [7] achieves worse result than the two motion compensated reconstructions and the effect is more obvious when we increase the gating window width.

The proposed method has a better result than method S [25] in the artifact level, visibility, details of the vessel structure. The proposed method achieves better quality with wider gating window.

The improvement of image quality at the stronger motion phase 𝑡 𝑠 is also obvious. This indicates that the proposed method can obtain a high 3D quality even at a non-optimal phase. b. Quantitative evaluation of the motion compensated reconstruction Rohkohl [START_REF] Rohkohl | CAVAREV-an open platform for evaluating 3D and 4D cardiac vasculature reconstruction[END_REF] 𝐷𝑠𝑐 is the Dice similarity coefficient defined as [START_REF] Crum W R | Generalized overlap measures for evaluation and validation in medical image analysis[END_REF]:

𝐷𝑠𝑐(𝑓 1 , 𝑓 2 ) = 2 ∑ 𝑓 1 (𝑥) • 𝑓 2 (𝑥) 𝑥 (∑ 𝑓 1 (𝑥) + 𝑓 2 (𝑥) 𝑥 ) ⁄ (3.53)
𝐷𝑠𝑐 ranges from 0, for no spatial overlap, to 1, for a perfect match. 𝑇(𝑓, 𝑎) is a threshold function to generate a binary volume of 𝑓, we have the definition as:

𝑇(𝑓, 𝑎)(𝒏) = { 1 𝑓(𝒏) ≥ 𝑎 0 𝑓(𝒏) < 𝑎 (3.54)
Then, the quality metric 𝑄 3𝐷 (𝑓) searches the image frame whose motion phase fits best to the reconstruction, it can be computed as:

𝑄 3𝐷 (𝑓) = 𝑚𝑎𝑥 𝑖∈*1,…,133+ 𝑄 𝑖 (𝑓) (3.55)
We show the quantitative evaluation metric value 𝑄 3𝐷 (𝑓) of the reconstruction results in Table 3.1. It can be seen that the ECG-gated reconstruction [7] has a lower metric value than the two motion compensated reconstruction. Its metric value decreases with the wider gating window width. The proposed method has a better metric value than method S. The reconstruction of phase 𝑡 𝑠 has the comparable metric values as phase 𝑡 𝑟 .

Compared with other methods in the ranking list of CAVAREV, 𝑄 3𝐷 (𝑓) value of the proposed method is at the top of the list (second place).

Reconstruction from improved calibration

As described in 3.4.1, the simulation dataset was generated by the calibrated projection matrices 𝑃 𝑖 ∈ 𝑅 3×4 , 𝑖 ∈ *1, … ,133+, which were obtained from an offline calibration of the C-arm system [46]. In the above section 3.4.1.3, we adopt 𝑃 𝑖 to evaluate the reconstruction algorithm conveniently. In this section, we come back and try to evaluate the contribution effect of calibration work in the reconstruction by the simulated data. We still use the cardiac motion-only dataset provided by CAVAREV.

We still adopt the proposed 2D motion compensated method to realize the reconstruction. Now, the reconstruction is done from two geometries: the calibrated geometries and the ideal geometries. The geometries can be applied in two steps of the 2D motion compensated reconstruction. One is the steps of the initial and compensated reconstruction. The acquisition geometry will affect the calculation of the system matrix 𝐴, which connects the 2D projection data and the 3D image volume. The other is the step of forward projection. We need to use the geometry at each acquisition angle to generate the forward projection.

We select the reconstruction phase 𝑡 𝑟 = 0.4, 𝜔 = 0.15 for the initial reconstruction, ω 𝑐 = 0.15,0.3,0.45,0.6 for the motion compensated reconstruction. We denote the reconstruction from calibrated geometries as 𝑓 𝐺 𝑐 , the reconstruction with ideal geometries as 𝑓 𝐺 𝑖 . All the reconstruction parameters are the same as in section 3.4.1.3.

a. Qualitative evaluation of the reconstruction

We plot the reconstruction results of 𝑓 𝐺 𝑐 and 𝑓 𝐺 𝑖 in Fig. 3.10. The volume rendering parameters [START_REF] Levoy | Display of surfaces from volume data[END_REF] are identical for each reconstruction.

As displayed in Fig. 3.10, 𝑓 𝐺 𝑖 has the obvious degenerated vessel structure and artifacts and the some parts of the vessel structure is not improved with the wider gating window. Another hand, 𝑓 𝐺 𝑐 has the complete vessel structure and the structure achieves less artifacts with the wider gating window. This shows that the calibration work is essential for the 2D motion compensated reconstruction.

Real Clinical Data

We adopt the real clinical dataset to evaluate the proposed method. The data was acquired on a GE C-arm CT system from Ponchaillou Hospital (Rennes, France). 117 projection images were acquired during a 155 ° rotation of the C-arm gantry.

The distance from source to detector (SDD) is 1050 mm, and the distance from source to iso-center (SAD) is 720 mm. The resolution of the projection image is 512×512, the isotropic pixel size is 0.2875 mm. No patient breathing motion exists over the acquisition and projection images have good contrast in the vessels. All the reconstructions were calculated in a volume of 256×256×256 voxels with isotropic Fig. 3.11. Cardiac motion and phase visualization voxel size 0.5 mm. We adopt the ideal geometries in our reconstruction and ignore the geometry deviation.

Cardiac motion phase selection

We adopt the described method in 3.3.2.1 to visualize the cardiac motion and corresponding phase. As shown in Fig. 3.11, there are 5 available cardiac cycles and 17 projection images of one cardiac cycle T. The two reconstructed phases are: rest phase 𝑡 𝑟 = 0.35 (end-systolic) and the stronger motion phase 𝑡 𝑠 = 0.65 .

For the rest phase 𝑡 𝑟 , we set 𝜔 = 0.2 for initial reconstruction and ω 𝑐 = 0.2,0.4, 0.6,0.8 for motion compensated reconstruction. For the stronger motion phase 𝑡 𝑠 , we set 𝜔 = 0.05 for initial reconstruction and ω 𝑐 = 0.05,0.1,0.2,0.4 for motion compensated reconstruction.

Segmentation and Registration

a. Segmentation of the original projections Figure 3.12 shows one typical example of the preprocessing of the original projection data. The projection image corresponds to the phase 𝑡 𝑟 of the first cardiac cycle. (a) is the original projection, (b) is the binary mask of the vessel structure. (c) is the estimated background, (d) is the extracted vessel structure. We observe that the vessel structure was extracted well from the original projection.

b. Image registration

In the clinical experiment, the parameters of the registration are set as: 𝑇 𝑟 =0.05, 𝜎 𝐹 = 𝜎 𝑀 =16, random subsampling of 2048 pixels, the order of B-spline transformation and interpolation is 3. Optimization parameters are 𝐴 𝑎 = 20, δ = 0.6, 𝑚 0 = 𝑚 1 = 0, 𝑁 𝜇 =10, κ = 0.1, K=500. The number of resolution level is 4, the row and columns of the images have the same sampling factor, the sampling factor

Motion compensated reconstruction evaluation

The ECG-gated 𝐿 0 prior iterative reconstruction method [7] and method S [25] are the reference methods. We apply 2 iteration cycles to realize the motion compensated reconstruction of method S and the proposed method.

a. Qualitative evaluation and comparison

We plot the reconstruction results at the two evaluation cardiac phases by volume rendering [START_REF] Levoy | Display of surfaces from volume data[END_REF] to evaluate qualitatively and compare these 3 methods from real clinical data. Fig. 3.14 shows the reconstruction at the rest phase 𝑡 𝑟 . Fig. 3.15 shows the reconstruction at phase 𝑡 𝑠 . All the volume rendering parameters are identical.

We compare the three reconstruction results from real clinical data with the same ω 𝑐 in Fig. 3.14. It can be seen that the two reconstruction results by motion compensated methods have an obvious improvement than the one by ECG-gated [7]. The proposed method has a better result than the reference method S [25] in the artifact level, visibility, contrast and details of the vessel structure.

We compare the three reconstructions with different ω 𝑐 in Fig. 3.14. The reconstruction by ECG-gated achieves more artifacts with the wider gating window 93 width. The motion compensated reconstruction by the proposed method achieves better quality with wider ω 𝑐 . The proposed motion compensated reconstruction with ω 𝑐 of 0.8 achieves the best result. We compare the three reconstruction results with the same ω 𝑐 in Fig. 3.15. It can be seen that the two reconstruction results by motion compensated methods have an obvious improvement than the one by ECG-gated [7]. The proposed method has a better result than the reference method S [25] in the artifact level, visibility and contrast of the vessel structure.

We compare the three reconstructions from clinical data with different ω 𝑐 in Fig. 3.15. The reconstruction without motion compensation achieves more artifacts with the wider gating window width. The reconstruction by the proposed method achieves the better quality with the wider gating windows. The compensated reconstruction with ω 𝑐 = 0.4 at phase 𝑡 𝑠 achieves an acceptable quality.

Comparing Fig. 3.14 and Fig. 3.15, the improvement of image quality at the stronger motion phase is even higher . The results illustrate that the proposed method can obtain an acceptable 3D quality even at the non-optimal phase with stronger residual motion.

Discussion

The evaluation results show that our proposed 2D motion compensation leads to an improved 3D coronary reconstruction both for the rest phase with less motion and the phase with stronger motion. In Fig. 3.8 from simulated data and Fig. 3.14 from clinical data, the contrast and detail of the reconstructions are enhanced by the motion compensation at the optimum rest phase. In Fig. 3.9 from simulated data and Fig. 3.15 from clinical data, the improvement in image quality is also obvious for the phase with stronger motion.

The 3D artifact increases strongly with increasing gating window width for ECG-gated iterative reconstruction [7]. This may results from two reasons. First, in a wider gating window, there exists the residual motion that will corrupt the ECG-gated iterative reconstruction. Second, this iterative reconstruction often needs a segmentation of the coronaries with high precision.

Another hand, the artifact is much lower when using motion compensation, which largely compensates the residual motion. As the gating window width increases, the overall image noise decreases and the contrast of the vessels improves. Motion compensation can obtain a better result even without a high precise segmentation of the coronary artery.

Motion compensation makes the use of a wider gating window possible, giving a better coronary reconstruction quality, while avoiding motion blur and a strong increase in anatomy errors. And a good trade-off between accuracy, computation time and visual image quality, it was decided to use a gating window of 0.6 with rectangular window of the rest phase, 0.4 of the stronger motion phase.

When we perform the ECG-gated reconstruction, we should estimate the optimal phase. However, the positions of the optimal reconstruction phases in the cardiac cycle are different due to individual condition, and the precise optimal phase is difficult to calculate. This may lead to the reconstruction artifacts directly by ECG-gated method. As shown in Fig. 3.9 and Fig. 3.15, we can obtain a good quality 3D visualization even with a non-optimum phase by the proposed method.

There are also some limitations of our approach. First, for the pre-processing of the original projection, we need to set the parameters manually. The segmentation result of the coronary arteries is not with a high precision. Second, for the registration part, we also need to tune the parameters manually. Sometimes, the method is not very effective for some tiny local deformation. Last, the iterative reconstruction algorithm adopted in the step of the initial and motion compensated reconstruction costs too much computation time.

Conclusion

The proposed 2D projection-based motion compensation reconstruction method has been evaluated and compare with two reference methods [7,25]. The whole procedure of the proposed method includes forward projection, preprocessing of the acquired projection, registration, 3D reconstruction. Each step will affect the final reconstruction. We try to improve each step in our algorithm. For the forward projection, we adopt the simplified distance driven projector 𝑃 𝑆𝐷𝐷 to generate a forward projection of the initial reconstruction. For the key problem of registration, we adopt the 2D registration algorithm based on the classical mutual information (MI) combined with a novel rigidity penalty. This cost criterion avoids the extraction of the centerline of the vessels and does the registration accurately and robustly. An advanced optimization method called adaptive stochastic gradient descent (ASGD) is applied to decrease the computation time of the registration part. For the 3D reconstruction algorithm, we combine the iterative reconstruction algorithm with the 2D motion compensation.

The two reference methods are ECG-gated iterative reconstruction [7], and the latest 2D motion compensated reconstruction by Schwemmer [25]. The experiment data include the simulated data and the clinical data. We have the qualitative evaluation for both simulated and clinical data, and quantitative evaluation for the simulated data.

The method improved the ECG-gated reconstruction results for both simulated and clinical data. Improvements in the visual image quality are apparent for both rest and stronger motion phases. As the 3D reconstruction with better quality at different phases can be achieved, a series of 3D reconstructions can be used to assess dynamic properties of the coronary arteries, like displacements and speeds. The considerable improvement in the image quality results from motion compensation increases the clinical usability of 3D coronary artery reconstruction. The combination of the calibration work and 2D motion compensated reconstruction will make the 3D reconstruction of coronary artery more accurate and fast. assumptions.

For our method, because the main deviation motion may occur around 𝛼 axis, we only consider the condition that secondary angle 𝛽 is zero. It is easy to extend our models when 𝛽 is not zero, see 4.2 in the following.

The second part is about 3D reconstruction of coronary arteries by 2D motion compensation. We propose a complete procedure, which corrects the residual motion of each projection to increase the number of available projection images. The whole procedure includes maximum intensity forward projection, preprocessing of the current projection, registration, initial and motion compensated reconstruction. For the forward projection, we adopt the simplified distance driven projector (𝑃 𝑆𝐷𝐷 ) to generate the maximum intensity forward projection, which is faster than distance driven (DD) and ray driven (RD) and have comparable accuracy. For the key part of iconic registration, we adopt the registration algorithm based on the classical mutual information (MI) with the novel rigidity penalty. This cost criterion can avoid the extraction of the centerline of the vessels and do the registration accurately and robustly. Another hand, we adopt the advanced adaptive stochastic gradient descent (ASGD) to do the nonlinear optimization of intensity-based image registration, which has less computation cost and is more robust than gradient descent (GD). For the 3D reconstruction algorithm from 2D X-ray rotational coronary angiography, we adopt the MAP iterative reconstruction. Comparing with the motion compensation method by Schwemmer et al. [2], the proposed 2D motion compensation algorithm has a better 3D reconstruction quality.

As our method adopt the 2D motion compensation combined with iterative reconstruction. This kind of iterative reconstruction algorithm can reduce artifact when the projection images are few. However, this reconstruction costs too much. Indeed, we calculate the computation time of the proposed reconstruction with clinical data (𝜔 𝑐 = 0.4, 𝑡 = 𝑡 𝑟 ). The whole computation includes the image registration processing and forward projection updating. We use the double core processor of Intel I5, 8G RAM memory. The intensity based image registration cost of each projection pair is about 40 s, the 3D reconstruction costs about 5 minutes, the maximum intensity forward projection only costs 0.26 s for each projection angle. The registration time can be reduced by parallel computing,but the iterative reconstruction still costs too much time.

As described in chapter 3, it is significant and necessary to do the calibration. The combination of the two improvements (calibration and 2D motion compensation) will make the 3D reconstruction of coronary artery more fast and accurate.

Conclusion

We have improved the 3D reconstruction of coronary arteries from rotational angiography effectively by data acquisition modeling of C-arm CBCT and 2D motion compensation. For the calibration work, we developed an analytic algorithm through the rotational projections of the classical helical phantom to estimate the geometry parameters. This physical phantom is easy to be obtained in the experiment environment. Then, we transfer the geometry to the C-arm system transformation. Last, we model the rotational movement based on the estimated posterior information by image matching, we consider three representative work positions: Head side, Left side, Table30 positions. The experiment results indicate that the proposed algorithm on the helical phantom is more accurate and robust than the original algorithm. We have the results that the focal length is constant, the smooth variation of the principle point should be considered in the models. The movement models have an acceptable precision to predict the acquisition parameters at any acquisition position. Experiment tests also indicate that Head position has the least deviation. Left position has minor deviation. The two positions can be modeled as a mean rigid motion. Table30 positions has the most severe deviation,it can be modeled as a mean rigid motion added with a residual translation. These models will simplify the calibration work greatly for the clinical application.

For the 3D reconstruction work of coronary arteries, we propose the complete 2D motion compensation method for the rotational X-ray angiography data. The proposed method has improved the 3D reconstruction quality of coronary arteries. Improvements in the visual image quality are apparent for both rest phase (estimated optimal phase) and other phase with stronger motion. The contrast and details of the reconstructions are enhanced by the motion compensation. The 3D artifact increases strongly with increasing gating window width for the uncompensated reconstructions. The artifact decreases when using a wider gating window of the motion compensated reconstruction. As a tradeoff consideration, a gating window of 40% width of the cardiac cycle is an appropriate choice in the clinical applications. The improvement in the image quality increased the visibility of vessels and clinical usability of 3D coronary artery. This will be helpful for a better interventional planning. As this method can be applied to different phases, a series of 3D reconstructions (4D) can be obtained to asses some dynamic properties of the coronary arteries. The combination of calibration and 2D motion compensation will improve the 3D reconstruction image further.

Suggesting new direction of study 4.2.1 Future Study of Zeego movement models

We propose the future research plan of the movement models of Zeego system in this section. We only consider the condition of 𝛽 = 0 in chapter 2 for the simplification. The complete movement models should consider the condition of 𝛽 ≠ 0, which is widely used in the clinical application. We describe the research plan in the following.

a. Data supplement

We should supplement the data as: We assume that 𝑓 is constant and 𝑢 𝑠 , 𝑣 𝑠 are modeled only as the polynomial of 𝛼. Extrinsic matrix is built as:

𝐸 = 𝐸 0 𝐷 𝛽 𝐷 𝛼 (4.2)
where 𝐸 0 = ,𝑅 0 |𝑇 0 -, 𝐷 𝛽 = ,𝑅 𝛽 |(𝐼 -𝑅 𝛽 )𝑃 𝛽with 𝑟 𝛽 = (0, -1,0) 𝑇 , 𝑃 𝛽 = (0,0,0) 𝑇 .

𝐷 𝛼 has the same expression as described in chapter 2.

The rough movement model of projection matrix 𝑃 is built as:

𝑃 = 𝐼𝐸 (4.
3)

The projection marker position can be estimated as:

𝑚 = 𝑃𝑀 (4.4)
𝑀 is the 3D coordinates of each marker. The nearest detected 2D coordinates to the predictive 2D coordinates 𝑚 is labeled as the projection of 𝑀.

b.2 Modeling the intrinsic parameters

For the supplementary data, we have the acquisition sequences corresponding to different nonzero 𝛽 values. We estimate the acquisition geometries of the sequences by the algorithm proposed in chapter 2.

We model the principal point 𝑢 𝑠 , 𝑣 𝑠 as the polynomial function of 𝛼, 𝛽:

𝑢 𝑠 (𝛼, 𝛽) = 𝑟 𝛽 is the rotation axis of 𝑅 𝛽 , 𝑃 𝛽 is the fixed point on the axis. They can be calibrated by the estimated geometries.

b.5 refinement

We use all the acquisition data to do the refinement. The models can be extended as:

𝑀 𝑒𝑎𝑑 , 𝑀 𝑙𝑒𝑓𝑡 0,1 ={𝜃 0 , 0 , 𝛾 0 , 𝑇 0 , 𝑟 𝛽 , 𝑃 𝛽 , 𝑟 𝛼 , 𝑃 𝛼 , 𝑎 𝑗 , 𝑏 𝑗 ; 𝛼, 𝛽, 𝑓} 𝑀 𝑡𝑎𝑏𝑙𝑒30 2 ={𝜃 0 , 0 , 𝛾 0 , 𝑇 0 , 𝑟 𝛽 , 𝑃 𝛽 , 𝑟 𝛼 , 𝑃 𝛼 , 𝑝𝑇 𝑥 , 𝑝𝑇 𝑦 , 𝑝𝑇 𝑧 , 𝑎 𝑗 , 𝑏 𝑗 ; 𝛼, 𝛽, 𝑓}

This work will be done in the future.

Future Study of 2D Motion Compensated Reconstruction

We propose the future research plans of the projection-based 2D motion compensated reconstruction in this section.

We propose to apply the move advanced multi-feature mutual information (α -MI) in the projection-based 2D motion compensation reconstruction. α -MI adopts not only image intensity, but also the features that describe local image structure. The registration error is improved compared to standard MI registration [3][4][5][6].

We derive the analytical derivative of 𝛼 -MI, which makes it possible to apply in the previous ASGD optimization method.

a. Definition of 𝛼 -MI

Define 𝑧(𝑢 𝑖 ) = ,𝑧 1 (𝑢 𝑖 ), … , 𝑧 𝑑 (𝑢 𝑖 )to be a vector of dimension d containing all feature values at point 𝑢 𝑖 . 𝑧 𝑓 (𝑢 𝑖 ) is the feature vector of the fixed image 𝐼 𝐹 at point 𝑢 𝑖 , and 𝑧 𝑚 .𝑇 𝜇 (𝑢 𝑖 )/ is that of the moving image 𝐼 𝑀 at the transformed point 𝑇 𝜇 (𝑢 𝑖 ). 𝑧 𝑓𝑚 .𝑢 𝑖 , 𝑇 𝜇 (𝑢 𝑖 )/ is the concatenation of the two feature vectors: 𝑧 𝑓𝑚 .𝑢 𝑖 , 𝑇 𝜇 (𝑢 𝑖 )/ = 0𝑧 𝑓 (𝑢 𝑖 ), 𝑧 𝑚 .𝑇 𝜇 (𝑢 𝑖 )/1 (4.9)

Three kNN graphs are constructed:

Τ 𝑖 𝑓 = 1 k p
 ‖𝑧 𝑓 (𝑢 𝑖 ) -𝑧 𝑓 (𝑢 𝑖𝑝 )‖ Τ 𝑖 𝑚 (𝜇) = Τ 𝑖 𝑚 (𝜇) can be derived similarly.

For the feature vector 𝑧(𝑢 𝑖 ), features that describe the local structure of images supply supplementary knowledge, which may improve the registration. All the features are invariant to rotation and translation. All the feature images can be derived from the original image, and can be computed before registration [6][7]. This work will be done in the future.

Motion Compensation Reconstruction with Calibration on Zeego System

In Chapter 3, we observe the contribution of the calibration work in the motion compensation reconstruction. In this section, we propose some experiment plan to observe the two improvements on Zeego system. The objects of the experiment are the dynamic cardiac phantom and the helical calibration phantom.

The procedure of our experiments includes two steps. First, as the calibration should be updated every six months or after some mechanical changes [8][9]. We should do the calibration to acquire three latest acquisition sequences of the helical phantom at Head side, Left side, Table 30 positions. Second, we do the acquisition of the dynamic cardiac phantom to obtain three sequences. The position, protocol, parameter setting of the two acquisition performances should be the same. We evaluate the contributions of the two improvements. We have three comparisons.

Comparison I, we do the 3D reconstruction with motion compensation and without motion compensation. The reconstruction is performed with calibrated geometry. This comparison can prove the contribution of the 2D motion compensation.

Comparison II, we do the 3D reconstruction by 2D motion compensation with the estimated geometry and with the ideal geometry. For the 3D reconstruction, the estimated geometry can be applied in two components of the whole procedure. One is the compensated reconstruction algorithm. The estimated geometry can correct the correspondence between the 3D space and the 2D projection. The other is the maximum intensity forward projection. We can generate the forward projection by the two geometries. This comparison can prove the contribution of the calibration work.

Comparison III, we combine the calibrated geometry and the proposed motion compensation reconstruction. The result will be compared with the DynaCT result. All of this designed experiments need to be done in the future.

II 3 . 4 . 2

 342 Real Clinical Data ................................................................................. 3.5 Discussion ....................................................................................................... 3.6 Conclusion ...................................................................................................... Reference .............................................................................................................. Chapter 4 Summary and Prospect .............................................................................. 4.1. Summary ...................................................................................................... 4.1.1 Discussion ........................................................................................... 4.1.2 Conclusion .......................................................................................... 4.2. Suggesting new direction of study ............................................................... 4.2.1 Future Study of Zeego movement models .......................................... 4.2.2 Future Study of 2D Motion Compensated Reconstruction ................. 4.2.3 Motion Compensation Reconstruction with Calibration on Zeego System .......................................................................................................... Reference ............................................................................................................ Appendix .................................................................................................................... Publication List ................................................................................................... Acknowledgement ...............................................................................................

Fig. 1 .

 1 Fig.1. Système angiographique rotationnel Siemens Artis Zeego (et son mouvement). Systèmes d'angiographie RX Rotationnel L'imagerie angiographique coronarienne R-X de nouvelle génération, l'imagerie rotationnelle, permet d'obtenir des vues multiples (une série d'images) des coronaires, sous différents angles de vues, lors d'une même injection de produit de contraste et lors de la rotation continue du système d'acquisition autour du patient. Une vision pseudo 3D est fournie au chirurgien permettant d'obtenir un diagnostic le plus fiable possible. Il doit cependant continuer à se représenter mentalement cette structure en 3D. Ce système rotationnel exploite des capteurs 2-D et une rotation limitée à 180 degrés autour du patient, effectuée en l'espace de 4 à 5 battements cardiaques pour acquérir une séquence de 80 à 150 projections. Contrairement aux systèmes mono ou bi-plans, les projections coniques acquises ne sont pas soustraites. Par ailleurs, le débordement du produit de contraste dans l'aorte introduit une zone sombre qui réduit fortement le contraste entre les coronaires et les tissus. Ce système rotationnel, se caractérise par des mouvements plus sophistiqués et diversifiés du couple émetteur/détecteur RX en combinant le mouvement d'arc en C avec celle du socle en L.

Fig. 2 .

 2 Fig.2. Amélioration des différents du calibrage de l'angiographe rotationnel.

Fig. 3 .

 3 Fig.3. Amélioration des différents étapes de la reconstruction 3D coronaire par compensation de mouvement.

Fig. 4 .

 4 Fig.4. Vue d'ensemble de l'approache proposée de calibration.

Fig. 5 .

 5 Fig.5. Vue d'ensemble de l'approche propose de reconstruction avec compensation de mouvement 2D

1

 1 Coronary Heart Disease and Coronary Atherosclerosis

Fig. 1 .

 1 Fig.1.1 coronary artery blockage caused by coronary atherosclerosis

  addition, Single Photon Emission Computed Tomography (SPECT), Positron Emission Tomography (PET), Magnetic Resonance Imaging (MRI) and echocardiogram are all important non-invasive examination technologies. b. Invasive examination b.1 Intravascular ultrasound Intravascular Ultrasound (IVUS) is one new interventional ultrasound diagnosis technology that has been developed recent years. It uses a micro-ultrasound probe mounted at the apex of the heart catheter to generate the ultrasonic signals.

Fig. 1 . 2

 12 Fig.1.2 Illustration of PTCA and stent implantation

Fig. 1 . 3

 13 Fig.1.3 Biplane angiography system

Fig. 1 . 4 .

 14 Fig.1.4. Acquisition procedure by a C-arm CT.

Fig. 1 . 5 .

 15 Fig.1.5. Traditional 3D rotational C-arm system with two axes (Artis-zee family)

  Fig.1.6 Illustration of the acquisition angles of C-arm system.

Fig. 1 . 7 . 1 ( 3 ) 2 ( 4 )( 5 )( 6 )

 17132456 Fig.1.7. Movement of Artis Zeego

  Left Side, Table rotated by 30 ° (f) Left Side, Table rotated by 15 ° (g) Park

Fig. 1 . 9 .

 19 Fig.1.9. Work Positions of Artis-Zeego System. (copyright by Siemens)

Fig. 1 .

 1 Fig.1.10. The angiography imaging system model.

Fig. 1 . 11 .

 111 Fig.1.11. Illustration of ECG-gating reconstruction. 5 cardiac cycles for this acquisition, the gating windows are given for a selected phase at 70% of the R-R interval.

Fig. 1 .

 1 Fig.1.12. The procedures of the improvements of the 3D reconstruction.

Fig. 2 . 1 .

 21 Fig.2.1. Acquisition geometry of 3D cone-beam imaging system.

Fig. 2 . 2 .

 22 Fig.2.2. Principle of the simplified distance-driven method proposed in [7].

Fig. 2 . 3 .

 23 Fig.2.3. Flowchart of the proposed data calibrating algorithm

Fig. 2 . 4 .

 24 Fig.2.4. Geometry of the calibration algorithm. One general position of the scanner, world coordinate system is built on the phantom.

Fig. 2 .

 2 Fig.2.5. 2D projection image, when the angle β = 0. The sequence of alternating large and small spheres represent a binary string with an 8 bit encoding, this established a unique 2D and 3D correspondence according to a projection angle and phantom position.

a. 1

 1 Fig.2.6. Illustration to calculate the projection points on the 𝑧 axis. a.1.1 Projection of 𝑧 axis

Fig. 2 . 7 .

 27 Fig.2.7. Illustration of the projection of the axes paralleled to 𝑥 axis. 𝑘 = 1, -1 a.1.2 Calculation of 𝑎 𝑧 , 𝑏 𝑧 , 𝑐 𝑧 Since 𝑛 ∈ (1, 𝑁𝑃 -3𝑁 2 ⁄ ), the projected points are enough to estimate 𝑝 𝑧 = (𝑎 𝑧 , 𝑏 𝑧 , 𝑐 𝑧 ).Given a set of 𝒓 𝒏 𝒁 = (0,0, 𝑍 𝑛 𝑍 ) and 𝑟 𝑛 𝑍′ = (𝑢 𝑛 𝑍′ , 𝑣 𝑛 𝑍′ ), equation (2.21) could be written as:

Fig. 2 . 8 .

 28 Fig.2.8. Examples to illustrate the projection of the axes paralleled to 𝑦 axis. 𝑘 = 1, 𝑛 = 0 in a.1 to solve the equation (2.35).

c. 39 )

 39 Intrinsic and extrinsic matricesThe definition of the intrinsic matrix 𝐼: 𝑘 𝑥 , 𝑘 𝑦 represent the pixel size of the detector plane.The extrinsic matrix 𝐸:𝑅 = ,𝒖 ̂, 𝒗 ̂, 𝒏 ̂-𝑇 𝑻 = -𝑅 • 𝒔 𝐸 = ,𝑅, 𝑻-(2.40)𝒖 ̂, 𝒗 ̂, 𝒏 ̂ mean the detector orientation. 𝒔 is the position of the X-ray source.

Fig. 2 . 9

 29 Fig.2.9 Dumay's C-arm geometric model

Fig. 2 .

 2 Fig.2.10 Gorges's C-arm geometric model 𝑇 = ,0,0, 𝑆𝐴𝐷-′

Fig. 2 .

 2 Fig.2.11 Li's C-arm geometric model. 𝐸 𝛼 = ,𝑅 𝛼 |(𝐼 -𝑅 𝛼 )𝑃 𝛼 + 𝜆𝑟 𝛼 -(2.42) 𝑅 𝛼 , 𝑅 𝛽 are the rotation around the two axes, 𝑃 𝛼 , 𝑃 𝛽 are the point on the axes, 𝑟 𝛼 is the rotation axis, 𝑢 𝑖 , 𝜆 𝑖 are the coefficients. The expression of the movement model is:

Fig. 2 .Where 𝑁 1 ,

 21 Fig.2.12 The relationship between the coordinate systems in Li's model

Fig. 2 .

 2 Fig. 2.13. Geometry of Zeego System. Work positions include Head, Left, Table 30. (a) is the Head position, (b) is the Left position, (c) is Table 30 position, (d) is the geometry of the trajectory.

  𝑃 𝑒𝑎𝑑 = *(𝛼, 𝛽)|𝛽 = 0 𝛼 ∈ *-160 °, -150 °, … ,0 °, 10 °, 20 °++ 𝑃 𝑙𝑒𝑓𝑡 = *(𝛼, 𝛽)|𝛽 = 0 𝛼 ∈ *-160 °, -150 °, … , 0 °, 10 °, 20 °++ 𝑃 𝑡𝑎𝑏𝑙𝑒30 = *(𝛼, 𝛽)|𝛽 = 0 𝛼 ∈ *-80 °, -70 °, … , 0 °, … , 90 °, 100 °++ The test dataset is defined as: 𝑃 𝑒𝑎𝑑 𝑇 = *(𝛼, 𝛽)|𝛽 = 0 𝛼 ∈ *-165 °, -155 °, … ,5 °, 15 °, 25 °++ 𝑃 𝑙𝑒𝑓𝑡 𝑇 = *(𝛼, 𝛽)|𝛽 = 0 𝛼 ∈ *-165 °, -155 °, … , 5 °, 15 °, 25 °++ 𝑃 𝑡𝑎𝑏𝑙𝑒30 𝑇 = *(𝛼, 𝛽)|𝛽 = 0 𝛼 ∈ *-75 °, -65 °, … , 85 °, 95 °++

Fig. 2 . 14 .

 214 Fig.2.14. Mean projection error of the two algorithms (a) Head side position (b) Left side position (c) Table30 position. Blue start is the result by our approach. Red cross is the result by Xu's algorithm.

Fig. 2 .

 2 Fig.2.15. Source trajectory under phantom and C-arm coordinate system (a) Head side position (b) Left side position (c) Table30 position.

  Fig.2.16 2D projection error in horizontal and vertical directions of the Zeego system 𝛼(°)

  , (a) are the reconstructed results of Head position. (b) are the reconstructed results of Left position. (c) are the reconstructed results of Table30 position. Left column is the reconstruction with ideal geometries, middle is the reconstruction by classical projection matrix method, right column is the reconstruction by the proposed movement model method.

  Fig.2.17. Reconstructed 3D spheres by the geometries

*

  𝐵 𝑖 + 𝑖 , {𝐵 𝑗 } 𝑗 , *𝐵 𝑘 + 𝑘 are the B-spline function bases along space coordinates, *𝐵 𝑙 + 𝑙 is the B-spline function basis along time coordinate and 𝑝 𝑖𝑗𝑘𝑙 is the 3D vector at control point given by indices i,j,k,l, belonging to the global parameter vector p. Φ(𝑝, 𝑋, 𝑡) provides the new position after motion of any 3D point at time 𝑡.

Fig. 3 . 1 .

 31 Fig.3.1. Workflow of the projection-based 2D motion compensated reconstruction

Fig. 3 . 2 .

 32 Fig.3.2. Workflow of the proposed 2D motion compensated reconstruction

Fig. 3 . 3 .

 33 Fig.3.3. Illustration of the reconstructed phase and gating window

Fig. 3 . 4 .

 34 Fig.3.4. Illustration of the registration process.

Fig. 3 . 5

 35 Fig.3.5 Cardiac motion and phase visualization

Figure 3 .

 3 6 shows one typical example of the preprocessing of the original projection data. The projection image corresponds to the phase 𝑡 = 0 of the first cardiac cycle. (a) is the original projection, (b) is the binary mask of the vessel structure. (c) is the estimated background. (d) is the extracted vessel structure. We observe that the vessel structure is extracted well from the original projection.

Fig. 3 . 6

 36 Fig.3.6 Preprocessing of the original projection

Fig. 3 . 7 .

 37 Fig.3.7. Checkerboard images showing the qualitative evaluation of image registration

4 Fig. 3 . 9 .

 439 Fig.3.9. Reconstructions of the simulation data at the stronger motion phase 𝑡 𝑠 with different gating window widths ω 𝑐 . Left column is the result by ECG-gated iterative reconstruction [7],middle column is the reconstruction by method S[25], right column is the reconstruction by the proposed method.

  Fig.3.10. Reconstructions with calibrated and ideal geometries.

  Fig.3.12. Preprocessing of original projection data

  Fig.3.13. Checkerboard images showing the qualitative evaluation of image registration

Fig. 3 .

 3 Fig.3.15. Reconstructions of the clinical data at the phase with stronger motion 𝑡 𝑠 with different gating window widths ω 𝑐 . (a) ω 𝑐 = 0.05 (b) ω 𝑐 = 0.1 (c) ω 𝑐 = 0.2 (d) ω 𝑐 = 0.4. Left column is the result by ECG-gated iterative reconstruction [7], middle column is the reconstruction by method S [25], right column is the reconstruction by the proposed method.

Fig. 4 .

 4 Fig.4.1 illustration of the projection of the helical phantom when 𝛽 is non-zero.

1

 1 𝑃 𝑒𝑎𝑑 = *(𝛼, 𝛽)|𝛽 ∈ *-40 °, -30 °, … ,0 °, 30 °, 40 °+ 𝛼 ∈ *-160 °, -150 °, … ,0 °, 10 °, 20 °++ 𝑃 𝑙𝑒𝑓𝑡 = *(𝛼, 𝛽)|𝛽 ∈ *-40 °, -30 °, … ,0 °, 30 °, 40 °+ 𝛼 ∈ *-160 °, -150 °, … , 0 °, 10 °, 20 °++ 𝑃 𝑡𝑎𝑏𝑙𝑒30 = *(𝛼, 𝛽)𝛽 ∈ *-40 °, -30 °, … ,0 °, 30 °, 40 °+ 𝛼 ∈ *-80 °, -70 °, … , 0 °, … , 90 °, 100 °++ b. Calibrating the complete movement models b.Labeling the projection markersAs shown in Fig.4.1, the projection of the helical curve is intersected in the case of 𝛽 ≠ 0. We extend the previous movement model in chapter 2 roughly with the ideal assumptions to tackle the label problem.Intrinsic matrix is built as: 𝛼) = 𝑝𝑢 1 𝛼 3 + 𝑝𝑢 2 𝛼 2 + 𝑝𝑢 3 𝛼 + 𝑝𝑢 4 𝑣 𝑠 (𝛼) = 𝑝𝑣 1 𝛼 3 + 𝑝𝑣 2 𝛼 2 + 𝑝𝑣 3 𝛼 + 𝑝𝑣 4 (4.1)

b. 3 4

 34 Modeling 𝑅 𝑝 𝑐 , 𝑇 𝑝 𝑐 , 𝐷 𝛼We adopt the acquisition sequence when 𝛽 is zero and perform the same procedure as chapter 2. All the following extrinsic parameters should be transformed to C-arm system. Modeling 𝑟 𝛽 , 𝑃 𝛽 We adopt the acquisition sequence when α = 0 and estimate the geometries around the rotation axis 𝛽. We model the movement around 𝛽 axis.𝐷 𝛽 = ,𝑅 𝛽 |(𝐼 -𝑅 𝛽 )𝑃 𝛽 -(4.7)



  ‖𝑧 𝑚 .𝑇 𝜇 (𝑢 𝑖 )/ -𝑧 𝑚 .𝑇 𝜇 (𝑢 𝑖𝑝 )‖𝑧 𝑓𝑚 .𝑢 𝑖 , 𝑇 𝜇 (𝑢 𝑖 )/ -𝑧 𝑓𝑚 .𝑢 𝑖𝑝 , 𝑇 𝜇 (𝑢 𝑖𝑝 )/‖ (4.10)A graph-based estimator for 𝛼 -MI is defined as[5]:𝛼 -𝑀𝐼 ̂(𝜇; 𝐼 𝐹 , 𝐼 𝑀 ) (1 -𝛼), 0 < 𝛼 < 1(user defined constant). N is the number of pixels in the overlapped domain, d is the length of the feature vector.b. Analytical derivative of 𝛼 -MIWe deduce the derivative 𝜕 𝜕𝜇 ⁄ 𝛼 -𝑀𝐼 ̂(𝜇; 𝐼 𝐹 , 𝐼 𝑀 ) in this section. For compact notation, we define 𝐺 𝑖 (𝜇) ) = 𝑧 𝑓𝑚 .𝑢 𝑖 , 𝑇 𝜇 (𝑢 𝑖 )/ -𝑧 𝑓𝑚 .𝑢 𝑖𝑝 , 𝑇 𝜇 (𝑢 𝑖𝑝 )/ (4.12)Then α -𝑀𝐼 ̂(𝜇; 𝐼 𝐹 , 𝐼 𝑀 )

(4. 13 )

 13 The derivative of 𝛼 -𝑀𝐼 ̂ equals:

  𝑝 11 𝑝 12 𝑝 13 𝑝 14 𝑝 21 𝑝 22 𝑝 23 𝑝 24 𝑝 31 𝑝 32 𝑝 33 𝑝 34

]

(2.7)

We combine

(2.5)

and (2.7) to obtain a new equation (2.8) as: 𝑈𝑝 = 0 (2.8) with 𝑝 = (𝑝 11 , 𝑝 12 , … 𝑝 34 ) 𝑇

  Since 𝑟 1 and 𝑟 2 are orthonormal, we have the constraints for each view:

			13)
	with
	𝐻 = 𝐼,𝑟 1 𝑟 2 𝑡-
	[ (𝑉 11 -𝑉 22 ) 𝑇 ] 𝑑 = 0 𝑇 𝑉 12	(2.14)
	Where 𝑉 𝑖𝑗 = [ 𝑖1 𝑗1 , 𝑖1 𝑗2 + 𝑖2 𝑗1 , 𝑖2 𝑗2 , 𝑖1 𝑗3 + 𝑖3 𝑗1 , 𝑖2 𝑗3 + 𝑖3 𝑗2 , 𝑖3 𝑗3 ]	𝑇
	𝐵 = 𝜆𝐼 -𝑇 𝐼 -1 ≡ 𝜆 [	𝐵 11 𝐵 12 𝐵 13 𝐵 12 𝐵 22 𝐵 23
		𝐵 13 𝐵 23 𝐵 33

]

Note that B is symmetric, we define it by a 6D vector:

𝑑 = ,𝐵

11 , 𝐵 12 , 𝐵 22 , 𝐵 13 , 𝐵 23 , 𝐵 33 -𝑇 (2.15) For n angles of the acquired images, we can stack equation (2.14) as:

  𝒓 𝑨 = (𝑅 cos 𝐴 , 𝑅 sin 𝐴 , 𝑏 0 * 𝐴 ) 𝒓 = (𝑅 cos 𝐵 , 𝑅 sin 𝐵 , 𝑏 0 * 𝐵 ) 𝒓 𝑪 = (𝑅 cos 𝐶 , 𝑅 sin 𝐶 , 𝑏 0 * 𝐶 ) 𝒓 𝑫 = (𝑅 cos 𝐷 , 𝑅 sin 𝐷 , 𝑏 0 * 𝐷 )

	𝐴 = 𝑛 , 𝐵 = 𝑛 + 2
	𝐶 = 𝑛 + 3 , 𝐷 = 𝑛 +	(2.28)

𝑏 (2.27)

𝑛 = 1, . . 𝑁𝑃, 𝑁𝑃 = 108, 𝑛 0 =55, N=40, 𝑅 = 70𝑚𝑚, 𝑏 = 1.3𝑚𝑚. See Fig.

2

.6, we try to find a group of parallelograms, each parallelogram is composed by 4 markers 𝑟 𝐴 , 𝑟 𝐵 , 𝑟 𝐶 , 𝑟 𝐷 :

  The intersection point of 𝑟 𝐴 𝑟 𝐶 and 𝑟 𝐵 𝑟 𝐷 is: Calculation of 𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥

		2.33)
	𝑘 = 1, (𝑛 𝐴 , 𝑛 𝐵 , 𝑛 𝐶 , 𝑛 𝐷 ) ∈ (1, 𝑁𝑃).
	𝑘 = 1, we seek 𝑖 1 ∈ ,max(1,2𝑛 0 + 𝑁 -𝑁𝑃, 𝑛 0 + 1) , min(𝑁𝑃 -𝑁, 2𝑛 0 -1)-,
	𝑖 1 ≠ 𝑛 0 + 𝑁/2.
	𝑘 = -1, we seek 𝑖 2 ∈ ,max(1 + 𝑁, 2𝑛 0 -𝑁𝑃) , min(𝑁𝑃, 2𝑛 0 -𝑁 -1, 𝑛 0 -
	1)-,𝑖 2 ≠ 𝑛 0 -𝑁/2.
	𝒓 𝑖 1 𝑿 = (𝑅 cos 𝑖 1 , 0, 𝑏 0 ) 𝒓 𝑖 2 𝑿 = (𝑅 cos 𝑖 2 , 0, -𝑏 0 )	(2.34)
	Their projection points 𝑟 𝑖 1 𝑋′ , 𝑟 𝑖 2
	Given a set of 𝒓 𝑖 1 𝑿 = (𝑋 𝑖 1 𝑋 , 0, 𝑍 𝑖 1 𝑋 ), 𝒓 𝑖 2 𝑿 = (𝑋 𝑖 2 𝑋 , 0, 𝑍 𝑖 2 𝑋 ) and 𝑟 𝑖 1 𝑋′ = (𝑢 𝑖 1 𝑋′ , 𝑣 𝑖 1 𝑋′ ), 𝑟 𝑖 2 𝑋′ = (𝑢 𝑖 2 𝑋′ , 𝑣 𝑖 2

𝑋′

compose the projection of the two axes paralleled to 𝑥 axis. They can be calculated by the intersection of lines 𝑟 𝐴′ 𝑟 𝐶′ and 𝑟 𝐵′ 𝑟 𝐷′ . a.2.2 𝑋′ ), we define 𝑝 𝑥 = (𝑎 𝑥 , 𝑏 𝑥 , 𝑐 𝑥 ). Equation (2.23) could be written as: 𝐴 𝑥 𝑝 𝑥 = 𝑏

(2.35)

  𝑅 sin 𝑗 , 𝑏 0 /2) Their projection points 𝑟 𝑗 1 𝑌 ′ , 𝑟 𝑗 2 𝑌 ′ , 𝑟 𝑗 3 𝑌 ′ , 𝑟 𝑗 4 𝑌 ′ compose the projection of the axes parallel to x axis, they can be calculated by the intersection of lines 𝑟 𝐴′ 𝑟 𝐶′ and 𝑟 𝐵′ 𝑟 𝐷′ . Calculation of 𝑎 𝑦 , 𝑏 𝑦 , 𝑐 𝑦

	𝒀 = (0, 𝑅 sin 𝑗 , -𝑏 0 /2) 𝒓 𝒋 𝟑 𝒓 𝒋 𝟒 𝒀 = (0, 𝑅 sin 𝑗 , -3𝑏 0 /2)	(2.38)
	a.3.2
	𝒓 𝒋 𝟐 𝒀 = (0, 𝑅 sin 𝑗 , 3𝑏 0 /2)

Table 2 .

 2 𝑡 𝑤 , 𝜃, 𝜖, ‖𝐶‖ 𝑑 𝑥 , 𝑑 𝑦 , 𝑅 𝑑 , 𝑡 𝑑 , 𝑡 𝑖 𝐸 0 , 𝑟 𝛼 , 𝑃 𝛼 , 𝑟 𝛽 , 𝑃 𝛽 , 𝜆 𝑖 𝑟 𝛼 , 𝑃 𝛼 , 𝐸 𝑎𝑠 , 𝐸 0 𝑟, 𝑓, 𝑎 𝑗 , 𝑏 𝑗

	22]

1. The comparisons of the C-arm movement models.

  𝜎 𝑣 𝑐 , 𝜎 𝑍 𝑠 are the standard deviation of 𝑣 𝑐 𝑐 , 𝑍 𝑆 𝑐 . 𝑍 𝑠 ̅̅̅ is the mean of 𝑍 𝑆 𝑐 . We use thePowell method[30] to realize the optimization. The objective function should be calculated over the calibration angles with minor deviation. We have the criterion to choose the calibration angles with minor deviation:Then, angles 𝛼, 𝛼 + 𝛼 𝑠 are chosen as the angles with minor deviation.𝑆 𝛼 is the source position at angle 𝛼, 𝑇 is the sampling interval angle between two adjacent calibration angles. 𝐿 𝑡 is the threshold. The initialization of 𝑹 𝒑 𝒄 can be estimated at the reference angle 𝛼 = 0.

	𝑹 𝑺 𝒄 = [	0 0 -1 0 0 -1 ]	(2.54)
			-1 0	0
	𝑺 is the estimated extrinsic matrix. We adopt the source trajectory of the chosen 𝑹 𝒑 calibration angles to compute the circle fitting, the circle center is 𝑷 𝒊𝒔𝒐 𝒄 = 𝟎, 𝒑 . 𝑷 𝒊𝒔𝒐 from (2.53), we have the initialization of 𝑻 𝒑 𝒄 :
		𝑻 𝒑 𝒄 = -𝑹 𝒑 𝒄 𝑷 𝒊𝒔𝒐 𝒑	(2.55)
	Last, we transform the extrinsic matrix:
		𝑸 = [	𝑹 𝒑 𝒄 0	𝒄 𝑻 𝒑 1	]
		𝑬 𝒄 = 𝑬 𝒑 𝑸 𝑻	(2.56)
				𝑍 𝑠 ̅̅̅ 2 + 𝜎 𝑍 𝑠	2	(2.51)
	𝐿 𝑆 𝛼 𝑆 𝛼+T 𝑠 = ‖𝑷 𝑺 𝜶 𝒑 -𝑷 𝑺 𝜶+𝑇 𝑠 𝒑	‖
		𝜃 = 𝛼 𝑠 /2 𝐿 = 2𝑆𝐴𝐷 𝑐 sin 𝜃
	If	|𝐿 𝑆 𝛼 𝑆 𝛼+1 -𝐿| ≤ 𝐿 𝑡	(2.52)
		𝑹 𝒑 𝒄 = 𝑹 𝑺 𝒄 𝑹 𝒑 𝑺	(2.53)
	with the ideal assumption

  𝑰 𝜶 𝑬 𝜶 𝑬 𝜶 = 𝑬 𝟎 𝑫 𝜶 𝑬 𝟎 = ,𝑹 𝟎 |𝑻 𝟎 -(2.57) 𝑰 𝜶 and 𝑬 𝜶 is the intrinsic and extrinsic matrix at angle 𝛼 . 𝑬 𝟎 (𝛼 = 0) is the reference matrix. 𝑫 𝜶 means the motion around axis 𝛼. In 2.4.1 and 2.4.2, we have the analyzed characteristics of 𝑰 𝜶 and 𝑫 𝜶 . In the following 2.4.3, we test many possible models. We model 𝑰 𝜶 by 𝑢 𝑠 , 𝑣 𝑠 , 𝐹. 𝐹 is constant, 𝑢 𝑠 , 𝑣 𝑠 can be modeled as polynomial functions:

  59)𝑹 𝜶 is the mean rotation matrix around axis 𝛼, 𝑷 𝜶 is the fixed point on the rotation axis 𝒓 𝜶 .𝑫 𝜶 can be modeled as a mean rotation motion with minor residual translation for the Table30 position:𝑫 𝜶 = ,𝑹 𝜶 |(𝑰 -𝑹 𝜶 )𝑷 𝜶 + 𝑻 𝜶 -𝑇 𝑥 (𝛼), 𝑇 𝑦 (𝛼), 𝑇 𝑧 (𝛼) are the components of 𝑻 𝜶 . 𝑝𝑇 𝑥 , 𝑝𝑇 𝑦 , 𝑝𝑇 𝑧 are the coefficients vectors of polynomial expansion of each component, M is the order.We have the following three representative models to build the projection matrix:𝑀 𝑒𝑎𝑑 , 𝑀 𝑙𝑒𝑓𝑡 0,1 =*𝜃 0 , 0 , 𝛾 0 , 𝑻 𝟎 , 𝒓 𝜶 , 𝑷 𝜶 , 𝑝𝑢, 𝑝𝑣; 𝛼, 𝐹+ (2.62) ={𝜃 0 , 0 , 𝛾 0 , 𝑻 𝟎 , 𝒓 𝜶 , 𝑷 𝜶 , 𝑝𝑇 𝑥 , 𝑝𝑇 𝑦 , 𝑝𝑇 𝑧 , 𝑝𝑢, 𝑝𝑣; 𝛼, 𝐹}

		(2.60)
	𝑻 𝜶 is the uncertain residual translation:
	𝑻 𝜶 = [𝑇 𝑥 (𝛼), 𝑇 𝑦 (𝛼), 𝑇 𝑧 (𝛼)] 𝑇
	𝑇 𝑥 (𝛼) = ∑ 𝑝𝑇 𝑥𝑖 𝛼 𝑀-𝑖 𝑀 𝑖=1 𝑇 𝑦 (𝛼) = ∑ 𝑝𝑇 𝑦𝑖 𝛼 𝑀-𝑖 𝑀 𝑖=1 𝑇 𝑧 (𝛼) = ∑ 𝑝𝑇 𝑧𝑖 𝛼 𝑀-𝑖 𝑀 𝑖=1	(2.61)
	𝑀 𝑡𝑎𝑏𝑙𝑒30
	2

Table 2 .

 2 3. Intrinsic parameter analysis calibration angles of the two algorithms in Fig.2.14. We calculate the mean computation time of the two algorithms in Table 2.2. From Fig.2.14 and Table

			Computation time (s)
			Head	Left	Table30
		Xu's	12.59	12.54		12.66
		algorithm
		proposed	10.29	10.40		10.42
		Table 2.2. Computation time of the two algorithms
		mean std		Mean fitting error (order=2,3,4)
	Head	Left	Table30	Head		Left	Table30
	𝐹 1196.67 2.9	1197.10 3.3	1196.89 3.	2.59/2.60/2.5	3.38/3.37/3.4	2.48/1.54/1.4
	9	1	2	7		0	0
	𝑢 𝑠 238.32 5.51 242.53 6.06 243.21 5.0	2.81/2.79/2.8	4.55/4.48/4.5	3.59/3.42/3.3
			1	2		4	6
	𝑣 𝑠 316.88 6.75 318.52 7.18 313.83 5.1	4.37/4.30/3.9	3.86/3.41/3.1	3.56/3.63/3.5
			9	5		4	2

  𝜃 0 , 0 , 𝛾 0 are three Euler angles of 𝑹 𝟎 . 𝑻 𝟎 is the translation vector of 𝑬 𝟎 . 𝑝𝑇 𝑥 , 𝑝𝑇 𝑦 , 𝑝𝑇 𝑧 are the coefficient vectors of polynomial of the residual translation vector 𝑻 𝜶 , see (2.61). 𝑝𝑢, 𝑝𝑣 are the coefficient vectors of polynomial of estimated 𝑢 𝑠 , 𝑣 𝑠 , see (2.58).

	𝛼, 𝐹, 𝑢 𝑠 , 𝑣 𝑠 , 𝑆𝐴𝐷+	(2.68)
	𝑚 0
	0 =*𝜃 0 , 0 , 𝛾 0 , 𝑻 𝟎 , 𝒓 𝜶 , 𝑷 𝜶 ; 𝛼, 𝐹, 𝑢 𝑠 , 𝑣 𝑠 +	(2.69)
	𝑚 1
	1 =*𝜃 0 , 0 , 𝛾 0 , 𝑻 𝟎 , 𝒓 𝜶 , 𝑷 𝜶 , 𝑝𝑢, 𝑝𝑣; 𝛼, 𝐹+	(2.70)
	𝑚 2
	2 ={𝜃 0 , 0 , 𝛾 0 , 𝑻 𝟎 , 𝒓 𝜶 , 𝑷 𝜶 , 𝑝𝑇 𝑥 , 𝑝𝑇 𝑦 , 𝑝𝑇 𝑧 , 𝑝𝑢, 𝑝𝑣; 𝛼, 𝐹}	(2.71)

𝛼, 𝐹, S D are read from sensor data. 𝑢 𝑠 , 𝑣 𝑠 are the center pixels of the detector.

Table 2 .

 2

		-70	-60	-50	-40	-30	30	40	50	60	70	80	90	100
	Table30 4.95 4.97 4.74 7.85 9.92 2.35 2.16 3.32 3.88 3.54 4.70 5.94 5.58
						(a) calibration dataset
	𝛼(°)	-105 -95	-85	-75	-65	-55	-45	-35	-25	-15	5	15	25
	Head	1.28 1.41 1.18 1.49 1.12 1.75 1.30 1.53 1.69 1.59 1.94 1.40 1.11
	Left	2.20 2.48 3.31 2.38 2.53 2.35 2.26 1.56 1.23 1.32 2.64 2.55 3.26
	𝛼(°)	-65	-55	-45	-35	-25	-15	-5	45	55	65	75	85	95
	Table30 4.93 4.37 4.96 9.27 9.18 7.52 4.09 5.11 2.93 3.44 4.24 3.96 6.07
						(b) test dataset

1.71 1.15 1.52 1.72 1.39 1.85 1.55 1.87 1.39 Left 2.35 2.38 2.44 2.48 2.45 2.22 2.14 2.18 1.22 1.24 1.43 1.70 2.20 𝛼(°) 6 2D projection error of Zeego system (in pixels)

Table 2 .

 2 17. Reconstructed 3D spheres by the geometries

		Ideal Geometry	Classical projection	Proposed movement
				matrix method	model method
		Head	Left Table30 Head Left Table30 Head Left Table30
	Computation	0.013 0.014 0.0095	64.61 63.54 64.55	26.79 26.30 27.97
	time(s)

7 Computation time of the algorithms

  𝑀,𝜇 (𝐼 𝐹 , 𝐼 𝑀 )𝑙𝑜𝑔 2 𝐼 𝐹 ,𝐼 𝑀 𝜌 𝐹,𝑀,𝜇 (𝐼 𝐹 ,𝐼 𝑀 ) 𝜌 𝐹,𝜇 (𝐼 𝐹 )𝜌 𝑀,𝜇 (𝐼 𝑀 )

  𝑅 and 𝑡 are the rotation matrix and translation vector respectively. 𝑟 𝑖𝑗 is the component of 𝑅. Three conditions on 𝑇 𝜇 (𝑥) can be derived.

					𝜇 (𝑥) = 𝑥 + 𝑇(𝑥) = 𝑅𝑥 + 𝑡	(3.32)
	Affinity condition:
					𝐶 𝑘𝑖𝑗 (𝑥) =	𝜕 2 𝑇 𝜇 𝑘 𝜕𝑥 𝑖 𝜕𝑥 𝑗 (𝑥) , for all 𝑘, 𝑖, 𝑗 = 1,2	(3.33)
	For the rigid transformation, 𝐶 𝑘𝑖𝑗 (𝑥) = 0.
	Orthonormality condition:
	O𝐶 𝑖𝑗 (𝑥) = ∑	2 𝑘=1	. 𝜕𝑇 𝜇 𝑘 𝜕𝑥 𝑖 (𝑥)	+ 𝛿 𝑘𝑖 / ( 𝜕𝑇 𝜇 𝑘 𝜕𝑥 𝑗 (𝑥)	+ 𝛿 𝑘𝑗 ) -𝛿 𝑖𝑗	, for all 𝑖, 𝑗 = 1,2 (3.34)
	For the rigid transformation, we have ∑	2 𝑘=1	𝑟 𝑘𝑖 𝑟 𝑘𝑗	= 𝛿 𝑖𝑗 and	𝜕𝑇 𝜇 𝑖 𝜕𝑥 𝑗 (𝑥)	= 𝑟 𝑖𝑗 -𝛿 𝑖𝑗 , then,
	O𝐶 𝑖𝑗 (𝑥) = 0.
	Properness:
						𝑃𝐶(𝑥) = det(𝑅) -1
						𝑟 𝑖𝑗 =	𝜕𝑇 𝜇 𝑖 𝜕𝑥 𝑗 (𝑥)	+ 𝛿 𝑖𝑗 .	(3.35)
	For the rigid transformation, 𝑃𝐶(𝑥) = 0.
	We have the definition of the rigidity penalty term Ψ(𝑇 𝜇 ) as the sum of all these
	conditions squared:
		Ψ(𝑇 𝜇 ; 𝐼 𝑀 ) =	1 ∑ 𝑐(𝑥+𝑇(𝑥)) 𝑥	∑ 𝑐(𝑥 + 𝑇(𝑥)) 𝑥
	× {𝑐 𝐴𝐶 ∑	𝑘,𝑖,𝑗	𝐴𝐶 𝑘𝑖𝑗 (𝑥) 2 + 𝑐 𝑂𝐶 ∑ 𝑂𝐶 𝑖𝑗 (𝑥) 2 + 𝑐 𝑃𝐶 𝑃𝐶 𝑖,𝑗	(𝑥) 2 }	(3.36)

  3

  ((𝑥 -𝑥 𝑘 )/𝜎) 𝑥 𝑘 are the control points. 𝛽 3 (𝑥) is the cubic multidimensional B-spline polynomial. 𝑙 𝑘 is the B-spline coefficient vectors. 𝜎 is the B-spline control point spacing, a small spacing of control points allows modeling of highly local non-rigid deformations. 𝑁 𝑥 is the set of all control points within the compact support of the B-spline at 𝑥, this means that the transformation of a point can be computed from only a couple of surrounding control points. The control points 𝑥 𝑘 are defined on a regular grid of the fixed image. The parameters 𝜇 are formed by the coefficients 𝑙 𝑘 . 𝑛 1 , 𝑛 2 are the number of control points of each direction of 2D image, 𝑠 1 , 𝑠 2 are the image size of each direction, then 𝑛 1 ≈ 𝑠 1 𝜎 1 ⁄ , 𝑛 2 ≈ 𝑠 2 𝜎 2 ⁄ . The number of parameters is 𝑛 1 × 𝑛 2 × 2.

	(3.38)

Table . 3

 . introduced the metric 𝑄 3𝐷 (𝑓) to evaluate the accuracy of the reconstructed cardiac vasculature morphology. The simulated dataset has 133 .1. 𝑄 3𝐷 (𝑓) values of the reconstruction results. m.c. means motion compensation. projection images and each projection corresponds to one motion state. CAVAREV platform provides the static binary reference volume 𝑓 𝑖 𝑟 representing the ground truth of the coronary arteries at each motion state i. 𝑓 is the 8 bit quantization reconstruction. The motion state based quality metric 𝑄 𝑖 (𝑓) can be computed by the following formulations:

	ω 𝑐 of the rest phase	0.15	0.3	0.45	0.6
	𝑡 𝑟
	ECG-gated iterative	0.741	0.733	0.702	0.692
	recon. (no m.c.) [7]
	Method S [25]	0.750	0.758	0.774	0.780
	Proposed method	0.762	0.767	0.792	0.813
	ω 𝑐 of the stronger	0.05	0.1	0.2	0.4
	motion phase 𝑡 𝑠
	ECG-gated iterative	0.718	0.715	0.707	0.705
	recon. (no m.c.) [7]
	Method S [25]	0.763	0.766	0.769	0.771
	Proposed method	0.770	0.773	0.778	0.781
	𝑄 𝑖 (𝑓) = 𝑚𝑎𝑥 𝑎∈*0,…,255+ 𝐷𝑠𝑐(𝑓 𝑖 𝑟 , 𝑇(𝑓, 𝑎))	(3.52)

  Quantitative evaluation of the reconstructionWe show the evaluation metric 𝑄 3𝐷 (𝑓)[START_REF] Rohkohl | CAVAREV-an open platform for evaluating 3D and 4D cardiac vasculature reconstruction[END_REF] of the two kinds of reconstructions in Table.3.2. From Table.3.2, we see that 𝑓 𝐺 𝑐 has a better metric value than 𝑓 𝐺 𝑖 . The calibration work has a significant contribution in the reconstruction result. The combination of the calibration work and 2D motion compensated reconstruction will make the 3D reconstruction of coronary artery more accurate and fast.

	ω 𝑐 of the rest phase	0.15	0.3	0.45	0.6
	𝑡 𝑟
	recon. 𝑓 𝐺 𝑖	0.688	0.690	0.702	0.706
	recon. 𝑓 𝐺 𝑐	0.762	0.767	0.792	0.813

Table.3.2. 𝑄 3𝐷 (𝑓) values of the reconstructions with calibrated and ideal geometries.

b.

  And the derivative of Τ 𝑖 𝑓𝑚 (𝜇) and 𝑇 𝑖 𝑚 (𝜇) are:

	𝜕 𝜕𝜇 𝑗	𝐺 𝑖 (𝜇) = (	𝜕 𝜕𝜇 𝑗	Τ 𝑖 𝑓𝑚 (𝜇) -	1 2	Τ 𝑖 𝑓𝑚 (𝜇)𝑇 𝑖 𝑚 (𝜇) -1 𝜕 𝜕𝜇 𝑗	𝑇 𝑖 𝑚 (𝜇)) √𝑇 𝑖 𝑓 𝑇 𝑖 𝑚 (𝜇) ⁄	(4.15)
						𝜕 𝜕𝜇 𝑗	Τ 𝑖 𝑓𝑚 (𝜇) =	1 p k 	𝜕 𝜕𝜇 𝑗	√𝑑 𝑖𝑝 𝑓𝑚 (𝜇) 𝑇 𝑑 𝑖𝑝 𝑓𝑚 (𝜇)
				=	1 p k 	1 2	‖𝑑 𝑖𝑝 𝑓𝑚 (𝜇)‖	-1 𝜕 𝜕𝜇 𝑗	.𝑑 𝑖𝑝 𝑓𝑚 (𝜇) 𝑇 𝑑 𝑖𝑝 𝑓𝑚 (𝜇)/
				=	1 p k 			𝑑 𝑖𝑝 𝑓𝑚 (𝜇) 𝑇 ‖𝑑 𝑖𝑝 𝑓𝑚 (𝜇)‖	𝜕 𝜕𝜇 𝑗	𝑑 𝑖𝑝 𝑓𝑚 (𝜇)
	=	1 p k 	𝑑 𝑖𝑝 𝑓𝑚 (𝜇) 𝑇 ‖𝑑 𝑖𝑝 𝑓𝑚 (𝜇)‖	[	𝜕𝑧 𝑚 𝜕𝑥 𝑖	.𝑇 𝜇 (𝑢 𝑖 )/	𝜕𝑇 𝜇 𝜕𝜇 𝑗	(𝑢 𝑖 ) -	𝜕𝑧 𝑚 𝜕𝑥 𝑖	.𝑇 𝜇 (𝑢 𝑖𝑝 )/	𝜕𝑇 𝜇 𝜕𝜇 𝑗	(𝑢 𝑖𝑝 )]
													(4.16)
	𝜕
	𝜕𝜇 𝑗
													 𝐺 𝑖 (𝜇) 2𝛾-1 𝜕 𝜕𝜇 𝑗	𝐺 𝑖 (𝜇)	(4.14)
	The derivative of 𝐺 𝑖 (𝜇) is written as:
													108

the most severe deviation and should be modeled with a residual translation part.

Acknowledgement

I would like to express my gratitude to all those people who helped me during the writing of this dissertation.

My deepest gratitude goes first to my supervisors, Professor Limin Luo, Christine Toumoulin, Jean-Claude Nunes. Thanks for their constant encouragement and guidance. Thank Jean-Claude, who helped me walk though all the stages of the writing of this dissertation. Thank Professor Luo and Chrristine, without their consistent and illuminating instruction, this dissertation could not have reached its present form.

Second, I would like to thank all the colleagues in the lab. Professor Huazhong Shu, Yang Chen, Miguel Castro, Ahmed Oukli, Helene Feuillatre, Jie Zhu. They all helped me a lot.

Thank my family, my close friends. Thanks for their considerations all through these years.

Last, I would like to thank China Scholarship Council (CSC). I had been in France for 18 months to finish the experiments with their financial support. Thank the hospital of Pontchaillou (Rennes, France) to provide all the clinical data in this dissertation. Thank myself for being persistent and patient. Thank this experience in my life.

(d) ω 𝑐 = 0.6 Fig. 3.8. Reconstructions of the simulation data at the rest phase 𝑡 𝑟 with different gating window widths ω 𝑐 . Left column is the result by ECG-gated iterative reconstruction [7], middle column is the reconstruction by method S [25], right column is the reconstruction by the proposed method. In this thesis, we have focused on the data acquisition modeling of Artis-zeego system and motion compensated reconstruction of coronary arteries from X-ray rotational angiography.

The first part about calibration work, we propose a new and complete data calibrating method of Artis-Zeego C-arm system. First, we proposed an extended analytic algorithm based on the classical helical phantom to estimate the geometry parameters. For the classical helical phantom, at beginning, it was designed for the classical projection matrix algorithm and parametric iterative calibration algorithms. It was convenient to be obtained in the experimental environment than other specific phantom. The two classical calibration algorithms can estimate the projection matrix accurately, however when we try to extract the intrinsic and extrinsic parameters, it is difficult to obtain a stable solution. For our method, we not only adopt the convenience of the helical phantom, but also can estimate the intrinsic and extrinsic parameters accurately and conveniently.

Second, in order to model the complex motion of Artis-Zeego, we transfer all the geometry parameters to the nominal C-arm system by minimizing the oscillation from the rotation displacement around the physical phantom. All the geometries are independent of the phantom placement, which is very important for some clinical function of this interventional system. This method is easier to perform than Ford [1] and quite suitable for our Zeego system.

Last, as Artis-Zeego system which offers new potentialities for image guided interventions is more complex and flexible than the traditional C-arm. There are two obvious differences between this novel interventional C-arm and the traditional one. One difference is that this advanced interventional system can achieve more work position results from its multi-axes, while the traditional C-arm achieves limited region with two rotation axes. The work positions of Artis-Zeego include Head side, Left side, Right side, Table15, Table30 position. The work position can't be set arbitrarily, instead, we just can choose one of them. In this dissertation, we analyze three representative work positions separately. The other, the motion of this 6 axes robotic system is more sophisticated, it is not the simple rigid motion or the rigid motion with some deviation along the rotation axis. To resolve this problem, we adopt the posterior information. We consider the intrinsic and extrinsic parameters separately and propose different movement models based on the estimated posterior information through image matching. The models are more convincing than any prior

Appendix Publication List

Si Li, Yang Chen, Limin Luo, Christine Toumoulin, Dictionary learning based sinogram inpainting for CT sparse reconstruction, Optik -International Journal for Light and Electron Optics 01/2014,125 (12):2862-2867 Helene Feuillatre, Jean Claude Nunes, Li Si, Christine Toumoulin, Coronary tree matching based on association graphs with artificial nodes, 6 th european conference of the international federation for medical and biological engineering, 7-11, September, 2014, Dubrovnik, Coatia