<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Quiz 3a - Implementation of the Harris Corner Detection
Eugenio Pacceli
Renato Oliveira
Brayan Acevedo
Step1:  First off, Harris detector computes a squared matrix M comprised basically of derivatives of image pixels on both x and y axis.
Then it calculates the corner responde number by subtracting the k*trace(M) from the determinant of M
Step3:  Afterwards, Harris detector proceeds with non-maximal suppression assuming a certain threshold (0.1 in this case)
Step4:  Compute Harris for ten levels of threshold
Keypoints are shown as red stars
Step5:  Quiz 3b - Implementation of the SIFT Feature Detector
Step6:  Execute custom SIFT for every given image in the current path then show keypoints as circles in the correspoding images.
SIFT algorithm begins detecting points that are invariant to scale changes. This is achieved by building a scale-space function named L(x,y,sigma) that consists of a convolution of a gaussian function and the image for every pixel. Following Lowe's paper (LOWE,2004) to efficiently detect keypoint locations in space scale we have proposed using a space scale extrema in the difference of Gaussian function convolved with the image.
In this implementation, the functions calcGaussianPyramid and getGaussianKernel compute the octaves and build a scape space for the image.
Afterwards, the difference of two consecutive guassian spaces separated by a constant factor k is carried out. This is implemented in the function calcDoG
Now, a local extrema detection step is executed. As Lowe stated each sample point is compared to each of its neighbors in the current image and nine of neighbors in the scale above and below.
Only the largest or the smallest of the points is selected as a candidate keypoint extrema.
To improve matching and stability a detailed fit to nearby data for location, scale and ratio of principal curvatures is performed. The function calcExtrema executes this step. This eliminates points with low contrast or that poorly localized along an edge.
And finally, the selected keypoints are drawn on the image.
Step7:  SIFT Feature Detector - Opencv Implementation
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import cv2
import matplotlib.pyplot as plt
import os
import glob
import random as rnd
from scipy.ndimage import filters
from PIL import Image
from numpy import *
from pylab import *
from pandas import *
np.seterr(divide='ignore', invalid='ignore')
Explanation: Quiz 3a - Implementation of the Harris Corner Detection
Eugenio Pacceli
Renato Oliveira
Brayan Acevedo
End of explanation
#Compute the Algorithm Harris corner detecion for implementation in grayscale image
def compute_harris_points(img, sigma=3):
    #compute derivates in the image 
    imx = np.zeros(img.size)
    imy = np.zeros(img.size)
    
    imx = filters.gaussian_filter(img, (sigma,sigma), (0,1))
    imy = filters.gaussian_filter(img, (sigma,sigma), (1,0))
    
    # compute the products of derivatives at every pixel
    Sxx = filters.gaussian_filter(imx*imx,sigma)
    Sxy = filters.gaussian_filter(imx*imy,sigma)
    Syy = filters.gaussian_filter(imy*imy,sigma)
    
    # determinant and trace
    Mdet = Sxx*Syy - Sxy**2
    Mtr = Sxx + Syy
    harris = np.divide(Mdet, Mtr)
    harris[np.isposinf(harris)] = 0
    harris[np.isnan(harris)] = 0
    return harris
