Automatic Exposure Correction And Local Contrast Setting For Diagnostic Viewing of Medical X-ray Images

Institutionen för medicinsk teknik

Department of Biomedical Engineering

Examensarbete

Automatic Exposure Correction And Local

Contrast Setting For Diagnostic Viewing of Medical

X-ray Images

Examensarbete utfört i medicinsk teknik vid Tekniska högskolan i Linköping

av

Ottar Pehrson Skidén

LiTH-IMT/MI30-A-EX--10/491--SE

Linköping 2010

Department of Biomedical Engineering Linköpings tekniska högskola Linköpings universitet Linköpings universitet SE-581 83 Linköping, Sweden 581 83 Linköping

Automatic Exposure Correction And Local

Contrast Setting For Diagnostic Viewing of Medical

X-ray Images

Examensarbete utfört i medicinsk teknik

vid Tekniska högskolan i Linköping

av

Ottar Pehrson Skidén

LiTH-IMT/MI30-A-EX--10/491--SE

Handledare: Mats Andersson

imt, Linköpings universitet

Hagen Spies

Sapheneia Commercial Products AB

Examinator: Hans Knutsson

imt, Linköpings universitet

Avdelning, Institution Division, Department

imt

Department of Biomedical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2010-05-10 Språk Language  Svenska/Swedish  Engelska/English   Rapporttyp Report category  Licentiatavhandling  Examensarbete  C-uppsats  D-uppsats  Övrig rapport  

URL för elektronisk version

http://www.imt.liu.se

ISBN — ISRN

LiTH-IMT/MI30-A-EX--10/491--SE Serietitel och serienummer

Title of series, numbering

ISSN —

Titel Title

Automatisk Exponeringskorrektion Och Lokal Kontrastinställning För Visning av Digitala Röntgenbilder

Automatic Exposure Correction And Local Contrast Setting For Diagnostic View-ing of Medical X-ray Images

Författare Author

Ottar Pehrson Skidén

Sammanfattning Abstract

To properly display digital X-ray images for visual diagnosis, a proper display range needs to be identified. This can be difficult when the image contains colli-mators or large background areas which can dominate the histograms. Also, when there are both underexposed and overexposed areas in the image it is difficult to display these properly at the same time. The purpose of this thesis is to find a way to solve these problems. A few different approaches are evaluated to find their strengths and weaknesses.

Based on Local Histogram Equalization, a new method is developed to put various constraints on the mapping. These include alternative ways to perform the histogram calculations and how to define the local histograms. The new method also includes collimator detection and background suppression to keep irrelevant parts of the image out of the calculations. Results show that the new method enables proper display of both underexposed and overexposed areas in the image simultaneously while maintaining the natural look of the image. More testing is required to find appropriate parameters for various image types.

Nyckelord

Keywords image processing, weight, histogram equalization, digital X-ray image, contrast setting, local

Abstract

To properly display digital X-ray images for visual diagnosis, a proper display range needs to be identified. This can be difficult when the image contains colli-mators or large background areas which can dominate the histograms. Also, when there are both underexposed and overexposed areas in the image it is difficult to display these properly at the same time. The purpose of this thesis is to find a way to solve these problems. A few different approaches are evaluated to find their strengths and weaknesses.

Based on Local Histogram Equalization, a new method is developed to put various constraints on the mapping. These include alternative ways to perform the histogram calculations and how to define the local histograms. The new method also includes collimator detection and background suppression to keep irrelevant parts of the image out of the calculations. Results show that the new method enables proper display of both underexposed and overexposed areas in the image simultaneously while maintaining the natural look of the image. More testing is required to find appropriate parameters for various image types.

Acknowledgments

I would like to thank Hagen Spies for the opportunity to do this work at Sapheneia and for providing me with useful thoughts and ideas, 1:st Res. Eng. Mats An-dersson for being my supervisor and Prof. Hans Knutsson for being my examiner. Great thanks to everyone at Sapheneia for making me feel welcome during my time there. Also, thank you Jenny for your encouragement and support.

Contents

1 Introduction 1

1.1 Problem . . . 1

1.2 Objectives . . . 2

1.3 Method . . . 2

1.4 Digital X-ray Images . . . 2

2 Tools 5 2.1 Orientation Estimation . . . 5

2.2 The Hough Transform . . . 7

3 Contrast Control 9 3.1 Global . . . 9 3.1.1 Histogram Equalization . . . 9 3.1.2 Gray-Level Grouping . . . 10 3.2 Local . . . 11 3.2.1 Block-wise Control . . . 12 3.2.2 Sliding Window . . . 14

3.3 Constraining the mapping . . . 16

3.3.1 Selective Approach . . . 16

3.3.2 Constraining the local mapping . . . 17

4 A New Approach 21 4.1 The Gradient Magnitude Weight . . . 21

4.2 Weighted Histogram Calculation . . . 22

4.2.1 Weighted Global Mapping . . . 23

4.2.2 Weighted Local Mapping . . . 24

4.3 Collimated Images . . . 24

4.4 Selective Approach With Weight . . . 27

4.4.1 Finding Breakpoints . . . 27

4.4.2 Local Mapping With Selective Approach . . . 28

4.5 Local Histograms With Global Information . . . 28

4.5.1 Using the Weight To Decide Local/Global-Ratio . . . 29

4.5.2 Constraining the Upper Limit . . . 30

4.6 Selecting the Weight . . . 31

4.6.1 The Histogram Weight . . . 31

x Contents

4.6.2 The Local/Global-Ratio Weight . . . 32 4.6.3 The Upper Limit Weight . . . 34 4.7 Multiple Constrained Local Histogram Equalization . . . 34

5 Results 37

6 Discussion 51

6.1 Suggestions For Further Work . . . 51 6.2 Conclusion . . . 52

Chapter 1

Introduction

This work has been conducted at Sapheneia Commercial Products AB (SCPAB) in Linköping, Sweden. SCPAB specializes in diagnostic image enhancement of X-ray, CT, MR and Ultrasound images to enable dose-reduction and improved image quality.

