Abstract:
A method of deriving a quantitative measure of a degree of calcification of a blood vessel such as an aorta by processing an image such as an X-ray image of at least a part of the blood vessel containing said calcification comprises: taking a starting set of digital data representative of an image of at least part of a blood vessel containing a calcification set against a background;
       estimating the boundary of the calcification;   using inpainting to replace digital data in said starting set representing the calcification with data extrapolating the boundary of the background to extend over the area of calcification, and so generating an inpainted set of digital data; and   computing the difference between the starting set of digital data and the inpainted set of digital data to obtain a quantitative measure of the degree of calcification of the blood vessel.

Description:
BACKGROUND OF THE INVENTION 
     The present invention relates to a method of deriving a quantitative measure of a degree of calcification of an aorta. 
     Inpainting is a technique that originates from retouching paintings where one wants to recreate lost or damaged structures in a legible way. Digital inpainting uses spatial or frequency information to restore partially damaged/removed images. 
     Various inpainting techniques are known that enable image restoration, in particular for photographs, videos and films. 
     It is known to detect and inpaint small regions in mammograms that possibly define a micro calcification to enable detection of calcified regions. Subsequently, features such as the average and the standard deviation of intensity values are extracted from both the pre- and the post-inpainting regions. A classifier is trained to distinguish between true micro calcifications and false positives based on the extracted features. The comparison between a region and its inpainting is used to enable detection. Thus a binary decision of whether a region is abnormal, i.e. different from its surroundings, is made. 
     There are, however, no methods available at present that use inpainting to give more than a basic indication of the presence of a calcification. In the present invention, it has been realised that such a method may be useful in the diagnosis of various diseases, for example, atherosclerosis. 
     Atherosclerosis is a process in which deposits of fatty substances, cholesterol, cellular waste products, calcium and other products build up in the inner lining of an artery. 
     Previous known methods of assessing atherosclerotic plaque include the assessment of the calcification index on lateral 2-D x-rays. In this method, aortic calcification is assessed at each vertebral segment. Calcific deposits are regarded as present if the densities are visible in an area parallel to the lumbar spine and anterior to the lower part of the spine. As the abdominal aorta in the L1-L4 region of the spine is occasionally situated lateral to the spine, aortic densities of the lumbar region can sometimes overlap the vertebrae. Densities overlapping the vertebrae are scored present only if they extend from or form a clear pattern with those of the lower part of the aorta. Calcific deposits in the abdominal aorta adjacent to each lumbar vertebra are assessed separately for the posterior and anterior wall of the aorta using the midpoint of the intervertebral space above and below the vertebrae as boundaries. 
     The lesions are graded according to specific categories: Lesions with small scattered deposits that extend along less than ⅓ of the longitudinal wall of the aorta; lesions where ⅓ or more, but less than ⅔ of the longitudinal wall is calcified and lesions where ⅔ or more of the longitudinal wall is calcified. 
     SUMMARY 
     The present invention seeks to provide a more accurate measure of the degree of calcification of the aorta or more generally any blood vessel using inpainting techniques. There is accordingly provided a method of deriving a quantitative measure of a degree of calcification of a blood vessel by processing an image of at least a part of the blood vessel containing said calcification, which method comprises: 
     taking a starting set of digital data representative of an image of at least part of a blood vessel containing an area of calcification, said area of calcification being set against a background area; 
     using an inpainting method to replace digital data in said starting set of digital data representing the calcification with data extrapolating the background to extend over the area of calcification, and so generating an inpainted set of digital data; and 
     computing the difference between the starting set of digital data and the inpainted set of digital data to obtain a result representative of the degree of calcification of the blood vessel. 
     The method may comprise the preliminary steps of estimating a boundary of the area of calcification, wherein the step of inpainting comprises replacing digital data within the boundary with data extrapolating the background outside the boundary to extend over the area within the boundary. Generally inpainting propagates data values in said boundary region into replacement data values for the selected area of calcification. Inpainting methods are described in numerous previous publications including US2004/0164996, U.S. Pat. No. 6,587,592, US2003/0012453 and US 2004/0161153, all of which are hereby incorporated by reference. 
     The method of the invention may further comprise: 
     calculating a signal to noise ratio along the boundary of the calcification between the inpainted set of digital data and the starting set of digital data;. 
     increasing the area within the boundary, and repeating the step of calculating the signal to noise ratio along the boundary; and 
     repeating the steps of increasing the area within the boundary and then calculating the signal to noise ratio between the inpainted set of digital data and the starting set of digital data along the boundary until a maximum signal to noise ratio is reached. 
     Such methods may further comprise calculating the average intensity difference between the inpainted set of digital data and the starting set of digital data, wherein the difference is-representative of the degree of calcification of the blood vessel. 
     Methods of inpainting suitable for use include all those currently known in the art or which may be developed hereafter including without limitation TV inpainting harmonic inpainting and average inpainting. 
     Such methods may be applied to any blood vessel but preferably to an artery, particularly the aorta. 
     The invention may be embodied in a preprogrammed computational device or an instruction set for such a device adapted to carry out the above method. 
    
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
       Embodiments of the present invention will hereinafter be described, by way of example, with reference to the accompanying drawings, in which: 
         FIG. 1  shows a spinal column with the lumbar region indicated; 
         FIG. 2  shows an illustration of divergence at P resulting in varying densities of particles as described in the present invention; 
         FIG. 3  shows an illustration of the algorithm of the present invention run on a test image; 
         FIG. 4A  shows the 20 sample area templates of the test image and  FIG. 4B  the regression lines in standard deviation against area in number of pixels for the three different inpainting methods; 
         FIGS. 5A and 5B  illustrate the bias introduced in the inpainting procedure, due to minor calcific deposits that were overlooked in the manual annotation process; 
         FIG. 6  is a table showing the SNR values and sum of the average differences between the original and inpainted images when the various inpainting methods described in the present invention are used on 10 test images; 
         FIGS. 7A and 7B  show comparative graphs of the difference inpainting methods of the present invention measured in standard procedure for detecting calcification of the aorta against the average intensity difference of individual images. 
     
    
    
     DESCRIPTION OF PREFERRED EMBODIMENTS 
     To generate a quantitative measure of the degree of calcification in the aorta the following steps can be taken: 
     1. Obtain lateral 2-D x-ray images of the spine. 
     2. Estimate the location of the aorta. 
     3. Detect calcified areas. 
     4. Obtain a quantitative measure of the degree of calcification. 
     5. Output the result. 
     The present invention is primarily concerned with Step 4, the actual process involved in obtaining a quantitative measure of the degree of calcification. 
     Methods of obtaining x-ray images are well known. This step will therefore not be described herein. Step 2 of estimating the location of the aorta is described in SPIE Medical Imaging: Image Processing 2004, L. A. Conrad-Hansen et al. “Prediction of the location of the lumbar aorta using the first four lumbar vertebrae as a predictor”, Edition 5370 (2), pgs 1271-1281. 
     In the described embodiment, lateral 2-D x-rays are used to provide a cost-efficient way of delivering reliable and easily reproducible data which may be used to extract information regarding the degree of calcification of the aorta. 
     The region of interest for which x-ray images are required is the lumbar region, denoted by L1-L4 in  FIG. 1 . As a result of the proximity of the thoracic region of the aorta to the rib cage, x-ray images taken from the thoracic region are generally contaminated with some noise. Therefore, the lumbar region is used. Furthermore, bifurcation of the aorta with the first parts of the common iliac arteries most commonly occurs at L4, which restricts the region of interest to the first four lumbar vertebrae. The bifurcation causes surges in the blood stream, which in turn contribute to the calcification process, so that calcific deposits primarily occur at L4 followed by L3. Therefore, if calcific deposits are found in the area of the aorta corresponding to L2, then the aortic region corresponding to the lower lying vertebrae is almost always even more calcificied. 
     Inpainting is a known technique that originates from retouching paintings to recreate lost or damaged structures in a legible way. The term “Inpainting” is derived from fusing the words “Image” and “Interpolation” together. Three different variations of inpainting are herein described, TV inpainting, harmonic inpainting and average inpainting. TV inpainting is used for its ability to preserve structures to some extent, whereas harmonic inpainting is used as it provides much smoother solutions. The general Bayes&#39; formulation of inpainting can be expressed as: 
               p   ⁡     (     u   ❘     u   0       )       =           p   ⁡     (       u   0     ❘   u     )       ⁢     p   ⁡     (   u   )           p   ⁡     (     u   0     )         ∝       p   ⁡     (       u   0     ❘   u     )       ⁢     p   ⁡     (   u   )                 
