GA-based multi-objective optimization technique for medical image denoising in wavelet domain

Abstract

Medical images that are acquired with reduced radiation exposure or following the administration of imaging agents with a low dose, are often known to experience problems by the noise stemming from acquisition hardware as well as psychological sources. This noise can adversely affect the quality of diagnosis, but also prevent practitioners from computing quantitative functional information. With a view to overcoming these challenges, the current paper puts forward optimization of multi-objective for denoising medical images within the wavelet domain. This proposed technique entails the use of genetic algorithm (GA) to get the threshold optimized within the denoising framework of wavelets. Two purposes are associated with this technique: First, its ability to adapt with different noise types of noise in the image without requiring prior information about the imaging process per se. In addition, it balances relevant diagnostic details’ preservation against the reduction of noise by considering the optimization of the error factor of Liu and SNR as the foundation of objective function. According to the implementation of this method on magnetic resonance (MR) and ultrasound (US) images of the brain, a better performance has been observed as compared to the existing wavelet-based denoising methods with regard to quantitative and qualitative metrics.

Keywords

Medical image denoising rician noise speckle noise wavelet thresholding threshold optimization optimization techniques multi-objective optimization

1 Introduction

Oftentimes, deterioration is observed in medical images due to the interference from various sources and other phenomena that impacts the processes of measurement in imaging and data acquisition systems. The small discrepancies between normal and abnormal tissues often get exacerbated by noise, which sometimes makes it difficult to directly analyse the images acquired [1]. Therefore, it is important to improve medical images’ quality not only to make an accurate diagnosis, but also to compute quantitative functional information. Traditional noise filtering techniques like median, mean and wiener filtering leverage pixel spatial information to remove the noise while preserving the feature information. The effectiveness of these filters rely on window size and its orientation; if not set appropriate results in blurred edges [2]. But, it is of paramount in medical images to safeguard the edge information to ascertain the difference between normal and pathogenic conditions. For this reason, noise filtering techniques that represent the signal in single-scale (frequency or time) are often unable to distinguish signals from noise [3].

In this context, methods based on multi-scale and wavelet have become increasingly popular to denoise images [4, 5]. Accordingly, noise filtering in wavelet domain is preferred compared to spatial filtering as they often blur low contrast edges. In many studies, it is proven that noise filtering in wavelet domain guarantees superior performance. Nevertheless, these methods involve threshold for noise filtering which should be selected in sync with the distribution of noise/signal component across various scale. This particularly holds true in medical images where it is equally important to depend upon resolution and SNR to ensure that the diagnosis is accurate. With a view to overcoming this problem, adaptive thresholding methods based on prior knowledge of noise are proposed in literature for different modalities of medical images. For example, In [6], the authors use generalized Nakagami distribution for modeling the speckle noise and estimate the threshold for suppressing speckles in wavelet domain. The success of these methods relies on the accuracy of noise prior being used. Unfortunately, different types of noises are associated with different modalities of medical images. There are also times when images from same modality suffer from different types of noises. To illustrate, speckle noise affects ultrasound images, whereas Rician noise corrupts the magnetic resonance images.

Towards this end, several adaptive techniques are proposed in literature to estimate threshold without prior knowledge of noise type. For example, The authors in [7] proposed a new learning approach to optimize both threshold and thresholding function simultaneously without priori knowledge of noise type. But, the method was for natural images. Whilst, the authors in [8, 9] used hidden Markov tree model (HMT) and Bayesian models to learn the correlation between wavelet coefficients to empirically estimate the noise level and determine threshold adaptively for each subband across different scales. Both these employed medical images to demonstrate the effectiveness of HMT and Bayesain model for threshold estimation. Though several denoising methods are proposed in literature for medical image denoising [10 –12] in recent years, yet the threshold estimation for wavelet denoising is still in its infancy. This paper aims to address the above described issue with heuristic optimization technique. To this end, the study first formulates the medical image noise suppression in wavelet domain as optimization problem and employs GA to search for the most optimal threshold for each wavelet subband across different scales. Later, the study defines appropriate GA operator that can best handle the optimization problem under the study. To achieve global optimal threshold value, multi-gene approach is adopted to encode the threshold value of all subbands and introduces mutli-objective optimization as fitness function integrating two conflicting function, SNR and Liu’s error factor with a view to ensuring feature preservation and noise reduction. Finally, two set of experiments were conducted to demonstrate the effectiveness of proposed GA-based wavelet method against state-of-the-art related works over medical images available in public domain. The results demonstrated the superior and consistent denoising performance of the proposed method on two different types of medical images.

2 Problem formulation

