Commits

Leszek committed 3c9bced

* Fixed Haar feature coefficients (this is what I get for sanitising research code at midnight)
* Fixed some gcc compilation issues
* Moved library to separate directory
* Added example command line program
* Added CMake support

  • Participants
  • Parent commits 421a97d

Comments (0)

Files changed (24)

+syntax: glob
+
+*.swp
+*.swo
+tags

File CMakeLists.txt

+project (pupiltracker)
+cmake_minimum_required(VERSION 2.8)
+
+add_subdirectory (lib)
+add_subdirectory (cmd)

File ConicSection.h

-#ifndef ConicSection_h__
-#define ConicSection_h__
-
-template<typename T>
-class ConicSection_
-{
-public:
-	T A,B,C,D,E,F;
-
-	ConicSection_(cv::RotatedRect r)
-	{
-		cv::Point_<T> axis((T)std::cos(CV_PI/180.0 * r.angle), (T)std::sin(CV_PI/180.0 * r.angle));
-		cv::Point_<T> centre(r.center);
-		T a = r.size.width/2;
-		T b = r.size.height/2;
-
-		initFromEllipse(axis, centre, a, b);
-	}
-
-	T algebraicDistance(cv::Point_<T> p)
-	{
-		return A*p.x*p.x + B*p.x*p.y + C*p.y*p.y + D*p.x + E*p.y + F;
-	}
-
-	T distance(cv::Point_<T> p)
-	{
-		//    dist
-		// -----------
-		// |grad|^0.45
-
-		T dist = algebraicDistance(p);
-		cv::Point_<T> grad = algebraicGradient(p);
-
-		T sqgrad = grad.dot(grad);
-
-		return dist / std::pow(sqgrad, T(0.45/2));
-	}
-
-	cv::Point_<T> algebraicGradient(cv::Point_<T> p)
-	{
-		return cv::Point_<T>(2*A*p.x + B*p.y + D, B*p.x + 2*C*p.y + E);
-	}
-
-	cv::Point_<T> algebraicGradientDir(cv::Point_<T> p)
-	{
-		cv::Point_<T> grad = algebraicGradient(p);
-		T len = std::sqrt(grad.ddot(grad));
-		grad.x /= len;
-		grad.y /= len;
-		return grad;
-	}
-
-protected:
-	void initFromEllipse(cv::Point_<T> axis, cv::Point_<T> centre, T a, T b)
-	{
-		T a2 = a * a;
-		T b2 = b * b;
-
-		A = axis.x*axis.x / a2 + axis.y*axis.y / b2;
-		B = 2*axis.x*axis.y / a2 - 2*axis.x*axis.y / b2;
-		C = axis.y*axis.y / a2 + axis.x*axis.x / b2;
-		D = (-2*axis.x*axis.y*centre.y - 2*axis.x*axis.x*centre.x) / a2
-			+ (2*axis.x*axis.y*centre.y - 2*axis.y*axis.y*centre.x) / b2;
-		E = (-2*axis.x*axis.y*centre.x - 2*axis.y*axis.y*centre.y) / a2
-			+ (2*axis.x*axis.y*centre.x - 2*axis.x*axis.x*centre.y) / b2;
-		F = (2*axis.x*axis.y*centre.x*centre.y + axis.x*axis.x*centre.x*centre.x + axis.y*axis.y*centre.y*centre.y) / a2
-			+ (-2*axis.x*axis.y*centre.x*centre.y + axis.y*axis.y*centre.x*centre.x + axis.x*axis.x*centre.y*centre.y) / b2
-			- 1;
-	}
-};
-typedef ConicSection_<float> ConicSection;
-
-#endif // ConicSection_h__

File PupilTracker.cpp

