A non-uniform quantization scheme for visualization of CT images

1.   Introduction

Modern medical science heavily depends on imaging technologies and post-processing for diagnosis and treatment planning. Imaging modalities such as ultrasounds, X-rays, computed tomography (CT), and magnetic resonance imaging (MRI) assist radiographers in clinical diagnosis. Although image acquisition techniques have become more advanced overtimes, the captured images still do not have satisfactory quality for an accurate medical interpretation. Environmental noises, acquisition techniques, sensor noise, emergency and lightening conditions, visual analysis scenarios ^[1], and various technical limitations of imaging tools are the reasons for low-quality images ^[2]. Over the past few decades, tremendous enhancement in information technology and medical devices has driven noteworthy improvements in the medical imaging domain. Some of these techniques include fuzzy gray level difference histogram equalization ^[3], intraindividual contrast enhancement for consistency in CT Images ^[4], F-shift transformation ^[5], adaptive anchored neighborhood regression ^[6], modified shark smell optimization algorithm ^[7], and Sushisen algorithm ^[8].

Computed tomography (CT) has become an essential domain in medical modalities to supplement medical ultrasonography and X-rays ^[9,10]. Researchers are mostly utilizing CT imaging datasets to create exact and productive strategies to fulfill the current need to determine illness location ^[11]. However, CT and other images obtained utilizing medical imaging techniques generally show insufficient illumination and contrast that cause weak boundaries among various soft and hard tissues ^[12,13]. Moreover, many electronic devices with limited computational potential cannot process large amounts of data and require systematic descriptors ^[14] and image enhancement techniques with lower memory consumption and higher efficiency. Therefore, there is a need for substantial efficiency and effective contrast enhancement of CT imaging data for accurate visualization and analysis.

Many algorithms have been proposed for the enhancement of CT images. Scale-space filtering works by embedding the original image into derived images for the enhancement of CT images. By performing the convolution of the source image with a Gaussian kernel of different variance, derived images are obtained ^[15]. However, edges are not appropriately preserved in scale-space filtering. To overcome the problem, Vincent et al. ^[16] employed the anisotropic diffusion filter, which can detect edges properly and provide better smoothening even in highly noisy scenarios. Kallel et al. ^[17] proposed an adaptive gamma correction algorithm using a discrete wavelet transform with a singular-value decomposition (DWT-SVD) technique for CT image enhancement.

For the energy range of medical imagining techniques, Compton scattering (CS), Rayleigh scattering (RS), and photoelectrical absorption (PE) are the events that occur when photons interact with matters ^[18]. Compton scattering information-dependent features may lower the contrast of CT images. To solve this, Chen et al. ^[18] proposed a physical parametric method based on intrinsic properties among CS and PE. Due to this association, better enhanced CT images are produced.

Different enhancement methods work for different types of medical images due to the strong impact of the characteristics of images upon enhancement measurement parameters. Agaian et al. ^[19] proposed a technique that works collectively on both CT and MRI image enhancement. In this technique, domain transform, shape-adaptive edge enhancement, and adaptive histogram equalization methods are used to enhance image quality. Another enhancement technique by Otsu et al. ^[20] sometimes decreases the contrast of essential regions and increases background noise by manipulating the gray-scale histogram values. To overcome the problems associated with histogram manipulation, Kaur et al. ^[21] proposed an unsharp masking filter-based technique for images with different types of noises. Subramani et al. ^[3] also proposed a novel adaptive fuzzy gray level difference histogram equalization technique to cope with the noise issues. After removing the ambiguities in the image using a fuzzy gray level histogram, a clipping step further restricts the unnecessary enhancement ratio to provide a balanced contrast for medical image enhancement. However, this method suffers from the data storage issue, and to overcome that, Li et al. ^[5] proposed an enhancement method that utilizes the contrast limited adaptive histogram equalization (CLAHE) technique in the F-shit transformation domain. Zhou et al. ^[7] proposed an optimization method to deal with both local and global enhancement of medical images.