since u 0  is known.
 
     The model term p(u|u 0 ) models the data formation process. It can be assumed that the observed image u 0  is obtained from a clean image u corrupted by additive zero-mean Gaussian noise spatially uncorrelated and stationary, with standard deviation σ outside the hole Ω, and that the data inside Ω has been totally lost. Therefore the model of acquisition degradation is given by 
               p   ⁡     (       u   0     ❘   u     )       =       c   1     ⁢     ⅇ     -       ∑     x   ∈     R   ⁢   \   ⁢   Ω         ⁢         (       u   ⁡     (   x   )       -       u   0     ⁡     (   x   )         )     2       2   ⁢     σ   2                       
where C1 is the normalising constant. The prior term p(u) is usually difficult to find since it models the probability that a given array of pixel values represents a meaningful image. In the absence of texture, some form of smoothness for images is assumed, i.e. the variation of pixel values around a given pixel location should be small. Therefore, a discrete gradient operator, ∇u(x) encoding the local variations of the image around a given pixel location can be introduced.
 
     The basics of TV inpainting will now be described. In order to enforce a reasonable smoothness, a Laplace distribution is assumed on all values, and the following prior is obtained: 
               p   ⁡     (   u   )       =       c   2     ⁢     ⅇ     -       ∑     x   ∈   R       ⁢            ∇     u   ⁡     (   x   )              μ                   
