MA3111: Mathematical Image Processing Intensity Transformations & Spatial Filtering

(1)

MA3111: Mathematical Image Processing Intensity Transformations & Spatial Filtering

Suh-Yuh Yang (楊肅煜)

Department of Mathematics, National Central University Jhongli District, Taoyuan City 32001, Taiwan

[email protected]

http://www.math.ncu.edu.tw/ ∼ syyang/

Suh-Yuh Yang (楊肅煜), Math. Dept., NCU, Taiwanc Intensity Transformations & Spatial Filtering – 1/41

(2)

Spatial domain and transform domain

The spatial domain approach and transform domain approach are two main categories in image processing:

Spatial domain:

refers to the image plane itself, and image processing methods in this category are based on direct manipulation of pixels in an image.

Transform domain:

involves first transforming an image into the transform domain, doing the processing there, and obtaining the inverse transform to bring the results back into spatial domain.

(3)

Outline of “intensity transformations & spatial filtering”

In this lecture, we will discuss a number of classical techniques for two principal categories of spatial domain processing:

Intensity transformations:

operate on single pixels of an image for tasks such as contrast manipulation and image thresholding.

Spatial filtering:

performs operations on the neighborhood of every pixel in an image. Examples of spatial filtering include image smoothing and sharpening.

The material of this lecture is based on Chapter 3 in [GW2018].

(4)

Spatial domain process

The spatial domain process is generally posed in the form:

g

(

x, y

) =

T

(

f

(

x, y

))

,

where f ( x, y ) is an input image, g ( x, y ) is the output image, and T is an operator on f defined over a neighborhood (typically a rectangle) of point ( x, y ) .

A 3 × 3 neighborhood about the point ( x 0 , y 0 ) . The neighborhood is moved from pixel to pixel in the image to generate the output image.

(5)

Spatial filtering and intensity transformation

A smoothing spatial filter T: suppose that the neighborhood is a square of size 3 × 3 and that operator T is defined as

compute the average intensity of the pixels in the neighborhood. Then T is a smoothing filter.

Consider an arbitrary location in an image f , say ( 100, 150 ) . Then

g ( 100, 150 ) = T f ( 100, 150 ) = ¹ 9

∑ 1 i=−1

∑ 1 j=−1

f ( 100 − i, 150 − j ) .

(A neighborhood processing technique)

Intensity transformation: The smallest possible neighborhood is of size 1 × 1. T becomes an intensity transformation of the form

g

(

x, y

) =

: s

=

T

(

r

)

:

=

T

(

f

(

x, y

))

. (A point processing technique)

(6)

Intensity transformation functions

Contrast stretching function:

“left figure”

produces an image of higher contrast than the original, by darkening the intensity levels below k and brightening the levels above k.

Thresholding function: In the limiting case shown in

“right figure,”

T ( r ) produces a two level (binary) image.

3.1

Background

121 the value of the output image g at those coordinates is equal to the result of applying T to the neighborhood with origin at ( , x y

₀ ₀

) in f. For example, suppose that the neighborhood is a square of size 3 3 × and that operator T is defined as “compute the average intensity of the pixels in the neighborhood.” Consider an arbitrary location in an image, say ( 100 150 The result at that location in the output image, , ).

g( 100 150 is the sum of f ( , ), 100 150 and its 8-neighbors, divided by 9. The center of , ) the neighborhood is then moved to the next adjacent location and the procedure is repeated to generate the next value of the output image g. Typically, the process starts at the top left of the input image and proceeds pixel by pixel in a horizontal (vertical) scan, one row (column) at a time. We will discuss this type of neighborhood processing beginning in Section 3.4.

The smallest possible neighborhood is of size 1 1 × . In this case, g depends only on the value of f at a single point ( , ) x y and T in Eq. (3-1) becomes an intensity (also called a gray-level, or mapping) transformation function of the form

s = ( ) T r (3-2)

where, for simplicity in notation, we use s and r to denote, respectively, the intensity of g and f at any point ( , ). x y For example, if T r ( ) has the form in Fig. 3.2(a), the result of applying the transformation to every pixel in f to generate the corresponding pixels in g would be to produce an image of higher contrast than the original, by darkening the intensity levels below k and brightening the levels above k. In this technique, sometimes called contrast stretching (see Section 3.2), values of r lower than k reduce (darken) the values of s, toward black. The opposite is true for values of r higher than k. Observe how an intensity value r

₀

is mapped to obtain the corresponding value s

₀

. In the limiting case shown in Fig. 3.2(b), T r ( ) produces a two- level (binary) image. A mapping of this form is called a thresholding function. Some fairly simple yet powerful processing approaches can be formulated with intensity transformation functions. In this chapter, we use intensity transformations princi- pally for image enhancement. In Chapter 10, we will use them for image segmenta- tion. Approaches whose results depend only on the intensity at a point sometimes are called point processing techniques, as opposed to the neighborhood processing techniques discussed in the previous paragraph.

Depending on the size of a neighborhood and its location, part of the neighborhood may lie outside the image. There are two solutions to this:

(1) to ignore the values outside the image, or (2) to pad image, as discussed in Section 3.4.

The second approach is preferred.

b a

FIGURE 3.2

Intensity transformation functions.

(a) Contrast stretching function.

(b) Thresholding function.