-#include "stdafx.h"
-#include "PupilTracker.h"
-#include "cvx.h"
-
-namespace 
-{
-	struct section_guard
-	{
-		std::string name;
-		tracker_log& log;
-		timer t;
-		section_guard(const std::string& name, tracker_log& log) : name(name), log(log), t() {  }
-		~section_guard() { log.add(name, t); }
-		operator bool() const {return false;}
-	};
-
-	inline section_guard make_section_guard(const std::string& name, tracker_log& log)
-	{
-		return section_guard(name,log);
-	}
-}
-
-#define SECTION(A,B) if (const auto& _section_guard_ = make_section_guard( A , B )) {} else
-
-
-
-
-class HaarSurroundFeature
-{
-public:
-	HaarSurroundFeature(int r1, int r2) : r_inner(r1), r_outer(r2)
-	{
-		//  _________________
-		// |        -ve      |
-		// |     _______     |
-		// |    |   +ve |    |
-		// |    |   .   |    |
-		// |    |_______|    |
-		// |         <r1>    |
-		// |_________<--r2-->|
-
-		// Number of pixels in each part of the kernel
-		int count_inner = r_inner*r_inner;
-		int count_outer = r_outer*r_outer - r_inner*r_inner;
-
-		/*
-		// Frobenius normalized values
-		//
-		// Want norm = 1 where norm = sqrt(sum(pixelvals^2)), so:
-		//  sqrt(count_inner*val_inner^2 + count_outer*val_outer^2) = 1
-		//
-		// Also want sum(pixelvals) = 0, so:
-		//  count_inner*val_inner + count_outer*val_outer = 0
-		//
-		// Solving both of these gives:
-		val_inner = std::sqrt( (double)count_outer/(count_inner*count_outer + sq(count_inner)) );
-		val_outer = -std::sqrt( (double)count_inner/(count_inner*count_outer + sq(count_outer)) );
-		*/
-
-		// Square radius normalised values
-		//
-		// Want the response to be scale-invariant, so scale it by the number of pixels inside it:
-		//  val_inner = 1/count = 1/r_outer^2
-		//
-		// Also want sum(pixelvals) = 0, so:
-		//  count_inner*val_inner + count_outer*val_outer = 0
-		//
-		// Hence:
-		val_inner = 1 / (r_outer*r_outer);
-		val_outer = val_inner*count_inner/count_outer;
-
-	}
-
-	double val_inner, val_outer;
-	int r_inner, r_outer;
-};
-
-cv::RotatedRect fitEllipse(const std::vector<PupilTrackerNative::EdgePoint>& edgePoints)
-{
-	std::vector<cv::Point2f> points;
-	points.reserve(edgePoints.size());
-
-	BOOST_FOREACH(const PupilTrackerNative::EdgePoint& e, edgePoints)
-		points.push_back(e.point);
-
-	return cv::fitEllipse(points);
-}
-
-
-bool PupilTrackerNative::findPupilEllipse(
-	const TrackerParams& params,
-	const cv::Mat& m,
-
-	PupilTrackerNative::findPupilEllipse_out& out,
-	tracker_log& log
-	)
-{
-	// --------------------
-	// Convert to greyscale
-	// --------------------
-
-	cv::Mat_<uchar> mEye;
-
-	SECTION("Grey and crop", log)
-	{
-		// Pick one channel if necessary, and crop it to get rid of borders
-		if (m.channels() == 1)
-		{
-			mEye = m;
-		}
-		else if (m.channels() == 3)
-		{
-			cv::cvtColor(m, mEye, CV_BGR2GRAY);
-		}
-		else if (m.channels() == 4)
-		{
-			cv::cvtColor(m, mEye, CV_BGRA2GRAY);
-		}
-		else
-		{
-			throw std::exception("Unsupported number of channels");
-		}
-	}
-
-	// -----------------------
-	// Find best haar response
-	// -----------------------
-
-	//             _____________________
-	//            |         Haar kernel |
-	//            |                     |
-	//  __________|______________       |
-	// | Image    |      |       |      |
-	// |    ______|______|___.-r-|--2r--|
-	// |   |      |      |___|___|      |
-	// |   |      |          |   |      |
-	// |   |      |          |   |      |
-	// |   |      |__________|___|______|
-	// |   |    Search       |   |
-	// |   |    region       |   |
-	// |   |                 |   |
-	// |   |_________________|   |
-	// |                         |
-	// |_________________________|
-	//
-
-	cv::Mat_<int32_t> mEyeIntegral;
-	int padding = 2*params.Radius_Max;
-
-	SECTION("Integral image", log)
-	{
-		cv::Mat mEyePad;
-		// Need to pad by an additional 1 to get bottom & right edges.
-		cv::copyMakeBorder(mEye, mEyePad, padding, padding, padding, padding, cv::BORDER_REPLICATE);
-		cv::integral(mEyePad, mEyeIntegral);
-	}
-
-	cv::Point2f pHaarPupil;
-	int haarRadius;
-
-	SECTION("Haar responses", log)
-	{
-		const int rstep = 2;
-		const int ystep = 4;
-		const int xstep = 4;
-
-		double minResponse = std::numeric_limits<double>::infinity();
-
-		for (int r = params.Radius_Min; r < params.Radius_Max; r+=rstep)
-		{   
-			// Get Haar feature
-			int r_inner = r;
-			int r_outer = 3*r;
-			HaarSurroundFeature f(r_inner, r_outer);
-
-			// Use TBB for rows
-			std::pair<double,cv::Point2f> minRadiusResponse = tbb::parallel_reduce(
-				tbb::blocked_range<int>(0, (mEye.rows-r - r - 1)/ystep + 1, ((mEye.rows-r - r - 1)/ystep + 1) / 8),
-				std::make_pair(0.0, PupilTrackerNative::UNKNOWN_POSITION),
-				[&] (tbb::blocked_range<int> range, const std::pair<double,cv::Point2f>& minValIn) -> std::pair<double,cv::Point2f>
-			{
-				std::pair<double,cv::Point2f> minValOut = minValIn;
-				for (int i = range.begin(), y = r + range.begin()*ystep; i < range.end(); i++, y += ystep)
-				{
-					//            �         �
-					// row1_outer.|         |  p00._____________________.p01
-					//            |         |     |         Haar kernel |
-					//            |         |     |                     |
-					// row1_inner.|         |     |   p00._______.p01   |
-					//            |-padding-|     |      |       |      |
-					//            |         |     |      | (x,y) |      |
-					// row2_inner.|         |     |      |_______|      |
-					//            |         |     |   p10'       'p11   |
-					//            |         |     |                     |
-					// row2_outer.|         |     |_____________________|
-					//            |         |  p10'                     'p11
-					//            �         �
-
-					int* row1_inner = mEyeIntegral[y+padding - r_inner];
-					int* row2_inner = mEyeIntegral[y+padding + r_inner + 1];
-					int* row1_outer = mEyeIntegral[y+padding - r_outer];
-					int* row2_outer = mEyeIntegral[y+padding + r_outer + 1];
-
-					int* p00_inner = row1_inner + r + padding - r_inner;
-					int* p01_inner = row1_inner + r + padding + r_inner + 1;
-					int* p10_inner = row2_inner + r + padding - r_inner;
-					int* p11_inner = row2_inner + r + padding + r_inner + 1;
-
-					int* p00_outer = row1_outer + r + padding - r_outer;
-					int* p01_outer = row1_outer + r + padding + r_outer + 1;
-					int* p10_outer = row2_outer + r + padding - r_outer;
-					int* p11_outer = row2_outer + r + padding + r_outer + 1;
-
-					for (int x = r; x < mEye.cols - r; x+=xstep)
-					{
-						int sumInner = *p00_inner + *p11_inner - *p01_inner - *p10_inner;
-						int sumOuter = *p00_outer + *p11_outer - *p01_outer - *p10_outer - sumInner;
-
-						double response = f.val_inner * sumInner + f.val_outer * sumOuter;
-
-						if (response < minValOut.first)
-						{
-							minValOut.first = response;
-							minValOut.second = cv::Point(x,y);
-						}
-
-						p00_inner += xstep;
-						p01_inner += xstep;
-						p10_inner += xstep;
-						p11_inner += xstep;
-
-						p00_outer += xstep;
-						p01_outer += xstep;
-						p10_outer += xstep;
-						p11_outer += xstep;
-					}
-				}
-				return minValOut;
-			},
-				[] (const std::pair<double,cv::Point2f>& x, const std::pair<double,cv::Point2f>& y) -> std::pair<double,cv::Point2f>
-			{
-				if (x.first < y.first)
-					return x;
-				else
-					return y;
-			}
-			);
-
-			if (minRadiusResponse.first < minResponse)
-			{
-				minResponse = minRadiusResponse.first;
-				// Set return values
-				pHaarPupil = minRadiusResponse.second;
-				haarRadius = r;
-			}
-			else
-				break;
-		}
-	}
-
-
-	// ---------------------------
-	// Pupil ROI around Haar point
-	// ---------------------------
-	cv::Rect roiHaarPupil = cvx::roiAround(cv::Point(pHaarPupil), haarRadius);
-	cv::Mat_<uchar> mHaarPupil;
-	cvx::getROI(mEye, mHaarPupil, roiHaarPupil);
-
-	out.roiHaarPupil = roiHaarPupil;
-	out.mHaarPupil = mHaarPupil;
-
-	// --------------------------------------------------
-	// Get histogram of pupil region, segment with KMeans
-	// --------------------------------------------------
-
-	const int bins = 256;
-
-	cv::Mat_<float> hist;
-	SECTION("Histogram", log)
-	{
-		int channels[] = {0};
-		int sizes[] = {bins};
-		float range[2] = {0, 256};
-		const float* ranges[] = {range};
-		cv::calcHist(&mHaarPupil, 1, channels, cv::Mat(), hist, 1, sizes, ranges);
-	}
-
-	out.histPupil = hist;
-
-	float threshold;
-	SECTION("KMeans", log)
-	{
-		// Try various candidate centres, return the one with minimal label distance
-		float candidate0[2] = {0, 0};
-		float candidate1[2] = {128, 255};
-		float bestDist = std::numeric_limits<float>::infinity();
-		float bestThreshold = std::numeric_limits<float>::quiet_NaN();
-
-		for (int i = 0; i < 2; i++)
-		{
-			cv::Mat_<uchar> labels;
-			float centres[2] = {candidate0[i], candidate1[i]};
-			float dist = cvx::histKmeans(hist, 0, 256, 2, centres, labels, cv::TermCriteria(cv::TermCriteria::COUNT, 50, 0.0));
-
-			float thisthreshold = (centres[0] + centres[1])/2;
-			if (dist < bestDist && boost::math::isnormal(thisthreshold))
-			{
-				bestDist = dist;
-				bestThreshold = thisthreshold;
-			}
-		}
-
-		if (!boost::math::isnormal(bestThreshold))
-		{
-			// If kmeans gives a degenerate solution, exit early
-			return false;
-		}
-
-		threshold = bestThreshold;
-	}
-
-	cv::Mat_<uchar> mPupilThresh;
-	SECTION("Threshold", log)
-	{
-		cv::threshold(mHaarPupil, mPupilThresh, threshold, 255, cv::THRESH_BINARY_INV);
-	}
-
-	out.threshold = threshold;
-	out.mPupilThresh = mPupilThresh;
-
-	// ---------------------------------------------
-	// Find best region in the segmented pupil image
-	// ---------------------------------------------
-
-	cv::Rect bbPupilThresh;
-	cv::RotatedRect elPupilThresh;
-
-	SECTION("Find best region", log)
-	{
-		cv::Mat_<uchar> mPupilContours = mPupilThresh.clone();
-		std::vector<std::vector<cv::Point> > contours;
-		cv::findContours(mPupilContours, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE);
-
-		if (contours.size() == 0)
-			return false;
-
-		std::vector<cv::Point>& maxContour = contours[0];
-		double maxContourArea = cv::contourArea(maxContour);
-		BOOST_FOREACH(std::vector<cv::Point>& c, contours)
-		{
-			double area = cv::contourArea(c);
-			if (area > maxContourArea)
-			{
-				maxContourArea = area;
-				maxContour = c;
-			}
-		}
-
-		cv::Moments momentsPupilThresh = cv::moments(maxContour);
-
-		bbPupilThresh = cv::boundingRect(maxContour);
-		elPupilThresh = cvx::fitEllipse(momentsPupilThresh);
-
-		// Shift best region into eye coords (instead of pupil region coords), and get ROI
-		bbPupilThresh.x += roiHaarPupil.x;
-		bbPupilThresh.y += roiHaarPupil.y;
-		elPupilThresh.center.x += roiHaarPupil.x;
-		elPupilThresh.center.y += roiHaarPupil.y;
-	}
-
-	out.bbPupilThresh = bbPupilThresh;
-	out.elPupilThresh = elPupilThresh;
-
-	// ------------------------------
-	// Find edges in new pupil region
-	// ------------------------------
-
-	cv::Mat_<uchar> mPupil, mPupilOpened, mPupilBlurred, mPupilEdges;
-	cv::Mat_<float> mPupilSobelX, mPupilSobelY;
-	cv::Rect bbPupil;
-	cv::Rect roiPupil = cvx::roiAround(cv::Point(elPupilThresh.center), haarRadius);
-	SECTION("Pupil preprocessing", log)
-	{
-		cv::Rect roiPadded(roiPupil.x-2, roiPupil.y-2, roiPupil.width+4, roiPupil.height+4);
-		// First get an ROI around the approximate pupil location
-		cvx::getROI(mEye, mPupil, roiPadded, cv::BORDER_REPLICATE);
-
-		cv::Mat morphologyDisk = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(5, 5));
-		cv::morphologyEx(mPupil, mPupilOpened, cv::MORPH_OPEN, morphologyDisk, cv::Point(-1,-1), 2);
-
-		if (params.CannyBlur > 0)
-		{
-			cv::GaussianBlur(mPupilOpened, mPupilBlurred, cv::Size(), params.CannyBlur);
-		}
-		else
-		{
-			mPupilBlurred = mPupilOpened;
-		}
-
-		cv::Sobel(mPupilBlurred, mPupilSobelX, CV_32F, 1, 0, 3);
-		cv::Sobel(mPupilBlurred, mPupilSobelY, CV_32F, 0, 1, 3);
-
-		cv::Canny(mPupilBlurred, mPupilEdges, params.CannyThreshold1, params.CannyThreshold2);
-
-		cv::Rect roiUnpadded(2,2,roiPupil.width,roiPupil.height);
-		mPupil = cv::Mat(mPupil, roiUnpadded);
-		mPupilOpened = cv::Mat(mPupilOpened, roiUnpadded);
-		mPupilBlurred = cv::Mat(mPupilBlurred, roiUnpadded);
-		mPupilSobelX = cv::Mat(mPupilSobelX, roiUnpadded);
-		mPupilSobelY = cv::Mat(mPupilSobelY, roiUnpadded);
-		mPupilEdges = cv::Mat(mPupilEdges, roiUnpadded);
-
-		bbPupil = cvx::boundingBox(mPupil);
-	}
-
-	out.roiPupil = roiPupil;
-	out.mPupil = mPupil;
-	out.mPupilOpened = mPupilOpened;
-	out.mPupilBlurred = mPupilBlurred;
-	out.mPupilSobelX = mPupilSobelX;
-	out.mPupilSobelY = mPupilSobelY;
-	out.mPupilEdges = mPupilEdges;
-
-	// -----------------------------------------------
-	// Get points on edges, optionally using starburst
-	// -----------------------------------------------
-
-	std::vector<cv::Point2f> edgePoints;
-
-	if (params.StarburstPoints > 0)
-	{
-		SECTION("Starburst", log)
-		{
-			// Starburst from initial pupil approximation, stopping when an edge is hit.
-			// Collect all edge points into a vector
-
-			// The initial pupil approximations are:
-			//    Centre of mass of thresholded region
-			//    Halfway along the major axis (calculated form second moments) in each direction
-
-			Concurrency::concurrent_vector<cv::Point2f> edgePointsConcurrent;
-
-			cv::Vec2f elPupil_majorAxis = cvx::majorAxis(elPupilThresh);
-			std::vector<cv::Point2f> centres;
-			centres.push_back(elPupilThresh.center - cv::Point2f(roiPupil.tl()));
-			centres.push_back(elPupilThresh.center - cv::Point2f(roiPupil.tl()) + cv::Point2f(elPupil_majorAxis));
-			centres.push_back(elPupilThresh.center - cv::Point2f(roiPupil.tl()) - cv::Point2f(elPupil_majorAxis));
-
-			BOOST_FOREACH(const cv::Point2f& centre, centres) {
-				Concurrency::parallel_for(0, params.StarburstPoints, [&] (int i) {
-					double theta = i * 2*PI/params.StarburstPoints;
-
-					// Initialise centre and direction vector
-					cv::Point2f pDir((float)std::cos(theta), (float)std::sin(theta));  
-
-					int t = 1;
-					cv::Point p = centre + (t * pDir);
-					while(p.inside(bbPupil))
-					{
-						uchar val = mPupilEdges(p);
-
-						if (val > 0)
-						{
-							float dx = mPupilSobelX(p);
-							float dy = mPupilSobelY(p);
-
-							float cdirx = p.x - (elPupilThresh.center.x - roiPupil.x);
-							float cdiry = p.y - (elPupilThresh.center.y - roiPupil.y);
-
-							// Check edge direction
-							double dirCheck = dx*cdirx + dy*cdiry;
-
-							if (dirCheck > 0)
-							{
-								// We've hit an edge
-								edgePointsConcurrent.push_back(cv::Point2f(p.x + 0.5f, p.y + 0.5f));
-								break;
-							}
-						}
-
-						++t;
-						p = centre + (t * pDir);
-					}
-				});
-			}
-
-			edgePoints = std::vector<cv::Point2f>(edgePointsConcurrent.begin(), edgePointsConcurrent.end());
-
-
-			// Remove duplicate edge points
-			std::sort(edgePoints.begin(), edgePoints.end(), [] (const cv::Point2f& p1, const cv::Point2f& p2) -> bool {
-				if (p1.x == p2.x)
-					return p1.y < p2.y;
-				else
-					return p1.x < p2.x;
-			});
-			edgePoints.erase( std::unique( edgePoints.begin(), edgePoints.end() ), edgePoints.end() );
-
-			if (edgePoints.size() < params.StarburstPoints/2)
-				return false;
-		}
-	}
-	else
-	{
-		SECTION("Non-zero value finder", log)
-		{
-			for(int y = 0; y < mPupilEdges.rows; y++)
-			{
-				uchar* val = mPupilEdges[y];
-				for(int x = 0; x < mPupilEdges.cols; x++, val++)
-				{
-					if(*val == 0)
-						continue;
-
-					edgePoints.push_back(cv::Point2f(x + 0.5f, y + 0.5f));
-				}
-			}
-		}
-	}
-
-
-	// ---------------------------
-	// Fit an ellipse to the edges
-	// ---------------------------
-
-	cv::RotatedRect elPupil;
-	std::vector<cv::Point2f> inliers;
-	SECTION("Ellipse fitting", log)
-	{
-		// Desired probability that only inliers are selected
-		const double p = 0.999;
-		// Probability that a point is an inlier
-		double w = params.PercentageInliers/100.0;
-		// Number of points needed for a model
-		const int n = 5;
-
-		if (params.PercentageInliers == 0)
-			return false;
-
-		if (edgePoints.size() >= n) // Minimum points for ellipse
-		{
-			// RANSAC!!!
-
-			double wToN = std::pow(w,n);
-			int k = static_cast<int>(std::log(1-p)/std::log(1 - wToN)  + 2*std::sqrt(1 - wToN)/wToN);
-
-			out.ransacIterations = k;
-
-			log.add("k", k);
-
-			//size_t threshold_inlierCount = std::max<size_t>(n, static_cast<size_t>(out.edgePoints.size() * 0.7));
-
-			// Use TBB for RANSAC
-			struct EllipseRansac_out {
-				std::vector<cv::Point2f> bestInliers;
-				cv::RotatedRect bestEllipse;
-				double bestEllipseGoodness;
-				int earlyRejections;
-				bool earlyTermination;
-
-				EllipseRansac_out() : bestEllipseGoodness(-std::numeric_limits<double>::infinity()), earlyTermination(false), earlyRejections(0) {}
-			};
-			struct EllipseRansac {
-				const TrackerParams& params;
-				const std::vector<cv::Point2f>& edgePoints;
-				int n;
-				const cv::Rect& bb;
-				const cv::Mat_<float>& mDX;
-				const cv::Mat_<float>& mDY;
-				int earlyRejections;
-				bool earlyTermination;
-
-				EllipseRansac_out out;
-
-				EllipseRansac(
-					const TrackerParams& params,
-					const std::vector<cv::Point2f>& edgePoints,
-					int n,
-					const cv::Rect& bb,
-					const cv::Mat_<float>& mDX,
-					const cv::Mat_<float>& mDY)
-					: params(params), edgePoints(edgePoints), n(n), bb(bb), mDX(mDX), mDY(mDY), earlyTermination(false), earlyRejections(0)
-				{
-				}
-
-				EllipseRansac(EllipseRansac& other, tbb::split)
-					: params(other.params), edgePoints(other.edgePoints), n(other.n), bb(other.bb), mDX(other.mDX), mDY(other.mDY), earlyTermination(other.earlyTermination), earlyRejections(other.earlyRejections)
-				{
-					//std::cout << "Ransac split" << std::endl;
-				}
-
-				void operator()(const tbb::blocked_range<size_t>& r)
-				{
-					if (out.earlyTermination)
-						return;
-					//std::cout << "Ransac start (" << (r.end()-r.begin()) << " elements)" << std::endl;
-					for( size_t i=r.begin(); i!=r.end(); ++i )
-					{
-						// Ransac Iteration
-						// ----------------
-						std::vector<cv::Point2f> sample;
-						if (params.Seed >= 0)
-							sample = randomSubset(edgePoints, n, static_cast<unsigned int>(i + params.Seed));
-						else
-							sample = randomSubset(edgePoints, n);
-
-						cv::RotatedRect ellipseSampleFit = fitEllipse(sample);
-						// Normalise ellipse to have width as the major axis.
-						if (ellipseSampleFit.size.height > ellipseSampleFit.size.width)
-						{
-							ellipseSampleFit.angle = std::fmod(ellipseSampleFit.angle + 90, 180);
-							std::swap(ellipseSampleFit.size.height, ellipseSampleFit.size.width);
-						}
-
-						cv::Size s = ellipseSampleFit.size;
-						// Discard useless ellipses early
-						if (!ellipseSampleFit.center.inside(bb)
-							|| s.height > params.Radius_Max*2
-							|| s.width > params.Radius_Max*2
-							|| s.height < params.Radius_Min*2 && s.width < params.Radius_Min*2
-							|| s.height > 4*s.width
-							|| s.width > 4*s.height
-							)
-						{
-							// Bad ellipse! Go to your room!
-							continue;
-						}
-
-						// Use conic section's algebraic distance as an error measure
-						ConicSection conicSampleFit(ellipseSampleFit);
-
-						// Check if sample's gradients are correctly oriented
-						if (params.EarlyRejection)
-						{
-							bool gradientCorrect = true;
-							BOOST_FOREACH(const cv::Point2f& p, sample)
-							{
-								cv::Point2f grad = conicSampleFit.algebraicGradientDir(p);
-								float dx = mDX(cv::Point(p));
-								float dy = mDY(cv::Point(p));
-
-								float dotProd = dx*grad.x + dy*grad.y;
-
-								gradientCorrect &= dotProd > 0;
-							}
-							if (!gradientCorrect)
-							{
-								out.earlyRejections++;
-								continue;
-							}
-						}
-
-						// Assume that the sample is the only inliers
-
-						cv::RotatedRect ellipseInlierFit = ellipseSampleFit;
-						ConicSection conicInlierFit = conicSampleFit;
-						std::vector<cv::Point2f> inliers, prevInliers;
-
-						// Iteratively find inliers, and re-fit the ellipse
-						for (int i = 0; i < params.InlierIterations; ++i)
-						{
-							// Get error scale for 1px out on the minor axis
-							cv::Point2f minorAxis(-std::sin(PI/180.0*ellipseInlierFit.angle), std::cos(PI/180.0*ellipseInlierFit.angle));
-							cv::Point2f minorAxisPlus1px = ellipseInlierFit.center + (ellipseInlierFit.size.height/2 + 1)*minorAxis;
-							float errOf1px = conicInlierFit.distance(minorAxisPlus1px);
-							float errorScale = 1.0f/errOf1px;
-
-							// Find inliers
-							inliers.reserve(edgePoints.size());
-							const float MAX_ERR = 2;
-							BOOST_FOREACH(const cv::Point2f& p, edgePoints)
-							{
-								float err = errorScale*conicInlierFit.distance(p);
-
-								if (err*err < MAX_ERR*MAX_ERR)
-									inliers.push_back(p);
-							}
-
-							if (inliers.size() < n) {
-								inliers.clear();
-								continue;
-							}
-
-							// Refit ellipse to inliers
-							ellipseInlierFit = fitEllipse(inliers);
-							conicInlierFit = ConicSection(ellipseInlierFit);
-
-							// Normalise ellipse to have width as the major axis.
-							if (ellipseInlierFit.size.height > ellipseInlierFit.size.width)
-							{
-								ellipseInlierFit.angle = std::fmod(ellipseInlierFit.angle + 90, 180);
-								std::swap(ellipseInlierFit.size.height, ellipseInlierFit.size.width);
-							}
-						}
-						if (inliers.empty())
-							continue;
-
-
-						// Calculate ellipse goodness
-						double ellipseGoodness = 0;
-						if (params.ImageAwareSupport)
-						{
-							BOOST_FOREACH(cv::Point2f& p, inliers)
-							{
-								cv::Point2f grad = conicInlierFit.algebraicGradientDir(p);
-								float dx = mDX(p);
-								float dy = mDY(p);
-
-								double edgeStrength = dx*grad.x + dy*grad.y;
-
-								ellipseGoodness += edgeStrength;
-							}
-						}
-						else
-						{
-							ellipseGoodness = inliers.size();
-						}
-
-						if (ellipseGoodness > out.bestEllipseGoodness)
-						{
-							std::swap(out.bestEllipseGoodness, ellipseGoodness);
-							std::swap(out.bestInliers, inliers);
-							std::swap(out.bestEllipse, ellipseInlierFit);
-
-							// Early termination, if 90% of points match
-							if (params.EarlyTerminationPercentage > 0 && out.bestInliers.size() > params.EarlyTerminationPercentage*edgePoints.size()/100)
-							{
-								earlyTermination = true;
-								break;
-							}
-						}
-
-					}
-					//std::cout << "Ransac end" << std::endl;
-				}
-
-				void join(EllipseRansac& other)
-				{
-					//std::cout << "Ransac join" << std::endl;
-					if (other.out.bestEllipseGoodness > out.bestEllipseGoodness)
-					{
-						std::swap(out.bestEllipseGoodness, other.out.bestEllipseGoodness);
-						std::swap(out.bestInliers, other.out.bestInliers);
-						std::swap(out.bestEllipse, other.out.bestEllipse);
-					}
-					out.earlyRejections += other.out.earlyRejections;
-					earlyTermination |= other.earlyTermination;
-
-					out.earlyTermination = earlyTermination;
-				}
-			};
-
-			EllipseRansac ransac(params, edgePoints, n, bbPupil, out.mPupilSobelX, out.mPupilSobelY);
-			try
-			{ 
-				tbb::parallel_reduce(tbb::blocked_range<size_t>(0,k,k/8), ransac);
-			}
-			catch (std::exception& e)
-			{
-				const char* c = e.what();
-				std::cerr << e.what() << std::endl;
-			}
-			inliers = ransac.out.bestInliers;
-			log.add("goodness", ransac.out.bestEllipseGoodness);
-
-			out.earlyRejections = ransac.out.earlyRejections;
-			out.earlyTermination = ransac.out.earlyTermination;
-
-
-			cv::RotatedRect ellipseBestFit = ransac.out.bestEllipse;
-			ConicSection conicBestFit(ellipseBestFit);
-			BOOST_FOREACH(const cv::Point2f& p, edgePoints)
-			{
-				cv::Point2f grad = conicBestFit.algebraicGradientDir(p);
-				float dx = out.mPupilSobelX(p);
-				float dy = out.mPupilSobelY(p);
-
-				out.edgePoints.push_back(EdgePoint(p, dx*grad.x + dy*grad.y));
-			}
-
-			elPupil = ellipseBestFit;
-			elPupil.center.x += roiPupil.x;
-			elPupil.center.y += roiPupil.y;
-		}
-
-		cv::Point2f pPupil = elPupil.center;
-
-		out.pPupil = pPupil;
-		out.elPupil = elPupil;
-		out.inliers = inliers;
-
-		return true;
-	}
-}

