基于sift的特征点匹配C语言源码

最新推荐文章于 2022-11-18 13:50:28 发布
最新推荐文章于 2022-11-18 13:50:28 发布
阅读量740 收藏 3
点赞数 2
CC 4.0 BY-SA版权
文章标签： sift c语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/yhtao923/article/details/99053726
本文分享了一段基于SIFT算法的特征点匹配源码，深入解析了从图像预处理到特征点检测与描述符生成的全过程。代码中详细展示了如何利用SIFT算法在不同尺度空间内检测关键点，并计算其描述符，最后通过比较描述符实现特征点匹配。适合计算机视觉和图像处理领域的研究者及开发者阅读。

摘要生成于 C知道 ，由 DeepSeek-R1 满血版支持， 前往体验 >

基于sift的特征点匹配源码，其中头文件在上传的附件中，本代码多年前本人所写，不足之处请批评指正。其中调用的部分头文件以及源码在之前上传的资源中。

#include <math.h>
#include <iostream> 
#include <string.h>
#include "matrix.h"
#include <time.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include "imageBase.h"
#include "image.h"
#include "keypoint.h"
void charImg2FloatImg(IMAGE_t *srcImg,IMAGE_t *dstImg)
{
	unsigned char *srcImgFData0;
	float *dstImgData0;
	int i,j;
	
	srcImgFData0 = (unsigned char *)srcImg->imageData;
	dstImgData0 = (float *)dstImg->imageData;
	for(i = 0;i < srcImg->height;i++)
	{
		for(j = 0;j < srcImg->width;j++)
		{
			dstImgData0[i * dstImg->width + j] = (float)((unsigned char)srcImgFData0[i * srcImg->widthStep + j]);

		}
	}
}

int getHaarKeyPoint(IMAGE_t *srcImg0,vect2f **srcPoint,int Yshift,int startPad)
{
	vect2f *pPoint0,tempPoint;
	IMAGE_t *edgeImg = NULL;
	IMAGE_t *srcFImg0,*edgeFImg;
	int i,j,p,q;
	float Dxy,Dxx,Dyy,Dyx,trace,Det;
	int x,y;
	int step = 8,widthStep;
	float *srcData;
	int pointNum = 0;
	float max = 0.0f;
	srcFImg0 = (IMAGE_t *)mvCreateImage(srcImg0->width,srcImg0->height,IPL_DEPTH_8U,1);
	edgeFImg = (IMAGE_t *)mvCreateImage(srcImg0->width,srcImg0->height,IPL_DEPTH_32F,1);

	charImg2FloatImg(srcImg0,edgeImg);
	widthStep = edgeFImg->width;
	srcData = (float *)edgeFImg->imageData;
	
	pPoint0 = (vect2f *)malloc(sizeof(vect2f) * 2000);
	for(i = 16;i < edgeImg->height - 16;i ++)
	{
		for(j = startPad + 16;j < edgeImg->width - 16;j ++)
		{
			max = 0.0f;
			int distOkFlag = 0;
			float distVal;
			for(p = 0;p < pointNum;p++)
			{
				tempPoint = pPoint0[p];
				distVal = (i - tempPoint.y) * (i - tempPoint.y) + (j - tempPoint.x) * (j - tempPoint.x);
				if(distVal > (512))
				{
					distOkFlag++;
				}
			}
			if(distOkFlag < pointNum)
				continue;

			for(p = -8;p < 8;p++)
			{
				for(q = -8;q < 8;q++)
				{
					x = j + q;
					y = i + p;
					Dxx = fabs((srcData[y * widthStep + x + step] - srcData[y * widthStep + x - 0])
						 - (srcData[y * widthStep + x + 0] + srcData[y * widthStep + x - step]));
					Dxy = fabs((srcData[(y - step) * widthStep + x + step] - srcData[y * widthStep + x - 0])
						 - (srcData[y * widthStep + x + 0] + srcData[(y + step) * widthStep + x - step]));
					Dyy = fabs((srcData[(y + step) * widthStep + x] - srcData[y * widthStep + x - 0])
						 - (srcData[y * widthStep + x + 0] + srcData[(y - step) * widthStep + x - 0]));
					Dyx = fabs((srcData[(y + step) * widthStep + x + step] - srcData[y * widthStep + x - 0])
						 - (srcData[y * widthStep + x + 0] + srcData[(y - step) * widthStep + x - step]));
					trace 	= Dxx + Dyy;
					Det 	= Dxx * Dyy - Dxy * Dyx;
					Det = fabs((trace * trace) - Det);
					if(Det > max)
					{
						max = Det;
						pPoint0[pointNum].x = (float)x;
						pPoint0[pointNum].y = (float)y;
					}
				}
			}

			
	
			//printf("Det: %f %f %f %f %f %d %d\n",Det,Dxx,Dyy,Dxy,Dyx,x,y);
			if(max > 5000.0f)
			{
				pointNum++;
			}
		}

	}
	
	*srcPoint = pPoint0;
	return pointNum;
}
void inter_hist_entry( float* hist, float rbin, double cbin,  
                               float obin, float mag, int d, int n )  
{  
    float d_r, d_c, d_o, v_r, v_c, v_o;  
    float * row, * h;  
    int r0, c0, o0, rb, cb, ob, r, c, o;  
  
    r0 = mvFloor( rbin );  
    c0 = mvFloor( cbin );  
    o0 = mvFloor( obin );  
    d_r = rbin - r0;  
    d_c = (float)cbin - c0;  
    d_o = obin - o0;  
  
    /* 
    The entry is distributed into up to 8 bins.  Each entry into a bin 
    is multiplied by a weight of 1 - d for each dimension, where d is the 
    distance from the center value of the bin measured in bin units. 
    */  
    for( r = 0; r <= 1; r++ )  
    {  
        rb = r0 + r;  
        if( rb >= 0  &&  rb < d )  
        {  
            v_r = mag * ( ( r == 0 )? 1.0f - d_r : d_r );  
            row = &hist[rb * d * n];  
            for( c = 0; c <= 1; c++ )  
            {  
                cb = c0 + c;  
                if( cb >= 0  &&  cb < d )  
                {  
                    v_c = v_r * ( ( c == 0 )? 1.0f - d_c : d_c );  
                    h = &row[cb * n];  
                    for( o = 0; o <= 1; o++ )  
                    {  
                        ob = ( o0 + o ) % n;  
                        v_o = v_c * ( ( o == 0 )? 1.0f - d_o : d_o );  
                        h[ob] += v_o;  
                    }  
                }  
            }  
        }  
    }  
}  

