Automatic Noise Reduction in Ultrasonic Computed Tomography Image for Adult Bone Fracture Detection

Abstract

Noise reduction in medical image analysis is still an interesting hot topic, especially in the field of ultrasonic images. Actually, a big concern has been given to automatically reducing noise in human-bone ultrasonic computed tomography (USCT) images. In this chapter, a new hardware prototype, called USCT, is used but images given by this device are noisy and difficult to interpret. Our approach aims to reinforce the peak signal-to-noise ratio (PSNR) in these images to perform an automatic segmentation for bone structures and pathology detection. First, we propose to improve USCT image quality by implementing the discrete wavelet transform algorithm. Second, we focus on a hybrid algorithm combining the k-means with the Otsu method, hence improving the PSNR. Our assessment of the performance shows that the algorithmic approach is comparable with recent methods. It outperforms most of them with its ability to enhance the PSNR to detect edges and pathologies in the USCT images. Our proposed algorithm can be generalized to any medical image to carry out automatic image diagnosis due to noise reduction, and then we have to overcome classical medical image analysis by achieving a short-time process.

1. Introduction

Bone disease exploration is assured by a variety of modalities of medical imaging. Bone mineral density is determined by standard scanners and x-rays. Although this technique delineates the bone structure, it remains an invasive method with qualitative information on bone structures. Ultrasound computed tomography (USCT) is the best to give us more details and a very interesting procedure for bone imaging associated with signal and image processing methods [1, 2]. A variety of medical imaging techniques such as x-rays and standard scanners determine bone mineral density. However, these modalities represent an ionizing method without giving quantitative results. USCT hardware has been propounded to solve this issue [1, 2]. USCT presents an important radiological technique because of its non-ionizing properties. However, it has been the subject of several studies due to the complexity of its visualization and the noisy USCT image quality [3]. In the light of this issue, noise reduction should be the first step to be interested in, to enhance the USCT image resolution and to detect bone diagnosis. Hence, many researchers have studied different techniques for ultrasound image noise reduction, and obtained results achieved the terms of quality improvement (37.14 dB of peak signal-to-noise ratio [PSNR]) [3]. But, these methods cannot be applied to USCT bone images due to its complexity. Actually, USCT bone imaging is a very difficult method that encounters many problems mainly related to high bone echogenicity. Noise reduction is one of the crucial topics in digital image processing and has been conducted in various fields such as ultrasound imaging [4]. Many obstacles still exist, but several results have been successful [5, 6]. Accordingly, we will try here to propose optimizations for image processing by implementing discrete Haar wavelet algorithm and a proposed hybrid algorithm combining k-means with the Otsu method. It is an important process for removing noise, as it produces promising results in terms of image resolution, noise removal and diagnosis detection. In this work, our objective aims to automatically remove noise levels generated during the rebound of ultrasonic waves against bony structures, from USCT images, improving the PSNR.

In this chapter, we divide the work into seven sections: Section 2 presents an overview of the medical history in this area, some physical considerations and the description of existing algorithms. Section 3 presents fundamental mathematical theorems that will be implemented in the next section. Section 4 describes the proposed hardware and software method to segment USCT. Section 5 shows the obtained results, the comparison of our work with previous work and the discussions. Section 6 concludes this chapter.

2. State of the art

In 2019, a Ram-Lak filter is implemented in the Radon domain to reduce the noise levels in the USCT images and facilitate their observation giving a PSNR with a value of 13.07 [3]. In 2018, different filters such as median and a high pass filter have been used to reduce noise in Magnetic Resonance Imaging (MRI) images as a pre-processing step and the first results are promising [7]. In addition, bilateral filter and trilateral filter are used for noise removal in [8] by considering the small structure as noise to be removed and conserving the large structure. Moreover, in [9] adaptive filter of Kuan, Frost or maximum medium median (MMM) have been produced to enhance ultrasound images resolution giving promising results. However, these methods fail to give encouraging results when applied to USCT images [3, 9].