The denoising problem aims is to estimate signal f from noisy observations f′ so as to preserve critical features. Multiplicative and additive noise are known to typically corrupt medical images. The modeling of the image f′ (x, y) corrupted with additive noise is given as below: $f^{'} (x, y) = f (x, y) + n (x, y)$ (1) Similarly, the modeling of the image corrupted with multiplicative noise f′ (x, y) is given as follows,

$f^{'} (x, y) = f (x, y) [1 + n (x, y)]$ (2) Here the noise function n (x, y) is used to generate random value with arbitrary distribution. For the past few decades, the application of wavelet transform (WT) for non-linear filtering called denoising, has received increased attention across researchers in the field of medical image analysis [12, 13]. This increase in interest is mainly owing to the following reasons: 1) simplicity of the approach, 2) computational complexity of the approach is reduced, and 3) approach enables image analysis at different frequency bands that are associated with different image feature details. In general, recursive decomposition is performed the discrete WT on lower frequency band obtained from previous decomposition. The recursive decomposition results in left recursive tree. The application of DWT to get low (approximation) and high (detail) frequency components at level j is written as follows mathematically [13],

$\begin{matrix} a^{j + 1} (x, y) = ((a^{j} * h) ↓ 2) (x, y) \\ d^{j + 1} (x, y) = ((a^{j} * g) ↓ 2) (x, y) \end{matrix}$ (3) Here g and h denote the highpass and lowpass filter. Also, the operators * and ↓ 2 denote discrete convolution and downsampling by 2. Using the same notation, the inverse wavelet transform can be written as follows,

$\begin{matrix} a^{j} (x, y) & = & ((a^{j + 1} ↑ 2) * h) \\ + ((d^{j + 1} ↑ 2) * g) (x, y) \end{matrix}$ (4) Within the context of denoising, the wavelet’s decorrelation property suggests to process the coefficients independent of each other; whereas the property of sparseness makes it possible to utilize thresholds for attenuating or removing the wavelet coefficients smaller than the level of noise. This process is termed as wavelet thresholding as it shrinks the given coefficients to zero. Wavelet coefficients can be shrunk by hard thresholding and soft thresholding as proposed by Donoho et al. [14]. Under the former, wavelet coefficients larger than the threshold λ is retained. On the other hand, the one which is smaller or equal to λ are set to zero as shown in Equation. 5.

$T_{hard} (w) = {\begin{matrix} w & | w | > λ \\ 0 & | w | \leq λ \end{matrix}$ (5) The soft thresholding scheme is formulated as follows,

$T_{soft} (w) = {\begin{matrix} w - λ & | w | > λ \\ 0 & | w | \leq λ \end{matrix}$ (6) A comparison of both schemes reveal that hard thresholding displays discontinuities at ±λ and could show greater instability or sensitivity toward small data changes. By contrast, soft thresholding tends to have more stability and fewer discontinuities. Despite the inherent simplicity and efficacy of thresholding, choosing a proper value of threshold is challenging since it has a direct impact on the denoising performance. As a case in point, selecting a massive threshold would shrink nearly all the coefficients, potentially leading to excessive smoothing as far as diagnostic features are concerned. Meanwhile, a smaller threshold may be unable to suppress artifacts or noise despite retaining sharp edges. In previous literature, Bayesian methods have been used to address wavelet thresholding with the use of prior probability distribution to formally model the relationship between DWT co-efficient, thus producing posterior distributions. It can then be possible to derive estimation rules using the conventional Bayesian approach, after specifying a function loss. Typically, Bayesian techniques have proved to gain better performance than state-of-the-art denoising methods in wavelet domain [8 , 15–17]. But, the reliance of these methods to a particular type of noise and the manual parameter settings of probability distribution function reduces their flexibility for processing medical images. This implies that as far as denoising medical images is concerned, application of intelligent optimization for threshold selection is of paramount importance.

3 GA-based optimal solution for threshold selection

GA is known to be one of the most frequently used techniques of artificial intelligence for purposes of optimization. GAs refer to stochastic searching algorithm premised on genetics and natural selection. Their efficacy, robustness, simplicity, and stability have been proven in complicated, non-linear, and ill-posed optimization problems in the context of image processing [18 –20]. GA has proven its ability in finding the approximate global maxima in a manner that is not predicated on a specific problem domain and also without requiring in-depth knowledge on the technique of processing. For example, [21], the authors have developed a GA optimized model-based tomographic reconstruction to detect and estimate multiple circular and elliptical objects of known intensities with limited number of noisy projections. Another work employs GA to develop an automatic medical image segmentation problem [22]. Likewise, near field effect of X-ray source was modelled for retrospective correction and the parameters of the model were estimated using GA [23]. Due to the success of these prior studies, the current research adopts the GA approach to undertake optimization of threshold in the context of wavelet denoising. The implementation of GA operators for the present study is described in the subsection following below,

3.1 Chromosome encoding

This study adopts a multi-gene approach for encoding the values of threshold that need optimization [18]. All chromosomes comprise 3 genes, each of which entails 3 sub-genes as illustrated in Fig. 1. It is evident that these 3 genes correspond to wavelet scales, with sub-genes denoting p_bl with 1 ≤ b ≤ 3 and 1 ≤ l ≤ 3. Here p_bl takes real value ranging between 0 and 1. Notably, the value denotes the optimum threshold that is assigned at scale l to subband b. A binary representation of string is used for encoding the float value of the sub-gene. Mathematically, the float value is computed as follows from the sub-gene binary string, $V = (- 1)^{s} . 2^{E} . (M)$ (7) Here, M and E represents mantissa part as fixed point real number and Exponent part as signed integer respectively. Also, S represents sign bit with 0 or 1 value.

Fig. 1

Chromosome encoding for GA-based wavelet thresholding.

3.2 Multi-objective fitness function

In real-world practice, there are class of problems that demand for simultaneous optimization of at least two conflicting objectives. These class of problems are known as multi-objective optimization problem (MoP). Unlike an optimizer of single objective optimization problem, an effective optimizer for MoP aims to find set of optimal solutions that can best reflect the tradeoff between all conflicting objective functions simultaneously. Over the past decades, several approaches have been proposed in literature to solve MoP. The well-known and highly suitable practical approach for solving MoP is weighted sum method. This method aggregates all objective functions to form single objective function and treats MoP as a composite of objective functions as follows, $MOO (x) = \sum_{i = 1}^{n} W_{i} F_{i} (x)$ (8) where F_i are different objective functions and W_i is the assigned weights to the objective functions. The key benefits of weighted sum method are its application simplicity and enabling to tune the significance of the objective functions according to the system requirement. Motivated by these benefits, the method is adopted in this study.

As discussed in introduction section, most of the medical image denoising algorithms reported in literature aims to suppress noise without emphasizing to preserve the diagnostic details. In real practice, during denoising process there is a tradeoff between noise removal and diagnostic feature preservation. To fill this gap, the fitness function of our GA algorithm is formulated as multi-objective optimization considering the two conflicting objective functions, SNR for noise removal and Liu error factor to preserve the essential diagnostic features. Leveraging the benefits of weighted-sum method, these two conflicting objective functions are aggregated to define the fitness function and is expressed as follows

$fit (p) = \max (SNR (f_{enh} (p)) - F (f_{enh} (p))$ (9) where

$\begin{matrix} SNR (f_{enh}) & = & \frac{{ROI}_{signal} (f_{enh})}{{ROI}_{noise} (f_{enh})} \\ F (f_{enh}) & = & \frac{1}{1000 (NM)} \sqrt{R} \sum_{i = 1}^{R} \frac{e_{i}^{2}}{\sqrt{A_{i}}} \\ f_{enh}^{[l, b]} & = & W^{- 1} (f_{wd}^{[l, b]}) \\ f_{wd}^{[l, b]} & = & T_{soft} (f_{w}^{[l, b]}, p^{[l, b]}) \\ f_{w}^{[l, b]} & = & W (f^{'}) \end{matrix}$

Here, f and f_enh represents the noisy input image and the denoised output image respectively. In the same way, the DWT and inverse DWT operators are denoted by W and W^-1 respectively; Also, the wavelet and thresholded coefficients in subband b at scale l is denoted by $f_{w}^{[l, b]}$ and $f_{wd}^{[l, b]}$ respectively. The soft thresholding operator on a wavelet coefficient $f_{w}^{[l, b]}$ using the threshold, p^[l,b] at subband b in the scale l is denoted by $T_{soft} (f_{w}^{[l, b]}, p^{[l, b]})$ .

Here, the visibility of ROI and the clarity of feature details is of paramount importance for medical diagnosis. Hence, the work here employs soft thresholding as it enables to achieve visual quality without abrupt artifacts, compared to hard thresholding.

3.3 Crossover

This operator entails picking chromosomes’ pairs randomly from the existing population to implement single-point crossover operation if the generated random number is smaller than the predefined crossover rate (rx) (crossover at the level of chromosome). As illustrated in Fig. 2, a random number is generated for each gene within the chromosome and the exchange of two genes takes place when the value of this random number is lower than a predefined crossover rate of genes (rgx). Following the crossover, all chromosomes are normalized for meeting the constraint of gene-value. Our experiments revealed that the proposed GA approach worked effectively when r_gx = r_x = 0.7 . Mathematically, the crossover procedure is represented as follows

$\begin{matrix} {Parent}_{1} & = [P_{m 1}, P_{m 2}, . . . ., P_{mc} . . . . P_{mN}], \\ {Parent}_{2} & = [P_{d 1}, P_{d 2}, . . . . P_{dc} . . . . P_{dN}], \\ P_{new 1} & = p_{mc} - rgx [p_{mc} - p_{dc}] \\ P_{new 2} & = p_{dc} + rgx [p_{mc} - p_{dc}] \\ {offspring}_{1} & = [P_{m 1}, P_{m 2}, . . . ., P_{new 1} . . . . P_{mN}] \\ {offspring}_{2} & = [P_{d 1}, P_{d 2}, . . . . P_{new 2} . . . . P_{dN}] \end{matrix}$ (10)

Fig. 2

Crossover operator of GA-based wavelet thresholding.

3.4 Mutation

After taking the aforementioned steps, the study subjects the chromosomes to mutation operation. All sub-genes were found to contain 32-bit floating point number in string array form as illustrated in Fig. 3. In the sub-genes, bits selected randomly were flipped. Two schemes of mutation were tested: two-points mutation and uniform mutation. Consideration was given to two parameters for both these schemes: mutation rate at the level of gene (r_gm) and mutation rate at the level of chromosome (r_m). if $P_{k 1}^{'}$ is mutated value, then mutation operation is given as follows, $P_{11}^{'} = (P_{11} + P_{k 1}^{'}) - P_{k 1}^{'}$ (11)

Fig. 3

Mutation operator of GA-based wavelet thresholding.

The application of a mutation operator was made only if the generated random number is lower than (r_m). The two-point mutation scheme randomly selects two genes but performs mutation over sub-gene in one of the gene by flipping each bits in the mantissa parts. Also, the value of sub-gene in another gene is revised. Whereas, under the scheme of uniform mutation, gene-level mutation of all genes within a chromosome was undertaken in a way that bears resemblance with the mutation scheme with two-points. Initially, attempts were made to mutate with probabilities of bit exchange ranging from 0.0 to 0.1. Improvements in convergence and better performance was observed for a two-point mutation whose probability was (r_m) = (r_gm) =0.003.

3.5 Stopping condition

This study employs the GA approach substantiating three conditions for stopping the optimization process. The first one occurs when the total iterations attain predetermined iterations; when not much change is observed in each generation’s fittest chromosome, and when the fitness values of all chromosomes remain the same, that is, after the convergence of the algorithm. The pseudocode of the proposed GA-approach for medical image denoising with optimal selection of wavelet thresholds for each subband across scales is given in Table-1.

Table 1
Pseudocode for GA optimization of wavelet thresholding

1. Parameter initialization of GA operators

1.a. mutation and crossover probability

1.b. create initial population randomly

2. Repeat

for each chromosome in the population do

Compute the chromosome fitness using Equation 7

Compute probability of the selection

Select parents and

Do crossover to create offspring

Mutate offspring

end for

Replace the old population by the new one

until ((SNR (f_enh) - F (f_enh)) is optimized)

3. Apply DWT on the image to obtain the four subbands

4. Apply Thresholding on all subbands except LL subband

using GA optimized threshold values

5. Obtain denoised image applying Inverse DWT

4 Dataset description

The performance and effectiveness of the new GA-based denoising method was compared with other state-of-the-art work for two types medical images, US and MR. For this purpose, we use two different benchmark datasets that contains real patient data namely, MICCAI BRATS [24] and Medpix [25].

As all the experiments designed in our study used these two datasets that are available in public domain, there was no need for ethical clearance. The US images used in this study were taken from Medpix online database. It is a radiology database archiving images from different medical imaging modalities, clinical case reports and teaching cases collected from 12,000 patients. Similarly, the MR brain scans were acquired from the benchmark dataset provided as part of MICCAI challenge. It contains multimodal MR brain scans for 285 cases collected from multiple institutions. For each case, the dataset includes four MR sequences namely, T1, T2, T1c, T2 and FLAIR with 259 high-grade and 76 low-grade Glimoas cases. For our experiment and comparative evaluation, we choose randomly two US and MR images from these datasets. The reference US and MR images used for our evaluation are shown in Fig. 4 (Row-1) and 4 (Row-2) respectively.

Fig. 4

The reference US (Row-1) and MR (Row-2) Images.

5 Results and discussion

This section presents and discuss the experimental results of the method proposed in this study on US and then with MR image in comparison against existing state-of-the-art works. To do this, we considered HMT approach of Crouse et al. [8] and Bayesian approach of Pizurica et al. [9]. The main reason for considered these works for comparison is that, to our knowledge, there are few works in the literature that leverage the benefits of wavelet coefficient correlations to learn the noise level in the underlying image, as we do in our work.

For comparison purposes, the implementation of these methods was premised on the authors’ algorithm given in their respective published papers. The achieved performance was also quantified by computing quality metrics. Another degree of flexibility of denoising in wavelet domain is the choices of the underlying wavelets. Therein, we apply Daubechies 4 (DB4) wavelet in all the evaluation experiments to preserve the fine details of image. The four vanishing moments in DB4 wavelet provides better approximation and frequency localization; thereby resulting in good performance.

All the experiments of the study were conducted on a personal computer with Windows 10 operating system and following specifications, Intel Core i7-8565H @ 1.8 GHz, 128 GB RAM. The algorithm of the proposed method was implemented using MatLab 2019a.

5.1 Application to US image

5.1.1 Speckle noise modeling

To validate the efficacy of the proposed method for denoising, low-dose modeled US image was simulated by corrupting the reference image shown in Fig. 4 (Row-1) with speckle noise. The spatially correlated speckle noises were generated according to the pseudocode given in Table-2, by lowpass filtering the reference image with unit-mean Gaussian random fields and considering the magnitude of the denoised image. The speckle’s correlation length is then controlled by determining the kernel’s size for the purpose of introducing correlation in the Gaussian noise. For modeling the reality, it is believed that a correlation derived with a size 3 kernel is typically sufficient. The variance of Gaussian random fields was changed to enable simulate speckle noise of different levels. Accordingly, several images were simulated varying the noise variance for experimental study. The US images speckled with noise variances of 0.03, 0.07, 0.1 and 0.2 were simulated as shown in Row-1 of Figs. 5 and 6 respectively. The noise variances were selected in such a manner that the resulting noisy image is analogous to the pSNR level experience in real environment practice.

Table 2
Pseudocode for speckle noise modeling [5]

1. Initialize noise parameter σ

2. Generate normalized random number

using built-in function randn()

3. Simulate speckle noise as defined below,

G(x,y) = randn(x,y)σ+sqrt(-1)randn(x,y)*σ

S(x,y) = lowpass_filter(G(x,y),kernel_size)

N(x,y)=normalize(abs(S(x,y)))

4. The reference image is speckled with noise

simulated in step 3 using Equation (2)

Fig. 5

Visual performance comparison of proposed method against related works on reference US image (left) in Fig. 4 (Row-1) for different speckle noise levels.

Fig. 6

Visual performance comparison of proposed method against related works on reference US image (right) in Fig. 4 (Row-1) for different speckle noise levels.

5.1.2 Comparative study

As informed earlier, the effectiveness of the new denoising technique is first tested with US images collected from Medpix database. The study considers two US images for comparative evaluation. These images as shown in Fig. 4 Row-1 depict coronal and sagittal view of US with dilated ventricles with evidence of hemorrhage in right lateral ventricle near caudothalamic groove. The first set of experiment aims to compare the denoising performance of the GA-based approach against HMT and Bayesian methods on US images. For this purpose, as detailed in previous subsection, the US images with simulated speckle noise for four different noise levels 0.03, 0.07, 0.1, 0.2 as shown in Row-1 of Figs. 5 and 6 were used.

Figs. 5 and 6 illustrates the visual performance of the proposed method for different noise levels column-wise against related works, HMT and Bayesian methods row-wise. The column-wise visual inspection of the results through different noise levels on two US images indicate that the noise suppression with diagnostic feature preservation decreases as data noise level increases, regardless of the denoising methods considered for comparison. This indicates that the increase in noise level affects the performance of all denoising methods. The row-wise comparison of the results through different denoising methods on two US images shows that HMT based wavelet denoising performs worse with distortion of ROIs compared to other two counterparts. On other hand, the Bayesian method show comparatively better results than that of HMT when it comes to noise suppression and feature preservation. Even though Bayesian method resolves ventricles and hemorrhage location, the edges are not delineated clearly.

The notable observation is that the proposed model shows reasonably better performs than the other wavelet denoising methods with more robust background noise suppression at the same time it retains sharp intensity variations in resolving the ventricles, and brings out the hemorrhage especially at higher noise level. This might be because the GA-based method enables to find the optimal threshold and suppress the noise adaptively with feature preservation in each band across different wavelet scales. On the whole, the visual comparison of results in Figs. 5 and 6, it is evident that the performance of wavelet denoising methods can be improved with diagnostic feature preservation especially when noise levels are higher if optimal thresholds can be determined for each wavelet band of different scales.

5.2 Application to MR image

5.2.1 Rician noise modeling

MR images obtained by Fourier reconstruction are shown to be sensitive to thermal noise. To corrupt the reference MR images, the additive white Gaussian noise is added to both real and imaginary channels according to the pseudocode given in Table-3. Noisy MR images required for our experiments were generated corrupting the reference images shown in Fig. 4 (Row-2) with complex zero mean white Gaussian noise. The variance of Gaussian random fields was changed to enable simulate rician noise of different levels. Accordingly, several MR images were simulated varying the noise level variance for experimental study. The MR images with noise variances of 0.07, 0.15, 0.30 and 0.5 were simulated as shown in Row-1 of Figs. 7 and 8 respectively.

Table 3
Pseudocode for rician noise modeling [5]

1. Initialize noise parameter σ

2. Generate normalized random number

using built-in function randn()

3. Simulate rician noise as defined below,

G(x,y) = randn(x,y)σ+sqrt(-1)randn(x,y)*σ

4. The reference image is corrupted with noise

simulated in step 3 using Equation (1)

Fig. 7

Visual performance comparison of proposed method against related works on reference MR image (left) in Fig. 4 (Row-2) for different rician noise levels.

Fig. 8

Visual performance comparison of proposed method against related works on reference MR image (right) in Fig. 4 (Row-2) for different rician noise levels.

5.2.2 Comparative study

The performance of the proposed method on US images encouraged us to evaluate its effectiveness on MR images. As described earlier, two MR images collected from MICCAI BRATS dataset were considered for comparative evaluation. Here, the first reference MR brain image was chosen with no evidence for hemorrhage but with extensive white matter lesions and no contrast enhancement. The second reference MR brain image was chosen with evidence for hemorrhage with enhancing mass lesion and associated edema. These two MR reference images are shown in Row-2 of Fig. 4. The second set of experiment aims to compare the denoising performance of the GA-based approach against HMT and Bayesian methods on MR images. For this purpose, the MR images with simulated rician noise for four different noise levels 0.07, 0.15, 0.3 and 0.5 as shown in Row-I of Figs. 7 and 8 were used.

The visual performance of the proposed method on MR images for different level of rician noise column-wise against the related works row-wise are presented in Figs. 7 and 8. The column-wise visual inspection of the results through different noise levels on two MR images are similar to that observed on US images indicating that the denoising performance degrades with increase in noise level regardless of the denoising methods considered for comparison. This confirms that high level noise greatly influences the denoising performance of a denoising methods. The row-wise comparison of the results through different denoising methods on two MR images also concurs with the results obtained on US images. For example, the resultant image from HMT based wavelet denoising reveals that HMT based wavelet methods fail to learn the noise type on MR images and shows either blurring on edges or ROI distortion. While, the Bayesian method removes noise and preserves the ROI better quality as compared to the HMT method. Yet, the homogenous regions show noise patches. The proposed method demonstrates to remove noise preserving the sharpness of image better compared to other competing methods. Also, the proposed method performs better in delineating the adjacent regions with similar intensity. This characteristic of the proposed approach might be due to the use of weighted sum method to optimize two conflicting parameters viz, SNR and resolution factor.

5.3 Quantitative performance evaluation

Quantitative or subjective analysis can be used to evaluate the enhanced image’s quality. The latter is premised and influenced by varying perception of humans. Unfortunately, mathematical formulae lack objectivity when it comes to realistically replicating the visual system of humans. As a result, the quality of enhancement hinges on factors like spatial compactness, homogeneity, and continuity, among others. This implies that it is highly improbably for one measure to meaningfully capture these factors. Thus, it is necessary to evaluate the goodness based on the quality of image enhancement in the field of interest application. However, mean square error (MSE) and Peak SNR (PSNR) are widely employed for measuring errors quantitatively.

In our work, three measures MSE, PSNR and structural similarity (SSIM) were calculated for the denoising results obtained from three different denoising methods, HMT, Bayesian and GA on two US and MRI images. The SSIM is considered to evaluate the effectiveness of denoising method for noise suppression against feature preservation. These measures are computed as follows, $MSE = \frac{1}{NxM} \sum_{i = 1}^{M} \sum_{j = 1}^{N} [f_{enh} (x, y) - f^{'} (x, y)]^{2}$ (12) Here, f_enh (x, y) and f′ (x, y) signify the enhanced and input noisy image of size N x N respectively; x and y denote the image pixel position along horizontal and vertical directions. The metric, PSNR was computed using, $PSNR = 10 \log [\frac{(\max (f^{'}))^{2}}{MSE}]$ (13) Here, max (f′) denotes the maximum pixel value in the image f′ because it is mapped to the maximum intensity in the acquired input image. $SSIM = \frac{(2 μ_{f} μ f^{'} + C_{1}) (2 σ_{f} σ f^{'} + C_{2})}{(μ_{f}^{2} μ {f^{'}}^{2} + C_{1}) (σ_{f}^{2} σ {f^{'}}^{2} + C_{2})}$ (14) Here, μ (f) andμf′ denotes the average of the image f′ and f. Similarly, σ (f) andσf′ denotes the standard deviation of the image f′ and f. The two variables C₁andC₂ are employed to stabilize the division with weak denominator;

In this step of analysis, the above discussed performance metrics such as MSE, PSNR and SSIM are computed for all resultant images presented in Figs. 5-8 and reported in Tables-4-7 respectively. Also, for clarity the best results for each dataset are highlighted in boldface.

Table 4

Quantitative Comparison of proposed method against related works for resultant US images in Fig. 5 in terms of (MSE/PSNR/SSIM)

Table 5

Quantitative Comparison of proposed method against related works for resultant US images in Fig. 6 in terms of (MSE/PSNR/SSIM)

Table 6

Quantitative Comparison of proposed method against related works for resultant MR images in Fig. 7 in terms of (MSE/PSNR/SSIM)

Table 7

Quantitative Comparison of proposed method against related works for resultant MR images in Fig. 8 in terms of (MSE/PSNR/SSIM)

From these results, it is worth noticing that HMT-based method demonstrates marginal performance under low noise level but degrades with increase in noise level displaying poor and discouraging performance on both US and MR images. In many applications, Bayesian-based methods have proved their outstanding performance, but it is surprise to note their average performance in this study. The reason might be because the success of these methods relies on prior knowledge about the noise model. In this study, Bayesian method was designed to model speckle and gaussian white noise for US and MR images respectively even then the assumption did not support Bayesian method to identify and suppress the noise with feature preservation under high noise level.

Inspection of the results remark that the proposed GA-based method yields best noise suppression performance with feature preservation ability on US and MR images with regard to all the compared metrics. Specifically, it can be observed that GA-based methods attain the best SSIM and higher than average performance in terms of PSNR and MSE on both US and MR datasets. Also, we can see the inconsistent performance of Bayesian-based methods across medical images. For instance, it shows comparably better performance on training set on US images. But, surprisingly it delivers the penalized performance on MR images.

On other hand, it is interesting to observe the proposed demonstrating consistent denoising performance on both US and MR images. The success of the proposed method may be due to the introduction multi-objective optimization involving SNR for noise suppression and Liu’s error factor for diagnostic feature preservation maintains the trade-off between the two conflicting objective function. This enables the proposed method to demonstrate consistent denoising performance finding the most optimal threshold for noise suppression with feature preservation without inducing blur.

The obtained quantitative metrics are also in agreement with qualitative images indicating that the proposed GA-based method is comparably efficient in achieving promising denoising results in terms of both PSNR and SSIM on both US and MR images. Thus, it is evident that the proposed GA-based method removes noise from medical images in an intelligent manner by learning the type of noise without necessitating any previous knowledge.

5.4 Computational efficiency

The time complexity of the proposed method was examined by calculating the arithmetic operations needed to implement the wavelet transform and GA operation. In this context, it is notable that the computational complexity of wavelet operation depends on image size, N and are denoted by O (N²). Whereas, the computational complexity of GA operation depends on population size, n and is denoted by O (n²). Thus, for a given image size, N and population size, n, the overall computational complexity of the proposed method approximates to O (N²)

For comparison purposes, the computation time taken for images of size (128,128), (256,256) and (512,512) by GA, Bayesian and HMT are shown in Fig. 9. It can be noted that the computation time taken by HMT and Bayesian increases with increase in image size. Therein, image of size (128,128) was considered for evaluation in all the experiments to ensure endurable computation cost among all the comparing methods. The computation time of HMT includes time taken for both training and denoising process. Therefore, the trained HMT-based method ensure to take less time for denoising. On contrary, both Bayesian- and GA-based methods process each image individually to understand the noise parameters of the image. For this reason, computation time of the GA- and Bayesian methods are higher as compared to the method based on HMT. Therein, future studies will aim to address GA’s computational expense by developing parallel algorithms that can be executed on high-speed computing resources in parallel.

Fig. 9

Time complexity of wavelet denoising methods.

6 Study limitations

The present study has several limitations that demands future investigations. First, the proposed noise reduction method is evaluated on medical images corrupted with simulated noise due to ethical issues. Therefore, evaluation of the proposed method on small animal data need to be verified to ensure its ability to process real low-dose patient data. Second, performance comparison has been conducted with state-of-the-art noise reduction methods. But, application of different heuristic optimization techniques for wavelet threshold determination has not been investigated. Third, the control parameter of GA were chosen in a trial and error manner which may vary for different modality images. Therefore, future work can explore on approach that can adaptively select the most optimal control parameters for GA.

7 Conclusions

This study introduces a novel, effective, and robust technique of denoising images in the wavelet domain based on GA. The study examines its performance using US and MR images before making comparisons with other approaches of wavelet premised on HMT and Bayesian models. According to the findings, the proposed method displays superior denoising performance with feature preservation compared to other existing works in terms of MSE, PSNR and SSIM even at high noise variance level of 0.5, thus reinforcing its ability to lower noise in not merely MR and ultrasound images, but also across other tomographic images. The interplay between resolution and SNR is another outcome of the study. This technique is unique in that it provides a plausible means of getting noise suppressed on different portions of an image, and no previous knowledge about the noise is required.

Footnotes

Acknowledgment

The authors are very grateful to thank their Institution, Prince Sattam Bin Abdulaziz University at AlKharj in Saudi Arabia for providing financial support and extend a special thanks to their Deanship of Scientific Research for technical resource support in accomplishing this work successfully.

References

Vaiyapuri

, Alaskar

, Whale optimization for wavelet-based unsupervised medical image segmentation: Application to ct and mr images, International Journal of Computational Intelligence Systems 13(1) (2020), 941–953.

Mateoa

J.L.

, Caballero

A.F.

, Finding out general tendencies in speckle noise reduction in ultrasound images, Expert Systems with Applications 36 (2009), 7786–7797.

Vaiyapuri

, A new automated approach for early lung cancer detection with improved diagnostic performance–a preliminary analysis, Indian Journal of Forensic Medicine & Toxicology 14(2), 2020.

Unser

, Aldroubi

, Laine

, Guest editorial: Wavelets in medical imaging, IEEE Trans on Medical Imaging 22 (2003), 285–288.

Thavavel

, Basha

J.J.

, Krishna

M.C.

, Murugesan

, Heuristic wavelet approach for low-dose epr tomographic reconstruction: An applicability analysis with phantom andimaging, Expert Systems with Applications 39(5) (2012), 5717–5726.

Gupta

, Kaura

, Chauhan

R.C.

, Saxena

S.C.

, A versatile technique for visual enhancement of medical ultrasound images, Digital Signal Processing 17 (2007), 542–560.

Gupta

, Kaura

, Chauhan

R.C.

, Saxena

S.C.

, Image denoising in the waveletdomain using a new adaptive thresholding function, Neurocomputing 72 (2009), 1012–1025.

Crouse

, Nowak

, Baraniuk

, Wavelet-based statistical signal processing using hidden markov models, IEEE Trans Signal Processing 46 (1998), 886–902.

Pizurica

, Philips

, Lemahieu

, Acheroy

, A versatile wavelet domain noise filtration technique for medical imaging, IEEE Trans on Med Imag 22 (2003), 323–331.

10.

Miri

, Sharifian

, Rashidi

, Ghods

, Medical image denoising based on 2d discrete cosine transform via ant colony optimization, Optik 156 (2018), 938–948.

11.

Elhoseny

, Shankar

, Optimal bilateral filter and convolutional neural network based denoising method of medical image measurements, Measurement 143 (2019), 125–135.

12.

Liu

, El Basha

M.D.

, Li

, Xiao

, Sanelli

P.C.

and Fang

, Deep evolutionary networks with expedited genetic algorithms for medical image denoising, Medical image analysis 54 (2019), 306–315.

13.

Subramanyam

V.P.

, Enhancement of ultrasound images using stochastic resonance-based wavelet transform, Computerized Medical Imaging and Graphics 32 (2008), 316–320.

14.

Donoho

D.L.

, Johnstone

I.M.

, De-noising by soft-thresholding, IEEE Trans Inform Theory 41 (1995), 613–627.

15.

Sengur

, Image denoising in the wavelet domain using a new adaptive thresholding function, Neurocomputing 72 (2009), 1012–1025.

16.

Rabbani

, Vafadoost

, Wavelet based image denoising based on a mixture of laplace distributions, Iran J Sci Technol Trans B Eng 30 (2006), 711–733.

17.

Achim

, Bezerianos

, Tsakalides

, Novel bayesian multiscale method for speckle removal in medical ultrasound images, IEEE Trans on Med Imag 20 (2001), 772–783.

18.

Yan

J.G.

, Li

B.Z.

, Yu

, Optimization design based on genetic algorithm of immunity, Journal of Machine Design 9 (2002), 14–17.

19.

Lukac

, Plataniotis

K.N.

, Smolka

, Venetsanopoulos

A.N.

, Color image filtering and enhancement based on genetic algorithms, In Proc IEEE ISCAS, pages 913–916, 2004.

20.

Shao

, Barner

K.E.

, Optimization of partition-based weighted sum filters and their application to image denoising, IEEE Trans Image Process 15 (2006), 1900–1915.

21.

Bichkar

R.S.

, Ray

A.K.

, Tomographic reconstruction of circular and elliptical objects using genetic algorithm, IEEE Signal Processing Letters 5 (1998), 248–251.

22.

Lai

, Chang

, A hierarchical evolutionary algorithm for automatic medical image segmentation, Expert Systems with Applications 36 (2009), 248–259.

23.

Korurek

, Yuksel

, Iscan

, Dokur

, Olmez

, Retrospective correction of near field effect of x-ray source in radiographic images by using genetic algorithms, Expert Systems with Applications, August 2009.

24.

Menze

B.H.

, Jakab

, Bauer

, Kalpathy-Cramer

, Farahani

, Kirby

, Burren

, Porz

, Slotboom

, Wiest

, et al., The multimodal brain tumor image segmentation benchmark (brats) , IEEE transactions on medical imaging 34(10) (2014), 1993–2024. URL https://www.med.upenn.edu/sbia/brats2017/data.html.

25.

American College of Radiology et al., Medpix™: medical image database, URL: https://medpix.nlm.nih.gov/case?id=03b244b0-7f36-4889-be14-c7589182cbfd, 2019.

1.	Parameter initialization of GA operators
	1.a. mutation and crossover probability
	1.b. create initial population randomly
2.	Repeat
	for each chromosome in the population do
	Compute the chromosome fitness using Equation 7
	Compute probability of the selection
	Select parents and
	Do crossover to create offspring
	Mutate offspring
	end for
	Replace the old population by the new one
	until ((SNR (f_enh) - F (f_enh)) is optimized)
3.	Apply DWT on the image to obtain the four subbands
4.	Apply Thresholding on all subbands except LL subband
	using GA optimized threshold values
5.	Obtain denoised image applying Inverse DWT

1.	Initialize noise parameter σ
2.	Generate normalized random number
	using built-in function randn()
3.	Simulate speckle noise as defined below,
	G(x,y) = randn(x,y)σ+sqrt(-1)randn(x,y)*σ
	S(x,y) = lowpass_filter(G(x,y),kernel_size)
	N(x,y)=normalize(abs(S(x,y)))
4.	The reference image is speckled with noise
	simulated in step 3 using Equation (2)