where c 2  is a normalising constant and μ√{square root over (2)} the standard deviation. Using Maximum A Posteriori (MAP) criterion, taking—log( ) of each term, and then shifting to a continuous formulation, μ can be found by minimizing the following energy expression:
 
     
       
         
           
             
               E 
               ⁡ 
               
                 ( 
                 u 
                 ) 
               
             
             = 
             
               
                 
                   1 
                   
                     2 
                     ⁢ 
                     
                       σ 
                       2 
                     
                   
                 
                 ⁢ 
                 
                   
                     ∫ 
                     
                       R 
                       / 
                       Ω 
                     
                   
                   ⁢ 
                   
                     
                       
                         ( 
                         
                           u 
                           - 
                           
                             u 
                             0 
                           
                         
                         ) 
                       
                       2 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ⅆ 
                       x 
                     
                   
                 
               
               + 
               
                 
                   1 
                   μ 
                 
                 ⁢ 
                 
                   
                     ∫ 
                     R 
                   
                   ⁢ 
                   
                     
                        
                       
                         ∇ 
                         u 
                       
                        
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ⅆ 
                       x 
                     
                   
                 
               
             
           
         
       
     
     To make the integration boundaries equal, if χ denotes the function χ(x)=0 if x εΩ, χ(x)=1 otherwise, and 
     
       
         
           
             λ 
             = 
             
               
                 
                   σ 
                   2 
                 
                 μ 
               
               . 
             
           
         
       
     
     The resulting energy expression can be written as
 
 E ( u )=∫ R χ( u−u   0 ) 2   dx+λ∫   R   |∇u|dx, 
 
     Having surpassed the problem with the integration boundaries, it is feasible to use the calculus of variations, so that a minimiser u for the energy is given by the following Euler-Lagrange equation 
                 ∂   E       ∂   u       =         ∂   F       ∂   u       -       ∂     ∂   x       ⁢       ∂   F         ∂   x     ⁢     ∂     u   x             -       ∂     ∂   y       ⁢       ∂   F         ∂   y     ⁢     ∂     u   y                     