File PupilTracker.h

-#pragma once
-
-#include <vector>
-#include <string>
-
-#include <opencv2/core/core.hpp>
-
-#include "timer.h"
-#include "ConicSection.h"
-
-struct tracker_log
-{
-public:
-	typedef std::vector<std::pair<std::string, std::string>>::iterator iterator;
-	typedef std::vector<std::pair<std::string, std::string>>::const_iterator const_iterator;
-
-	template<typename T>
-	void add(const std::string& key, const T& val)
-	{
-		m_log.push_back(std::make_pair(key, boost::lexical_cast<std::string>(val)));
-	}
-	template<>
-	void add<timer>(const std::string& key, const timer& val)
-	{
-		std::stringstream ss;
-		ss.precision(2);
-		ss.setf(std::ios::fixed);
-		ss << (val.elapsed()*1000.0) << "ms";
-		m_log.push_back(std::make_pair(key, ss.str()));
-	}
-
-	iterator begin() { return m_log.begin(); }
-	const_iterator begin() const { return m_log.begin(); }
-	iterator end() { return m_log.end(); }
-	const_iterator end() const { return m_log.end(); }
-
-private:
-	std::vector<std::pair<std::string, std::string>> m_log;
-};
-
-namespace PupilTrackerNative {
-	
-	struct TrackerParams
-	{
-		int Radius_Min;
-		int Radius_Max;
-
-		double CannyBlur;
-		double CannyThreshold1;
-		double CannyThreshold2;
-		int StarburstPoints;
-
-		int PercentageInliers;
-		int InlierIterations;
-		bool ImageAwareSupport;
-		int EarlyTerminationPercentage;
-		bool EarlyRejection;
-		int Seed;
-	};
-
-	const cv::Point2f UNKNOWN_POSITION = cv::Point2f(-1,-1);
-
-	struct EdgePoint
-	{
-		cv::Point2f point;
-		double edgeStrength;
-
-		EdgePoint(const cv::Point2f& p, double s) : point(p), edgeStrength(s) {}
-		EdgePoint(float x, float y, double s) : point(x,y), edgeStrength(s) {}
-
-		bool operator== (const EdgePoint& other)
-		{
-			return point == other.point;
-		}
-	};
-
-	struct findPupilEllipse_out {
-		cv::Rect roiHaarPupil;
-		cv::Mat_<uchar> mHaarPupil;
-		
-		cv::Mat_<float> histPupil;
-		double threshold;
-		cv::Mat_<uchar> mPupilThresh;
-
-		cv::Rect bbPupilThresh;
-		cv::RotatedRect elPupilThresh;
-
-		cv::Rect roiPupil;
-		cv::Mat_<uchar> mPupil;
-		cv::Mat_<uchar> mPupilOpened;
-		cv::Mat_<uchar> mPupilBlurred;
-		cv::Mat_<uchar> mPupilEdges;
-		cv::Mat_<float> mPupilSobelX;
-		cv::Mat_<float> mPupilSobelY;
-
-		std::vector<EdgePoint> edgePoints;
-		std::vector<cv::Point2f> inliers;
-		int ransacIterations;
-		int earlyRejections;
-		bool earlyTermination;
-
-		cv::Point2f pPupil;
-		cv::RotatedRect elPupil;
-
-		findPupilEllipse_out() : pPupil(UNKNOWN_POSITION), threshold(-1) {}
-	};
-	bool findPupilEllipse(
-		const TrackerParams& params,
-		const cv::Mat& m,
-
-		findPupilEllipse_out& out,
-		tracker_log& log
-		);
-
-}//PupilTrackerNative

