Abstract:
Method of determining (a) symmetry line(s) in a radiation image. 
     A method of determining (a) symmetry line(s) in horizontal (vertical) direction in an image that is represented by a matrix of digital pixel values by (i) calculating a symmetry point for each column (row) of the image, and (ii) interconnecting calculated symmetry points.

Description:
The application claims the benefit of U.S. Provisional Application No. 60/058,274 filed Sep. 9, 1997. 
    
    
     FIELD OF THE INVENTION 
     The invention relates to an image processing method to find the symmetry line and orientation angle of a radiation image, such as a medical radiation image. 
     BACKGROUND OF THE INVENTION 
     One of the basic properties of shapes and objects is symmetry. The human body shows symmetry along an axis oriented from top to bottom. 
     Medical images originating from image acquisition modalities such as computed projection radiography show symmetry to a large extent for examination types of body parts which are essentially symmetric. Examples of such examination types are e.g. frontal skull, frontal spine, frontal thorax, frontal abdomen, frontal pelvis. 
     In U.S. Pat. No. 5,533,143 a method is disclosed for determining whether an image is symmetric or asymmetric. Image processing conditions are based on the judgment that the image is a frontal (symmetric) object image or a lateral (asymmetric) image. 
     The disclosed method is based on the assumption that the hypothetical symmetry axis is always fixed and coincides with a vertical center line passing through the center of the image. The method does not discriminate between vertical or horizontal symmetry. 
     By studying a large mixture of all symmetric examination types, it is observed that (1) the symmetry line in most cases does not coincide with the center line of the image, nor is its angle with respect to a perfectly vertical (or perfectly horizontal) line equal to zero and (2) horizontal symmetry occurs equally well as vertical symmetry. 
     The above-mentioned observations limit the accuracy and applicability of the prior art method disclosed in the above-mentioned patent. 
     In U.S. Pat. No. 4,903,310 and U.S. Pat. No. 4,951,201 various histogram based methods are proposed to discriminate between a frontal and a lateral body posture. 
     It is however not possible to decide by applying these methods whether the object is symmetrical and whether the object is positioned horizontally or vertically. 
     The above methods are thus not adequate. 
     Experience has learnt that radiologists prefer to examine radiographic images in standard viewing conditions, i.e. with the image displayed so that a symmetry axis in the image is substantially horizontally or vertically oriented. 
     OBJECT OF THE INVENTION 
     It is an object of the invention to provide a method of determining a symmetry line of a radiation image. 
     It is a further object of the invention to provide a method of determining an orientation angle associated with such a symmetry line. 
     It is still a further object to provide a method by means of which images can be provided having a symmetry axis directed in a predetermined orientation. 
     Further objects of the present will become apparent from the description hereafter. 
     SUMMARY OF THE INVENTION 
     The objects of the present invention are achieved by a method of determining a symmetry line in a radiation image represented by a raster of pixel values comprising the steps of 
     (i) calculating a symmetry point for each column (row) of said image, 
     (ii) interconnecting symmetry points calculated for each column (row). 
     The way in which (a) symmetry point(s) is calculated, is described further on. 
     The method of the present invention can be applied to a digital signal representation of a radiation image such as a medical X-ray image. But it can also be applied to other kinds of radiation images such as images originating from non-destructive testing (NDT) of objects by irradiation of these objects. 
     Collimation and direct exposure areas in the image need to be excluded from consideration since they are not necessarily present in the image in a symmetrical way. In European patent application 97201683.6 an ensemble of methods is disclosed to perform a raw segmentation of a radiation image into collimation areas (foreground), direct exposure areas (background) and diagnostic areas (body parts). 
     The method of the present invention can also be applied to each of the sub-images of a multiply exposed image. 
     The method of the present invention permits to compute a symmetry index for a horizontal and vertical direction in the image, and to classify the image as either horizontally or vertically oriented based on the horizontal and vertical symmetry index referred to further on. 
     The matrix of pixel values representing a radiation image can be obtained by means of a wide variety of image acquisition techniques. 
     In one example of an image acquisition technique, a radiation image which is temporarily stored in a photostimulable phosphor screen is read out by scanning the screen with stimulating radiation, detecting light emitted by the screen upon radiation and converting detected light into a digital signal representation. 
     However, other acquisition techniques may be envisaged. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     Further details of the present invention as well as preferred embodiments hereof will be explained hereafter with reference to the drawings in which 
     FIG. 1 shows a thorax image in the presence of collimation areas and direct exposure areas. The automatically computed vertical symmetry axis is superimposed. 
     FIG. 2 shows the fitted symmetry axis of the thorax image. 
     FIG. 3 shows a pelvis image in the presence of collimation areas and direct exposure areas. The automatically computed horizontal symmetry axis is superimposed. 
     FIG. 4 shows the fitted symmetry axis of the pelvis image. 
    
    
     DETAILED DESCRIPTION 
     Computation of the Symmetry Line 
     The method of the present invention is described in the following for application to the columns j of an image. An analogous description applies to rows i. 
     Given an image represented by a raster of M×N pixel values g(i,j) with M rows ranging from (0, . . . , M−1) and N columns ranging from (0, . . . , N−1), and with larger gray values g representing less exposed image pixels. 
     A symmetry measure is computed for each pixel on row i, assuming that the pixel be the most symmetrical point on the line. 
     The basic consideration is that, in the assumption that the central pixel be the symmetry point for the interval of which it is the center, the gray values at symmetrical positions at either side of the symmetry point will nearly be equal, hence their difference will nearly be zero. 
     The symmetry measure can thus be expressed as a sum S of squared or absolute differences of gray values taken symmetrically around a hypothesized symmetry center c:          S        (       i   =   c     ,   j     )       =       ∑     k   =   1       m        (   i   )                [       g        (       c   +   k     ,   j     )       -     g        (       c   -   k     ,   j     )         ]     2                 c   =   T     ,   …              ,     M   -   T   -   1                            
     in which g(i,j) denotes the gray value at row i and column j, c is the row index of the current hypothesized symmetry center, and T is an offset parameter denoting the number of pixels to skip at the start of the column resp. to omit at the end of the column. 
     In the start and end interval T, it is assumed that too few points are available to compute S reliably. 
     The measure S will be minimal for the most symmetrical point in the column j. The parameter m denotes the maximum number of difference terms at row location c and is given by 
     
       
         m(i)=min(i,M−i−1) 
       
     
     Each sum for a row index i needs proper normalization in order to make correct comparisons possible. A straightforward normalization is to divide S by the number of terms m involved in its computation:                S   _     1          (     c   ,   j     )       =           S        (     c   ,   j     )       m                   c     =   T       ,   …              ,     M   -   T   -   1                            
     A further refinement takes into account particularities of geometry and exposure level of the examination. It has been observed that (1) the symmetry line has higher probability to lie nearby the center of the (sub-)image than nearby the borders of the image, and (2) the symmetry line runs through more exposed body parts. 
     Therefore, the refined normalization first takes the number m of terms into account, and second the gray value of the hypothesized center c itself:                S   _     2          (     c   ,   j     )       =           S        (     c   ,   j     )             g   u          (     c   ,   j     )       ·     m   v                       c     =   T       ,   …              ,     M   -   T   -   1.                            
     In the context of the present invention, u=2 and v=3. The use of exponents larger than 1 has been found to have a stabilizing effect on the location of symmetry points for subsequent rows. Clearly, the measure {overscore (S)} 2  favours potential symmetry centers having sufficiently large gray value g, that is, being more attenuated, and which have more number of terms m involved in their computation. Intuitively stated, symmetry lines nearby the image center line (but not necessarily coinciding with it) and more attenuating body parts such as bone are better candidates for being the actual symmetry line than any other lines. Examination types such as thorax, spine and skull all adhere to these underlying assumptions. In the sequel, it is assumed that the refined normalization is used, i.e. {overscore (S)}={overscore (S)} 2 . 
     Finally, the symmetry center C h (j) of column j is assigned to row c for which {overscore (S)}(c,j) is minimal: 
     
       
         C h (j)=min({overscore (S)}(c,j), c=T, . . . , M−T−1) 
       
     
     The sequence 
     
       
         C h (j) j=0, . . . , N−1 
       
     
     is defined as the horizontal symmetry line of the image. 
     An analog description holds for deriving a vertical symmetry line C v (i). Similarly, a diagonal symmetry line can be derived for objects having a principal axis along the main diagonal or secondary diagonal of the image. 
     In general, for linear symmetry of two-dimensional objects, i.e. symmetry expressed by a symmetry line running along a one-dimensional curve, at any point along the symmetry curve, the difference accumulation measure {overscore (S)} reaches a minimum along a line locally perpendicular to the symmetry curve. 
     Hence, for symmetry lines slightly oblique with respect to the perfectly horizontal resp. vertical direction, the method presented here yields an approximation to the ideal algorithm in which the measure {overscore (S)} is computed along several trajectories through the hypothesized symmetry center, one trajectory of which is perpendicular to the yet unknown symmetry curve and thus having the smaller value of {overscore (S)} associated with it. 
     The computation of a symmetry point becomes unreliable when the gray value profile is constant. Indeed, according to definition S 1  of the normalized sum, every point along a constant profile is an equally potential symmetry center since the symmetry measure is constant zero everywhere. According to definition S 2 , the center point of the line will a priori be symmetry point. Therefore, a quality measure is needed to express the reliability of the symmetry point. Obviously, when the line shows a large degree of signal energy variation, the most symmetrical point is also considered more reliable since, despite the variability of gray value data, all difference terms around the true symmetry center amount to a small value. The inverse of the standard deviation of data around the gray value of the hypothesized symmetry center is therefore a good candidate to measure the reliability of the hypothesized symmetry center. This reliability measure can thus be expressed as          R        (     c   ,   j     )       =     1   -     1           σ   2          (     c   ,   j     )       +   1                 with                   σ   2          (     c   ,   j     )       =                  1   m          (         ∑     k   =   1     m            [       g        (       c   -   k     ,   j     )       -     g        (     c   ,   j     )         ]     2       +                                        ∑     k   =   1     m            [       g        (       c   +   k     ,   j     )       -     g        (     c   ,   j     )         ]     2       )       .                                
     The additive term 1 in the denominator avoids zero division for perfectly fiat gray value profiles. obviously, 0&lt;R(c,j)≦1. The associated sequence of reliability measures for each of the mos t symmetrical points in a column j is defined as 
     
       
         R h (j)=R(i=C h (j),j) j=0, . . . , N−1 
       
     
     The computations of {overscore (S)}(c,j) and R(c,j) are further refined in two ways: 
     points (c,j) belonging to either the collimation area or direct exposure area are excluded from consideration since potential symmetry points must lie in diagnostic areas. The availability of the collimation areas and direct exposure areas as bitmap images is segmented on the original x-ray as disclosed in European patent application 97201683.6. These points are disregarded by putting {overscore (S)}(c,j) and R(c,j) to their corresponding extreme values 
     
       
         {overscore (S)}(c,j)=LARGE_FLOAT and R(c,j)=0.0 
       
     
     terms in the sums of {overscore (S)}(c,j) and R(c,j) are disregarded whenever they exceed an upper limit on the difference value. This threshold is empirically determined as a fraction of the full scale range L: 
     T diff =L/4 
     This has the effect of ignoring large valued terms due to gray values belonging to direct exposure areas. 
     At this point it is clear that the concept of minimizing sums of difference values along trajectories perpendicular to the symmetry line searched for and running through the hypothesized symmetry center is applicable for other kinds of symmetry as well. 
     Rotational symmetry, in which a two-dimensional object shows symmetry around its center point (mostly coinciding with the object&#39;s center), can be detected by finding the point at which the sum of difference values between points lying symmetrically at either side of a hypothesized symmetry center, is minimal. An alternate form of rotational symmetry is presented by the N-degree rotational symmetry in which an object is identical to its rotation by N degrees around its symmetry center. Here, differences between rotated points are accumulated and minimized. 
     Rotational symmetry frequently occurs in radiation images of man-made objects such as rings, wheels, gears etc. 
     Fitting a Polynomial to the Symmetry Line 
     The symmetry line computation is completed by idealizing the empirically found sequence C h (j) or C v (i) to a polynomial equation. For most examination types (such as thorax, pelvis, skull, abdomen), the symmetry line is a straight line. The regression model is the well-known linear combination of 1 and x (namely a+bx). More general, the model solved for can be a linear combination of any K specified functions of x. Here x stands for the running abscissa coordinate along rows i for the vertical symmetry axis fitting or columns j for the horizontal symmetry axis fitting. Prior art such as described in “W. H. Press, Numerical Recipes in C, The art of scientific computing, Cambridge University Press, 1st ed. 1988, Chapter 14, Modeling of Data” is used to solve a general linear least squares fitting or regression. For example, the functions could be 1,x,x 2 , . . . ,{circumflex over (x)} K−1 , in which case the general linear combination 
     
       
         y(x)=a 1 +a 2 x+a 3 x 2 + . . . a K x K−1   
       
     
     is a polynomial of degree K−1. The general form of this kind of model is          y                   (   x   )       =       ∑     k   =   1     K            a   k            X   k          (   x   )                                  
     where X 1 (x), . . . , X K (x) are arbitrary fixed functions of x, called the basis functions. The X k (x) can be nonlinear functions of x. The term linear only refers to the model&#39;s linear dependence on the parameters a k . As best parameters a k , those that minimize the chi-square figure of merit are computed            χ   2     =       ∑     j   =   1     N            [           C   h          (   j   )       -       ∑     k   =   1     K            a   k            X   k          (   j   )               1   -       R   h          (   j   )           ]     2         ,                          
     where the notation is given for the horizontal symmetry line, with C h (j) the location in each column j of the most symmetrical point and 1−R h (j) its associated measurement error. The influence on the fitting result of the incorporation of measurement errors is such that unreliable symmetry centers will be taken less into account. Hence, the fit result will have the effect to only ‘follow’ reliable data points, and ‘ignore’ the others. This property is particular advantageous for examination types having portions of relatively flat gray value regions where it is difficult to assign a true symmetry center. The symmetry center points of those regions are disregarded in the computation and do not influence the correct symmetry line. 
     Standardization of Image Viewing Conditions 
     The most symmetrical orientation, either horizontally or vertically, can be derived on the basis of the sequence of symmetry measures 
     S h (j) j=0, . . . , N−1 
     S v (i) i=0, . . . , M−1 
     by computing their associated symmetry index          SI   h     =       1   N            ∑     j   =   0       N   -   1              S   h          (   j   )                                  
     resp.          SI   v     =       1   M            ∑     i   =   0       M   -   1              S   v          (   i   )                                  
     and taking the smaller of the two as the symmetry index and associated main orientation of the image, i.e. 
      SI=min(SI h ,SI v ) 
     The orientation angle α is derived from the symmetry line and the image is rotated by an angle −α to make the orientation coincide with the perfectly horizontal or perfectly vertical orientation depending on the applicable case. Prior art is used to rotate a digital image artifact-free over a given angle, such as described in “G. Wolberg, Digital Image Warping, IEEE Computer Society Press Monograph, 1990, Chapter 7.3. Rotation”. The rotated image is subsequently sent to a soft-copy or hard-copy device for display, review or diagnosis. 
     Stepwise Summary of Operations 
     A stepwise description of operations for performing the image standardization algorithm is given in the sequel. 
     Step 1: for each sub-image of a possibly partitioned image 
     Step 2: for each column j 
     Step 3: for each pixel cat position i in column j, T≦c≦M−T−1, excluding pixels in collimation areas or direct exposure areas 
     Step 4: Compute 
     the maximum range for the subsequent sums 
     
       
         m=min(i,M−i−1) 
       
     
     the symmetrical difference sums around hypothesized symmetry center c          S        (       i   =   c     ,   j     )       =       ∑     k   =   1     m            [       g        (       c   +   k     ,   j     )       -     g        (       c   -   k     ,   j     )         ]     2                              
     the normalized symmetry measure            S   _          (     c   ,   j     )       =       S        (     c   ,   j     )             g   u          (     c   ,   j     )       ·     m   v                                
     the associated reliability measure          R        (     c   ,   j     )       =     1   -     1           σ   2          (     c   ,   j     )       +   1                                  
     Step 2′: repeat for the next pixel in the column 
     Step 4: From the sequences of symmetry measure and associated reliability measure 
     find the most symmetric point in the column 
     
       
         C h (j)=min({overscore (S)}(c,j), c=T, . . . , M−T−1) 
       
     
     find the associated reliability 
     
       
         R h (j)=R(i=C h (j),j) 
       
     
     Step 3′: repeat for the next column j+1 
     Step 5: Form the sequences of symmetry centers and associated reliability measure for all columns 
     C h (j) j=0, . . . , N−1 
     R h (j) j=0, . . . , N−1 
     Step 6: Compute the linear least squares fitting on the data tuples (j,C j (j))and associated measurement error 1−R h (j) to find the horizontal symmetry line SL h  of the examination body part. 
     Step 7: Compute the symmetry index SI h  corresponding to the horizontal symmetry line 
     Step 8: Repeat steps 3 through 7 to obtain the vertical symmetry line SL v  and associated symmetry index SI v . 
     Step 9: Retain the smaller of the symmetry indexes (SI h ,SI v ) for the horizontal and vertical hypothesized symmetry lines (SL h ,SL v ) as the final image&#39;s main symmetry line SL and associated symmetry index SI. 
     Step 10: Derive the orientation α of the retained symmetry line SL. 
     Step 11: Rotate the image over an angle −α to its standard viewing position. 
     Step 1′: Repeat steps 2 through 11 for each sub-image of a possibly multiply exposed image