1 Introduction
The reconstruction of 3D attenuation coefficient images from 2D projections is more efficient than the reconstruction of slices from a line detector array 〚1, 2, 3〛. The X-ray energy is much more efficiently used on 2D detectors. Moreover 3D homogeneous reconstruction can be obtained avoiding the difficult problem of inter slice interpolation encountered in 3D imaging from a set of 2D reconstruction. A wide range of work have been devoted to 3D reconstruction by a very active community 〚1, 3–7〛. The development of rotational angiography, the use of radiology in operating rooms for acquiring 3D information 〚8〛, the development of new digital X-ray sensors technology 〚9〛, increase the interest for 3D imaging. Imaging systems need to be calibrated for obtaining quantitative information 〚10〛. Moreover, distortions due to magnetic fields occur on image intensifiers. These distortions depend on the position and orientation of the sensor and on the magnetic field perturbations during the measurement. Thus, the image intensifier distortion of mobile imaging systems need to be dynamically estimated and corrected.
A lot of calibration devices are based on a mesh of opaque markers of known geometry. The mesh is positioned between the detector and the X-ray source 〚8, 10〛. Markers are detected in each X-ray projection. The acquisition geometry of the imaging system and the distortion of the projection can be then estimated from the deformation of the marker mesh geometry. Unfortunately, because of the marker opacity, the information hidden by the marker is lost on the 2D projection.
This drawback can be overcome in the context of 3D reconstruction from 2D projections. Lost information can be hopefully interpolated from other data in the 3D set of successive projections. Interpolation formulas can be obtained by Fourier series based on sampling conditions of the data (here the parallel-beam X-ray transform or the cone-beam X-ray transform). Standard schemes are 2D rectangular sampling on the detector combined by equidistant angular sampling. In the next section we show that the X-ray transform can be more efficiently sampled. Efficient sampling schemes are subsets of standard schemes. Thus standard schemes are redundant. The main idea of the third section is to exploit this redundancy in order to design a new calibration marker mesh. With the use of this mesh, essentially no information is lost in the sense of Shannon. We propose some numerical experiments and real experiments on dry bones. In the fourth section we present a Fourier interpolation scheme for retrieving the lacking data. For simplicity, the idea is developed on the X-ray transform (parallel geometry) but it can be extended numerically to the cone-beam transform.
2 Efficient sampling
The original idea of efficient interlaced sampling in 2D tomography is due to Cormack 〚11〛. The approach of Rattey and Lindgren 〚12〛 is based on Fourier analysis: they showed that the essential support K of the Fourier transform of the Radon transform of an essentially b-band-limited function has a bow-type shape. This geometrical property can be used in order to satisfy the generalized Shannon non-overlapping conditions of the sets associated to sampling on the lattice generated by the 2×2 non-singular matrix W (see 〚13〛). Indeed, if the Shannon non-overlapping conditions is satisfied for a non-singular matrix W and a set K in the Fourier space, then the Fourier analysis yields the following interpolation formula for a function g ( sufficiently regular on ) from the sampling points
1 |
Thus if is the essential support of i.e., can be negligible, then the interpolation error is low. In 2D tomography, because of the bow-tie shape of the interlaced sampling 〚14〛 leads to a more compact non-overlapping scheme of sets than the most compact standard sampling scheme generated by a standard diagonal matrix (2D orthogonal grid).
In 〚15〛, the previous approach is generalised to establish the sampling conditions of the 3D X-ray transform of an essentially b-band-limited function
As g is periodic in its first variable, its Fourier transform is defined by:
The essential support of the Fourier transform of is the set shown in Fig. 1. It is precisely shown 〚15〛 that if f is essentially b-band-limited ( negligible) then is negligible (the sampling error is driven by this term, see 〚14〛).
The minimal standard 3D orthogonal grid satisfying the non-overlapping conditions of the sets is generated by the matrix WS, see (4). The interlaced scheme satisfies also the non-overlapping condition of the sets with twice less data than the best standard grid, see 〚16〛
4 |
3 3D virtually invisible marker geometry
The main idea of our work is to make use of the redundancy of standard grids. Opaque markers are introduced in 3D (set of 2D projections) on such a way that, theoretically, the measured data are sufficient to interpolate the data which are not available because of the markers. We first remark that:
As 0<ϑ′<ϑ<1, the matrix satisfies the non-overlapping conditions of the set In practice ϑ is chosen close to 1, thus almost the same efficiency is obtained by as by because is then close to 1. Thus, if half of the standard points located in the translated scheme are covered by opaque markers, then the half remaining scheme is enough to recover the full information. The same remark could be done with an hexagonal grid associated with the scheme 〚5〛. This idea is illustrated in Fig. 2 for an interlaced scheme. We solve the problem by reconstructing from the interlaced data. We then compute the full projections from the solution in order to show that the lacking data (supposed to be covered by markers) are well interpolated from the available interlaced data.
In practice, the calibration grids are sparse. However, in order to ensure that no information is lost, we interlace the calibration markers between two projections according to the results of previous section. The even projection are first measured. Afterwards we shift the calibration grid of one sampling unit (here the size of the ball) in the direction s. Then the odd projections can be processed. A simple algebraic approach is performed for the reconstruction from these data. We use a conjugate gradient on the regularized least squares problems
In Figs. 3 and 4 we show reconstructions from both simulated and real data. In both cases, the X-ray source trajectory is half a circle around the measured object. Measurement geometry with image intensifiers are actually divergent. The distance of the source to the axis of rotation is approximatively seven times the diameter of the cylinder of reconstruction, thus the geometry is in fact almost parallel. In order to avoid a fast oscillating calibration marker grid, we can first acquire the even projection angles, then we can shift the calibration marker grid only once and finally we can acquire the odd projection angles.
4 Interpolation on a standard scheme from interlaced data
Reconstruction from an irregular grid, such as in the presence of opaque markers, can also be achieved by interpolating the data on the complete standard grid and then using standard algorithms, as proposed in 〚19〛. In the following, we consider the Fourier interpolation on a standard scheme from an interlaced scheme in 2D tomography. This could be easily generalized to the 3D X-ray transform.
Let us suppose that half of the complete standard data are covered by markers. The geometry of the marker positions in the sinogram has to be interlaced such that the complementary interlaced data are still available. In this case, no information is lacking from the Shannon theory point of view. We propose here a Fast Fourier interpolation on a complete scheme from the data on an interlaced scheme.
In the general setting, the interlaced matrix WI corresponding to a standard scheme
From the Poisson formula associated to the interlaced scheme, we have:
8 |
In order to use the Fast Fourier Transform, we want to compute according to (8) with:
9 |
With the previous notations, we have:
Proof. The group is a subgroup of
10 |
We now remark that:
Using now
Thus can be computed from two FFTs of dimension based on the grid points and respectively. Indeed we have
The value of are then obtained by inverse FFT of the 2D grids points
In 2D tomography we can use where p is the number of projections on and h2=2/q (the reconstruction radius is normalised to 1). The sampling condition is that p must be slightly larger than (the band limit of f), see 〚14〛. In Fig. 5 we show a numerical example of the Fourier interpolation based on the previous development. This could be easily generalised to 3D problems.
5 Conclusion
In this work we have proposed interlaced geometries of markers for the dynamic calibration of image intensifiers. The design of these geometries is based on spectral properties of the 3D X-ray transform. For 2D or 3D reconstruction, these interlaced geometries of markers in the projections have almost no effect on the reconstruction. From the Shannon theoretical point of view, essentially no information is lost. We have illustrated this result with 3D reconstructions from both simulated and real data. We have given a Fast Fourier Transform based method for the interpolation on standard scheme from interlaced scheme in 2D. This interpolation scheme can be easily generalized to 3D. Interpolating on standards schemes from the available data first, allows then to use classical algorithms for the reconstruction.
Acknowledgements
We would like to thank Dr. Pierre Bessou from the neuroradiology department of the Grenoble University Hospital for providing the real data. This work have been supported by a grant of the Région Rhône-Alpes within the project “Santé et HP: de la Vision au Pilotage” and the project “ADéMo”.
Version abrégée
L'objectif de ce travail est de proposer de nouveaux schémas de grille de calibrage pour les amplificateurs de brillance dans le cadre de leur utilisation pour la reconstruction 3D en imagerie interventionnelle. Depuis plus de dix ans, la chirurgie assistée par ordinateur a connu de nombreux succès et a permis d'améliorer la précision de certains gestes pointus, tels que le vissage pédiculaire dans la chirurgie de la vertèbre. Ces progrès sont en partie dus à l'utilisation quantitative des images médicales lors des interventions chirurgicales. Les amplificateurs de brillance fournissent des images radiologiques en salle d'opération. Les distorsions dues au champ magnétique doivent être mesurées et corrigées. La technique habituelle de correction de ces distorsions est basée sur l'utilisation de mires de calibrage constituées d'une grille de marqueurs radio-opaques, qui permet d'estimer la distorsion (par la déformation de la grille), puis de corriger l'image. Dans le contexte de la reconstruction 3D, un grand nombre de projections radiographiques est transformé par un algorithme de reconstruction en une image 3D des coefficients d'atténuation de la région radiographiée. Nous nous sommes fixé l'objectif de construire une mire de calibrage furtive, dans le sens où cette mire serait visible dans chacune des projections et permettrait de corriger la distorsion dynamiquement sur chaque projection radiographique, mais ne supprimerait essentiellement pas d'information au sens de Shannon dans l'ensemble des projections. Comme nous le verrons, cet objectif peut être atteint par l'analyse des conditions d'échantillonnage du modèle de la mesure.
Après une correction de gain et d'offset, suivie d'une transformation logarithmique des données (selon la loi de Lambert–Beer, l'atténuation est exponentielle), le modèle utilisé pour la reconstruction en tomographie 3D est celui de la transformée en rayons X. Dans notre travail, nous limitons l'exposé à la géométrie parallèle, mais la démarche pourrait être généralisée aux géométries coniques, au moins numériquement. Dans le cas de la géométrie parallèle, si nous restreignons la trajectoire de la source et du détecteur à un cercle, nous pouvons mesurer la fonction d'atténuation à reconstruire par :
Dans un premier temps, nous rappelons les conditions d'échantillonnage de g lorsque f est supposée essentiellement b-bande limitée Elles sont basées sur l'identification du support essentiel de la transformée de Fourier de g. On peut montrer que les schémas standards (équidistants en φ, s et t), c'est-à-dire les grilles rectangulaires dans l'espace des paramètres, sont redondants. Le meilleur schéma standard (ensemble des points ) est presque deux fois moins performant que le schéma entrelacé (ensemble des points ) avec:
Après quelques expérimentations numériques, nous proposons dans la dernière partie une méthode d'interpolation de Fourier rapide d'un schéma entrelacé vers un schéma standard. Ceci permet de montrer numériquement un exemple simulé en tomographie 2D que la suppression d'une moitié des points d'échantillonnage (supposés recouverts de marqueurs de calibrage) n'enlève essentiellement pas d'information au sens de Shannon et que cette information peut-être retrouvée efficacement par des techniques de transformée de Fourier rapide.