1.1

Problem

Because of their huge dynamic range, raw digital X-ray images are not suitable for direct visual diagnosis by human observers. To prepare such images for display an identification of a relevant display method is necessary. The histogram of the image can be used to determine the intensity distribution and to set a suitable display range. This can be difficult due to the presence of collimated areas and background areas, which even though they are usually irrelevant can dominate the histograms.

Many X-ray images are acquired using collimators. Collimators are used to shield parts of the body that are not to be examined from being unnecessarily exposed to the X-ray radiation. This result in edges and parts of the acquired image that are irrelevant to diagnostic purposes, these areas will interfere with the desired histogram of the image. Figure 1.1 (b) shows a collimated image.

Also, the image can contain underexposed and overexposed areas simultane-ously which are difficult to present in good quality at the same time since giving one part of the grayscale more room tends to lead to suppression of the other parts. It can be argued that an examination is conducted to look at a specific part of the body, so the display range only needs to be set for that part of the image. However, it is always good to get as much information as possible, to be able to see other parts well. Figure 1.2 illustrates the difficulty of displaying different areas simultaneously. In figure 1.2 (b), the mapping was performed to see the screws and the lower part of the image well. This, however, makes the lungs and hands almost disappear. In figure 1.2 (c), the lungs are displayed properly but the lower part of the image is very low in contrast.

2 Introduction

1.2

Objectives

The goal of this thesis is not to perform image enhancement, it is finding a robust way to set the display range of the image prior to, or after, the actual enhancement. The objectives include:

• Finding a way to suppress irrelevant areas in the image.

• Enabling a proper display of under- and overexposed areas simultaneously. • Doing the above without disturbing the original intensity relationships too

much, i.e. maintaining the natural look of the image.

1.3

Method

Using a variety of images of different anatomies and from different manufacturers, a few methods will be explored, combined and further developed too meet the requirements of diagnostic image display.

1.4

Digital X-ray Images

The amount of X-rays that is absorbed by the body (or any object) differs de-pending on the kind of tissue (bone, skin, cartilage and so on), its density, as well as the thickness of the region being exposed. This is what makes it possible to use X-rays for medical imaging, to distinguish tissues from one another. But to find an adequate exposure limit can be a bit tricky if the exposed object has a large density variation, e.g. if the entire body is being examined it is difficult to expose the lungs and the hips simultaneously in a good way. Some areas will be underexposed (too few photons detected) and some will be overexposed (too many photons detected). This leads to a low SNR in these areas. Digital X-ray images enables post-processing of the acquired images so that enhancement algorithms can be applied. Figure 1.1 shows two of the images that were used in the evalua-tion. Notice the artifacts that form a border around figure 1.1 (b). Artifacts like these can cause misleading information, for example spikes in the histograms and irrelevant structures when performing gradient based operations.

There are two different kinds of digital X-ray systems, Computed Radiogra-phy (CR) and Direct RadiograRadiogra-phy (DR), they differ in the way they register the photons [6]. In CR, the X-rays are indirectly conversed to a digital image via a latent image. To do this a phosphor-based imaging plate is used instead of the traditional film. When the X-rays come in to this layer they excite electrons which are then trapped in the film. When the film is later scanned using a laser system, the electrons relax and emit visible light which is registered and digitized, a digital image is obtained.

DR is a newer technique and is a way to to obtain a digital image without having to scan a film after exposure. There are various ways of doing this but one way is to use large flat panel detectors with thin-film transistor arrays. The array

1.4 Digital X-ray Images 3

is covered with a scintillator material, e.g. thallium doped cesium iodide which fluoresces when the X-rays strike. This light is then registered and turned into an electrical signal by photodiodes.

(a) (b)

4 Introduction

(a)

(b) (c)

Figure 1.2: Original image (a), image mapped to display the lower part of the image well (b) and image mapped to display the lung region well (c).

Chapter 2

Tools

This chapter provides background theory on the methods that have been used in the thesis.

2.1

Orientation Estimation

Everything interesting in an image has structure in the form of lines and edges, meaning that there will be more or less sudden intensity differences in these areas. Finding these is therefore of interest in most image processing tasks. Sudden inten-sity differences means that a point on an edge will have relatively large derivative in one direction or the other. Filtering the image with a derivative kernel will pick up these derivatives and thereby detect the edge. To find the edges in all directions, two kernels with perpendicular directions are applied and then combined:

G =qG2

x+ G2y (2.1)

This will result in a gradient magnitude image in which the intensity encodes the edge strength, see figure 2.1. Examples of kernels that can be used are the derivative of a Gaussian or the Sobel operator. Here, quadrature filters were used [2].

Quadrature filters consist of an even real kernel (line detector), fe, and an

odd imaginary kernel (edge detector), fo, that have a relative phase shift of π/2.

The kernels are convolved with the image and combined to obtain the magnitude image, q, in the kernel direction:

f = fe+ ifo (2.2) qe= fe∗ image qo= fo∗ image (2.3) q =pq2 e+ q2o (2.4) 5

6 Tools

(a) (b)

Figure 2.1: Original image (a) and edge magnitude image (b, square root for a better view).

The argument of the quadrature filter response yields the phase in the filter di-rection, a dominant even filter response represents a positive or negative line with

θ = 0 and θ = π respectively.

θ = arg(qe+ iqo) (2.5)

The direction of an edge can lie between [0 2π]. If an edge is rotated by 180◦ it will change sign, if a line is rotated 180◦ it will remain the same. To avoid the orientation being 180◦different depending on if the sign is positive or negative, four quadrature filters evenly distributed with directions ϕk in the following manner:

ϕk = (k − 1) π 4 k = 1...4 (2.6) z = 4 X k=1 qk  cos(2ϕk) sin(2ϕk)  =  q1− q3 q2− q4  (2.7) z = z1+ iz2 (2.8) θ = arg(z) (2.9)

