diff --git a/src/papillon.cpp b/src/papillon.cpp
index 5062328..626d510 100644
--- a/src/papillon.cpp
+++ b/src/papillon.cpp
@@ -17,11 +17,17 @@ class Traite_image {
 		const static int THRESHOLD_DETECT_SENSITIVITY = 10;
 		const static int BLUR_SIZE = 5;
 		const static int THRESHOLD_MOV = 5;
+		const static int MOVEMENT_THRES = 0.1;
 
 
 		Mat prev;
 		Mat last_T;
 		bool first = true;
+
+		// Features vectors
+		vector <Point2f> prev_ftr, cur_ftr;
+
+		// Downsize factor
 		int resize_f = 2;
 
 		int theObject[2] = {0,0};
@@ -72,7 +78,8 @@ class Traite_image {
             Rect myROI(next_stab.size().width/8, next_stab.size().height/8, next_stab.size().width*3/4, next_stab.size().height*3/4);
             Mat next_stab_cropped = next_stab(myROI);
             Mat prev_cropped = prev(myROI);
-			searchForMovementOptFlow(prev_cropped, next_stab_cropped, output);
+            trackingOptFlow(prev_cropped, next_stab_cropped, output);
+			//searchForMovementOptFlow(prev_cropped, next_stab_cropped, output);
 
 
 			pub_img.publish(cv_bridge::CvImage(msg->header, "rgb8", output).toImageMsg());
@@ -109,14 +116,16 @@ class Traite_image {
 			calcOpticalFlowPyrLK(prev_grey, cur_grey, prev_corner, cur_corner, status, err);
 
 			// weed out bad matches
+			prev_ftr.clear();
+			cur_ftr.clear();
 			for(size_t i=0; i < status.size(); i++) {
 				if(status[i]) {
-					prev_corner2.push_back(prev_corner[i]);
-					cur_corner2.push_back(cur_corner[i]);
+					prev_ftr.push_back(prev_corner[i]);
+					cur_ftr.push_back(cur_corner[i]);
 				}
 			}
 
-			Mat T = estimateRigidTransform(prev_corner2, cur_corner2, true); // false = rigid transform, no scaling/shearing
+			Mat T = estimateRigidTransform(prev_ftr, cur_ftr, true); // false = rigid transform, no scaling/shearing
 
 			if(T.data == NULL) {
 				last_T.copyTo(T);
@@ -239,59 +248,40 @@ class Traite_image {
 
 		}
 
-		void trackingOptFlow(Mat prev, Mat cur, Mat &output) {
-			Mat cur_grey, prev_grey;
-			cur.copyTo(output);
-			cvtColor(prev, prev_grey, COLOR_BGR2GRAY);
-			cvtColor(cur, cur_grey, COLOR_BGR2GRAY);
+		void warpPoints(vector<Point2f> p, vector<Point2f> &p_warp, Mat T, bool invert=false) {
+			if(invert)
+				invertAffineTransform(T, T);
 
-			Mat flow;
-			calcOpticalFlowFarneback(prev_grey, cur_grey, flow, 0.5, 3, 15, 3, 5, 1.2, 0);
-			vector<Mat> flow_coord(2);
-			Mat flow_norm, angle;
-			split(flow, flow_coord);
-			cartToPolar(flow_coord[0], flow_coord[1], flow_norm, angle);
+			p_warp.clear();
+			for(size_t i=0; i < p.size(); ++i) {
+				Mat src(3/*rows*/,1 /* cols */,CV_64F);
 
-			//threshold(flow_norm, flow_norm, THRESHOLD_DETECT_SENSITIVITY, 255, THRESH_BINARY);
-			// Blur to eliminate noise
-			blur(flow_norm, flow_norm, Size(BLUR_SIZE, BLUR_SIZE));
-			threshold(flow_norm, flow_norm, THRESHOLD_DETECT_SENSITIVITY, 255, THRESH_BINARY);
-			flow_norm.convertTo(flow_norm, CV_8U);
+				src.at<double>(0,0)=p[i].x;
+				src.at<double>(1,0)=p[i].y;
+				src.at<double>(2,0)=1.0;
 
-			bool objectDetected = false;
-			Mat temp;
-			flow_norm.copyTo(temp);
-			//these two vectors needed for output of findContours
-			vector< vector<Point> > contours;
-			vector<Vec4i> hierarchy;
-			//find contours of filtered image using openCV findContours function
-			//findContours(temp,contours,hierarchy,CV_RETR_CCOMP,CV_CHAIN_APPROX_SIMPLE );// retrieves all contours
-			findContours(temp,contours,hierarchy,CV_RETR_EXTERNAL,CV_CHAIN_APPROX_SIMPLE );// retrieves external contours
-
-			//if contours vector is not empty, we have found some objects
-			if(contours.size()>0)objectDetected=true;
-			else objectDetected = false;
-
-			if(objectDetected){
-				//the largest contour is found at the end of the contours vector
-				//we will simply assume that the biggest contour is the object we are looking for.
-				vector< vector<Point> > largestContourVec;
-				largestContourVec.push_back(contours.at(contours.size()-1));
-				//make a bounding rectangle around the largest contour then find its centroid
-				//this will be the object's final estimated position.
-				objectBoundingRectangle = boundingRect(largestContourVec.at(0));
+				Mat dst = T*src; //USE MATRIX ALGEBRA
+				p_warp.push_back(Point2f(dst.at<double>(0,0),dst.at<double>(1,0)));
 			}
-			//make some temp x and y variables so we dont have to type out so much
-			int x = objectBoundingRectangle.x;
-			int y = objectBoundingRectangle.y;
-			int width = objectBoundingRectangle.width;
-			int height = objectBoundingRectangle.height;
+		}
 
-			//draw a rectangle around the object
-			rectangle(output, Point(x,y), Point(x+width, y+height), Scalar(0, 255, 0), 2);
+		void trackingOptFlow(Mat prev, Mat cur, Mat &output) {
+			vector <Point2f> curc_stab;
+			Mat T = estimateRigidTransform(prev_ftr, cur_ftr, true); // false = rigid transform, no scaling/shearing
+			ROS_INFO("ready to warp");
+			warpPoints(cur_ftr, curc_stab, T, true);
+			ROS_INFO("warped");
 
-			//write the position of the object to the screen
-			putText(output,"Tracking object at (" + intToString(x)+","+intToString(y)+")",Point(x,y),1,1,Scalar(255,0,0),2);
+			vector <Point2f> objects;
+			for(size_t i=0; i < prev_ftr.size(); ++i) {
+				float flow_norm = norm(prev_ftr[i] - cur_ftr[i]) / prev.size().height;
+				if(flow_norm > MOVEMENT_THRES)
+					objects.push_back(cur_ftr[i]);
+			}
+
+			for(size_t i=0; i < objects.size(); ++i) {
+				circle(output, objects[i], 5, Scalar(0, 200, 0), 1);
+			}
 		}
 
 		inline bool isFlowCorrect(Point2f u)