diff --git a/Descriptor.h b/Descriptor.h
index b3e523b..cad0da9 100644
--- a/Descriptor.h
+++ b/Descriptor.h
@@ -7,22 +7,19 @@
  * this code :)
  ******************************************************/
 
+#pragma once
 #ifndef DESCRIPTOR_H
 #define DESCRIPTOR_H
-
+#include<vector>
 using namespace std;
 
 class Descriptor
 {
 public:
 	float			xi, yi;		// The location
-	vector<double>	fv;			// The feature vector
-
-	Descriptor()
-	{
-	}
+	vector<float>	fv;			// The feature vector
 
-	Descriptor(float x, float y, vector<double> const& f)
+	Descriptor(float x, float y, vector<float> const& f)
 	{
 		xi = x;
 		yi = y;
@@ -30,4 +27,4 @@ class Descriptor
 	}
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/KeyPoint.h b/KeyPoint.h
index 745b385..9d6053b 100644
--- a/KeyPoint.h
+++ b/KeyPoint.h
@@ -7,6 +7,7 @@
  * this code :)
  ******************************************************/
 
+#pragma once
 #ifndef _KEYPOINT_H
 #define _KEYPOINT_H
 
@@ -20,13 +21,12 @@ class Keypoint
 public:
 	float			xi;
 	float			yi;	// It's location
-	vector<double>	mag;	// The list of magnitudes at this point
-	vector<double>	orien;	// The list of orientations detected
-	unsigned int	scale;	// The scale where this was detected
-
-	Keypoint() { }
-	Keypoint(float x, float y) { xi=x; yi=y; }
-	Keypoint(float x, float y, vector<double> const& m, vector<double> const& o, unsigned int s)
+	vector<float>	mag;	// The list of magnitudes at this point
+	vector<float>	orien;	// The list of orientations detected
+	int	scale;	// The scale where this was detected
+	
+	//defining the constructor
+	Keypoint(float x, float y, vector<float>m, vector<float>o, int s)
 	{
 		xi = x;
 		yi = y;
@@ -36,4 +36,4 @@ class Keypoint
 	}
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/MySIFT.cpp b/MySIFT.cpp
index 2ece459..03d7c2e 100644
--- a/MySIFT.cpp
+++ b/MySIFT.cpp
@@ -10,26 +10,18 @@
 // MySIFT.cpp : Defines the entry point for the console application.
 //
 
-#include "stdafx.h"
+#include<iostream>
+#include<opencv2/opencv.hpp>
+#include "mySIFT.h"
+using namespace cv;
+using namespace std;
 
-#include "SIFT.h"
-
-#include <cv.h>
-#include <highgui.h>
-
-// The main function!
 int main()
 {
-	// Create an instance of SIFT
-	SIFT *sift = new SIFT("C:\\house.jpg", 4, 2);
-
-	sift->DoSift();				// Find keypoints
-	sift->ShowAbsSigma();		// Display the sigma table
-	sift->ShowKeypoints();		// Show the keypoints
-	cvWaitKey(0);				// Wait for a keypress
-
-	// Cleanup and exit
-	delete sift;
+	Mat img = imread("car.jpg");
+	meSIFT x(img,4,2);
+	x.DoSift();
+	x.ShowKeypoints();
 	return 0;
 }
 
diff --git a/SIFT.cpp b/SIFT.cpp
index ae60532..f82c935 100644
--- a/SIFT.cpp
+++ b/SIFT.cpp
@@ -7,10 +7,11 @@
  * this code :)
  ******************************************************/
 
-#include "stdafx.h"
-
 #include "SIFT.h"
-
+#include "stdafx.h"
+#include <climits>
+//This file containes all the method of 
+//the class meSIFT
 #define SIGMA_ANTIALIAS			0.5
 #define SIGMA_PREBLUR			1.0
 #define CURVATURE_THRESHOLD		5.0
@@ -23,1010 +24,813 @@
 #define FVSIZE					128
 #define	FV_THRESHOLD			0.2
 
-// SaveFloatingPointImage()
-// The standard HighGUI functions can save only 8bit images. This
-// function converts a floating point image (with values 0..1) into
-// a normal 8bit image (with values 0..255), and then saves it.
-void SaveFloatingPointImage(const char *filename, IplImage* img)
-{
-	IplImage* dup = cvCreateImage(cvGetSize(img), 8, 1);
-	cvCvtScale(img, dup, 255.0);
-
-	cvSaveImage(filename, dup);
-	cvReleaseImage(&dup);
-}
+using namespace cv;
 
-// Constructor: Give a preloaded image and provide the desired number
-// of octaves and intervals
-SIFT::SIFT(IplImage* img, int octaves, int intervals)
+//intialising the the meSIFT class
+//with a constructor,taking 3 arguments
+//the src image, the number of octaves in the image
+//and the number of intervals
+meSIFT::meSIFT(Mat src,int octaves,int intervals)
 {
-	// Store the image internally
-	m_srcImage = cvCloneImage(img);
-
-	// Set the number of octaves and intervals
+	srcImage = src.clone();
 	m_numOctaves = octaves;
 	m_numIntervals = intervals;
-
-	// Proceed to initialize the algorithm
-	GenerateLists();
-}
-
-// Constructor: Give a filename and the desired number of octaves
-// and intervals
-SIFT::SIFT(const char* filename, int octaves, int intervals)
-{
-	// Load the image
-	m_srcImage = cvLoadImage(filename);
-
-	// Set the number of octaves and intervals
-	m_numOctaves = octaves;
-	m_numIntervals = intervals;
-
-	// Proceed to initialize the algorithm
-	GenerateLists();
-}
-
-// AllocateMemory()
-// This function allocates memory for the scale space (the multiple gaussian images,
-// the difference of guassian images)
-void SIFT::GenerateLists()
-{
-	// A variable for the loops
-	unsigned int i=0;
-
-	// Create a 2D array of gaussian blurred images
-	m_gList = new IplImage**[m_numOctaves];
-	for(i=0;i<m_numOctaves;i++)
-		m_gList[i] = new IplImage*[m_numIntervals+3];
-
-	// Create a 2D array to store images generated after the
-	// DoG operation
-	m_dogList = new IplImage**[m_numOctaves];
-	for(i=0;i<m_numOctaves;i++)
-		m_dogList[i] = new IplImage*[m_numIntervals+2];
-
-	// Create a 2D array that will hold if a particular point
-	// is an extrema or not
-	m_extrema = new IplImage**[m_numOctaves];
-	for(i=0;i<m_numOctaves;i++)
-		m_extrema[i] = new IplImage*[m_numIntervals];
-
-	// Create a 2D array of decimal numbers. It holds the sigma
-	// used to blur the gaussian images.
-	m_absSigma = new double*[m_numOctaves];
-	for(i=0;i<m_numOctaves;i++)
-		m_absSigma[i] = new double[m_numIntervals+3];
-}
-
-// Destructor
-// Cleanup after you're done
-SIFT::~SIFT()
-{
-	unsigned int i, j;
-	for(i=0;i<m_numOctaves;i++)
-	{
-		// Release all images in that particular octave
-		for(j=0;j<m_numIntervals+3;j++)	cvReleaseImage(&m_gList[i][j]);
-		for(j=0;j<m_numIntervals+2;j++)	cvReleaseImage(&m_dogList[i][j]);
-		for(j=0;j<m_numIntervals;j++)	cvReleaseImage(&m_extrema[i][j]);
-
-		// Delete memory for that array
-		delete [] m_gList[i];
-		delete [] m_dogList[i];
-		delete [] m_extrema[i];
-		delete [] m_absSigma[i];
-	}
-
-	// Delete the 2D arrays
-    delete [] m_gList;
-    delete [] m_dogList;
-    delete [] m_extrema;
-    delete [] m_absSigma;
 }
 
 // DoSift()
 // This function does everything in sequence.
-void SIFT::DoSift()
+void meSIFT::DoSift()
 {
 	BuildScaleSpace();
 	DetectExtrema();
 	AssignOrientations();
-	ExtractKeypointDescriptors();
+	//ExtractKeypointDescriptors();
+	return;
 }
 
 // BuildScaleSpace()
-// This function generates all the blurred out images for each octave
+// This function generates al`l the blurred out images for each octave
 // and also the DoG images