result in an angle of the orientation vector, z, that is twice the orientation of the structure, orthogonally oriented structures will have maximally different ori-entation vectors. This is an attractive feature if the neighbourhood only has one

2.2 The Hough Transform 7

dominating orientation, as opposite direction angles will now be mapped next to each other instead of opposite each other. However, if the neighbourhood is more complicated, e.g. if there is a T-junction where two orthogonal lines meet, the magnitude responses will be suppressed since they will cancel each other out. More on this can be found in [2].

2.2

The Hough Transform

Locating and extracting the strongest edges in an image can often be of interest in image processing. Each edge point in an image fulfills the equation:

y = kx + m (2.10)

for some k and m. If the equation is instead written as:

k = y/x − m/x (2.11) the image is said to be transformed into the parameter space. This transform is called the Hough transform and can be used to locate lines in an image. A point in the parameter space represents a line in the image space. Since the slope in equation 2.11 can become infinte, another defintion is suitable:

ρ = xcos(θ) + ysin(θ) (2.12) where ρ is the distance of the line from the coordinate center and θ is the ori-entation of the line. Using the local oriori-entation estimation from section 2.1 will reduce the computational burden since both the positions of the lines as well as their orientations are already known. This way, each point in the image will only yield one point in the parameter space [4]. Also, the contribution of each point in the image to the parameter space can be weighted using the edge magnitude. The points in the transformed image that have the highest intensity correspond to the strongest lines in the original image. Figure 2.2 shows a hough transform with the intermediate steps. In the orientation image, the hue of the colour represents the orientation and the brightness represents the certainty measure determined by the gradient magnitude.

8 Tools

(a) (b)

(c) (d)

(e)

Figure 2.2: Test image with noise (a), gradient magnitude image (b), orientation angle (c), resulting hough transform (d) and the detected lines superimposed on the test image (e).

Chapter 3

Contrast Control

This chapter introduces the methods that were used to create transformation func-tions, namely, conventional Histogram Equalization (HE) and Gray-Level Group-ing (GLG). These are later used to create a new approach, presented in the next chapter. They were also developed further to overcome some of the limitations they exhibit. Different variations of the methods, including global, local, selective, contrast-limited variations and mixtures of these were used.

3.1

Global

There are many different contrast control techniques to choose from depending on the application, these will of course differ greatly in complexity and output. Many of the techniques available are generalizations of conventional techniques, like the many different alterations of the well known Histogram Equalization. Here, we will consider Histogram Equalization and a newer method called Gray-Level Grouping.

3.1.1

Histogram Equalization

Histogram Equalization (HE) is a very simple method of adjusting the contrast of an image. Given an image with L intensity values r = [0 L − 1], the probability of a pixel having a certain intensity is

pdf (r) = H(r)/(M × N ) (3.1) which is called the probability density function of the image. H is the histogram of the image and M and N is the width and height of the image. Basically, the pdf is the histogram normalized to [0 1].

To obtain a mapping function that will adjust the contrast of the image, the cumulative sum of the pdf is used. This is the conventional histogram equalization, figure 3.1 (a). T (i) = (L − 1) i X j=0 pdf (j) i = 0..L − 1 (3.2) 9

10 Contrast Control

This procedure seeks to create a flat histogram, i.e. spreading out the intensity values evenly over the whole dynamic range of the image. The main limitation of HE is that it can not produce good results when there are large histogram compo-nents at some locations and small at others, the smaller compocompo-nents will always be suppressed. If, for example, an image has a relatively large dark background, the equalized image will have a washed out appearance due to the background being mapped over a large part of the grayscale. The equalization procedure can also cause over-contrast that will look unnatural. Figure 3.2 (a) shows a histogram equalized image.

(a)

Figure 3.1: HE of the histogram corresponding to the test image (a).

As mentioned earlier, there are many different variations of HE. Bi-histogram equalization is a common method that splits the histogram into two parts based on its mean and then performs HE on the two parts separately. Usually this yields a better result but since HE is still applied it has the same limitations, though to a smaller extent. The same limitations apply to other HE variations. Other ways of constraining the equalization procedure will be dealt with later in this chapter.

3.1.2

Gray-Level Grouping

Gray-Level Grouping (GLG) is a different way of obtaining a mapping function. It is unlike many of the histogram equalization variations in its algorithm and was developed to utilize the entire grayscale more efficiently. The basic idea is to group the histogram components into a number of bins that are then redistributed uniformly over the grayscale.

First, the nonzero histogram components are assigned to bins and the gray-level interval in each bin is recorded. The bins are then grouped based on their size,

3.2 Local 11

(a) (b)

Figure 3.2: Original image (a), HE image (b).

in each iteration the smallest component is grouped with its smallest neighbour and the gray-level intervals are updated. This process can be repeated until there are only two bins or a specified number of bins left. Each bin will now contain a gray-level interval from the input image.

The remaining bins are then distributed evenly over the entire grayscale and within each bin, the corresponding gray-level interval is linearly distributed. The number of bins that the original histogram is grouped into can be set (Fast-GLG) or every combination can be performed and compared using a contrast measure to obtain the best result (GLG). The latter however is more computationally intensive and it is difficult to find a suitable measure. Figure 3.3 (b) shows a Fast-GLG image. The full algorithm can be found in [1].

3.2

Local

The main limitations of the global methods include the fact that they will suppress the relatively small histogram components and that they are not able to solve the problem of displaying under- and overexposed areas in the image simultaneously. Using local information makes it possible to overcome these limitations, objects in low contrast areas can be discerned more readily compared to the global methods, see figure 3.4.

12 Contrast Control

(a) (b)

Figure 3.3: Original image (a), Fast-GLG image (b).

in a global sense, but in a local. This means that two pixels with the same intensity in the original image can have different intensities in the resulting image. Also, with only local information, the mapping may spread out the local intensities over too large a range if no constraints are applied. Care has to be taken so the resulting image does not have an unnatural look. This can be a bit difficult, as we will see in the next sections.

