fvscanny.cpp

本系统是一个指纹识别系统,可以用来识别用户的指纹,可以用在安全性高的行业或部门.

// Taken from the canny.c file
// Robot Vision Group
// Dept. of Artifical Intelligence
// University of Edinburgh
//
// Authors: Bob Fisher, Dave Craft, A Fitzgibbion
// orginal date September 86
// program: canny.c
//
// purpose: to apply a simple symmetric canny operator to any size image
// 
// Revised September 2002
// By Shivang Patel
// Revisions made to work with Fingerprint Verification System (FVS)

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "fvsimage.h"
#include "vector.h"
#include "fvstypes.h"


/*
double hypotenuse(double x, double y);
double gaussian(double x, double s);
FvsError_t canny_core(double stddev, FvsImage<FVS_IMAGE_TYPE> & img, Vector<unsigned char> & derivative_mag, 
					  Vector<unsigned char> & magnitude, Vector <unsigned char> & orientation);
int follow(int i, int j, int low, FvsImage<FVS_IMAGE_TYPE> & img, Vector<unsigned char> & magnitude, Vector<unsigned char> & orientation);
*/

// Scales angles in radians to fit within the 0-255 range of unsigned
// char variables
#define ORIENT_SCALE 40.0

// hypot can't cope when both it's arguments are zero, hence this hack...
double hypotenuse(double x, double y)
{
	if (x == 0.0 && y == 0.0) return 0.0;
	else return hypot(x,y);
}

// gaussian function (centered at the origin and ignoring the factor of 
// 1/(s * sqrt(2 *PI))) 
double gaussian (double x, double s)
{
	return (exp((-x*x)/(2*s*s)));
}