where F(u,u x ,u y )=χ(u−u 0 ) 2 +λ√{square root over (u x   2 +u y   2 )}, so that:
 
     
       
         
           
             
               
                 ∂ 
                 E 
               
               
                 ∂ 
                 u 
               
             
             = 
             
               
                 
                   2 
                   ⁢ 
                   
                     χ 
                     ⁡ 
                     
                       ( 
                       
                         u 
                         - 
                         
                           u 
                           0 
                         
                       
                       ) 
                     
                   
                 
                 - 
                 
                   λ 
                   ⁢ 
                   
                     
                       
                         
                           u 
                           xx 
                         
                         ⁢ 
                         
                           u 
                           y 
                           2 
                         
                       
                       - 
                       
                         2 
                         ⁢ 
                         
                           u 
                           x 
                         
                         ⁢ 
                         
                           u 
                           xy 
                         
                       
                       + 
                       
                         
                           u 
                           yy 
                         
                         ⁢ 
                         
                           u 
                           x 
                           2 
                         
                       
                     
                     
                       
                         ( 
                         
                           
                             u 
                             x 
                             2 
                           
                           + 
                           
                             u 
                             y 
                             2 
                           
                         
                         ) 
                       
                       
                         3 
                         2 
                       
                     
                   
                 
               
               = 
               0 
             
           
         
       
     
     This scheme can be made discrete using gradient descent. However, it cannot be proven to yield a unique solution. A more elegant solution can be found by introducing the divergence of a vector field shown in  FIGS. 2A ,  2 B and  2 C. A vector field on R 2  is a map that attaches at every point P(x,y) (shown in figure x) a vector 
     
       
         
           
             
               
                 υ 
                 _ 
               
               ⁡ 
               
                 ( 
                 P 
                 ) 
               
             
             = 
             
               
                 υ 
                 ⁡ 
                 
                   ( 
                   
                     x 
                     , 
                     y 
                   
                   ) 
                 
               
               = 
               
                 ( 
                 
                   
                     
                       
                         
                           υ 
                           1 
                         
                         ⁡ 
                         
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                     
                   
                   
                     
                       
                         
                           υ 
                           2 
                         
                         ⁡ 
                         
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                     
                   
                 
                 ) 
               
             
           
         
       
     
     The divergence of  υ  is given by: 
               div   ⁢     υ   _       =       ∇     υ   _       =         ∂     υ   1         ∂   x       +       ∂     υ   2         ∂   y                 
and expresses at a given point P, the local behaviour of the field.
 
     As shown in  FIG. 2A , where there is positive divergence of P, the vector field diverges in the neighbourhood of P and the density of particles in the field decreases. Similarly, as shown in  FIG. 2B , where there is negative divergence, the vector field converges in the neighbourhood of P and the density of particles increases. Finally, where there is no divergence at P, as shown in  FIG. 2C , the density of particles remains constant. 
     A natural divergence term therefore appears as: 
                   ∂     ∂   x       ⁢       ∂   F       ∂     u   x           +       ∂     ∂   y       ⁢       ∂   F       ∂     u   y             =     div   ⁡     (             ∂   F       ∂     u   x                     ∂   F       ∂     u   y               )             
which in the case of TV inpainting corresponds to the term:
 
                 div   (             u   x           u   x   2     +     u   y   2                       u   y           u   x   2     +     u   y   2                 )     =       div   (       ∇   u            ∇   u            )     =   0       ,     λ   =       σ   2     μ             