3.2.1

Block-wise Control

Dividing the image into non-overlapping sub-images, or blocks, of arbitrary size enables a mapping based on local histograms. Each block is processed separately and with a suitable method. To avoid too much variation between the mappings of the different blocks, the mapping of one block can be slightly weighted with the mappings of neighbouring blocks. Finally the blocks are fused together at the block borders using bilinear interpolation. If the blocks are simply stitched together, there will be strong edge effects at the borders. This method will not be truly local since only the pixels in the middle of each block will be mapped based on a symmetric neighbourhood, the others will use a mapping that is interpolated. However, the block method will be able to resolve low-contrast areas much better than the global methods, see figure 3.5.

3.2 Local 13

(a)

(b)

(c)

Figure 3.4: Original image (a), global HE image (b), local HE image (c). The large bright area suppress the small low contrast dark area in the global HE.

14 Contrast Control

(a) (b)

Figure 3.5: Image blocks HE-processed without interpolation (a) and with interpolation (b).

3.2.2

Sliding Window

To perform a truly local process, a window of suitable dimensions is defined around the starting pixel in the image and the method of choice, often local histogram equalization (LHE), is applied for only that pixel, using its local information. Then, the window is moved to next pixel and the process is repeated for all pixels in the image. This is a significantly more computationally intensive procedure than the others. GLG is not very well suited for sliding window operations as it is requires too much computation and is therefore very slow when applied to all the pixels in an image.

The size of the window determines how strong the local effect will be, a small window will be able to extract low contrast areas better than a large one but the noise amplification will be larger as well. If the window is too small there will be too few intensity values spread out leading too a noisy look. The mapping will be more and more like the global mapping as the size of the window is increased, figure 3.6.

3.2 Local 15

(a) (b)

(c) (d)

Figure 3.6: LHE with windows of different sizes, 101 × 101 (a), 201 × 201 (b), 401 × 401 (c), 801 × 801 (d) pixels respectively and an image size of 2795 × 1764

16 Contrast Control

3.3

Constraining the mapping

3.3.1

Selective Approach

When an image contains irrelevant parts, such a large background or is collimated, as mentioned earlier, these areas will tend to dominate the histograms and they will largely dominate the resulting mapping as well. Finding a way around this is not necessarily simple. However, if a suitable breakpoint can be identified these parts of the histograms can be suppressed by limiting irrelevant intensity values to be mapped into a smaller range than they would be if no restrictions were introduced. Each segment determined by the breakpoints is processed separately within their new limits, which are set to some suitable values. The problem is how to find the breakpoints. One way to do this will be presented in section 4.4.1.

The relevant intensities will be allowed to occupy a larger part of the dynamic range in the global mappings and thereby improve the result, often significantly. Figure 3.7 shows an image acquired using Fast-GLG without restrictions and an image acquired using selective Fast-GLG, note that the background is suppressed much better in the latter. If the intensities of the irrelevant parts of the image lie in the same parts of the histogram as the relevant ones however, this technique is of course useless and other approaches will have to be used. This is one of the things that the new approach presented in chapter 4 tries to solve by eliminating the influence of the irrelevant areas.

(a) (b)

3.3 Constraining the mapping 17

3.3.2

Constraining the local mapping

The above mentioned local methods, in their original form, obviously produce results that are not good enough for medical images. The result should maintain some form of natural look, e.g. the lung should not be much brighter than the hips in the resulting image. To prevent this, some way of constraining the mapping function must be introduced.

Different ways of constraining the local histogram equalization have been pro-posed. Pizer, [5], introduced a contrast limitation where the histogram is clipped at a certain level and re-normalized. The maximum slope of the resulting map-ping function will thereby be smaller compared to the mapmap-ping function without a clip-limit. This way, drastic changes between neighboring intensities in the map-ping function will be decreased, preventing over-contrast. Histogram clipmap-ping can also be applied to global operations. Instead of just clipping the histogram above a certain threshold, a soft function like arctan can be used to suppress the high peaks in the histogram, this way the original shape of the histogram is preserved, see equation 3.3 where a lower ψ suppresses the histogram peaks more. Figure 3.8 shows histograms modified with a clip limit (a) and with arctan (b).

H0 = atan(H/ψ) (3.3)

(a) (b)

Figure 3.8: Clipped histogram (a) and histogram modified with arctan (b).

Iyad and Hao, [3], used a calculation of the HE transformation function that introduces a second term which uses the mean of the image to constrain the map-ping and thereby avoiding noise enhancement and maintaining the mean image brightness. Zhu, [7], presented a method, part of which uses an alternative way of defining the local histogram. Instead of using only the histogram of the neigh-bourhood around a pixel, the global histogram is allowed to influence as well. The way this is done can be seen in equation 3.4 where 0 ≤ α ≤ 1, Hwindow is the

normalized histogram of the window and Hrestis the normalized histogram of the

rest of the image. With this definition, the resulting local/global histogram will be normalized no matter which α is used.

H = αHwindow+ (1 − α)Hrest (3.4)

This means that local histograms will have all the input gray levels in them but to a lesser extent so the truly local histogram is still emphasized, figure 3.9.

18 Contrast Control

When the mapping function is calculated this will keep the local values from being spread out over the entire available range. Figure 3.10 shows local histogram equalization images using arctan, local/global histograms and a combination of the two constraints. When comparing these images it can be seen that using only arctan, the mapping is still allowed to use the entire grayscale since the local histograms contain the same components as before. The local/global approach produces good results in some parts of the image but over-emphasizes the local histograms in others, the combination of the two results in an image with the same problem. This will be dealt with in the upcoming chapter.

(a) (b)

Figure 3.9: Local histogram of an area in the test image (a) and alternative local histogram (b).

3.3 Constraining the mapping 19

(a) (b)

(c) (d)

