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