FvsError_t canny_core(double stddev, FvsImage<FVS_IMAGE_TYPE> & img, Vector<unsigned char> & derivative_mag, 
					  Vector<unsigned char> & magnitude, Vector <unsigned char> & orientation)
{
	int filter_width; // length of 1-D gaussian mask
	Vector<float> gsmooth_x, gsmooth_y;
	Vector<float> derivative_x, derivative_y;
	int i, j, k, n; // counter
	int t; // temp grad magnitude variable
	double a, b, c, d, g0; // mask generation intermediate vars
	double ux, uy;
	double t1, t2;
	double G[20], dG[20], D2G[20]; // gaussian & derivative filter masks
	double gc, gn, gs, gw, ge, gnw, gne, gsw, gse; 
	int jlimit, jstart;
	int ilimit;
	int jfactor;
	int kfactor1, kfactor2;
	int kfactor;
	int cindex, nindex, sindex, windex, eindex, nwindex, neindex, swindex, seindex;
	int low = 1, high = 255; // tracker hysteresis parameters
	int col_plus, col_minus; // col + 1 and col-1 respectiably
	int picsize;
	int mag_overflow_count = 0; // used to measure how often mag array overflows
	int man_underflow_count = 0; // used to measure how often mag array underflow

	picsize = img.width() * img.height();

	// calculate coefficents for 1-dimensional G, dG/dn and delta squared G filters
	for (n = 0; n < 20; n++)
	{
		a = gaussian(((double)n), stddev);
		if (a > 0.005 || n<2)
		{
			b = gaussian((double)n-0.5, stddev);
			c = gaussian((double)n+0.5, stddev);
			d = gaussian((double)n, stddev);

			G[n] = (a + b + c)/3/(6.283185*stddev*stddev);
			dG[n] = c - b;
			D2G[n] = 1.6 * d - a; // DOG
		}
		else break;
	} // for n

	filter_width = n;
	
	// allocate space for gaussian smoothing arrays
	gsmooth_x.resize(picsize);
	gsmooth_y.resize(picsize);

	// produce x and y convolutions with gaussian
	ilimit = img.width() - (filter_width - 1);
	jstart = img.width() - (filter_width - 1);
	jlimit = img.width() - (img.height() - (filter_width - 1));

	for (i = filter_width - 1; i < ilimit; ++i)
	{
		for (jfactor = jstart; jstart < jlimit; jfactor += img.width())
		{
			cindex = i + jfactor;
			t1 = img[i][cindex] * G[0];
			t2 = t1;

			for (k = 1, kfactor = cindex - img.width(), kfactor2 = cindex + img.width();
			k < filter_width; k++, kfactor1 -= img.width(), kfactor2 += img.width())
			{
				t1 += G[k] * (img[kfactor1 - (kfactor1/img.width())][kfactor1/img.width()] +
					img[kfactor2 - (kfactor2/img.width())][kfactor2/img.width()]);
				t2 += G[k] * (img[(cindex - k) - ((cindex - k)/img.width())][((cindex - k)/
					img.width())] + img[(cindex + k) - ((cindex + k)/img.width())][((cindex +k)/img.width())]);
			} // for k

			gsmooth_x[cindex] = t1;
			gsmooth_y[cindex] = t2;

		} // for jfactor
	} // for i

	// allocate space for gradient arrays 
	derivative_x.resize(picsize);
	derivative_y.resize(picsize);

	// produce x and y convolutions with derivative of gaussian
	for (i = filter_width - 1; i < ilimit; ++i)
	{
		for (jfactor = jstart; jfactor < jlimit; jfactor += img.width())
		{
			t1 = 0;
			cindex = i + jfactor;
			for (k = 1; k < filter_width; ++k)
				t1 += dG[k] * (gsmooth_x[cindex - k] - gsmooth_x[cindex + k]);
			derivative_x[cindex] = t1;
		} // for jfactor
	} // for i

	for (i = n; i < img.width() - n; ++i)
	{
		for (jfactor = jstart; jfactor < jlimit; jfactor += img.width())
		{
			t2 = 0; 
			cindex = i + jfactor;
			for (k = 1, kfactor = img.width(); k < filter_width; k++, kfactor += img.width())
				t2 += dG[k] * (gsmooth_y[cindex - kfactor] - gsmooth_y[cindex + kfactor]);
			derivative_y[cindex] = t2;
		} // for jfactor 
	} // for i

	// non-maximum suppression (4 cases for orientation of line of max slope)
	ilimit = img.width() - filter_width;
	jstart = img.width() - filter_width;
	jlimit = img.width() * (img.height() - filter_width);
	
	for (i = filter_width; i < ilimit; ++i)
	{
		for (jfactor = jstart; jfactor < jlimit; jfactor += img.width())
		{
			// calculate current indeces
			cindex = i + jfactor;
			nindex = cindex - img.width();
			sindex = cindex + img.width();
			windex = cindex - 1;
			eindex = cindex + 1;
			nwindex = nindex - 1;
			neindex = nindex + 1;
			swindex = sindex - 1;
			seindex = sindex + 1;
			ux = derivative_x[cindex];
			uy = derivative_y[cindex];
			gc = hypotenuse(ux, uy);
			
			nwindex = (nwindex<0)?0:nwindex;

			// scale fv to fit into an unsighed change arry
			t = gc * 20.0;
			
			derivative_mag[cindex] = (t < 256?t:255);
			gn = hypotenuse(derivative_x[nindex], derivative_y[nindex]);
			gs = hypotenuse(derivative_x[sindex], derivative_y[sindex]);
			gw = hypotenuse(derivative_x[windex], derivative_y[windex]);
			ge = hypotenuse(derivative_x[eindex], derivative_y[eindex]);
			gne = hypotenuse(derivative_x[neindex], derivative_y[neindex]);
			gnw = hypotenuse(derivative_x[nwindex], derivative_y[nwindex]);
			gsw = hypotenuse(derivative_x[swindex], derivative_y[swindex]);
			gse = hypotenuse(derivative_x[seindex], derivative_y[seindex]);

			if (ux * uy > 0)
			{
				if (fabs(ux) < fabs(uy))
				{
					if ((g0 = fabs(uy * gc)) < fabs(ux * gse + (uy - ux) * gs) || 
						g0 <= fabs(ux * gnw + (uy - ux)*gn))
						continue;
				} else {
					if ((g0 = fabs(ux * gc)) < fabs(uy * gse + (ux - uy) * gc) || 
						g0 <= fabs(uy * gnw + (ux - uy) * gw))
						continue;
				} 
			} else {
				if (fabs(ux) < fabs(uy))
				{
					if ((g0 = fabs(uy * gc)) < fabs(ux * gne - (uy + ux) * gn) || 
						g0 <= fabs(ux * gsw - (uy + ux) * gs))
						continue;
				} else {
					if ((g0 = fabs(ux * gc)) < fabs(uy * gne - (ux + uy) * gc) ||
						g0 <= fabs(uy * gsw - (ux + uy) * gw))
						continue;
				}
			}

			// seems to be a good scale factor
			magnitude [cindex] = derivative_mag[cindex];

			// pi * 40 ~= 128 - direction is (thought of as ) a sign byte
			orientation[cindex] = atan2(ux, uy) * ORIENT_SCALE;
		}
	}
	
	//SP: vector class destructor will free all memory when function ends
	return FvsOk;
}