Figure 3.10: LHE with local histograms modified with arctan (a), LHE with global/local histograms (b), LHE with both constraints (c) and LHE without constraints (d).

Chapter 4

A New Approach

In this chapter, a new approach to contrast equalization is presented. This builds on the calculations and constraints on the methods presented in the preceding chapter.

4.1

The Gradient Magnitude Weight

To find a way to map an image in an acceptable way for diagnostic viewing, the main issue that needs to be solved is that irrelevant parts of the image must not interfere with the mappings. For example if the image is collimated or has a large background, these areas need to be suppressed in the map calculations. This is needed not only for the global methods but also for the local regions as will be explained in section 4.5.

For the local transformations, the following additional issues need to be ad-dressed:

• The local mappings occupy too much of the available grayscale, amplifying noise and causing over-contrast.

• The resulting image does not maintain a natural brightness balance over the entire image.

Since certain areas in an image can sometimes dominate the histogram, a way of focusing on the interesting parts of the image is desirable. This would decrease the influence of the irrelevant parts of the image when performing the map calculations. Areas such as background and other ’flat’ areas in the image are often relatively homogenous compared to the medically interesting areas. Therefore, we propose an approach that uses weights when performing for example histogram calculations. Such weight images can also be used for other operations as will be explained in the upcoming sections.

22 A New Approach

(a)

(b)

Figure 4.1: The test image and gradient magnitude weight used in the histogram cal-culation (a), traditional histogram (left) and histogram calculated using the gradient magnitude weight (right) (b).

4.2

Weighted Histogram Calculation

To give the areas with higher energy more room in the histogram, the gradient magnitude of the image can be used as a weight in the histogram calculation. A pixel with an intensity value r in the original image will contribute to the histogram with the value at the same pixel position in the weight image. The ordinary histogram calculation can be thought of as using a flat weight image that has the same value in every pixel. Here, every pixel will provide an equal count to the histogram. If instead the gradient magnitude is used as the weight, the contribution to a specific bin from each pixel will be different, depending on whether the pixel is in a region where the gradient magnitude is high or low. Figure 4.1 shows an example of using a flat weight (the traditional histogram) and using the edge magnitude as the weight. The weighted histogram calculation is

4.2 Weighted Histogram Calculation 23

described by equation 4.1 where (x, y) is the pixel position in the image.

Hweighted(r) =

Hweighted(r) + Whistogram(x, y) image(x, y) = r

Hweighted(r) image(x, y) 6= r

(4.1)

4.2.1

Weighted Global Mapping

In the global mappings, the weight will shift the focus from homogenous areas to the high-energy areas. This is, as explained earlier, of interest when the back-ground is large. However, when the object of interest contains large under- or over-exposed parts, these will sometimes be suppressed even more than the background as the background often contains noise and artifacts from the image acquisition. The weight has to be used with care and adjusted to the type of input image. Notice in figure 4.2 (a) that the lungs and the hands are low in contrast when no weight is used, this is due to the large bright area dominating the histogram. In figure 4.2 (b), a weight is used to put focus on the high energy areas and the contrast is significantly increased in those parts of the image, also see figure 4.3. However, the suppression of the low energy area is also rather significant.

(a) (b)

Figure 4.2: Selective fast-GLG with no weight (a) and selective fast-GLG with gradient magnitude weight (b).

24 A New Approach

(a) (b)

Figure 4.3: Zoomed in parts of figure 4.2.

4.2.2

Weighted Local Mapping

Using the weight in the local mappings will not suppress the homogenous areas to the same extent as in the global mappings if no other constraints are used. If the local window is centered in a homogenous area, the normalized local histogram is not much different than without the weight since the gradient magnitude will be relatively flat in that particular window. Where there is a high variation in the magnitude image however, such as the boundaries between background and the object, the weight will influence the mapping. Depending on the constraints used, the weight can have a larger influence in the homogenous areas as well, section 4.5. In figure 4.4 an LHE-image without weight and one with weight is shown.

4.3

Collimated Images

Collimators are frequently used in conventional X-ray imaging and the collimated area of an image is often dominating. In the ideal case, the collimator completely blocks out the X-rays so that no photons reach the detector, resulting in a large area with no energy in the gradient magnitude image. When this is the case, the weighted approach in the previous section can be applied since the collimated area

4.3 Collimated Images 25

(a) (b)

Figure 4.4: LHE with no weight (a) and LHE with gradient magnitude weight (b).

will not contribute to the histogram, except for the edges where the collimator and the area of interest meet.

The ideal case however, is unfortunately very rare, the collimated areas are usually low in energy but there is still enough to influence the histogram signifi-cantly since the area is often very large. Figure 4.5 illustrates this problem and figure 4.6 (b) shows the desired histogram. As can be seen, the histogram is very much influenced even by the low energy collimator.

To completely get rid of the of the collimator’s influence on the histogram, finding a mask is necessary, see figure 4.6 (a). Since the strongest edges in the image are most often the collimator edges, the gradient magnitude image can be used in the Hough Transform to find the corresponding lines. The magnitude image is used to weight the contributions in the Hough Transform as explained in section 2.2. The mask is then constructed by using the area that is spanned between the intersections of the lines. If the collimator is only located at one side of the image, the image boundary on the other side is used as the line on that side. Artifacts like those in figure 1.1 (b) will be picked up by the gradient filters and cause misleading edges in the Hough Transform. Ignoring the input image pixels close to the borders usually solves this problem.

26 A New Approach

(a) (b)

Figure 4.5: Gradient magnitude weight for a collimated image (a) and the corresponding histogram (b).

(a)

(b)

Figure 4.6: Mask obtained by the Hough Transform and the masked gradient magnitude weight (a), the corresponding histogram using the masked weight (b).

4.4 Selective Approach With Weight 27

4.4

Selective Approach With Weight

4.4.1

Finding Breakpoints