The Fourier transform with its decomposition of a trigonometric series has been widely introduced in the field of signal processing as well as in medical image processing, but it is insufficient for giving a piece of complete information in both time and frequency areas simultaneously. Indeed, using the classical Born approximation and the spatial Fourier transform, the result is a poor contrast-to-noise ratio (CNR) image. Some previous work to improve the CNR has been introduced while handling signal and image processing [2, 5]. Although image analysis is an important process for noise removal, it cannot produce the intended result in some digital image processing [4]. Actually, a flow of image processing should be applied to improve the USCT image quality and to achieve automatic diagnosis detection. The active contour method was introduced in the process of USCT image segmentation. Its use avoided the issue of noise in USCT images [10]. In [11], authors appended this algorithm on USCT paired-bone images. However, the results were not good enough and the detection of the distances between the two bone forms (tibia and fibula) was not possible.

The method, called the “Wavelet-based Coded Excitation” (WCE) method [5] is based on the wavelet decomposition of the signal and on a suitable transmitted incident wave correlated with the experimental set-up. Indeed, the contrast with reports of noise was improved, but the detection of edges and areas of child-matched bone by ultrasound with USCT was impossible [5, 11]. However, this method has remained very interesting. During the last two decades, a lot of research has involved the use of the transformed wavelet for removing noise through its energy compaction and its multiresolution parameter proprieties [12, 13]. Thanks to these parameters, we can obtain different versions (dilated or compressed, and translated or shifted) from the same mother wavelet.

The k-means algorithm is known for its simplicity in clustering a database into clusters. Nevertheless, this number of clusters has to be identified. It has been used for resolving the issue of data clustering. k-means followed by fuzzy c-means were introduced by Alan Jose et al. in [14, 15] to detect a brain tumor, and the result of the segmented image remained for the feature extraction and the image resolution was improved. Image segmentation using the k-means algorithm is a process of separating images into different regions with high resolution. The purpose of such segmentation is the detection of regions of interest simultaneously with noise removal in an image. Moreover, image segmentation has been an increasingly expanded issue especially in the area of medical imaging and more specifically in USCT considering the inhomogeneity of pixels and complex anatomical topology [6].

Otsu method has been used for a long time to medical image analysis for image resolution enhancement [16] obtaining the right diagnostic. In [7], a proposed Otsu method has been used for MRI images for tumor brain detection.

3. Mathematical theoretical

3.1 Wavelet transform

The wavelet transform is a mathematical function that allows image compression and signal processing. It resolves the problem of the Fourier analyses. It is defined by the following equation:

$$\Psi \left(a,b\right)=\frac{1}{\sqrt{\mathrm{a}}}\Psi \left(t-\frac{\mathrm{b}}{\mathrm{a}}\right)$$E1

where a is the scale parameter or the expansion factor, and b is the translation or proposition parameter. The bigger a is, the more the wavelet is dilated.

3.1.1 Principle of wavelet transform

The wavelet transform is a multiresolution analysis tool able to transfer accurate temporal and spatial information. In the literature, various noise reduction techniques concerning wavelet approaches have been put forward [14]. From an original image, low frequencies are analyzed by the application of a low-pass filter. Then, high frequencies are analyzed by the application of a high-pass filter. It allows dividing the information of the image in an approximation and detail as depicted in Figure 1. We will use the Haar wavelet to keep the edge detection of an image as follows:

• On the lines, estimate the medium of each pixel of paired data.

• On the lines, calculate the variance between each datum and its respective average.

• Point the averages in the first half of the matrix data which corresponds to the L approximation image.

• Point the average in one-half of the data corresponding to the H detailed image matrix.

• Rehearse the process on the first half of the L data, but on columns so we get two images corresponding to the approximation of LL (means of dimension matrix) and LH (matrix of differences in dimension) image details.

• Renew the process on the LL image approximation to have a decomposition level greater than 1.

