harrislike.cc

/*

    This file is part of the FAST-ER machine learning system.
    Copyright (C) 2008  Edward Rosten and Los Alamos National Laboratory

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
#include <cvd/image.h>
#include <cvd/convolution.h>
#include <gvars3/instances.h>
#include <vector>

#include "harrislike.h"

using namespace std;
using namespace CVD;
using namespace GVars3;



////////////////////////////////////////////////////////////////////////////////
//
// Harris-like corners.
// 

///\cond never
template<class C> inline C sq(const C& c)
{
	return c*c;
}


struct HarrisScore
{
	static float Compute(float xx, float xy, float yy)
	{
		return (xx * yy - xy * xy) - 0.04 * sq(xx + yy);
	}
};

struct ShiTomasiScore
{
	static float Compute(float xx, float xy, float yy)
	{
		float l1 = xx + yy + sqrt(sq(xx - yy)+4.0*xy*xy);   
		float l2 = xx + yy - sqrt(sq(xx - yy)+4.0*xy*xy);
		return min(abs(l1), abs(l2));
	}
};

struct PosInserter
{
	static void insert(vector<ImageRef>& i, const pair<float, ImageRef>& p)
	{
		i.push_back(p.second);
	}
};

struct PairInserter
{
	static void insert(vector<pair<float, ImageRef> >& i, const pair<float, ImageRef>& p)
	{
		i.push_back(p);
	}
};
///\endcond

template<class Score, class Inserter, class C, class B> void harris_like(const Image<B>& i, C& c, unsigned int N, float blur, float sigmas)
{

	Image<float> xx(i.size()), xy(i.size()), yy(i.size());

	zeroBorders(xx);
	zeroBorders(xy);
	zeroBorders(yy);

	typedef typename Pixel::traits<B>::wider_type gType;
	
	//Compute gradients
	for(int y=1; y < i.size().y - 1; y++)
		for(int x=1; x < i.size().x - 1; x++)
		{
			gType gx = (gType)i[y][x-1] - i[y][x+1];
			gType gy = (gType)i[y-1][x] - i[y+1][x];
			
			//Compute the gradient moments
			xx[y][x] = gx * gx;
			xy[y][x] = gx * gy;
			yy[y][x] = gy * gy;
		}

	convolveGaussian_fir(xx, xx, blur, sigmas);
	convolveGaussian_fir(xy, xy, blur, sigmas);
	convolveGaussian_fir(yy, yy, blur, sigmas);
	
	//Avoid computing the score along the image borders where the
	//result of the convolution is not valid.
	int	kspread = (int)ceil(sigmas * blur);

	//Compute harris score
	for(int y=kspread; y < i.size().y-kspread; y++)
		for(int x=kspread; x <i.size().x-kspread; x++)
			xx[y][x] = Score::Compute(xx[y][x], xy[y][x], yy[y][x]);

	vector<pair<float, ImageRef> > corners;
	corners.reserve(10000);

	//Find local maxima
	for(int y=kspread; y < i.size().y-kspread; y++)
		for(int x=kspread; x <i.size().x-kspread; x++)
		{
			float c = xx[y][x];

			if( c > xx[y-1][x-1]  &&
				c > xx[y-1][x+0]  &&
				c > xx[y-1][x+1]  &&
				c > xx[y-0][x-1]  &&
				c > xx[y-0][x+1]  &&
				c > xx[y+1][x-1]  &&
				c > xx[y+1][x+0]  &&
				c > xx[y+1][x+1])
			{
				corners.push_back(make_pair(-c, ImageRef(x, y)));
			}
		}

	if(corners.size() > N)
	{
		nth_element(corners.begin(), corners.begin() + N, corners.end());
		corners.resize(N);
	}


	for(unsigned int i=0; i < corners.size(); i++)
		Inserter::insert(c, corners[i]);
}


void HarrisDetector(const Image<float>& i, vector<pair<float, ImageRef> >& c, unsigned int N, float blur, float sigmas)
{
	harris_like<HarrisScore, PairInserter>(i, c, N, blur, sigmas);
}


void HarrisDetect::operator()(const CVD::Image<CVD::byte>& i, std::vector<CVD::ImageRef>& c, unsigned int N) const 
{
	float blur = GV3::get<float>("harris.blur", 2.5, 1);
	float sigmas = GV3::get<float>("harris.sigmas", 2.0, 1);
	harris_like<HarrisScore,PosInserter>(i, c, N, blur, sigmas);
}

void ShiTomasiDetect::operator()(const CVD::Image<CVD::byte>& i, std::vector<CVD::ImageRef>& c, unsigned int N) const
{
	float blur = GV3::get<float>("shitomasi.blur", 2.5, 1);
	float sigmas = GV3::get<float>("shitomasi.sigmas", 2.0, 1);
	harris_like<ShiTomasiScore, PosInserter>(i, c, N, blur, sigmas);
}