Magnetic resonance cholangiopancreatography image enhancement for automatic disease detection

Abstract

AIM: To sufficiently improve magnetic resonance cholangiopancreatography (MRCP) quality to enable reliable computer-aided diagnosis (CAD).

METHODS: A set of image enhancement strategies that included filters (i.e. Gaussian, median, Wiener and Perona-Malik), wavelets (i.e. contourlet, ridgelet and a non-orthogonal noise compensation implementation), graph-cut approaches using lazy-snapping and Phase Unwrapping MAxflow, and binary thresholding using a fixed threshold and dynamic thresholding via histogram analysis were implemented to overcome the adverse characteristics of MRCP images such as acquisition noise, artifacts, partial volume effect and large inter- and intra-patient image intensity variations, all of which pose problems in application development. Subjective evaluation of several popular pre-processing techniques was undertaken to improve the quality of the 2D MRCP images and enhance the detection of the significant biliary structures within them, with the purpose of biliary disease detection.

RESULTS: The results varied as expected since each algorithm capitalized on different characteristics of the images. For denoising, the Perona-Malik and contourlet approaches were found to be the most suitable. In terms of extraction of the significant biliary structures and removal of background, the thresholding approaches performed well. The interactive scheme performed the best, especially by using the strengths of the graph-cut algorithm enhanced by user-friendly lazy-snapping for foreground and background marker selection.

CONCLUSION: Tests show promising results for some techniques, but not others, as viable image enhancement modules for automatic CAD systems for biliary and liver diseases.

Key Words: Bile ducts; Liver diseases; Image enhancement; Structure detection; Magnetic resonance cholangiopancreatography

INTRODUCTION

Image enhancement is the first step in most image processing applications. Collectively, it may consist of several smaller tasks to serve one or more of the following purposes: (1) noise reduction; (2) background suppression; (3) size, intensity, color, brightness and contrast normalization; (4) region-of-interest (ROI) enhancement; (5) artifact removal; (6) quality and resolution improvement; and (7) pre-segmentation. The effectiveness of the image enhancement step greatly influences the subsequent steps to be undertaken. Ideally, this pre-processing step would enable elimination of all image inconsistencies and result in only the ROI being presented for evaluation, thus empowering concentration on the task at hand (e.g. labelling and diagnosis) without any hindrances of unwanted characteristics in the image.

In modern medical diagnosis of diseases affecting the bile ducts, a sequence in magnetic resonance imaging (MRI), called magnetic resonance cholangiopancreatography (MRCP), is used. Like MRI, MRCP is non-ionizing (no exposure to potentially harmful radiation), non-invasive (does not require any surgical procedures, and in many cases is even undertaken without contrast medium), flexible (enables internal organs to be visualized in many orientations and with different parameter settings) and usually requires no hospitalization and minimal preparation before examination. As powerful electromagnets (1.5 Tesla or higher) are used, this procedure is not applicable for patients with metallic and electronic implants.

This modality, although producing better quality images than other modalities (e.g. ultrasound), suffers from a number of shortcomings. These include: acquisition noise from the equipment, ambient noise from the environment, the presence of background tissue, other organs and anatomical influences such as body fat, and breathing motion. The amount of fluid in the body, diseases, deformities, artifacts, orientation, partial volume effect (loss of detailed structural information as a volume is represented in 2D) and acquisition settings can all significantly influence the image and make it very difficult for automated processing of the images. The human mind and visual system, along with specialised training, allows for medical practitioners to comprehend and compensate for such issues when examining the images. In certain cases, the human specialist uses these characteristics to orientate the image and better understand the anatomical situation presented in the image. The computer, however, especially when presented only with a single 2D MRCP image, is often defeated in such tasks, unless reliable pre-processing is undertaken to minimize unwanted influences.

Over the years, a multitude of algorithms, techniques and approaches have been used in image enhancement pre-processing. Varying in complexity, performance and even in the targeted data, the objective has been to produce enhancement of the images either as an intermediate or end result in a multi-stage processing scheme. Some of the more popular techniques for medical image processing are evaluated in this work, considering their applicability and performance on 2D MRCP images, with the objective of studying their suitability for the enhancement of MRCP images for use in preliminary detection of biliary diseases in medical computer-aided diagnosis (CAD) systems.