-void SIFT::BuildScaleSpace()
+void meSIFT::BuildScaleSpace()
 {
-	printf("Generating scale space...\n");
-	// For loops
-	unsigned int i,j;
-
-	// floating point grayscale image
-	IplImage* imgGray = cvCreateImage(cvGetSize(m_srcImage), IPL_DEPTH_32F , 1);
-	IplImage* imgTemp = cvCreateImage(cvGetSize(m_srcImage), 8 , 1);
-	
-	// Create a duplicate. We don't want to mess the original
-	// If the image is colour, it is converted to grayscale
-	if(m_srcImage->nChannels==3)
-	{
-		cvCvtColor(m_srcImage, imgTemp, CV_BGR2GRAY);
-	}
-	else
-	{
-		cvCopy(m_srcImage, imgTemp);
-	}
 
-	// Finally, generate the floating point image... convert 0..255 range into 0..1
-	for(int x=0;x<imgTemp->width;x++)
-	{
-		for(int y=0;y<imgTemp->height;y++)
-		{
-			cvSetReal2D(imgGray, y, x, cvGetReal2D(imgTemp, y, x)/255.0);
-		}
-	}
+	Mat src_gray;
+	//intialising a floating point image
+	Mat imgFloat;
+	//we need a gray scale image 
+	cvtColor(srcImage,src_gray, COLOR_BGR2GRAY);
+	//we set each and every pixel of the floating point image
+	//converting from range(0 to 255) to (0 to 1)
+	src_gray.convertTo(imgFloat,CV_32F,1.0/255,0);
 
 	// Lowe claims blur the image with a sigma of 0.5 and double it's dimensions
 	// to increase the number of stable keypoints
-	cvSmooth(imgGray, imgGray, CV_GAUSSIAN, 0, 0, SIGMA_ANTIALIAS);
-
-	// Create an image double the dimensions, resize imgGray and store it in m_gList[0][0]
-	m_gList[0][0] = cvCreateImage(cvSize(imgGray->width*2, imgGray->height*2), IPL_DEPTH_32F , 1);
-	cvPyrUp(imgGray, m_gList[0][0]);
-
-	// Preblur this base image
-	cvSmooth(m_gList[0][0], m_gList[0][0], CV_GAUSSIAN, 0, 0, SIGMA_PREBLUR);
-
-	// SaveFloatingPointImage("C:\\SIFT Test\\Gaussian\\g_octave_0_scale_0.jpg", m_gList[0][0]);
-
+	//so we apply Gaussian blur
+	GaussianBlur(imgFloat, imgFloat, Size(0, 0), SIGMA_ANTIALIAS);               //CHANGEABLE
+	resize(imgFloat,imgFloat,Size(imgFloat.size()*2));
+	m_glist.push_back(vector<Mat>());
+	m_glist[0].push_back(imgFloat.clone());
+	pyrUp(imgFloat, m_glist[0][0]);
+	GaussianBlur(m_glist[0][0],m_glist[0][0],Size(0,0),SIGMA_PREBLUR);
+	//imshow("0", m_glist[0][0]);
+	//waitKey(0);
+
+	//keeping a track of the sigmas
 	double initSigma = sqrt(2.0f);
+	m_absSigma.push_back(vector<double>());
+	m_absSigma[0].push_back(initSigma * 0.5);
 
-	// Keep a track of the sigmas
-	m_absSigma[0][0] = initSigma * 0.5;
-
-	// Now for the actual image generation
-	for(i=0;i<m_numOctaves;i++)
+	//now for the actual image generation
+	//we run this 2D loop
+	int i, j;
+	for (i = 0; i < m_numOctaves; i++)
 	{
 		// Reset sigma for each octave
 		double sigma = initSigma;
-		CvSize currentSize = cvGetSize(m_gList[i][0]);
-
-		for(j=1;j<m_numIntervals+3;j++)
+		Size s = m_glist[i][0].size();
+		m_dogList.push_back(vector<Mat>());
+		/*if (i > 0)
+		{
+			m_absSigma.push_back(vector<double>());
+			m_glist.push_back(vector<Mat>());
+		}*/
+		for (j = 1; j<int(m_numIntervals+3); j++)
 		{
-			// Allocate memory
-			m_gList[i][j] = cvCreateImage(currentSize, 32, 1);
-
 			// Calculate a sigma to blur the current image to get the next one
-			double sigma_f = sqrt(pow(2.0,2.0/m_numIntervals)-1) * sigma;
-            sigma = pow(2.0,1.0/m_numIntervals) * sigma;
-
+			double sigma_f = sqrt(pow(2.0, 2.0 / m_numIntervals) - 1) * sigma;
+			sigma = pow(2.0, 1.0 / m_numIntervals) * sigma;
 			// Store sigma values (to be used later on)
-			m_absSigma[i][j] = sigma * 0.5 * pow(2.0f, (float)i);
-
-			// Apply gaussian smoothing)
-			cvSmooth(m_gList[i][j-1], m_gList[i][j], CV_GAUSSIAN, 0, 0, sigma_f);
-
-			// Calculate the DoG image
-			m_dogList[i][j-1] = cvCreateImage(currentSize, 32, 1);
-			cvSub(m_gList[i][j-1], m_gList[i][j], m_dogList[i][j-1]);
-
-			// Save the images generated for fun :)
-			/*char* filename = new char[200];
-			sprintf(filename, "C:\\SIFT Test\\Gaussian\\g_octave_%d_scale_%d.jpg", i, j);
-			SaveFloatingPointImage(filename, m_gList[i][j]);
-
-			sprintf(filename, "C:\\SIFT TEST\\DOG\\dog_octave_%d_scale_%d.jpg", i, j-1);
-			SaveFloatingPointImage(filename, m_dogList[i][j-1]);*/
+			m_absSigma[i].push_back(sigma * 0.5 * pow(2.0f, (float)i));
+
+			//applying Gaussian blur
+			Mat temp;
+			GaussianBlur(m_glist[i][j-1], temp, Size(0, 0), sigma_f);
+			m_glist[i].push_back(temp);
+
+			//finding and allocating the DoG image using the vector
+			Mat temp2;
+			subtract(m_glist[i][j-1],m_glist[i][j],temp2,Mat());
+			m_dogList[i].push_back(temp2);
+			//imshow(to_string(10*j+i),m_glist[i][j]);
+			//waitKey(0);
+			//imshow(to_string(10*j+i),m_dogList[i][j-1]);
+			//waitKey(0);
 		}
-
-		// If we're not at the last octave
-		if(i<m_numOctaves-1)
+		if (i < m_numOctaves - 1)
 		{
-			// Reduce size to half
-			currentSize.width/=2;
-			currentSize.height/=2;
-
-			// Allocate memory and resample the image
-			m_gList[i+1][0] = cvCreateImage(currentSize, 32, 1);
-			cvPyrDown(m_gList[i][0], m_gList[i+1][0]);
-			m_absSigma[i+1][0] = m_absSigma[i][m_numIntervals];
-
-			// Store the image
-			/*char* filename = new char[200];
-			sprintf(filename, "C:\\SIFT Test\\Gaussian\\g_octave_%d_scale_0.jpg", i+1);
-			SaveFloatingPointImage(filename, m_gList[i+1][0]);*/
+			Mat temp3;
+			pyrDown(m_glist[i][0], temp3, s/2);
+			m_absSigma.push_back(vector<double>());
+			m_glist.push_back(vector<Mat>());
+			m_glist[i+1].push_back(temp3);
+			m_absSigma[i+1].push_back(m_absSigma[i][m_numIntervals]);
 		}
 	}
 }
 
 // DetectExtrema()
 // Locates extreme points (maxima and minima)