File cmd/CMakeLists.txt

+find_package (OpenCV REQUIRED)
+
+add_executable (
+    pupiltrackercmd
+
+    PupilTrackerCmd.cpp
+)
+
+target_link_libraries(
+    pupiltrackercmd
+
+    pupiltracker
+    ${OpenCV_LIBS}
+    tbb
+)

File cmd/PupilTrackerCmd.cpp

+#include <iostream>
+
+#include <opencv2/highgui/highgui.hpp>
+
+#include "../lib/PupilTracker.h"
+#include "../lib/cvx.h"
+
+void imshowscale(const std::string& name, cv::Mat& m, double scale)
+{
+    cv::Mat res;
+    cv::resize(m, res, cv::Size(), scale, scale, cv::INTER_NEAREST);
+    cv::imshow(name, res);
+}
+
+int main(int argc, char* argv[])
+{
+    if (argc < 2) {
+        std::cerr << "Need filename" << std::endl;
+        return 1;
+    }
+
+    std::cout << "Opening " << argv[1] << std::endl;
+    cv::VideoCapture vc(argv[1]);
+    if (!vc.isOpened()) {
+        std::cerr << "Could not open " << argv[1] << std::endl;
+        return 2;
+    }
+
+
+    cv::Mat m;
+    while (true)
+    {
+        vc >> m;
+        if (m.empty())
+        {
+            vc.open(argv[1]);
+            if (!vc.isOpened()) {
+                std::cerr << "Could not open " << argv[1] << std::endl;
+                return 2;
+            }
+            vc >> m;
+            if (m.empty()) {
+                return 1;
+            }
+        }
+
+
+        PupilTracker::TrackerParams params;
+		params.Radius_Min = 10;
+		params.Radius_Max = 60;
+
+		params.CannyBlur = 1.6;
+		params.CannyThreshold1 = 30;
+		params.CannyThreshold2 = 50;
+		params.StarburstPoints = 0;
+
+		params.PercentageInliers = 40;
+		params.InlierIterations = 2;
+		params.ImageAwareSupport = true;
+		params.EarlyTerminationPercentage = 95;
+		params.EarlyRejection = true;
+		params.Seed = -1;
+
+        PupilTracker::findPupilEllipse_out out;
+        tracker_log log;
+        PupilTracker::findPupilEllipse(params, m, out, log); 
+
+        cvx::cross(m, out.pPupil, 5, CV_RGB(255,255,0));
+        cv::ellipse(m, out.elPupil, CV_RGB(255,0,255));
+
+
+        imshowscale("Haar Pupil", out.mHaarPupil, 3);
+        imshowscale("Pupil", out.mPupil, 3);
+        //imshowscale("Thresh Pupil", out.mPupilThresh, 3);
+        imshowscale("Edges", out.mPupilEdges, 3);
+        cv::imshow("Result", m);
+
+        if (cv::waitKey(10) != -1)
+            break;
+    }
+}

File cvx.cpp