MATERIALS AND METHODS

Image enhancement algorithms

The scope of this work is limited to enhancing MRCP images such that noise is minimized and the clarity of the significant bile ducts is maximized. The term significant here is used to describe dilated bile ducts. Bile ducts produce the enzyme bile, which is used in the digestion of fat, absorption of oil-based vitamins and minerals, and in the removal of fat-soluble waste from the body. Diseases affecting the biliary structures, either from within the structure itself or from surrounding organs, cause blockage of the bile ducts. This, in turn, causes draining of bile from the ducts in the liver to the small intestines to be restricted. The build-up of bile causes the ducts to swell, and are thus dilated. In most cases, the main clue to the presence of diseases affecting the bile ducts in an MRCP image is the presence of the dilated ducts. As the targeted application of this work is for medical CAD systems for biliary diseases, the tests are performed accordingly. Note that the sensitivity of the algorithms may be tweaked by adjusting one or more coefficients, in order to visualize the less significant ducts as well, if so desired.

This section describes the popular techniques developed for the pre-processing of common images, which have then been successfully applied to medical images. Some customized techniques in the recent literature are also considered and described below. The results obtained from evaluating their suitability for MRCP pre-processing is given in the next section.

MRCP produces greyscale images, lacking color components, thus the information is represented by intensities. Technically, MRI machines measure the signal of the energy release that takes place when energized hydrogen protons that have been forced in a certain direction by a radio frequency (RF) pulse from the MRI machine return to their normal spin (precess) and low energy state once the RF pulse is turned off[1]. The mathematical information received is converted using a Fourier transform and are represented as a greyscale series of 2D slice images that are used in diagnosis.

Intensity thresholding

Thresholding is possibly the most primitive and least resource consuming of the pre-processing techniques. Used for binary image processing, it is known to be an effective and very popular technique in most applications, and is present in one form or another at various levels of complexity. For the purposes of this section, only simple thresholding is discussed. In a greyscale image such as a MRCP image, thresholding would be undertaken by specifying an intensity threshold, T, and removing all parts of the image, I, that fall below that intensity threshold. The process for a pixel at the coordinates (x,y) satisfies (1).

The problem is in determining the appropriate value for T. Experimentally, it was found that the intensity distribution in MRCP images differed among images. Normalization would be required before an automatic selection of T could be performed. Even then, it was found that standard normalization techniques, such as scaling, were ineffective with MRCP images. In this work, the threshold was set at intensity T = 150. This was determined experimentally and supported by prior experience with MRCP images to reasonably remove most of the background tissue. The threshold is applied on images after normalization by scaling them between 0-255, which is suitable for most computer monitors that support an 8 bits per pixel (bpp) format. Original MRCP images are stored in the Digital Imaging and Communications in Medicine file format and may be up to 16 bpp (usually 12 bpp).

More accurate results could be obtained through interactive selection of T on an image-by-image basis. This would be tedious but the interactive threshold scheme, if applied to an appropriate training set of MRCP images, could also be used to collect statistics that may aid in T approximation for future use. To provide a basis for comparison, an interactive threshold implementation using just a mouse click on the ROI is also implemented and evaluated in this work.

Common filters

Filtering is a popular method for removing undesired parts of images. The types of filters vary in the technique employed and the feature or characteristic affected. For instance, bandpass filters eliminate parts outside (above and below) the set bandwidth while preserving the parts within the bandwidth.

Gaussian: Smoothing is commonly undertaken using linear filters such as the Gaussian function (the kernel is based on the normal distribution curve), which tends to produce good results in reducing the influence of noise with respect to the image. The 1D and 2D Gaussian distributions with standard deviation for a data point (x) and pixel (x,y), are given by (2) and (3), respectively[2].

The kernel could be extended to further dimensions as well. For an image, the 2D Gaussian distribution is used to provide a point-spread; i.e. blurring over neighboring pixels. This is implemented on every pixel in the image using the convolution operation. The degree of blurring is controlled by the sigma or blurring coefficient, as well as the size of the kernel used (squares with an odd number of pixels; e.g. 3 × 3, 5 × 5 pixels, so that the pixel being acted upon is in the middle). The processing can be speeded up by implementing the filtering in the frequency rather than spatial domain, especially for the slower convolution operation (which is implemented as the faster multiplication operation in the frequency domain).

