Automatic processing of digital X-ray bilateral filtration method medical images by

Abstract

The method of bilateral filtering of digital X-ray medical images has been improved by automatic selection of the filter kernel parameters, namely the standard deviations of the bilateral filter kernel in the spatial region and in the region of brightness. The parameters of the filter kernel are calculated using the noise level in the image and the average spatial period of the image. Its standard deviation is used as the noise level. The mean spatial period of an image is calculated by its Fourier spectrum. Bilateral filtering provided a reduction of Gaussian and pulsed noises while the image contours remained clear. To improve the visual quality of images, the local contrast increasing and gamma correction were also performed after bilateral filtering. The proposed methods of image processing are implemented as software in the Matlab system. The results of processing experimental X-ray medical images by the proposed method showed a significant increase of their signal-to-noise ratio and visual quality.

Keywords

• 1 Digital X-ray medical images
• noise level
• bilateral filtering
• image processing
• gamma correction

1. Introduction

Processing of experimental digital X-ray medical images is complicated due to the presence of significant noise levels and low contrast of individual areas of images [ 1-3 ]. Such image imperfections lead to decreasing their signal-to-noise ratio (SNR) and visual quality. As a result, it reduces the accuracy of detecting the details of the studied objects. The peculiarity of X-ray medical images is that their imperfections are difficult to eliminate during image formation in sensors because X-ray dose for patients is restricted. For this reason, the urgent task is to reduce the noise level in medical images by digital filtering. However, known filtering methods using, for example, a Gaussian or median or wavelet filters, lead to blurring of image contours and reduce their detail [ 4-5 ]. Wavelet filters provide less contour blur compared to a Gaussian filter, but require a fairly complex procedure for selecting a family of wavelets and thresholds for wavelet coefficients.

Therefore, in this paper we propose to remove noise in images using a bilateral filter, which provides a significant increase of SNR while the image contours remain clear [ 6-8 ]. Due to the development of the bilateral filtering method, the automatic selection of the filter kernel parameters is provided, namely the standard deviations of the kernel in the spatial region and in the brightness region. At the same time, automating the selection of filter kernel parameters not only increases the speed of the method, but also reduces the influence of subjective factors on the result of image filtering, which are peculiar to manual processing. For better visualization of medical images, it is also advisable to increase their local contrast and perform gamma correction of images. This processing of X-ray medical images can notably increase their visual quality, which significantly increases the efficiency of subsequent image analysis, in particular, the accuracy of diagnosis of patients.

2. Improving the visual quality of X-ray medical images

High visual quality of medical images is a prerequisite for an objective assessment of patients. The paper proposes to improve the visual quality of images by the method of bilateral filtering, because such filtering reduces noise while the image details important for their analysis remain clear. 2.1.

Bilateral image filtering algorithm

Bilateral image filtering is designed to reduce the noise level on experimental X-ray medical images and is performed according to the following algorithm (Fig. 1). The first step is reading the initial image f, which contains whole image of the object under study or single part of it. Digital images f are processed as rectangular matrices f = f (i, k), where i = 1, ..., M; k = 1, ..., N; M is the height of the image (in pixels), N is the width of the image (in pixels) [ 4-5 ]. The brightness of the images is normalized in the range from 0 to 1.

The parameters of the bilateral filter kernel depend on the noise level σNE in the image and the spatial period TS of the image, so the next steps of the algorithm correspond to calculation of the image parameters σNE and TS. Medical X-ray images mainly contain Gaussian and impulse noise, the total level of which can be described by the standard deviation σNE. The calculation of the noise level is performed by a high-precision method using low-frequency filtering when selection of the noise component and taking into account the Region Of Interest (ROI), in which the noise component is dominated [ 9-11 ].

To find the spatial period TS of the image f, at first its Fourier spectrum F calculated using twodimensional Discrete Fast Fourier Transform (DFFT) [ 3, 5, 12 ]

=1 =1 ( , ) = ( , ) ∙ − ∙ 2

( − 1)

( − 1) + ,

An important parameter of the radial distribution is its average radial spatial frequency vS, calculated by the formula [ 3, 13 ] where

m is frequency number (index) in height; n is frequency number in width; m = 1, 2, ..., M; n = 1, 2, ..., N; M is the height of the digital image (in pixels); N is the width of the digital image (in pixels); j is an imaginary unit.