void imgDeviaOfGauss(IMAGE_t *Img0,IMAGE_t *Img1,IMAGE_t *ImgDst)
{
	int p,q,c;
	int img0Width;
	int ch = Img1->nChannels;
	float *imgData0,*imgData1,*imgData2 = (float *)malloc(Img0->width * Img0->height * sizeof(float));
	
	imgData0 = (float *)Img0->imageData;
	imgData1 = (float *)Img1->imageData;
	img0Width = Img1->width;
	
	for(p = 0;p < Img1->height;p++)
	{
		for(q = 0;q < Img1->width;q++)
		{
			uint idx = (p * img0Width + q) * ch;
			for(c = 0;c < ch;c++)
			{
				imgData2[idx + c] = (imgData1[idx + c] - imgData0[idx + c]);
				//printf("%f\n",imgData2[idx + c] );
			}
		}
	}
	if(ImgDst->depth == IPL_DEPTH_8U)
		mvImagef32Norm(imgData2,ImgDst->imageData,Img1->width,Img1->height,ch,0);
	else if(ImgDst->depth == IPL_DEPTH_32F)
		memcpy((uchar *)ImgDst->imageData,imgData2,Img0->width * Img0->height * sizeof(float));

	free(imgData2);
	
}

int getSiftKeyPoint(IMAGE_t *srcImg,SIFT_KEYPOINTS_FEATER_t * keyPointFeater,
	int gaussSize,int scaleSpace,float sigma,int layer)
{
	
	int i,j,p,q;
	int layerNum = 0;
	//IplImage *cvImage0,*cvImage1;
//	CvPoint point;	
	float **kernel32fPtr;
	//cvShowImage("leftTempImg",leftTempImg);

	
	if(layer < 0 || layer > 8 || gaussSize > 7 || srcImg == NULL)
		return NULL;
	if(gaussSize != 3 && gaussSize != 5 && gaussSize != 7)
		return NULL;
	SIFT_KEYPOINTS_t *keyPointInf;

	int pointIndex = 0,keyPointIndex = 0,keyPointIndex1 = 0,maxPointNum = ((srcImg->width * srcImg->height) >> 4);
	IMAGE_t *_srcImg;
	IMAGE_t **gaussImage = (IMAGE_t **)malloc(scaleSpace * layer * sizeof(IMAGE_t *));
	IMAGE_t **dogImage = (IMAGE_t **)malloc(scaleSpace * layer * sizeof(IMAGE_t *));
	kernel32fPtr = (float **)malloc(scaleSpace * layer * sizeof(float *));
	keyPointInf = (SIFT_KEYPOINTS_t *)malloc(maxPointNum * sizeof(SIFT_KEYPOINTS_t));
	


	if(srcImg->depth != IPL_DEPTH_8U)
	{
		return NULL;
	}
	IMAGE_t *srcImgeCh1;
	if(srcImg->nChannels == 3)
	{
		srcImgeCh1 = mvCreateImage(srcImg->width,srcImg->height,IPL_DEPTH_8U,1);
		mvImageBGR2Gray(srcImg,srcImgeCh1);
		
	}
	else
	{
		srcImgeCh1 = srcImg;
	}

	_srcImg = mvCreateImage(srcImgeCh1->width,srcImgeCh1->height,IPL_DEPTH_32F,1);
	float *srcImageData = (float *)_srcImg->imageData;
	for(i = 0;i < srcImgeCh1->width * srcImgeCh1->height;i++)
	{
		srcImageData[i] = (float)srcImgeCh1->imageData[i];
	}

	for(i = 0;i < layer;i++)
	{
		for(j = 0;j < scaleSpace;j++)
		{
			float *kernel32f = (float *)malloc((2 * (j + 1) + 1) * sizeof(float));
			initGaussKernel(NULL,kernel32f,sigma * pow(2.0f,(float)((float)i + (float)j/scaleSpace)),(2 * (j + 1) + 1));
			gaussImage[i * scaleSpace + j] = getGaussDownPyramidsf32(_srcImg,kernel32f,(2 * (j + 1) + 1),i);
			kernel32fPtr[i * scaleSpace + j] = (float *)kernel32f;
			//cvImage0 = mvImage2CvImage(gaussImage[i * scaleSpace + j]);
		
			//cvShowImage("image0",cvImage0);
			//cvWaitKey(0);
		}
	}
	for(i = 0;i < layer;i++)
	{
		for(j = 0;j < scaleSpace - 1;j++)
		{
			int index = i * scaleSpace + j;
			dogImage[i * (scaleSpace - 1) + j] 
				= mvCreateImage(_srcImg->width >> i,_srcImg->height >> i,IPL_DEPTH_8U,1);
			imgDeviaOfGauss(gaussImage[index],gaussImage[index + 1],dogImage[i * (scaleSpace - 1) + j]);
			//cvImage0 = mvImage2CvImage(dogImage[i * (scaleSpace - 1) + j]);
		
			//cvShowImage("image0",cvImage0);
			//cvWaitKey(0);
		}
	}
	uchar **pDogImgData = (uchar **)malloc(layer * (scaleSpace - 1) * sizeof(uchar *));
	for(layerNum = 0;layerNum < layer;layerNum++)
	{
		for(p = 0;p < scaleSpace - 1;p++)
		{
			pDogImgData[layerNum * (scaleSpace - 1) + p] = dogImage[layerNum * (scaleSpace - 1) + p]->imageData;
		}
	}
	uchar calcMaxPix[27];
	int height,width;
	uchar maxFlag = 0;
	for(layerNum = 0;layerNum < layer;layerNum++)
	{
		height = dogImage[layerNum * (scaleSpace - 1)]->height;
		width = dogImage[layerNum * (scaleSpace - 1)]->width;
		int winSize = 32 >> layerNum;
		for(i = 4;i < height - 4;i++)
		{
			for(j = 4;j < width - 4;j++)
			{
				for(p = 1;p < scaleSpace - 2;p++)
				{
					uchar *pCalcMaxPix = calcMaxPix;
					uchar *currDogImgData = pDogImgData[layerNum * (scaleSpace - 1) + p - 0];
					uchar *lastDogImgData = pDogImgData[layerNum * (scaleSpace - 1) + p - 1];
					uchar *nextDogImgData = pDogImgData[layerNum * (scaleSpace - 1) + p + 1];
					uchar currPix =  currDogImgData[(i - 0) * width + j - 0];
					*pCalcMaxPix++ = currDogImgData[(i - 0) * width + j - 0];
					
					*pCalcMaxPix++ = lastDogImgData[(i - 1) * width + j - 1];
					*pCalcMaxPix++ = lastDogImgData[(i - 1) * width + j - 0];
					*pCalcMaxPix++ = lastDogImgData[(i - 1) * width + j + 1];
					*pCalcMaxPix++ = lastDogImgData[(i - 0) * width + j - 1];
					*pCalcMaxPix++ = lastDogImgData[(i - 0) * width + j - 0];
					*pCalcMaxPix++ = lastDogImgData[(i - 0) * width + j + 1];
					*pCalcMaxPix++ = lastDogImgData[(i + 1) * width + j - 1];
					*pCalcMaxPix++ = lastDogImgData[(i + 1) * width + j - 0];
					*pCalcMaxPix++ = lastDogImgData[(i + 1) * width + j + 1];

					*pCalcMaxPix++ = currDogImgData[(i - 1) * width + j - 1];
					*pCalcMaxPix++ = currDogImgData[(i - 1) * width + j - 0];
					*pCalcMaxPix++ = currDogImgData[(i - 1) * width + j + 1];
					*pCalcMaxPix++ = currDogImgData[(i - 0) * width + j - 1];
					*pCalcMaxPix++ = currDogImgData[(i - 0) * width + j + 1];
					*pCalcMaxPix++ = currDogImgData[(i + 1) * width + j - 1];
					*pCalcMaxPix++ = currDogImgData[(i + 1) * width + j - 0];
					*pCalcMaxPix++ = currDogImgData[(i + 1) * width + j + 1];

					*pCalcMaxPix++ = nextDogImgData[(i - 1) * width + j - 1];
					*pCalcMaxPix++ = nextDogImgData[(i - 1) * width + j - 0];
					*pCalcMaxPix++ = nextDogImgData[(i - 1) * width + j + 1];
					*pCalcMaxPix++ = nextDogImgData[(i - 0) * width + j - 1];
					*pCalcMaxPix++ = nextDogImgData[(i - 0) * width + j - 0];
					*pCalcMaxPix++ = nextDogImgData[(i - 0) * width + j + 1];
					*pCalcMaxPix++ = nextDogImgData[(i + 1) * width + j - 1];
					*pCalcMaxPix++ = nextDogImgData[(i + 1) * width + j - 0];
					*pCalcMaxPix++ = nextDogImgData[(i + 1) * width + j + 1];
					maxFlag = 0;
					for(q = 1;q < 27;q++)
					{
						if(currPix < calcMaxPix[q])
						{
							maxFlag = 1;
							q = 26;
						}
					}
					if(maxFlag == 0)
					{
						if(pointIndex < maxPointNum)
						{
							/*
							int k,g;
							for(k = 0;k < 3;k++)
							{
								for(g = 0;g < 3;g++)
								{
									printf("%d ",calcMaxPix[k * 3 + g]);
								}
								printf("\n");
							}
							printf("\n");
							*/
							keyPointInf[pointIndex].dogData = currPix;
							keyPointInf[pointIndex].x = (float)(j);
							keyPointInf[pointIndex].y = (float)(i);
							keyPointInf[pointIndex].scl_octv = sigma * pow(2.0f,(float)((float)layerNum + (float)p/scaleSpace));
							keyPointInf[pointIndex].scl_octv_num = p;
							keyPointInf[pointIndex].layer = layerNum;
							keyPointInf[pointIndex].valid = 1;
							
						}
						pointIndex++;
	
					}
					
				}	
				
			}
		}
	}
	float belt,coff,amod;
	float dyy,dxx,dxy,dyx;
	float deltx = 0.0f,delty = 0.0f,delts = 0.0f;
	int interp;
	float mag,ang;
	int bin,radius;
	float histSigma;
	float exp_denom;
	float dx,dy,ds,d2x,d2y,d2s,dxs,dys;
	int imageIndex,dogIndex;
	float hist[MV_SIFT_ORI_HIST_BINS];
	double prev, tmp, h0;
	double omax;  
	int maxbin; 
	float  w;
	float hist_width; 
	int n = 8;
	int d = 4;
	float sum = 0.0f;
	float scale;
	int pad;
	float scl_octv_num = 0.0f;
	for(i = 0;i < pointIndex;i++)
	{
		maxFlag = 0;
		imageIndex = keyPointInf[i].layer * (scaleSpace) + (keyPointInf[i].scl_octv_num);
		dogIndex = keyPointInf[i].layer * (scaleSpace - 1) + (keyPointInf[i].scl_octv_num);
		height = dogImage[dogIndex]->height;
		width = dogImage[dogIndex]->width;
		uchar *pTmpImageData = (uchar *)pDogImgData[dogIndex]+ mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);
		dyy = (float)(pTmpImageData[width] - pTmpImageData[-width]);
		dxx = (float)(pTmpImageData[1] - pTmpImageData[-1]);
		dxy = (float)(pTmpImageData[width - 1] - pTmpImageData[-width + 1]);
		dyx = (float)(pTmpImageData[width + 1] - pTmpImageData[-width - 1]);
		belt = (dxx * dyy - dxy * dyx);
		coff = 100.0f;
		if(fabs(belt) > 0.001f)
		{
			coff = (dxx + dyy) * (dxx + dyy)/belt;
		}
		if(fabs(coff) > 12.1f)
		{
			keyPointInf[i].valid = 0;
			continue;
		}
		hist_width = MV_SIFT_DESCR_SCL_FCTR * keyPointInf[i].scl_octv;
		radius = (int)(hist_width *  ( (float)d + 1.0f ) + 0.5f);
		//radius = MV_SIFT_PAD_SIZE * 2;
		if(keyPointInf[i].y < radius || keyPointInf[i].x < radius
			|| keyPointInf[i].y > (height - radius) || keyPointInf[i].x > (width - radius))
			continue;
		
		uchar *pTmpImageDataPrev = (uchar *)pDogImgData[dogIndex - 1] + mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);
		uchar *pTmpImageDataNext = (uchar *)pDogImgData[dogIndex + 1] + mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);

		for(interp = 0;interp < 5;interp++)
		{
			dx = dxx/2.0f;
			dy = dyy/2.0f;
			ds = (pTmpImageData[0] - pTmpImageDataPrev[0])/2.0f;

			d2x = (float)(pTmpImageData[1] + pTmpImageData[-1]) - 2 * pTmpImageData[0];
			d2y = (float)(pTmpImageData[width] + pTmpImageData[-width]) - 2 * pTmpImageData[0];
			d2s = (float)(pTmpImageDataPrev[0] + pTmpImageDataNext[0]) - 2 * pTmpImageData[0];
			dxy = (float)((pTmpImageData[-width - 1] + pTmpImageData[width + 1]) 
				- (pTmpImageData[-width + 1] + pTmpImageData[width - 1]))/4;
			dxs = (float)((pTmpImageDataNext[-1] + pTmpImageDataPrev[1])
				- (pTmpImageDataPrev[-1] + pTmpImageDataNext[1]))/4;
			dys = (float)((pTmpImageDataNext[width] - pTmpImageDataPrev[-width])
				- (pTmpImageDataNext[-width] + pTmpImageDataPrev[width]))/4;
	
			float a[3][3] = {{d2x,	dxy,	dxs},
							{dxy,	d2y,	dys},
							{dxs,	dys,	d2s}};
			float ia[3][3];
			amod = a[0][0] * a[1][1] * a[2][2] + a[0][1] * a[1][2] * a[2][0] + a[0][2] * a[1][0] * a[2][1]
				 - a[0][0] * a[1][2] * a[2][1] - a[0][1] * a[1][0] * a[2][2] - a[0][2] * a[1][1] * a[2][0];
			if(fabs(amod) > 0.001f)
			{
				ia[0][0] = (a[1][1] * a[2][2] - a[2][1] * a[1][2]);
				ia[0][1] = -(a[0][1] * a[2][2] - a[2][1] * a[0][2]);
				ia[0][2] = a[0][1] * a[1][2] - a[0][2] * a[1][1];
				ia[1][0] = a[1][2] * a[2][0] - a[2][2] * a[1][0];
				ia[1][1] = -a[0][2] * a[2][0] - a[2][2] * a[0][0];
				ia[1][2] = a[0][2] * a[1][0] - a[0][0] * a[1][2];
				ia[2][0] = a[1][0] * a[2][1] - a[2][0] * a[1][1];
				ia[2][1] = -a[0][0] * a[2][1] - a[2][0] * a[0][1];
				ia[2][2] = a[0][0] * a[1][1] - a[0][1] * a[1][0];
				for(p = 0;p < 3;p++)
				{
					for(q = 0;q < 3;q++)
					{
						ia[p][q] = ia[0][0]/amod;
					}
				}
				deltx = ia[0][0] * dx + ia[0][1] * dy + ia[0][2] * ds;
				delty = ia[1][0] * dx + ia[1][1] * dy + ia[1][2] * ds;
				delts = ia[2][0] * dx + ia[2][1] * dy + ia[2][2] * ds;
			}
			else
			{
				maxFlag = 1;
				break;
			}
			if(fabs(deltx) < 0.05f && fabs(delty) < 0.05f && fabs(delts) < 0.05f)
			{
				break;
			}
			if(fabs(deltx) > 1.00f || fabs(delty) > 1.00f || fabs(delts) > 1.00f)
			{
				maxFlag = 1;
				break;
			}
			keyPointInf[i].x += deltx;
			keyPointInf[i].y += delty;
			scl_octv_num += delts;
			keyPointInf[i].scl_octv_num = mvRound(scl_octv_num);
			if(keyPointInf[i].x < 8
				|| keyPointInf[i].y < 8
				|| keyPointInf[i].x > (width - 8)
				|| keyPointInf[i].y > (height - 8))
			{
				maxFlag = 1;
				break;
			}
			pTmpImageData = (uchar *)pDogImgData[dogIndex] + mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);
			pTmpImageDataPrev = (uchar *)pDogImgData[dogIndex - 1] + mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);
			pTmpImageDataNext = (uchar *)pDogImgData[dogIndex + 1] + mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);
			dx = (pTmpImageData[1] - pTmpImageData[-1])/2.0f;
			dy = (pTmpImageData[width] - pTmpImageData[-width])/2.0f;
			ds = (pTmpImageData[0] - pTmpImageDataPrev[0])/2.0f;
			float contr = dx * keyPointInf[i].x + 
				dy * keyPointInf[i].x + ds * keyPointInf[i].scl_octv;
			if(fabs(contr) < (keyPointInf[i].scl_octv/scaleSpace))
			{
				maxFlag = 1;
				break;
			}
		}
		if(maxFlag)
		{
			keyPointInf[i].valid = 0;
			continue;
		}
		memset(hist,0,sizeof(hist));
		histSigma = keyPointInf[i].scl_octv * MV_SIFT_ORI_SIG_FCTR;
		exp_denom = 2.0f * histSigma * histSigma;
		radius = mvRound(MV_SIFT_ORI_RADIUS * keyPointInf[i].scl_octv);
		imageIndex = keyPointInf[i].layer * scaleSpace + mvRound(keyPointInf[i].scl_octv_num);
		//float *pF32TmpImageData = (float *)gaussImage[imageIndex]->imageData + mvRound(keyPointInf[i].y) * width + mvRound(keyPointInf[i].x);
		//int pad = radius + MV_SIFT_PAD_SIZE - mvRound(keyPointInf[i].x);
		//if(pad < 0)
			pad = 0;
		for(p = -radius + pad;p < radius + pad;p++)
		{
			for(q = -radius + pad;q < radius + pad;q++)
			{
				dyy = (float)(pTmpImageData[width * (p - 1) + q] - pTmpImageData[width * (p + 1) + q]);
				dxx = (float)(pTmpImageData[width * p + q + 1] - pTmpImageData[width * p + q - 1]);
				mag = sqrt((float)(dxx * dxx + dyy * dyy));
				ang = atan2((float)dyy,(float)dxx);
				w = exp( -( p*p + q*q ) / exp_denom );  
                bin = mvRound( MV_SIFT_ORI_RADIUS * ( ang + PI ) / (2.0f * PI) );  
                bin = ( bin < MV_SIFT_ORI_HIST_BINS )? bin : 0;  
                hist[bin] += w * mag;  
				
			}
		}
		
		h0 = hist[0];
		prev = hist[MV_SIFT_ORI_HIST_BINS - 1];  
	    for( p = 0; p < MV_SIFT_ORI_HIST_BINS; p++ )  
	    {  
	        tmp = hist[p];  
	        hist[p] = 0.25f * (float)prev + 0.5f * hist[p] +   
	            0.25f * ( ( p+1 == MV_SIFT_ORI_HIST_BINS )? (float)h0 : hist[p+1] );  
	        prev = tmp;  
	    }  
	 

		omax = hist[0];  
		maxbin = 0;  
		for( p = 1; p < MV_SIFT_ORI_HIST_BINS; p++ )  
		{
			if( hist[p] > omax )  
			{  
			    omax = hist[p];  
			    maxbin = p;  
			}  
		}
		keyPointFeater[keyPointIndex].ang = (float)maxbin * MV_SIFT_ORI_HIST_BINS;
		keyPointFeater[keyPointIndex].scl_octv = keyPointInf[i].scl_octv;
		keyPointFeater[keyPointIndex].x = keyPointInf[i].x;
		keyPointFeater[keyPointIndex].y = keyPointInf[i].y;
		keyPointFeater[keyPointIndex].layer = keyPointInf[i].layer;
		keyPointFeater[keyPointIndex].scl_octv_num = keyPointInf[i].scl_octv_num;
		for( p = 1; p < MV_SIFT_ORI_HIST_BINS; p++ )
		{
			if( hist[p] > hist[p - 1] && hist[p] > hist[p + 1] && hist[p] > (maxbin * MV_SIFT_ORI_PEAK_RATIO))
			{
				keyPointFeater[keyPointIndex].ang = (float)p * MV_SIFT_ORI_HIST_BINS;
				keyPointFeater[keyPointIndex].scl_octv = keyPointInf[i].scl_octv;
				keyPointFeater[keyPointIndex].x = keyPointInf[i].x;
				keyPointFeater[keyPointIndex].y = keyPointInf[i].y;
				keyPointFeater[keyPointIndex].layer = keyPointInf[i].layer;
				keyPointFeater[keyPointIndex].scl_octv_num = keyPointInf[i].scl_octv_num;
				keyPointIndex++;
				break;
			}
		}
		
		if(imageIndex == 5)
		{
			//point.x = mvRound(keyPointFeater[keyPointIndex].x);
			//point.y = mvRound(keyPointFeater[keyPointIndex].y);
			//cvCircle(cvImage0,point,4,cvScalar(255u,255u,255u));
		}
		
	}
	
	double r_rot, c_rot, grad_ori = 0.0f, rbin, cbin, obin, PI2 = 2.0 * PI;  
	
	int count = 0;
	int row,col;
	float cos_t;  
    float sin_t;
	float bins_per_rad;  
    float r_dot,c_dot;
    float pixVal,x00,x01,x10,x11;
	float descF32[128];

