Using numpy and scipy on images, applying filters, etc

misc_image_manipulations_and_filters.py

import sys
from time import time
from PIL import Image
import numpy as np
from scipy.ndimage.filters import generic_filter, gaussian_filter


def trace(fn):
   """A decorator to time your functions"""
   def trace_func(*args, **kwargs):
       print fn.__name__ + '...',
       sys.stdout.flush()
       beg = time()
       ret = fn(*args, **kwargs)
       tot = time() - beg
       print '%.3f' % tot
       return ret
   return trace_func

def matlab_style_gauss2D(shape=(3,3),sigma=0.5):
    """
    2D gaussian mask - should give the same result as MATLAB's
    fspecial('gaussian',[shape],[sigma])
    http://stackoverflow.com/questions/17190649/how-to-obtain-a-gaussian-filter-in-python
    """
    m,n = [(ss-1.)/2. for ss in shape]
    y,x = np.ogrid[-m:m+1,-n:n+1]
    h = np.exp( -(x*x + y*y) / (2.*sigma*sigma) )
    h[ h < np.finfo(h.dtype).eps*h.max() ] = 0
    sumh = h.sum()
    if sumh != 0:
        h /= sumh
    return h
   
@trace
def my_image_filter2(img, filter):
    img = Image.open(img)
    img.convert('L')
    array = np.array(img)
    flt_filter = filter.flatten()
    
    def fnc(buffer):
        return (flt_filter * buffer).sum()
    
    filtered_array = generic_filter(array, fnc, size=filter.shape, mode='constant', cval=0)
    return Image.fromarray(filtered_array,'L')

@trace
def my_image_filter_gaussian(img):
    img = Image.open(img)
    img.convert('L')
    array = np.array(img)
    
    filtered_array = gaussian_filter(array, sigma=3, mode='constant', cval=0)
    return Image.fromarray(filtered_array,'L')

@trace
def median_filter(img, shape):
    img = Image.open(img)
    img.convert('L')
    array = np.array(img)
    filtered_array = np.zeros(array.shape, np.uint8)
    
    x_win_width = shape[0]/2
    y_win_width = shape[1]/2
    
    for y in range(array.shape[1]):
        #import pdb;pdb.set_trace()
        for x in range(array.shape[0]):
            #print x,y
            #print "looking at neigbhorhood",x-x_win_width,x+x_win_width+1, y-y_win_width,y+y_win_width+1
            neighborhood = array[x-x_win_width:x+x_win_width+1, y-y_win_width:y+y_win_width+1]
            if neighborhood.shape == shape:
                #import pdb;pdb.set_trace()
                vals = neighborhood.flatten()
                vals.sort()
                filtered_array[x,y] = vals[len(vals)/2]
            else:
                pass #Leave as 0 for borders
    return Image.fromarray(filtered_array,'L')

def normalize(array, min_range, max_range):
    array = (min_range + ((array-min_range)*(max_range-min_range)))/(array.max() - array.min())
    return array
    

@trace
def gradient(img):
    img = Image.open(img)
    img.convert('L')
    array = np.array(img)
    gx, gy = np.gradient(array)
    mag = np.sqrt(gx**2 + gy**2)
    #stretch over 0-255
    mag = normalize(mag, 0, 255)
    direction = np.nan_to_num(np.arctan(gx/gy))
    print direction
    direction = normalize(direction,0,255)
    print np.histogram(mag, bins=255)[0]
    print np.histogram(direction, bins=255)
    return Image.fromarray(mag,'L')
    
@trace
def edges(img):
    from skimage import feature
    img = Image.open(img)
    img.convert('L')
    array = np.array(img)
    print repr(array)
    out = np.uint8(feature.canny(array, sigma=3, ) * 255)
    print repr(out)
    return Image.fromarray(out,mode='L')
    
'''
@trace
def edges(img):
    from skimage import feature
    from skimage import io
    img = io.imread('Q_3.jpg')
    print repr(img)
    #img.convert('L')
    #array = np.array(img)
    #print repr(array)
    out = feature.canny(img, sigma=1, )
    print repr(out)
    return Image.fromarray(out,'L')    
'''

@trace
def subtract_images(img1,img2):
    img1 = Image.open(img1)
    img1=img1.convert('L')
    img2 = Image.open(img2)
    img2 = img2.convert('L')
    array1 = np.array(img1)
    array2 = np.array(img2)
    result = array1 - array2
    print result
    #why no negative values, why doesn't normalize work?
    result = normalize(result, 0, 255)
    print result
    return Image.fromarray(result,'L')
            
@trace
def my_image_filter(img, filter):
    print "Applying filter", filter
    img = Image.open(img)
    img.convert('L')
    array = np.array(img)
    filtered_array = np.zeros(array.shape, np.uint8)
    
    x_win_width = filter.shape[0]/2
    y_win_width = filter.shape[1]/2
    
    for y in range(array.shape[1]):
        #import pdb;pdb.set_trace()
        for x in range(array.shape[0]):
            #print x,y
            #print "looking at neigbhorhood",x-x_win_width,x+x_win_width+1, y-y_win_width,y+y_win_width+1
            neighborhood = array[x-x_win_width:x+x_win_width+1, y-y_win_width:y+y_win_width+1]
            if neighborhood.shape == filter.shape:
                #import pdb;pdb.set_trace()
                filtered_array[x,y] = np.sum(filter * neighborhood)
            else:
                pass #Leave as 0 for borders
    return Image.fromarray(filtered_array,'L')


if __name__ == '__main__':
    edges('Q_3.jpg').save('Q_3_edges.jpg')
    #subtract_images('walk_1.jpg','walk_2.jpg').save('walk_new.jpg')
    #gradient('Q_3.jpg').save('Q_3_gradient.jpg')
    sys.exit()
    #my_image_filter_gaussian('House1.jpg').save('House1_scipyguass.jpg')
    median_filter('Noisyimage1.jpg', (5,5)).save('Noisyimage1_median.jpg')
    median_filter('Noisyimage2.jpg', (5,5)).save('Noisyimage2_median.jpg')        
    my_image_filter('Noisyimage1.jpg', np.ones((5,5)) * (1/25.)).save('Noisyimage1_avg.jpg')
    my_image_filter('Noisyimage2.jpg', np.ones((5,5)) * (1/25.)).save('Noisyimage2_avg.jpg')
    sys.exit()
    filters = dict(
        avg_kernal_3_3 = np.ones((3,3)) * (1/9.),
        avg_kernal_5_5 = np.ones((5,5)) * (1/25.),
        sobel_3_3_x = np.array(
        [
        [-1,0,1],
        [-2,0,2],
        [-1,0,1],
        ]),
        sobel_3_3_y = np.array(
        [
        [-1,-2,-1],
        [0,0,0],
        [1,2,1],
        ]),
        prewitt_3_3_x = np.array(
        [
        [-1,0,1],
        [-1,0,1],
        [-1,0,1],
        ]),
        prewitt_3_3_y = np.array(
        [
        [1,1,1],
        [0,0,0],
        [-1,-1,-1],
        ]),        
        guass_s1 = matlab_style_gauss2D(shape=(3,3), sigma=1),
        guass_s2 = matlab_style_gauss2D(shape=(5,5), sigma=2),
        guass_s3 = matlab_style_gauss2D(shape=(7,7), sigma=3),
    )
    for fname, filter in sorted(filters.items()):
        my_image_filter('House1.jpg', filter).save('House1_'+fname+'.jpg')
        my_image_filter('House2.jpg', filter).save('House2_'+fname+'.jpg')