• A low-pass filter application remains to get an L image which is compressed and the application of a high-pass filter leads to obtain an H image which introduces image details. This process is done by Eqs. (2) and (3) as follows:

Figure 1.

Principle of wavelet application.

$$\mathit{YH}\left[k\right]=\sum _{n}x\left[n\right]G\left[2k-n\right]$$E2

$$\mathit{Yl}\left[k\right]=\sum _{n}x\left[n\right]H\left[2k-n\right]$$E3

3.2 k-means algorithm

It is based on grouping similar data points into clusters. There is no prediction involved. Its algorithm is illustrated by these steps [18].

• Fix the k number of cluster values.

• Identify the k cluster centers.

• Determine the cluster center.

• Determine the pixel distance for each cluster center.

• If the distance is close to the center value, budge to that cluster.

• Otherwise, move to the next cluster.

• Re-identify the center.

3.3 k-means algorithm combined with Otsu algorithm + morphologic algorithm

The proposed hybrid algorithm uses the combination of k-means with the Otsu method and the morphologic algorithm. The k-means algorithm implementation is important. It works well in a large number of cases and it is a powerful tool to have in the closet point. Unfortunately, in medical image processing, it is not sufficient for region detection and edge detection. k-Means combined with the morphologic and Otsu algorithms give us sufficient results.

3.3.1 Morphologic algorithm

Morphological filters are a valuable aid in the segmentation and noise removal process. Morphological filtering is based on mathematical morphological operations, well applied to binary images as well as monochrome (grayscaled) images. To be limited to binary image morphological filtering, morphological operations aim to perfect the improper structure of an image [17]. The morphological erosion of X by B is defined by the principle of duality where X is the set of points described in the space and B is the structuring element. Its equation is written under the following form:

$$\mathrm{\epsilon }\mathrm{B}\overline{\left(\mathrm{x}\right)}=\overline{\mathrm{\delta B}\left(\mathrm{x}\right)}$$E4

$$\mathrm{\epsilon }\mathrm{B}\left(\mathrm{x}\right)=\overline{\mathrm{\delta }\mathrm{B}\overline{\left(\mathrm{x}\right)}}=\overline{\overline{\mathrm{x}}\oplus \overline{\mathrm{B}}}=\mathrm{x}⊝\overline{\mathrm{B}}$$E5

3.3.2 Otsu algorithm

The Otsu method of the threshold is the most powerful and global threshold method. It performs image binarization based on the histogram image shape. It assumes that the image for binarization contains the only foreground and background pixels [19, 20]. Using the simple formula in the Otsu algorithm, we get:

$${\mathrm{\sigma }}^2=\mathrm{\Psi A}{\left(\mathrm{u}\mathrm{A}-\mathrm{u}\right)}^2+\mathrm{\Psi B}{\left(\mathrm{u}\mathrm{B}-\mathrm{u}\right)}^2$$E6

where σ² is the variance between both clusters, Ψ_A is the probability of class A, Ψ_B is the probability of class A, u_A is the average gray of class A, u _B is the average gray of class B, and u is the threshold value which divides the image into two classes A and B. The best threshold u maximizes the variance between both clusters. It computes the optimal threshold by minimizing the intra-class variance that separates the foreground pixels from background pixels [21]. The main purpose of the Otsu method is to find the threshold values where the sum of the values of the foreground and background pixels has to be minimal [22].

3.4 Peak-signal -to noise ratio (PSNR)

The PSNR is defined by the ratio of the power of the signal and the corrupting noise presented in the image.

4. Materials and methods

4.1 Materials

The used hardware device is illustrated in Ref. [5] in section material and method.

4.2 Methods

Our approach is based on a hardware and software co-design to enhance the signal-to-noise ratio (SNR), to automatically segment the USCT images and finally to detect fractures in bones. It is described in Figure 2.

Figure 2.

Synoptic flow of our approach.

4.3 Hardware method