A centered Fourier spectrum FC is obtained from the Fourier spectrum F, where the centre of the frequency rectangle corresponds to zero frequencies. Next, the power spectrum PS (or power spectral density) [ 5 ] of the image calculated, which is equal to the square of the modulus FC:

To obtain the power spectrum PS, its radial distribution PR (d) is calculated, where d = 0, 1,…, NR, NR = [N / 2], d are integer values of the distance from the element of the spectrum (m, n) to its center. Each distance d for the radial distribution PR (d) corresponds to the value of the radial spatial frequency

The spatial period TS of the image is calculated based on the frequency vS

= | |2. ( ) = / .

= 1 =1 =1 = ( ) ( ) ( ).

(1) (2) (3) (4) (5)

After that, the parameters σSB and σB of the kernel wB of the bilateral filter are calculated on the basis of predefined image parameters σNE and TS. The kernel of the filter wB = (wB (m, n)) is described by the formula: (6) (7) where m = 1, ..., Mw; n = 1, ..., Nw; m is the row number of the kernel elements; n is the column number; Mw, Nw are the sizes of the filter kernel in height and width, respectively; σSB is standard deviation of the bilateral filter kernel in the spatial region; σB is standard deviation of the of the bilateral filter kernel in the brightness region; mc and nc are coordinates of the kernel center of the filter wB in height and width, respectively; fw (m, n) is the brightness of the image pixel, which corresponds to the kernel element with numbers (m, n); fC is the brightness of the image pixel that corresponds to the kernel center.

The value of the elements sum of the filter kernel wB is normalized to 1.

According to formula (6), the kernel of the bilateral filter wB can be described by the elemental product of two kernels wSB and wBB: ( , ) = exp −( − )2 + ( − )2 2 2 where m = 1,..., Mw; n = 1,..., Nw.

The method of bilateral noise filtering [ 6-8 ] allows keeping the clarity of the contours, because such a method applies spatial weighted averaging of the image brightness. In the method of bilateral filtering, two Gaussian filters are used: one filter performs processing in the spatial region (with standard deviation σSB), and the other – in the brightness region (with standard deviation σB).

The dimensions of the filter kernel wB are calculated taking into account the rule 3σ for the twodimensional Gaussian distribution by the formula:

= [6 ∙ ], = [6 ∙ ]. (9)

In the case of studied X-ray medical images, the average spatial period is in the range from 32 to 64 pixels. Therefore the dependences of σSB (σNE, TS) are investigated by image filtering (previously the image brightness is modulated by the sum of two sinusoids with period TS in width and height and Gaussian noise with σNE level was artificially added) for period values of 32 and 64 pixels. When TS = 32 pixels, the dependence of σSB (σNE) is described by the empirical formula: where с1 = 0.16; с2 = 6.5.

When TS = 64 pixels, the dependence of σSB (σNE) is described by the empirical formula: 32 = 1 + 2 ∙ 64 = 3 + 4 ∙ , where с3 = 0.16; с4 = 10.2.

Therefore, taking into account formulas (10) and (11) at 32 ≤ TS ≤ 64, the dependence of σSB (σNE, TS) is described by the empirical formula:

64 − − 32 (12) = 32 ∙

+ 64 ∙ 32 32 .

(8) (10) (11) (13) (14) Formula (12) can also be written as = с1 + 2 ∙ 64 − 32 + 4 ∙ − 32 32 , where с1 = 0.16; с2 = 6.5; с4=10.2.

According to the rule 3σ for the normal distribution the standard deviation of bilateral filter kernel in the brightness region is calculated by the formula:

= 3 ∙

Further, values of the elements of the kernel wB, are calculated in cycles with counters i and k (Fig. 1); then convolution of these values with the corresponding fragment of the image f gave the value gB (i, k) for the pixel of the filtered image. Thus, as a result of bilateral filtering the image gB is calculated. 2.2.

Software implementation of the method of bilateral filtration

The software "p_Bilateral_Filter" for bilateral image filtering is created on the basis of the algorithm (Fig. 1) by means of Matlab [ 5 ] as a set of m-files. The program provides reading images in various graphic formats (including tiff and jpg formats), as well as in DICOM format (Digital Imaging and Communications in Medicine). In this case, the image fragment f with the most important areas is selected from the complete initial image fA and further processing is performed over this fragment.

