diff --git a/__pycache__/convolution.cpython-37.pyc b/__pycache__/convolution.cpython-37.pyc
new file mode 100644
index 0000000..0a28000
Binary files /dev/null and b/__pycache__/convolution.cpython-37.pyc differ
diff --git a/__pycache__/harris.cpython-37.pyc b/__pycache__/harris.cpython-37.pyc
new file mode 100644
index 0000000..5ff47a5
Binary files /dev/null and b/__pycache__/harris.cpython-37.pyc differ
diff --git a/__pycache__/keypoint_detector.cpython-37.pyc b/__pycache__/keypoint_detector.cpython-37.pyc
new file mode 100644
index 0000000..ee2de8d
Binary files /dev/null and b/__pycache__/keypoint_detector.cpython-37.pyc differ
diff --git a/chessboard.png b/chessboard.png
new file mode 100644
index 0000000..29e51fe
Binary files /dev/null and b/chessboard.png differ
diff --git a/class_photo1.jpg b/class_photo1.jpg
new file mode 100644
index 0000000..40cb998
Binary files /dev/null and b/class_photo1.jpg differ
diff --git a/class_photo1re.jpg b/class_photo1re.jpg
new file mode 100644
index 0000000..b8df5e4
Binary files /dev/null and b/class_photo1re.jpg differ
diff --git a/class_photo2.jpg b/class_photo2.jpg
new file mode 100644
index 0000000..106eb15
Binary files /dev/null and b/class_photo2.jpg differ
diff --git a/class_photo2re.jpg b/class_photo2re.jpg
new file mode 100644
index 0000000..78188c7
Binary files /dev/null and b/class_photo2re.jpg differ
diff --git a/convolution.py b/convolution.py
new file mode 100644
index 0000000..c87a641
--- /dev/null
+++ b/convolution.py
@@ -0,0 +1,158 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from PIL import Image
+import sys
+
+# convolve image g with kernel h
+def convolve(g, h, std=False, color_channel=False):
+    g = g.astype(dtype='float32')
+
+    # image dimensions and color channels (add color channel if grayscale)
+    g_shape = None
+    if (len(g.shape) == 2):
+        g_shape = g.shape
+        height, width = g.shape
+        g = g.reshape((height, width, 1))
+    height, width, colors = g.shape
+    g_h = np.zeros_like(g)
+
+    # size of neighborhood
+    size = (int)((h.shape[0]-1)/2)
+
+    # iterate each color band
+    for color in range(colors):
+        # iterate each (u,v) pixel (offset 1 from edges)
+        for u in range(size, height-size):
+            for v in range(size, width-size):
+                g_h[u,v,color] = convolve_pixel(g, h, size, u, v, color)
+
+    if std == True:
+        standardize_color_channel(g, width, height, colors, size)
+
+    # reshape to grayscale if necessary
+    if (g_shape != None):
+        g_h = g_h.reshape(g_shape)
+
+    return g_h
+
+# find the convolution of (g*h)(u,v) at pixel [u,v,color]
+def convolve_pixel(g, h, size, u, v, color):
+    g_neighborhood = g[u-size:u+size+1, v-size:v+size+1, color]
+    g_h = 0
+
+    # convolve over 3x3 neighborhood
+    g_h = np.multiply(g_neighborhood,h)
+    g_h = g_h.sum()
+
+    return g_h
+
+# standardize color channels over range (0,255)
+def standardize_color_channel(g, width, height, colors, size):
+    # highest and lowest intensity
+    I_max = None
+    I_min = None
+    for color in range(colors):
+        for u in range(size, height-size):
+            for v in range(size, width-size):
+                if I_min == None or g[u,v,color] < I_min:
+                    I_min = g[u,v,color]
+                if  I_max == None or g[u,v,color] > I_max:
+                    I_max = g[u,v,color]
+    I_range = I_max - I_min
+    for color in range(colors):
+        for u in range(1, height-1):
+            for v in range(1, width-1):
+                g[u,v,color] = (255 / I_range) * (g[u,v,color] + I_min)
+    return g
+
+# Linear (Box) Blur kernel
+def linear_kernel(size):
+    n = 1+2*size
+    h = np.zeros((n, n))
+
+    for j in range(-size, size+1):
+        for k in range(-size, size+1):
+            h[j+size,k+size] = 1
+
+    # Calculate Unity constant
+    Z = 1/np.sum(h)
+    h = Z * h
+
+    return h
+
+# Gaussian Blur kernel
+def gaussian_kernel(sigma, size):
+    n = 1+2*size
+    h = np.zeros((n, n))
+
+    for j in range(-size, size+1):
+        for k in range(-size, size+1):
+            h[j+size,k+size] = np.exp(-1 * (np.power(j,2) + np.power(k,2)) / (2 * np.power(sigma,2)))
+
+    # Calculate Unity constant
+    Z = 1/np.sum(h)
+    h = Z * h
+
+    return h
+
+# Sobel Operator kernel (u), for Edge Detection
+def sobel_u_kernel():
+    h = np.array([[-1, 0, 1],
+                  [-2, 0, 2],
+                  [-1, 0, 1]])
+    return h
+
+# Sobel Operator kernel (v), for Edge Detection
+def sobel_v_kernel():
+    h = np.array([[-1, -2, -1],
+                  [ 0,  0,  0],
+                  [ 1,  2,  1]])
+    return h
+
+# convert image to grayscale using "perceptual luminance-preserving formula"
+def rgb2gray(rgb):
+    return np.dot(rgb[...,:3], [0.299, 0.587, 0.114])
+
+#############################################################################
+##############################       MAIN     ###############################
+#############################################################################
+
+# # import image
+# img_color = plt.imread('noisy_big_chief.jpeg')
+# img_gray = img_color.mean(axis=2)
+#
+# # show both image
+# # fig, ax = plt.subplots(1, figsize=(12,8))
+# # plt.imshow(img_color)
+# # plt.show()
+# # plt.imshow(img_gray, cmap='gray')
+# # plt.show()
+# # plt.imsave('bnc_gray.jpeg', img_gray, cmap='gray')
+#
+# # Linear Blur/Smoothing
+# h_linear = linear_kernel(5)
+# img_linear = convolve(img_color, h_linear).astype(dtype='uint8')
+# plt.imshow(img_linear)
+# plt.show()
+# plt.imsave('bnc_linear.jpeg', img_linear)
+#
+# # Gaussian Blur/Smoothing
+# h_gaussian = gaussian_kernel(3, 5)
+# img_gaussian = convolve(img_color, h_gaussian).astype(dtype='uint8')
+# plt.imshow(img_gaussian)
+# plt.show()
+# plt.imsave('bnc_gaussian.jpeg', img_gaussian)
+#
+# # Sobel-u Edge Detection
+# h_sobelu = sobel_u_kernel()
+# img_sobelu = convolve(img_gray, h_sobelu)
+# plt.imshow(img_sobelu, cmap='gray')
+# plt.show()
+# plt.imsave('bnc_sobelu.jpeg', img_sobelu, cmap='gray')
+#
+# # Sobel-v Edge Detection
+# h_sobelv = sobel_v_kernel()
+# img_sobelv = convolve(img_gray, h_sobelv)
+# plt.imshow(img_sobelv, cmap='gray')
+# plt.show()
+# plt.imsave('bnc_sobelv.jpeg', img_sobelv, cmap='gray')
diff --git a/keypoint_detector.py b/keypoint_detector.py
new file mode 100644
index 0000000..2a37914
--- /dev/null
+++ b/keypoint_detector.py
@@ -0,0 +1,37 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from PIL import Image
+import sys
+
+# local imports
+from convolution import *
+
+# find local maxima of Harris Response matrix
+def local_maxima(g, h, size):
+    return 0
+
+# Check if a given pixel is a local maxima
+def check_if_local_maxima(H, size, width, height, u, v):
+    # test pixel Harris score
+    px = H[u,v]
+    # iterate over neighborhood
+    for j in range(u-size, u+size+1):
+        for k in range(v-size, v+size+1):
+            if j<0 or k<0 or j>height-1 or k>width-1:
+                continue
+            # skip itself
+            if j==u and k==v:
+                continue
+            # if pixel of greater score found, not maxima
+            elif px < H[j,k]:
+                return False
+    # else, is maxima
+    return True
+
+def sum_sq_error(p1, p2):
+	sum = 0
+
+	for i, j in zip(p1, p2):
+		sum += (i - j)**2
+
+	return (sum)
diff --git a/noisy_big_chief.jpeg b/noisy_big_chief.jpeg
new file mode 100644
index 0000000..621b299
Binary files /dev/null and b/noisy_big_chief.jpeg differ
diff --git a/project_description.ipynb b/project_description.ipynb
index 421c6d0..be28190 100644
--- a/project_description.ipynb
+++ b/project_description.ipynb
@@ -90,7 +90,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.5"
+   "version": "3.7.1"
   }
  },
  "nbformat": 4,