k k r0

s0 & T(r0)

Dark Light

r r

s & T(r) s & T(r)

T(r) T(r)

DIP4E_GLOBAL_Print_Ready.indb 121 6/16/2017 2:03:05 PM

www.EBooksWorld.ir

g

(

x, y

) =

: s

=

T

(

r

)

:

=

T

(

f

(

x, y

))

(7)

Some examples: g ( x, y ) = : s = T ( r ) : = T ( f ( x, y ))

Negative transformation: The negative of an image with intensity levels in the range [ 0, L − 1 ] is obtained by s = L − 1 − r.

Log transformation: s = c log ( 1 + r ) , where c > 0 is a constant.

Power-law (gamma) transformation: s = cr ^γ or s = c ( r + ε ) ^γ , where c and γ are positive constants.

Piecewise linear transformation

3.2 Some Basic Intensity Transformation Functions 123

Reversing the intensity levels of a digital image in this manner produces the equivalent of a photographic negative. This type of processing is used, for example, in enhancing white or gray detail embedded in dark regions of an image, especially when the black areas are dominant in size. Figure 3.4 shows an example. The original image is a digital mammogram showing a small lesion. Despite the fact that the visual content is the same in both images, some viewers find it easier to analyze the fine details of the breast tissue using the negative image.

Identity

0 L/4 L/2 3L/4 L ! 1

Input intensity levels, r 0

L/4 L/2 3L/4 L ! 1

Output intensity levels, s

Log Negative

nth power nth root

Inverse log (exponential) FIGURE 3.3

Some basic intensity transformation functions. Each curve was scaled independently so that all curves would fit in the same graph. Our interest here is on the shapes of the curves, not on their relative values.

b a FIGURE 3.4 (a) A digital mammogram.

(b) Negative image obtained using Eq. (3-3).

(Image (a) Courtesy of General Electric Medical Systems.)

www.EBooksWorld.ir

3.2 Some Basic Intensity Transformation Functions 125 interval [ ,0 255 and showing the spectrum in the same 8-bit display. The level of ] detail visible in this image as compared to an unmodified display of the spectrum is evident from these two images. Most of the Fourier spectra in image processing publications, including this book, have been scaled in this manner.

POWER-LAW (GAMMA) TRANSFORMATIONS Power-law transformations have the form

s=cr^g (3-5)

where c and g are positive constants. Sometimes Eq. (3-5) is written as s=c r(+^e)^g to account for offsets (that is, a measurable output when the input is zero). However, offsets typically are an issue of display calibration, and as a result they are normally ignored in Eq. (3-5). Figure 3.6 shows plots of s as a function of r for various values of g. As with log transformations, power-law curves with fractional values of g map a narrow range of dark input values into a wider range of output values, with the opposite being true for higher values of input levels. Note also in Fig. 3.6 that a fam- ily of transformations can be obtained simply by varying g. Curves generated with values of g> 1 have exactly the opposite effect as those generated with values of g< 1. When c = =g 1 Eq. (3-5) reduces to the identity transformation.

The response of many devices used for image capture, printing, and display obey a power law. By convention, the exponent in a power-law equation is referred to as gamma [hence our use of this symbol in Eq. (3-5)]. The process used to correct these power-law response phenomena is called gamma correction or gamma encoding.

For example, cathode ray tube (CRT) devices have an intensity-to-voltage response that is a power function, with exponents varying from approximately 1.8 to 2.5. As the curve for g= 2 5. in Fig. 3.6 shows, such display systems would tend to produce

g & 0.04 g & 0.10

g & 0.20 g & 0.40

g & 0.67 g & 1

g & 1.5 g & 2.5

g & 5.0 g & 10.0

g & 25.0

0 L/4 L/2 3L/4 L ! 1

Input intensity levels,r 0

L/4 L/2 3L/4 L ! 1

FIGURE 3.6 Plots of the gamma equation s=c r^g for various values of g (c = 1 in all cases). Each curve was scaled independently so that all curves would fit in the same graph. Our interest here is on the shapes of the curves, not on their relative values.

www.EBooksWorld.ir

3.2 Some Basic Intensity Transformation Functions 129 Contrast Stretching

Low-contrast images can result from poor illumination, lack of dynamic range in the imaging sensor, or even the wrong setting of a lens aperture during image acquisi- tion. Contrast stretching expands the range of intensity levels in an image so that it spans the ideal full intensity range of the recording medium or display device.

Figure 3.10(a) shows a typical transformation used for contrast stretching. The locations of points ( , )r s1 1 and ( , )r s2 2 control the shape of the transformation function.

If r₁= and rs₁ ₂= the transformation is a linear function that produces no changes s₂ in intensity. If r₁= , sr₂ ₁= , and s0 ₂=L− the transformation becomes a threshold-1 ing function that creates a binary image [see Fig. 3.2(b)]. Intermediate values of ( , )r s1 1

and ( , )s r₂ ₂ produce various degrees of spread in the intensity levels of the output image, thus affecting its contrast. In general, r₁≤ and sr₂ ₁≤s₂ is assumed so that the function is single valued and monotonically increasing. This preserves the order of intensity levels, thus preventing the creation of intensity artifacts. Figure 3.10(b) shows an 8-bit image with low contrast. Figure 3.10(c) shows the result of contrast stretching, obtained by setting ( , )r s₁ ₁ =(r_min, )0 and ( ,r s₂ ₂)=(r_max,L−1), where r_min and r_max denote the minimum and maximum intensity levels in the input image,