The DICOM format [ 14 ] describes the medical industry standard for the creation, storage, transmission and visualization of digital medical images and patient examination documents. An important attribute of the DICOM file is "BitDepth", i.e. the number of bits n that describe the brightness or color of the image pixel. The Discrete Fast Fourier Transform is performed using the built-in "fft2" function in the Matlab system. Processed images are stored in DICOM files or in *.tiff image files (with a color depth of 8 and 16 bits).

2.3. Results of X-ray medical image # 1 processing by bilateral and wavelet filtration

Consider an example of processing a typical X-ray medical image of the lungs (image # 1). The complete initial image fA (Fig. 2) is read from a DICOM file (grayscale with a color depth n = 14). The image f is obtained as a fragment of the image fA, which contains an informative (diagnostically important) central area (Fig. 3). The image f shows the roots of the lungs, however, such an image contains a significant level of noise, which reduces its visual quality.

Using the software "p_Bilateral_Filter" developed on the base of the bilateral filtering method the noise level on the image is reduced in the following sequence. At first, power spectrum PS and average spatial period (TS = 40,627 pixels) as well as the noise level (σNE = 0.0210) were calculated for the initial image f. Next, the standard deviations of the bilateral filter kernel were calculated on the base of the of TS and σNE values, and the initial image was filtered.

As a result, a filtered gB image was obtained (Fig. 4) with a noise level reduced by more than an order of magnitude, which significantly increases the visual quality of image, its signal-to-noise ratio (by more than an order of magnitude) and the efficiency of subsequent analysis. After noise reduction, clearer visualization of small elements of the medical image (in particular, calcifications) is provided.

In order to improve the visual quality, additional processing of the filtered gB image was performed by increasing the local contrast and gamma correction.

For the low-noise gB image, the local contrast increasing is necessary to provide a high-quality gC image for all its areas (Fig. 5). For this purpose, a fast-acting method of increasing the local image contrast using the enveloping of minimum and maximum values of the image brightness was applied within the local windows [ 15 ]. In the case of excessive local image contrast increasing the artifacts (e.g., halos) can appear, so the maximum level of local contrast increasing is limited by a fixed value of Scale_Max (for example, Scale_Max = 3). Increased local contrast is especially effective in the study of bone structure.

The results of the implemented methods of bilateral filtration and increase of local contrast are especially noticeable in the analysis of profiles [ 16 ] of processed images (Fig. 6). Analysis of the image profile gB (Fig. 6c) shows that bilateral filtering significantly reduces the noise level (compared to the profile of the original image in Fig. 6b), but the contours of objects (e.g. edges) did not lose clarity. Analysis of the gC image profile (Fig. 6d) shows that after increasing the local contrast, the range of changes in the brightness z of the image increases significantly.

For comparison, wavelet filtering of the studied image was also performed using Daubechies wavelets [ 5 ]. Automatic wavelet filtering of the image by Wavelet Toolbox (Matlab), namely by tool "Stationary Wavelet Transform Denoising 2-D", provided partial noise removal only. This noise reduction is almost invisible in the image visually, and the image profile is slightly smoothed (Fig. 6e). Wavelet filtering can reduce the noise level in the image more, but it requires a rather complex procedure for selecting the wavelet family, the order of the wavelet in the family and the thresholds for the wavelet coefficients. Therefore, the method of bilateral filtration is more advisable to use for automatically filtering X-ray medical images.

In order to improve visualization of dark or light areas for the filtered image gB, gamma correction with the correction parameter γ was performed, and resulting image gG was calculated. Due to gamma correction, in particular, the structure of the lung roots is visualized with better quality (Fig. 7). d) e) Figure 6: Profiles of z(r) images: a) initial image f (Fig. 3); b) profile z(r) for the image f; c) image gB profile after bilateral filtration (Fig. 4); d) image gC profile after increasing local contrast (Fig. 5); e) image profile after Daubechies wavelet filtering (order 5); z is image brightness, r is profile length (pixels); Qp is the number of profile points b) Figure 7: Images gG after gamma correction, calculated on the basis of the image gB (Fig. 4): a) γ = 0.5; b) γ = 1.5