so that the following Euler-Lagrange Partial Differential Equation (PDE) can be solved:
 
     
       
         
           
             
               
                 
                   χ 
                   ⁡ 
                   
                     ( 
                     
                       u 
                       - 
                       
                         u 
                         0 
                       
                     
                     ) 
                   
                 
                 - 
                 
                   λ 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     div 
                     ( 
                     
                       
                         ∇ 
                         u 
                       
                       
                          
                         
                           ∇ 
                           u 
                         
                          
                       
                     
                     ) 
                   
                 
               
               = 
               0 
             
             , 
             
               λ 
               = 
               
                 
                   σ 
                   2 
                 
                 μ 
               
             
           
         
       
     
     This is a non linear elliptic equation, not well defined when |∇u|=0. This problem is usually fixed by replacing |∇u| by φ(|∇u| 2 ) where φ(s 2 )=√{square root over (s 2 +ε 2 )}, ε&gt;0 small. The new equation to solve is then:
 
χ( u−u   0 )−2λdiv(φ′(|∇ u|   2 )∇ u )=0
 
     In order to cope with the non linearity, a fixed-point strategy is used for estimating the φ′(|∇u| 2 ). This leads to the following resolution algorithm: 
     
       
         
               
               
             
           
               
                   
                   
               
             
             
               
                   
                 1. Set β 0  (x) ≡ 1 
               
               
                   
                 2. For i = 0...N do 
               
               
                   
                  * Compute the solution u i  of the linear PDE 
               
               
                   
                   χ(u i  − u 0 ) − 2λdiv(β i ∇u i ) = 0   (*) 
               
               
                   
                  * Set β i+1  = φ ′ (|∇u i | 2 ) 
               
               
                   
                 3. Output u N   
               
               
                   
                   
               
             
          
         
       
     
     Each linear equation (*) is made discrete and solved directly as follows. Given a pixel location s ε R, denote by v(s) its 4-neighbourhood. If r ε v(s), then 
               r   +   s     2         
will denote the mid-point between r and s. The discretization of the divergence is performed using central differences with half grid steps. Then, skipping the superscripts i, at pixel location s equation (*) is made discrete as:
 
     
       
         
           
             
               
                 
                   χ 
                   s 
                 
                 ⁡ 
                 
                   ( 
                   
                     
                       u 
                       s 
                     
                     - 
                     
                       u 
                       
                         0 
                         ⁢ 
                         s 
                       
                     
                   
                   ) 
                 
               
               - 
               
                 2 
                 ⁢ 
                 λ 
                 ⁢ 
                 
                   
                     ∑ 
                     
                       r 
                       ∈ 
                       
                         v 
                         ⁡ 
                         
                           ( 
                           s 
                           ) 
                         
                       
                     
                     
                         
                     
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     β 
                     ⁢ 
                     
                       
                         r 
                         + 
                         s 
                       
                       2 
                     
                     ⁢ 
                     
                       ( 
                       
                         
                           u 
                           r 
                         
                         - 
                         
                           u 
                           s 
                         
                       
                       ) 
                     
                   
                 
               
             
             = 
             0 
           
         
       
     
     with boundary conditions incorporated as u r :=u s  whenever r falls out of R. The resulting system (for all s ε R) is solved by a Gauss-Seidel iterative scheme, which, at a given iteration reads 
     For all s ε R do 
               u   s     ←           χ   s     ⁢     u     0   ⁢   s         +     2   ⁢   λ   ⁢       ∑     r   ∈     v   ⁡     (   s   )                   ⁢           ⁢       β       r   +   s     2       ⁢     u   r                 χ   s     +     2   ⁢   λ   ⁢       ∑     r   ∈     v   ⁡     (   s   )           ⁢           ⁢     β       r   +   s     2                     
where the symbol←means “replace the old value by the new one”. The iterative process is stopped after a convergence criteria, which is reached when the change in u s  falls beneath a certain threshold, has been fulfilled or after a fixed number of iterations. Note that it can be seen as a local lowpass filtering process—the filter coefficients are positive and their sum is 1—biased toward the original value u 05  when χ s =1.
 
     The coefficient β r+s/2   i  can be made discrete as follows: 
     If s=(j,k),r=(j−1,k), it is required from the definition of β i  to compute |∇u i−1 | j−1/2,k   2 . The x-derivative is approximated by a standard half grid step central difference, centered at 
               j   -     1   2       ,     k   .           