-// Relatively simple stuff
-void SIFT::DetectExtrema()
+void meSIFT::DetectExtrema()
 {
-	printf("Detecting extrema...\n");
-
-	// Looping variables
-	unsigned int i, j, xi, yi;
-
-	// Some variables we'll use later on
 	double curvature_ratio, curvature_threshold;
-	IplImage *middle, *up, *down;
 	int scale;
 	double dxx, dyy, dxy, trH, detH;
-
-	unsigned int num=0;				// Number of keypoins detected
-	unsigned int numRemoved=0;		// The number of key points rejected because they failed a test
-
-	curvature_threshold = (CURVATURE_THRESHOLD+1)*(CURVATURE_THRESHOLD+1)/CURVATURE_THRESHOLD;
-
-	// Detect extrema in the DoG images
-	for(i=0;i<m_numOctaves;i++)
+	unsigned int num = 0;				// Number of keypoints detected
+	unsigned int numRemoved = 0;		// The number of key points rejected because they failed a test
+	Mat middle, up, down;
+	curvature_threshold = (CURVATURE_THRESHOLD + 1) * (CURVATURE_THRESHOLD + 1) / CURVATURE_THRESHOLD;
+
+	//Detect extrema in the DoG image
+	//using almost four loops
+	for (int i = 0; i < m_numOctaves; i++)
 	{
-		scale = (int)pow(2.0, (double)i);
-		
-		for(j=1;j<m_numIntervals+1;j++)
+		scale= (int)pow(2.0, (double)i);
+		m_extrema.push_back(vector<Mat>());
+		for (int j = 1; j < m_numIntervals + 1; j++)
 		{
 			// Allocate memory and set all points to zero ("not key point")
-			m_extrema[i][j-1] = cvCreateImage(cvGetSize(m_dogList[i][0]), 8, 1);
-			cvZero(m_extrema[i][j-1]);
+			Mat temp = Mat::zeros(m_dogList[i][0].size(), CV_8UC1);
+			m_extrema[i].push_back(temp);
 
-			// Images just above and below, in the current octave
+			//finding the images just above and below
+			//for the current octave 
 			middle = m_dogList[i][j];
-			up = m_dogList[i][j+1];
-			down = m_dogList[i][j-1];
+			up = m_dogList[i][j + 1];
+			down = m_dogList[i][j - 1];
 
-			for(xi=1;xi<m_dogList[i][j]->width-1;xi++)
+			//intializing certain variables
+			int xi, yi;
+			for (xi = 1; xi < m_dogList[i][j].rows - 1; xi++)
 			{
-				for(yi=1;yi<m_dogList[i][j]->height-1;yi++)
+				for (yi = 1; yi < m_dogList[i][j].cols - 1; yi++)
 				{
 					// true if a keypoint is a maxima/minima
 					// but needs to be tested for contrast/edge thingy
 					bool justSet = false;
-
-					double currentPixel = cvGetReal2D(middle, yi, xi);
-
-					// Check for a maximum
-					if (currentPixel > cvGetReal2D(middle, yi-1, xi  )	&&
-                        currentPixel > cvGetReal2D(middle, yi+1, xi  )  &&
-                        currentPixel > cvGetReal2D(middle, yi  , xi-1)  &&
-                        currentPixel > cvGetReal2D(middle, yi  , xi+1)  &&
-                        currentPixel > cvGetReal2D(middle, yi-1, xi-1)	&&
-                        currentPixel > cvGetReal2D(middle, yi-1, xi+1)	&&
-                        currentPixel > cvGetReal2D(middle, yi+1, xi+1)	&&
-                        currentPixel > cvGetReal2D(middle, yi+1, xi-1)	&&
-                        currentPixel > cvGetReal2D(up, yi  , xi  )      &&
-                        currentPixel > cvGetReal2D(up, yi-1, xi  )      &&
-                        currentPixel > cvGetReal2D(up, yi+1, xi  )      &&
-                        currentPixel > cvGetReal2D(up, yi  , xi-1)      &&
-                        currentPixel > cvGetReal2D(up, yi  , xi+1)      &&
-                        currentPixel > cvGetReal2D(up, yi-1, xi-1)		&&
-                        currentPixel > cvGetReal2D(up, yi-1, xi+1)		&&
-                        currentPixel > cvGetReal2D(up, yi+1, xi+1)		&&
-                        currentPixel > cvGetReal2D(up, yi+1, xi-1)		&&
-                        currentPixel > cvGetReal2D(down, yi  , xi  )    &&
-                        currentPixel > cvGetReal2D(down, yi-1, xi  )    &&
-                        currentPixel > cvGetReal2D(down, yi+1, xi  )    &&
-                        currentPixel > cvGetReal2D(down, yi  , xi-1)    &&
-                        currentPixel > cvGetReal2D(down, yi  , xi+1)    &&
-                        currentPixel > cvGetReal2D(down, yi-1, xi-1)	&&
-                        currentPixel > cvGetReal2D(down, yi-1, xi+1)	&&
-                        currentPixel > cvGetReal2D(down, yi+1, xi+1)	&&
-                        currentPixel > cvGetReal2D(down, yi+1, xi-1)   )
+					float curr = middle.at<float>(Point(yi, xi));
+					if (curr > middle.at<float>(Point(yi + 1, xi)) &&
+						curr > middle.at<float>(Point(yi - 1, xi)) &&
+						curr > middle.at<float>(Point(yi, xi + 1)) &&
+						curr > middle.at<float>(Point(yi, xi - 1)) &&
+						curr > middle.at<float>(Point(yi + 1, xi + 1)) &&
+						curr > middle.at<float>(Point(yi - 1, xi - 1)) &&
+						curr > middle.at<float>(Point(yi - 1, xi + 1)) &&
+						curr > middle.at<float>(Point(yi + 1, xi - 1)) &&
+						curr > up.at<float>(Point(yi, xi)) &&
+						curr > up.at<float>(Point(yi + 1, xi)) &&
+						curr > up.at<float>(Point(yi - 1, xi)) &&
+						curr > up.at<float>(Point(yi, xi + 1)) &&
+						curr > up.at<float>(Point(yi, xi - 1)) &&
+						curr > up.at<float>(Point(yi + 1, xi + 1)) &&
+						curr > up.at<float>(Point(yi - 1, xi - 1)) &&
+						curr > up.at<float>(Point(yi - 1, xi + 1)) &&
+						curr > up.at<float>(Point(yi + 1, xi - 1)) &&
+						curr > down.at<float>(Point(yi, xi)) &&
+						curr > down.at<float>(Point(yi + 1, xi)) &&
+						curr > down.at<float>(Point(yi - 1, xi)) &&
+						curr > down.at<float>(Point(yi, xi + 1)) &&
+						curr > down.at<float>(Point(yi, xi - 1)) &&
+						curr > down.at<float>(Point(yi + 1, xi + 1)) &&
+						curr > down.at<float>(Point(yi - 1, xi - 1)) &&
+						curr > down.at<float>(Point(yi - 1, xi + 1)) &&
+						curr > down.at<float>(Point(yi + 1, xi - 1))
+						)
 					{
-						cvSetReal2D(m_extrema[i][j-1], yi, xi, 255);
 						num++;
 						justSet = true;
+						m_extrema[i][j - 1].at<uchar>(Point(yi, xi)) = 255;
 					}
-					// Check if it's a minimum
-					else if (currentPixel < cvGetReal2D(middle, yi-1, xi  )	&&
-                        currentPixel < cvGetReal2D(middle, yi+1, xi  )  &&
-                        currentPixel < cvGetReal2D(middle, yi  , xi-1)  &&
-                        currentPixel < cvGetReal2D(middle, yi  , xi+1)  &&
-                        currentPixel < cvGetReal2D(middle, yi-1, xi-1)	&&
-                        currentPixel < cvGetReal2D(middle, yi-1, xi+1)	&&
-                        currentPixel < cvGetReal2D(middle, yi+1, xi+1)	&&
-                        currentPixel < cvGetReal2D(middle, yi+1, xi-1)	&&
-                        currentPixel < cvGetReal2D(up, yi  , xi  )      &&
-                        currentPixel < cvGetReal2D(up, yi-1, xi  )      &&
-                        currentPixel < cvGetReal2D(up, yi+1, xi  )      &&
-                        currentPixel < cvGetReal2D(up, yi  , xi-1)      &&
-                        currentPixel < cvGetReal2D(up, yi  , xi+1)      &&
-                        currentPixel < cvGetReal2D(up, yi-1, xi-1)		&&
-                        currentPixel < cvGetReal2D(up, yi-1, xi+1)		&&
-                        currentPixel < cvGetReal2D(up, yi+1, xi+1)		&&
-                        currentPixel < cvGetReal2D(up, yi+1, xi-1)		&&
-                        currentPixel < cvGetReal2D(down, yi  , xi  )    &&
-                        currentPixel < cvGetReal2D(down, yi-1, xi  )    &&
-                        currentPixel < cvGetReal2D(down, yi+1, xi  )    &&
-                        currentPixel < cvGetReal2D(down, yi  , xi-1)    &&
-                        currentPixel < cvGetReal2D(down, yi  , xi+1)    &&
-                        currentPixel < cvGetReal2D(down, yi-1, xi-1)	&&
-                        currentPixel < cvGetReal2D(down, yi-1, xi+1)	&&
-                        currentPixel < cvGetReal2D(down, yi+1, xi+1)	&&
-                        currentPixel < cvGetReal2D(down, yi+1, xi-1)   )
+
+					//checking if it's a minima
+					if (curr < middle.at<float>(Point(yi + 1, xi)) &&
+						curr < middle.at<float>(Point(yi - 1, xi)) &&
+						curr < middle.at<float>(Point(yi, xi + 1)) &&
+						curr < middle.at<float>(Point(yi, xi - 1)) &&
+						curr < middle.at<float>(Point(yi + 1, xi + 1)) &&
+						curr < middle.at<float>(Point(yi - 1, xi - 1)) &&
+						curr < middle.at<float>(Point(yi - 1, xi + 1)) &&
+						curr < middle.at<float>(Point(yi + 1, xi - 1)) &&
+						curr < up.at<float>(Point(yi, xi)) &&
+						curr < up.at<float>(Point(yi + 1, xi)) &&
+						curr < up.at<float>(Point(yi - 1, xi)) &&
+						curr < up.at<float>(Point(yi, xi + 1)) &&
+						curr < up.at<float>(Point(yi, xi - 1)) &&
+						curr < up.at<float>(Point(yi + 1, xi + 1)) &&
+						curr < up.at<float>(Point(yi - 1, xi - 1)) &&
+						curr < up.at<float>(Point(yi - 1, xi + 1)) &&
+						curr < up.at<float>(Point(yi + 1, xi - 1)) &&
+						curr < down.at<float>(Point(yi, xi)) &&
+						curr < down.at<float>(Point(yi + 1, xi)) &&
+						curr < down.at<float>(Point(yi - 1, xi)) &&
+						curr < down.at<float>(Point(yi, xi + 1)) &&
+						curr < down.at<float>(Point(yi, xi - 1)) &&
+						curr < down.at<float>(Point(yi + 1, xi + 1)) &&
+						curr < down.at<float>(Point(yi - 1, xi - 1)) &&
+						curr < down.at<float>(Point(yi - 1, xi + 1)) &&
+						curr < down.at<float>(Point(yi + 1, xi - 1))
+						)
 					{
-						cvSetReal2D(m_extrema[i][j-1], yi, xi, 255);
 						num++;
 						justSet = true;
+						m_extrema[i][j - 1].at<uchar>(Point(yi, xi)) = 255;
 					}
 
-					// The contrast check
-					if(justSet && fabs(cvGetReal2D(middle, yi, xi)) < CONTRAST_THRESHOLD)
+					//contrast check
+					if (justSet && fabs(middle.at<float>(Point(yi, xi))) < CONTRAST_THRESHOLD)
 					{
-						cvSetReal2D(m_extrema[i][j-1], yi, xi, 0);
 						num--;
 						numRemoved++;
-
-						justSet=false;
+						m_extrema[i][j - 1].at<uchar>(Point(yi, xi)) = 0;
+						justSet = false;
 					}
-
-					// The edge check
-					if(justSet)
+					//edge check
+					if (justSet)
 					{
-						dxx = (cvGetReal2D(middle, yi-1, xi) +
-							  cvGetReal2D(middle, yi+1, xi) -
-							  2.0*cvGetReal2D(middle, yi, xi));
-
-						dyy = (cvGetReal2D(middle, yi, xi-1) +
-							  cvGetReal2D(middle, yi, xi+1) -
-							  2.0*cvGetReal2D(middle, yi, xi));
+						dxx = (middle.at<float>(Point(yi + 1, xi)) + middle.at<float>(Point(yi - 1, xi)) - 2.0 * middle.at<float>(Point(yi, xi)));
+						dyy = (middle.at<float>(Point(yi, xi + 1)) + middle.at<float>(Point(yi, xi - 1)) - 2.0 * middle.at<float>(Point(yi, xi)));
+						dxy = (middle.at<float>(Point(yi + 1, xi + 1)) + middle.at<float>(Point(yi - 1, xi - 1)) - middle.at<float>(Point(yi - 1, xi + 1)) - middle.at<float>(Point(yi + 1, xi - 1))) / 4.0;
 
-						dxy = (cvGetReal2D(middle, yi-1, xi-1) +
-							  cvGetReal2D(middle, yi+1, xi+1) -
-							  cvGetReal2D(middle, yi+1, xi-1) - 
-							  cvGetReal2D(middle, yi-1, xi+1)) / 4.0;
 
 						trH = dxx + dyy;
-						detH = dxx*dyy - dxy*dxy;
+						detH = dxx * dyy - dxy * dxy;
 
-						curvature_ratio = trH*trH/detH;
+						curvature_ratio = trH * trH / detH;
 						//printf("Threshold: %f - Ratio: %f\n", curvature_threshold, curvature_ratio);
-						if(detH<0 || curvature_ratio>curvature_threshold)
+						if (detH<0 || curvature_ratio>curvature_threshold)
 						{
-							cvSetReal2D(m_extrema[i][j-1], yi, xi, 0);
 							num--;
 							numRemoved++;
-
-							justSet=false;
+							m_extrema[i][j - 1].at<uchar>(Point(yi, xi)) = 0;
+							justSet = false;
 						}
 					}
 				}
 			}
 
-			// Save the image
-			/*char* filename = new char[200];
-			sprintf(filename, "C:\\SIFT Test\\Extrema\\extrema_oct_%d_scale_%d.jpg", i, j-1);
-			cvSaveImage(filename, m_extrema[i][j-1]);*/
+			//imshow(to_string(10*j+i),m_extrema[i][j-1]);
+			//waitKey(0);
 		}
 	}
 
 	m_numKeypoints = num;
-	printf("Found %d keypoints\n", num);
-	printf("Rejected %d keypoints\n", numRemoved);
+	//printf("Found %d keypoints\n", num);
+	//printf("Rejected %d keypoints\n", numRemoved);
+}
+
+// GetKernelSize()
+// Returns the size of the kernal for the Gaussian blur given the sigma and
+// cutoff value.
+int meSIFT::GetKernelSize(double sigma, double cut_off)
+{
+	int i;
+	for (i = 0; i < MAX_KERNEL_SIZE; i++)
+		if (exp(-((double)(i * i)) / (2.0 * sigma * sigma)) < cut_off)
+			break;
+	int size = 2 * i - 1;
+	return size;
 }
 