0 L/4 L/2 3L/4 L ! 1

Input intensities, r 0

L/4 L/2 3L/4 L ! 1

Output intensities, s

(r2, s2)

(r1, s1) T(r) b

a d c FIGURE 3.10 Contrast stretching.

(a) Piecewise linear transformation function. (b) A low- contrast electron microscope image of pollen, magnified 700 times.

(c) Result of contrast stretching.

(d) Result of thresholding.

(Original image courtesy of Dr.

Roger Heady, Research School of Biological Sciences, Australian National University, Canberra, Australia.)

www.EBooksWorld.ir Suh-Yuh Yang (楊肅煜), Math. Dept., NCU, Taiwanc Intensity Transformations & Spatial Filtering – 7/41

(8)

Negative images and log images

3.2 Some Basic Intensity Transformation Functions 123

Identity

0 L/4 L/2 3L/4 L ! 1

L/4 L/2 3L/4 L ! 1

Log Negative

www.EBooksWorld.ir

Identity

0 L/4 L/2 3L/4 L ! 1

L/4 L/2 3L/4 L ! 1

Log Negative

www.EBooksWorld.ir

Negative image of a digital mammogram (X-ray)

124 Chapter3 Intensity Transformations and Spatial Filtering LOG TRANSFORMATIONS

The general form of the log transformation in Fig. 3.3 is

s=clog(1+r) (3-4)

where c is a constant and it is assumed that r≥ 0. The shape of the log curve in Fig. 3.3 shows that this transformation maps a narrow range of low intensity values in the input into a wider range of output levels. For example, note how input levels in the range [ ,0L4] map to output levels to the range [ ,0 3L4]. Conversely, higher values of input levels are mapped to a narrower range in the output. We use a transformation of this type to expand the values of dark pixels in an image, while compressing the higher-level values. The opposite is true of the inverse log (exponential) transformation.

Any curve having the general shape of the log function shown in Fig. 3.3 would accomplish this spreading/compressing of intensity levels in an image, but the power-law transformations discussed in the next section are much more versatile for this purpose. The log function has the important characteristic that it compresses the dynamic range of pixel values. An example in which pixel values have a large dynamic range is the Fourier spectrum, which we will discuss in Chapter 4. It is not unusual to encounter spectrum values that range from 0 to 10⁶ or higher. Processing numbers such as these presents no problems for a computer, but image displays can- not reproduce faithfully such a wide range of values. The net effect is that intensity detail can be lost in the display of a typical Fourier spectrum.

Figure 3.5(a) shows a Fourier spectrum with values in the range 0 to 1 5 10. × ⁶. When these values are scaled linearly for display in an 8-bit system, the brightest pixels dominate the display, at the expense of lower (and just as important) values of the spectrum. The effect of this dominance is illustrated vividly by the relatively small area of the image in Fig. 3.5(a) that is not perceived as black. If, instead of displaying the values in this manner, we first apply Eq. (3-4) (with c= 1 in this case) to the spectrum values, then the range of values of the result becomes 0 to 6.2. Trans- forming values in this way enables a greater range of intensities to be shown on the display. Figure 3.5(b) shows the result of scaling the intensity range linearly to the b

a FIGURE 3.5 (a) Fourier spectrum displayed as a grayscale image.

(b) Result of applying the log transformation in Eq. (3-4) with c= 1. Both images are scaled to the range [0, 255].

www.EBooksWorld.ir

124 Chapter3 Intensity Transformations and Spatial Filtering LOG TRANSFORMATIONS

The general form of the log transformation in Fig. 3.3 is

s=clog(1+r) (3-4)

where c is a constant and it is assumed that r≥ 0. The shape of the log curve in Fig. 3.3 shows that this transformation maps a narrow range of low intensity values in the input into a wider range of output levels. For example, note how input levels in the range [ ,0L4] map to output levels to the range [ ,0 3L4]. Conversely, higher values of input levels are mapped to a narrower range in the output. We use a transformation of this type to expand the values of dark pixels in an image, while compressing the higher-level values. The opposite is true of the inverse log (exponential) transformation.

Any curve having the general shape of the log function shown in Fig. 3.3 would accomplish this spreading/compressing of intensity levels in an image, but the power-law transformations discussed in the next section are much more versatile for this purpose. The log function has the important characteristic that it compresses the dynamic range of pixel values. An example in which pixel values have a large dynamic range is the Fourier spectrum, which we will discuss in Chapter 4. It is not unusual to encounter spectrum values that range from 0 to 10⁶ or higher. Processing numbers such as these presents no problems for a computer, but image displays can- not reproduce faithfully such a wide range of values. The net effect is that intensity detail can be lost in the display of a typical Fourier spectrum.

Figure 3.5(a) shows a Fourier spectrum with values in the range 0 to 1 5 10. × ⁶. When these values are scaled linearly for display in an 8-bit system, the brightest pixels dominate the display, at the expense of lower (and just as important) values of the spectrum. The effect of this dominance is illustrated vividly by the relatively small area of the image in Fig. 3.5(a) that is not perceived as black. If, instead of displaying the values in this manner, we first apply Eq. (3-4) (with c= 1 in this case) to the spectrum values, then the range of values of the result becomes 0 to 6.2. Trans- forming values in this way enables a greater range of intensities to be shown on the display. Figure 3.5(b) shows the result of scaling the intensity range linearly to the b