The y-derivative is interpolated using the central differences at (j,k) an at (j−1,k), so that
 
     
       
         
           
             
               
                 
                    
                   
                     ∇ 
                     
                       u 
                       
                         i 
                         - 
                         1 
                       
                     
                   
                    
                 
                 
                   
                     j 
                     - 
                     
                       1 
                       2 
                     
                   
                   , 
                   k 
                 
               
               2 
             
             ≈ 
             
               
                 
                   ( 
                   
                     
                       u 
                       jk 
                     
                     - 
                     
                       u 
                       
                         j 
                         - 
                         
                           1 
                           ⁢ 
                           k 
                         
                       
                     
                   
                   ) 
                 
                 2 
               
               + 
               
                 
                   1 
                   16 
                 
                 ⁢ 
                 
                   
                     ( 
                     
                       
                         u 
                         
                           jk 
                           + 
                           1 
                         
                       
                       - 
                       
                         u 
                         
                           jk 
                           - 
                           1 
                         
                       
                       + 
                       
                         u 
                         
                           j 
                           - 
                           
                             1 
                             ⁢ 
                             k 
                           
                           + 
                           1 
                         
                       
                       - 
                       
                         u 
                         
                           j 
                           - 
                           
                             1 
                             ⁢ 
                             k 
                           
                           - 
                           1 
                         
                       
                     
                     ) 
                   
                   2 
                 
               
             
           
         
       
     
     Although it can be proven that a minimiser ū 0  always exists under reasonable conditions, it might fail to be unique. It is therefore necessary to specify a good starting value for ū 0  for the inpainting algorithm as it allows for faster convergence. Outside the inpainting domain Ω, ū 0  (x, y)=u 0  (x, y) and inside Ω, ū 0  (x, y) is drawn from a Gaussian distribution with mean and variance corresponding to the mean and variance of u 0  outside Ω. 
     For harmonic inpainting, a Gaussian distribution is assumed for the prior probability 
               ??   ⁡     (   ??   )       =       c   2     ⁢     ⅇ     -       ∑     x   ∈   R               ⁢           ⁢              ∇     u   ⁡     (   x   )              2       2   ⁢     μ   2                       
where c 2  again is a normalising constant and μ is the standard deviation. The energy term can then be expressed as
   E ( u )=∫ R χ( u−u   0 ) 2   dx+λ∫   R   |∇u|   2   dx,   
where χdenotes the function χ(x)=0 if x ε Ω, χ(x)=1 otherwise, and
 
     
       
         
           
             λ 
             = 
             
               
                 
                   σ 
                   2 
                 
                 
                   μ 
                   2 
                 
               
               . 
             
           
         
       
     
     A third form of inpainting, average inpainting, constitutes the simplest form of inpainting, where Ω is filled homogeneously with the value S resulting from averaging over the immediate boundary of Ω according to 
               S   =       1   n     ⁢       ∑     i   =   1     n     ⁢           ⁢     t   i           ,         
where n is the number of boundary pixels and t the respective pixel value.
 
     The signal-to-noise ratio (SNR) is calculated for each pixel along the immediate boundary of the calcified area according to 
               SNr   =         I   inp     -     I   orig         std   ⨯   area         ,         