Using the selective approach described in section 3.3.1 involves deciding where on the grayscale to put the breakpoints. It is not always easy to do this by looking only at the histogram of the image. However, if a histogram without weight is compared to a weighted histogram the difference between them is higher in the homogenous areas such as the background. Taking the difference between the two histograms helps to identify a suitable breakpoint that will suppress the background, see equation 4.2 and figure 4.7. The breakpoint is most often set to suppress the background, when there is a foreground, for example a collimator, that has intensity values that are higher than the relevant parts of the image, the selective approach can be used to suppress these as well.

Hdif f erence= Hweighted− Horiginal (4.2)

(a) (b)

(c)

Figure 4.7: Histogram with no weight (a), histogram with weight raised to the power of 10 (b) and the difference between the two (c).

The difference is positive when there is a higher probability that a pixel has a specific intensity value in the weighted histogram, usually the background will be on the negative side. The difference can either be looked at manually to decide where to put the breakpoint or it can be determined automatically by for exam-ple taking the first intensity value that has a large enough positive probability difference, say 40% of the maximum positive difference. To really put focus on the strong edges between background and object, raising the weight to the power of 5-10 can help. For collimated images that still has some background in it, the mask and the masked histogram weight are used to obtain Horiginaland Hweighted.

28 A New Approach

4.4.2

Local Mapping With Selective Approach

The selective approach can be of great use also in the local methods. It not only helps to suppress the background, it also introduces a constraint on how far the local mapping will be allowed to go in the relevant areas of the image. A major problem however, occurs when the local histogram lies on both sides of a breakpoint. For example, if the histogram has a small part on one side of the breakpoint, those values will be mapped onto a relatively larger part of the grayscale than the neighbouring values on the other side of the breakpoint that can be suppressed by the larger histogram components on that side. This leads to an abrupt change in intensity in that particular area which looks very unnatural. Figure 4.8 (a) illustrates this problem which will be addressed in the next section.

(a) (b)

Figure 4.8: Zoomed in area of an selective-LHE image with a window size of 401 × 401 (a) and same area using local/global histograms with α = 0.4 (b).

4.5

Local Histograms With Global Information

Using selective approach solves the problem of a large background influencing the mapping and using gradient magnitude weight when calculating the histograms puts focus on the high-energy areas in the image. There is still the problem of noise amplification and over-contrast in the transformed image since the mappings use too much of the available grayscale. Also, the resulting image still has an unnatural look. Figure 4.9 shows two images of the same part of the evaluation image. The first is an LHE image with no constraints, the mapping uses the entire grayscale, which causes over-contrast and the noise amplification is large. The second is an LHE image with local/global histograms, the mapping is not allowed to use the entire grayscale and the over-contrast and noise amplification is not as significant.

4.5 Local Histograms With Global Information 29

(a) (b)

Figure 4.9: Zoomed in area of an LHE image with a window size of 401 × 401 (a) and same area using local/global histograms (b).

The local histogram definition in equation 3.4 makes it possible to control how much of the grayscale the local mappings are allowed to occupy. Setting α in equation 3.4 according to equation 4.3 results in a local histogram that is the same as the global histogram and when using equation 4.4, the local histogram is emphasized. α = X window /X rest (4.3) α > X window /X rest (4.4)

Another important effect of using this formulation is that the original bright-ness relationships between the different parts of the image are largely maintained, giving the transformed image a more natural look. The unwanted effect when a local histogram is on both sides of a breakpoint and causes abrupt changes in intensity is also decreased since there are now global histogram information that will make the ’tail’ of the local histogram less significant. Figure 4.8 (b) is an example where the global histogram components decreases this effect.

4.5.1

Using the Weight To Decide Local/Global-Ratio

To display different parts of the image with enough contrast, some parts may require more focus on the local histogram than others. Introducing another weight that specifies an α-value for each pixel makes this possible. Low-contrast areas that are under- or overexposed, but where there are structures of interest are examples of when this is needed. The goal is to be able to display as much

30 A New Approach

information as possible but without increasing too much noise in areas with a low energy contribution while avoiding over-contrast in regions with a large energy contribution. Figure 4.10 (a) shows an α-weight, obtained using equation 4.8, that determines the local/global-ratio for each pixel and figure 4.11 (a) shows the resulting image.

(a) (b)

Figure 4.10: α-weight (a) and upper limit weight (b).

4.5.2

Constraining the Upper Limit

The resulting image after using the α-weight will naturally span the entire grayscale. Even though the contrast has been increased in a bright area, the range might lie a bit too high up on the grayscale. It is easier to discern details at lower brightness, or luminance, values and compressing the dynamic range slightly often gives more natural look to the image. Here, we suggest doing this by introducing another weight. This weight is also based on the gradient magnitude of the image and it puts a restriction on the upper limit that the mapping of each pixel can achieve. In the test image in figure 4.11 (a), the lower part of the image is a bit too high up on the grayscale. When the upper limit weight in figure 4.10 (b) is multiplied with the result after using the α-weight, the range is slightly moved down the grayscale in the low energy areas and largely kept the same in the high energy areas which

4.6 Selecting the Weight 31

(a) (b)

Figure 4.11: Result after using α-weight (a) and the same image with the applied upper-limits (b).

results in the image in figure 4.11 (b).

4.6

Selecting the Weight

In this section, a few things to consider when choosing the different weights are mentioned.

4.6.1

The Histogram Weight

To find a suitable weight for the histogram calculation, the type of image should be considered. If it is an image with no background and the energy profile is fairly homogenous, it might not be necessary to use a weight at all. If the image contains large, underexposed areas where there is still relevant information present, the weight should not focus too much on the strongest edges in the image, for example the background and object border. This tends to lead to too much suppression of the underexposed areas. One way to define the weight is presented in equation 4.5 where G is the gradient magnitude image normalized to [0 1], δ decides how

32 A New Approach

much emphasis is put on the stronger edges. Two different histogram weights can be seen in figure 4.12

WH= Gδ (4.5)

Figure 4.12: Two different histogram weights with δ = 1 and δ = 0.5.

