CS1001.py

Extended Introduction to Computer Science with Python, Tel-Aviv University, Spring 2013

Recitation 13 - 10-13.6.2013

Last update: 10.6.2013

General image processing

This module is available here:


In [5]:
from IPMatrix import Matrix

These functions are used to convert bmp/jpg files to the matrix file format, bitmap:


In [7]:
import numpy
from PIL import Image

def bitmap2image(file):
    ''' convert format .bitmap image to .jpg/.bmp '''
    im = Matrix.load(file)
    n,m = im.dim()
    ar = numpy.zeros((n,m)) #an nXm array of zeros
    for i in range(n):
        for j in range(m):
            ar[i,j] = im[i,j]
    image = Image.fromarray(numpy.uint8(ar))
    image.save(file.split(".")[0] + ".bmp")

def image2bitmap(file):
    ''' convert .jpg/.bmp file to format .bitmap
        which can be loaded into class Matrix using Matrix.load() ''' 
    image = Image.open(file).convert("L") # converts to 8 bit black and white image
    im = numpy.asarray(image)             # stores image in an array
    im.setflags(write=1)                  # makes array writeable
    n,m = im.shape[0], im.shape[1]        # image dimensions

    new_file_name = file.split(".")[0] + ".bitmap"
    new_file = open(new_file_name, 'w')

    print(n,m, file=new_file)
    for i in range(n):
        for j in range(m):
            print(int(im[i,j]), end=" ", file = new_file)
        print("", file=new_file)#newline
    new_file.close()

Lets load an image and start working with it:


In [8]:
im = Matrix.load("eiffel.bitmap")
im


Out[8]:
<IPMatrix.Matrix at 0x86c8e48>

In [9]:
def mult(im,k):
    return im * k

In [10]:
mult(im, 2)


Out[10]:
<IPMatrix.Matrix at 0x86b78d0>

In [8]:
def add(im, k):
    n,m = im.dim()
    im2 = Matrix(n,m)
    for i in range(n):
        for j in range(m):
            im2[i,j] = (im[i,j] + k) % 256
    return im2

In [9]:
add(im,50)


Out[9]:
<IPMatrix.Matrix at 0x7b86320>

In [10]:
add(im,-50)


Out[10]:
<IPMatrix.Matrix at 0x7b79780>

High-order functions!


In [11]:
def pixelwise_operation(im, op):
    n,m = im.dim()
    im2 = Matrix(n,m)
    for i in range(n):
        for j in range(m):
            im2[i,j] = op(im, i, j)
    return im2

In [12]:
def add(im, k):
    op = lambda im,i,j: im[i,j] + k
    return pixelwise_operation(im, op)

In [13]:
add(im, 100)


Out[13]:
<IPMatrix.Matrix at 0x8832828>

In [17]:
secret(im) # SHHHHH, the code is below


Out[17]:
<IPMatrix.Matrix at 0x9960470>

In [18]:
def secret(im):
    n,m = im.dim()
    im2 = Matrix(n,m)
    for i in range(n):
        im2.rows[i] = sorted(im.rows[i])
    return im2

In [19]:
def negate(im):
    op = lambda im,i,j: 255 - im[i,j]
    return pixelwise_operation(im, op)

In [20]:
negate(im)


Out[20]:
<IPMatrix.Matrix at 0x995d898>

In [21]:
def shift(im, v_shift, h_shift):
    op = lambda im,i,j: im[(i - v_shift) % im.dim()[0], (j - h_shift) % im.dim()[1]]
    return pixelwise_operation(im, op)

In [22]:
shift(im, 100, 0)


Out[22]:
<IPMatrix.Matrix at 0x995d550>

In [23]:
shift(im, 0, 100)


Out[23]:
<IPMatrix.Matrix at 0x995d828>

In [24]:
shift(im, 50, 50)


Out[24]:
<IPMatrix.Matrix at 0x995da58>

In [25]:
def upside_down(im):
    op = lambda im,i,j: im[im.dim()[0] - i - 1, j]
    return pixelwise_operation(im, op)

In [26]:
upside_down(im)


Out[26]:
<IPMatrix.Matrix at 0x99740f0>

In [27]:
def concatenate(im1, im2):
    n1,m1 = im1.dim()
    n2,m2 = im2.dim()
    n = max(n1,n2)
    im3 = Matrix(n,m1+m2)
    im3[:n1,:m1] = im1
    im3[:n2, m1:m1+m2] = im2
    return im3

In [28]:
concatenate(im,im)


Out[28]:
<IPMatrix.Matrix at 0xa8b0358>

In [29]:
def strech(im, ratio=1):
    n,m = im.dim()
    im2 = Matrix(n * ratio, m)
    for i in range(n):
        for j in range(m):
            for k in range(ratio):
                im2[ratio * i + k , j] = im[i,j]
    return im2

In [30]:
strech(im,3)


Out[30]:
<IPMatrix.Matrix at 0xa8bc0f0>

This strech just multiplies every pixel.

A better solution would be to insert the average between adjacent pixels.

The same goes for squeeze.

Denoising

In class we saw two types o noises:

  • Gaussian noise (also called white noise) - reduced well by local means
  • Salt and pepper - reduced well by local medians

