Software Patent: Method for correcting artefacts in a digital image
Table of contents:
Bibliography
|
Title |
|
|
Title - german |
|
|
Title - french |
|
|
Granted. Opposition window closed |
|
|
EP19980203786 19981110 |
|
|
|
|
|
JP2000172839-A, EP1001370-B1 |
|
|
ep1001370 |
|
|
EP19980203786 19981110 |
Abstract
||<#ffd0d0>[[patent( A method of correcting artefacts in a processed digital image represented by a multi-scale gradient representation being obtained by non-linear processing of a multi-scale gradient representation of an original image. A modified gradient representation is generated by modifying gradient images of the multi-scale gradient representation of the processed image so that the directions of the gradient images are pixel-wise equal to those of the gradient images of the multi-scale gradient representation of the original image. A reconstruction algorithm is then applied to the modified gradient representation. )]] ||
Description
[[patent( Description of <strong>EP1001370</strong><br/> <br/> <br/> Field of the invention<br/>
<br/> [0001] This invention relates to digital image processing.<br/>
<br/> <br/> Background of the invention<br/>
<br/> [0002] In imaging systems where the final output image has to be evaluated by a human observer, a problem arises if the original image, as obtained from an image sensing device, contains detail information at various degrees of coarseness, within a wide amplitude range. This situation is typical when the sensor has a good signal to noise ratio over a large dynamic range, which is the case with CR and CT. In common working environments such as a hospital there is no time for selecting the optimal window/level setting (if any) for each individual image, so there is a need to display a single image on film or monitor screen, which reveals all the relevant diagnostic details with an acceptable contrast over the whole dynamic range. This problem has been recognized in the field of digital radiography, see:<br/>
<br/> [0003] Many methods for contrast enhancement have been developed, such as the commonly known techniques of unsharp masking (Cocklin M.L., Gourlay A.R., Jackson P.H., Kaye G., Kerr I.H., and Lams P., Digital processing of chest radiographs, Imange and Vision Computing, Vol. 1, No. 2, 67-78 (1983)) and adaptive histogram modification (Pizer S.M., Amburn E.P., Austin J.D., Cromartie R., Geselowitz A., Greer T., Ter Haar Romery B., Zimmerman J.B., and Zuiderveld K., Adaptive histogram equalisation and its variations, Computer Vision, Graphics, and Image Processing, Vol. 39, 355-368 (1987)). These traditional methods suffer from two problems: <br/>
<br/> [0004] To overcome the first problem, multiscale methods for contrast enhancement have been developed in the past few years, see Vuylsteke P., Schoeters E., Multiscale Image Contrast Amplification (further referred to as MUSICA algorithm), Proceedings of SPIE, Vol. 2167, 5510560 (1994) and Lu J., Healy D.M. Jr., Contrast enhancement via multiscale gradient transformation, In Proceedings of SPIE: Wavelet applications, Orlando, FL (1994);<br/>
<br/> [0005] While effectively curing this first problem, these methods also seem to perform well with regard to the second problem. The MUSICA algorithm is being used routinely in several hospitals around the world and there are no complaints about artefacts in CR.<br/>
<br/> <br/> Objects of the invention<br/>
<br/> [0006] It is an object of the present invention to provide a method which corrects for artefacts originating from nonlinear treatment of components (coefficients) in a multi-scale representation. <br/> <br/> [0007] Further objects of the present invention will become apparent from the description hereafter. <br/> <br/> <br/> Description of the invention<br/>
<br/> <br/> Statement of the invention<br/>
<br/> [0008] The objects of the present invention are achieved by a method as set forth in claim 1. <br/> <br/> [0009] The method of this invention is based on a multiscale gradient representation such as described in: Mallat S., Zhong S.,<br/>
<br/> [0010] The method has two basic steps. The first step is the nonlinear processing step, e.g. a contrast enhancement step (e.g. using an algorithm similar to that of Lu and Healy (see above)).<br/>
<br/> <br/> Detailed description of the present invention<br/>
<br/> <br/> Multiscale gradient representation<br/>
<br/> [0011] A digital image f is a map from a square grid of pixels P to the real numbers <br/> EMI3.1<br/> . The multiscale gradient representation is a decomposition of the image into a series (say K, K=1,2,3, ...) of vector-valued images on P (i.e. maps from P to <br/> EMI3.2<br/> ) W0f, ..., WK-1f and a residual image SKf , each having the size of the original image (no subsampling takes place). Let S0f = f and for j = 1, ... , K, Sjf is obtained from Sj-1f by applying a smoothing filter, then Wjf is defined to be the (discrete) gradient of Sjf . We will refer to W0f , ..., WK-1f , SKf as the multiscale gradient representation (MGR) of digital image f. <br/> <br/> [0012] The smoothing and gradient filtering are obtained by two-dimensional convolution.<br/>
<br/> [0013] The filters at multiple scales are generated from a basic filter at the smallest scale. The corresponding filters at subsequent larger scales j=1, ...K, are obtained by inserting 2<j>-1 zero-valued coefficients between any two coefficients of the basic filter. <br/> <br/> [0014] In a preferred embodiment, a basic smoothing filter h0 is used with coefficients (1/16, 4/16, 6/16, 4/16, 1/16).<br/>
Sjf = (hj-1 (x)hj-1) * Sj-1f , for 0<j≤K<br/> Wjf = ((gj (x) delta 0)*Sjf , ( delta 0 <br/>
<br/> [0015] Clearly Sjf is a smoothed version of Sj-1f , while Wjf is the (discrete) gradient of Sjf at scale j.<br/>
<br/> [0016] The decomposition is schematized in Figure 1 ( in this Figure, superscripts (1) and (2) refer to the horizontal and vertical component of the gradients respectively).<br/>
<br/> [0017] Filters pj and qj at larger scales j are generated by inserting 2<j>-1 zero coefficients between successive taps, as described above. <br/> <br/> [0018] The IMG reconstruction is carried out as explained in Zong, S, "Edge representation from wavelet transform maxima", Ph.D. Thesis, New York University, 1990 as follows: <br/>
<br/> [0019] The image is recovered from the low-resolution image SKf and the gradients W0f,...,WK-1f by repeated use of this equation for j=K-1,...,0. (see Figure 2). <br/> <br/> [0020] The filters p0 and q0 should be sufficiently smooth to avoid discontinuities when reconstructing the images from a set of nonlinearly modified transform coefficients (as will always be the case when doing contrast enhancement). <br/> <br/> [0021] In a preferred embodiment p0 = (1/4,0,1/2,0/1/4) and q0=(-1/32,-5/32,-13/32,-25/32,25/32,13/32,5/32,1/32).<br/>
<br/> [0022] In a first embodiment, the number of pixels remains constant at each scale. Such a representation is commonly called a multiscale stack. <br/> <br/> [0023] In a second embodiment, the numbers of pixels is lowered in each step of the decomposition, by subsampling. Such a representation is commonly denoted by a multiscale pyramid. <br/> <br/> [0024] The MGR decomposition is then defined as : <br/> EMI6.1<br/>
<br/> [0025] Here &darr& denotes downsampling by a factor of two, keeping only the pixels with even coordinates (2x,2y), where (x,y) refers to the pixel coordinates in the current grid. <br/> <br/> [0026] The pixels of the resulting image are renamed from (2x,2y) to (x,y). The number of pixels is halved in both dimensions, in each stage of the pyramidal decomposition. The reconstruction from the MGR pyramid is done using the equation <br/>
<br/> [0027] The reconstruction procedure (5) is equivalent to subsampling the result of the reconstruction (3), taking into account that p'_0 (x) p'_0) &darr& f = (p0 (x) p0)f, for any arbitrary image f, if p0 has zero coefficients between subsequent taps. The latter condition is fulfilled if p0 = (1/4,0,1,2,0,1/4) <br/> <br/> [0028] Equation (5) tells us how to obtain a subsampled version &darr& S'_jf of S'_jf. The full-size reconstructed image S'_jf is next calculated from &darr& S'_jf and W'_jf (W'_jfbeing the gradient of S'_jf) in the following way.<br/>
<br/> [0029] The pixels of the full-size reconstructed image at positions (2x+1,2y) lying in between the two sample pixels (2x,2y) and (2x+2,2y) are computed by: <br/> EMI7.1<br/>
<br/> [0030] The pixels with coordinates (2x,2y+1) lying in between the two sample pixels (2x,2y) and (2x,2y+2) are computed by : <br/> EMI7.2<br/>
<br/> [0031] The pixels with coordinates (2x + 1,2y + 1) in between the four sample pixels(2x,2y),(2x+2,2y),(2x,2y+2) and (2x+2,2y+2) are computed by: <br/> EMI7.3<br/>
<br/> [0032] We will refer to the latter procedure (described by equations 6-8) as gradient based interpolation (GBI). <br/> <br/> [0033] The pyramidal reconstruction scheme is depicted in Figure 4. <br/> <br/> <br/> Contrast enhancement step<br/>
<br/> [0034] First the multiscale gradient representation W0f,...,WK-1f,SKf (or W'_0f,...,W'_K-1f,S'_Kf in case of the pyramid) of the image f is calculated. <br/> <br/> [0035] To enhance the contrast in the image we want to modify the transform coefficients in such a way that details of small amplitude are enlarged at the expense of the ones with larger amplitude, and this uniformly at all scales. Therefore the transform coefficients must be modified in a nonlinear, monotonic way. <br/> <br/> [0036] In preferred embodiment, the following contrast enhancement function is used: <br/>
y(x) = <FENCE>m DIVIDED x</FENCE><p> if c≤x≤m <br/>
<br/> [0037] The function y is continuous and monotonically decreasing, has a constant value <br/>
<br/> [0038] The modified multiscale representation V0,..., VK-1, SKf (alternatively V'_0,...,V'_K-1,S'_Kf in case of a pyramidal decomposition) is defined as <br/> EMI9.1<br/>
<br/> [0039] This has the following effect: gradient vectors with small magnitude are multiplied by a larger factor than gradient vectors with large magnitude, so details of small amplitude are favoured at the expense of more striking ones. <br/> <br/> [0040] The enhanced output image f is then calculated by carrying out the IMG on V0, ..., VK-1, SKf (or V'_0,...,V'_K-1,V'_Kf).<br/>
<br/> [0041] Note that in this modification, only the gradient magnitudes are affected, not their directions. This is done to prevent artefacts upon reconstruction. However, this goal is only partly achieved, because the method treats the transform coefficients as if they were independent of each other.<br/>
<br/> [0042] The distortion of edge gradient orientation due to independent modification of MGR pixels at subsequent scales, is an unwanted effect, and may give rise to visible artefacts. <br/> <br/> <br/> Artefact correction step<br/>
<br/> [0043] The above problem is solved if the reconstruction is applied to a set of gradient images X0,X1,...,XK-1 instead of V0,...,VK-1 (or V'_0,...,V'_K-1) which comply to the following constraints <br/>
EMI10.1<br/>
<br/> [0044] In D3, the term 'more equalized' is to be understood in the sense of equation (10) : Xj is called more equalized than Wjf when, if<br/>
<br/> [0045] In other words, the artefact correction method aims at redirecting (D2) the gradient vectors of the enhanced (D3) image to the gradient directions of the original image in a way consistent with the multiscale gradient representation structure (D1). <br/> <br/> [0046] We describe below two alternative embodiments of methods that construct such a set X0,X1,...,XK-1. <br/> <br/> <br/> - Alternating convex projections,<br/>
<br/> [0047] It is well known that projecting alternatively onto two convex nondisjoint subsets of a linear space is a convergent scheme that yields a point in the intersection of these two sets (see Youla D., Webb H., Image restoration by the method of convex projections, IEEE Trans. on Medical Imaging 1, 81-101 (1982)). <br/> <br/> [0048] This knowledge is used in a first embodiment to construct gradient images satisfying the above demands D1,D2 and D3. <br/> <br/> [0049] The method starts from the modified representation, V0,...,VK-1,SKf (or from V'_0,...,V'_.-1,S'_Kf) defined in (10) and iterates an alternating sequence of projection operations onto the two vector spaces as follows : <br/>
EMI11.1<br/>
EMI11.2<br/>
EMI11.3<br/>
<br/> [0050] The result of the first projection is fed into the second projection. The resulting vector-images are again fed into the first projection stage, and so on. <br/> <br/> [0051] The first projection onto the space W of multiscale gradients is done by performing on a set of vector-images the multiscale gradient reconstruction procedure (equations (3) or (4), (5), (6), (7) and subsequently applying to the resulting image the multiscale gradient decomposition procedure (equation (1) or (4)). This projection operation is denoted by PW. <br/> <br/> [0052] The second projection operation onto the vector space V PARALLEL of images Xj of which the vector-valued pixel values are parallel to the corresponding gradient pixels of the multiscale gradient representation Wjf, is done by taking the orthogonal projection onto Wjf(i) (Wjf(i) for every i and j. The orthogonal projection of a vector v onto the vector w is equal to <br/>
<br/> [0053] The output image is obtained by applying the IMG reconstruction process to the result of iterating the above alternation of projection operations <br/> EMI12.1<br/>
<br/> [0054] We need approximately 10 iterations for convergence, i.e. we need to do the operations PW and PV PARALLEL 10 times. This algorithm results in a multiscale gradient representation that satisfies all three above constraints D1.,D2. and D3. <br/> <br/> <br/> - Multiscale gradient projection<br/>
<br/> [0055] A second embodiment of the artefact correction process is described next, which is based on the multiscale stack embodiment of the MGR; subsequently a third embodiment of artefact correction based on the pyramidal representation will be disclosed. <br/> <br/> [0056] The second embodiment is based on the fact that for every j, 0 ≤ j < K -1 there exists a linear operator Fj on the space of vectorimages such that for every image f <br/> EMI12.2<br/>
<br/> <br/> Indeed,<br/>
<br/>
EMI13.1<br/>
<br/> [0058] One of the advantages of the z-transform is that it converts a convolution into a product, i.e. for two filters a and b the z-transform of a*b is A(z)B(z).<br/>
<br/> [0059] This yields <br/>
<br/> [0060] The artefact correction method according to a second embodiment starts by decomposing the enhanced image f into its multiscale gradient representation, yielding U0,..., UK-1,SK. <br/> <br/> [0061] In accordance with the criterion D2, these vectors have to be redirected in the direction of the multiscale gradients of the original image. <br/> <br/> [0062] For j=O,...,K-1, let Pj denote the pixelwise orthogonal projection operator onto Wjf. 1-Pj is then the pixelwise orthogonal rejection w.r.t. Wjf. The orthogonal rejection v ORTHOGONAL of a vector v w.r.t. the vector w is given by the formula <br/>
<br/> [0063] As a result of projecting U0 onto W0f according to P0 some information is lost and this can lead to a loss of contrast upon reconstruction.<br/>
EMI14.1<br/>
EMI14.2<br/>
EMI14.3<br/>
<br/> [0064] The factor 1/2 comes in because the filter b0 used in the operation F0 has gain (=sum of all filter taps) equal to 2, so by adding just <br/> EMI15.1<br/>
<br/> [0065] By repeating this procedure the vector image at smallest scale X0 is defined by: <br/> EMI15.2<br/>
<br/> [0066] This will be the full contribution to scale 0.<br/>
EMI15.3<br/>
EMI15.4<br/>
<br/> [0067] The final output image is then reconstructed by applying the IMG to <br/> EMI15.5<br/>
<br/> [0068] According to our findings, <br/> EMI15.6<br/>
<br/> [0069] This algorithm is faster than the alternating convex projection algorithm and preserves better the contrast gained in the contrast enhancement step. <br/> <br/> [0070] We now indicate what changes should be made when using the pyramid decomposition scheme, according to a third embodiment of artefact correction. <br/> <br/> [0071] There still exists a linear operator F'_0 such that for every image f <br/>
<br/> [0072] The filter b0 is the same as defined before, but the operator F'_0 also involves subsampling, denoted by &darr& . The discarded components of projection (1-Pj) are added back in a way similar to the second embodiment, but the effect of subsampling has to be encompassed by proper interpolation. <br/> <br/> [0073] In our preferred embodiment bilinear interpolation was chosen, denoted by the operator I. This operator is applied to each of the contributions at larger scales, in order the match the image dimension of the next smaller scale. <br/> <br/> [0074] The corrected gradient images Xj are then obtained by : <br/>
<br/> [0075] Of course, when projecting a gradient field at some larger scale onto a gradient field at a smaller scale, we have to identify the corresponding pixels, i.e. if Y is k scales larger than Z the pixel (x,y) of Y should be identified with the pixel (x2<k>,y2<k>) of Z . )]] |
Claims
The claims shown below can differ from the final claims filed EPO. Please read the original
documents, if such a difference is of importance.
[[patent( Claims of <strong>EP1001370</strong><br/> <br/> <br/> 1. A method of correcting artefacts in a processed digital image f represented by a multi-scale gradient representation <br/> EMI18.1<br/>
EMI18.2<br/>
<br/> 2. A method according to claim 1 wherein <br/>
EMI18.3<br/>
EMI18.4<br/>
EMI18.5<br/>
<br/> 3. A method according to claims 1 or 2, wherein the multiscale gradient representation of the original image is obtained by (i)applying one-dimensional gradient filters independently to the rows and columns of a digital image representation of the original image thereby yielding the horizontal and vertical components of the multiscale gradient image at the smallest scale, and applying a two-dimensional low-pass filter to said digital representation of the original image,thereby yielding an approximation of the original image at a next larger scale, (ii)identifying the above operations as the first step of an iterative loop, and performing additional iterations using the approximation image resulting from a previous iteration instead of the original image, with the filters being adapted to the current scale, thereby yielding gradient and approximation images of the original image at the successively larger scales. <br/> <br/> 4. A method according to claim 3, where said gradient filters have filter coefficients (0.5, -0.5) or a multiple thereof and where the low-pass filtering is obtained by applying a horizontal one-dimensional low-pass filter to the rows of an image, followed by a vertical one-dimensional low-pass filter applied to the columnns of the resulting image of the horizontal filter, both low-pass filters having coefficients selected from the following set:(0.5, 0.5), (0.25, 0.5, 0.25), (0.0625, 0.25, 0.375, 0.25, 0.0625). <br/> <br/> 5. A method according to claim 2, wherein constraints (1) to (3) are fulfilled by <br/>
<br/> 6. A method according to claim 2 wherein constraints (1) to (3) are fulfilled by : <br/>
EMI20.1<br/>
EMI20.2<br/>
EMI21.1<br/>
<br/> 7. A method according to claim 1 wherein said non-linear processing is a contrast enhancing processing. )]] |
Comment
Sources
The data shown above was obtained from:
Espacenet bibliography
v3.espacenet PDF source
l2.espacenet bibliography
EPOline register
Epoline dossier
The data in this page was last updated . Read Legal Disclaimer.