// follow a connect edge 
int follow(int i, int j, int low, FvsImage<FVS_IMAGE_TYPE> & img, Vector<unsigned char> & magnitude, Vector<unsigned char> & orientation)
{
	int k, l; // counters
	int i_plus_1, i_minus_1, j_plus_1, j_minus_1;
	long index, kindex;
	char break_flag;

	i_plus_1 = i + 1;
	i_minus_1 = i - 1;
	j_plus_1 = j + 1;
	j_minus_1 = j - 1;
	index = i + j * img.width();

	if (j_plus_1 >= img.height()) j_plus_1 = img.height() - 1;
	if (j_minus_1 < 0) j_minus_1 = 0;
	if (i_plus_1 >= img.height()) i_plus_1 = img.height() - 1;
	if (i_minus_1 < 0) i_minus_1 = 0;

	if (!img[i][j])
	{
		// current point must be added to the list of tracked points
		img[i][j] = magnitude[j * img.width() + i];

		// now we check immediately adhacent points to see if they too 
		// could be added to the track
		break_flag = 0;
		for (k = i_minus_1; k <= i_plus_1; k++)
		{
			for (l = j_minus_1; l <= j_plus_1; ++l)
			{
				kindex = k+l * img.width();
				// SP: Threshold low value should be replaced by 1
				if (!(l == j && k == i) && magnitude[kindex] >= 1) 
				{
					if (follow(k, l, low, img,magnitude,orientation)) 
					{
						break_flag = 1;
						break;
					} // if follow
				} // if (!(l==j && k==i)
			} // for l
			if (break_flag) break;
		} // for k
		return 1;
	} 
	else return 0;
}

// Track the edges in the magnitude away, starting at points that exceed the "high"
// threshold and continuing to track as long as point values don't duck below the "low"
// threshold (yes it's hystercsis... I'm sure that hyster means "womb" (as in hysterical),
// but I don't know what its doin in a common engineering term)
void thresholding_tracker(int high, int low, FvsImage<FVS_IMAGE_TYPE> & img, Vector<unsigned char> & magnitude,
						  Vector<unsigned char> & orientation)
{
	int i, j, k; // counters

	// clear data array before tracking
/*	for (i = 0; i < img.height(); i++)
	{
		for (j = 0; j < img.width(); j++)
		{
			img[i][j] = 0;
		}
	}
*/
	// threshold image with hystersis: follow from each strong edge point
	for (i = 0; i < img.width(); i++)
	{
		for (j = 0; j < img.height(); j++)
		{
			if (magnitude[i + img.width() * j] >= high)
				follow(i, j, low, img, magnitude, orientation);
		}
	}
}


int FVSCannyEdgeDect(FvsImage<FVS_IMAGE_TYPE> & img, int low=1, int high=255)
{
	double s = 1.0; // standard deviation of gaussian
	Vector<unsigned char> derivative_mag; // mag of del G before non-maximum suppression
	Vector<unsigned char> magnitude, orientation; // mag & orintation agter non-max suppression

	// resize vectors
	derivative_mag.resize(img.width() * img.height());
	magnitude.resize(img.width() * img.height());
	orientation.resize(img.width() * img.height());

	// call canny core routines - these perform a gaussian smoothing, calculate a gradient
	// array, and suppress non-maximal points in this array
	canny_core(s, img, derivative_mag, magnitude, orientation);



	// here write orientation and magnitude to a file to see how they look

	

	// track edges 
	thresholding_tracker(high, low, img, magnitude, orientation);

	// to see the image output img to a file
	return 0;
}