Commit f6e1ceff authored by Thiago Santini's avatar Thiago Santini

Adds PuRe

parent 0e1ddbaa
/*
* Copyright (c) 2018, Thiago Santini
*
* Permission to use, copy, modify, and distribute this software and its
* documentation for non-commercial purposes, without fee, and without a written
* agreement is hereby granted, provided that:
*
* 1) the above copyright notice, this permission notice, and the subsequent
* bibliographic references be included in all copies or substantial portions of
* the software
*
* 2) the appropriate bibliographic references be made on related publications
*
* In this context, non-commercial means not intended for use towards commercial
* advantage (e.g., as complement to or part of a product) or monetary
* compensation. The copyright holder reserves the right to decide whether a
* certain use classifies as commercial or not. For commercial use, please contact
* the copyright holders.
*
* REFERENCES:
*
* Thiago Santini, Wolfgang Fuhl, Enkelejda Kasneci, PuRe: Robust pupil detection
* for real-time pervasive eye tracking, Computer Vision and Image Understanding,
* 2018, ISSN 1077-3142, https://doi.org/10.1016/j.cviu.2018.02.002.
*
*
* IN NO EVENT SHALL THE AUTHORS BE LIABLE TO ANY PARTY FOR DIRECT,
* INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
* ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
* THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* THE AUTHORS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE AUTHORS
* HAVE NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
* MODIFICATIONS.
*/
#include "PuRe.h"
#include <climits>
#include <iostream>
#include <QDebug>
#include <QElapsedTimer>
#include <opencv2/highgui.hpp>
//#define SAVE_ILLUSTRATION
using namespace std;
using namespace cv;
string PuRe::desc = "PuRe (Santini et. al 2018)";
PuRe::PuRe() :
baseSize(320,240),
expectedFrameSize(-1,-1),
outlineBias(5)
{
mDesc = desc;
/*
* 1) Canthi:
* Using measurements from white men
* Mean intercanthal distance 32.7 (2.4) mm
* Mean palpebral fissure width 27.6 (1.9) mm
* Jayanth Kunjur, T. Sabesan, V. Ilankovan
* Anthropometric analysis of eyebrows and eyelids:
* An inter-racial study
*/
meanCanthiDistanceMM = 27.6f;
//meanCanthiDistanceMM = 32.7f;
/*
* 2) Pupil:
* 2 to 4 mm in diameter in bright light to 4 to 8 mm in the dark
* Clinical Methods: The History, Physical, and Laboratory Examinations. 3rd edition.
* Chapter 58The Pupils
* Robert H. Spector.
*/
maxPupilDiameterMM = 8.0f;
minPupilDiameterMM = 2.0f;
}
PuRe::~PuRe()
{
}
void PuRe::estimateParameters(int rows, int cols)
{
/*
* Assumptions:
* 1) The image contains at least both eye corners
* 2) The image contains a maximum of 5cm of the face (i.e., ~= 2x canthi distance)
*/
float d = sqrt( pow(rows,2) + pow(cols,2) );
maxCanthiDistancePx = d;
minCanthiDistancePx = 2*d/3.0;
maxPupilDiameterPx = maxCanthiDistancePx*(maxPupilDiameterMM/meanCanthiDistanceMM);
minPupilDiameterPx = minCanthiDistancePx*(minPupilDiameterMM/meanCanthiDistanceMM);
}
void PuRe::init(const Mat &frame)
{
if (expectedFrameSize == Size(frame.cols, frame.rows))
return;
expectedFrameSize = Size(frame.cols, frame.rows);
float rw = baseSize.width / (float) frame.cols;
float rh = baseSize.height / (float) frame.rows;
scalingRatio = min<float>( min<float>(rw, rh) , 1.0 );
}
Mat PuRe::canny(const Mat &in, bool blurImage, bool useL2, int bins, float nonEdgePixelsRatio, float lowHighThresholdRatio)
{
(void) useL2;
/*
* Smoothing and directional derivatives
* TODO: adapt sizes to image size
*/
Mat blurred;
if (blurImage) {
Size blurSize(5,5);
GaussianBlur(in, blurred, blurSize, 1.5, 1.5, BORDER_REPLICATE);
} else
blurred = in;
Sobel(blurred, dx, dx.type(), 1, 0, 7, 1, BORDER_REPLICATE);
Sobel(blurred, dy, dy.type(), 0, 1, 7, 1, BORDER_REPLICATE);
/*
* Magnitude
*/
double minMag = 0;
double maxMag = 0;
float *p_res;
float *p_x, *p_y; // result, x, y
cv::magnitude(dx, dy, magnitude);
cv::minMaxLoc(magnitude, &minMag, &maxMag);
/*
* Threshold selection based on the magnitude histogram
*/
float low_th = 0;
float high_th = 0;
// Normalization
magnitude = magnitude / maxMag;
// Histogram
int *histogram = new int[bins]();
Mat res_idx = (bins-1) * magnitude;
res_idx.convertTo(res_idx, CV_16U);
short *p_res_idx=0;
for(int i=0; i<res_idx.rows; i++){
p_res_idx = res_idx.ptr<short>(i);
for(int j=0; j<res_idx.cols; j++)
histogram[ p_res_idx[j] ]++;
}
// Ratio
int sum=0;
int nonEdgePixels = nonEdgePixelsRatio * in.rows * in.cols;
for(int i=0; i<bins; i++){
sum += histogram[i];
if( sum > nonEdgePixels ){
high_th = float(i+1) / bins ;
break;
}
}
low_th = lowHighThresholdRatio*high_th;
delete histogram;
/*
* Non maximum supression
*/
const float tg22_5 = 0.4142135623730950488016887242097f;
const float tg67_5 = 2.4142135623730950488016887242097f;
uchar *_edgeType;
float *p_res_b, *p_res_t;
edgeType.setTo(0);
for(int i=1; i<magnitude.rows-1; i++) {
_edgeType = edgeType.ptr<uchar>(i);
p_res=magnitude.ptr<float>(i);
p_res_t=magnitude.ptr<float>(i-1);
p_res_b=magnitude.ptr<float>(i+1);
p_x=dx.ptr<float>(i);
p_y=dy.ptr<float>(i);
for(int j=1; j<magnitude.cols-1; j++){
float m = p_res[j];
if (m < low_th)
continue;
float iy = p_y[j];
float ix = p_x[j];
float y = abs( (double) iy );
float x = abs( (double) ix );
uchar val = p_res[j] > high_th ? 255 : 128;
float tg22_5x = tg22_5 * x;
if (y < tg22_5x) {
if (m > p_res[j-1] && m >= p_res[j+1])
_edgeType[j] = val;
} else {
float tg67_5x = tg67_5 * x;
if (y > tg67_5x) {
if (m > p_res_b[j] && m >= p_res_t[j])
_edgeType[j] = val;
} else {
if ( (iy<=0) == (ix<=0) ) {
if ( m > p_res_t[j-1] && m >= p_res_b[j+1])
_edgeType[j] = val;
} else {
if ( m > p_res_b[j-1] && m >= p_res_t[j+1])
_edgeType[j] = val;
}
}
}
}
}
/*
* Hystheresis
*/
int pic_x=edgeType.cols;
int pic_y=edgeType.rows;
int area = pic_x*pic_y;
int lines_idx=0;
int idx=0;
vector<int> lines;
edge.setTo(0);
for(int i=1;i<pic_y-1;i++){
for(int j=1;j<pic_x-1;j++){
if( edgeType.data[idx+j] != 255 || edge.data[idx+j] != 0 )
continue;
edge.data[idx+j] = 255;
lines_idx = 1;
lines.clear();
lines.push_back(idx+j);
int akt_idx = 0;
while(akt_idx<lines_idx){
int akt_pos=lines[akt_idx];
akt_idx++;
if( akt_pos-pic_x-1 < 0 || akt_pos+pic_x+1 >= area )
continue;
for(int k1=-1;k1<2;k1++)
for(int k2=-1;k2<2;k2++){
if(edge.data[(akt_pos+(k1*pic_x))+k2]!=0 || edgeType.data[(akt_pos+(k1*pic_x))+k2]==0)
continue;
edge.data[(akt_pos+(k1*pic_x))+k2] = 255;
lines.push_back((akt_pos+(k1*pic_x))+k2);
lines_idx++;
}
}
}
idx+=pic_x;
}
return edge;
}
void PuRe::filterEdges(cv::Mat &edges)
{
// TODO: there is room for improvement here; however, it is prone to small
// mistakes; will be done when we have time
int start_x = 5;
int start_y = 5;
int end_x = edges.cols - 5;
int end_y = edges.rows - 5;
for(int j=start_y; j<end_y; j++)
for(int i=start_x; i<end_x; i++){
uchar box[9];
box[4]=(uchar)edges.data[(edges.cols*(j))+(i)];
if(box[4]){
box[1]=(uchar)edges.data[(edges.cols*(j-1))+(i)];
box[3]=(uchar)edges.data[(edges.cols*(j))+(i-1)];
box[5]=(uchar)edges.data[(edges.cols*(j))+(i+1)];
box[7]=(uchar)edges.data[(edges.cols*(j+1))+(i)];
if((box[5] && box[7])) edges.data[(edges.cols*(j))+(i)]=0;
if((box[5] && box[1])) edges.data[(edges.cols*(j))+(i)]=0;
if((box[3] && box[7])) edges.data[(edges.cols*(j))+(i)]=0;
if((box[3] && box[1])) edges.data[(edges.cols*(j))+(i)]=0;
}
}
//too many neigbours
for(int j=start_y; j<end_y; j++)
for(int i=start_x; i<end_x; i++){
uchar neig=0;
for(int k1=-1;k1<2;k1++)
for(int k2=-1;k2<2;k2++){
if(edges.data[(edges.cols*(j+k1))+(i+k2)]>0)
neig++;
}
if(neig>3)
edges.data[(edges.cols*(j))+(i)]=0;
}
for(int j=start_y; j<end_y; j++)
for(int i=start_x; i<end_x; i++){
uchar box[17];
box[4]=(uchar)edges.data[(edges.cols*(j))+(i)];
if(box[4]){
box[0]=(uchar)edges.data[(edges.cols*(j-1))+(i-1)];
box[1]=(uchar)edges.data[(edges.cols*(j-1))+(i)];
box[2]=(uchar)edges.data[(edges.cols*(j-1))+(i+1)];
box[3]=(uchar)edges.data[(edges.cols*(j))+(i-1)];
box[5]=(uchar)edges.data[(edges.cols*(j))+(i+1)];
box[6]=(uchar)edges.data[(edges.cols*(j+1))+(i-1)];
box[7]=(uchar)edges.data[(edges.cols*(j+1))+(i)];
box[8]=(uchar)edges.data[(edges.cols*(j+1))+(i+1)];
//external
box[9]=(uchar)edges.data[(edges.cols*(j))+(i+2)];
box[10]=(uchar)edges.data[(edges.cols*(j+2))+(i)];
box[11]=(uchar)edges.data[(edges.cols*(j))+(i+3)];
box[12]=(uchar)edges.data[(edges.cols*(j-1))+(i+2)];
box[13]=(uchar)edges.data[(edges.cols*(j+1))+(i+2)];
box[14]=(uchar)edges.data[(edges.cols*(j+3))+(i)];
box[15]=(uchar)edges.data[(edges.cols*(j+2))+(i-1)];
box[16]=(uchar)edges.data[(edges.cols*(j+2))+(i+1)];
if( (box[10] && !box[7]) && (box[8] || box[6]) ){
edges.data[(edges.cols*(j+1))+(i-1)]=0;
edges.data[(edges.cols*(j+1))+(i+1)]=0;
edges.data[(edges.cols*(j+1))+(i)]=255;
}
if( (box[14] && !box[7] && !box[10]) && ( (box[8] || box[6]) && (box[16] || box[15]) ) ){
edges.data[(edges.cols*(j+1))+(i+1)]=0;
edges.data[(edges.cols*(j+1))+(i-1)]=0;
edges.data[(edges.cols*(j+2))+(i+1)]=0;
edges.data[(edges.cols*(j+2))+(i-1)]=0;
edges.data[(edges.cols*(j+1))+(i)]=255;
edges.data[(edges.cols*(j+2))+(i)]=255;
}
if( (box[9] && !box[5]) && (box[8] || box[2]) ){
edges.data[(edges.cols*(j+1))+(i+1)]=0;
edges.data[(edges.cols*(j-1))+(i+1)]=0;
edges.data[(edges.cols*(j))+(i+1)]=255;
}
if( (box[11] && !box[5] && !box[9]) && ( (box[8] || box[2]) && (box[13] || box[12]) ) ){
edges.data[(edges.cols*(j+1))+(i+1)]=0;
edges.data[(edges.cols*(j-1))+(i+1)]=0;
edges.data[(edges.cols*(j+1))+(i+2)]=0;
edges.data[(edges.cols*(j-1))+(i+2)]=0;
edges.data[(edges.cols*(j))+(i+1)]=255;
edges.data[(edges.cols*(j))+(i+2)]=255;
}
}
}
for(int j=start_y; j<end_y; j++)
for(int i=start_x; i<end_x; i++){
uchar box[33];
box[4]=(uchar)edges.data[(edges.cols*(j))+(i)];
if(box[4]){
box[0]=(uchar)edges.data[(edges.cols*(j-1))+(i-1)];
box[1]=(uchar)edges.data[(edges.cols*(j-1))+(i)];
box[2]=(uchar)edges.data[(edges.cols*(j-1))+(i+1)];
box[3]=(uchar)edges.data[(edges.cols*(j))+(i-1)];
box[5]=(uchar)edges.data[(edges.cols*(j))+(i+1)];
box[6]=(uchar)edges.data[(edges.cols*(j+1))+(i-1)];
box[7]=(uchar)edges.data[(edges.cols*(j+1))+(i)];
box[8]=(uchar)edges.data[(edges.cols*(j+1))+(i+1)];
box[9]=(uchar)edges.data[(edges.cols*(j-1))+(i+2)];
box[10]=(uchar)edges.data[(edges.cols*(j-1))+(i-2)];
box[11]=(uchar)edges.data[(edges.cols*(j+1))+(i+2)];
box[12]=(uchar)edges.data[(edges.cols*(j+1))+(i-2)];
box[13]=(uchar)edges.data[(edges.cols*(j-2))+(i-1)];
box[14]=(uchar)edges.data[(edges.cols*(j-2))+(i+1)];
box[15]=(uchar)edges.data[(edges.cols*(j+2))+(i-1)];
box[16]=(uchar)edges.data[(edges.cols*(j+2))+(i+1)];
box[17]=(uchar)edges.data[(edges.cols*(j-3))+(i-1)];
box[18]=(uchar)edges.data[(edges.cols*(j-3))+(i+1)];
box[19]=(uchar)edges.data[(edges.cols*(j+3))+(i-1)];
box[20]=(uchar)edges.data[(edges.cols*(j+3))+(i+1)];
box[21]=(uchar)edges.data[(edges.cols*(j+1))+(i+3)];
box[22]=(uchar)edges.data[(edges.cols*(j+1))+(i-3)];
box[23]=(uchar)edges.data[(edges.cols*(j-1))+(i+3)];
box[24]=(uchar)edges.data[(edges.cols*(j-1))+(i-3)];
box[25]=(uchar)edges.data[(edges.cols*(j-2))+(i-2)];
box[26]=(uchar)edges.data[(edges.cols*(j+2))+(i+2)];
box[27]=(uchar)edges.data[(edges.cols*(j-2))+(i+2)];
box[28]=(uchar)edges.data[(edges.cols*(j+2))+(i-2)];
box[29]=(uchar)edges.data[(edges.cols*(j-3))+(i-3)];
box[30]=(uchar)edges.data[(edges.cols*(j+3))+(i+3)];
box[31]=(uchar)edges.data[(edges.cols*(j-3))+(i+3)];
box[32]=(uchar)edges.data[(edges.cols*(j+3))+(i-3)];
if( box[7] && box[2] && box[9] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[7] && box[0] && box[10] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[1] && box[8] && box[11] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[1] && box[6] && box[12] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[0] && box[13] && box[17] && box[8] && box[11] && box[21] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[2] && box[14] && box[18] && box[6] && box[12] && box[22] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[6] && box[15] && box[19] && box[2] && box[9] && box[23] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[8] && box[16] && box[20] && box[0] && box[10] && box[24] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[0] && box[25] && box[2] && box[27] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[0] && box[25] && box[6] && box[28] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[8] && box[26] && box[2] && box[27] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box[8] && box[26] && box[6] && box[28] )
edges.data[(edges.cols*(j))+(i)]=0;
uchar box2[18];
box2[1]=(uchar)edges.data[(edges.cols*(j))+(i-1)];
box2[2]=(uchar)edges.data[(edges.cols*(j-1))+(i-2)];
box2[3]=(uchar)edges.data[(edges.cols*(j-2))+(i-3)];
box2[4]=(uchar)edges.data[(edges.cols*(j-1))+(i+1)];
box2[5]=(uchar)edges.data[(edges.cols*(j-2))+(i+2)];
box2[6]=(uchar)edges.data[(edges.cols*(j+1))+(i-2)];
box2[7]=(uchar)edges.data[(edges.cols*(j+2))+(i-3)];
box2[8]=(uchar)edges.data[(edges.cols*(j+1))+(i+1)];
box2[9]=(uchar)edges.data[(edges.cols*(j+2))+(i+2)];
box2[10]=(uchar)edges.data[(edges.cols*(j+1))+(i)];
box2[15]=(uchar)edges.data[(edges.cols*(j-1))+(i-1)];
box2[16]=(uchar)edges.data[(edges.cols*(j-2))+(i-2)];
box2[11]=(uchar)edges.data[(edges.cols*(j+2))+(i+1)];
box2[12]=(uchar)edges.data[(edges.cols*(j+3))+(i+2)];
box2[13]=(uchar)edges.data[(edges.cols*(j+2))+(i-1)];
box2[14]=(uchar)edges.data[(edges.cols*(j+3))+(i-2)];
if( box2[1] && box2[2] && box2[3] && box2[4] && box2[5] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box2[1] && box2[6] && box2[7] && box2[8] && box2[9] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box2[10] && box2[11] && box2[12] && box2[4] && box2[5] )
edges.data[(edges.cols*(j))+(i)]=0;
if( box2[10] && box2[13] && box2[14] && box2[15] && box2[16] )
edges.data[(edges.cols*(j))+(i)]=0;
}
}
}
void PuRe::findPupilEdgeCandidates(const Mat &intensityImage, Mat &edge, vector<PupilCandidate> &candidates)
{
/* Find all lines
* Small note here: using anchor points tends to result in better ellipse fitting later!
* It's also faster than doing connected components and collecting the labels
*/
vector<Vec4i> hierarchy;
vector<vector<Point> > curves;
findContours( edge, curves, hierarchy, CV_RETR_LIST, CV_CHAIN_APPROX_TC89_KCOS );
removeDuplicates(curves, edge.cols);
// Create valid candidates
for (size_t i=curves.size(); i-->0;) {
PupilCandidate candidate(curves[i]);
if (candidate.isValid(intensityImage, minPupilDiameterPx, maxPupilDiameterPx, outlineBias))
candidates.push_back( candidate );
}
}
void PuRe::combineEdgeCandidates(const cv::Mat &intensityImage, cv::Mat &edge, std::vector<PupilCandidate> &candidates)
{
(void) edge;
if (candidates.size() <= 1)
return;
vector<PupilCandidate> mergedCandidates;
for (auto pc=candidates.begin(); pc!=candidates.end(); pc++) {
for (auto pc2=pc+1; pc2!=candidates.end(); pc2++) {
Rect intersection = pc->combinationRegion & pc2->combinationRegion;
if (intersection.area() < 1)
continue; // no intersection
//#define DBG_EDGE_COMBINATION
#ifdef DBG_EDGE_COMBINATION
Mat tmp;
cvtColor(intensityImage, tmp, CV_GRAY2BGR);
rectangle(tmp, pc->combinationRegion, pc->color);
for (unsigned int i=0; i<pc->points.size(); i++)
cv::circle(tmp, pc->points[i], 1, pc->color, -1);
rectangle(tmp, pc2->combinationRegion, pc2->color);
for (unsigned int i=0; i<pc2->points.size(); i++)
cv::circle(tmp, pc2->points[i], 1, pc2->color, -1);
imshow("Combined edges", tmp);
imwrite("combined.png", tmp);
//waitKey(0);
#endif
if (intersection.area() >= min<int>(pc->combinationRegion.area(),pc2->combinationRegion.area()))
continue;
vector<Point> mergedPoints = pc->points;
mergedPoints.insert(mergedPoints.end(), pc2->points.begin(), pc2->points.end());
PupilCandidate candidate( mergedPoints );
if (!candidate.isValid(intensityImage, minPupilDiameterPx, maxPupilDiameterPx, outlineBias))
continue;
if (candidate.outlineContrast < pc->outlineContrast || candidate.outlineContrast < pc2->outlineContrast)
continue;
mergedCandidates.push_back( candidate );
}