+
 // AssignOrientations()
 // For all the key points, generate an orientation.
-void SIFT::AssignOrientations()
+
+void meSIFT::AssignOrientations()
 {
-	printf("Assigning orientations...\n");
-	unsigned int i, j, k, xi, yi;
+	int i, j, k, xi, yi;
 	int kk, tt;
 
-	// These images hold the magnitude and direction of gradient 
-	// for all blurred out images
-	IplImage*** magnitude = new IplImage**[m_numOctaves];
-	IplImage*** orientation = new IplImage**[m_numOctaves];
+	//the 2-D vectors for Magnitude and orientation
+	vector<vector<Mat> >magnitude;
+	vector<vector<Mat> >orientation;
 
-	// Allocate some memory
-	for(i=0;i<m_numOctaves;i++)
+	//iterate through all the octaves
+	for (i = 0; i < m_numOctaves; i++)
 	{
-		magnitude[i] = new IplImage*[m_numIntervals];
-		orientation[i] = new IplImage*[m_numIntervals];
-	}	
-
-	// These two loops are to calculate the magnitude and orientation of gradients
-	// through all octaces once and for all. We don't run around calculating things
-	// again and again that way.
-
-	// Iterate through all octaves
-	for(i=0;i<m_numOctaves;i++)
-	{
-		// Iterate through all scales
-		for(j=1;j<m_numIntervals+1;j++)
+		magnitude.push_back(vector<Mat>());
+		orientation.push_back(vector<Mat>());
+		for (j = 1; j < m_numIntervals + 1; j++)
 		{
-			magnitude[i][j-1] = cvCreateImage(cvGetSize(m_gList[i][j]), 32, 1);
-			orientation[i][j-1] = cvCreateImage(cvGetSize(m_gList[i][j]), 32, 1);
-
-			cvZero(magnitude[i][j-1]);
-			cvZero(orientation[i][j-1]);
+			Mat temp1 = Mat::zeros(m_glist[i][j].size(),CV_32FC1);
+			Mat temp2 = Mat::zeros(m_glist[i][j].size(), CV_32FC1);
+			magnitude[i].push_back(temp1);
+			orientation[i].push_back(temp2);
 
 			// Iterate over the gaussian image with the current octave and interval
-			for(xi=1;xi<m_gList[i][j]->width-1;xi++)
+			for (xi = 1; xi < m_glist[i][j].rows-1; xi++)
 			{
-				for(yi=1;yi<m_gList[i][j]->height-1;yi++)
+				for (yi = 1; yi < m_glist[i][j].cols-1; yi++)
 				{
-					// Calculate gradient
-					double dx = cvGetReal2D(m_gList[i][j], yi, xi+1) - cvGetReal2D(m_gList[i][j], yi, xi-1);
-					double dy = cvGetReal2D(m_gList[i][j], yi+1, xi) - cvGetReal2D(m_gList[i][j], yi-1, xi);
+					float dx = m_glist[i][j].at<float>(Point(yi, xi + 1)) - m_glist[i][j].at<float>(Point(yi, xi - 1));
+					float dy=  m_glist[i][j].at<float>(Point(yi+1, xi)) - m_glist[i][j].at<float>(Point(yi-1, xi));
 
-					// Store magnitude
-					cvSetReal2D(magnitude[i][j-1], yi, xi, sqrt(dx*dx + dy*dy));
+					//storing the magnitude in the new magnitude grid
+					magnitude[i][j - 1].at<float>(Point(yi, xi)) = sqrt(dx*dx+dy*dy);
 
-					// Store orientation as radians
-					double ori=atan(dy/dx);					
-					cvSet2D(orientation[i][j-1], yi, xi, cvScalar(ori));
+					//calculating and storing the orientation
+					float ori = atan(dy / dx);
+					orientation[i][j - 1].at<float>(Point(yi, xi)) = ori;
 				}
 			}
-
-			// Save these images for fun
-			/*char* filename = new char[200];
-			sprintf(filename, "C:\\SIFT Test\\Mag\\mag_oct_%d_scl_%d.jpg", i, j-1);
-			cvSaveImage(filename, magnitude[i][j-1]);
-
-			sprintf(filename, "C:\\SIFT Test\\Ori\\ori_oct_%d_scl_%d.jpg", i, j-1);
-			cvSaveImage(filename, orientation[i][j-1]);*/
+			//imshow(to_string(10 * j + i), magnitude[i][j-1]);
+			//waitKey(0);
+			//imshow(to_string(10 * j + i), orientation[i][j-1]);
+			//waitKey(0);
 		}
 	}
 
 	// The histogram with 8 bins
-	double* hist_orient = new double[NUM_BINS];
-
-	// Go through all octaves
-	for(i=0;i<m_numOctaves;i++)
+	float* hist_orient = new float[NUM_BINS];
+	//iterate thru all the octaves
+	for (int i = 0; i < m_numOctaves; i++)
 	{
-		// Store current scale, width and height
-		unsigned int scale = (unsigned int)pow(2.0, (double)i);
-		unsigned int width = m_gList[i][0]->width;
-		unsigned int height= m_gList[i][0]->height;
-
-		// Go through all intervals in the current scale
-		for(j=1;j<m_numIntervals+1;j++)
+		int scale = int(pow(2.0,double(i)));
+		int width = m_glist[i][0].cols;
+		int height = m_glist[i][0].rows;
+		//iterate thru all the intervals in the current scale
+		for (int j = 1; j < m_numIntervals + 1; j++)
 		{
 			double abs_sigma = m_absSigma[i][j];
-
-			// This is used for magnitudes
-			IplImage* imgWeight = cvCreateImage(cvSize(width, height), 32, 1);
-			cvSmooth(magnitude[i][j-1], imgWeight, CV_GAUSSIAN, 0, 0, 1.5*abs_sigma);
+			Mat img_weight=Mat::zeros(magnitude[i][j-1].size(),CV_32FC1);
+			GaussianBlur(magnitude[i][j-1], img_weight, Size(0, 0), 1.5 * abs_sigma);
 
 			// Get the kernel size for the Guassian blur
-			int hfsz = GetKernelSize(1.5*abs_sigma)/2;
-
-			// Temporarily used to generate a mask of region used to calculate 
-			// the orientations
-			IplImage* imgMask = cvCreateImage(cvSize(width, height), 8, 1);
-			cvZero(imgMask);
-
-			// Iterate through all points at this octave and interval
-			for(xi=0;xi<width;xi++)
+			int h = meSIFT::GetKernelSize(1.5 * abs_sigma) / 2;
+			//temporarily creating a mask of region
+			//to calculate the orientation
+			Mat img_mask = Mat::zeros(Size(width,height),CV_8UC1);
+
+			//iterate thru all the points in the octave 
+			//and the interval
+			for (xi = 1; xi < height-1; xi++)
 			{
-				for(yi=0;yi<height;yi++)
+				for (yi = 1; yi < width-1; yi++)
 				{
-					// We're at a keypoint
-					if(cvGetReal2D(m_extrema[i][j-1], yi, xi)!=0)
+					int val = (int)m_extrema[i][j - 1].at<uchar>(Point(yi, xi));
+					if (val!=0)
 					{
 						// Reset the histogram thingy
-						for(k=0;k<NUM_BINS;k++)
-							hist_orient[k]=0.0;
-
+						for (k = 0; k < NUM_BINS; k++)
+							hist_orient[k] = 0.0;
 						// Go through all pixels in the window around the extrema
-						for(kk=-hfsz;kk<=hfsz;kk++)
+						for (kk = -h; kk <= h; kk++)
 						{
-							for(tt=-hfsz;tt<=hfsz;tt++)
+							for (tt = -h; tt <= h; tt++)
 							{
 								// Ensure we're within the image
-								if(xi+kk<0 || xi+kk>=width || yi+tt<0 || yi+tt>=height)
+								if (xi + kk < 0 || xi + kk >= height || yi + tt < 0 || yi + tt >= width)
 									continue;
 
- 								double sampleOrient = cvGetReal2D(orientation[i][j-1], yi+tt, xi+kk);
+								float sample_orient = orientation[i][j-1].at<float>(Point(yi + tt, xi + kk));
 
-								if(sampleOrient <=-M_PI || sampleOrient>M_PI)
-									printf("Bad Orientation: %f\n", sampleOrient);
-								
-								sampleOrient+=M_PI;
+								//if (sampleOrient <= -M_PI || sampleOrient > M_PI)
+									//printf("Bad Orientation: %f\n", sampleOrient);      //CHANGEABLE
+
+								sample_orient+=(float)M_PI;
+
+								int sample_orient_degrees = int(sample_orient * 180 / M_PI);
+								hist_orient[(int)sample_orient_degrees / (360 / NUM_BINS)] += img_weight.at<float>(Point(yi + tt, xi + kk));
+								img_mask.at<uchar>(Point(yi + tt, xi + kk))=255;
 
-								// Convert to degrees
-								unsigned int sampleOrientDegrees = sampleOrient * 180 / M_PI;
-								hist_orient[(int)sampleOrientDegrees / (360/NUM_BINS)] += cvGetReal2D(imgWeight, yi+tt, xi+kk);
-								cvSetReal2D(imgMask, yi+tt, xi+kk, 255);
 							}
 						}
-
+						
 						// We've computed the histogram. Now check for the maximum
-						double max_peak = hist_orient[0];
-						unsigned int max_peak_index = 0;
-						for(k=1;k<NUM_BINS;k++)
+						float max_peak = hist_orient[0];
+						int max_peak_index = 0;
+						for (k = 1; k < NUM_BINS; k++)
 						{
-							if(hist_orient[k]>max_peak)
+							if (hist_orient[k] > max_peak)
 							{
-								max_peak=hist_orient[k];
+								max_peak = hist_orient[k];
 								max_peak_index = k;
 							}
 						}
 
-						// List of magnitudes and orientations at the current extrema
-						vector<double> orien;
-						vector<double> mag;
-						for(k=0;k<NUM_BINS;k++)
+						//list of magnitudes and orientation at the current extrema
+						vector<float> orien;
+						vector<float> mag;
+
+						for (int k = 0; k < NUM_BINS; k++)
 						{
 							// Do we have a good peak?
-							if(hist_orient[k]> 0.8*max_peak)
+							if (hist_orient[k] > 0.8 * max_peak)
 							{
 								// Three points. (x2,y2) is the peak and (x1,y1)
 								// and (x3,y3) are the neigbours to the left and right.
 								// If the peak occurs at the extreme left, the "left
 								// neighbour" is equal to the right most. Similarly for
 								// the other case (peak is rightmost)
-								double x1 = k-1;
-								double y1;
-								double x2 = k;
-								double y2 = hist_orient[k];
-								double x3 = k+1;
-								double y3;
-
-								if(k==0)
+								float x1 = (float)(k - 1);
+								float y1;
+								float x2 = (float)k;
+								float y2 = hist_orient[k];
+								float x3 = (float)(k + 1);
+								float y3;
+								//cout << x1 << " " << x2 << " " << x3 << endl;
+								if  (k ==0)
 								{
-									y1 = hist_orient[NUM_BINS-1];
+									y1 = hist_orient[NUM_BINS - 1];
 									y3 = hist_orient[1];
 								}
-								else if(k==NUM_BINS-1)
+								else if (k == NUM_BINS - 1)
 								{
-									y1 = hist_orient[NUM_BINS-1];
+									y1 = hist_orient[NUM_BINS - 2];             ///CHANGED
 									y3 = hist_orient[0];
 								}
 								else
 								{
-									y1 = hist_orient[k-1];
-									y3 = hist_orient[k+1];
+									y1 = hist_orient[k - 1];
+									y3 = hist_orient[k + 1];
 								}
 
+
+
 								// Next we fit a downward parabola aound
 								// these three points for better accuracy
 
 								// A downward parabola has the general form
 								//
-                                // y = a * x^2 + bx + c
-                                // Now the three equations stem from the three points
-                                // (x1,y1) (x2,y2) (x3.y3) are
+								// y = a * x^2 + bx + c
+								// Now the three equations stem from the three points
+								// (x1,y1) (x2,y2) (x3.y3) are
 								//
-                                // y1 = a * x1^2 + b * x1 + c
-                                // y2 = a * x2^2 + b * x2 + c
-                                // y3 = a * x3^2 + b * x3 + c
+								// y1 = a * x1^2 + b * x1 + c
+								// y2 = a * x2^2 + b * x2 + c
+								// y3 = a * x3^2 + b * x3 + c
 								//
-                                // in Matrix notation, this is y = Xb, where
-                                // y = (y1 y2 y3)' b = (a b c)' and
-                                // 
-                                //     x1^2 x1 1
-                                // X = x2^2 x2 1
-                                //     x3^2 x3 1
-                                //
-                                // OK, we need to solve this equation for b
-                                // this is done by inverse the matrix X
-                                //
-                                // b = inv(X) Y
-
-								double *b = new double[3];
-								CvMat *X = cvCreateMat(3, 3, CV_32FC1);
-								CvMat *matInv = cvCreateMat(3, 3, CV_32FC1);
-
-								cvSetReal2D(X, 0, 0, x1*x1);
-								cvSetReal2D(X, 1, 0, x1);
-								cvSetReal2D(X, 2, 0, 1);
-
-								cvSetReal2D(X, 0, 1, x2*x2);
-								cvSetReal2D(X, 1, 1, x2);
-								cvSetReal2D(X, 2, 1, 1);
-
-								cvSetReal2D(X, 0, 2, x3*x3);
-								cvSetReal2D(X, 1, 2, x3);
-								cvSetReal2D(X, 2, 2, 1);
-
-								// Invert the matrix
-								cvInv(X, matInv);
-
-								b[0] = cvGetReal2D(matInv, 0, 0)*y1 + cvGetReal2D(matInv, 1, 0)*y2 + cvGetReal2D(matInv, 2, 0)*y3;
-								b[1] = cvGetReal2D(matInv, 0, 1)*y1 + cvGetReal2D(matInv, 1, 1)*y2 + cvGetReal2D(matInv, 2, 1)*y3;
-								b[2] = cvGetReal2D(matInv, 0, 2)*y1 + cvGetReal2D(matInv, 1, 2)*y2 + cvGetReal2D(matInv, 2, 2)*y3;
-
-								double x0 = -b[1]/(2*b[0]);
-
+								// in Matrix notation, this is y = Xb, where
+								// y = (y1 y2 y3)' b = (a b c)' and
+								// 
+								//     x1^2 x1 1
+								// X = x2^2 x2 1
+								//     x3^2 x3 1
+								//
+								// OK, we need to solve this equation for b
+								// this is done by inverse the matrix X
+								//
+								// b = inv(X) Y
+
+								float* b = new float[3];
+								Mat X = Mat(Size(3,3),CV_32FC1);
+								Mat matInv = Mat(Size(3,3),CV_32FC1);
+								X.at<float>(Point(0,0))=x1*x1;
+								X.at<float>(Point(1,0))=x1;
+								X.at<float>(Point(2,0))=1;
+								X.at<float>(Point(0,1))=x2*x2;
+								X.at<float>(Point(1,1))=x2;
+								X.at<float>(Point(2,1))=1;
+								X.at<float>(Point(0,2))=x3*x3;
+								X.at<float>(Point(1,2))=x3;
+								X.at<float>(Point(2, 2)) =1;
+
+								//finding and storing the inverse
+								invert(X,matInv);
+								/*for (int i = 0; i < 3; i++)
+								{
+									for (int j = 0; j < 3; j++)
+									{
+										cout << matInv.at<float>(Point(j,i)) << " ";
+									}
+									cout << endl;
+								}*/
+								
+								b[0] = matInv.at<float>(Point(0,0))*y1+matInv.at<float>(Point(1,0))*y2+matInv.at<float>(Point(2,0))*y3;
+								b[1] = matInv.at<float>(Point(0,1))*y1+matInv.at<float>(Point(1,1))*y2+matInv.at<float>(Point(2,1))*y3;
+								b[2] = matInv.at<float>(Point(0, 2))*y1 + matInv.at<float>(Point(1, 2))*y2 + matInv.at<float>(Point(2, 2))*y3;
+								
+								float x0 = -b[1] / (2 * b[0]);
+								//
+								//cout << b[0] << endl;
 								// Anomalous situation
-								if(fabs(x0)>2*NUM_BINS)
-									x0=x2;
-
-								while(x0<0)
+								if (fabs(x0) > 2 * NUM_BINS)
+									x0 = x2;
+								while (x0 < 0)
 									x0 += NUM_BINS;
-								while(x0>= NUM_BINS)
-									x0-= NUM_BINS;
+								while (x0 >= NUM_BINS)
+									x0 -= NUM_BINS;
 
 								// Normalize it
-								double x0_n = x0*(2*M_PI/NUM_BINS);
-
-								assert(x0_n>=0 && x0_n<2*M_PI);
-								x0_n -= M_PI;
-								assert(x0_n>=-M_PI && x0_n<M_PI);
+								float s = (float)M_PI / NUM_BINS;
+								float x0_n = x0 * (2 * s);
+								//cout << x0_n << endl;
+								assert(x0_n >= 0 && x0_n < 2 * M_PI);
+								x0_n -= (float)M_PI;
+								assert(x0_n >= -M_PI && x0_n < M_PI);
 
 								orien.push_back(x0_n);
+								//cout << orien.back() << " ";
 								mag.push_back(hist_orient[k]);
+
+
+							
 							}
 						}
-
+						//cout << endl;
 						// Save this keypoint into the list
-						m_keyPoints.push_back(Keypoint(xi*scale/2, yi*scale/2, mag, orien, i*m_numIntervals+j-1));
+						float r = float(xi * scale / 2);
+						float s = float(yi * scale / 2);
+						//cout << orien[0] << endl;
+						//if (orien.empty()) cout << "EMpty!" << endl;
+						Keypoint k(r, s, mag, orien, i* m_numIntervals + j - 1);
+						m_keyPoints.push_back(k);
 					}
 				}
 			}
 
-			// Save the regions!
-			/*char* filename = new char[200];
-			sprintf(filename, "C:\\SIFT Test\\Orientation Region\\ori_region_oct_%d_scl_%d.jpg", i, j-1);
-			cvSaveImage(filename, imgMask);*/
-			cvReleaseImage(&imgMask);
+			//imshow(to_string(10 * j + i), img_mask);
+			//waitKey(0);
 		}
 	}
 
 	// Finally, we're done with all the magnitude and orientation images.
-	// Erase them from RAM
-	assert(m_keyPoints.size() == m_numKeypoints);
-	for(i=0;i<m_numOctaves;i++)
-	{
-		for(j=1;j<m_numIntervals+1;j++)
-		{
-			cvReleaseImage(&magnitude[i][j-1]);
-			cvReleaseImage(&orientation[i][j-1]);
-		}
-
-		delete [] magnitude[i];
-		delete [] orientation[i];
-	}
+	//assert(m_keyPoints.size() == m_numKeypoints);
 
-	delete [] magnitude;
-	delete [] orientation;
 }
 