Wu et al. ^[22] proposed an advanced contrast enhancement reversible data hiding technique for medical imagining. In this method, gray-scale values in the image background are recognized as a region of interest (ROI). After labeling histogram bins of gray-level values, the un-labeled bins enhance ROI, which leads to better contrast enhancement and visually appealing results. Malik et al. ^[23] proposed to use a median filter and adaptive histogram processing for better smoothing, noise removal, and enhancement of CT images. In soft tissues and background areas, traditional enhancements techniques result in overly enhanced contrast along with noise. For this, Yangdai et al. ^[24] proposed a weighted retinex algorithm to adjust the detail enhancement in the ROI. The noise around ROI areas results in poorly enhanced images. To overcome this issue, An et al. ^[25] proposed an image fusion technique based on Haar Wavelet Transform (HWT) for enhancement and filtering the noise around ROI. The decomposition of input images into multi-scaled edges in various wavelet areas using HWT is termed fusion. This method results in low noise enhanced CT images.

Most histogram-based methods are not helpful to enhance the soft tissues of CT brain images. However, by maintaining background brightness and stretching the matching grey level of soft tissues, the method proposed by Tan et al. ^[26] helps to enhance the contrast of the soft regions in medical images. To solve the edge degradation problem in enhanced images, Bhandari et al. ^[27] proposed a two-step algorithm. Initially, a discrete wavelet transform is applied to the image, while in the second step, singular-value decomposition (SVD) is used for processing the low-low (LL) sub-band. To overcome the problem of contrast over enhancement, Sundaram et al. ^[28] proposed a histogram-modified local contrast enhancement (HMLCE) technique. This method provides resultant images with better contrast enhancement and preserves local details for appropriate interpretation of mammogram images. A seed-dependent adaptive region growing method is described by Ganesan et al. ^[29] for enhancement. This technique provides better results in terms of efficiency and effectiveness. Al-Juboori et al. ^[30] proposed a method to enhance mammogram images by merging the adaptive histogram equalization technique with the retinex method.

Due to specific image acquisition environments, digital imaging and communication in medicine (DICOM) images usually suffer from luminance and contrast problems ^[31]. To obtain clear and accurate medical images, a luminance-level modulation and gradient modulation (LM & GM) method is proposed by Zhao et al. ^[32]. However, this technique suffers from noise problems. Tang et al. ^[33] developed a method with a pre-processing step utilizing stacked generative adversarial networks (SGAN) to overcome contrast and high noise issues. After overcoming the noise issue, the resultant image is produced by stacked generative adversarial networks (SGAN) having enhanced boundaries and better contrast. Chen et al. ^[34] proposed a network for improving low-dose CT image quality using various down-sampling techniques. A technique by Naseem et al. ^[35], Cross Modality Guided Enhancement (CMGE), aims to solve the visualization problem of CT images.

Image enhancement methods primarily work using pixel-by-pixel operations, which causes slower processing of larger images. Ameen et al. ^[36] proposed a spatial domain image size-dependent normalization method to produce the enhanced images efficiently. Schaap et al. ^[37] proposed an efficient novel anisotropic diffusion method that increases the anisotropic method's smoothness for better visualization to mitigate the problem of over smoothness. Gandhamal et al. ^[37] suggested a local and global gray-level S-curve transformation technique for better enhancement of medical images.

CT images are mostly in DICOM format that encodes the image pixels values as 16-bit integers. Many image enhancement techniques proposed for 8-bit images can be utilized for 16-bit CT images too; however, the results can face over-illumination and extended contrast enhancement ^[37]. Another issue related to 16-bits CT images is that they cannot be shown directly on most of the existing displays, which support 8-bit images only. Therefore, for visualization purposes, these images need to be transformed into 8-bit images. The compression of the dynamic range of images to make them displayable on standard displays is termed tone mapping ^[38]. Many tone-mapping techniques developed for synthetic and natural high dynamic range (HDR) images can also be used for CT image visualization and enhancement. However, tone-mapping techniques are generally designed to preserve the naturalness and realism of the scene, but for medical images, the target is to enhance the contrast to distinguish the abnormalities visibly. Therefore, general tone-mapping techniques may not produce the desired enhancement of medical images.