diff --git a/stitcher.py b/stitcher.py
new file mode 100644
index 0000000..5681edd
--- /dev/null
+++ b/stitcher.py
@@ -0,0 +1,148 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from keypoint_detector import *
+
+class Stitcher(object):
+	def __init__(self,image_1,image_2):
+		self.images = [image_1,image_2]
+		self.images[0] = self.images[0].mean(axis=2)
+		self.images[1] = self.images[1].mean(axis=2)
+
+	def find_keypoints(self, img):
+		"""
+		Step 1: This method locates features that are "good" for matching.  To do this we will implement the Harris 
+		corner detector
+		"""
+		# import image and convert to grayscale
+		I = img
+
+		# detect_keypoints(I)
+		g = I
+
+		h_sobelu = sobel_u_kernel()
+		h_sobelv = sobel_v_kernel()
+		h_gaussian = gaussian_kernel(2, 2)
+
+		g_u = convolve(g, h_sobelu)
+		g_u2 = np.multiply(g_u, g_u)
+		g_u2 = convolve(g_u2, h_gaussian)
+
+		g_v = convolve(g, h_sobelv)
+		g_v2 = np.multiply(g_v, g_v)
+		g_v2 = convolve(g_v2, h_gaussian)
+
+		g_uv = np.multiply(g_u, g_v)
+		g_uv = convolve(g_uv, h_gaussian)
+
+		# Harris Response matrix
+		H = np.zeros_like(g)
+		# Determinant(g)
+		H_det = np.multiply(g_u2, g_v2) - np.multiply(g_uv, g_uv)
+		# Trace(g) (small number to prevent division by zero)
+		H_tr = g_u2 + g_v2 + 1e-10
+		# H = det(g)/tr(g)
+		H = np.divide(H_det, H_tr)
+		
+		# Find local maxima of Harris matrix scores
+		height, width = H.shape
+		max = []
+		# size of neighborhood
+		size = 50
+
+		# iterate over pixels
+		for u in range(height):
+			for v in range(width):
+				# add to list if local maxima
+				if check_if_local_maxima(H, size, width, height, u, v):
+					max.append([v,u])
+		
+		for i in range(len(max)):
+			plt.scatter(x=max[i][0], y=max[i][1], c='r')
+
+		plt.imshow(I,cmap=plt.cm.gray)
+		plt.show()
+		
+		return max
+
+	def generate_descriptors(self,img,points, l):
+		"""
+		Step 2: After identifying relevant keypoints, we need to come up with a quantitative description of the 
+		neighborhood of that keypoint, so that we can match it to keypoints in other images.
+		"""
+		l = int(l/2)
+		
+		descriptors = []
+		for i in range(len(points)):
+			x1 = points[i][0]-l
+			x2 = points[i][0]+l+1
+			y1 = points[i][1]-l
+			y2 = points[i][1]+l+1
+			
+			if(x1 < 0):
+				x1 = 0
+			if(y1 < 0):
+				y1 = 0
+			if(x2 > img.shape[0]):
+				x2 = img.shape[0]
+			if(y2 > img.shape[1]):
+				y2 = img.shape[1]
+			
+			descriptors.append([img[x1:x2:1, y1:y2:1], points[i]]) #Returns an array with [actual img data, img location (u,v)]
+		return descriptors
+	
+	def match_keypoints(self):
+		"""
+		Step 3: Compare keypoint descriptions between images, identify potential matches, and filter likely
+		mismatches
+		"""
+		matches = []
+		m1 = myStitcher.find_keypoints(self.images[0])
+		m2 = myStitcher.find_keypoints(self.images[1])
+		desc1 = myStitcher.generate_descriptors(self.images[0], m1, 21)
+		desc2 = myStitcher.generate_descriptors(self.images[1], m2, 21)
+		
+		for p1 in desc1:
+			x1, y1 = p1[1][0], p1[1][1]
+			errors = []
+			for p2 in desc2:
+				x2, y2 = p2[1][0], p2[1][1]
+				errors.append([[x2, y2], sum_sq_error(p1[0].flatten(), p2[0].flatten())])
+			errors.sort(key=lambda x: x[1])
+			matches.append([[x1, y1], errors[0][0]])
+		return matches
+
+	def find_homography(self):
+		"""
+		Step 4: Find a linear transformation (of various complexities) that maps pixels from the second image to 
+		pixels in the first image
+		"""
+
+	def stitch(self):
+		"""
+		Step 5: Transform second image into local coordinate system of first image, and (perhaps) perform blending
+		to avoid obvious seams between images.
+		"""
+
+im1 = plt.imread('class_photo1re.jpg')
+im2 = plt.imread('class_photo2re.jpg')
+myStitcher = Stitcher(im1, im2)
+matches = myStitcher.match_keypoints()
+
+I = np.concatenate((im1, im2), axis=1)
+
+plt.imshow(I, cmap=plt.cm.gray)
+
+offset = im1.shape[1]
+lengthOfMatches = len(matches)
+
+x1 = [matches[i][0][0] for i in range(lengthOfMatches)]
+y1 = [matches[i][0][1] for i in range(lengthOfMatches)]
+
+x2 = [matches[i][1][0] + offset for i in range(lengthOfMatches)]
+y2 = [matches[i][1][1] for i in range(lengthOfMatches)]
+
+for i in range(lengthOfMatches):
+	plt.plot([x1[i], x2[i]], [y1[i], y2[i]])
+
+plt.show()