-// ExtractKeypointDescriptors()
-// Generates a unique descriptor for each keypoint descriptor
-void SIFT::ExtractKeypointDescriptors()
+
+
+// BuildInterpolatedGaussianTable()
+// This function actually generates the bell curve like image for the weighted
+// addition earlier.
+Mat meSIFT::BuildInterpolatedGaussianTable(int size, double sigma)
 {
-	printf("Extract keypoint descriptors...\n");
+	int i, j;
+	double half_kernel_size = size / 2 - 0.5;
+
+	double sog = 0;
+	Mat ret=Mat::zeros(size,size,CV_32FC1);
 
-	// For loops
-	unsigned int i, j;
+	assert(size % 2 == 0);
 
-	// Interpolated thingy. We're dealing with "inbetween" gradient
-	// magnitudes and orientations
-	IplImage*** imgInterpolatedMagnitude = new IplImage**[m_numOctaves];
-	IplImage*** imgInterpolatedOrientation = new IplImage**[m_numOctaves];
-	for(i=0;i<m_numOctaves;i++)
+	double temp = 0;
+	for (i = 0; i < size; i++)
 	{
-		imgInterpolatedMagnitude[i] = new IplImage*[m_numIntervals];
-		imgInterpolatedOrientation[i] = new IplImage*[m_numIntervals];
+		for (j = 0; j < size; j++)
+		{
+			temp = meSIFT::gaussian2D(i - half_kernel_size, j - half_kernel_size, sigma);
+			ret.at<float>(Point(j, i)) = (float)temp;
+			sog += (float)temp;
+		}
 	}
 
-	// These two loops calculate the interpolated thingy for all octaves
-	// and subimages
-	for(i=0;i<m_numOctaves;i++)
+	for (i = 0; i < size; i++)
 	{
-		for(j=1;j<m_numIntervals+1;j++)
+		for (j = 0; j < size; j++)
 		{
-			unsigned int width = m_gList[i][j]->width;
-			unsigned int height =m_gList[i][j]->height;
+			float val = ret.at<float>(Point(j, i));
+			ret.at<float>(Point(j, i)) = (float)(1.0 / (sog * val));
+		}
+	}
 
-			// Create an image and zero it out.
-			IplImage* imgTemp = cvCreateImage(cvSize(width*2, height*2), 32, 1);
-			cvZero(imgTemp); 
+	return ret;
+}
 