a FIGURE 3.5 (a) Fourier spectrum displayed as a grayscale image.

(b) Result of applying the log transformation in Eq. (3-4) with c= 1. Both images are scaled to the range [0, 255].

www.EBooksWorld.ir

Log transformation of Fourier spectrum with c = 1

(9)

Images of gamma transformation

Many devices used for image capture, printing, and display obey a power law, e.g., cathode ray tube (陰極射線管) ¹²⁶

^Chapter³Intensity Transformations and Spatial Filtering

images that are darker than intended. Figure 3.7 illustrates this effect. Figure 3.7(a) is an image of an intensity ramp displayed in a monitor with a gamma of 2.5. As expected, the output of the monitor appears darker than the input, as Fig. 3.7(b) shows.

In this case, gamma correction consists of using the transformation s r= ^{1 2 5}^. =r^{0 4}^. to preprocess the image before inputting it into the monitor. Figure 3.7(c) is the result.

When input into the same monitor, the gamma-corrected image produces an output that is close in appearance to the original image, as Fig. 3.7(d) shows. A similar analysis as above would apply to other imaging devices, such as scanners and printers, the dif- ference being the device-dependent value of gamma (Poynton [1996]).

EXAMPLE 3.1 : Contrast enhancement using power-law intensity transformations.

In addition to gamma correction, power-law transformations are useful for general-purpose contrast manipulation. Figure 3.8(a) shows a magnetic resonance image (MRI) of a human upper thoracic spine with a fracture dislocation. The fracture is visible in the region highlighted by the circle. Because the image is predominantly dark, an expansion of intensity levels is desirable. This can be accomplished using a power-law transformation with a fractional exponent. The other images shown in the figure were obtained by processing Fig. 3.8(a) with the power-law transformation function of Eq. (3-5). The values

Sometimes, a higher gamma makes the displayed image look better to viewers than the original because of an increase in contrast.

However, the objective of gamma correction is to produce a faithful display of an input image.

b ad c FIGURE 3.7 (a) Intensity ramp image. (b) Image as viewed on a simulated monitor with a gamma of 2.5. (c) Gamma- corrected image.

(d) Corrected image as viewed on the same monitor. Compare (d) and (a).

Original image as viewed on a monitor with a gamma of 2.5

Original image Gamma Correction

Gamma-corrected image Gamma-corrected image as viewed on the same monitor

www.EBooksWorld.ir

Intensity ramp images with c = _{1, γ} = 2.5 and correction s = r ^1/(2.5)

(10)

Gamma transformation: MRI of a fractured human spine

of gamma corresponding to images (b) through (d) are 0.6, 0.4, and 0.3, respectively (c= 1 in all cases).

Observe that as gamma decreased from 0.6 to 0.4, more detail became visible. A further decrease of gamma to 0.3 enhanced a little more detail in the background, but began to reduce contrast to the point where the image started to have a very slight “washed-out” appearance, especially in the background.

The best enhancement in terms of contrast and discernible detail was obtained with g= 0 4. . A value of g= 0 3. is an approximate limit below which contrast in this particular image would be reduced to an unacceptable level.

EXAMPLE 3.2 : Another illustration of power-law transformations.

Figure 3.9(a) shows the opposite problem of that presented in Fig. 3.8(a). The image to be processed b

a d c FIGURE 3.8 (a) Magnetic resonance image (MRI) of a fractured human spine (the region of the fracture is enclosed by the circle).

(b)–(d) Results of applying the transformation in Eq. (3-5) with c= 1 and g= 0 6. , 0.4, and 0.3, respectively.

David R. Pickens, Department of Radiology and Radiological Sciences, Vanderbilt University Medical Center.)

www.EBooksWorld.ir

Region of the fracture is enclosed by the circle: c = 1, γ = 0.6, 0.4, 0.3

(11)

Gamma transformation: aerial images (空空空拍拍拍影影影像像像)

128

Chapter3 Intensity Transformations and Spatial Filtering b

ad c FIGURE 3.9 (a) Aerial image.

(b)–(d) Results of applying the transformation in Eq. (3-5) with g= 3 0. , 4.0, and 5.0, respectively.

(c= 1 in all cases.) (Original image courtesy of NASA.)

now has a washed-out appearance, indicating that a compression of intensity levels is desirable. This can be accomplished with Eq. (3-5) using values of g greater than 1. The results of processing Fig. 3.9(a) with g= 3 0. , 4.0, and 5.0 are shown in Figs. 3.9(b) through (d), respectively. Suitable results were obtained using gamma values of 3.0 and 4.0. The latter result has a slightly more appealing appearance because it has higher contrast. This is true also of the result obtained with g= 5 0. . For example, the airport runways near the middle of the image appears clearer in Fig. 3.9(d) than in any of the other three images.

PIECEWISE LINEAR TRANSFORMATION FUNCTIONS

An approach complementary to the methods discussed in the previous three sec- tions is to use piecewise linear functions. The advantage of these functions over those discussed thus far is that the form of piecewise functions can be arbitrarily complex.