This paper proposes a new clustering-based algorithm that is very effective for contrast enhancement of CT images, thereby leading to better visualization. The algorithm is computationally very efficient and works much faster than the existing contrast enhancement techniques. Moreover, the proposed algorithm converts 16-bit images to 8-bit for visualization on standard displays.

K-means clustering is a well-known technique for optimal partitioning data into a given number of bins. However, the k-means algorithm requires an initial guess that heavily influences the results as an algorithm with a strong tendency to converge to locally optimal solutions. K-means and its several invariants are slow in convergence when data size and the number of the clusters are large, as is the case with CT images, which have millions of pixels and need to be clustered into 256 bins (assuming a standard 8-bit format). We propose a novel divisive clustering algorithm that can efficiently handle this large dimension problem and has several advantages over k-means and its variants. For example, it does not require initialization, and therefore the results are consistent across multiple runs, unlike k-means, whose results depend on the initial solution. Experimental results show the proposed algorithm achieves nearly optimal performance and is highly effective in enhancing the contrast of CT images. It enhances both low and high-intensity regions and preserves the edges. We present extensive studies to validate the performance over a large number of test images.

The rest of the paper is organized as follows. Section 2 gives a detailed description of the proposed clustering method. In section 3, we apply it to 16-bit CT images to tone-map them to 8-bit images. We compare the results with the outputs of several existing algorithms in terms of visual quality and execution time. A subjective study conducted to evaluate the enhanced images' visual quality is also reported in this section. Finally, the paper is concluded in Section 4.

2.   Methodology

This section describes a detailed discussion of the proposed divisive clustering technique. Starting with entire image pixels in a single cluster, we keep adding a new cluster iteratively by splitting one of the existing clusters into two until 256 clusters are obtained. The clusters are then labeled in [0,255] range to transform the 16-bit medical images to contrast-enhanced standard 8-bit images. The flow of the complete algorithm is illustrated below in Figure 1, and a detailed explanation of the steps is given below.

2.1.   Splitting algorithm

To explain the proposed clustering technique, we first describe a novel technique of splitting a cluster into two. Let us assume that a set of data points ${S}_{k}, \mathrm{ }1\le k\le N$ is contained in a single cluster, and the data samples are arranged in ascending order. If every point in this cluster is replaced by the mean value of data, $\stackrel{-}{S},$ then the collective error induced in the samples in the cluster defined as the sum of the squared differences can be written as:

Now let us divide the cluster into two and again represent their respective samples by their mean values. This would increase the precision as two candidate values could be used instead of one to represent the data points. A good strategy would be to pick the split point that would reduce the overall combined error of two sub-clusters. With $N$ data points in the parent cluster, there are $N-1$ possible ways to split it. Let's assume that the parent cluster is split into two child clusters so that the first bin has first $M$ elements and the second bin has remaining $N-M$ elements. Characterizing ${\overline{S}}_{1}$ and ${\overline{S}}_{2}$ as mean values of both sub-clusters, the total error of these sub-clusters can be calculated as followed:

The split point $M$ is the boundary between the newly created sub-clusters. Its selection is important as it determines the errors given by (4). The optimal point with minimum combined error e can be picked after calculating the combined error of two sub-clusters for each of the N – 1 candidate split points. However, this brute-force technique is not requisite. As the data is sorted, by moving $M$ from left to right, the error in the left sub-cluster ${e}_{1}$ will increase while there will be a decrease in ${e}_{2}$, the right sub-cluster's error. Therefore, the combined error e given by (4) is a parabolic function with a single minimum point, as illustrated below in Figure 2. This point can be obtained by applying a binary search method having complexity O (log₂ N).

Assume that the selected cluster to be split is defined by the interval $\left[a, b\right]$. We calculate the slope of the error function e in the middle of this interval. The slope is computed by calculating the error at two adjacent points i.e., midpoint $m = (a+b)/2$ and $m+1$ and taking their difference. Denoting the error at these two points as ${e}_{m}$ and ${e}_{m+1}$, the slope can be written as:

