图形图像处理

开发平台：
Visual C++

ScaleSpace.cpp：源码内容
							// ScaleSpace.cpp: implementation of the CScaleSpace class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "TestSIFT.h"
#include "ScaleSpace.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
// Variable s is from explanation in Lowe paper
CScaleSpace::CScaleSpace(int s,double basePixScale)
{
	m_scaleLevel = s;
	m_dogCount = m_scaleLevel+2;	
	m_gaussCount = m_scaleLevel+3;
	m_gaussImgs = NULL;
	m_dogImgs = NULL;
	m_magnitudes = NULL;
	m_directions = NULL;
	m_basePixScale = basePixScale;
}
CScaleSpace::~CScaleSpace()
{
	if (m_dogImgs != NULL)
	{
		for (int i = 0; i < m_dogCount; i++)
			cvReleaseImage (&m_dogImgs[i]);
		delete m_dogImgs;
	}
	if (m_gaussImgs != NULL)
	{
		for (int i = 0; i < m_gaussCount; i++)
			cvReleaseImage (&m_gaussImgs[i]);
		delete m_gaussImgs;
	}
}
void CScaleSpace::buildDoG(IplImage *img,double startsigma,double firstScale)
{
	m_gaussImgs = new IplImage*[m_gaussCount];
	for (int i = 0; i < m_gaussCount; i++)
		m_gaussImgs[i] = cvCreateImage (cvGetSize(img),IPL_DEPTH_32F,1);
	if (img->depth == IPL_DEPTH_32F)
		cvCopy(img,m_gaussImgs[0]);
	else
	{
		for (int y = 0; y < cvGetSize(img).height; y++)
			for (int x = 0; x < cvGetSize(img).width; x++)
				cvSet2D (m_gaussImgs[0],y,x,cvGet2D(img,y,x));
	}
	m_dogImgs = new IplImage*[m_dogCount];
	for (i = 0; i < m_dogCount; i++)
		m_dogImgs[i] = cvCreateImage (cvGetSize(img),IPL_DEPTH_32F,1);
	
	/* This explaination is from Sebastian Nowozin's C# code
	 *
	 * Gaussian G(sigma), with relation
	 * G(sigma_1) * G(sigma_2) = G(sqrt(sigma_1^2 + * sigma_2^2))
	 *
	 * Then, we have:
	 *
	 * G(k^{p+1}) = G(k^p) * G(sigma),
	 * and our goal is to compute every iterations sigma value so this
	 * equation iteratively produces the next level.  Hence:
	 *
	 * sigma = sqrt{left(k^{p+1}right)^2 - left(k^pright)^2}
	 *   = sqrt{k^{2p+2} - k^{2p}}
	 *   = sqrt{k^2p (k^2 - 1)}
	 *   = k^p sqrt{k^2 - 1}
	 *
	 * In the code below, 'w' is the running k^p, where p increases by one
	 * each iteration.  kTerm is the constant sqrt{k^2 - 1} term.
	 */
	 // SToK (scales) = 2^(1/scales) == k
	double k = exp (1.0/m_scaleLevel*log(2));
	double kTerm = sqrt (k*k - 1.0);
	double w = startsigma;
	for (i = 1 ; i < m_gaussCount ; i++) {
		cvSmooth( m_gaussImgs[i-1], m_gaussImgs[i],CV_GAUSSIAN,0,0,w);
		cvAbsDiff (m_gaussImgs[i-1],m_gaussImgs[i],m_dogImgs[i-1]);
		w *= k;
	}		
}
IplImage* CScaleSpace::getDoG(int i)
{
	return m_dogImgs[i];
}
int CScaleSpace::getDoGCount()
{
	return m_dogCount;
}
void CScaleSpace::findPeak(double r,double dogThresh,std::vector<CScalePoint>* vec)
{
	// dogThresh is in range [0,1]
	dogThresh *= 255; 
	// {x,y,scale offset}
	int dir[27][3] = 
	{{-1,-1,0},{0,-1,0},{1,-1,0},{-1,0,0},{1,0,0},{-1,1,0},{0,1,0},{1,1,0},
	{-1,-1,-1},{0,-1,-1},{1,-1,-1},{-1,0,-1},{0,0,-1},{1,0,-1},{-1,1,-1},{0,1,-1},{1,1,-1},
	{-1,-1,1},{0,-1,1},{1,-1,1},{-1,0,1},{0,0,1},{1,0,1},{-1,1,1},{0,1,1},{1,1,1}};
	for (int i = 1; i < m_dogCount-1; i++)
	{
		for (int y = 1; y < cvGetSize (m_dogImgs[i]).height-1; y++)
			for (int x = 1; x < cvGetSize (m_dogImgs[i]).width-1; x++)
			{
				double val = cvGet2D (m_dogImgs[i],y,x).val[0];
				bool maxima = true;
				for (int k = 0; k < 27; k++)
				{
					if (val < dogThresh || val < cvGet2D (m_dogImgs[i+dir[k][2]],y+dir[k][1],x+dir[k][0]).val[0])
					{
						maxima = false;
						break;
					}
				}
				if (maxima)
				{
					// Eliminate edge response
					double dxx = cvGet2D(m_dogImgs[i],y,x+1).val[0] +
								 cvGet2D(m_dogImgs[i],y,x-1).val[0] -
								 2 * cvGet2D(m_dogImgs[i],y,x).val[0];
					double dyy = cvGet2D(m_dogImgs[i],y+1,x).val[0] +
								 cvGet2D(m_dogImgs[i],y-1,x).val[0] -
								 2 * cvGet2D(m_dogImgs[i],y,x).val[0];
					double dxy = 0.25 * (cvGet2D(m_dogImgs[i],y+1,x+1).val[0] -
										 cvGet2D(m_dogImgs[i],y-1,x+1).val[0] -
										 cvGet2D(m_dogImgs[i],y+1,x-1).val[0] +
										 cvGet2D(m_dogImgs[i],y-1,x-1).val[0]);					
					double tr = dxx+dyy; // trace 
					double det = dxx*dyy - dxy*dxy; //det					
					double rsq = r+1;
					rsq *= rsq;
					if (tr*tr/(det+1e-20) > (rsq)/r)
						maxima = false;
				}
				if (maxima)
				{
					CScalePoint p (x,y,i);
					vec->push_back(p);
				}
				//else
				//	cvSet2D (dogImg[i],y,x,cvScalar(0));
			}		
	}
}
void CScaleSpace::filterAndLocalizePeaks(std::vector<CScalePoint>& peaks,std::vector<CScalePoint>* filtered, 
										 double dValueLoThresh, double scaleAdjustThresh, int relocationMaximum)
{
	dValueLoThresh *= 255;
	CvMat *processedMap = cvCreateMat (cvGetSize(m_dogImgs[0]).height,cvGetSize(m_dogImgs[0]).width,CV_8UC1);
	cvSetZero (processedMap);
	for (int i = 0; i < peaks.size(); i++) 
	{
		// Kuas : already reject edge like in findPeak
		//if (IsTooEdgelike (spaces[peak.Level], peak.X, peak.Y, edgeRatio))
		//	continue;
		// When the localization hits some problem, i.e. while moving the
		// point a border is reached, then skip this point.
		if (LocalizeIsWeak (peaks[i], relocationMaximum, processedMap))
			continue;
		
		if (abs (peaks[i].m_fineS) > scaleAdjustThresh)
			continue;
		// Additional local pixel information is now available, threshhold
		// the D(^x)
		// Console.WriteLine ("{0} {1} {2} # == DVALUE", peak.Y, peak.X, peak.Local.DValue);
		if (abs (peaks[i].m_DValue) <= dValueLoThresh)
			continue;
	
		// its edgy enough, add it
		filtered->push_back (peaks[i]);
	}
	cvReleaseMat (&processedMap);
	//return (filtered);
}
bool CScaleSpace::LocalizeIsWeak(CScalePoint& peak, int steps, CvMat* processed)
{
	bool needToAdjust = true;
	int adjusted = steps;
	while (needToAdjust) 
	{
		int x = peak.m_x;
		int y = peak.m_y;
		// Points we cannot say anything about, as they lie on the border
		// of the scale space
		// Kuas : seem to be not happen because in FindPeaks not get these peaks
		//if (point.Level <= 0 || point.Level >= (spaces.Length - 1))
		//	return (true);
		IplImage* dog = m_dogImgs[peak.m_level];
		if (x <= 0 || x >= (cvGetSize(dog).width-1))
			return (true);
		if (y <= 0 || y >= (cvGetSize(dog).height-1))
			return (true);
		double dp;
		CvMat *adj = cvCreateMat (3,1,CV_32F);
		GetAdjustment (peak, peak.m_level, x, y, &dp, adj);
		// Get adjustments and check if we require further adjustments due
		// to pixel level moves. If so, turn the adjustments into real
		// changes and continue the loop. Do not adjust the plane, as we
		// are usually quite low on planes in thie space and could not do
		// further adjustments from the top/bottom planes.
		double adjS = cvGet2D (adj,0,0).val[0];//adj[0, 0];
		double adjY = cvGet2D (adj,1,0).val[0];//adj[1, 0];
		double adjX = cvGet2D (adj,2,0).val[0];//adj[2, 0];
		cvReleaseMat (&adj);
		if (abs (adjX) > 0.5 || abs (adjY) > 0.5) 
		{
			// Already adjusted the last time, give up
			if (adjusted == 0) {
				//Console.WriteLine ("too many adjustments, returning");
				return (true);
			}
			adjusted -= 1;
			// Check that just one pixel step is needed, otherwise discard
			// the point
			double distSq = adjX * adjX + adjY * adjY;
			if (distSq > 2.0)
				return (true);
			peak.m_x = (int) (peak.m_x + adjX + 0.5);
			peak.m_y = (int) (peak.m_y + adjY + 0.5);
			//point.Level = (int) (point.Level + adjS + 0.5);
			//TRACE ("moved point by (%f,%f: %f) to (%d,%d: %d)",
			//	adjX, adjY, adjS, peak.m_x, peak.m_y, peak.m_level);
			continue;
		}
		// Check if we already have a keypoint within this octave for this
		// pixel position in order to avoid dupes. (Maybe we can move this
		// check earlier after any adjustment, so we catch dupes earlier).
		// If its not in there, mark it for later searches.
		//
		// FIXME: check why there does not seem to be a dupe at all
		if (cvGet2D (processed,peak.m_y,peak.m_x).val[0] != 0)
			return (true);
		cvSet2D (processed,peak.m_y,peak.m_x,cvScalar(1));
		// Save final sub-pixel adjustments.
		peak.m_fineS = adjS;
		peak.m_fineX = adjX;
		peak.m_fineY = adjY;
		peak.m_DValue = cvGet2D (dog,peak.m_y,peak.m_x).val[0] + 0.5 * dp;
		needToAdjust = false;
	}
	return false;
}
void CScaleSpace::GetAdjustment(CScalePoint& point, int level, int x, int y, double* dp,CvMat* adj)
{
	*dp = 0.0;
	//if (point.Level <= 0 || point.Level >= (spaces.Length - 1))
	//	throw (new ArgumentException ("point.Level is not within [bottom-1;top-1] range"));
	IplImage* below = m_dogImgs[level - 1];
	IplImage* current = m_dogImgs[level];
	IplImage* above = m_dogImgs[level + 1];
	CvMat* H = cvCreateMat( 3, 3, CV_32F);
	//H[0, 0] = below[x, y] - 2 * current[x, y] + above[x, y];
	cvSet2D (H,0,0,cvScalar (cvGet2D (below,y,x).val[0] - 2 * cvGet2D(current,y,x).val[0] + cvGet2D(above,y,x).val[0]) );
	cvSet2D (H,0,1,cvScalar(0.25 * (cvGet2D(above,y+1,x).val[0] - cvGet2D(above,y-1,x).val[0] - (cvGet2D(below,y+1,x).val[0] - cvGet2D(below,y-1,x).val[0]))));
	cvSet2D (H,1,0,cvGet2D (H,0,1));
	cvSet2D (H,0,2,cvScalar(0.25 * (cvGet2D(above,y,x+1).val[0] - cvGet2D(above,y,x-1).val[0] - (cvGet2D(below,y,x+1).val[0] - cvGet2D(below,y,x-1).val[0]))));
	cvSet2D (H,2,0,cvGet2D(H,0,2)); 
	cvSet2D (H,1,1,cvScalar(cvGet2D(current,y-1,x).val[0] - 2 * cvGet2D(current,y,x).val[0] + cvGet2D(current,y+1,x).val[0]));
	cvSet2D (H,1,2,cvScalar(0.25 * (cvGet2D(current,y+1,x+1).val[0] - cvGet2D(current,y+1,x-1).val[0] -(cvGet2D(current,y-1,x+1).val[0] - cvGet2D(current,y-1,x-1).val[0]))));
	cvSet2D (H,2,1,cvGet2D(H,1,2));
	cvSet2D (H,2,2,cvScalar(cvGet2D(current,y,x-1).val[0] - 2 * cvGet2D(current,y,x).val[0] + cvGet2D(current,y,x+1).val[0]));
	CvMat* d = cvCreateMat(3, 1, CV_32F);
	cvSet2D (d,0,0,cvScalar(0.5 * (cvGet2D(above,y,x).val[0] - cvGet2D(below,y,x).val[0])));
	cvSet2D (d,1,0,cvScalar(0.5 * (cvGet2D(current,y+1,x).val[0] - cvGet2D(current,y-1,x).val[0])));
	cvSet2D (d,2,0,cvScalar(0.5 * (cvGet2D(current,y,x+1).val[0] - cvGet2D(current,y,x-1).val[0])));
	CvMat* b = cvCloneMat(d);
	//b.Negate ();
	cvSet2D (b,0,0,cvScalar(-cvGet2D(b,0,0).val[0]));
	cvSet2D (b,1,0,cvScalar(-cvGet2D(b,1,0).val[0]));
	cvSet2D (b,2,0,cvScalar(-cvGet2D(b,2,0).val[0]));
	// Solve: A x = b
	//H.SolveLinear (b);
	// adj ::= x
	cvSolve (H,b,adj);
	*dp = cvDotProduct(adj,d);
	
	cvReleaseMat (&H);
	cvReleaseMat (&d);
	cvReleaseMat (&b);
}
#define sqr(x) (x*x)
void CScaleSpace::GenMagnitudeAndDirectionMaps()
{
	// We leave the first entry to null, and ommit the last. This way, the
	// magnitudes and directions maps have the same index as the
	// imgScaled maps they belong to.
	m_magnitudes = new IplImage*[m_dogCount - 1];
	m_directions = new IplImage*[m_dogCount - 1];
	// Build the maps, omitting the border pixels, as we cannot build
	// gradient information there.
	for (int s = 1 ; s < (m_dogCount - 1); s++) 
	{
		m_magnitudes[s] = cvCreateImage (cvGetSize(m_dogImgs[0]),IPL_DEPTH_32F,1);
		m_directions[s] = cvCreateImage (cvGetSize(m_dogImgs[0]),IPL_DEPTH_32F,1);
		for (int y = 1 ; y < cvGetSize(m_gaussImgs[s]).height - 1 ; y++) 
		{
			for (int x = 1 ; x < cvGetSize(m_gaussImgs[s]).width - 1 ; x++) 
			{
				// gradient magnitude m
				double dx = cvGet2D(m_gaussImgs[s],y,x+1).val[0] - 
							cvGet2D(m_gaussImgs[s],y,x-1).val[0];
				double dy = cvGet2D(m_gaussImgs[s],y+1,x).val[0] - 
							cvGet2D(m_gaussImgs[s],y-1,x).val[0];
				cvSet2D (m_magnitudes[s],y,x,cvScalar( sqrt (dx*dx + dy*dy) ) );
				// gradient direction theta
				cvSet2D (m_directions[s],y,x,cvScalar(atan2(dy,dx) ));
			}
		}
	}	
}
void CScaleSpace::ClearMagnitudeAndDirectionMaps()
{
	for (int i = 1; i < m_dogCount-1; i++)
	{
		cvReleaseImage (&m_magnitudes[i]);
		cvReleaseImage (&m_directions[i]);
	}
	delete m_magnitudes;
	delete m_directions;
}
void CScaleSpace::GenerateKeypoints (std::vector<CScalePoint>& localizedPeaks, std::vector<CKeyPoint*>* keypoints,int scaleCount, double octaveSigma)
{
	for (int i = 0; i < localizedPeaks.size(); i++) 
	{
		std::vector<CKeyPoint*> thisKeyPoints;
		// Generate zero or more keypoints from the scale point locations.
		// TODO: make the values configurable
		// Songkran : make 36,0.8 be constant
		GenerateKeypointSingle (m_basePixScale, localizedPeaks[i], 36, 0.8, scaleCount, octaveSigma, &thisKeyPoints);
		
		// Generate the feature descriptor.
		CreateDescriptors (&thisKeyPoints,m_magnitudes[localizedPeaks[i].m_level], m_directions[localizedPeaks[i].m_level], 2.0, 4, 8, 0.2);
		
		// Only copy over those keypoints that have been successfully
		// assigned a descriptor (feature vector).
		// Songkran : all key must be already assign feature vector i think
		for (int j = 0; j < thisKeyPoints.size(); j++) 
		{
			//if (kp.HasFV == false)
			//	throw (new Exception ("should not happen"));
			// Transform the this level image relative coordinate system
			// to the original image coordinates by multiplying with the
			// current img scale (which starts with either 0.5 or 1.0 and
			// then always doubles: 2.0, 4.0, ..)
			// Note that the kp coordinates are not used for processing by
			// the detection methods and this has to be the last step.
			// Also transform the relative-to-image scale to an
			// absolute-to-original-image scale.
			// Songkran : 
			// m_kpScale = octaveSigma * exp (((point.m_level + point.m_fineS) / scaleCount) * log(2));
			CKeyPoint* kp = thisKeyPoints[j];
			if (kp->m_featureVec == NULL)
				TRACE ("Songkran : Something wrong with feature point descriptorn");
			kp->m_x *= kp->m_imgScale;
			kp->m_y *= kp->m_imgScale;
			kp->m_kpScale *= kp->m_imgScale;
			keypoints->push_back (kp);
		}
	}
}
// Average the content of the direction bins.
void CScaleSpace::AverageWeakBins (double* bins, int binCount)
{
	// TODO: make some tests what number of passes is the best. (its clear
	// one is not enough, as we may have something like
	// ( 0.4, 0.4, 0.3, 0.4, 0.4 ))
	for (int sn = 0 ; sn < 4 ; ++sn) 
	{
		double firstE = bins[0];
		double last = bins[binCount - 1];
		for (int sw = 0 ; sw < binCount ; ++sw) 
		{
			double cur = bins[sw];
			double next = (sw == (binCount - 1)) ? firstE : bins[(sw + 1) % binCount];
			bins[sw] = (last + cur + next) / 3.0;
			last = cur;
		}
	}
}
// Fit a parabol to the three points (-1.0 ; left), (0.0 ; middle) and
// (1.0 ; right).
//
// Formulas:
// f(x) = a (x - c)^2 + b
//
// c is the peak offset (where f'(x) is zero), b is the peak value.
//
// In case there is an error false is returned, otherwise a correction
// value between [-1 ; 1] is returned in 'degreeCorrection', where -1
// means the peak is located completely at the left vector, and -0.5 just
// in the middle between left and middle and > 0 to the right side. In
// 'peakValue' the maximum estimated peak value is stored.
bool CScaleSpace::InterpolateOrientation (double left, double middle,
							 double right, double* degreeCorrection, double* peakValue)
{
	// Kuas:
	// if we set f(x) = 0 then
	// a ( x^2 -2cx + c^2) + b= 0
	// x1+x2 = 2c/a ,so we get
	// a = 2c/(x1+x2)
	// and so .. i don't under stand this code
	double a = ((left + right) - 2.0 * middle) / 2.0;
	// Not a parabola
	if (a == 0.0)
		return (false);
	*peakValue = 0; //Double.NaN;
	*degreeCorrection = 0; //Double.NaN;
	double c = (((left - middle) / a) - 1.0) / 2.0;
	double b = middle - c * c * a;
	if (c < -0.5 || c > 0.5)
	{
		TRACE ("KUAS : ERROR IN InterpolateOrientationn");
	}
	*degreeCorrection = c;
	*peakValue = b;
	return (true);
}
// Create the descriptor vector for a list of keypoints.
//
// keypoints: The list of keypoints to be processed. Everything but the
//     descriptor must be filled in already.
// magnitude/direction: The precomputed gradient magnitude and direction
//     maps.
// considerScaleFactor: The downscale factor, which describes the amount
//     of pixels in the circular region relative to the keypoint scale.
//     Low values means few pixels considered, large values extend the
//     range. (Use values between 1.0 and 6.0)
// descDim: The dimension size of the output descriptor. There will be
//     descDim * descDim * directionCount elements in the feature vector.
// directionCount: The dimensionality of the low level gradient vectors.
// fvGradHicap: The feature vector gradient length hi-cap threshhold.
//     (Should be: 0.2)
//
// Some parts modelled after Alexandre Jenny's Matlab implementation.
//
// Songkran : I think this is not implements according to Lowe's paper
//			  but the idia behind is quite the same
//			  in this implementation he use variable size of window (change by scale of DoG img)
void CScaleSpace::CreateDescriptors (std::vector<CKeyPoint*>* keypoints,
	IplImage* magnitude, IplImage* direction,
	double considerScaleFactor, int descDim, int directionCount,
	double fvGradHicap)