4.6.2

The Local/Global-Ratio Weight

The local/global-ratio weight decides how much local information is to be used and thereby how much of the grayscale the local mapping is allowed to use. A suitable initial weight to start from is the weight that gives the global histogram in all points. Based on the histogram weight, this image can easily be calculated according to equation 4.6, where h is a matrix the size of the window to be used, containing only ones.

αi=

P WH− (h ∗ WH)

P WH

(4.6)

Since αi gives the global histogram, using it without any other constraints

will give the global HE. Based on αi, the weight can be modified according to

the characteristics of the input image type. To allow underexposed areas in the image more room on the grayscale, a higher α for these areas than others might be necessary. When choosing the α-weight it needs to be considered how much focus has been put on the high-energy areas in the image. Equation 4.7 gives an example of how the weight can be defined, β < 1 results in a higher increase of local histogram emphasis for low energy areas than for high energy areas and γ increases the local emphasis linearly for all pixels, see figure 4.13. If β = 1 and

γ = 1, the initial weight is obtained.

4.6 Selecting the Weight 33

Figure 4.13: Initial α-weight with window size 401 × 401 and α-weight with β = 0.1 and

γ = 5.

Another, slightly more complicated, way to obtaining α and to give the low energy areas even more focus than the high energy ones is to use a negative of the gradient magnitude together with αi. Equation 4.8 shows how this can be done,

αi is added to make sure every pixel at least receives a ratio that results in the

global histogram. See figure 4.14 for an example.

αneg = abs(αi− max(αi))β× max(αi) × γ

α = αneg+ αi

(4.8)

Figure 4.14: Initial α-weight with window size 401 × 401 and α-weight with negative of gradient magnitude and β = 5 and γ = 5.

Even though global information is used in the alternative definition of the local histogram, the window size will still be important. The width of the histogram

34 A New Approach

within the window will generally be smaller for a small window than for a large one. Using a large amount of a small window histogram will more often lead to over-contrast than for a large window histogram.

4.6.3

The Upper Limit Weight

The upper limit weight can be constructed based on the gradient magnitude image as well. However, it needs to be low-pass filtered to avoid sharp edges since these would make the edges in the resulting image too strong, also, the upper limit cannot vary too much as the natural look has to be maintained. The parts of the image that might require an upper limit constraint to be displayed better are the ones that are low in energy and will therefore receive the lowest upper limit but are also the ones that are mapped to the higher intensities. This makes it possible to put a stronger constraint on these areas without disturbing the natural look of the image. Equation 4.9 is a suggestion for how the upper limit weight can be obtained,  controls how low the lowest upper limit will be. WH is normalized to

[0 1].

Wu= (WH∗ lpf ilter) (4.9)

Figure 4.15: Upper limit weight with  = 0.05.

4.7

Multiple Constrained Local Histogram

Equal-ization

Combining the different weights and constraints results in a new variation of the Local Histogram Equalization and will be referred to as Multiple Constrained

4.7 Multiple Constrained Local Histogram Equalization 35

Local Histogram Equalization (MCLHE). It can be summarized as LHE with the following constraints:

• A gradient magnitude weight, WH, that puts focus on the high energy areas

in the histogram.

• A selective option if there is a background that needs to be suppressed. • Collimation detection and removal using the Hough Transform or if it is

possible, the selective approach.

• A local/global combination of histograms and a weight to decide α for each pixel. This maintains a natural look and limits the local mapping from using too much of the grayscale.

• An upper limit weight that puts a maximum intensity limit that the mapping can achieve for each pixel, this acts to compress the dynamic range of the image.

• Further contrast limitation in the form of an arctan function that prevents large spikes in the histogram and thereby prevents the slope of the mapping from being too large.

Figure 4.16 shows three images with the same constraints except for ψ in equation 3.3. In the next chapter, a few different images, the results of performing MCLHE on them and the final result after processing by Sapheneias software will be presented.

36 A New Approach

(a) (b)

(c) (d)

Figure 4.16: Original thorax image (a), MCLHE with ψ = 0.001 (b), ψ = 0.01 (c) and

Chapter 5

Results

This chapter will present examples of a few different types of images and the results obtained by applying MCLHE on them. The results will be compared with a global transformation and Matlab’s inbuilt function CLAHE (Contrast Limited Adaptive Histogram Equalization) which uses the block-wise approach for local operations. Also, the result after processing by Sapheneia’s software will be presented.

Figure 5.1 (a) shows a profile thorax image that contains large underexposed areas which almost completely dominates the histogram of the image, along with the background. In figure 5.1 (b-d), a global method, Matlab’s CLAHE algorithm and a hand made LUT is shown. The MCLHE image with corresponding weights is shown in figure 5.2 which also shows a processed image. It can be seen that MCLHE retains more of the details in the upper parts of the image compared to the other methods and still displays the structures in the underexposed lower parts.

A similar type of image that illustrates the problem of having a large variation in tissue density is a full body examination. In figure 5.3 and 5.4 a full body X-ray image is shown which cannot be displayed properly using a global mapping if the different areas are to be displayed simultaneously. MCLHE retains the natural intensity composition and increases the contrast without suppressing the small histogram components, see figure 5.5 where a processed image is also shown. A good example of a small structure, with small histogram components, that is almost completely lost in the global methods and very much flattened in the CLAHE image is the thumb in the upper left part of the image, see figure 5.6.

The spine image in figure 5.7 has a more homogenous energy distribution than the previous examples. This means that the histogram weight will not affect the look of the histogram as much. Since the other weights are based on the gradient magnitude as well, they will also have a less significant effect than in the previous examples compared to if a flat weight and a fixed α-value was used. The CLAHE image illustrates the problem of using non-overlapping blocks and fusing them with interpolation, the dark bands on both sides of the spine arise from the difference in mappings between the blocks. MCLHE still produces an image with increased local contrast without suppressing the smaller histogram components, see figure

38 Results

5.8.