A positive slope means that the minimum is present in the left half of the interval, and we can eliminate the right half from future searches. Likewise, a negative slope shows that we can exclude the left half of the cluster from future searches. Therefore, if the slope is negative, then left boundary a is moved to the midpoint m else b is moved to the midpoint. After a few iterations, i.e., ${\mathrm{log}}_{2}N$ iterations, the interval size becomes one so no further reduction would be possible. We end this process when this condition is reached, i.e., the interval size is reduced to one, or a point with the slope zero is obtained. In terms of minimizing the mean square error, the obtained point using this method would be the optimal split-point $M$ of the cluster.

The process of dividing a cluster into two can be repeated iteratively until the desired number of clusters is obtained. In each iteration, the cluster that has the maximum error can be selected for splitting, and this strategy would result in a maximum reduction in overall error. After each iteration, the size of the cluster to be split reduces, and the splitting procedure becomes even faster.

2.2.   Divisive clustering of CT images

Using the algorithm of splitting a cluster of data samples into two, we segment the input CT image data into 256 clusters. The complete algorithm for doing this is outlined below:

1.Convert two-dimensional image data to the one-dimensional array and sort in ascending order. Place all the data in a single cluster.

2.Initialize a few additional arrays – an array of cluster edges initialized with the first and last point of the sorted data and sorted in ascending order; an array to store mean values of the clusters; the third array to store the error in each cluster initialized with the value calculated using (1).

3.Split the cluster into two using the procedure explained above. Add the split point M to the array of edges and update the array of cluster errors accordingly.

4.Find the cluster that has the maximum error among the existing clusters.

5.Repeat steps 3 and 4 until 256 clusters are formed.

After clustering, the next step, as illustrated in Figure 1, is to label the clustered data. Every pixel contained in a cluster is allocated an integer value equal to the index number of the cluster. This assigns 8-bit integer values to each pixel. This labeled data is brought to the original order before sorting. This can be done efficiently by storing the original index numbers of the data before sorting. Finally, the data is converted to a two-dimensional image structure for visualization.

Note that our algorithm selects the split points such that the overall error is reduced. This results in the grouping of pixels with similar values in each sub-cluster, which in turn ensures high inter-cluster similarity and high intra-cluster dissimilarity. By assigning a different label to each cluster, we ensure that a good contrast is maintained between dissimilar pixels and regions.

As mentioned above, we construct an array of mean values of clusters as well. This 256x1 array of 16-bit values is stored with the labeled data and can be used to recover the original dynamic range of the data if needed. To do this, we can simply replace the labels with the corresponding cluster's mean value.

3.   Experimental evaluations

3.1.   Materials