-			// Scale it up. This will give us "access" to in betweens
-			cvPyrUp(m_gList[i][j], imgTemp);
+// gaussian2D
+// Returns the value of the bell curve at a (x,y) for a given sigma
+double meSIFT::gaussian2D(double x, double y, double sigma)
+{
+	double ret = 1.0 / (2 * M_PI * sigma * sigma) * exp(-(x * x + y * y) / (2.0 * sigma * sigma));
 
-			// Allocate memory and zero them
-			imgInterpolatedMagnitude[i][j-1] = cvCreateImage(cvSize(width+1, height+1), 32, 1);
-			imgInterpolatedOrientation[i][j-1] = cvCreateImage(cvSize(width+1, height+1), 32, 1);
-			cvZero(imgInterpolatedMagnitude[i][j-1]);
-			cvZero(imgInterpolatedOrientation[i][j-1]);
 
+	return ret;
+}
+
+// ExtractKeypointDescriptors()
+// Generates a unique descriptor for each keypoint descriptor
+/*
+void meSIFT::ExtractKeypointDescriptors()
+{
+	vector<vector<Mat> >imginterpolatedMagnitude;
+	vector<vector<Mat> >imginterpolatedOrientation;
+	for (int i = 0; i < m_numOctaves; i++)
+	{
+		imginterpolatedMagnitude.push_back(vector<Mat>());
+		imginterpolatedOrientation.push_back(vector<Mat>());
+		for (int j = 1; j < m_numIntervals + 1; j++)
+		{
+			Mat temp = Mat::zeros(m_glist[i][j].size()*2,CV_32FC1);
+			// Scale it up. This will give us "access" to in betweens
+			pyrUp(m_glist[i][j], temp);
+			int width = m_glist[i][j].cols;
+			int height = m_glist[i][j].rows;
+			Mat temp1 = Mat::zeros(Size(width+1,height+1), CV_32FC1);
+			Mat temp2 = Mat::zeros(Size(width + 1, height + 1), CV_32FC1);
+			imginterpolatedMagnitude[i].push_back(temp1);
+			imginterpolatedOrientation[i].push_back(temp2);
 			// Do the calculations
-			for(float ii=1.5;ii<width-1.5;ii++)
+			for (int ii = 1; ii < height - 1; ii++)
 			{
-				for(float jj=1.5;jj<height-1.5;jj++)
+				for (int jj = 1; jj < width - 1; jj++)
 				{
 					// "inbetween" change
-					double dx = (cvGetReal2D(m_gList[i][j], jj, ii+1.5) + cvGetReal2D(m_gList[i][j], jj, ii+0.5))/2 - (cvGetReal2D(m_gList[i][j], jj, ii-1.5) + cvGetReal2D(m_gList[i][j], jj, ii-0.5))/2;
-					double dy = (cvGetReal2D(m_gList[i][j], jj+1.5, ii) + cvGetReal2D(m_gList[i][j], jj+0.5, ii))/2 - (cvGetReal2D(m_gList[i][j], jj-1.5, ii) + cvGetReal2D(m_gList[i][j], jj-0.5, ii))/2;
-
-					unsigned int iii = ii+1;
-					unsigned int jjj = jj+1;
-					assert(iii<=width && jjj<=height);
-
-					// Set the magnitude and orientation
-					cvSetReal2D(imgInterpolatedMagnitude[i][j-1], jjj, iii, sqrt(dx*dx + dy*dy));
-					cvSetReal2D(imgInterpolatedOrientation[i][j-1], jjj, iii, (atan2(dy,dx)==M_PI)? -M_PI:atan2(dy,dx) );
+					float dx = (float)(m_glist[i][j].at<float>(Point(jj,ii+1))+m_glist[i][j].at<float>(Point(jj,ii)))/2.0-(float)(m_glist[i][j].at<float>(Point(jj,ii))+m_glist[i][j].at<float>(Point(jj,ii-1)))/2.0;
+					float dy = (float)(m_glist[i][j].at<float>(Point(jj,ii+1))+m_glist[i][j].at<float>(Point(jj,ii)))/2.0-(float)(m_glist[i][j].at<float>(Point(jj,ii))+m_glist[i][j].at<float>(Point(jj,ii-1)))/2.0;
+					int iii = ii + 1;
+					int jjj = jj + 1;
+					if (iii <= height && jjj <= width)
+					{
+						imginterpolatedMagnitude[i][j - 1].at<float>(Point(jjj, iii)) = sqrt(dx*dx+dy*dy);
+						imginterpolatedOrientation[i][j-1].at<float>(Point(jjj,iii))= (atan2(dy, dx) == M_PI) ? float(-M_PI) : float(atan2(dy, dx));
+					}
 				}
 			}
-
 			// Pad the edges with zeros
-			for(unsigned int iii=0;iii<width+1;iii++)
+			for (int iii = 0; iii < height+1; iii++)
 			{
-				cvSetReal2D(imgInterpolatedMagnitude[i][j-1], 0, iii, 0);
-				cvSetReal2D(imgInterpolatedMagnitude[i][j-1], height, iii, 0);
-				cvSetReal2D(imgInterpolatedOrientation[i][j-1], 0, iii, 0);
-				cvSetReal2D(imgInterpolatedOrientation[i][j-1], height, iii, 0);
+				imginterpolatedMagnitude[i][j-1].at<float>(Point(0,iii))= 0.0;
+				imginterpolatedMagnitude[i][j-1].at<float>(Point(width,iii))= 0.0;
+				imginterpolatedOrientation[i][j-1].at<float>(Point(0,iii))= 0.0;
+				imginterpolatedOrientation[i][j-1].at<float>(Point(width,iii))= 0.0;
 			}
-
-			for(unsigned int jjj=0;jjj<height+1;jjj++)
+			for (int jjj = 0; jjj < width+1; jjj++)
 			{
-				cvSetReal2D(imgInterpolatedMagnitude[i][j-1], jjj, 0, 0);
-				cvSetReal2D(imgInterpolatedMagnitude[i][j-1], jjj, width, 0);
-				cvSetReal2D(imgInterpolatedOrientation[i][j-1], jjj, 0, 0);
-				cvSetReal2D(imgInterpolatedOrientation[i][j-1], jjj, width, 0);
+				imginterpolatedMagnitude[i][j - 1].at<float>(Point(jjj, 0)) = 0.0;
+				imginterpolatedMagnitude[i][j - 1].at<float>(Point(jjj,height)) = 0.0;
+				imginterpolatedOrientation[i][j - 1].at<float>(Point(jjj,0)) = 0.0;
+				imginterpolatedOrientation[i][j - 1].at<float>(Point(jjj,height)) = 0.0;
 			}
-
-			// Now we have the imgInterpolated* ready. Store and get started
-			/*char* filename = new char[200];
-			sprintf(filename, "C:\\SIFT Test\\Interpolated Mag\\intmag_oct_%d_scl_%d.jpg", i, j-1);
-			cvSaveImage(filename, imgInterpolatedMagnitude[i][j-1]);
-
-			sprintf(filename, "C:\\SIFT Test\\Interpolated Ori\\intori_oct_%d_scl_%d.jpg", i, j-1);
-			cvSaveImage(filename, imgInterpolatedOrientation[i][j-1]);*/
-
-			cvReleaseImage(&imgTemp);
-
+			//imshow(to_string(10 * j + i),imginterpolatedMagnitude[i][j-1]);
+			//waitKey(0);
+			//imshow(to_string(10 * j + i),imginterpolatedOrientation[i][j-1]);
+			//waitKey(0);
 		}
 	}
-
 	// Build an Interpolated Gaussian Table of size FEATURE_WINDOW_SIZE
 	// Lowe suggests sigma should be half the window size
 	// This is used to construct the "circular gaussian window" to weight 
 	// magnitudes