In fact, as you will see shortly, a practical implementation of some important transformations can be formulated only as piecewise linear functions. The main disadvan- tage of these functions is that their specification requires considerable user input.

www.EBooksWorld.ir

c = 1, γ = 3.0, 4.0, 5.0

(12)

Piecewise linear transformation: contrast stretching

3.2 Some Basic Intensity Transformation Functions 129 Contrast Stretching

Figure 3.10(a) shows a typical transformation used for contrast stretching. The locations of points ( , )r s₁₁ and ( , )r s₂ ₂ control the shape of the transformation function.

If r1= and rs1 2= the transformation is a linear function that produces no changes s2 in intensity. If r₁= , sr₂ 1= , and s0 ₂=L− the transformation becomes a threshold-1 ing function that creates a binary image [see Fig. 3.2(b)]. Intermediate values of ( , )r s₁ ₁ and ( , )s r2 2 produce various degrees of spread in the intensity levels of the output image, thus affecting its contrast. In general, r₁≤ and sr₂ ₁≤s₂ is assumed so that the function is single valued and monotonically increasing. This preserves the order of intensity levels, thus preventing the creation of intensity artifacts. Figure 3.10(b) shows an 8-bit image with low contrast. Figure 3.10(c) shows the result of contrast stretching, obtained by setting ( , ) (r s₁ ₁ =r_min, )0 and ( , ) (r s₂ ₂ =r_max,L−1), where rmin and rmax denote the minimum and maximum intensity levels in the input image,

0 L/4 L/2 3L/4 L ! 1

L/4 L/2 3L/4 L ! 1

(r₂, s₂)

(r₁, s₁) T(r) b

ad c FIGURE 3.10 Contrast stretching.

www.EBooksWorld.ir

Figure 3.10(a) shows a typical transformation used for contrast stretching. The locations of points ( , )r s₁ ₁ and ( , )r s₂ ₂ control the shape of the transformation function.

If r₁= and rs₁ ₂= the transformation is a linear function that produces no changes s₂ in intensity. If r1= , sr2 1= , and s0 2=L− the transformation becomes a threshold-1 ing function that creates a binary image [see Fig. 3.2(b)]. Intermediate values of ( , )r s₁ ₁ and ( , )s r₂ ₂ produce various degrees of spread in the intensity levels of the output image, thus affecting its contrast. In general, r1≤ and sr2 1≤s2 is assumed so that the function is single valued and monotonically increasing. This preserves the order of intensity levels, thus preventing the creation of intensity artifacts. Figure 3.10(b) shows an 8-bit image with low contrast. Figure 3.10(c) shows the result of contrast stretching, obtained by setting ( , ) (r s₁ ₁=r_min, )0 and ( , ) (r s₂ ₂ =r_max,L−1), where rmin and rmax denote the minimum and maximum intensity levels in the input image,

0 L/4 L/2 3L/4 L ! 1

L/4 L/2 3L/4 L ! 1

(r₂, s₂)

(r₁, s₁) T(r) b

www.EBooksWorld.ir

Figure 3.10(a) shows a typical transformation used for contrast stretching. The locations of points ( , )r s₁ ₁ and ( , )r s₂ ₂ control the shape of the transformation function.

If r₁= and rs₁ ₂= the transformation is a linear function that produces no changes s₂ in intensity. If r1= , sr2 1= , and s0 2=L− the transformation becomes a threshold-1 ing function that creates a binary image [see Fig. 3.2(b)]. Intermediate values of ( , )r s₁ ₁ and ( , )s r₂ ₂ produce various degrees of spread in the intensity levels of the output image, thus affecting its contrast. In general, r1≤ and sr2 1≤s2 is assumed so that the function is single valued and monotonically increasing. This preserves the order of intensity levels, thus preventing the creation of intensity artifacts. Figure 3.10(b) shows an 8-bit image with low contrast. Figure 3.10(c) shows the result of contrast stretching, obtained by setting ( , ) (r s₁ ₁ =r_min, )0 and ( , ) (r s₂ ₂ =r_max,L−1), where rmin and rmax denote the minimum and maximum intensity levels in the input image,

0 L/4 L/2 3L/4 L ! 1

L/4 L/2 3L/4 L ! 1

(r₂, s₂)

(r₁, s₁) T(r) b

www.EBooksWorld.ir

A low-contrast electron microscope image of pollen; Result of contrast stretching; Result of thresholding

3.2 Some Basic Intensity Transformation Functions 129 Contrast Stretching

Figure 3.10(a) shows a typical transformation used for contrast stretching. The locations of points ( , )r s1 1 and ( , )r s2 2 control the shape of the transformation function.

If r₁= and rs₁ ₂= the transformation is a linear function that produces no changes s₂ in intensity. If r₁= , sr₂ 1= , and s0 ₂=L− the transformation becomes a threshold-1 ing function that creates a binary image [see Fig. 3.2(b)]. Intermediate values of ( , )r s1 1