Explanation: First off, Harris detector computes a squared matrix M comprised basically of derivatives of image pixels on both x and y axis.
Then it calculates the corner responde number by subtracting the k*trace(M) from the determinant of M
End of explanation
def doHarrisNonMaxSupression(harrisim,min_dist=10,threshold=0.1):
    #Return corners from a Harris response image
    #min_dist is the minimum number of pixels separating
    #corners and image boundary. 
    global t 
    global dist
    dist=min_dist
    t=threshold
    #print(t)
    # find top corner candidates above a threshold
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # get coordinates of candidates
    coords = array(harrisim_t.nonzero()).T
    # ...and their values
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # sort candidates
    index = argsort(candidate_values)
    # store allowed point locations in array
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # select the best points taking min_distance into account
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
            (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
    return filtered_coords
def plot_harris_points(image,filtered_coords):
    # Plots corners found in image. 
    plt.figure(figsize=(20,12))
    gray()
    plt.imshow(image)
    plt.title('Harris corner detection, dist=%s and threshold=%s'%(dist,t))
    plt.plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*',color = 'r')
    plt.axis('off')
    plt.show()
Explanation: Afterwards, Harris detector proceeds with non-maximal suppression assuming a certain threshold (0.1 in this case)
End of explanation
im = Image.open('boat_images/img1.pgm')
plt.figure(figsize=(20,12))
gray()
plt.imshow(im, cmap = 'gray')
harrisim = compute_harris_points(im)
for i in range(1, 11,1):
        xx=i* 0.01
        j=10
        filtered_coords = doHarrisNonMaxSupression(harrisim,j,xx)
        plot_harris_points(im, filtered_coords)
Explanation: Compute Harris for ten levels of threshold
Keypoints are shown as red stars
End of explanation
def getGaussianKernel(sigma, kernelHeight=51, kernelWidth=51):
    assert(kernelHeight % 2 == 1 and kernelWidth % 2 == 1)
    yOffset = (kernelHeight - 1) / 2
    xOffset = (kernelWidth - 1) / 2
    kernel = np.ndarray((kernelHeight, kernelWidth), np.float64)
    for y in range(-yOffset, yOffset+1, 1):
        for x in range(-xOffset, xOffset+1, 1):
            kernel[y+yOffset][x+xOffset] = (1. / (2.*np.pi*sigma**2)) * np.exp(-(x**2 + y**2) / (2 * sigma**2))
    kernel /= kernel.sum()
    return kernel
def calcGaussianPyramid(org_img):
    img = org_img.copy()
    bluredImg = img.copy()
    sigma = 1.6
    octaveCount = 7
    sigmaCount = 4
    gp = np.ndarray(shape=(octaveCount,), dtype=np.ndarray)
    for o in range(0, octaveCount):
        gp[o] = np.ndarray(shape=(sigmaCount+1, img.shape[0], img.shape[1]), dtype=np.float64)
        gp[o][0] = bluredImg.copy()
        for s in range(1, sigmaCount + 1):
            #k = 2**(float(s)/float(sigmaCount))
            k = np.sqrt(2.0)**s
            kernel = getGaussianKernel(k*sigma)
            bluredImg = cv2.filter2D(img, -1, kernel)
            gp[o][s] = bluredImg.copy()
        if (o < octaveCount-1):
            img = downscale(img)
            bluredImg = downscale(bluredImg)
    return gp
def calcDifference(img0, img1, threshold = 0):
    assert(img0.shape == img1.shape)
    diffImg = np.ndarray(img0.shape, np.float64)
    for y in range(diffImg.shape[0]):
        for x in range(diffImg.shape[1]):
            difference = abs(img1[y][x] - img0[y][x])
            if difference > threshold:
                diffImg[y][x] = difference
            else:
                diffImg[y][x] = 0
    return diffImg
def calcDoG(gp):
 
    DoG = np.ndarray(shape=gp.shape, dtype=np.ndarray)
    for o in range(DoG.shape[0]):
        DoG[o] = np.ndarray(shape=(gp[o].shape[0]-1, gp[o].shape[1], gp[o].shape[2]), dtype=np.float64)
        for s in range(DoG[o].shape[0]):
            DoG[o][s] = calcDifference(gp[o][s], gp[o][s+1])
    return DoG
def getNeighbourhood(octave, s, y, x, radius=1):
    neighbourhood = octave[s-radius:s+radius+1, y-radius:y+radius+1, x-radius:x+radius+1]
    return neighbourhood
def calcExtrema(DoG, threshold=0.3, radius=1):
    keypoints = np.ndarray(shape=DoG.shape, dtype=np.ndarray)
    sigma = 1.6
    sigmaCount = DoG[0].shape[0]
    for o in range(DoG.shape[0]):
        keypoints[o] = np.ndarray(shape=(DoG[o].shape[0]-(2*radius),), dtype=list)
        for s in range(radius, DoG[o].shape[0]-radius):
            keypoints[o][s-radius] = []
            k = 2**(float(s)/float(sigmaCount))
            for y in range(radius, DoG[o].shape[1]-radius):
                for x in range(radius, DoG[o].shape[2]-radius):
                    value = DoG[o][s, y, x]
                    neighbourhood = getNeighbourhood(DoG[o], s, y, x, radius=radius).flatten()
                    neighbourhood.sort()
                    min2 = neighbourhood[1]
                    max2 = neighbourhood[-2]
                    if value < min2 or (value > threshold and value > max2):
                        scale = 2**o
                        keypoints[o][s-radius].append((scale * y + scale/2, scale * x + scale/2, scale * k*sigma))
    return keypoints
def normalize(img):
    normImg = np.ndarray(shape=img.shape, dtype=np.float64)
    max_val = img.max()
    if max > 0:
        normImg = img/float(max_val)
        normImg *= 255.
    else:
        return img.copy()
    return normImg.astype(np.uint8)
def scale(img, factor=2):
    assert(len(img.shape) == 2)
    rows, cols = img.shape
    scaledImg = np.ndarray((rows*factor, cols*factor), np.float64)
    for y in range(0, scaledImg.shape[0]):
        for x in range(0, scaledImg.shape[1]):
            scaledImg[y][x] = img[y/factor][x/factor]
    return scaledImg
def downscale(img):
    assert(len(img.shape) == 2)
    rows, cols = img.shape
    scaledImg = np.ndarray((rows/2, cols/2), np.float64)
    for y in range(0, scaledImg.shape[0]):
        for x in range(0, scaledImg.shape[1]):
            scaledImg[y][x] = img[2*y][2*x]
    return scaledImg
def drawKeypoints(img, kp):
    if (len(img.shape) < 3 or img.shape[2] == 1):
        kpImg = cv2.cvtColor(img, cv2.COLOR_GRAY2BGR)
    else:
        kpImg = img.copy()
    for y, x, scale in kp:
        r = rnd.randrange(0,255)
        g = rnd.randrange(0,255)
        b = rnd.randrange(0,255)
        cv2.circle(kpImg, (int(x), int(y)), int(scale), (r, g, b))
    
    return kpImg
def plotImage(title,image):
    plt.figure(figsize=(20,12))
    gray()
    plt.imshow(image)
    plt.title(title)
    plt.axis('off')
    plt.show()
Explanation: Quiz 3b - Implementation of the SIFT Feature Detector
End of explanation
images = glob.glob('boat_images/img_*.png')
print('Images Loaded!')    
for filename in images:
    img = cv2.imread(filename, 0)
    gp = calcGaussianPyramid(img)
    DoG = calcDoG(gp)
    radius = 1
    keypoints = calcExtrema(DoG, radius=radius)
    kpImg = img.copy()
    for o in range(keypoints.shape[0]):
        for s in range(radius, DoG[o].shape[0]-radius):
            kp = keypoints[o][s-radius]
            kpImg = drawKeypoints(kpImg, kp)
    plotImage("Custom SIFT "+ filename, kpImg)    
    
        
Explanation: Execute custom SIFT for every given image in the current path then show keypoints as circles in the correspoding images.
SIFT algorithm begins detecting points that are invariant to scale changes. This is achieved by building a scale-space function named L(x,y,sigma) that consists of a convolution of a gaussian function and the image for every pixel. Following Lowe's paper (LOWE,2004) to efficiently detect keypoint locations in space scale we have proposed using a space scale extrema in the difference of Gaussian function convolved with the image.
In this implementation, the functions calcGaussianPyramid and getGaussianKernel compute the octaves and build a scape space for the image.
Afterwards, the difference of two consecutive guassian spaces separated by a constant factor k is carried out. This is implemented in the function calcDoG
Now, a local extrema detection step is executed. As Lowe stated each sample point is compared to each of its neighbors in the current image and nine of neighbors in the scale above and below.
Only the largest or the smallest of the points is selected as a candidate keypoint extrema.
To improve matching and stability a detailed fit to nearby data for location, scale and ratio of principal curvatures is performed. The function calcExtrema executes this step. This eliminates points with low contrast or that poorly localized along an edge.
And finally, the selected keypoints are drawn on the image.
End of explanation
#img = (Image.open('boat_images/img1.pgm').convert('L'))
#img.save('boat_images/img_1.png')
#img = (Image.open('boat_images/img2.pgm').convert('L'))
#img.save('boat_images/img_2.png')
#img = (Image.open('boat_images/img3.pgm').convert('L'))
#img.save('boat_images/img_3.png')
#img = (Image.open('boat_images/img4.pgm').convert('L'))
#img.save('boat_images/img_4.png')
#img = (Image.open('boat_images/img5.pgm').convert('L'))
#img.save('boat_images/img_5.png')
#img = (Image.open('boat_images/img6.pgm').convert('L'))
#img.save('boat_images/img_6.png')
images = glob.glob('boat_images//img*.pgm')
for filename in images:
    img = cv2.imread(filename, 0)
    sift = cv2.SIFT()
    kp, desc = sift.detectAndCompute(img, None)
    imgfinal=cv2.drawKeypoints(img,kp,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
    plotImage("SIFT Opencv  "+filename, imgfinal)   
## matching
sift = cv2.SIFT()
imgA = cv2.imread('boat_images//img1.pgm', 0)
imgB = cv2.imread('boat_images//img2.pgm', 0)
kpA, desA = sift.detectAndCompute(imgA,None)
kpB, desB = sift.detectAndCompute(imgB,None)
bf = cv2.BFMatcher()
matches = bf.knnMatch(desA,desB, k=2)
#img3 = cv2.drawMatchesKnn(imgA,kpA,imgB,kpB,good,flags=2)
#plt.imshow(img3)
#plt.show()
Explanation: SIFT Feature Detector - Opencv Implementation
End of explanation
<END_TASK>