where I inp  denotes the inpainted image, I orig  the original image, std the standard deviation of the estimated noise level, and area the number of pixels in Ω. The mechanism is illustrated in the below algorithm and the results are shown on  FIG. 3  on a test image.
 
     
       
         
               
             
               
               
             
           
               
                   
               
               
                 Algorithm 1 Pixelwise expansion scheme: 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 1.  let N = number of calcified areas in image 
               
               
                   
                 2.  let B = number of pixels in an isocurve 
               
               
                   
                 3.  let p denote a pixel 
               
               
                   
                 4.  For i = 1...N do 
               
               
                   
                    inpaint Ω 
               
               
                   
                    calculate SNR 
               
               
                   
                    while SNR not max do 
               
               
                   
                     calculate outer distance map 
               
               
                   
                     find next isocurve 
               
               
                   
                     For j= 1...B do 
               
               
                   
                      expand Ω by p j   
               
               
                   
                      inpaint and calculate SNR 
               
               
                   
                      if SNR is increased, include p j   
               
               
                   
                    inpaint new Ω 
               
               
                   
                    calculate new SNR 
               
               
                   
                 5.  output SNR. 
               
               
                   
                   
               
             
          
         
       
     
     As shown in  FIG. 3 , the algorithm starts with the original input image with faulty annotation (in this case, manual annotation). The area with the annotated boundary is inpainted, then, in accordance with the algorithm, the SNR ratio for each pixel along the boundary is calculated, and depending on the result of this, the boundary is redefined to include the calcified pixels that might have been undetected. Again, the area within the boundary is inpainted until the maximum signal to noise ratio is output. 
     The std denotes, as explained above, the standard deviation of the estimated noise level, and is computed by estimating the non-calcified background and taking the difference of the estimated and the actual background. The data set for the background estimation contains ninety hand annotated x-rays of the lumbar region with pixel intensities ranging from 0 to 2,048 with displays of various degrees of calcification. 
     In order to estimate how well the individual inpainting techniques perfom on the x-ray images, 20 templates of calcification shape were chosen at random from manual annotations (see  FIG. 4A ). Each template was placed at 80 randomly chosen non-calcified aorta segments in different images and the intensity differences betweeen the inpainted area and the original area were computed. For each area template the distance map was calculated. The standard deviations for the error estimates were calculated for each individual area bounded by the isocurves (i.e. curves that have the same distance to an initial reference curve along the entire boundary) resulting from the distance map.  FIG. 4B  shows regression lines through the calculated standard deviations. The regression lines express the standard deviations in the total intensity difference for the three inpainting methods as a function of areasize. 
     
       
         
               
               
               
               
             
               
               
               
               
               
             
           
               
                   
                   
               
               
                   
                 TV 
                 Harmonic 
                 Average 
               
               
                   
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 std 
                 9.8 
                 11.6 
                 10.3 
               
               
                   
                   
               
             
          
         
       
     
     The above table shows the standard deviations of the pixelwise error for the three inpainting methods. From this it can be seen that TV inpainting produces the least error per pixel whereas harmonic inpainting performs worst. 
     As shown in  FIGS. 5A and 5B , the inpainting procedure is biased by minute calcific deposits just outside the annotated areas. These calcific rests are missed by the physician and became apparent only when zooming into pixel level. Since the inpainting methods rely entirely on boundary information, we had to expand the annotated areas in order to avoid faulty results as much as possible. Thus, the next step is to iterate the inpainting process for each calcified area of an image, so that each area is expanded in order to yield the maximum SNR (see Algorithm 1). 
     Following this, ten images that showed light to moderate calcifications were chosen based on the criteria that if anything decisive on images that are not heavily calcified can be concluded, then more severe calcifications can be easily quantified. 
     As shown in  FIG. 6 , the ten selected images show light to moderate calcifications and contain 1-8 calcified areas. The number area denotes how many calcifications an image contains, Dif total denotes the sum of the average differences between the original and the inpainted image areas and SNR total the sum of the individual signal-to-noise ratios. The final three rows show the respective scores of the standard procedure. 
     The quality of the method described is shown in  FIGS. 7A and 7B , where all the differences for the respective areas for each image are summed up and plotted in sorted order against the respective scoring systems of the standard procedure.  FIG. 7A  shows the result for TV inpainting, and  FIG. 7B  shows the result for harmonic inpainting. From the graphs, it is apparent that each inpainting method offers more possibility for discerning the different stages of plaque development than previous known procedures. 
     It will be appreciated that modifications to or variations of the embodiments described and illustrated may be made within the scope of this application as set out in the appended claims.