and ( , )s r₂ ₂ produce various degrees of spread in the intensity levels of the output image, thus affecting its contrast. In general, r1≤ and sr2 1≤s2 is assumed so that the function is single valued and monotonically increasing. This preserves the order of intensity levels, thus preventing the creation of intensity artifacts. Figure 3.10(b) shows an 8-bit image with low contrast. Figure 3.10(c) shows the result of contrast stretching, obtained by setting ( , )r s₁ ₁ =(r_min, )0 and ( ,r s₂ ₂)=(r_max,L−1), where r_min and r_max denote the minimum and maximum intensity levels in the input image,

0 L/4 L/2 3L/4 L ! 1

L/4 L/2 3L/4 L ! 1

(r2, s2)

(r1, s1) T(r) b

ad c

FIGURE 3.10 Contrast stretching.

www.EBooksWorld.ir

3.1

Background

121 the value of the output image g at those coordinates is equal to the result of applying T to the neighborhood with origin at (

x y₀

,

₀

) in f. For example, suppose that the neighborhood is a square of size 3 × and that operator T is defined as “com- 3 pute the average intensity of the pixels in the neighborhood.” Consider an arbitrary location in an image, say ( 100 150 The result at that location in the output image, , ).

g(

100 150 is the sum of f ( , ), 100 150 and its 8-neighbors, divided by 9. The center of , ) the neighborhood is then moved to the next adjacent location and the procedure is repeated to generate the next value of the output image g. Typically, the process starts at the top left of the input image and proceeds pixel by pixel in a horizontal (vertical) scan, one row (column) at a time. We will discuss this type of neighborhood processing beginning in Section 3.4.

The smallest possible neighborhood is of size 1 1 × . In this case, g depends only on the value of f at a single point ( , )

x y

and T in Eq. (3-1) becomes an intensity (also called a gray-level, or mapping) transformation function of the form

s

= ( )

T r (3-2)

where, for simplicity in notation, we use s and r to denote, respectively, the intensity of g and f at any point ( , ).

x y

For example, if T r ( ) has the form in Fig. 3.2(a), the result of applying the transformation to every pixel in f to generate the corresponding pixels in g would be to produce an image of higher contrast than the original, by darkening the intensity levels below k and brightening the levels above k. In this technique, sometimes called contrast stretching (see Section 3.2), values of r lower than k reduce (darken) the values of s, toward black. The opposite is true for values of r higher than k. Observe how an intensity value r

₀

is mapped to obtain the corresponding value s

₀

. In the limiting case shown in Fig. 3.2(b), T r ( ) produces a two- level (binary) image. A mapping of this form is called a thresholding function. Some fairly simple yet powerful processing approaches can be formulated with intensity transformation functions. In this chapter, we use intensity transformations princi- pally for image enhancement. In Chapter 10, we will use them for image segmenta- tion. Approaches whose results depend only on the intensity at a point sometimes are called point processing techniques, as opposed to the neighborhood processing techniques discussed in the previous paragraph.

Depending on the size of a neighborhood and its location, part of the neighborhood may lie outside the image. There are two solutions to this:

(1) to ignore the values outside the image, or (2) to pad image, as discussed in Section 3.4.

The second approach is preferred.

b a

FIGURE 3.2 Intensity transformation functions.

(a) Contrast stretching function.

(b) Thresholding function.

k k r0

s0 & T (r0)

DarkLight

Dark Light

r r

s & T (r) s & T (r)

T(r) T(r)

www.EBooksWorld.ir

Thresholding function: r 1 = r 2 = k, s 1 = 0, s 2 = L − 1

(13)

Intensity-level slicing (強強強度度度準準準位位位切切切片片片)

Intensity-level slicing is to highlight a specific range of intensities in an image,

e.g., enhancing features in satellite imagery such as masses of water, and enhancing flaws in X-ray images.

One approach is to display in one value (say, white) all the values in the range of interest and in another (say, black) all other intensities, i.e., produces a binary image.

The second approach brightens (or darkens) the desired range of intensities, but leaves all other intensity levels in the image unchanged.

(L) first approach; (R) second approach

(14)

Examples of the intensity-level slicing

(L) aortic angiogram; (M) first approach; (R) second approach, with the selected range set near black

(15)

Histogram (直直直方方方圖圖圖)

Let r k , for k = 0, 1, · · · , L − 1, denote the intensities of an L-level image f ( x, y ) . The unnormalized histogram of f is defined as

h

(

r

_k ) =

n

_k

, k

=

0, 1,

· · ·

, L

−

1,

where n _k is the number of pixels in f with intensity r _k .

The normalized histogram (image histogram) of f is defined as

p

(

r

k ) =

^h

(

r

k )

MN

=

ⁿ

^k

MN,

where f is an M × N image. That is, p ( r _k ) is the probability of intensity level r _k occurring in an image. Then

L−1 ∑

k=0

p

(

r

_k ) =

1.

Histograms are simple to compute and are also suitable for fast hardware implementations, thus making histogram-based techniques a popular tool for real-time image processing.

(16)

Four image types and their corresponding histograms

The horizontal axis of the histograms are values of r _k and the vertical axis are values of p ( r _k )

(17)

Intensity transformation for histogram

Let the variable r denote the intensities of an image to be processed.

Assume that r ∈ [ 0, L − 1 ] with r = 0 representing black and r = L − 1 representing white. We consider the intensity transformation

s

=

T

(

r

)

, 0

≤

r

≤

L

−

1.

For a given intensity value r in the input image, T produces an output intensity value s. We assume that

T

(

r

)