-#include "stdafx.h"
-
-#include "cvx.h"
-
-
-void cvx::blurAndDecimate(const cv::Mat& src, cv::Mat& dst, double sigma, int& octave, int borderType)
-{
-	cv::Mat tmp;
-
-	double sigma_rem;
-	octave = sigma2octave(sigma, sigma_rem);
-
-	// Blur and decimate image down to the correct octave
-	tmp = src;
-	for (int i = 0; i < octave; ++i)
-	{
-		cv::pyrDown(tmp, dst, cv::Size(), borderType);
-		tmp = dst;
-	}
-
-	// Scale sigma_rem to octave scale, and do the remaining blur
-	sigma_rem /= pow2(octave);
-	cv::GaussianBlur(tmp, dst, cv::Size(), sigma_rem, sigma_rem, borderType);
-}
-
-
-void cvx::getROI(const cv::Mat& src, cv::Mat& dst, const cv::Rect& roi, int borderType)
-{
-	cv::Rect bbSrc = boundingBox(src);
-	cv::Rect validROI = roi & bbSrc;
-	if (validROI == roi)
-	{
-		dst = cv::Mat(src, validROI);
-	}
-	else
-	{
-		// Figure out how much to add on for top, left, right and bottom
-		cv::Point tl = roi.tl() - bbSrc.tl();
-		cv::Point br = roi.br() - bbSrc.br();
-
-		int top = std::max(-tl.y, 0);  // Top and left are negated because adding a border
-		int left = std::max(-tl.x, 0); // goes "the wrong way"
-		int right = std::max(br.x, 0);
-		int bottom = std::max(br.y, 0);
-
-		cv::Mat tmp(src, validROI);
-		cv::copyMakeBorder(tmp, dst, top, bottom, left, right, borderType);
-	}
-}
-
-
-void cvx::hessianRatio(const cv::Mat& src, cv::Mat& dst)
-{
-	int type = src.type();
-
-	cv::Mat src_x, src_y;
-	cv::Sobel(src, src_x, type, 1, 0);
-	cv::Sobel(src, src_y, type, 0, 1);
-	cv::Mat src_xx, src_xy, src_yy;
-	cv::Sobel(src_x, src_xx, type, 1, 0);
-	cv::Sobel(src_x, src_xy, type, 0, 1);
-	cv::Sobel(src_y, src_yy, type, 0, 1);
-
-	cv::Mat src_h_trace = src_xx + src_yy;
-	cv::Mat src_h_trace2;
-	cv::multiply(src_h_trace, src_h_trace, src_h_trace2);
-
-	cv::Mat src_xx_yy;
-	cv::multiply(src_xx, src_yy, src_xx_yy);
-	cv::Mat src_xy2;
-	cv::multiply(src_xy, src_xy, src_xy);
-
-	cv::Mat src_h_det = src_xx_yy - src_xy2;
-
-	cv::divide(src_h_trace2, src_h_det, dst);
-}
-
-float cvx::histKmeans(const cv::Mat_<float>& hist, int bin_min, int bin_max, int K, float init_centres[], cv::Mat_<uchar>& labels, cv::TermCriteria termCriteria)
-{
-	CV_Assert( hist.rows == 1 || hist.cols == 1 && K > 0 );
-
-	labels = cv::Mat_<uchar>::zeros(hist.size());
-	int nbins = hist.total();
-	float binWidth = (bin_max - bin_min)/nbins;
-	float binStart = bin_min + binWidth/2;
-
-	cv::Mat_<float> centres(K, 1, init_centres, 4);
-
-	int iters = 0;
-	bool finalRun = false;
-	while (true)
-	{
-		++iters;
-		cv::Mat_<float> old_centres = centres.clone();
-
-		int i_bin;
-		cv::Mat_<float>::const_iterator i_hist;
-		cv::Mat_<uchar>::iterator i_labels;
-		cv::Mat_<float>::iterator i_centres;
-		uchar label;
-
-		float sumDist = 0;
-		int movedCount = 0;
-
-		// Step 1. Assign each element a label
-		for (i_bin = 0, i_labels = labels.begin(), i_hist = hist.begin();
-			 i_bin < nbins;
-			 ++i_bin, ++i_labels, ++i_hist)
-		{
-			float bin_val = binStart + i_bin*binWidth;
-			float minDist = sq(bin_val - centres(*i_labels));
-			int curLabel = *i_labels;
-
-			for (label = 0; label < K; ++label)
-			{
-				float dist = sq(bin_val - centres(label));
-				if (dist < minDist)
-				{
-					minDist = dist;
-					*i_labels = label;
-				}
-			}
-
-			if (*i_labels != curLabel)
-				movedCount++;
-
-			sumDist += (*i_hist) * std::sqrt(minDist);
-		}
-
-		if (finalRun)
-			return sumDist;
-
-		// Step 2. Recalculate centres
-		cv::Mat_<float> counts(K, 1, 0.0f);
-		for (i_bin = 0, i_labels = labels.begin(), i_hist = hist.begin();
-			 i_bin < nbins;
-			 ++i_bin, ++i_labels, ++i_hist)
-		{
-			float bin_val = binStart + i_bin*binWidth;
-
-			centres(*i_labels) += (*i_hist) * bin_val;
-			counts(*i_labels) += *i_hist;
-		}
-		for (label = 0; label < K; ++label)
-		{
-			if (counts(label) == 0)
-				return std::numeric_limits<float>::infinity();
-
-			centres(label) /= counts(label);
-		}
-
-		// Step 3. Detect termination criteria
-		if (movedCount == 0)
-			finalRun = true;
-		else if (termCriteria.type | cv::TermCriteria::COUNT && iters >= termCriteria.maxCount)
-			finalRun = true;
-		else if (termCriteria.type | cv::TermCriteria::EPS)
-		{
-			float max_movement = 0;
-			for (label = 0; label < K; ++label)
-			{
-				max_movement = std::max(max_movement, sq(centres(label) - old_centres(label)));
-			}
-			if (sqrt(max_movement) < termCriteria.epsilon)
-				finalRun = true;
-		}
-	}
-	return std::numeric_limits<float>::infinity();
-}
-
-cv::RotatedRect cvx::fitEllipse(const cv::Moments& m)
-{
-	cv::RotatedRect ret;
-
-	ret.center.x = m.m10/m.m00;
-	ret.center.y = m.m01/m.m00;
-
-	double mu20 = m.m20/m.m00 - ret.center.x*ret.center.x;
-	double mu02 = m.m02/m.m00 - ret.center.y*ret.center.y;
-	double mu11 = m.m11/m.m00 - ret.center.x*ret.center.y;
-
-	double common = std::sqrt(sq(mu20 - mu02) + 4*sq(mu11));
-
-	ret.size.width = std::sqrt(2*(mu20 + mu02 + common));
-	ret.size.height = std::sqrt(2*(mu20 + mu02 - common));
-
-	double num, den;
-	if (mu02 > mu20) {
-		num = mu02 - mu20 + common;
-		den = 2*mu11;
-	}
-	else {
-		num = 2*mu11;
-		den = mu20 - mu02 + common;
-	}
-
-	if (num == 0 && den == 0)
-		ret.angle = 0;
-	else
-		ret.angle = (180/PI) * std::atan2(num,den);
-
-	return ret;
-}
-cv::Vec2f cvx::majorAxis(const cv::RotatedRect& ellipse)
-{
-	return cv::Vec2f(ellipse.size.width*std::cos(PI/180*ellipse.angle), ellipse.size.width*std::sin(PI/180*ellipse.angle));
-}

File cvx.h

-#ifndef cvx_h__
-#define cvx_h__
-
-#include <opencv2/core/core.hpp>
-#include <opencv2/imgproc/imgproc.hpp>
-
-#include "utils.h"
-
-namespace cvx
-{
-	// Octave<->sigma conversion functions
-	// Used Matlab to minimise difference between Gaussian of sigma 's'
-	// versus pyramid image of level 'n'.
-	// Got
-	//   n |  0   1   2   3   4   5   6   7   8 
-	//   s |  0
-	// Modelled as:
-	//   s = sqrt( (a^n - 1) / b )
-	// where
-	//   a = 3.9512
-	//   b = 1.4182
-	// and
-	//   n = log_a(b * s^2 + 1)
-
-	namespace
-	{
-		const double a = 3.9512;
-		const double b = 1.4182;
-		const double loga = std::log(a);
-	}
-
-	inline double octave2sigmaSq(int octave, double sigma_rem = 0.0)
-	{
-		return (std::pow(a, octave) - 1)/b + sq(sigma_rem);
-	}
-	inline double octave2sigma(int octave, double sigma_rem = 0.0)
-	{
-		return std::sqrt(octave2sigmaSq(octave, sigma_rem));
-	}
-	inline int sigma2octave(double sigma)
-	{
-		return (int)std::floor( std::log( b * sq(sigma) + 1) / loga );
-	}
-	inline int sigma2octave(double sigma, double& sigma_rem)
-	{
-		int octave = sigma2octave(sigma);
-		sigma_rem = std::sqrt( sq(sigma) - octave2sigmaSq(octave));
-		return octave;
-	}
-
-
-	template<typename T>
-	inline cv::Rect_<T> roiAround(const cv::Point_<T>& centre, T radius)
-	{
-		return roiAround(centre.x, centre.y, radius);
-	}
-	template<typename T>
-	inline cv::Rect_<T> roiAround(T x, T y, T radius)
-	{
-		return cv::Rect_<T>(x - radius, y - radius, 2*radius + 1, 2*radius + 1);
-	}
-
-
-
-	inline void cross(cv::Mat& img, cv::Point centre, int radius, const cv::Scalar& colour, int thickness = 1, int lineType = 8, int shift = 0)
-	{
-		cv::line(img, centre + cv::Point(-radius, -radius), centre + cv::Point(radius, radius), colour, thickness, lineType, shift = 0);
-		cv::line(img, centre + cv::Point(-radius, radius), centre + cv::Point(radius, -radius), colour, thickness, lineType, shift = 0);
-	}
-	inline void plus(cv::Mat& img, cv::Point centre, int radius, const cv::Scalar& colour, int thickness = 1, int lineType = 8, int shift = 0)
-	{
-		cv::line(img, centre + cv::Point(0, -radius), centre + cv::Point(0, radius), colour, thickness, lineType, shift = 0);
-		cv::line(img, centre + cv::Point(-radius, 0), centre + cv::Point(radius, 0), colour, thickness, lineType, shift = 0);
-	}
-
-	inline cv::Rect boundingBox(const cv::Mat& img)
-	{
-		return cv::Rect(0,0,img.cols,img.rows);
-	}
-
-	void blurAndDecimate(const cv::Mat& src, cv::Mat& dst, double sigma, int& octave, int borderType = cv::BORDER_REPLICATE);
-	void getROI(const cv::Mat& src, cv::Mat& dst, const cv::Rect& roi, int borderType = cv::BORDER_REPLICATE);
-	void hessianRatio(const cv::Mat& src, cv::Mat& dst);
-
-	float histKmeans(const cv::Mat_<float>& hist, int bin_min, int bin_max, int K, float init_centres[], cv::Mat_<uchar>& labels, cv::TermCriteria termCriteria);
-
-	cv::RotatedRect fitEllipse(const cv::Moments& m);
-	cv::Vec2f majorAxis(const cv::RotatedRect& ellipse);
-
-
-
-
-}
-
-#endif // cvx_h__