2.4. Results of processing the X-ray medical image # 2 using the method of bilateral filtration

Here we consider an example of processing a typical X-ray medical image with implants (image # 2). The complete initial image fA (Fig. 8) is read from a DICOM file (grayscale with a color depth of n = 14). The image f is obtained as a fragment of the image fA, which contains diagnostically important area (Fig. 9). The image f shows an implant, however, this image contains a high noise level, which reduces its visual quality.

Applying developed software "p_Bilateral_Filter" and method of bilateral filtration to the image # 2 (Fig. 9) the noise level is reduced (similarly to the image # 1). At first, the power spectrum and the average spatial period (TS = 47,391 pixels), as well as the noise level (σNE = 0.0125) were calculated for the initial image. Next, based on the values of TS and σNE, the standard deviation of the bilateral filter kernel was calculated and the initial image was filtered. As a result, a filtered gB image was obtained (Fig. 10) with a noise level reduced by more than an order of magnitude, which significantly increases the visual quality of image and the efficiency of its subsequent analysis. After reducing the noise level, a clearer visualization of both bone tissue and implant is provided.

To improve the visual quality, additional processing of the filtered gB image was performed by increasing the local contrast and gamma correction.

On the filtered image gB the local contrast was increased [ 14 ], resulting in an image gC with high visual quality for all its parts (Fig. 11). To prevent the appearance of artefacts (e.g. halos), the maximum increase level in local contrast is limited to a fixed value (Scale_Max = 3). Increasing the local contrast of the image was effective in the study of both bone structure and implant.

Analysis of the studied image profiles (Fig. 12) shows that as a result of bilateral filtering (Fig. 12c) the noise level is significantly reduced (compared to the profile of the original image in Fig. 12b), but the contours of the objects remain clear. After increasing the local contrast, the range of changes in the brightness z of the image increases (Fig. 12d). d) Figure 12: Profiles z(r) of images: a) initial image f (Fig. 9); b) profile z(r) of image f; c) profile of image gB after bilateral filtration (Fig. 10); d) profile of image gC after increasing local contrast (Fig. 11); z is image brightness, r is profile length (pixels); Qp is the number of profile points

The better visualization of dark or light areas for the filtered image gB was obtained by gamma correction of image (Fig. 13). b) Figure 13: Image gG after gamma correction, calculated on the basis of image gB (Fig. 10): a) γ = 0.5; b) γ = 1.5

Brightness inversion was also performed for the filtered image (Fig. 14), because details in dark areas are visible clearer.

A series of other X-ray medical images were processed similarly. In all cases, a significant increase in the visual quality of the images was obtained. 2.5.

Conclusions

The method of bilateral filtering has been improved by automatically calculating the parameters of the filter kernel based on the noise level of the image and its average spatial period. The noise level is calculated as its standard deviation by a method based on low-frequency image filtering. The spatial period is calculated based on the Fourier power spectrum of the image. An empirical formula (12) for calculating standard deviation of the bilateral filter kernel has been developed.

The method of bilateral filtering is software implemented in the Matlab system, the initial images are read from files in DICOM format, as well as jpg and tiff. Due to the automatic determination of filtering parameters, the developed software can replace the human operator in the pre-processing of X-ray medical images.

Experimental testing of an improved method of bilateral filtering has shown that it can significantly reduce the noise level (more than an order of magnitude) keeping the clarity of the contours. This significantly improves the visual quality of images and the accuracy of their subsequent analysis. The advantage of the method of bilateral filtration in comparison with automatic wavelet filtration is shown.

For filtered images, an increase their local contrast, gamma correction and inversion were additionally implemented. This image processing allows you to distinguish even inconspicuous details of objects that located in too dark or too light areas.

3. Acknowledgements

We are grateful to Vasyl V. Lagoda (Vinnytsia Pirogov Regional Clinical Hospital, Vinnytsya, Ukraine) and Tetyana M. Sidorchuk (National Pirogov Memorial Medical University, Vinnytsya, Ukraine) for providing the X-ray medical images used in our study, as well as for useful advice on research topics.

4. References

References