The acquisition of the signal is done in retro diffusion, transmission and diffraction. The frequencies of measurements are performed at 1 MHz. Each transducer plays the role of a transmitter and a receiver. A transmitted signal is sent and diffracted on the object to be imaged, reaching the other seven receivers. There are three configurations per zone. The wave is transmitted, or back-propagated or diffracted as illustrated in Ref. [2].

4.4 Software method

Using an adult bone defect as depicted in Figure 3 (a). We have a USCT noisy bone image in Figure 3 (b). USCT image processing recognition is kept by the application of some algorithms such as the Haar wavelet transformation and the improved k-means, as depicted by the synoptic flow in Figure 4.

Figure 3.

Ultrasonic bone tomography with USCT: (a) adult cortical bone defect (copyright/rights reserved, CNRS-LMA-Marseille), (b) USCT image defect.

Figure 4.

Hybrıd algorithm combining k-means with Otsu method.

4.4.1 Implemented Haar wavelet algorithm

Algorithm 1: Direct Haar wavelets transform

 1: Download the source which contains the converted binary image

 2: Download f-L and f-H filter file.

 3: Treatment following lines.

 4: Treatment following columns.

 5: //Sgl is the converted binary image.

 6: while.

 7: {

 8: Bf [0] and Bf [1] take the sgl values.

 9: CVh = 0 compute the mean $\mathrm{CVh}=\frac{\mathit{Bf}\left(0\right)+\mathit{Bf}\left(1\right)}{2.}$

10: CVl = 0 compute the difference $\mathrm{CVl}=(\mathit{Bf}\left(0\right)-\mathit{CVh}$)

11: }

12: Display the values of CVh, CVl

13: Display LL, HL, LH and HH image values

4.4.1.1 Proposed k-means algorithm

Our proposed k-means algorithm is described in Figure 4.

5. Results

5.1 Fracture detection

Figure 5 shows an example of how to identify an artificial defect on a human femur. The defect is through an incision at the end of the sample whose width is reduced in depth. Three heights are analyzed: one (H1) in the very fine portion of the crack, a second (H2) in the intermediate zone, and the third (H3) in the very open upper zone. We can see in Figures 6 and 7 that the implemented Haar wavelet detects the external edge. Hence, it is very clear, but the fracture is not detected and the problem still exists. It is because of the noise, the inhomogeneity of pixel intensities and the difficulty of separating anatomic regions. We do not have any idea about the intensity of the three regions. The problem is how to detect the fracture and how to obtain the different regions of interest. For this fact, a hybrid algorithm is proposed in the second step.

Figure 5.

Adult bone defect.

Figure 6.

External edge detection by wavelet transform with FMC interface, (a1): LH1, (b1): HL, (c1): HH1.

Figure 7.

External edge detection by wavelet transformation: (a) LL: PSNR = 16.18, (b) LH: PSNR = 24.65, (c) HL: PSNR = 25.15, (d) HH: PSNR = 22.98.

An example of the application of the k-means algorithm in the identification of a bone defect is discussed as an interesting breakthrough work of USCT bone imaging. Despite the ability of automatic fracture detection with k-means, as shown in Figure 8(a), (b) and (c), the results are not sufficient because we need to have an excellent quality of images and to extract the background that noises the ultrasonic image. With k-means, we detect regions and fractures. We need to extract the background and detect the edges. For this fact, we apply the morphologic algorithm. As shown in Figure 9(d), (e) and (f), with the morphologic algorithm, edges are well detected and we can calculate the distances with more precision. We can clearly distinguish the different zones, especially in the H1 part where the rack is not very open. The PSNRs are (H1) 13.08, (H2) 13.11 and (H3) 13.14, respectively. As seen in Figure 10 (g), (h) and (i), due to the application of the Otsu method, we obtain the image histogram value. After that, we extract the background. Finally, we obtain every region independent of the other. We are interested in the internal region which is the medullary cavity with a diameter of 0.8 cm. These results are similar to the reality measures where the cortical bone width D2 is about 3 mm, as given in Figure 8 (b) and the cancellous bone represents a D3 width of 1 mm as depicted in Figure 9(d).