//	2.0, 4, 8,
//	0.2);
{
	if (keypoints->size() <= 0)
		return;
	considerScaleFactor *= (*keypoints)[0]->m_kpScale;
	double dDim05 = ((double) descDim) / 2.0;
	// Now calculate the radius: We consider pixels in a square with
	// dimension 'descDim' plus 0.5 in each direction. As the feature
	// vector elements at the diagonal borders are most distant from the
	// center pixel we have scale up with sqrt(2).
	// Songkran : I think in Lowe paper he suggest to use radius = (descDim*4)/2 * sqrt(2)
	int radius = (int) (((descDim + 1.0) / 2) *
		sqrt (2.0) * considerScaleFactor + 0.5);
	// Precompute the sigma for the "center-most, border-less" gaussian
	// weighting.
	// (We are operating to dDim05, CV book tells us G(x), x > 3 sigma
	//  negligible, but this range seems much shorter!?)
	//
	// In Lowe03, page 15 it says "A Gaussian weighting function with
	// sigma equal to one half the width of the descriptor window is
	// used", so we just use his advice.
	// Songkran : should the line below recorrect to 2.0 * (r/2) * (r/2)
	double sigma2Sq = 2.0 * dDim05 * dDim05;
	for (int i = 0; i < keypoints->size(); i++)
	{
		CKeyPoint* kp = (*keypoints)[i];
		// The angle to rotate with: negate the orientation.
		double angle = -kp->m_orientation;
		kp->CreateVector (descDim, descDim, directionCount);
		//Console.WriteLine ("  FV allocated");
		for (int y = -radius ; y < radius ; ++y) 
		{
			for (int x = -radius ; x < radius ; ++x) 
			{
				// Rotate and scale
				double yR = sin (angle) * x + cos (angle) * y;
				double xR = cos (angle) * x - sin (angle) * y;
				yR /= considerScaleFactor;
				xR /= considerScaleFactor;
				// Now consider all (xR, yR) that are anchored within
				// (- descDim/2 - 0.5 ; -descDim/2 - 0.5) to
				//    (descDim/2 + 0.5 ; descDim/2 + 0.5),
				// as only those can influence the FV.
				if (yR >= (dDim05 + 0.5) || xR >= (dDim05 + 0.5) ||
					xR <= -(dDim05 + 0.5) || yR <= -(dDim05 + 0.5))
					continue;
				int currentX = (int) (x + kp->m_x + 0.5);
				int currentY = (int) (y + kp->m_y + 0.5);
				if (currentX < 1 || currentX >= (cvGetSize(magnitude).width - 1) ||
					currentY < 1 || currentY >= (cvGetSize(magnitude).height - 1))
					continue;
				// Weight the magnitude relative to the center of the
				// whole FV. We do not need a normalizing factor now, as
				// we normalize the whole FV later anyway (see below).
				// xR, yR are each in -(dDim05 + 0.5) to (dDim05 + 0.5)
				// range
				double magW = exp (-(xR * xR + yR * yR) / sigma2Sq) * 
					cvGet2D (magnitude,currentY,currentX).val[0];
				// Anchor to (-1.0, -1.0)-(dDim + 1.0, dDim + 1.0), where
				// the FV points are located at (x, y)
				// Songkran : i think this should be yR += dDim05;
				// Songkran : and range is from (-0.5,-0.5) to (dDim+0.5,dDim+0.5)
				//yR += dDim05 - 0.5;
				//xR += dDim05 - 0.5;
				yR += dDim05;
				xR += dDim05;
				// Build linear interpolation weights:
				// A B
				// C D
				//
				// The keypoint is located between A, B, C and D.
				// Songkran : this is trilinear interpolation (3 variable x,y,dir)
				int xIdx[2];
				int yIdx[2];
				int dirIdx[2];
				double xWeight[2];
				double yWeight[2];
				double dirWeight[2];
				bool flagx[2],flagy[2],flagdir[2];
				memset (flagx,false,2*sizeof(bool));
				memset (flagy,false,2*sizeof(bool));
				memset (flagdir,false,2*sizeof(bool));
				// Songkran : have to check xR < descDim also ??
				if (xR >= 0 && xR < descDim) 
				{
					xIdx[0] = (int) xR;
					xWeight[0] = (1.0 - (xR - xIdx[0]));
					flagx[0] = true;
				}
				if (yR >= 0 && yR < descDim) 
				{
					yIdx[0] = (int) yR;
					yWeight[0] = (1.0 - (yR - yIdx[0]));
					flagy[0] = true;
				}
				if (xR+1 < descDim) 
				{
					xIdx[1] = (int) (xR + 1.0);
					//	xWeight[0] + xWeight[1] = 1
					//  1 - xR + xIdx[0] + xR - xIdx[0] - 1 + 1
					xWeight[1] = xR - xIdx[1] + 1.0;
					flagx[1] = true;
				}
				if (yR+1 < descDim) 
				{
					//	yWeight[0] + yWeight[1] = 1
					//  1 - yR + yIdy[0] + yR - yIdy[0] - 1 + 1
					yIdx[1] = (int) (yR + 1.0);
					yWeight[1] = yR - yIdx[1] + 1.0;
					flagy[1] = true;
				}
				// Rotate the gradient direction by the keypoint
				// orientation, then normalize to [-pi ; pi] range.
				// Songkran : should correct to dir += 2*Math.PI ?
				double dir = cvGet2D(direction,currentY,currentX).val[0] - kp->m_orientation;
				if (dir <= -PI)
					//dir += Math.PI;
					dir += 2*PI;
				if (dir > PI)
					//dir -= Math.PI;
					dir -= 2*PI;
				double idxDir = (dir * directionCount) / (2.0 * PI);
				if (idxDir < 0.0)
					idxDir += directionCount;
				dirIdx[0] = (int) idxDir;
				dirIdx[1] = (dirIdx[0] + 1) % directionCount;
				dirWeight[0] = 1.0 - (idxDir - dirIdx[0]);
				dirWeight[1] = idxDir - dirIdx[0];
				for (int iy = 0 ; iy < 2 ; ++iy) 
					for (int ix = 0 ; ix < 2 ; ++ix) 
						for (int id = 0 ; id < 2 ; ++id) 
							if (flagx[ix] && flagy[iy])
								kp->FVSet (xIdx[ix], yIdx[iy], dirIdx[id],
								kp->FVGet (xIdx[ix], yIdx[iy], dirIdx[id]) +
								xWeight[ix] * yWeight[iy] * dirWeight[id] * magW);
			}
		}
		// Normalize and hicap the feature vector, as recommended on page
		// 16 in Lowe03.
		kp->CapAndNormalizeFV (fvGradHicap);
	}
}
// Assign each feature point one or more standardized orientations.
// (section 5 in Lowe's paper)
//
// We use an orientation histogram with 36 bins, 10 degrees each. For
// this, every pixel (x,y) lieing in a circle of 'squareDim' diameter
// within in a 'squareDim' sized field within the image L ('gaussImg') is
// examined and two measures calculated:
//
//    m = sqrt{ (L_{x+1,y} - L_{x-1,y})^2 + (L_{x,y+1} - L_{x,y-1})^2 }
//    theta = tan^{-1} ( frac{ L_{x,y+1} - L_{x,y-1} }
//                { L_{x+1,y} - L_{x-1,y} } )
//
// Where m is the gradient magnitude around the pixel and theta is the
// gradient orientation. The 'imgScale' value is the octave scale,
// starting with 1.0 at the finest-detail octave, and doubling every
// octave. The gradient orientations are discreetized to 'binCount'
// directions (should be: 36). For every peak orientation that lies within
// 'peakRelThresh' of the maximum peak value, a keypoint location is
// added (should be: 0.8).
//
// Note that 'space' is the gaussian smoothed original image, not the
// difference-of-gaussian one used for peak-search.
// Kuas : imgScale is basePixScale (start scale (relative dimension with input img) of that octave)
void CScaleSpace::GenerateKeypointSingle (double imgScale, CScalePoint point,
	int binCount, double peakRelThresh, int scaleCount,
	double octaveSigma, std::vector<CKeyPoint*>*vec)
{
	// The relative estimated keypoint scale. The actual absolute keypoint
	// scale to the original image is yielded by the product of imgScale.
	// But as we operate in the current octave, the size relative to the
	// anchoring images is missing the imgScale factor.
	//double kpScale = octaveSigma * exp (((point.m_level + point.m_fineS) / scaleCount) * log(2));
double kpScale = octaveSigma * (exp (((point.m_level) / scaleCount) * log(2))+ point.m_fineS);
	// Lowe03, "A gaussian-weighted circular window with a sigma three
	// times that of the scale of the keypoint".
	//
	// With sigma = 3.0 * kpScale, the square dimension we have to
	// consider is (3 * sigma) (until the weight becomes very small).
	// Kuas : Lowe04 say sigma = 1.5 * kpScale
	double sigma = 1.5 * kpScale;
	int radius = (int) (3.0 * sigma / 2.0 + 0.5);
	int radiusSq = radius * radius;
	IplImage* magnitude = m_magnitudes[point.m_level];
	IplImage* direction = m_directions[point.m_level];
	// As the point may lie near the border, build the rectangle
	// coordinates we can still reach, minus the border pixels, for which
	// we do not have gradient information available.
	int xMin = point.m_x - radius > 1 ? point.m_x - radius : 1;
	int xMax = point.m_x + radius < cvGetSize(magnitude).width-1 ? point.m_x + radius : cvGetSize(magnitude).width-1;
	int yMin = point.m_y - radius > 1 ? point.m_y - radius : 1;
	int yMax = point.m_y + radius < cvGetSize(magnitude).height-1 ? point.m_y + radius : cvGetSize(magnitude).height-1;
	// Precompute 1D gaussian divisor (2 sigma^2) in:
	// G(r) = e^{-frac{r^2}{2 sigma^2}}
	double gaussianSigmaFactor = 2.0 * sigma * sigma;
	double *bins = new double[binCount];
	memset (bins,0,sizeof(double)*binCount);
	// Build the direction histogram
	for (int y = yMin ; y < yMax ; ++y) {
		for (int x = xMin ; x < xMax ; ++x) {
			// Only consider pixels in the circle, else we might skew the
			// orientation histogram by considering more pixels into the
			// corner directions
			int relX = x - point.m_x;
			int relY = y - point.m_y;
			if (sqr(relX) + sqr(relY) > radiusSq)
				continue;
			// The gaussian weight factor.
			double gaussianWeight = exp (- ((relX * relX + relY * relY) / gaussianSigmaFactor));
			// find the closest bin and add the direction
			//int binIdx = FindClosestRotationBin (binCount, direction[x, y]);
			double angle = cvGet2D (direction,y,x).val[0];
			angle += PI;
			angle /= 2.0 * PI;
			// calculate the aligned bin
			angle *= binCount;
			int binIdx = (int) angle;
			if (binIdx == binCount)
				binIdx = 0;
			bins[binIdx] += cvGet2D(magnitude,y,x).val[0] * gaussianWeight;
		}
	}
	// As there may be succeeding histogram entries like this:
	// ( ..., 0.4, 0.3, 0.4, ... ) where the real peak is located at the
	// middle of this three entries, we can improve the distinctiveness of
	// the bins by applying an averaging pass.
	//
	// TODO: is this really the best method? (we also loose a bit of
	// information. Maybe there is a one-step method that conserves more)
	AverageWeakBins (bins, binCount);
	// find the maximum peak in gradient orientation
	double maxGrad = 0.0;
	int maxBin = 0;
	for (int b = 0 ; b < binCount ; ++b) {
		if (bins[b] > maxGrad) {
			maxGrad = bins[b];
			maxBin = b;
		}
	}
	// First determine the real interpolated peak high at the maximum bin
	// position, which is guaranteed to be an absolute peak.
	//
	// XXX: should we use the estimated peak value as reference for the
	//   0.8 check or the original bin-value?
	double maxPeakValue, maxDegreeCorrection;
	InterpolateOrientation (bins[maxBin == 0 ? (binCount - 1) : (maxBin - 1)],
		bins[maxBin], bins[(maxBin + 1) % binCount],
		&maxDegreeCorrection, &maxPeakValue);
	// Now that we know the maximum peak value, we can find other keypoint
	// orientations, which have to fulfill two criterias:
	//
	//  1. They must be a local peak themselves. Else we might add a very
	//     similar keypoint orientation twice (imagine for example the
	//     values: 0.4 1.0 0.8, if 1.0 is maximum peak, 0.8 is still added
	//     with the default threshhold, but the maximum peak orientation
	//     was already added).
	//  2. They must have at least peakRelThresh times the maximum peak
	//     value.
	bool* binIsKeypoint = new bool[binCount];
	for (b = 0 ; b < binCount ; ++b) 
	{
		binIsKeypoint[b] = false;
		// The maximum peak of course is
		if (b == maxBin) 
		{
			binIsKeypoint[b] = true;
			continue;
		}
		// Local peaks are, too, in case they fulfill the threshhold
		if (bins[b] < (peakRelThresh * maxPeakValue))
			continue;
		int leftI = (b == 0) ? (binCount - 1) : (b - 1);
		int rightI = (b + 1) % binCount;
		if (bins[b] <= bins[leftI] || bins[b] <= bins[rightI])
			continue;	// not local peak
		binIsKeypoint[b] = true;
	}
	// All the valid keypoint bins are now marked in binIsKeypoint, now
	// build them.
	// find other possible locations
	double oneBinRad = (2.0 * PI) / binCount;
	for (b = 0 ; b < binCount ; ++b) 
	{
		if (binIsKeypoint[b] == false)
			continue;
		int bLeft = (b == 0) ? (binCount - 1) : (b - 1);
		int bRight = (b + 1) % binCount;
		// Get an interpolated peak direction and value guess.
		double peakValue;
		double degreeCorrection;
		if (InterpolateOrientation (bins[bLeft], bins[b], bins[bRight],
			&degreeCorrection, &peakValue) == false)
		{
			//throw (new InvalidOperationException ("BUG: Parabola fitting broken"));
			TRACE ("BUG: Parabola fitting brokenn");
		}
		// [-1.0 ; 1.0] -> [0 ; binrange], and add the fixed absolute bin
		// position.
		// We subtract PI because bin 0 refers to 0, binCount-1 bin refers
		// to a bin just below 2PI, so -> [-PI ; PI]. Note that at this
		// point we determine the canonical descriptor anchor angle. It
		// does not matter where we set it relative to the peak degree,
		// but it has to be constant. Also, if the output of this
		// implementation is to be matched with other implementations it
		// must be the same constant angle (here: -PI).
		double degree = (b + degreeCorrection) * oneBinRad - PI;
		if (degree < -PI)
			degree += 2.0 * PI;
		else if (degree > PI)
			degree -= 2.0 * PI;
		CKeyPoint* kp = new CKeyPoint (m_gaussImgs[point.m_level],
			point.m_x + point.m_fineX,
			point.m_y + point.m_fineY,
			imgScale, kpScale, degree);
		vec->push_back (kp);
	}
	delete binIsKeypoint;
	delete bins;
}

						