-	CvMat *G = BuildInterpolatedGaussianTable(FEATURE_WINDOW_SIZE, 0.5*FEATURE_WINDOW_SIZE);
-	
-	vector<double> hist(DESC_NUM_BINS);
-
-	// Loop over all keypoints
-	for(unsigned int ikp = 0;ikp<m_numKeypoints;ikp++)
+	Mat G = meSIFT::BuildInterpolatedGaussianTable(FEATURE_WINDOW_SIZE, 0.5 * FEATURE_WINDOW_SIZE);
+	vector<float> hist(DESC_NUM_BINS);
+	for (int ikp = 0; ikp < m_numKeypoints; ikp++)
 	{
-		unsigned int scale = m_keyPoints[ikp].scale;
+		int scale = m_keyPoints[ikp].scale;
 		float kpxi = m_keyPoints[ikp].xi;
 		float kpyi = m_keyPoints[ikp].yi;
-
 		float descxi = kpxi;
 		float descyi = kpyi;
-
-		unsigned int ii = (unsigned int)(kpxi*2) / (unsigned int)(pow(2.0, (double)scale/m_numIntervals));
-		unsigned int jj = (unsigned int)(kpyi*2) / (unsigned int)(pow(2.0, (double)scale/m_numIntervals));
-
-		unsigned int width = m_gList[scale/m_numIntervals][0]->width;
-		unsigned int height = m_gList[scale/m_numIntervals][0]->height;
-
-		vector<double> orien = m_keyPoints[ikp].orien;
-		vector<double> mag = m_keyPoints[ikp].mag;
-
+		int ii = (int)(kpxi * 2) / (int)(pow(2.0, (double)scale / m_numIntervals));
+		int jj = (int)(kpyi * 2) / (int)(pow(2.0, (double)scale / m_numIntervals));
+		int width = m_glist[scale / m_numIntervals][0].cols;
+		int height = m_glist[scale / m_numIntervals][0].rows;
+		vector<float> orien = m_keyPoints[ikp].orien;
+		vector<float> magn = m_keyPoints[ikp].mag;
+		//cout << orien[2] << endl;
+		//cout << magn[0] << endl;
 		// Find the orientation and magnitude that have the maximum impact
 		// on the feature
-		double main_mag = mag[0];
-		double main_orien = orien[0];
-		for(unsigned int orient_count=1;orient_count<mag.size();orient_count++)
+		float main_magn = INT_MIN;
+		float main_orien = INT_MIN;
+		for (int orient_count = 1; orient_count < magn.size(); orient_count++)
 		{
-			if(mag[orient_count]>main_mag)
+			if (magn[orient_count] > main_magn)
 			{
 				main_orien = orien[orient_count];
-				main_mag = mag[orient_count];
+				main_magn = magn[orient_count];
 			}
 		}
-
-		unsigned int hfsz = FEATURE_WINDOW_SIZE/2;
-		CvMat *weight = cvCreateMat(FEATURE_WINDOW_SIZE, FEATURE_WINDOW_SIZE, CV_32FC1);
-		vector<double> fv(FVSIZE);
-
-		for(i=0;i<FEATURE_WINDOW_SIZE;i++)
+		int hfsz = FEATURE_WINDOW_SIZE / 2;
+		Mat weight = Mat(Size(FEATURE_WINDOW_SIZE, FEATURE_WINDOW_SIZE), CV_32FC1);
+		vector<float> fv(FVSIZE);
+		int i, j;
+		for (i = 0; i < FEATURE_WINDOW_SIZE; i++)
 		{
-			for(j=0;j<FEATURE_WINDOW_SIZE;j++)
+			for (j = 0; j < FEATURE_WINDOW_SIZE; j++)
 			{
-				if(ii+i+1<hfsz || ii+i+1>width+hfsz || jj+j+1<hfsz || jj+j+1>height+hfsz)
-                    cvSetReal2D(weight, j, i, 0);
+				if (ii + i + 1 < hfsz || ii + i + 1 > width + hfsz || jj + j + 1 < hfsz || jj + j + 1 > height + hfsz)
+					weight.at<float>(Point(j, i)) = 0;
 				else
-					cvSetReal2D(weight, j, i, cvGetReal2D(G, j, i)*cvGetReal2D(imgInterpolatedMagnitude[scale/m_numIntervals][scale%m_numIntervals], jj+j+1-hfsz, ii+i+1-hfsz));
+				{
+					float val1 = G.at<float>(Point(j, i));
+					float val2 = imginterpolatedMagnitude[scale / m_numIntervals][scale % m_numIntervals].at<float>(Point(jj + j + 1 - hfsz, ii + i + 1 - hfsz));
+					weight.at<float>(Point(j, i)) = val1 * val2;
+				}
 			}
 		}
-
 		// Now that we've weighted the required magnitudes, we proceed to generating
 		// the feature vector
-
 		// The next two two loops are for splitting the 16x16 window
 		// into sixteen 4x4 blocks
-		for(i=0;i<FEATURE_WINDOW_SIZE/4;i++)			// 4x4 thingy
+		for (i = 0; i < (int)FEATURE_WINDOW_SIZE / 4; i++)			// 4x4 thingy
 		{
-			for(j=0;j<FEATURE_WINDOW_SIZE/4;j++)
+			for (j = 0; j < (int)FEATURE_WINDOW_SIZE / 4; j++)
 			{
 				// Clear the histograms
-				for(unsigned int t=0;t<DESC_NUM_BINS;t++)
-					hist[t]=0.0;
-
+				for (int t = 0; t < DESC_NUM_BINS; t++)
+					hist[t] = 0.0;
 				// Calculate the coordinates of the 4x4 block
-				int starti = (int)ii - (int)hfsz + 1 + (int)(hfsz/2*i);
-				int startj = (int)jj - (int)hfsz + 1 + (int)(hfsz/2*j);
-				int limiti = (int)ii + (int)(hfsz/2)*((int)(i)-1);
-				int limitj = (int)jj + (int)(hfsz/2)*((int)(j)-1);
-
+				int starti = (int)ii - (int)hfsz + 1 + (int)(hfsz / 2 * i);
+				int startj = (int)jj - (int)hfsz + 1 + (int)(hfsz / 2 * j);
+				int limiti = (int)ii + (int)(hfsz / 2) * ((int)(i)-1);
+				int limitj = (int)jj + (int)(hfsz / 2) * ((int)(j)-1);
 				// Go though this 4x4 block and do the thingy :D
-				for(int k=starti;k<=limiti;k++)
+				for (int k = starti; k <= limiti; k++)
 				{
-					for(int t=startj;t<=limitj;t++)
+					for (int t = startj; t <= limitj; t++)
 					{
-						if(k<0 || k>width || t<0 || t>height)
+						if (k<0 || k>width || t<0 || t>height)
 							continue;
-
 						// This is where rotation invariance is done
-						double sample_orien = cvGetReal2D(imgInterpolatedOrientation[scale/m_numIntervals][scale%m_numIntervals], t, k);
+						float sample_orien = imginterpolatedOrientation[scale / m_numIntervals][scale % m_numIntervals].at<float>(Point(t,k));
 						sample_orien -= main_orien;
-
-						while(sample_orien<0)
-							sample_orien+=2*M_PI;
-
-						while(sample_orien>2*M_PI)
-							sample_orien-=2*M_PI;
-
+						if (sample_orien < 0)
+						{
+							while (sample_orien < 0)
+								sample_orien +=(float) 2 * M_PI;
+						}
+							if (sample_orien > 2 * M_PI)
+							{
+								sample_orien = sample_orien - float(2 * M_PI);
+							}
+							int sample_orien_d=0;
+							int bin=0;
+							float bin_f=0.0;
+							float vxls=0.0;
 						// This should never happen
-						if(!(sample_orien>=0 && sample_orien<2*M_PI))
-							printf("BAD: %f\n", sample_orien);
-						assert(sample_orien>=0 && sample_orien<2*M_PI);
-
-						unsigned int sample_orien_d = sample_orien*180/M_PI;
-						assert(sample_orien_d<360);
-
-						unsigned int bin = sample_orien_d/(360/DESC_NUM_BINS);					// The bin
-						double bin_f = (double)sample_orien_d/(double)(360/DESC_NUM_BINS);		// The actual entry
-
-						assert(bin<DESC_NUM_BINS);
-						assert(k+hfsz-1-ii<FEATURE_WINDOW_SIZE && t+hfsz-1-jj<FEATURE_WINDOW_SIZE);
-
-						// Add to the bin
-						hist[bin]+=(1-fabs(bin_f-(bin+0.5))) * cvGetReal2D(weight, t+hfsz-1-jj, k+hfsz-1-ii);
+						//if (!(sample_orien >= 0 && sample_orien < 2 * M_PI))
+							//printf("BAD: %f\n", sample_orien);
+							if (sample_orien >= 0 && sample_orien < 2 * M_PI)
+							{
+								sample_orien_d = sample_orien * 180 / M_PI;
+								if (sample_orien_d < 360)
+								{
+									bin = sample_orien_d / (360 / DESC_NUM_BINS);
+									bin_f= (float)sample_orien_d / (float)(360 / DESC_NUM_BINS);		// The actual entry// The bin
+								}
+							}
+							if (bin < DESC_NUM_BINS && k + hfsz - 1 - ii < FEATURE_WINDOW_SIZE && t + hfsz - 1 - jj < FEATURE_WINDOW_SIZE)
+							{
+								// Add to the bin
+								vxls = weight.at<float>(Point(t + hfsz - 1 - jj, k + hfsz - 1 - ii));
+								hist[bin] += (1 - fabs(bin_f - (bin + 0.5))) * vxls;
+							}
 					}
 				}
-
+				
 				// Keep adding these numbers to the feature vector
-				for(unsigned int t=0;t<DESC_NUM_BINS;t++)
+				for (int t = 0; t < DESC_NUM_BINS; t++)
 				{
-					fv[(i*FEATURE_WINDOW_SIZE/4+j)*DESC_NUM_BINS+t] = hist[t];
+					fv[(i * FEATURE_WINDOW_SIZE / 4 + j) * DESC_NUM_BINS + t] = hist[t];
 				}
 			}
 		}
-
 		// Now, normalize the feature vector to ensure illumination independence
-		double norm=0;
-		for(unsigned int t=0;t<FVSIZE;t++)
-			norm+=pow(fv[t], 2.0);
+		float norm = 0;
+		for (int t = 0; t < FVSIZE; t++)
+			norm +=(float) pow(fv[t], 2.0);
 		norm = sqrt(norm);
-
-		for(unsigned int t=0;t<FVSIZE;t++)
-			fv[t]/=norm;
-
+		for (int t = 0; t < FVSIZE; t++)
+			fv[t] /= norm;
 		// Now, threshold the vector
-		for(unsigned int t=0;t<FVSIZE;t++)
-			if(fv[t]>FV_THRESHOLD)
+		for (int t = 0; t < FVSIZE; t++)
+			if (fv[t] > FV_THRESHOLD)
 				fv[t] = FV_THRESHOLD;
-
 		// Normalize yet again
-		norm=0;
-		for(unsigned int t=0;t<FVSIZE;t++)
-			norm+=pow(fv[t], 2.0);
+		norm = 0;
+		for (int t = 0; t < FVSIZE; t++)
+			norm += pow(fv[t], 2.0);
 		norm = sqrt(norm);
-
-		for(unsigned int t=0;t<FVSIZE;t++)
-			fv[t]/=norm;
-
+		for (int t = 0; t < FVSIZE; t++)
+			fv[t] /= norm;
 		// We're done with this descriptor. Store it into a list
-		m_keyDescs.push_back(Descriptor(descxi, descyi, fv));
-	}
-
-	assert(m_keyDescs.size()==m_numKeypoints);
-
-	// Get rid of memory we don't need anylonger
-	for(i=0;i<m_numOctaves;i++)
-	{
-		for(j=1;j<m_numIntervals+1;j++)
-		{
-			cvReleaseImage(&imgInterpolatedMagnitude[i][j-1]);
-			cvReleaseImage(&imgInterpolatedOrientation[i][j-1]);
-		}
-
-		delete [] imgInterpolatedMagnitude[i];
-		delete [] imgInterpolatedOrientation[i];
-	}
-
-	delete [] imgInterpolatedMagnitude;
-	delete [] imgInterpolatedOrientation;
-}
-
-// GetKernelSize()
-// Returns the size of the kernal for the Gaussian blur given the sigma and
-// cutoff value.
-unsigned int SIFT::GetKernelSize(double sigma, double cut_off)
-{
-    unsigned int i;
-    for (i=0;i<MAX_KERNEL_SIZE;i++)
-        if (exp(-((double)(i*i))/(2.0*sigma*sigma))<cut_off)
-            break;
-    unsigned int size = 2*i-1;
-    return size;
-}
-
-// BuildInterpolatedGaussianTable()
-// This function actually generates the bell curve like image for the weighted
-// addition earlier.
-CvMat* SIFT::BuildInterpolatedGaussianTable(unsigned int size, double sigma)
-{
-	unsigned int i, j;
-	double half_kernel_size = size/2 - 0.5;
-
-	double sog=0;
-	CvMat* ret = cvCreateMat(size, size, CV_32FC1);
-
-	assert(size%2==0);
-
-	double temp=0;
-	for(i=0;i<size;i++)
-	{
-		for(j=0;j<size;j++)
-		{
-			temp = gaussian2D(i-half_kernel_size, j-half_kernel_size, sigma);
-			cvSetReal2D(ret, j, i, temp);
-			sog+=temp;
-		}
-	}
-
-	for(i=0;i<size;i++)
-	{
-		for(j=0;j<size;j++)
-		{
-			cvSetReal2D(ret, j, i, 1.0/sog * cvGetReal2D(ret, j, i));
-		}
+		Descriptor k(descxi, descyi, fv);
+		m_keyDescs.push_back(k);
 	}
-
-	return ret;
-}
-
-// gaussian2D
-// Returns the value of the bell curve at a (x,y) for a given sigma
-double SIFT::gaussian2D(double x, double y, double sigma)
-{
-	double ret = 1.0/(2*M_PI*sigma*sigma) * exp(-(x*x+y*y)/(2.0*sigma*sigma));
-
-	
-	return ret;
+	assert(m_keyDescs.size() == m_numKeypoints);
 }