For evaluation of the proposed algorithm, we compare the effectiveness (in terms of accuracy and visual quality) and efficiency of our algorithm with some existing methods. We conducted the visual assessment and subjective evaluation using five different 16-bit CT test images taken from publicly available datasets – head ^[39], chest ^[39], lungs1 ^[40], arms, and abdominal (https://wiki.idoimaging.com), and lungs2 ^[40].

For comparison, we picked two algorithms particularly designed for CT images enhancement Ameen et al. ^[36] and Gandhamal et al. ^[37] and six well-known tone-mapping operators: Kim and Kautz ^[41], Khan et al. ^[42], Reinhard and Delvin ^[43], Shibata et al. ^[44], Li et al. ^[38] and Liang et al. ^[45]. For ^{[36,37,38,42,44,45]}, we used the codes provided by authors themselves. While the codes used for remaining tone mapping operators are publicly available in the High Dynamic Range (HDR) Toolbox provided by Franceso Banterle (https://github.com/banterle/HDR_Toolbox). For all of them, we used default parameters defined by the authors. All the methods, including our algorithm, are executed using MATLAB. A Dell computer with 4 GB RAM and Intel Core i-3 CPU at 1.90 GHz processing speed is utilized for simulation.

3.2.   Visual assessment

The enhanced output images of the head CT image using the above-mentioned methods and our proposed algorithm are illustrated below in Figure 3.

In Reinhard and Delvin ^[43], edges information is preserved to some extent but suffers from loss of local contrast and details. Kim and Kautz ^[41] keep the detail of edges, though it suffers from low contrast overall. Khan et al.'s technique ^[42] represents better visibility of both higher and lower intensity areas, though it over-enhanced the edges information. Li et al. ^[38] and Liang et al. ^[45] preserve the edges information; however, they suffer from information loss due to over enhancement. Shibata et al. ^[44] represent information lost in the regions where abrupt intensity change occurred, whereas the areas with no intensity changes maintained their information. Al-Ameen et al. ^[36] and Gandhamal et al. ^[37] methods were purposely designed for CT images. Al-Ameen's method ^[36] provides better brightness level and contrast, but it also results in loss of local texture details in lower intensity regions. The information in higher and lower intensity range areas is lost using Gandhamal et al. ^[37] method. Compared to all other algorithms, our method provides a better balance of contrast among different regions and achieves a visually pleasing result.

The results of the other images are shown in Figures 4-8, and they all show the superiority of our algorithm over others. In all the results, our method well preserves the edges and sufficiently enhances the information of both higher and lower intensity areas, thereby leading to a clear view of image details. This is because our algorithm splits the highest error cluster at the point of minimum error to keep the similar pixels in one group. This results in assigning the same label to similar pixels. A pixel having different intensity values is likely to be grouped in different clusters with different labels. This whole process leads to good contrast preservation in the whole image. The proposed method has a consistent performance across all images, whereas other methods that perform well on some images have unsatisfactory performance others. To evaluate the overall performance for all images, we conducted a subjective study and present its results next.

3.3.   Subjective evaluation

The enhanced CT images of the mentioned datasets, as shown in Figures 3-8, were ranked by a number of participants in a subjective evaluation. The random arrangement of images for each dataset was shown to each participant on a 32-inch screen. A group of 30 participants was randomly selected from our software engineering department, most of them having basic to advanced knowledge of the image processing domain. All the nine enhanced 8-bit CT scan results of each dataset were simultaneously shown to each participant, who was requested to give a score from the range of 1 to 10, where 1 represents poor visual quality and 10 means excellent visual quality. The participants were allowed to observe the images as long as they wanted.

Table 1 shows the average rating scores given by participants to each of the images. The methods listed in descending order based on average scores given by participants are – the proposed method, Khan et al. ^[42], Al-Ameen et al. ^[36], Shibata et al. ^[44], Kim and Kautz ^[41], Reinhard and Delvin ^[43], Liang et al. ^[45], Li et al. ^[38] and Gandhamal et al. ^[37]. Hence, this user study validates that our proposed algorithm has the best visual quality as compared to all the other tested approaches. It also shows that the tone-mapping operators designed for natural and synthetic HDR images do not perform well on CT images. The method specially designed for CT images Al-Ameen et al. ^[36] was ranked at the third position after the proposed method and Khan et al. ^[42]. The winning score for each image is shown in bold font in the table, and the runner-up is shown in italic bold font. It can be seen that the proposed method is the winner in all cases, whereas no other method is constantly at the second position. This shows consistency of the performance of our method across a variety of images.

3.4.   Execution time

We also compared the computational efficiency of our method with the rest. The average time taken by each method to enhance the test images is mentioned below in Figure 9. The results showed that our algorithm has the shortest execution time among all the enhancement methods. We note that, compared with all the other methods, our proposed technique provides the best performance in terms of both execution time and image quality, thereby making the proposed algorithm a feasible medical image enhancement technique.

4.   Conclusion

A divisive clustering-based technique has been proposed for better enhancement and visualization of computed tomography images. The technique works efficiently on very large data, unlike the existing clustering algorithms that fail to converge in the same cases. Therefore, the proposed clustering algorithm is suitable for images which usually come with a very large number of pixels. Experimental results were presented using several test images, which showed that the proposed approach performed better than the existing methods in terms of visual quality and speed of execution. A user study was conducted in which the participants assigned the highest scores to the results produced by our method, further validating the performance.

Acknowledgments

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number MoE-IF-G-20-11.

Conflict of interest

The authors have no conflict of interest.