is a monotonic increasing function in the interval

[

0, L

−

1

]

. T

(

r

) ∈ [

0, L

−

1

]

for all r

∈ [

0, L

−

1

]

.

If we need to use the inverse r

=

T

⁻¹ (

s

)

, s

∈

Range

(

T

)

, then we assume T

(

r

)

is a strictly monotonic increasing function in

[

0, L

−

1

]

.

(18)

Histogram equalization (HE): g ( x, y ) = : s = T ( r ) : = T ( f ( x, y ))

We are given a grayscale image f : Ω → [ 0, 1 ] . The cumulative histogram

(cumulative distribution function)

T is defined by considering f as a random variable: for η ∈ [ 0, 1 ] , we define

T

( η )

:

=

Prob

(

f

≤ η )

=

¹

|

_Ω

|

{(

x, y

) ∈

_{Ω : f}

(

x, y

) ≤ η }

_.

Then T : [ 0, 1 ] → [ 0, 1 ] is a monotonic increasing function.

The histogram equalized image g : Ω → [ 0, 1 ] is obtained by defining

g

(

x, y

)

:

=

T

(

f

(

x, y

))

.

(19)

Histogram equalized image g ∼ U ( 0, 1 ) if T is invertible

If T is strictly increasing, then T is invertible and the

cumulative distribution function

of the histogram equalized image g is

Prob

(

g

≤ η ) =

Prob

(

T

(

f

) ≤ η ) =

Prob

(

f

≤

T

⁻¹ ( η ))

=

T

(

T

⁻¹ ( η )) = η.

Hence, the

probability density function

of g is

p

(

t

) =

1 for 0

≤

t

≤

1, 0 elsewhere.

Therefore, g has a uniform distribution, i.e.,

g

∼ U (

0, 1

)

_.

Remark: Let X be a random variable and p ( t ) the probability density function (pdf) of X. The cumulative distribution function (cdf) of X is

F

( η )

:

=

Prob

(

X

≤ η ) = Z _η

− ∞

p

(

t

)

dt.

(20)

Example of histogram equalized image

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

0 50 100 150 200 250

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

0 50 100 150 200 250

Histogram equalization of 400 × 600 image: (top) before; (bottom) after; and the corresponding histograms

Matlab commands: imhist(A) , histeq(A)

(21)

Histogram-equalized images

3.3 Histogram Processing 141

FIGURE 3.20 Left column: Images from Fig. 3.16. Center column: Corresponding histogram-equalized images. Right column: histograms of the images in the center column (compare with the histograms in Fig. 3.16).

www.EBooksWorld.ir

Histogram-equalized images and the corresponding normalized histograms

(22)

Histogram-equalized images (cont’d)

3.3 Histogram Processing 141

FIGURE 3.20 Left column: Images from Fig. 3.16. Center column: Corresponding histogram-equalized images. Right column: histograms of the images in the center column (compare with the histograms in Fig. 3.16).

www.EBooksWorld.ir

Histogram-equalized images and the corresponding normalized histograms

(23)

Spatial filter (空空空間間間濾濾濾波波波)

Spatial filtering modifies an image by replacing the value of each pixel by a function of the values of the pixel and its neighbors.

(discrete!) If the operation performed on the image pixels is linear, then the filter is called a linear spatial filter. Otherwise, the filter is a nonlinear spatial filter.

A linear spatial filter performs a sum-of-products operation between an image f and a filter kernel w.

The kernel is an array whose size defines the neighborhood of operation, and whose coefficients (entries) determine the nature of the filter.

Other terms used to refer to a spatial filter kernel are

mask, template,

and

window.

We use the term

“filter kernel”

or simply

“kernel.”

(24)

Filter kernel of a linear spatial filter

Linear spatial filtering using a 3 × 3 kernel

(25)

Linear spatial filtering

3 × 3 kernel: at any point ( x, y ) in the image f , the response g ( x, y ) of the filter is the sum of products of the kernel entries (coefficients) and the image pixels:

g

(

x, y

) =

w

(−

1,

−

1

)

f

(

x

−

1, y

−

1

) +

w

(−

1, 0

)

f

(

x

−

1, y

) + · · · +

w

(

0, 0

)

f

(

x, y

)

| {z }

+ · · · +

w

(

1, 1

)

f

(

x

+

1, y

+

1

)

.

As x and y are varied, the

center (origin) of the kernel

moves from pixel to pixel, generating the filtered image g.

m × n kernel: Assume that m = 2a + 1 and n = 2b + 1. Then

g

(

x, y

) =

∑ a i=−a

∑ b j=−b

w

(

i, j

)

f

(

x

+

i, y

+

j

)

.

This is talking about the

(discrete) spatial correlation (空間相關性).

The mechanics of

(discrete) spatial convolution (空間卷積,

∗

or

~

)

are the same, except that the kernel is rotated by 180 ^◦ counterclockwise.

(26)

Convolution of two functions

Let f and g be two integrable real-valued functions with compact supports in R. Then the convolution of f and g is defined as a function in variable t,

(

f

∗

g

)(

t

)

:

= Z _∞

− ∞

f

( τ )

g

(

t

− τ )

dτ, t

∈ _R.

The operation can be described as

a weighted average of the input f according to the weighting (or kernel) g at each point in time t.