File high_resolution_timer.hpp

-// Copyright (c) 2005 Hartmut Kaiser
-// Copyright (c) 2005 Christopher Diggins
-//
-// Distributed under the Boost Software License, Version 1.0.
-// (See accompanying file LICENSE_1_0.txt or copy at 
-// http://www.boost.org/LICENSE_1_0.txt)
-//
-// Disclaimer: Not a Boost library
-
-#if !defined(BOOST_HIGH_RESOLUTION_TIMER_HPP)
-#define BOOST_HIGH_RESOLUTION_TIMER_HPP
-
-#include <boost/config.hpp>
-#include <boost/timer.hpp>
-
-#if !defined(BOOST_WINDOWS)
-
-//  For platforms other than Windows, simply fall back to boost::timer
-namespace boost {
-	typedef boost::timer high_resolution_timer;
-}
-
-#else
-
-#include <stdexcept>
-#include <limits>
-#include <windows.h>
-
-namespace boost {
-
-	///////////////////////////////////////////////////////////////////////////////
-	//
-	//  high_resolution_timer 
-	//      A timer object measures elapsed time.
-	//      CAUTION: Windows only!
-	//
-	///////////////////////////////////////////////////////////////////////////////
-	class high_resolution_timer
-	{
-	public:
-		// ctor
-		high_resolution_timer() 
-		{
-			start_time.QuadPart = 0;
-			frequency.QuadPart = 0;
-
-			if (!QueryPerformanceFrequency(&frequency))
-				throw std::runtime_error("Couldn't acquire frequency");
-
-			restart(); 
-		} 
-
-		// restart timer
-		void restart() 
-		{ 
-			t.restart();
-			if (!QueryPerformanceCounter(&start_time))
-				throw std::runtime_error("Couldn't initialize start_time");
-		} 
-
-		// return elapsed time in seconds
-		double elapsed() const                  
-		{ 
-			LARGE_INTEGER now;
-			if (!QueryPerformanceCounter(&now))
-				throw std::runtime_error("Couldn't get current time");
-
-			// QueryPerformanceCounter() workaround
-			// http://support.microsoft.com/default.aspx?scid=kb;EN-US;q274323
-			double d1 = double(now.QuadPart - start_time.QuadPart) / frequency.QuadPart;
-			double d2 = t.elapsed();
-			return ((d1 - d2) > 0.5) ? d2 : d1;
-		}
-
-		// return estimated maximum value for elapsed()
-		double elapsed_max() const   
-		{
-			return (double((std::numeric_limits<LONGLONG>::max)())
-				- double(start_time.QuadPart)) / double(frequency.QuadPart); 
-		}
-
-		// return minimum value for elapsed()
-		double elapsed_min() const            
-		{ 
-			return 1.0 / frequency.QuadPart; 
-		}
-
-	private:
-		timer t; // backup in case of QueryPerformanceCounter() bug
-		LARGE_INTEGER start_time;
-		LARGE_INTEGER frequency;
-	}; 
-
-} // namespace boost
-
-#endif  // !defined(BOOST_WINDOWS)
-
-#endif  // !defined(BOOST_HIGH_RESOLUTION_TIMER_HPP)

File lib/CMakeLists.txt

+add_library (
+    pupiltracker
+
+    PupilTracker.cpp
+    cvx.cpp
+    utils.cpp
+)

File lib/ConicSection.h

+#ifndef __CONIC_SECTION_H__
+#define __CONIC_SECTION_H__
+
+template<typename T>
+class ConicSection_
+{
+public:
+	T A,B,C,D,E,F;
+
+	ConicSection_(cv::RotatedRect r)
+	{
+		cv::Point_<T> axis((T)std::cos(CV_PI/180.0 * r.angle), (T)std::sin(CV_PI/180.0 * r.angle));
+		cv::Point_<T> centre(r.center);
+		T a = r.size.width/2;
+		T b = r.size.height/2;
+
+		initFromEllipse(axis, centre, a, b);
+	}
+
+	T algebraicDistance(cv::Point_<T> p)
+	{
+		return A*p.x*p.x + B*p.x*p.y + C*p.y*p.y + D*p.x + E*p.y + F;
+	}
+
+	T distance(cv::Point_<T> p)
+	{
+		//    dist
+		// -----------
+		// |grad|^0.45
+
+		T dist = algebraicDistance(p);
+		cv::Point_<T> grad = algebraicGradient(p);
+
+		T sqgrad = grad.dot(grad);
+
+		return dist / std::pow(sqgrad, T(0.45/2));
+	}
+
+	cv::Point_<T> algebraicGradient(cv::Point_<T> p)
+	{
+		return cv::Point_<T>(2*A*p.x + B*p.y + D, B*p.x + 2*C*p.y + E);
+	}
+
+	cv::Point_<T> algebraicGradientDir(cv::Point_<T> p)
+	{
+		cv::Point_<T> grad = algebraicGradient(p);
+		T len = std::sqrt(grad.ddot(grad));
+		grad.x /= len;
+		grad.y /= len;
+		return grad;
+	}
+
+protected:
+	void initFromEllipse(cv::Point_<T> axis, cv::Point_<T> centre, T a, T b)
+	{
+		T a2 = a * a;
+		T b2 = b * b;
+
+		A = axis.x*axis.x / a2 + axis.y*axis.y / b2;
+		B = 2*axis.x*axis.y / a2 - 2*axis.x*axis.y / b2;
+		C = axis.y*axis.y / a2 + axis.x*axis.x / b2;
+		D = (-2*axis.x*axis.y*centre.y - 2*axis.x*axis.x*centre.x) / a2
+			+ (2*axis.x*axis.y*centre.y - 2*axis.y*axis.y*centre.x) / b2;
+		E = (-2*axis.x*axis.y*centre.x - 2*axis.y*axis.y*centre.y) / a2
+			+ (2*axis.x*axis.y*centre.x - 2*axis.x*axis.x*centre.y) / b2;
+		F = (2*axis.x*axis.y*centre.x*centre.y + axis.x*axis.x*centre.x*centre.x + axis.y*axis.y*centre.y*centre.y) / a2
+			+ (-2*axis.x*axis.y*centre.x*centre.y + axis.y*axis.y*centre.x*centre.x + axis.x*axis.x*centre.y*centre.y) / b2
+			- 1;
+	}
+};
+typedef ConicSection_<float> ConicSection;
+
+#endif // __CONIC_SECTION_H__

File lib/PupilTracker.cpp