Wiener: The Wiener filter[3] may also be used for smoothing. This filter is the mean squares error-optimal stationary linear filter for images degraded by additive noise and blurring[4]. It is usually applied in the frequency domain (by taking the Fourier transform) and can be represented as (4).

where

H(u,v) = Fourier transform of the point spread function

P_s(u,v) = Power spectrum of the signal process, obtained by taking the Fourier transform of the signal autocorrelation

P_n(u,v) = Power spectrum of the noise process, obtained by taking the Fourier transform of the noise autocorrelation

It should be noted that there are some known limitations to Wiener filters. They are able to suppress frequency components that have been degraded by noise but do not reconstruct them. Wiener filters are also unable to undo blurring caused by bandlimiting of H(u,v), which is a phenomenon in real-world imaging systems[4].

Median: A popular non-linear noise-reduction filter is the median filter, which has been shown to be good at removing salt-and-pepper noise in images[5]. Sometimes known as a rank filter, this spatial filter suppresses isolated noise by replacing each pixel’s intensity by the median of the intensities of the pixels in its neighbourhood. It is widely used in denoising and image smoothing applications. Median filters exhibit edge-preserving characteristics (cf. linear methods such as average filtering tends to blur edges), which is very desirable for many image processing applications as edges contain important information for segmenting, labelling and preserving detail in images. This filter may be represented by (5).

where

w_F = w x w filter window with pixel (u,v) as its middle

Perona-Malik: Isotropic filters provide generalized diffusion to an image. When edges in the image are to be preserved, anisotrophic diffusion is required. The most popular filter for this is the Perona-Malik[6], which has widespread use. It is commonly believed that the Perona-Malik equation provides a potential algorithm for image segmentation, noise removing, edge detection, and image enhancement[7]. The basic idea behind the Perona-Malik algorithm is to evolve an original image under an edge-controlled diffusion operator. The equation may be represented as (6)[8].

where

g(s) = the decreasing smooth function, g(0) = 1; 0 < g(s) → 0 for s→∞

The original Perona-Malik equation, although numerically stable, was thought to be ill-posed by many. Subsequently, much research has gone into improving the algorithm, mainly by regularizing g. Also known as the Perona-Malik function, g is selected to be in (7)[8], making it “non-linear scale-space”.

Interactive graph-cut

Graph-cut has gained popularity in recent years as an effective object extraction technique. Essentially, often interactively, the user will select parts of the foreground and the background of the image. Graph-cut then uses a maximum flow/minimum cut algorithm to automatically identify the foreground objects in the image. There have been many variations of this technique used with varying degrees of success; e.g. Phase Unwrapping MAxflow (PUMA) via graph-cuts[9]. Combined with techniques such as Lazy-Snapping[10], the graph-cut method allows for good pre-processing through minimal interactive input by the user.

An example of the lazy-snapping graph-cut implementation used in this work is given in Figure 1, where Figure 1B shows the result of the interactive selection of the foreground and background (one line each) and the corresponding image boundaries identified by the algorithm. A silhouette of the resulting identified foreground object (biliary structures) is given in Figure 1C for clarity.

Figure 1 Lazy-snapping graph-cut identification of the significant biliary structures. A: Original magnetic resonance cholangiopancreatography thin slice image; B: Interactive identification of foreground and background, with the corresponding object boundaries identified by the algorithm; C: Silhouette of the identified objects.

In this work, three implementations of graph-cut were tested. The PUMA graph-cut implementation was used for denoising. The second implementation takes the PUMA result and thresholds it at intensity 150 for object extraction. The third was object extraction using the interactive lazy-snapping graph-cut implementation. The results obtained will be discussed later in this paper.

Histogram analysis