-
+*/
 // ShowKeypoints()
 // Displays keypoints as calculated by the algorithm
-void SIFT::ShowKeypoints()
+
+void meSIFT::ShowKeypoints()
 {
-	IplImage* img = cvCloneImage(m_srcImage);
+	Mat img = srcImage.clone();
 
-	for(int i=0;i<m_numKeypoints;i++)
+	for (int i = 0; i < m_keyPoints.size(); i++)
 	{
 		Keypoint kp = m_keyPoints[i];
-
-		cvLine(img, cvPoint(kp.xi, kp.yi), cvPoint(kp.xi, kp.yi), CV_RGB(255,255,255), 3);
-		cvLine(img, cvPoint(kp.xi, kp.yi), cvPoint(kp.xi+10*cos(kp.orien[0]), kp.yi+10*sin((double)kp.orien[0])), CV_RGB(255,255,255), 1);
+		circle(img, Point(kp.xi, kp.yi),6,Scalar(0,255,0),1);
+		//line(img, Point(kp.xi, kp.yi), Point(kp.xi, kp.yi), CV_RGB(255, 255, 255), 3);
+		if (kp.orien.empty()) continue;
+		else
+		{
+			line(img, Point(kp.xi, kp.yi), Point(kp.xi+10 *cos(kp.orien[0]), kp.yi + 10 * sin((float)kp.orien[0])), CV_RGB(255, 255, 255), 1);
+		}
 	}
 
-	cvNamedWindow("Keypoints");
-	cvShowImage("Keypoints", img);
+	imshow("Keypoints", img);
+	waitKey(0);
 }
 
 // ShowAbsSigma()
 // This function shows the sigmas used for various images.
-void SIFT::ShowAbsSigma()
+/*void meSIFT::ShowAbsSigma()
 {
-	for(int i=0;i<m_numOctaves;i++)
+	for (int i = 0; i < m_numOctaves; i++)
 	{
-		for(int j=1;j<m_numIntervals+4;j++)
+		for (int j = 1; j < m_numIntervals + 4; j++)
 		{
-			printf("%f\t", m_absSigma[i][j-1]);
+			printf("%f\t", m_absSigma[i][j - 1]);
 		}
 		printf("\n");
 	}
-}
\ No newline at end of file
+}*/
diff --git a/SIFT.h b/SIFT.h
index 86758c6..4c4676a 100644
--- a/SIFT.h
+++ b/SIFT.h
@@ -1,3 +1,6 @@
+#pragma once
+#include <opencv2/opencv.hpp>
+#include <opencv2/highgui.hpp>
 /****************************************************** 
  *	Code by Utkarsh Sinha
  *	Based on JIFT by Jun Liu
@@ -6,48 +9,39 @@
  * Use, reuse, modify, hack, kick. Do whatever you want with
  * this code :)
  ******************************************************/
-
-#include <cv.h>
-#include <highgui.h>
-
+#include <vector>
+#include <opencv2/imgproc.hpp>
 #include "keypoint.h"
 #include "descriptor.h"
-
-class SIFT
+using namespace cv;
+class meSIFT
 {
-public:
-	SIFT(IplImage* img, int octaves, int intervals);
-	SIFT(const char* filename, int octaves, int intervals);
-	~SIFT();
-
-	void DoSift();
-
-	void ShowKeypoints();
-	void ShowAbsSigma();
-
 private:
-	void GenerateLists();
+	vector<vector<Mat> >m_glist;    // A 2D array to hold the different gaussian blurred images 
+	vector<vector<Mat> >m_dogList;   // A 2D array to hold the different DoG images
+	vector<vector<Mat> >m_extrema;    // A 2D array to hold binary images. In the binary image, 1 = extrema, 0 = not extrema
+	vector<vector<double> >m_absSigma;		// A 2D array to hold the sigma used to blur a particular image
+	vector<Keypoint> m_keyPoints;	// Holds each keypoint's basic info
+	vector<Descriptor> m_keyDescs;	// Holds each keypoint's descriptor
+	Mat srcImage;			// The image we're working on
+	int m_numOctaves;		// The desired number of octaves
+	int m_numIntervals;	// The desired number of intervals
+	int m_numKeypoints;	// The number of keypoints detected
+
+	//void GenerateLists();
 	void BuildScaleSpace();
 	void DetectExtrema();
 	void AssignOrientations();
 	void ExtractKeypointDescriptors();
-
-	unsigned int GetKernelSize(double sigma, double cut_off=0.001);
-	CvMat* BuildInterpolatedGaussianTable(unsigned int size, double sigma);
+	int GetKernelSize(double sigma, double cut_off = 0.001);
+	Mat BuildInterpolatedGaussianTable(int size, double sigma);
 	double gaussian2D(double x, double y, double sigma);
 
+public:
+	meSIFT(Mat img, int octaves, int intervals);
 
-private:
-	IplImage* m_srcImage;			// The image we're working on
-	unsigned int m_numOctaves;		// The desired number of octaves
-	unsigned int m_numIntervals;	// The desired number of intervals
-	unsigned int m_numKeypoints;	// The number of keypoints detected
-
-	IplImage***	m_gList;		// A 2D array to hold the different gaussian blurred images
-	IplImage*** m_dogList;		// A 2D array to hold the different DoG images
-	IplImage*** m_extrema;		// A 2D array to hold binary images. In the binary image, 1 = extrema, 0 = not extrema
-	double**	m_absSigma;		// A 2D array to hold the sigma used to blur a particular image
+	void DoSift();
+	void ShowKeypoints();
+	//void ShowAbsSigma();
 
-	vector<Keypoint> m_keyPoints;	// Holds each keypoint's basic info
-	vector<Descriptor> m_keyDescs;	// Holds each keypoint's descriptor
-};
\ No newline at end of file
+};