+#include "PupilTracker.h"
+
+#include <iostream>
+
+#include <boost/foreach.hpp>
+
+#include <opencv2/core/core.hpp>
+#include <opencv2/imgproc/imgproc.hpp>
+#include <opencv2/highgui/highgui.hpp>
+
+#include <tbb/tbb.h>
+
+#include "cvx.h"
+
+namespace 
+{
+	struct section_guard
+	{
+		std::string name;
+		tracker_log& log;
+		timer t;
+		section_guard(const std::string& name, tracker_log& log) : name(name), log(log), t() {  }
+		~section_guard() { log.add(name, t); }
+		operator bool() const {return false;}
+	};
+
+	inline section_guard make_section_guard(const std::string& name, tracker_log& log)
+	{
+		return section_guard(name,log);
+	}
+}
+
+#define SECTION(A,B) if (const section_guard& _section_guard_ = make_section_guard( A , B )) {} else
+
+
+
+
+class HaarSurroundFeature
+{
+public:
+	HaarSurroundFeature(int r1, int r2) : r_inner(r1), r_outer(r2)
+	{
+		//  _________________
+		// |        -ve      |
+		// |     _______     |
+		// |    |   +ve |    |
+		// |    |   .   |    |
+		// |    |_______|    |
+		// |         <r1>    |
+		// |_________<--r2-->|
+
+		// Number of pixels in each part of the kernel
+		int count_inner = r_inner*r_inner;
+		int count_outer = r_outer*r_outer - r_inner*r_inner;
+
+		// Frobenius normalized values
+		//
+		// Want norm = 1 where norm = sqrt(sum(pixelvals^2)), so:
+		//  sqrt(count_inner*val_inner^2 + count_outer*val_outer^2) = 1
+		//
+		// Also want sum(pixelvals) = 0, so:
+		//  count_inner*val_inner + count_outer*val_outer = 0
+		//
+		// Solving both of these gives:
+		//val_inner = std::sqrt( (double)count_outer/(count_inner*count_outer + sq(count_inner)) );
+		//val_outer = -std::sqrt( (double)count_inner/(count_inner*count_outer + sq(count_outer)) );
+
+		// Square radius normalised values
+		//
+		// Want the response to be scale-invariant, so scale it by the number of pixels inside it:
+		//  val_inner = 1/count = 1/r_outer^2
+		//
+		// Also want sum(pixelvals) = 0, so:
+		//  count_inner*val_inner + count_outer*val_outer = 0
+		//
+		// Hence:
+		val_inner = 1.0 / (r_inner*r_inner);
+		val_outer = -val_inner*count_inner/count_outer;
+
+	}
+
+	double val_inner, val_outer;
+	int r_inner, r_outer;
+};
+
+cv::RotatedRect fitEllipse(const std::vector<PupilTracker::EdgePoint>& edgePoints)
+{
+	std::vector<cv::Point2f> points;
+	points.reserve(edgePoints.size());
+
+	BOOST_FOREACH(const PupilTracker::EdgePoint& e, edgePoints)
+		points.push_back(e.point);
+
+	return cv::fitEllipse(points);
+}
+
+
+bool PupilTracker::findPupilEllipse(
+	const TrackerParams& params,
+	const cv::Mat& m,
+
+	PupilTracker::findPupilEllipse_out& out,
+	tracker_log& log
+	)
+{
+	// --------------------
+	// Convert to greyscale
+	// --------------------
+
+	cv::Mat_<uchar> mEye;
+
+	SECTION("Grey and crop", log)
+	{
+		// Pick one channel if necessary, and crop it to get rid of borders
+		if (m.channels() == 1)
+		{
+			mEye = m;
+		}
+		else if (m.channels() == 3)
+		{
+			cv::cvtColor(m, mEye, CV_BGR2GRAY);
+		}
+		else if (m.channels() == 4)
+		{
+			cv::cvtColor(m, mEye, CV_BGRA2GRAY);
+		}
+		else
+		{
+			throw std::runtime_error("Unsupported number of channels");
+		}
+	}
+
+	// -----------------------
+	// Find best haar response
+	// -----------------------
+
+	//             _____________________
+	//            |         Haar kernel |
+	//            |                     |
+	//  __________|______________       |
+	// | Image    |      |       |      |
+	// |    ______|______|___.-r-|--2r--|
+	// |   |      |      |___|___|      |
+	// |   |      |          |   |      |
+	// |   |      |          |   |      |
+	// |   |      |__________|___|______|
+	// |   |    Search       |   |
+	// |   |    region       |   |
+	// |   |                 |   |
+	// |   |_________________|   |
+	// |                         |
+	// |_________________________|
+	//
+
+	cv::Mat_<int32_t> mEyeIntegral;
+	int padding = 2*params.Radius_Max;
+
+	SECTION("Integral image", log)
+	{
+		cv::Mat mEyePad;
+		// Need to pad by an additional 1 to get bottom & right edges.
+		cv::copyMakeBorder(mEye, mEyePad, padding, padding, padding, padding, cv::BORDER_REPLICATE);
+		cv::integral(mEyePad, mEyeIntegral);
+	}
+
+	cv::Point2f pHaarPupil;
+	int haarRadius;
+
+	SECTION("Haar responses", log)
+	{
+		const int rstep = 2;
+		const int ystep = 4;
+		const int xstep = 4;
+
+		double minResponse = std::numeric_limits<double>::infinity();
+
+		for (int r = params.Radius_Min; r < params.Radius_Max; r+=rstep)
+		{   
+			// Get Haar feature
+			int r_inner = r;
+			int r_outer = 3*r;
+			HaarSurroundFeature f(r_inner, r_outer);
+
+			// Use TBB for rows
+			std::pair<double,cv::Point2f> minRadiusResponse = tbb::parallel_reduce(
+				tbb::blocked_range<int>(0, (mEye.rows-r - r - 1)/ystep + 1, ((mEye.rows-r - r - 1)/ystep + 1) / 8),
+				std::make_pair(std::numeric_limits<double>::infinity(), PupilTracker::UNKNOWN_POSITION),
+				[&] (tbb::blocked_range<int> range, const std::pair<double,cv::Point2f>& minValIn) -> std::pair<double,cv::Point2f>
+			{
+				std::pair<double,cv::Point2f> minValOut = minValIn;
+				for (int i = range.begin(), y = r + range.begin()*ystep; i < range.end(); i++, y += ystep)
+				{
+					//            �         �
+					// row1_outer.|         |  p00._____________________.p01
+					//            |         |     |         Haar kernel |
+					//            |         |     |                     |
+					// row1_inner.|         |     |   p00._______.p01   |
+					//            |-padding-|     |      |       |      |
+					//            |         |     |      | (x,y) |      |
+					// row2_inner.|         |     |      |_______|      |
+					//            |         |     |   p10'       'p11   |
+					//            |         |     |                     |
+					// row2_outer.|         |     |_____________________|
+					//            |         |  p10'                     'p11
+					//            �         �
+
+					int* row1_inner = mEyeIntegral[y+padding - r_inner];
+					int* row2_inner = mEyeIntegral[y+padding + r_inner + 1];
+					int* row1_outer = mEyeIntegral[y+padding - r_outer];
+					int* row2_outer = mEyeIntegral[y+padding + r_outer + 1];
+
+					int* p00_inner = row1_inner + r + padding - r_inner;
+					int* p01_inner = row1_inner + r + padding + r_inner + 1;
+					int* p10_inner = row2_inner + r + padding - r_inner;
+					int* p11_inner = row2_inner + r + padding + r_inner + 1;
+
+					int* p00_outer = row1_outer + r + padding - r_outer;
+					int* p01_outer = row1_outer + r + padding + r_outer + 1;
+					int* p10_outer = row2_outer + r + padding - r_outer;
+					int* p11_outer = row2_outer + r + padding + r_outer + 1;
+
+					for (int x = r; x < mEye.cols - r; x+=xstep)
+					{
+						int sumInner = *p00_inner + *p11_inner - *p01_inner - *p10_inner;
+						int sumOuter = *p00_outer + *p11_outer - *p01_outer - *p10_outer - sumInner;
+
+						double response = f.val_inner * sumInner + f.val_outer * sumOuter;
+
+						if (response < minValOut.first)
+						{
+							minValOut.first = response;
+							minValOut.second = cv::Point(x,y);
+						}
+
+						p00_inner += xstep;
+						p01_inner += xstep;
+						p10_inner += xstep;
+						p11_inner += xstep;
+
+						p00_outer += xstep;
+						p01_outer += xstep;
+						p10_outer += xstep;
+						p11_outer += xstep;
+					}
+				}
+				return minValOut;
+			},
+				[] (const std::pair<double,cv::Point2f>& x, const std::pair<double,cv::Point2f>& y) -> std::pair<double,cv::Point2f>
+			{
+				if (x.first < y.first)
+					return x;
+				else
+					return y;
+			}
+			);
+
+			if (minRadiusResponse.first < minResponse)
+			{
+				minResponse = minRadiusResponse.first;
+				// Set return values
+				pHaarPupil = minRadiusResponse.second;
+				haarRadius = r;
+			}
+		}
+	}
+    // Paradoxically, a good Haar fit won't catch the entire pupil, so expand it a bit
+    haarRadius = (int)(haarRadius * SQRT_2);
+
+	// ---------------------------
+	// Pupil ROI around Haar point
+	// ---------------------------
+	cv::Rect roiHaarPupil = cvx::roiAround(cv::Point(pHaarPupil.x, pHaarPupil.y), haarRadius);
+	cv::Mat_<uchar> mHaarPupil;
+	cvx::getROI(mEye, mHaarPupil, roiHaarPupil);
+
+	out.roiHaarPupil = roiHaarPupil;
+	out.mHaarPupil = mHaarPupil;
+
+	// --------------------------------------------------
+	// Get histogram of pupil region, segment with KMeans
+	// --------------------------------------------------
+
+	const int bins = 256;
+
+	cv::Mat_<float> hist;
+	SECTION("Histogram", log)
+	{
+		int channels[] = {0};
+		int sizes[] = {bins};
+		float range[2] = {0, 256};
+		const float* ranges[] = {range};
+		cv::calcHist(&mHaarPupil, 1, channels, cv::Mat(), hist, 1, sizes, ranges);
+	}
+
+	out.histPupil = hist;
+
+	float threshold;
+	SECTION("KMeans", log)
+	{
+		// Try various candidate centres, return the one with minimal label distance
+		float candidate0[2] = {0, 0};
+		float candidate1[2] = {128, 255};
+		float bestDist = std::numeric_limits<float>::infinity();
+		float bestThreshold = std::numeric_limits<float>::quiet_NaN();
+
+		for (int i = 0; i < 2; i++)
+		{
+			cv::Mat_<uchar> labels;
+			float centres[2] = {candidate0[i], candidate1[i]};
+			float dist = cvx::histKmeans(hist, 0, 256, 2, centres, labels, cv::TermCriteria(cv::TermCriteria::COUNT, 50, 0.0));
+
+			float thisthreshold = (centres[0] + centres[1])/2;
+			if (dist < bestDist && boost::math::isnormal(thisthreshold))
+			{
+				bestDist = dist;
+				bestThreshold = thisthreshold;
+			}
+		}
+
+		if (!boost::math::isnormal(bestThreshold))
+		{
+			// If kmeans gives a degenerate solution, exit early
+			return false;
+		}
+
+		threshold = bestThreshold;
+	}
+
+	cv::Mat_<uchar> mPupilThresh;
+	SECTION("Threshold", log)
+	{
+		cv::threshold(mHaarPupil, mPupilThresh, threshold, 255, cv::THRESH_BINARY_INV);
+	}
+
+	out.threshold = threshold;
+	out.mPupilThresh = mPupilThresh;
+
+	// ---------------------------------------------
+	// Find best region in the segmented pupil image
+	// ---------------------------------------------
+
+	cv::Rect bbPupilThresh;
+	cv::RotatedRect elPupilThresh;
+
+	SECTION("Find best region", log)
+	{
+		cv::Mat_<uchar> mPupilContours = mPupilThresh.clone();
+		std::vector<std::vector<cv::Point> > contours;
+		cv::findContours(mPupilContours, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE);
+
+		if (contours.size() == 0)
+			return false;
+
+		std::vector<cv::Point>& maxContour = contours[0];
+		double maxContourArea = cv::contourArea(maxContour);
+		BOOST_FOREACH(std::vector<cv::Point>& c, contours)
+		{
+			double area = cv::contourArea(c);
+			if (area > maxContourArea)
+			{
+				maxContourArea = area;
+				maxContour = c;
+			}
+		}
+
+		cv::Moments momentsPupilThresh = cv::moments(maxContour);
+
+		bbPupilThresh = cv::boundingRect(maxContour);
+		elPupilThresh = cvx::fitEllipse(momentsPupilThresh);
+
+		// Shift best region into eye coords (instead of pupil region coords), and get ROI
+		bbPupilThresh.x += roiHaarPupil.x;
+		bbPupilThresh.y += roiHaarPupil.y;
+		elPupilThresh.center.x += roiHaarPupil.x;
+		elPupilThresh.center.y += roiHaarPupil.y;
+	}
+
+	out.bbPupilThresh = bbPupilThresh;
+	out.elPupilThresh = elPupilThresh;
+
+	// ------------------------------
+	// Find edges in new pupil region
+	// ------------------------------
+
+	cv::Mat_<uchar> mPupil, mPupilOpened, mPupilBlurred, mPupilEdges;
+	cv::Mat_<float> mPupilSobelX, mPupilSobelY;
+	cv::Rect bbPupil;
+	cv::Rect roiPupil = cvx::roiAround(cv::Point(elPupilThresh.center.x, elPupilThresh.center.y), haarRadius);
+	SECTION("Pupil preprocessing", log)
+	{
+        const int padding = 3;
+
+		cv::Rect roiPadded(roiPupil.x-padding, roiPupil.y-padding, roiPupil.width+2*padding, roiPupil.height+2*padding);
+		// First get an ROI around the approximate pupil location
+		cvx::getROI(mEye, mPupil, roiPadded, cv::BORDER_REPLICATE);
+
+		cv::Mat morphologyDisk = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(5, 5));
+		cv::morphologyEx(mPupil, mPupilOpened, cv::MORPH_OPEN, morphologyDisk, cv::Point(-1,-1), 2);
+
+		if (params.CannyBlur > 0)
+		{
+			cv::GaussianBlur(mPupilOpened, mPupilBlurred, cv::Size(), params.CannyBlur);
+		}
+		else
+		{
+			mPupilBlurred = mPupilOpened;
+		}
+
+		cv::Sobel(mPupilBlurred, mPupilSobelX, CV_32F, 1, 0, 3);
+		cv::Sobel(mPupilBlurred, mPupilSobelY, CV_32F, 0, 1, 3);
+
+		cv::Canny(mPupilBlurred, mPupilEdges, params.CannyThreshold1, params.CannyThreshold2);
+
+		cv::Rect roiUnpadded(padding,padding,roiPupil.width,roiPupil.height);
+		mPupil = cv::Mat(mPupil, roiUnpadded);
+		mPupilOpened = cv::Mat(mPupilOpened, roiUnpadded);
+		mPupilBlurred = cv::Mat(mPupilBlurred, roiUnpadded);
+		mPupilSobelX = cv::Mat(mPupilSobelX, roiUnpadded);
+		mPupilSobelY = cv::Mat(mPupilSobelY, roiUnpadded);
+		mPupilEdges = cv::Mat(mPupilEdges, roiUnpadded);
+
+		bbPupil = cvx::boundingBox(mPupil);
+	}
+
+	out.roiPupil = roiPupil;
+	out.mPupil = mPupil;
+	out.mPupilOpened = mPupilOpened;
+	out.mPupilBlurred = mPupilBlurred;
+	out.mPupilSobelX = mPupilSobelX;
+	out.mPupilSobelY = mPupilSobelY;
+	out.mPupilEdges = mPupilEdges;
+
+	// -----------------------------------------------
+	// Get points on edges, optionally using starburst
+	// -----------------------------------------------
+
+	std::vector<cv::Point2f> edgePoints;
+
+	if (params.StarburstPoints > 0)
+	{
+		SECTION("Starburst", log)
+		{
+			// Starburst from initial pupil approximation, stopping when an edge is hit.
+			// Collect all edge points into a vector
+
+			// The initial pupil approximations are:
+			//    Centre of mass of thresholded region
+			//    Halfway along the major axis (calculated form second moments) in each direction
+
+			tbb::concurrent_vector<cv::Point2f> edgePointsConcurrent;
+
+			cv::Vec2f elPupil_majorAxis = cvx::majorAxis(elPupilThresh);
+			std::vector<cv::Point2f> centres;
+			centres.push_back(elPupilThresh.center - cv::Point2f(roiPupil.tl().x, roiPupil.tl().y));
+			centres.push_back(elPupilThresh.center - cv::Point2f(roiPupil.tl().x, roiPupil.tl().y) + cv::Point2f(elPupil_majorAxis));
+			centres.push_back(elPupilThresh.center - cv::Point2f(roiPupil.tl().x, roiPupil.tl().y) - cv::Point2f(elPupil_majorAxis));
+
+			BOOST_FOREACH(const cv::Point2f& centre, centres) {
+				tbb::parallel_for(0, params.StarburstPoints, [&] (int i) {
+					double theta = i * 2*PI/params.StarburstPoints;
+
+					// Initialise centre and direction vector
+					cv::Point2f pDir((float)std::cos(theta), (float)std::sin(theta));  
+
+					int t = 1;
+					cv::Point p = centre + (t * pDir);
+					while(p.inside(bbPupil))
+					{
+						uchar val = mPupilEdges(p);
+
+						if (val > 0)
+						{
+							float dx = mPupilSobelX(p);
+							float dy = mPupilSobelY(p);
+
+							float cdirx = p.x - (elPupilThresh.center.x - roiPupil.x);
+							float cdiry = p.y - (elPupilThresh.center.y - roiPupil.y);
+
+							// Check edge direction
+							double dirCheck = dx*cdirx + dy*cdiry;
+
+							if (dirCheck > 0)
+							{
+								// We've hit an edge
+								edgePointsConcurrent.push_back(cv::Point2f(p.x + 0.5f, p.y + 0.5f));
+								break;
+							}
+						}
+
+						++t;
+						p = centre + (t * pDir);
+					}
+				});
+			}
+
+			edgePoints = std::vector<cv::Point2f>(edgePointsConcurrent.begin(), edgePointsConcurrent.end());
+
+
+			// Remove duplicate edge points
+			std::sort(edgePoints.begin(), edgePoints.end(), [] (const cv::Point2f& p1, const cv::Point2f& p2) -> bool {
+				if (p1.x == p2.x)
+					return p1.y < p2.y;
+				else
+					return p1.x < p2.x;
+			});
+			edgePoints.erase( std::unique( edgePoints.begin(), edgePoints.end() ), edgePoints.end() );
+
+			if (edgePoints.size() < params.StarburstPoints/2)
+				return false;
+		}
+	}
+	else
+	{
+		SECTION("Non-zero value finder", log)
+		{
+			for(int y = 0; y < mPupilEdges.rows; y++)
+			{
+				uchar* val = mPupilEdges[y];
+				for(int x = 0; x < mPupilEdges.cols; x++, val++)
+				{
+					if(*val == 0)
+						continue;
+
+					edgePoints.push_back(cv::Point2f(x + 0.5f, y + 0.5f));
+				}
+			}
+		}
+	}
+
+
+	// ---------------------------
+	// Fit an ellipse to the edges
+	// ---------------------------
+
+	cv::RotatedRect elPupil;
+	std::vector<cv::Point2f> inliers;
+	SECTION("Ellipse fitting", log)
+	{
+		// Desired probability that only inliers are selected
+		const double p = 0.999;
+		// Probability that a point is an inlier
+		double w = params.PercentageInliers/100.0;
+		// Number of points needed for a model
+		const int n = 5;
+
+		if (params.PercentageInliers == 0)
+			return false;
+
+		if (edgePoints.size() >= n) // Minimum points for ellipse
+		{
+			// RANSAC!!!
+
+			double wToN = std::pow(w,n);
+			int k = static_cast<int>(std::log(1-p)/std::log(1 - wToN)  + 2*std::sqrt(1 - wToN)/wToN);
+
+			out.ransacIterations = k;
+
+			log.add("k", k);
+
+			//size_t threshold_inlierCount = std::max<size_t>(n, static_cast<size_t>(out.edgePoints.size() * 0.7));
+
+			// Use TBB for RANSAC
+			struct EllipseRansac_out {
+				std::vector<cv::Point2f> bestInliers;
+				cv::RotatedRect bestEllipse;
+				double bestEllipseGoodness;
+				int earlyRejections;
+				bool earlyTermination;
+
+				EllipseRansac_out() : bestEllipseGoodness(-std::numeric_limits<double>::infinity()), earlyTermination(false), earlyRejections(0) {}
+			};
+			struct EllipseRansac {
+				const TrackerParams& params;
+				const std::vector<cv::Point2f>& edgePoints;
+				int n;
+				const cv::Rect& bb;
+				const cv::Mat_<float>& mDX;
+				const cv::Mat_<float>& mDY;
+				int earlyRejections;
+				bool earlyTermination;
+
+				EllipseRansac_out out;
+
+				EllipseRansac(
+					const TrackerParams& params,
+					const std::vector<cv::Point2f>& edgePoints,
+					int n,
+					const cv::Rect& bb,
+					const cv::Mat_<float>& mDX,
+					const cv::Mat_<float>& mDY)
+					: params(params), edgePoints(edgePoints), n(n), bb(bb), mDX(mDX), mDY(mDY), earlyTermination(false), earlyRejections(0)
+				{
+				}
+
+				EllipseRansac(EllipseRa