//	cvImage0 = mvImage2CvImage(gaussImage[scaleSpace]);
	for(p = 0;p < keyPointIndex;p++)
	{
		memset(descF32,0,sizeof(descF32));
		imageIndex = keyPointFeater[p].layer * scaleSpace + mvRound(keyPointFeater[p].scl_octv_num);
		height = gaussImage[imageIndex]->height;
		width = gaussImage[imageIndex]->width;
		float *pTmpImageData = (float *)gaussImage[imageIndex]->imageData + mvRound(keyPointFeater[p].y) * width + mvRound(keyPointFeater[p].x);

		cos_t = cos( keyPointFeater[p].ang );
		sin_t = sin( keyPointFeater[p].ang );
		bins_per_rad = (float)n / (float)PI2;
		exp_denom = d * d * 0.5f;
		hist_width = MV_SIFT_DESCR_SCL_FCTR * keyPointFeater[p].scl_octv;
		radius = (int)((float)hist_width * sqrt(2.0f) * ( d + 1.0f ) * 0.5f + 0.5f);
		float *pImgRng = (float *)malloc(4 * radius * radius * sizeof(float));
		// = radius + MV_SIFT_PAD_SIZE - mvRound(keyPointInf[i].x);
		//if(pad < 0)
		pad = 0;
		for(i = -radius + pad;i < radius + pad;i++)
		{
			for(j = -radius + pad;j < radius + pad;j++)
			{
				c_rot = ( (float)j * cos_t - (float)i * sin_t );  
	            r_rot = ( (float)j * sin_t + (float)i * cos_t );  

				row = (int)r_rot;
				col = (int)c_rot;
				x00 = pTmpImageData[width * (row + 1) + col + 0];
				x01 = pTmpImageData[width * (row + 1) + col + 1];
				x10 = pTmpImageData[width * row + col + 0];
				x11 = pTmpImageData[width * row + col + 1];

				r_dot = (float)fabs(r_rot - row);
				c_dot = (float)fabs(c_rot - col);
				
				pixVal = x00 * (1.0f - c_dot) * (1.0f - r_dot) + x01 * (1.0f - c_dot) * r_dot 
					+ x11 * c_dot * r_dot + x10 * c_dot * (1.0f - r_dot);
				pImgRng[(i + radius) * 2 * radius + j + radius] = pixVal;
	
			}
		}
		for(i = -radius;i < radius;i++)
		{
			for(j = -radius;j < radius;j++)
			{
				c_rot = ( (float)j * cos_t - (float)i * sin_t ) / (float)hist_width;  
	            r_rot = ( (float)j * sin_t + (float)i * cos_t ) / (float)hist_width;  
	            rbin = r_rot + d / 2 - 0.5;  
	            cbin = c_rot + d / 2 - 0.5;  
	
  
	            if( rbin > -1.0  &&  rbin < d  &&  cbin > -1.0  &&  cbin < d )  
	            {
					if( rbin > 0  &&  rbin < height - 1  &&  cbin > 0  &&  cbin < width - 1 )  
					{
						dyy = pImgRng[(i + radius + 1) * radius * 2 + j + radius] - pImgRng[(i + radius - 1) * radius * 2 + j + radius];
						dxx = pImgRng[(i + radius) * radius * 2 + j + radius + 1] - pImgRng[(i + radius) * radius * 2 + j + radius - 1];
						mag = sqrt((float)(dxx * dxx + dyy * dyy));
						ang = atan2((float)dyy, (float)dxx);
	               
	                    while( ang < 0.0f )  
	                        ang += (float)PI2;  
	                    while( ang >= PI2 )  
	                        ang -= (float)PI2;  

	                    obin = ang * bins_per_rad;  
	                    w = (float)exp( -(c_rot * c_rot + r_rot * r_rot) / exp_denom );  
						//histFeat[rbin][cbin][obin] += mag * w;
	
	                    inter_hist_entry( (float *)&descF32[0], (float)rbin, cbin, (float)obin, mag * w, d, n ); 
						count++;
					}
	                
	            }
				
			}
		}
		keyPointIndex1++;
		
		sum = 0.0f;
		for(i = 0;i < 128;i++)
		{
			sum += descF32[i];
		}
		scale = 1.0f/(sum);
		for(i = 0;i < 128;i++)
		{
			descF32[i] = descF32[i] * scale;
			if(descF32[i] > MV_SIFT_DESCR_MAG_THR)
				descF32[i] = MV_SIFT_DESCR_MAG_THR;
		}
		sum = 0.0f;
		for(i = 0;i < 128;i++)
		{
			sum += descF32[i];
		}
		scale = 1.0f/(sum);
		for(i = 0;i < 128;i++)
		{
			descF32[i] = descF32[i] * scale;
		}

		for( i = 0; i < 128; i++ )  
	    {  
	        keyPointFeater[p].desc[i] = descF32[i];  
	    } 
	#if 0
	/*
		for(i = 0;i < 4;i++)
		{
			for(j = 0;j < 4;j++)
			{
				for(q = 0;q < 8;q++)
				{
					printf("%f ",descF32[i * 32 + j * 8 + q]);
				}
				printf("\n");
			}
		}
		printf("\n");
	
	if(keyPointFeater[p].layer == 0)
	{
		point.x = mvRound(keyPointFeater[p].x) >> 1;
		point.y = mvRound(keyPointFeater[p].y) >> 1;
	}
	else
	{
		point.x = mvRound(keyPointFeater[p].x) << (keyPointFeater[p].layer - 1);
		point.y = mvRound(keyPointFeater[p].y) << (keyPointFeater[p].layer - 1);
	}
		*/
		//cvCircle(cvImage0,point,4,cvScalar(255u,255u,255u));
			
	#endif
		free(pImgRng);
		
	}
	//cvShowImage("image0",cvImage0);
		//cvWaitKey(0);
	for(i = 0;i < layer;i++)
	{
		for(j = 0;j < scaleSpace;j++)
		{
			free(kernel32fPtr[i * scaleSpace + j]);
		}
	}
	free(gaussImage);
	free(dogImage);
	free(kernel32fPtr);
	free(keyPointInf);
	free(pDogImgData);
	
	return keyPointIndex1;
}
void siftTest(void)
{
	CvPoint point0,point1;
	vect2f  image_points0[500];
	vect2f  image_points1[500];
	IplImage *cvImage0,*cvImageDst;
	int pointNum0,pointNum1,matchNum = 0;

	IMAGE_t * testImage0 = mvLoadImage("D:/vc++/imagelib/imagelibtest/image/1.bmp",3);
	IMAGE_t * reSizeImage0 = mvCreateImage(1280,960,8,testImage0->nChannels);
	mvImageResize(testImage0,reSizeImage0);
	SIFT_KEYPOINTS_FEATER_t *keyPointFeater0 = (SIFT_KEYPOINTS_FEATER_t *)malloc(500 * sizeof(SIFT_KEYPOINTS_FEATER_t));
	pointNum0 = getSiftKeyPoint(reSizeImage0,keyPointFeater0,7,4,1.5,4);

	IMAGE_t * testImage1 = mvLoadImage("D:/vc++/imagelib/imagelibtest/image/2.bmp",3);
	IMAGE_t * reSizeImage1 = mvCreateImage(1280,960,8,testImage1->nChannels);
	mvImageResize(testImage1,reSizeImage1);
	SIFT_KEYPOINTS_FEATER_t *keyPointFeater1 = (SIFT_KEYPOINTS_FEATER_t *)malloc(500 * sizeof(SIFT_KEYPOINTS_FEATER_t));
	pointNum1 = getSiftKeyPoint(reSizeImage1,keyPointFeater1,7,4,1.5,4);

	//matchNum = keyPointMatch(keyPointFeater0,keyPointFeater1,image_points0,image_points1,pointNum0,pointNum1);

	IMAGE_t *srcImgeCh3 = mvCreateImage(640,480,IPL_DEPTH_8U,3);
	IMAGE_t *srcImgeCh1 = mvCreateImage(640,480,IPL_DEPTH_8U,1);
	mvImageResize(reSizeImage0,srcImgeCh3);
	mvImageBGR2Gray(srcImgeCh3,srcImgeCh1);
	IMAGE_t *dstImgeCh3 = mvCreateImage(640,480,IPL_DEPTH_8U,3);
	IMAGE_t *dstImgeCh1 = mvCreateImage(640,480,IPL_DEPTH_8U,1);
	mvImageResize(reSizeImage1,dstImgeCh3);
	mvImageBGR2Gray(dstImgeCh3,dstImgeCh1);
	
	int i,j,k;
	int idx;

	for(i = 0;i < pointNum0;i++)
	{
		float min = 100000.0f;

		for(k = 0;k < pointNum1;k++)
		{
			float sum = 0.0f;
			float *desc0 = &keyPointFeater0[i].desc[0];
			float *desc1 = &keyPointFeater1[k].desc[0];
			for(j = 0;j < 128;j++)
			{
				sum += (desc0[j] - desc1[j]) * (desc0[j] - desc1[j]);
			}
			
			if(sum < min)
			{
				idx = k;
				min = sum;
			}
		}
		keyPointFeater0[i].valid = 0;
		float *desc0 = &keyPointFeater0[i].desc[0];
		float *desc1 = &keyPointFeater1[idx].desc[0];
		if(fabs(min) < 0.003f)
		{
			/*
			for(g = 0;g < 4;g++)
			{
				for(h = 0;h < 4;h++)
				{
					for(p = 0;p < 8;p++)
					{
						printf("%f ",fabs(desc0[g * 32 + h * 8 + p] - desc1[g * 32 + h * 8 + p]));
					}
					printf("\n");
				}
			}
			*/
			printf("\n");
			
			if(keyPointFeater0[i].layer > 0)
			{
				image_points0[matchNum].x = keyPointFeater0[i].x * (float)(0x01 << (keyPointFeater0[i].layer - 1));
				image_points0[matchNum].y = keyPointFeater0[i].y * (float)(0x01 << (keyPointFeater0[i].layer - 1));
			}
			else
			{
				image_points0[matchNum].x = keyPointFeater0[i].x * 0.5f;
				image_points0[matchNum].y = keyPointFeater0[i].y * 0.5f;
			}
			if(keyPointFeater1[idx].layer > 0)
			{
				image_points1[matchNum].x = keyPointFeater1[idx].x * (float)(0x01 << (keyPointFeater1[idx].layer - 1));
				image_points1[matchNum].y = keyPointFeater1[idx].y * (float)(0x01 << (keyPointFeater1[idx].layer - 1));
			}
			else
			{
				image_points1[matchNum].x = keyPointFeater1[idx].x * 0.5f;
				image_points1[matchNum].y = keyPointFeater1[idx].y * 0.5f;
			}
			printf("%d %f %f %f %f %f\n",matchNum,min,image_points0[matchNum].x,image_points0[matchNum].y,image_points1[matchNum].x,image_points1[matchNum].y);
			printf("%f %f\n",image_points1[matchNum].x - image_points0[matchNum].x,image_points1[matchNum].y - image_points0[matchNum].y);
			cvImage0 = mvImage2CvImage(srcImgeCh1,0);
			cvImageDst = mvImage2CvImage(dstImgeCh1,0);
			
			point0.x = mvRound(image_points0[matchNum].x);
			point0.y = mvRound(image_points0[matchNum].y);
			cvCircle(cvImage0,point0,10,cvScalar(255u,255u,255u));
			
			point1.x = mvRound(image_points1[matchNum].x);
			point1.y = mvRound(image_points1[matchNum].y);
			cvCircle(cvImageDst,point1,10,cvScalar(255u,255u,255u));
			cvShowImage("image0",cvImage0);
			cvShowImage("image1",cvImageDst);
			cvWaitKey(0);

			cvReleaseImage(&cvImage0);
			cvReleaseImage(&cvImageDst);
			//printf("%f %f\n",src[matchNum].x,src[matchNum].y);
			//printf("%f %f\n",dst[matchNum].x,dst[matchNum].y);
			matchNum++;
		}

	}	
}
yhtao923
关注
用C语言实现的SIFT图像匹配算法
08-28
用C语言实现的SIFT特征匹配算法，不用配置GSL库，会C语言的人都能看懂，SIFT的函数用C语言重新实现，每一步SIFT步骤都有图像输出演示，很好的东东！
SURF特征匹配 C语言版
08-27
opencv移植，纯C语言代码（C语言实现SURF特征检测），全部开源，程序是vs2017，仅打开图片用到OpenCV。对比过OpenCV情况，匹配情况基本相同，可以供大家参考借鉴
sift特征匹配_纯c语言实现
06-14
SIFT特征匹配纯C语言代码，不调用opencv！不调用opencv！不调用opencv！结果存为txt，vs2013开发，里面有图片直接运行。
sift算法c语言实现
hjwang1的专栏
05-12 1181
转自：http://blog.csdn.net/lsh_2013/article/details/46826141 前段时间在做三维测量方面的研究，需要得到物体表面三维数据，sift算法是立体匹配中的经典算法，下面是对RobHess的SIFT源代码的注释。部分内容参考网上，在这里向各位大神表示感谢！ 头文件及函数声明 #include "sift.h" #include "imgfeatures.h" #include "utils.h" #include <cxcore.h> #inc
sift特征匹配C源码
09-14
sift算法，使用C语言实现的源码，在windows上使用vs2010上可以运行，需要配置OpenCV使用
C语言基于OpenCV的SIFT特征匹配源码解析
综上所述，该压缩包提供的是一套基于OpenCV和SIFT算法的C语言特征匹配源代码，适用于需要进行图像处理和特征提取的开发者。通过SIFT算法实现的特征匹配具有很强的抗干扰能力，能够帮助开发者在复杂的图像处理任务中...
C语言实现SIFT特征匹配源码解析
综上所述，SIFT特征匹配C源码涉及到了计算机视觉算法实现、C语言编程、OpenCV库使用和Windows平台下Visual Studio开发环境配置等多方面知识。对这些知识点的掌握有助于开发者在图像处理和计算机视觉领域中进行更深入...
SIFT特征提取与区域匹配的C语言源代码分享
标签中的“sift 特征提取 匹配 c源码”概括了本内容的三个关键点：SIFT算法、特征提取与匹配的过程、以及源代码的编程语言（C语言）。 文件名称列表中的“SIFT 分区域_区域匹配”暗示了源代码中可能包含将图像划分...
SIFT算法实战：C语言程序设计源码解析
5. 特征匹配：在两幅图像之间匹配SIFT特征点，通常使用欧氏距离或其他距离度量方法，选择距离最近的一组特征点作为匹配对。为了提高匹配的准确性，可以采用RANSAC算法对匹配点进行几何一致性检验。 C语言程序设计...
sift:使用纯C语言实现sift特征匹配(on going)
05-18
sift 本程序在Linux下编写 编译方法: 安装cmake 在终端下输入 cmake . 在终端下输入 make
sift特征点提取c
07-06
sift算法实现图像特征点提取，这是用c语言版本实现基本的提取功能。
SIFT特征提取+匹配
10-20
该程序调用opencv库函数实现的SIFT特征点提取及匹配，其中有两组图片以供测试。程序配置库为opencv2.4.9+vs2013。
opencv实现的sift算法源码，包含了图像的SIFT特征提取，详细讲解SIFT关键点的提取过程，图像间的基于SIFT特征的图像匹配
11-06
opencv实现的SIFT算法源码，包含图像的SIFT特征提取算法，以及图像之间的基于SIFT特征的匹配算法
opencv249+C语言实现SIFT影像匹配
09-13
Opencv249+C语言实现sift影像匹配、特征提取，武大数字摄影测量编程实习，内部包含可执行打包文件，可直接使用，也有源码。
SIFT算法的C++实现源代码
04-27
用C++编写的SIFT算法，非常合理的阈值，效果很好的sift实现。
在VS2010中应用SIFT（C）源码
weixin_30619101的博客
04-28 247
PS: 很久没做CV的事情了，这是很早以前刚入门时候的一篇，以后再有CV相关工作会发布在新的个人站点：http://my.phirobot.com/blog/category/cv.htmlCV分类下。 update @ 2013-03-17 15:29 from [FreedomShe] 后来才发现Rob Hess修改过源码了，下载链接貌似也改变了，现在提供以前版本源码（本人修改过后的...
基于surf的特征点匹配C语言源码
yhtao923的博客
08-10 836
代码中调用的头文件以及部分源码在之前上传的资源中。 bool SameSide(vect3f A, vect3f B, vect3f C, vect3f P) { vect3f AB = vector3Sub(B , A) ; vect3f AC = vector3Sub(C , A) ; vect3f AP = vector3Sub(P , A) ; ...
C++实现SIFT匹配
qq_42655403的博客
11-18 1677
【代码】C++实现SIFT匹配。
SITF 特征匹配 代码
hanlinger_的博客
11-01 2414
#coding=utf-8 import cv2 import scipy as sp img1 = cv2.imread('0.jpg',0) # queryImage img2 = cv2.imread('1.jpg',0) # trainImage # Initiate SIFT detector sift = cv2.xfeatures2d.SIFT_create() # find...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值