The basic scheme for denoising is that we change the value of each pixel wrt to the value of its neighbours using some function. The neighbors of a pixel is pixels in a square around the pixel.

methodprosconsworks well with
local meansfast
good on smooth areas
blurs edges
affected by outliers
Gaussian noise
local medianspreserves edges
not affected by outliers
slow (median calculation)
eliminates fine details (contours)
salt & pepper noise

Examples

This image has a XxX white spot:


In [5]:
X = 5
im[50:50 + X,50:50 + X] = Matrix(X, X, 255)
im


Out[5]:
<IPMatrix.Matrix at 0x85a18d0>

Which denoising method would you use to denoise it, with what neighborhood size?

To minimize the effect of the denoising on the rest of the image, we will use the local medians method.

We start by implementing a function that returns the neighborhood :


In [13]:
def items(mat):
  lst = []
  n,m = mat.dim()
  for i in range(n):
    lst = lst + [mat[i,j] for j in range(m)]
  return lst

def median(lst):
    return sorted(lst)[len(lst) // 2]

def local_medians(A, k=1):
  n,m = A.dim()
  res = A.copy()
  for i in range(k,n-k):
    for j in range(k,m-k):
      neighborhood = items(A[i-k:i+k+1,j-k:j+k+1])
      res[i,j] = median(neighborhood)
  return res

In [8]:
local_medians(im)


Out[8]:
<IPMatrix.Matrix at 0x90c34a8>

In [23]:
local_medians(im, k=3)


Out[23]:
<IPMatrix.Matrix at 0x90d9358>

In [15]:
local_medians(im, k=5)


Out[15]:
<IPMatrix.Matrix at 0x90d9128>

When using a small value for k the spot is not cleaned, when using a large value the image is blurred.

Solution - bounded local medians

When darkness/brightness is extreme (close to 0 or 255) we can do better:

  • fix only extreme values
  • neighborhoods only include non-extreme values

In [11]:
def local_medians_bounded(A, k=1, black=10, white=250):
  n,m = A.dim()
  res1 = A.copy()
  #fix white
  for i in range(k,n-k):
    for j in range(k,m-k):
      neighborhood = items(A[i-k:i+k+1,j-k:j+k+1])
      res1[i,j] = fix_white(neighborhood, white)
  #fix black        
  res2 = res1.copy()
  for i in range(k,n-k):
    for j in range(k,m-k):
      neighborhood = items(res1[i-k:i+k+1,j-k:j+k+1])
      res2[i,j] = fix_black(neighborhood, black)

  return res2

def fix_white(lst, white_thr):
  center = len(lst)//2
  if lst[center] > white_thr:
    new_lst = [i for i in lst if i <= white_thr]
    if new_lst != []:
      return median(new_lst)
  return lst[center]  #either center was not extreme
                      #or all its neighbors are extreme

def fix_black(lst, black_thr):
  center = len(lst)//2
  if lst[center] < black_thr:
    new_lst = [i for i in lst if i >= black_thr]
    if new_lst != []:
      return median(new_lst)
  return lst[center]  #either center was not extreme
                      #or all its neighbors are extreme

In [14]:
local_medians_bounded(im, k=3)


Out[14]:
<IPMatrix.Matrix at 0x95bd4e0>

Guess who?

First, how do we add Salt and Pepper noise?


In [15]:
def add_SP(mat, p=0.01):
  ''' Generates salt and pepper noise: Each pixel is "hit" indep.
      with prob. p. If hit, it has fifty fifty chance of becoming
      white or black. '''
  n,m = mat.dim()
  new = Matrix(n,m)
  for i in range (n):
      for j in range (m):
          rand = random.randint(0,20000)
          if rand < p * 20000:
              if rand % 2 == 1:
                  new[i,j] = 0
              elif rand % 2 == 0:
                  new[i,j] = 255
          else: 
              new[i,j] = mat[i,j]
  return new

Now we load a picture of a famous person, noise it and try to denoise it:


In [16]:
person = Matrix.load("person.bitmap")
person_noise = add_SP(person, p=0.7)
person_noise


Out[16]:
<IPMatrix.Matrix at 0x95ce9e8>

In [17]:
local_medians(person_noise)


Out[17]:
<IPMatrix.Matrix at 0x95bd668>

In [18]:
person_denoised = local_medians_bounded(person_noise)
person_denoised


Out[18]:
<IPMatrix.Matrix at 0x95c30f0>

In [19]:
local_medians_bounded(person_denoised)


Out[19]:
<IPMatrix.Matrix at 0x95bd898>

Here is the original:


In [22]:
person


Out[22]:
<IPMatrix.Matrix at 0x90d94e0>

Fin

This notebook is part of the Extended introduction to computer science course at Tel-Aviv University.

The notebook was written using Python 3.2 and IPython 0.13.1 with Pillow 2.0.0 and NumPy 1.7.1rc1.

The code is available at https://raw.github.com/yoavram/CS1001.py/master/recitation13.ipynb.

The notebook can be viewed online at http://nbviewer.ipython.org/urls/raw.github.com/yoavram/CS1001.py/master/recitation13.ipynb.

This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.