Add gabor filter to filters in image processing

Image_Processing/src/Filters/src/gabor filter/frequency.py

@@ -0,0 +1,123 @@
+# Metody biometryczne
+# Przemyslaw Pastuszka
+
+from PIL import Image, ImageDraw
+import utils
+import argparse
+import math
+
+
+def points_on_line(line, W):
+    im = Image.new("L", (W, 3 * W), 100)
+    draw = ImageDraw.Draw(im)
+    draw.line([(0, line(0) + W), (W, line(W) + W)], fill=10)
+    im_load = im.load()
+
+    points = []
+    for x in range(0, W):
+        for y in range(0, 3 * W):
+            if im_load[x, y] == 10:
+               points.append((x, y - W))
+
+    del draw
+    del im
+
+    dist = lambda x, y: (x - W / 2) ** 2 + (y - W / 2) ** 2
+
+    return sorted(points)[:W]
+
+def vec_and_step(tang, W):
+    (begin, end) = utils.get_line_ends(0, 0, W, tang)
+    (x_vec, y_vec) = (end[0] - begin[0], end[1] - begin[1])
+    length = math.hypot(x_vec, y_vec)
+    (x_norm, y_norm) = (x_vec / length, y_vec / length)
+    step = length / W
+
+    return (x_norm, y_norm, step)
+
+def block_frequency(i, j, W, angle, im_load):
+    tang = math.tan(angle)
+    ortho_tang = -1 / tang
+
+    (x_norm, y_norm, step) = vec_and_step(tang, W)
+    (x_corner, y_corner) = (0 if x_norm >= 0 else W, 0 if y_norm >= 0 else W)
+
+    grey_levels = []
+
+    for k in range(0, W):
+        line = lambda x: (x - x_norm * k * step - x_corner) * ortho_tang + y_norm * k * step + y_corner
+        points = points_on_line(line, W)
+        level = 0
+        for point in points:
+            level += im_load[point[0] + i * W, point[1] + j * W]
+        grey_levels.append(level)
+
+    treshold = 100
+    upward = False
+    last_level = 0
+    last_bottom = 0
+    count = 0.0
+    spaces = len(grey_levels)
+    for level in grey_levels:
+        if level < last_bottom:
+            last_bottom = level
+        if upward and level < last_level:
+            upward = False
+            if last_bottom + treshold < last_level:
+                count += 1
+                last_bottom = last_level
+        if level > last_level:
+            upward = True
+        last_level = level
+
+    return count / spaces if spaces > 0 else 0
+
+def freq(im, W, angles):
+    (x, y) = im.size
+    im_load = im.load()
+    freqs = [[0] for i in range (0, int(x / W))]
+
+    for i in range(1, int(x / W - 1)):
+        for j in range(1, int(y / W - 1)):
+            freq = block_frequency(i, j, W, angles[i][j], im_load)
+            freqs[i].append(freq)
+        freqs[i].append(0)
+
+    freqs[0] = freqs[-1] = [0 for i in range(0, y / W)]
+
+    return freqs
+
+def freq_img(im, W, angles):
+    (x, y) = im.size
+    freqs = freq(im, W, angles)
+    freq_img = im.copy()
+
+    for i in range(1, x / W - 1):
+        for j in range(1, y / W - 1):
+            box = (i * W, j * W, min(i * W + W, x), min(j * W + W, y))
+            freq_img.paste(freqs[i][j] * 255.0 * 1.2, box)
+
+    return freq_img
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Image frequency")
+    parser.add_argument("image", nargs=1, help = "Path to image")
+    parser.add_argument("block_size", nargs=1, help = "Block size")
+    parser.add_argument('--smooth', "-s", action='store_true', help = "Use Gauss for smoothing")
+    args = parser.parse_args()
+
+    im = Image.open(args.image[0])
+    im = im.convert("L")  # covert to grayscale
+    im.show()
+
+    W = int(args.block_size[0])
+
+    f = lambda x, y: 2 * x * y
+    g = lambda x, y: x ** 2 - y ** 2
+
+    angles = utils.calculate_angles(im, W, f, g)
+    if args.smooth:
+        angles = utils.smooth_angles(angles)
+
+    freq_img = freq_img(im, W, angles)
+    freq_img.show()

Image_Processing/src/Filters/src/gabor filter/frequest.py