The collimated image in figure 5.9 shows a leg with low exposure that also has a fairly homogenous energy distribution profile. Again, CLAHE produces dark bands in the image. This is not an acceptable effect for medical images. The MCLHE image can be seen in figure 5.10. Figure 5.11 shows a cropped version of the collimated veterinarian X-ray image in figure 1.1 (b) along with a global mapping and figure 5.12 shows the MCLHE image along with the processed image.

39

(a) (b)

(c) (d)

Figure 5.1: Original image (a), Global Selective Fast-GLG (b), Matlab’s CLAHE with 2 × 2 tiles and ’exponential’ set to 2 (c) and manual mapping (d) .

40 Results

(a)

(b) (c)

Figure 5.2: Histogram weight, α-weight and upper limit weight (a) used to obtain the MCLHE image (b) and MCLHE image processed by Sapheneia’s commercial software (c).

41

(a) (b)

42 Results

(a) (b)

Figure 5.4: Matlab’s CLAHE with 2 × 2 tiles and ’exponential’ set to 2 (a) and manual mapping (b).

43

(a) (b)

Figure 5.5: MCLHE image (a) and MCLHE image processed by Sapheneia’s commercial software (b).

44 Results

(a) (b) (c)

Figure 5.6: Zoomed in parts of the thumb in the body image. MCLHE (a), Global Selective Fast-GLG (b) and Matlab’s CLAHE (c).

45

(a) (b)

(c) (d)

Figure 5.7: Spine image (a), Global Selective Fast-GLG (b), Matlab’s CLAHE (c) and Manual mapping (d).

46 Results

(a)

(b)

Figure 5.8: MCLHE image (a) and MCLHE image processed by Sapheneia’s commercial software (b).

47

(a) (b)

(c) (d)

Figure 5.9: Leg image (a), Global Selective Fast-GLG (b), Matlab’s CLAHE (c) and Manual mapping (d).

48 Results

(a)

(b)

Figure 5.10: MCLHE image (a) and MCLHE image processed by Sapheneia’s commer-cial software (b).

49

(a)

(b)

50 Results

(a)

(b)

Figure 5.12: MCLHE image (a) and MCLHE image processed by Sapheneia’s commer-cial software (b).

Chapter 6

Discussion

To be able to obtain a proper contrast setting for all of the different regions in an image, a local approach is necessary. The block-wise approach is most often not a good candidate since it is difficult to find suitable block sizes and the interpolation can cause artifacts that are unacceptable for diagnostic viewing. The sliding window approach is the only method that makes all pixels the center pixel in their respective surrounding, at the cost of increased computation time.

The proposed method solves the problem of displaying under- and overexposed areas simultaneously while maintaining the natural composition of the image as much as possible. However, it is difficult to find proper weights and constraints without having to tune them for a specific image type. The more parameters need-ing tunneed-ing the less robust the method will be. Findneed-ing a good set of parameters for specific image types seems to be necessary.

For collimated images it is important to be able to find a mask that excludes the irrelevant areas. This is fairly easy when the images contain distinct collimator edges without artifacts that cause misleading edges which can be picked up in the hough transform. In some images however, the collimator edges can be difficult to detect since they are not always as dominating in the gradient magnitude images as in the examples presented here.

6.1

Suggestions For Further Work

To improve the proposed method there a few things to focus on.

• To get a good result, the weight parameters need to be tuned for the specific image or image type. Identifying a good set of parameters for each image type would be desirable. To do this, the method needs to be tested on a wider set of images for each examination type, e.g more full body examination images, more spine images, more extremity images and so on.

• The square window that is used to obtain the local histograms might not be the best shape since it is not an ideal representation of a pixel’s surrounding.

52 Discussion

A Gaussian window is probably a better candidate, as it would put more weight on a pixel’s closest neighbours. However, the attractive feature of the square window, where the local histogram can simply be updated when moving from one pixel to another instead of having to calculate the entire histogram again, will be lost. A possible approach is to use two or more windows of different sizes, which can be weighted differently along with the global histogram, to obtain the local histograms.

6.2

Conclusion

We have presented a modification of local histogram equalization to make it suit-able for preparing digital X-ray images for diagnostic viewing. By using the gra-dient magnitude of the image as a weight in the histogram calculations and in deciding how much of the local histogram to be used in the equalization procedure we enable a local contrast setting within a suitable range. This means that the contrast setting in one region of the image does not influence the contrast setting in another while at the same time maintaining the natural look of the image as a whole. Also, a way of eliminating the influence of collimators using the hough transform has been presented. The results have been compared to existing global and local methods and show promising results.

Bibliography

[1] Zhi Yu Chen, B.R. Abidi, D.L. Page, and M.A. Abidi. A generalized and automatic image contrast enhancement using gray level grouping. In Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on, volume 2, pages II–II, May 2006.

[2] Gösta H. Granlund and Hans Knutsson. Signal Processing for Computer Vi-sion. Kluwer Academic Publishers, 1995.

[3] Iyad Jafar and Hao Ying. New algorithms for contrast enhancement in grayscale images based on the variational definition of histogram equalization. Integr. Comput.-Aided Eng., 15(2):131–147, 2008.

[4] Bernd Jähne. Digital Image Processing. Springer-Verlag Berlin Heidelberg, 2005.

[5] Stephen M. Pizer, E. Philip Amburn, John D. Austin, Robert Cromartie, Ari Geselowitz, Trey Greer, Bart Ter Haar Romeny, and John B. Zimmerman. Adaptive histogram equalization and its variations. Comput. Vision Graph. Image Process., 39(3):355–368, 1987.

[6] M J Yaffe and J A Rowlands. X-ray detectors for digital radiography. Physics in Medicine and Biology, 42(1):1–39, 1997.

[7] Hui Zhu, Francis H. Y. Chan, and F. K. Lam. Image contrast enhancement by constrained local histogram equalization. Comput. Vis. Image Underst., 73(2):281–290, 1999.