When a fully automatic system is desired, all stages of processing must be conducted without requiring user intervention. For this, understanding of the image intensity frequency distribution may aid in better understanding of the significant parts of the image. Recent literature on signal processing of MRCP images has identified observable patterns in the intensity histogram of those images. As described in[11], the peaks in the histogram tend to correspond to different characteristics of the anatomy and could be used to eliminate most of the noise, which would aid in enhancing the parts related to the bile ducts. Capitalizing on the thresholding scheme, dynamic thresholding via automatic selection of T could be performed in this way.

From analysis of the MRCP test images obtained in this work, it is noted that different types of MRCP images exist, with differences in their histogram distributions. The common types of MRCP images used in this work include thin slices, thick slabs and projected images. Thin slices are the essential 2D images showing the slices of the abdominal areas, each approximately 4-8 mm thick. A thick slab is usually much thicker, at about 50 mm, and has better signal strength (in terms of signal to noise ratio, SNR), thus providing a clearer image of the biliary structures (albeit the partial volume effect). A projected image is generally constructed using thick slab images at different angles (orientations) to approximate a rotational view of the area enclosed between the slabs. The projected images have a black border around them. The first peak in the histogram often represents the air in the body of a typical MRCP image, but this would be the second peak in projected images as the first peak would be the black border.

From the experimental testing conducted in this work, it was found that, in most cases, removing the entire first (second, in the case of projected images) peak until the minimum trough before the histogram rises again, results in extraction of the significant bile ducts. An example of this is given in Figure 2. The histogram in Figure 2B shows several peaks. The peaks corresponding to the border, air and soft tissues are removed at the point indicated on the slider. The result is then rescaled to fill the 0-255 intensity range, shown in Figure 2C.

Figure 2 Dynamic thresholding through histogram analysis to identify significant biliary structures. A: Original projection-type magnetic resonance cholangiopancreatography image; B: Intensity histogram of the image; C: Resulting thresholded image, rescaled to a 0-255 intensity level.

Multiresolution analysis

Impulse noise or abnormal intensity pixels in an image can severely affect the accuracy of many algorithms. One way of overcoming this phenomenon is through analysis at various resolutions, often by performing generalization akin to lossy compression to eliminate less significant information in the image, to emphasize the more significant regions. This allows for recognition of important aspects of the image, unhindered by the details (a lot of which includes background tissue and artifacts).

Wavelets: The most popular multiresolution tool is wavelets. There is much literature on the use of wavelets in noise reduction, such as wavelet thresholding in[12] and[13], and for medical images such as ultrasound[14] and recently on MRI[15]. There are many types of wavelets to choose from, each created for their specific characteristics. Examples include Haar, Daubechies, Symlets, Coiflets, Biorthogonal, Meyer, Gaussian, Mexican hat, Morlet, Complex Gaussian, Shannon, and Frequency B-Spline. A summary of some of the wavelet families and associated properties can be found in[16].

A non-orthogonal wavelet-based denoising implementation by[17] is used in this work. This technique, unlike the conventional wavelet denoising implementations, is reported to preserve the phase information, which is vital to human perception.

Directional wavelets: Conventional wavelets are restricted to horizontal and vertical directional limitations. Directional wavelets have been created to overcome this restriction. Among the most popular of these include curvelets[18], ridgelets[18], contourlets[19] and the directional wavelet transform[20]. The contourlet is also known as the pyramidal directional filter bank. They have been applied successfully to many types of image processing problems including edge detection[21], image enhancement[22], compression[23], texture retrieval[24] and adaptive denoising[25]. The contourlet and ridgelet toolboxes used in this paper were from[26]. A popular curvelet library is available at[27].

Scale-space analysis

Related to the multiresolution approach, scale-space analysis examines an image at varying levels of scales in order to attach significance to the different parts and objects in the image. The changing of scale is affected by using blurring filters such as the Gaussian kernel. However, isotropic filters, such as the Gaussian derivatives, tend to also blur the edges of structures, which in the case of MRCP images, may destroy the structural information. The solution to this problem is to replace it with anisotropic filters, such as the Perona-Malik[6] or Euclidean Shortening Flow[28] methods, which have better edge-preserving properties.

RESULTS