[1] J. Jan , Medical Image Processing, Reconstruction and Analysis. Concepts and Methods , New York, CRC Press, 2019 . doi: 10 .1201/b22391.

[2] E.R. Ranschaert . S. Morozov, P.R. Algra , Artificial Intelligence in Medical Imaging. Opportunities, Applications and Risks , Springer Nature Switzerland, 2019 . doi: 10 .1007/978-3- 319 -94878-2.

[3] T. Olson , Applied Fourier Analysis From Signal Processing to Medical Imaging , Springer Science+Business Media, LLC, 2017 . doi: 10 .1007/978-1- 4939 -7393-4.

[4] S. Krigg , Computer Vision Metrics. Survey, Taxonomy, and Analysis , Apress, Berkeley, CA, 2014 . doi: 10 .1007/978-1- 4302 -5930-5.

[5] R. Gonzalez , R. Woods , Digital image processing, 4th edidion , Pearson/ Prentice Hall, NY, 2018 .

[6] K. Sugimoto , S.I. Kamata , Compressive bilateral filtering , IEEE Transactions on Image Processing 24 ( 2015 ) 3357 - 3369 . doi: 10 .1109/TIP. 2015 . 2442916 .

[7] K.N. Chaudhury , S.D. Dabhade , Fast and provably accurate bilateral filtering , IEEE Transactions on Image Processing 25 ( 2016 ) 2519 - 2528 . doi: 10 .1109/TIP. 2016 . 2548363 .

[8] V. Anoop , P. R. Bipin , Medical Image Enhancement by a Bilateral Filter Using Optimization Technique , Journal of Medical Systems 43 / 240 ( 2019 ) 1 - 12 . doi: 10 .1007/s10916-019-1370-x.

[9] S.V. Balovsyak , Kh. S. Odaiska , Automatic Determination of the Gaussian Noise Level on Digital Images by High-Pass Filtering for Regions of Interest , Cybernetics and Systems Analysis 54 ( 4 ) ( 2018 ) 662 - 670 . doi: 10 .1007/s10559-018-0067-3.

[10] O. Berezsky , S. Verbovyy , O. Pitsun , Hybrid Intelligent Information Technology for Biomedical Image Processing. IEEE 13th International Scientific and Technical Conference on Computer Sciences and Information Technologies (CSIT) , Lviv ( 2018 ): 420 - 423 . doi: 10 .1109/STCCSIT. 2018 . 8526711 .

[11] A.B. Lozynskyy , I.M. Romanyshyn , B.P. Rusyn , Intensity estimation of noise-like signal in presence of uncorrelated pulse interferences , Radioelectronics and Communigations Systems 62 ( 5 ) ( 2019 ) 214 - 222 .

[12] P. Nair , A. Popli , K.N. Chaudhury , A Fast Approximation of the Bilateral Filter using the Discrete Fourier Transform , Image Processing On Line (IPOL) , 7 ( 2017 ) 115 - 130 . doi: 10 .5201/ipol. 2017 . 184 .

[13] A. Phinyomark , S. Thongpanja , H. Hu , P. Phukpattranont , Ch. Limsakul, The usefulness of Mean and Median Frequencies in Electromyography Analysis , in G.R. Naik (Eds.), Computation Intelligence in Electromyography Analysis, Chapter 8 , A Perspective Current Applications and Future Challenges, InTech, 2012 , pp. 195 - 220 . doi: 10 .5772/50639.

[14] O.S. Pianykh , Digital Imaging and Communications in Medicine (DICOM). A Practical Introduction and Survival Guide, Springer-Verlag, 2012 .

[15] S.V. Balovsyak , I.M. Fodchuk , Yu.M. Solovay , Ia.V. Lutsyk , Multilevel method of local contrast increase and images heterogeneous background removal, Cybernetics and Computer Engineering 182 ( 2015 ) 15 - 26 . doi: 10 .15407/kvt182. 02 . 015 . (in Russian)

[16] S.V. Balovsyak , O.V. Derevyanchuk , I.M. Fodchuk , Method of calculation of averaged digital image profiles by envelopes as the conic sections , Advances in Intelligent Systems and Computing (AISC) 754 ( 2019 ) 204 - 212 . doi: 10 .1007/978-3- 319 -91008-6_ 21 .