Figure 8.

Denoised USCT image with k-means algorithm: (a) adult bone USCT image, (b) segmented USCT adult bone image with the k-means algorithm, k = 3, (c) segmented USCT adult bone image with the k-means algorithm, k = 2.

Figure 9.

Denoised USCT with k-means algorithm combined with the morphologic algorithm: (d) H1 fracture defect detection in USCT image, (e) H2 fracture defect detection in USCT image, (f) H3 fracture defect detection in USCT image.

Figure 10.

Denoised USCT with k-means algorithm combined with morphologic algorithm and Otsu algorithms: (g) cortical bone detection, (h) cortical bone detection with enhanced PSNR, (i) medullary cavity detection.

5.2 PSNR results

The PSNR results provided in Table 1 demonstrate the best result with a value of 25.15 with an implemented Haar wavelet. However, a value of 13.08 is recorded with k-means combined with the morphologic algorithms and Otsu method. In general, image areas with higher PSNR and CNR estimates indicate having better contrast resolution [23]. The PSNR is important as it is a good measure of image quality. In detecting lesions in the body, however, a high PSNR alone will not guarantee that sufficient contrast exists to make the lesion detectable. The CNR between the lesion and background is important as it serves as a quantitative metric for low-contrast lesion detection: The higher the CNR between the lesion and the background, the more likely the lesion detection [24, 25].

Algorithm	Image	Size	PSNR
—	Source	256*256	17.87
Haar Wavelet (Figures 6, and 7)	Image LL	64*64	16.18
Proposed Algorithm (Figures 8–10)	Image LH	64*64	24.65
	Image HL	64*64	25.15
	Image HH; Image LH1; Image HL1; Image HH1	64*64; 128*128; 128*128; 128*128	22.98; 22.04; 21.38; 19.15
	Image (a); Image (b); Image (c); Image (d); Image (e); Image (f); Image (g); Image (h)	256*256; 256*256; 256*256; 256*256; 256*256; 256*256; 256*256; 256*256	10.62; 10.94; 11.32; 13.08; 13.11; 13.14; 11.36; 11.42

Table 1.

PSNR results.

5.3 Time process

As shown in Figure 11, the time execution process is about 33 s while time execution with k-means combined with the Otsu and morphologic algorithms is about 1.8 s.

Figure 11.

Time process.

6. Discussions

6.1 Comparative study with related work

A comparative study with the state of the art is depicted in Table 2. Our implemented Haar wavelet shows a PSNR enhancement of 7.27 dB compared to the original image. Compared to that found by [2, 5], we have gained an enhancement with a result of 43%, outperforming the related work. In the discrete wavelet transform, one approximation can be further split, as it will remove the noise from images, but for higher noise, it will lose the details and create irregularities in the edges [26, 27]. In this context, we have used the 2D Haar wavelet decomposition to detect osteopathologies, we have to detect edges for distance measurements between bones. As a result, we have improved the PSNR and we have detected external edges, but the fracture has not been detected. With our hybrid proposed algorithms, we have detected it from the bone image. Consequently, we have outperformed [2, 5] detecting pathologies from USCT images. In 2019 [3], the authors achieved a USCT PSNR value of 21.17 dB. We outperform them with a PSNR enhancement value of 4 dB. In addition, our denoised images represent a high USCT image quality and we have produced a free USCT database, given the difficulty to have a big data of USCT images.

Parameters	Our Work	[5] 2017	[3] 2019
PSNR	25.15	low	21.17
Edge detection; Region detection; Diagnostic detection	++++; ++++; + ++++	------; -----; ----; -----	+++ +; +++; ++; ++

Table 2.

Comparative study with related works.

6.1.1 Edge detection