The purpose of this work was to examine the performance of popular pre-processing techniques applied for noise reduction and image enhancement of the ROI; i.e. the significant dilated biliary structures of the MRCP images. Certain assumptions are used to form the basis of this work. First, the test data is robust, with data from different patients with different builds and disease characteristics, taken from different sources by various radiographers using MRI machines of different makes and models with different orientations and parameter settings. The data was sourced from the internet to enable availability and the mentioned robust characteristics. Through collaboration with medical experts, the types of test images used were confirmed to be typical examples of clinical MRCP images used for diagnosis. Noise is assumed to be unknown, thus equal treatment is given to all images.

As the algorithms to be tested differed in their characteristics, they were grouped into several test sets. The milder algorithms for denoising of the images were tested first. This was split into two sub-tests; first for the MRCP images containing the presence of a significant amount of background tissue, and the second set for images containing numerous objects (multiple parts of the bile ducts and other organs). The images in these tests were then reused for the next test using the algorithms implementing more severe filtering, with the target of highlighting only the significant biliary structures in the images.

Denoising results

The experiment for denoising was conducted using the Gaussian, Wiener, median and Perona-Malik filters, PUMA graph-cut and the directional wavelets (contourlet and ridgelet). Some of the results obtained for the first sub-test for MRCP images, with heavy influence of noise and background tissue, are shown in Figure 3. It was observed that the performance of the different algorithms varied significantly. The Gaussian filter with a 5 × 5 kernel and σ = 5 produced considerable blurring of the images, reducing the background tissue detail [obvious for images (a) and (b)]. The resulting lesser influence of the background, relative to the biliary structures, eases the background noise and tissue removal process. However, it was noticed that when the background tissue intensity was high, as in (d), the average intensity difference between the biliary structures and the background was reduced. This complicates subsequent processing. The isotropic filter also smoothed the edges of the biliary structures, which is undesirable as it makes accurate structure identification difficult.

Figure 3 Denoising results for images with heavy tissue background.

The Wiener filter also blurred the image, and caused higher smoothened background intensity [see the bottom left of (a) and mid-left of (b)]. A similar effect is noticed with the median filter as well. As expected, the anisotropic Perona-Malik filter achieved blurring of the background that is comparable to the Gaussian filter, but preserved the edges better [see (a) and (c)]. The PUMA graph-cut enhanced the images with phase information, increasing the detail. However, this was counter-productive as the background and noise was also more prominent. The contourlet achieved slight noise reduction and diffusion of detail, but the ridgelet caused emphasis of the background. Overall, the anisotropic Perona-Malik filter was found to be the most suitable for noise reduction in 2D MRCP images with heavy background tissue present.

Results of the sub-test for denoising MRCP images with multiple objects are given in Figure 4. Generally, all the blurring filters performed better on these images as there was less background tissue present. The blurring improved the composition of the objects by compensating for inter-pixel inconsistencies, such as holes and fluctuations in intensities due to impulse noise, improving object detection. Gaussian analysis performed well in reducing the influence of the tissue when less background was present, while PUMA and ridgelet analyses enhanced the (unwanted) detail [see (b) and (c)].

Figure 4 Denoising results for images with multiple objects.

Object extraction results

The main objective of this work was to identify appropriate image enhancement for the detection and extraction of the significant biliary structures, by reducing or eliminating the background. The more aggressive of the discussed algorithms were tested for this purpose, using the same sets of images.

Figure 5 gives the results for the images with heavy background influence. Binary hard-thresholding at intensity 150 removed the background tissue in all of the images. Unfortunately, some parts of the biliary ducts were also eliminated, especially in (b) where large parts of the bile ducts were represented by lower intensities. Using the PUMA graph-cut enhancement prior to thresholding did not improve the results. The result was the same as those of Threshold 150 and so are not shown separately. Interactive thresholding produced satisfactorily accurate results. The histogram-based dynamic threshold scheme also performed well. The non-orthogonal noise compensation wavelet implementation produced good results, but left behind a lot of tissue in images with very heavy backgrounds [e.g. (d)]. The interactive lazy-snapping graph-cut was the most accurate, even with only a marker each for the foreground and background [except for (d), which required two background markers].

Figure 5 Structure extraction for images with heavy tissue background.

