Programming for Engineers in Python

Recitation 12
Image Processing
 Plan: Image Processing with numpy
 Binary segmentation

 Image gradient

 Image brightening

 Morphological operators
 Erosion
 Dilation
 Smoothing

 Denoising

 Ternary segmentation

2
 A 2D table of values
(pixels), each in
Grayscale Image 0..255:
• 0 = Black
• 255 = White

105 114 116 116 160 121 97 90 124 119 188 144 112 116 78 32 19 61 40

3
 Image processing – basic functions
 Reading an image from disk:
from scipy import misc
im = misc.imread('C:/Koala.jpg')
 Creating an “empty” image matrix:
im = numpy.zeros( (height,width), dtype=numpy.uint8 )

 Important: each pixel is in the range 0-255 (type np.uint8)

o Numerical operations cause overflow, e.g.:
numpy.uint8(200) + numpy.uint8(100) = 44 = 300 mod 256
o Therefore, before doing numerical operations on image pixels, convert the
whole image or a specific pixel to int 32-bit using numpy.int_ :
a = numpy.uint8(200) ; b = numpy.uint8(100)
numpy.int_(a) + numpy.int_(b) = 300
4
 Image processing – more functions
>>> A = numpy.array([[1,7,2],[6,9,3]])
>>> A
array([[1, 7, 2],
[6, 9, 3]])
 Find the maximal / minimal value:
numpy.min(A) → 1
numpy.max(A) → 9
 Find the average / median value:
numpy.mean(A) → 4.66
numpy.median(A) → 4.5

5
 Image processing – more functions
 We will need the following numpy function later:
numpy.concatenate( (A,B), axis=0)
array([[1, 2, 3],
>>> A = numpy.array([[1,2,3],[4,5,6]]) [4, 5, 6]])

>>> B = numpy.array([[7,8,9],[1,1,1]]) array([[7, 8, 9], array([[1, 2, 3],

[1, 1, 1]])
[4, 5, 6],
>>> numpy.concatenate((A,B), axis=0) [7, 8, 9],
[1, 1, 1]])
>>> C = numpy.array([[50,50]]) array([[50, 50]])
array([[50],
>>> C = C.T [50]])
array([[ 1, 2, 3, 50],
>>> numpy.concatenate((A,C), axis=1) [ 4, 5, 6, 50]])

6
 Binary Segmentation
 Goal: grayscale image  black and white image
 pixel > threshold  change to white (255)
 pixel <= threshold  change to black (0)

 Motivation: for example, save space

Threshold = 50 Threshold = 128 Threshold = 200

7
 Binary Segmentation: Code
>>> import matplotlib.pyplot as plt
>>> threshold = 128
>>> from scipy import misc
>>> binIm = (im > threshold) * 255
>>> im = misc.imread('C:/Koala.jpg')
>>> plt.figure()
>>> type(im)
>>> plt.imshow(im, cmap=plt.cm.gray)
<type 'numpy.ndarray'>
>>> plt.figure()
>>> im.shape, im.dtype
>>> plt.imshow(binIm, cmap=plt.cm.gray)
((598, 848), dtype('uint8'))
>>> plt.show()
>>> np.min(im), np.max(im)
(0, 255)

8
 Image gradient
 Given an image we may want to locate points of major change in
intensity.
 In the gradient image, every pixel is the difference between the
corresponding pixel in the original image and a neighboring pixel.
 We can compute the gradient wrt each of the image dimensions.

9
 Image vertical gradient: Code
 The idea: shift up and subtract images

zeroRow = np.zeros((1,im.shape[1]), dtype=np.uint8)

upShiftedIm = np.concatenate((im[1:], zeroRow), axis=0)
diff = np.int_(upShiftedIm) - np.int_(im)
imDy = np.abs(diff) array([[1, 2, 3],
[4, 5, 6],
plt.imshow(imDy, cmap=plt.cm.gray) [7, 8, 9]])
plt.show() =>
array([[4, 5, 6],
[7, 8, 9],
 Can also be computed directly using for loops ! [0, 0, 0]])

10
 Brightening an image
 We can hardly see the gradient image.
o Solution: let’s brighten it.

>>> brImDy = np.minimum(np.int_(imDy) * 5, 255)

11
 Image processing
 Many more operations can be performed on images, in particular:
o Geometric transformations: shifting, rotation, flipping
o Filtering: replace a pixel by a function of the neighboring pixels
 The neighborhood can be defined by a structuring element, e.g. square

 A good survey with numpy/scipy examples can be found in:

http://scipy-lectures.github.io/advanced/image_processing/

12
 Morph by Neighbors
 All the operators we’ll see have the same basic idea:
Change a pixel according to its
neighborhood
The ‘neighborhood’ of a pixel in
the input image

nx,ny=3 nx,ny=1 nx,ny=2

13
 http://luispedro.org/software/pymorph

Morphological Operators
 Morphology is a technique for the analysis and processing
of geometrical structures
 Erosion:
 Take binary image ; produce more black, less white.
 Edges of white areas (“foreground”) erode away.
 So black areas (“background”) grow.

14
 Morphological Operators
 Dilation – inverse of erosion.
 Boundaries of black areas erode away.
 So white areas grow.

15
 Morph – Common Code
The specific
Framework: morph. operator

def morphological(im, operator=np.min, nx=5, ny=5):

newIm = np.zeros(im.shape)
for x in range(im.shape[0]): current pixel
location
for y in range(im.shape[1]):
nlst = neighbours(im, x, y, nx, ny)
newIm[x,y] = operator(nlst)
return newIm
Size of
neighborhood
16
 Get Neighbors – Code