If we assimilate the method of the Haar wavelet transform with visual C++ with the sliding window algorithm using MATLB [11, 12], we can say that we have made a net automatic edge detection and we have improved the resolution. The increase in the USCT imaging process utilizing only signal processing cannot be isolated. The way to enhance the resolution of the image [5], we can detect some bone abnormalities. First, the associated signals are treated before processing the image reconstruction. Afterward, the automatic image processing tool and precision are used. This will be the right way for the method perspective.

6.1.2 Time execution with visual C++ and python

We can say that we have reduced the time of execution of the image to 1.8 seconds. The implemented Haar wavelet algorithm is faster than the one proposed by Lasaygues in 2006. Comparing our time execution with the run time using Matlab. We have saved and gained with about 5.45 times. However, with k-means, Python runtime is about 1.8 seconds. Thus, our proposed k-means algorithm is the best method to save time compared to that found with the algorithm of the Haar wavelet. On the other hand, we must think of a deep learning algorithm that combines the neural network and the genetic algorithm to help accelerate the time execution by its implementation on a Graphic Processor Unit (GPU) system.

6.1.3 Region detection

Comparing our work done with the suggested algorithm and the method of classification with the Laboratory of Mechanic and Acousttics (LMA) work [5, 12] done with the sliding window, the present error to determine the bone cortical area is roughly between 4% and 20%. It is very high. However, our proposed algorithm aims to detect different regions of interest, as depicted in Figures 9 and 10.

6.1.4 Results validations

A comparative study with ground truth: To validate our results, a comparative study has been done with an echographic image of human bones. As illustrated in Figure 12, multiple structures of an ultrasonic bone have not been detected because of the issue of echographic frequencies to perforate bones.

Figure 12.

Echographic bone image segmentation results: (a) Echographic bone image, (b) segmented image.

In the ground truth, the used bone to be imaged is a real adult Padilla bone with a fracture of three levels (H1, H2 and H3) as depicted in Figure 5. The lengths of these fracture diameters are approximately similar to those lengths found with our proposed algorithm. Moreover, the width of the cancellous and cortical bone as well as the diameter of the medullary cavity is similar to the ground truth width with an error of + − 1 mm. Indeed, this similarity of diameter measurements indicates the performance of our used method.

Our implemented Haar wavelet using the FMC interface and Visual C++ can be applied to various tomographic images, such as cerebral computed tomography (CT) and retinal vessel computed tomography for edge and pathology detection. Figure 13 shows the performance of our implemented Haar wavelet algorithm.

Figure 13.

Vessels retinal detection in CT image.

Our k-means combined with the morphologic and the Otsu algorithms have resolved the problem of adult bone region detection and outperformed the wavelet signal processing method [2, 5]. Furthermore, it has detected osteopathologies such as fractures with its highest ability to remove noise and save time.

7. Conclusion

In this chapter, an improved hardware/software design has been used. In the beginning, we have presented the USCT prototype to avoid any x-ray exposure to human beings. Thus, we have had USCT noisy images. As a solution, we have implemented the Haar wavelet and improved k-means combined with the morphologic and Otsu algorithms. It is a very crucial method that has solved the problem of noise in USCT images. Therefore, with this approach, the resolution has been improved, the time execution has been reduced, the noise has been removed and we have gained automatic diagnosis detection. In fact, we have to say that compared to what has been already published, the method presents more resolved results. With unsupervised developed k-means, the USCT image is segmented into three anatomic regions of interest, namely, the cortical and cancellous bone and the medullary cavity. These regions will be the labeled features of bone tomographic images. The next step will be devoted to the combination of unsupervised deep learning and supervised deep learning to automatically classify our image data into two classes (pathologic bone images and healthy bone images) and for the extraction of the region of interest from big data of bone tomographic images. Besides, interest will be given to the reconstruction of 3D images.

Funding

This study was funded by the Ministry of High Education and Scientific Research in Tunisia and the Laboratory of Mechanics and Acoustics in Marseille.

Conflict of interest

Authors declare that they have no conflict of interest.

Ethical Approval

This chapter does not contain any studies with human participants performed by any of the authors.

Informed consent

This chapter does not contain patient data.