@@ -0,0 +1,108 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Apr 22 02:51:53 2016
+
+@author: utkarsh
+"""
+
+
+
+# FREQEST - Estimate fingerprint ridge frequency within image block
+#
+# Function to estimate the fingerprint ridge frequency within a small block
+# of a fingerprint image.  This function is used by RIDGEFREQ
+#
+# Usage:
+#  freqim =  freqest(im, orientim, windsze, minWaveLength, maxWaveLength)
+#
+# Arguments:
+#         im       - Image block to be processed.
+#         orientim - Ridge orientation image of image block.
+#         windsze  - Window length used to identify peaks. This should be
+#                    an odd integer, say 3 or 5.
+#         minWaveLength,  maxWaveLength - Minimum and maximum ridge
+#                     wavelengths, in pixels, considered acceptable.
+# 
+# Returns:
+#         freqim    - An image block the same size as im with all values
+#                     set to the estimated ridge spatial frequency.  If a
+#                     ridge frequency cannot be found, or cannot be found
+#                     within the limits set by min and max Wavlength
+#                     freqim is set to zeros.
+#
+# Suggested parameters for a 500dpi fingerprint image
+#   freqim = freqest(im,orientim, 5, 5, 15);
+#
+# See also:  RIDGEFREQ, RIDGEORIENT, RIDGESEGMENT
+
+### REFERENCES
+
+# Peter Kovesi 
+# School of Computer Science & Software Engineering
+# The University of Western Australia
+# pk at csse uwa edu au
+# http://www.csse.uwa.edu.au/~pk
+
+
+import numpy as np
+import math
+import scipy.ndimage
+#import cv2
+def frequest(im,orientim,windsze,minWaveLength,maxWaveLength):
+    rows,cols = np.shape(im);
+
+    # Find mean orientation within the block. This is done by averaging the
+    # sines and cosines of the doubled angles before reconstructing the
+    # angle again.  This avoids wraparound problems at the origin.
+
+
+    cosorient = np.mean(np.cos(2*orientim));
+    sinorient = np.mean(np.sin(2*orientim));    
+    orient = math.atan2(sinorient,cosorient)/2;
+
+    # Rotate the image block so that the ridges are vertical    
+
+    #ROT_mat = cv2.getRotationMatrix2D((cols/2,rows/2),orient/np.pi*180 + 90,1)    
+    #rotim = cv2.warpAffine(im,ROT_mat,(cols,rows))
+    rotim = scipy.ndimage.rotate(im,orient/np.pi*180 + 90,axes=(1,0),reshape = False,order = 3,mode = 'nearest');
+
+    # Now crop the image so that the rotated image does not contain any
+    # invalid regions.  This prevents the projection down the columns
+    # from being mucked up.
+
+    cropsze = int(np.fix(rows/np.sqrt(2)));
+    offset = int(np.fix((rows-cropsze)/2));
+    rotim = rotim[offset:offset+cropsze][:,offset:offset+cropsze];
+
+    # Sum down the columns to get a projection of the grey values down
+    # the ridges.
+
+    proj = np.sum(rotim,axis = 0);
+    dilation = scipy.ndimage.grey_dilation(proj, windsze,structure=np.ones(windsze));
+
+    temp = np.abs(dilation - proj);
+
+    peak_thresh = 2;    
+
+    maxpts = (temp<peak_thresh) & (proj > np.mean(proj));
+    maxind = np.where(maxpts);
+
+    rows_maxind,cols_maxind = np.shape(maxind);
+
+    # Determine the spatial frequency of the ridges by divinding the
+    # distance between the 1st and last peaks by the (No of peaks-1). If no
+    # peaks are detected, or the wavelength is outside the allowed bounds,
+    # the frequency image is set to 0    
+
+    if(cols_maxind<2):
+        freqim = np.zeros(im.shape);
+    else:
+        NoOfPeaks = cols_maxind;
+        waveLength = (maxind[0][cols_maxind-1] - maxind[0][0])/(NoOfPeaks - 1);
+        if waveLength>=minWaveLength and waveLength<=maxWaveLength:
+            freqim = 1/np.double(waveLength) * np.ones(im.shape);
+        else:
+            freqim = np.zeros(im.shape);
+
+    return(freqim);
+

Image_Processing/src/Filters/src/gabor filter/gabor.py

@@ -0,0 +1,77 @@
+# Metody biometryczne
+# Przemyslaw Pastuszka
+
+from PIL import Image, ImageDraw
+import utils
+import argparse
+import math
+import frequency
+import os
+
+def gabor_kernel(W, angle, freq):
+    cos = math.cos(angle)
+    sin = math.sin(angle)
+
+    yangle = lambda x, y: x * cos + y * sin
+    xangle = lambda x, y: -x * sin + y * cos
+
+    xsigma = ysigma = 4
+
+    return utils.kernel_from_function(W, lambda x, y:
+        math.exp(-(
+            (xangle(x, y) ** 2) / (xsigma ** 2) +
+            (yangle(x, y) ** 2) / (ysigma ** 2)) / 2) *
+        math.cos(2 * math.pi * freq * xangle(x, y)))
+
+def gabor(im, W, angles):
+    (x, y) = im.size
+    im_load = im.load()
+
+    freqs = frequency.freq(im, W, angles)
+    print ("computing local ridge frequency done")
+
+    gauss = utils.gauss_kernel(3)
+    utils.apply_kernel(freqs, gauss)
+
+    for i in range(1, x / W - 1):
+        for j in range(1, y / W - 1):
+            kernel = gabor_kernel(W, angles[i][j], freqs[i][j])
+            for k in range(0, W):
+                for l in range(0, W):
+                    im_load[i * W + k, j * W + l] = utils.apply_kernel_at(
+                        lambda x, y: im_load[x, y],
+                        kernel,
+                        i * W + k,
+                        j * W + l)
+
+    return im
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Gabor filter applied")
+    parser.add_argument("-t","--image", nargs=1, help = "Path to image")
+    parser.add_argument("-l","--block_size", nargs=1, help = "Block size")
+    parser.add_argument("--save", action='store_true', help = "Save result image as src_image_enhanced.gif")
+    args = parser.parse_args()
+
+    im = Image.open(args.image[0])
+    im = im.convert("L")  # covert to grayscale
+    im.show()
+
+    W = int(args.block_size[0])
+
+    f = lambda x, y: 2 * x * y
+    g = lambda x, y: x ** 2 - y ** 2
+
+    angles = utils.calculate_angles(im, W, f, g)
+    print(angles)
+    print ("calculating orientation done")
+
+    angles = utils.smooth_angles(angles)
+    print("smoothing angles done")
+
+    result = gabor(im, W, angles)
+    result.show()
+
+    if args.save:
+        base_image_name = os.path.splitext(os.path.basename(args.image[0]))[0]
+        im.save(base_image_name + "_enhanced.gif", "GIF")

Image_Processing/src/Filters/src/gabor filter/gabor_enhancement.py

@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Apr 18 11:42:58 2016
+
+@author: utkarsh
+"""
+
+import numpy as np
+#import cv2
+#import numpy as np;
+import matplotlib.pylab as plt;
+import scipy.ndimage
+import sys
+import cv2
+
+from image_enhance import image_enhance
+
+
+if(len(sys.argv)<2):
+    print('loading sample image');
+    img_name = '01.jpg'
+    img = scipy.ndimage.imread('images/' + img_name);
+
+elif(len(sys.argv) >= 2):
+    img_name = sys.argv[1];
+    img = scipy.ndimage.imread(sys.argv[1]);
+
+if(len(img.shape)>2):
+    # img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
+    img = np.dot(img[...,:3], [0.299, 0.587, 0.114]);
+
+
+print(img.shape)
+
+rows,cols = np.shape(img);
+aspect_ratio = np.double(rows)/np.double(cols);
+
+new_rows = 350;             # randomly selected number
+new_cols = new_rows/aspect_ratio;
+
+#img = cv2.resize(img,(new_rows,new_cols));
+img = scipy.misc.imresize(img,(np.int(new_rows),np.int(new_cols)));
+
+enhanced_img = image_enhance(img);
+enhanced_img = 255*np.uint8(enhanced_img)
+kernel = np.ones((5,5),np.uint8)
+# closing = cv2.morphologyEx(enhanced_img, cv2.MORPH_OPEN, kernel)
+erosion = cv2.erode(enhanced_img,kernel,iterations = 1)
+
+cv2.imshow('output',enhanced_img)
+cv2.waitKey(0)
+cv2.imshow('outclosed',erosion)
+cv2.waitKey(0)
+
+
+if(1):
+    print('saving the image')
+    scipy.misc.imsave('../enhanced/' + img_name, enhanced_img)
+else:
+    plt.imshow(enhanced_img,cmap = 'Greys_r');