Let f and g be two integrable real-valued functions with compact supports in R ^d . Then the convolution of f and g is defined as a function in variable x,

(

f

∗

g

)( x )

:

= Z

R

^df

( y )

g

( x − y )

dy,

x ∈ R ^d

.

Commutativity: f ∗ g = g ∗ f

Associativity: ( f ∗ g ) ∗ h = f ∗ ( g ∗ h ) Distributivity: f ∗ ( g + h ) = ( f ∗ g ) + ( f ∗ h )

(27)

Commutativity: f ∗ g = g ∗ f

∵ f and g are two integrable functions with compact supports in R.

∴ ∃ L > 0 such that f ( _t ) = ₀ = _g ( _t ) _{for t} 6∈ [− L, L ] _.

∴ ( f ∗ g )( t ) : = Z _∞

− ∞ f ( τ ) g ( t − τ ) dτ = Z _L

−L f ( τ ) g ( t − τ ) dτ, ∀ t ∈ _R Let η = −( τ − t ) . Then τ = t − η and dη = − dτ, and we have Z L

−L f ( τ ) g ( t − τ ) dτ = Z t−L

t+L f ( t − η ) g ( η )(− dη ) = Z t+L

t−L f ( t − η ) g ( η ) dη.

If t ≥ 0, then Z _t+L

t−L f ( t − η ) g ( η ) dη = Z _L

−L f ( t − η ) g ( η ) dη = ( g ∗ f )( t ) . If t < 0, then

Z t+L

t−L f ( t − η ) g ( η ) dη = Z L

−L f ( t − η ) g ( η ) dη = ( g ∗ f )( t ) .

∴ (

f

∗

g

)(

t

) = (

g

∗

f

)(

t

)

,

∀

t

∈ R

(28)

Convolution of two 1-D functions

(

f

∗

g

)(

t

)

:

= Z _∞

− ∞

f

( τ )

g

(

t

− τ )

dτ, t

∈ R

Wikipedia: https://en.wikipedia.org/wiki/Convolution

(29)

Discrete convolution: 1-D

The discrete convolution of input (signal) f and kernel g is defined by

(

f

∗

g

)(

t

)

:

=

∑ ∞ τ=−∞,τ∈Z

f

( τ )

g

(

t

− τ )

, t

∈ _Z.

When f and g have finite supports, a finite summation may be used.

f and g can be viewed as piecewise constant functions in each unit integer interval.

(30)

Correlation vs. convolution: 1-D example

(31)

1-D discrete full convolution

Let u = [ u 1 , · · · , u n ] ^> ∈ R ⁿ and v = [ v 1 , · · · , v m ] ^> ∈ R ^m . The convolution of u and v is defined as

u ∗ v : =







u

1

v

1

u

1

v

2 +

u

2

v

1

u

1

v

3 +

u

2

v

2 +

u

3

v

1

...

u

n−2

v

m +

u

n−1

v

m−1 +

u

n

v

m−2

u

n−1

v

m +

u

n

v

m−1

u

n

v

m







∈ R ^m+n−1

.

(32)

Convolution: 1-D example in M ATLAB

u=[1 1 1];

v=[1 1 0 0 0 1 1];

—————- w1=conv(u,v) w1=

1 2 2 1 0 1 2 2 1

—————-

w2=conv(u,v,’same’) w2=

1 0 1

—————-

w3=conv(u,v,’valid’) w3=

1 × 0 empty · · ·

—————-

w4=conv(v,u,’valid’) w4=

2 1 0 1 2

(33)

Correlation vs. convolution: 2-D example

(34)

2-D discrete convolution: conv2(f,K,’valid’)

no padding, stride 1

In M ATLAB : conv2(f,K,’valid’)

Full convolution: conv2(f,K) =⇒ ( 7 + 3 − 1 ) × ( 7 + 3 − 1 ) matrix!

(35)

Stride (步步步長長長)

(correction: 34 34)

Convolution of a 6 × 6 matrix and a 3 × 3 filter with stride 3, no padding

=⇒ 2 × 2 matrix

(36)

Padding (填填填補補補)

Convolution of a 6 × 6 matrix with zero-padding 1 and a 3 × 3 filter with stride 1

In M ATLAB : conv2(f,K,’same’)

(37)

A Matlab file for convolution and correlation

clear all clc

m=5;%image size

w=3;%window size of convolution I=reshape(1:mˆ2,m,m)

K=reshape(1:wˆ2,w,w)

%convolution

conv2 output=conv2(I,K,’valid’)

%manual implementation C=zeros(m-w+1,m-w+1);

for i=1:m-w+1 for j=1:m-w+1

C(i,j)=sum(sum(I(i:i+w-1,j:j+w-1). ∗ rot90(K,2)));

end end C

(38)

A Matlab file for convolution and correlation (cont’d)

%correlation

corr output=filter2(K,I,’valid’)

% manual implementation D=zeros(m-w+1,m-w+1);

for i=1:m-w+1 for j=1:m-w+1

D(i,j)=sum(sum(I(i:i+w-1,j:j+w-1). ∗ K));

end end D

—————————————-

% function ‘imfilter’ is provided in the MATLAB toolbox imfilter conv output=imfilter(I,K,’conv’,’same’)

MA3111: Mathematical Image Processing Intensity Transformations & Spatial Filtering