In the case of less background tissue but more objects, the results in Figure 6 show that thresholding at intensity 150 removed all the background tissue and significant amounts of unwanted objects. Parts of the bile ducts lost in this scheme may be recovered through post-processing with a scheme such as region-growing seeded with the detected bile ducts. Interactive threshold selection managed to reduce the amount of lost biliary branches, as seen in (b). The interactive lazy-snapping graph-cut required three markers for the background in order to obtain the results, although only one marker had to be used for the foreground. Overall, the significant bile ducts were successfully detected and extracted, albeit with the lower intensity parts also being removed. Unfortunately, none of the algorithms were able to remove background objects that were of similar high intensity as the significant biliary structures. Post-processing utilizing higher-level algorithms, such as shape or texture analysis, would be required for proper labelling and removal of these objects.

Figure 6 Structure extraction for images with multiple objects.

Scale-space is implemented as an isotropic approach using Gaussian filtering, or anisotropically using Perona-Malik or ESF, at different blurring levels. The results are not shown here as they are essentially what have been shown in the individual Gaussian and Perona-Malik implementations above. The strength of scale-space will be more apparent when information derived in each scale (level) of blurring is mapped from the largest to the lowest, allowing for hierarchical information among the objects or parts thereof to be identified. An example implementation for bile duct hierarchical structure detection using MRCP images can be found in[29].

DISCUSSION

This paper evaluated several well-known and proven image enhancement approaches, specifically on a robust test set of MRCP images used for the diagnosis of diseases affecting the bile ducts. Tests were conducted for denoising, as well as object detection and extraction, in tests sets with heavy background noise and tissue influence as well as with multiple objects. The tested algorithms included filters (e.g. Gaussian, median, Wiener and Perona-Malik), wavelets (e.g. contourlet, ridgelet and a non-orthogonal noise compensation implementation), graph-cut approaches using lazy-snapping and PUMA, and binary thresholding using a fixed threshold and dynamic thresholding via histogram analysis.

The results varied as expected since each algorithm capitalized on different characteristics of the images. For denoising, the Perona-Malik and contourlet approaches were found to be the most suitable. In terms of extraction of the significant biliary structures and removal of background, the thresholding approaches performed well. The interactive scheme performed the best, especially by using the strengths of the graph-cut algorithm enhanced by the user-friendly lazy-snapping for foreground and background marker selection. Improving accuracy in labelling and extraction would require further post-processing, via higher-level strategies such as shape, texture and statistical analysis, or through self-learning and optimizing algorithms such as neural networks or support vector machines. The selected algorithms above could be considered as suitable for preliminary disease detection affecting bile ducts in CAD systems that may be developed in the future.

COMMENTS

Background

The noisy, dynamic intensity range and artifact susceptible nature of magnetic resonance cholangiopancreatography (MRCP) images makes it very difficult to develop computer-aided diagnosis systems that can aid in the preliminary screening of the images, as the developed algorithms are often defeated by the inter-image variances. Automated screening is an increasing need with the ever increasing amounts of digital image data generated by the various sequences of MRCP examinations.

Research frontiers

Identification of appropriate strategies to overcome the weaknesses of the images is necessary for the development of reliable automated systems for MRCP. Specifically, for the purposes of screening, techniques that enhance the biliary structures while minimizing the image noise and unimportant structures.

Innovations and breakthroughs

The field of image processing has introduced a large number of image processing schemes for the purposes of pre-processing images. These differ in terms of strategies and abilities. However, an optimum technique has yet to be identified among the various schemes due to lack of performance comparison on the effects of those schemes on MRCP images.

Applications

The characteristics of the various algorithms are studied and applied experimentally to the images. The optimum parameter settings are determined to achieve the best results of the respective algorithm on the MRCP image test database. The outcomes are analyzed.

Terminology

Denoising is the process of removal of noise present in the signal and can be in the form of environmental ambient radiation, equipment noise due to the movements of the magnets and electromagnetic interference in the magnetic resonance imaging equipment, motion noise due to minor movements of the patient or internal body fluids, and various artifacts. Image enhancement in this work relates to improving the quality of the image by highlighting areas/structures of interest while suppressing area/structures that make identification difficult.

Peer review

The paper is well-written paper but lacks applicability in every-day radiology practice. However, the comments regarding the improvement of image quality is note-worthy.