An algorithm for correcting X-ray image distortions caused by central projection A.V. Ustinov2 , N.Yu. Ilyasova1,2, N.S. Demin1 1 Samara National Research University, 34 Moskovskoye shosse, 443086, Samara, Russia 2 Image Processing Systems Unstitute - Branch of the Federal Scientific Research Centre "Crystallography and Photonics" of Russian Academy of Sciences, 151 Molodogvardeyskaya street, 443001, Samara, Russia Abstract We propose an algorithm for correcting X-ray image distortions caused by central projection. Relationships between coordinates of X-ray image points obtained using parallel and central projection are derived. We describe two correction techniques using which the original image can be made similar to the one based on parallel projection. In this way, the three-dimensional heart vessel model can be essentially simplified and diagnostic parameters can be assessed with higher accuracy, resulting in a more accurate early disease diagnosis. Key words: X-ray image; distortion correction; central projection 1. Introduction Roentgenology is an extensively used and dynamic branch of medicine that uses X-ray imagery for the diagnosis of a variety of diseases. By way of illustration, X-ray angiography is utilized for diagnosis of cardiovascular diseases [1]. As a rule, the diagnosis is made based on visual assessment of angiograms, however, its accuracy essentially depends on the projection angle. A three-dimensional model of heart vessels serves to visualize three-dimensional geometric and topological information, thus enabling the diagnosis to be made with higher accuracy [2-6]. The initial data is a sequence of DICOM frames [7, 8]. The projection procedure is affected by a number of technical limitations, resulting in the imagery characterized by a variety of distortions. Various techniques for distortion correction were described in Ref. [9,10,13]. The X-ray imagery is obtained using specialized clinical equipment, such as an operating-room X- ray unit C-ARM. Examples of such equipment are illustrated in Figure 1. The unit is composed of an X-ray source and receiver connected by an arc-shaped holder freely moving on a support. With such a design, the X-ray camera has two degrees of freedom. The camera can also move relative to the holder in the longitudinal direction, enabling the image to be scaled. The spatial position of the camera is described by two angles: primary and secondary angles of the camera position. a) b) c) d) Fig. 1. Specialized clinical equipment: a) AXIOM Multista, b) AXIOM Artis BC, c) primary angle, d) secondary angle. The primary angle (denoted as α) is provided by rotating the camera holder and the support as a whole relative to the mount beam. The secondary angle β (analogous to the geographic latitude) is provided by sliding the arc-shaped holder on the support guide, with the camera and X-ray source moving on a circular arc. The imaging is done by a diverging X-ray, with the divergence angle defined by technical characteristics of the scanning device and usually found within 10-12°. As a result, the X- 3rd International conference “Information Technology and Nanotechnology 2017” 10 Computer Optics and Nanophotonics / A.V. Ustinov , N.Yu. Ilyasova, N.S. Demin ray image is observed in the central projection. Based on such projections, it is not possible to reconstruct a model of the original object without additional information concerning imaging conditions, which is not available. At the same time, the original image can be reconstructed from parallel projections without use of any additional information [7,8]. Our idea is that effects caused by central projection can partially be compensated for by making the X-ray image look as if it were built using parallel projections. By correcting X-ray image distortions caused by central projection, the three-dimensional heart vessel model can be essentially simplified [12] and diagnostic parameters can be assessed with higher accuracy [5], resulting in a more accurate early disease diagnosis. 2. Mathematical model Because the refractive index for X-rays for all substances is very close to one, it is impossible to make a collimator that converts a divergent beam from a source close to a point beam into a parallel beam. This leads to an object image deformation even at the planar object. At straight rays falling (the primary and secondary angles are zero), the distortion reduces to a change in scale - the image remains similar to the image obtained by parallel projection. At inclined falling an additional distortion appears: the scale becomes different on image area, which leads to a violation of similarity (a change in the proportions of the object in the case of oblique incidence also occurs in parallel projection). We propose an algorithm for compensating of central projection influence (divergent rays). It should be noted that full compensation is possible only for a planar object, because otherwise it is principally impossible to indicate precise values of some parameters. Nevertheless, the achieved compensation is enough for more precision of three-dimensional trace process described in Ref. [7, 12]. In order that to compensate the central projection influence we are need to get a relation of image point coordinates at parallel and central projection. For this purpose a calculation of world and planar coordinates of projection point is required. Let a point S is a light source and distance OS  H is given. A direction to the light source is defined as OS  nH . Point A  ( x0 ; y0 ; h) is an intersection of projection plane with straight-line having a directing vector SA, where h is distance from image plane to the object at zero angles. Now we obtain the projection of point A  ( x; y; z ). Planar coordinates of point are u  OAu  xux  yu y  zuz , v  OAv  xvx  yvy  zvz , where u , v are basic vectors of planar coordinates system. If the plane is defined by an equation Ax  By  Cz  D  0 and straight-line is defined by formulas x  x0  lt ; y  y0  mt ; z  z 0  nt , then coordinates of point A are obtained at substitution of value t derived from equation ( Al  Bm  Cn )t  Ax0  By0  Cz 0  D  0 . In our case: D  0; ( A; B ; C )  ( n x ; n y ; n z ); ( x 0 ; y 0 ; z 0 )  ( x0 ; y 0 ; h ); ( l ; m ; n )  AS  nH  ( x0 ; y 0 ; h )  ( Hn x  x0 ; Hn y  y 0 ; Hn z  z 0 ). At the parallel projection: (l ; m; n )  n  ( A; B; C ) . Taking in account these values we get the projection point coordinates: x0 ( B 2  C 2 )  y0 AB  hAC x , 1  ( Ax0  By0  Ch) / H  x0 AB  y0 ( A2  C 2 )  hBC y , (1) 1  ( Ax0  By0  Ch) / H  x0 AC  y0 BC  h( A2  B 2 ) z . 1  ( Ax0  By0  Ch) / H Fig. 2. Effect of divergent beam. Left - straight falling (α=β=0), right - inclined (α=45°, β=0).The second frame has larger crosshairs on left part in image than in right part. 3rd International conference “Information Technology and Nanotechnology 2017” 11 Computer Optics and Nanophotonics / A.V. Ustinov , N.Yu. Ilyasova, N.S. Demin Numerators are projection point coordinates at the parallel projection and L  Ax0  By0  Ch in denominators is signed distance from point A to the projection plane. Taking in account formulas for planar coordinates we get: u  u par 1  L / H  , v  v par 1  L / H  . In some cases value of H we can get from DICOM-file header, else we can evaluate it with aid of mira image - special etalon image (grid in fig.2). In item 3 an algorithm invented for evaluation of H is described. About value of h: if the object is not planar then h can not be obtain in principle, since a reconstruction of three-dimensional object is need to this. For approximate solution of correction task we propose three approaches: 1) If object length in OZ axis direction compare with H is small and object is near to receiver then we use h=0. 2) If object length is small but object is not near to receiver there we use some average value, for example, a heart size. 3) Assuming that distortions are not large (three-dimensional tracing process is not failed) we build three-dimensional tree ignoring the central projection influence. That h equal to z-coordinate of corresponding tree point. After correction the three- dimensional tree building is repeated. But this variant increases calculation volume. 3. Determination of distance from the source to the receiver with usage of mira For the determination of distance from the source to receiver with usage of mira, we preliminarily consider an auxiliary case. Let we have straight falling rays (at neglecting distortion of central projection) and a task is planar: source, object and coordinates origin lie in one plane perpendicular to the plane of the receiver. The distance from point source to receiver plane (the origin) is H, the distance from point object to the origin is h. Then a following relation holds: xc   H H  h  x p , it is connects coordinates of the object image on the receiver plane with the parallel projection and the central one. Now we consider our case. The task is still planar - the source and the non-point object passing through the origin lie in one plane perpendicular to the plane of the receiver. Let he object is segment of length 2a lying horizontally in this plane and the origin is its center, i.e. h=0 for segment center. A left point will be the first, and a right point will be the second. Rays fall bias at angle φ to vertical (the angle is counter on clockwise). The distance of the point source to the origin is H. Under this scheme of observation the central point of segment-object retains its zero coordinate on segment-image, but it will not be a center of segment-image. For the parallel projection, image stays symmetrical: if we use a length (instead of coordinate with sign) then x p1  x p 2  a cos  . Because of the segment is not lie in in receiver plane, value h  0 for its ending points. A module of this value is equal for the first and the second points: h1  h2  a sin  . But for the first (left) point it has minus sign. For absolute values we have final result: H xc1  a cos  ; H  a sin  (2) H xc 2  a cos  . H  a sin  This means that if we know angle φ and the lengths of segments measured on the image which are equal on the object, then we can calculate distance H and the length of segment by formulas: 2 xc1 xc 2 a ; ( xc1  xc 2 ) cos  (3) 2 xc1 xc 2 H  tg . xc 2  xc1 A really used mira image contains many segments. Since mira grid on object is rectangular then it is desirable to use of images photographed when one of camera orientation angle is zero. This give two advantages: angle φ is calculated easily (it is equal to modulus of the second orientation angle) and lines of equal scale stay be straight and parallel to image sides. If β=0, these lines are vertical, if α=0, ones are horizontal. For definiteness, we shall consider case of β=0. At this condition we have the following algorithm. Step 1. We do search of mira grid crosshair nearest to the image center. Step 2. We do step back on some squares to left (not necessary on line containing center). A difference of horizontal coordinates of central and new crosshair is value xc1 . 3rd International conference “Information Technology and Nanotechnology 2017” 12 Computer Optics and Nanophotonics / A.V. Ustinov , N.Yu. Ilyasova, N.S. Demin Step 3. We do step back on some squares to right (not necessary on line containing center). A difference of horizontal coordinates of new and central crosshair is value xc 2 . Step 4. The sought distance H is calculated by the second formula (3). If it is negative then modulus is used - minus sign means that right segment is shorter of left one that takes place of rays falling left from vertical. At obtaining of formulas (2) we assumed rays falling right from vertical. 4. Correction distortion of central projection with cycle on source image The first method correction distortion of central projection consists in following: for points on source image we search corresponding points on corrected image. If x0 , y0 are found then correction is implemented by formulas:  Ax  By0  Ch  u par  u  1  0   H  .  Ax0  By0  Ch  v par  v  1    H  Solution of equation set for values x0 , y0 is following: a (u  b1h )  a2 (v  b2 h) x0  4 a1a4  a2 a3 . (4) a (v  b2 h )  a3 (u  b1 h ) y0  1 a1a4  a2 a3 Knowing expressions for vectors dependence on primary and secondary angles n  ( sin  ;sin  cos  ;cos  cos  ) u  (cos  ;sin  sin  ;cos  sin  ) v  (0;cos  ;  sin  ) we get expressions for coefficients a1 , a2 , a3 , a4 , b1 , b2 : a1    u H  sin   cos  , a2   u H  sin  cos   sin  sin  , a3    v H  sin  , (5) a4   v H  sin  cos   cos  , b1    u H  cos  cos   cos  sin  , b2    u H  cos  cos   sin  . But if we want to form the correct image this method is undesirable: because of discretization part of output image points (in area of stretching) will not have a prototype. There is a grid from points without the prototype on fig. 3 that usually is not suited for us. Fig. 3. Correction with cycle on source image. Left – source image, right – corrected one. A white grid is points without prototype. 3rd International conference “Information Technology and Nanotechnology 2017” 13 Computer Optics and Nanophotonics / A.V. Ustinov , N.Yu. Ilyasova, N.S. Demin 5. Correction distortion of central projection with cycle on corrected image For avoidance of empty points we employ the second method correction effects of central projection. Here the standard approach of image spatial transformation (for example, rotation or reflecting in non-planar mirror) is used - a cycle employs on output image and the prototype of current point is found on the source image. If x0 , y0 are found then the prototype point has coordinates:  Ax  By0  Ch  u  u par /  1  0 ,  H   Ax  By0  Ch  v  v par /  1  0 .  H  Coordinates of point x0 , y0 is calculated by formulas: a4 (u par  b1h)  a2 (v par  b2 h) x0  , a1 a4  a2 a3 (6) a1 (v par  b2 h)  a3 (u par  b1h) y0  . a1a4  a2 a3 Coefficients are equal to a1  cos  ; a2  sin  sin  ; a3  0; a4  cos  ; b1   cos  sin  ; b2  sin  , where  ,  are primary and secondary angles of camera rotation, and h is distance from the object to the projection plane. After a substitution of these values we find simpler formulas: a4 (u par  b1h)  a2 (v par  b2 h) x0  , a1a4 (7) v par  b2 h y0  . a4 On fig. 4 we see that the grid from points without prototype is absent, i.e. the correction is more precise. Fig. 4. Correction with cycle on corrected image. Left – source image, right – corrected one. 6. Conclusion In this work, we propose an algorithm enabling X-ray image distortions caused by central projection to be corrected for. Considering that it is not possible to reconstruct the original object from a three-dimensional heart vessel model without additional information concerning the imaging conditions, which is not available in most cases, the process gets more complicated. At the same time, the reconstruction based on parallel projections requires no additional information. We propose a relation of image point coordinates at parallel and central projection. We describe two correction techniques using which the original image can be made similar to the one based on parallel projection. In this way, the three-dimensional heart vessel model can be essentially simplified and diagnostic parameters can be assessed with higher accuracy, resulting in a more accurate early disease diagnosis. 3rd International conference “Information Technology and Nanotechnology 2017” 14 Computer Optics and Nanophotonics / A.V. Ustinov , N.Yu. Ilyasova, N.S. Demin Acknowledgements The work was carried out with the partial support of the Ministry of Education and Science of the Russian Federation within the framework of the activities of the Program for Enhancing the Competitiveness of the SSAU among the world's leading scientific and educational centers for 2013-2020; Grants of the Russian Foundation for Basic Research No. 14-07- 97040, No. 15-29- 03823, No. 15-29- 07077, No. 16-57-48006, No. 16-41-630761; Program № 6 of fundamental research Department of Nanotechnologies and Information Technologies of the Russian Academy of Sciences "Bioinformatics, modern information technologies and mathematical methods in medicine" 2016 -2017. References [1] Chandra T, Pukenas B, Mohan S, Melhem E. Contrast-Enhanced Magnetic Resonance Angiography. Magnetic Resonance Imaging Clinics of North America 2012; 20(4): 687–698. [2] Ilyasova N. Computer Systems for Geometrical Analysis of Blood Vessels Diagnostic Images. Optical Memory and Neural Networks (Information Optics) 2014; 23(4): 278–286. [3] Ilyasova N. Methods to evaluate the three-dimensional features of blood vessels. Optical Memory and Neural Networks (Information Optics) 2015; 24(1): 36–47. [4] Ilyasova N. Evaluation of Geometric Characteristics of the Spatial Structure of Vessels. Pattern Recognition and Image Analysis 2015; 25(4): 621–625. [5] Ilyasova NYu. Estimating the geometric features of a 3D vascular structure. Computer Optics 2014; 38(3): 529–538. [6] Ilyasova NYu. Methods for digital analysis of human vascular system. Literature review. Computer Optics 2013; 37(4): 511–535. [7] Ilyasova NYu, Kupriyanov AV, Khramov AG. Information technologies of image analysis in medical diagnostics. Radio and communication, 2012. [8] Soifer VA. Computer Image Processing. Part II: Methods and algorithms: Appendix A2. Biomedi-cal Images Processing. VDM Verlag, 2009. [9] Hong C, Lee D-H, Han BS. Characteristics of geometric distortion correction with increasing field-of-view in open-configuration MRI. Magnetic Resonance Imaging 2014; 32(6): 786–790. [10] Cardoso PL, Dymerska B, Bachratá B, Fischmeister FPhS, Mahr N, Matt E, Trattnig S, Beisteiner R, Robinson SD. The clinical relevance of distortion correction in presurgical fMRI.NeuroImage 2016. [11] Menga C, Zhang J, Zhou F, Wang T. New method for geometric calibration and distortion correction of conventional C-arm. Computers in Biology and Medicine 2014; 52: 49–56. [12] Ilyasova NYu, Kazanskiy NL, Korepanov AO, Kupriyanov AV, Ustinov AV, Khramov AG, Computer technology for reconstructing the 3D structure of coronary arteries from angiographic projections. Computer Optics 2009; 33(3): 281–318. [13] Moravec J, Hub M. Automatic correction of barrel distorted images using a cascaded evolutionary estimator. Information Sciences 2016; 366: 70–98. 3rd International conference “Information Technology and Nanotechnology 2017” 15