def neighbours(im, x, y, nx=1, ny=1):
return im[max(x-nx, 0) : min(x+nx+1, im.shape[0]), \
max(y-ny, 0) : min(y+ny+1, im.shape[1])]

Default Neighborhood
17
 Erosion and dilation using morph
def erosion(im, nx=5, ny=5):
return morphological(im, np.min, nx, ny)

• Black neighbor  change pixel to black

def dilation(im, nx=5, ny=5):

return morphological(im, np.max, nx, ny)

• White neighbor  change pixel to white

18
 Morphological Operators - Run

erosion

dilation

19
 Detect Edges with Dilation
 Idea:
 Thicken white edges a bit, using dilation
 Compute difference from original image

Dilated Original Difference

im1 = misc.imread('square1.bmp')
im2 = dilation(im1, 1, 1)
imDiff = im2 – im1
20
 Denoising
 We want to “clean” these pictures from noise
 Suggestions?

21
 Denoising - Mean
Take 1: Smoothing
Use the morphological framework with mean (average):
def denoise_mean(im, nx=5, ny=5):
return morphological(im, np.mean, nx, ny)

22
 Denoising - Median
Take 2: Median
def denoise_median(im, nx=5, ny=5):
return morphological(im, np.median, nx, ny)

23
 Denoising - Median
 Can we do better?
 We did not exploit what we know about the ‘noise’
 The ‘noise’ can only be
 Almost white
 Almost black

 So change only very bright or very dark pixels!

 What is very bright?
 Pixel>W
 What is very dark?
 Pixel<B
24
 Denoising - operation
 Pixels not too bright or dark remain with no change.
 Otherwise, examine the pixel’s neighborhood nlst.
 too dark:
 Get the median of nlst, but ignore ‘too dark’ pixels
 too bright:
 Get the median of nlst, but ignore ‘too bright’ pixels

25
 Denoising – Bounded Median
def my_median(nbrs, min_val=0, max_val=255): def fix_white(nbrs, W=200):
goodNbrs = \ center = (nbrs.shape[0]/2, nbrs.shape[1]/2)
nbrs[(nbrs >= min_val) & (nbrs <= max_val)] if nbrs[center] > W: # too bright
if goodNbrs.shape[0] > 0: return my_median(nbrs, max_val=W)
return np.median(goodNbrs) else:
else: return nbrs[center]
center = (nbrs.shape[0]/2, nbrs.shape[1]/2)
return nbrs[center]
def fix_black(nbrs, B=50):
center = (nbrs.shape[0]/2, nbrs.shape[1]/2)
if nbrs[center] < B: # too dark
return my_median(nbrs, min_val=B)
else:
return nbrs[center]

26
 Denoising – Bounded Median
def denoise_bounded_median(im, nx=3, ny=3):
num_repeat = 3
for i in range(num_repeat):
im = morphological(im, fix_white, nx, ny)
im = morphological(im, fix_black, nx, ny)
return im

27
 Denoising – Bounded Median
Iteration
number
1 2

3 4

28
 Ternary segmentation


29
 Example

Before After
30
 Ternary segmentation
 Example: If the values of the image pixels are 1 to 90 then:
first threshold = 30; second threshold = 60
 Possible algorithm:
 Sort all n*m pixel values in the image into sorted_pxls
(one-dimensional list)
 Set first threshold = sorted_pxls[n*m/3]
 Set second threshold = sorted_pxls[2*n*m/3]
 Complexity: O(nm*log(nm))

 Let’s try to solve in O(nm)

31
 Ternary segmentation
 We will break the problem into smaller problems:
1. Create a histogram.
2. Create a cumulative histogram.
3. Find 2 thresholds from cumulative histogram.
4. Transform the image using the thresholds.

 Step by step…
 Modularize the solution

32
 Image histogram
 pic.png contains a 20x20 pixels image
 Each pixel is a gray level in 0..255

 Q1: Implement the function build_hist():

 Input: grayscale image.
 Output: histogram of gray levels in image.

 So our histogram will be a list with 256 values:

 Index = gray level
 Value = #pixels with that gray level
33
 Image Histogram - Code
def build_hist(im):
hist = [0]*256
for x in range(im.shape[0]):
for y in range(im.shape[1]):
gray_level = im[x,y]
hist[gray_level] += 1
return hist

34
 Cumulative histogram
 Q2: Implement the function cum_sum(lst):
 Input: list of numbers
 Output: list of the cumulative sums of lst

12 43 17 22

12 55 72 94

35
 Cumulative sum – Code
def cum_sum(numbers):
res = []
for i in range(len(numbers)):
prevSum = res[i-1] if i > 0 else 0
res.append(prevSum + numbers[i])
return res

36
 Ternary segmentation
 Q3: Given an image, find two thresholds th1, th2 such that:

37
 Finding the thresholds
def find_thrds(im):
hist = build_hist(im)
c_lst = cum_sum(hist)
i=0
while float(c_lst[i]) / im.size < 1./3:
i += 1
th1 = i
while float(c_lst[i]) / im.size < 2./3:
i += 1
th2 = i
return (th1, th2)
38
 Ternary segmentation
 Q4: Implement the ternary segmentation function itself:
o Input: a grayscale image
o Output: a new image where every pixel has one of 3 gray levels:
0, 128 or 255

def ternary_sementation(im):
th1, th2 = find_thrds(im)
newIm = np.zeros(im.shape)
newIm = newIm + ((im > th1) & (im <= th2)) * 128
newIm = newIm + (im > th2) * 255
return newIm

39