学步园

                      教你一步一步用c语言实现sift算法、下

作者：July、二零一一年三月十二日
出处：http://blog.csdn.net/v_JULY_v。
参考：Rob Hess维护的sift 库
环境：windows xp+vc6.0
条件：c语言实现。
说明：本BLOG内会陆续一一实现所有经典算法。
------------------------

本文接上，教你一步一步用c语言实现sift算法、上，而来：
函数编写
    ok，接上文，咱们一个一个的来编写main函数中所涉及到所有函数，这也是本文的关键部分：
//下采样原来的图像，返回缩小2倍尺寸的图像<br /> CvMat * halfSizeImage(CvMat * im)<br /> {<br /> unsigned int i,j;<br /> int w = im->cols/2;<br /> int h = im->rows/2;<br /> CvMat *imnew = cvCreateMat(h, w, CV_32FC1);</p> <p>#define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]<br /> #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]<br /> for ( j = 0; j < h; j++)<br /> for ( i = 0; i < w; i++)<br /> Imnew(j,i)=Im(j*2, i*2);<br /> return imnew;<br /> }</p> <p>//上采样原来的图像，返回放大2倍尺寸的图像<br /> CvMat * doubleSizeImage(CvMat * im)<br /> {<br /> unsigned int i,j;<br /> int w = im->cols*2;<br /> int h = im->rows*2;<br /> CvMat *imnew = cvCreateMat(h, w, CV_32FC1);</p> <p>#define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]<br /> #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]</p> <p> for ( j = 0; j < h; j++)<br /> for ( i = 0; i < w; i++)<br /> Imnew(j,i)=Im(j/2, i/2);</p> <p> return imnew;<br /> }</p> <p>//上采样原来的图像，返回放大2倍尺寸的线性插值图像<br /> CvMat * doubleSizeImage2(CvMat * im)<br /> {<br /> unsigned int i,j;<br /> int w = im->cols*2;<br /> int h = im->rows*2;<br /> CvMat *imnew = cvCreateMat(h, w, CV_32FC1);</p> <p>#define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]<br /> #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]</p> <p> // fill every pixel so we don't have to worry about skipping pixels later<br /> for ( j = 0; j < h; j++)<br /> {<br /> for ( i = 0; i < w; i++)<br /> {<br /> Imnew(j,i)=Im(j/2, i/2);<br /> }<br /> }<br /> /*<br /> A B C<br /> E F G<br /> H I J<br /> pixels A C H J are pixels from original image<br /> pixels B E G I F are interpolated pixels<br /> */<br /> // interpolate pixels B and I<br /> for ( j = 0; j < h; j += 2)<br /> for ( i = 1; i < w - 1; i += 2)<br /> Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2, i/2+1));<br /> // interpolate pixels E and G<br /> for ( j = 1; j < h - 1; j += 2)<br /> for ( i = 0; i < w; i += 2)<br /> Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2+1, i/2));<br /> // interpolate pixel F<br /> for ( j = 1; j < h - 1; j += 2)<br /> for ( i = 1; i < w - 1; i += 2)<br /> Imnew(j,i)=0.25*(Im(j/2, i/2)+Im(j/2+1, i/2)+Im(j/2, i/2+1)+Im(j/2+1, i/2+1));<br /> return imnew;<br /> }</p> <p>//双线性插值，返回像素间的灰度值<br /> float getPixelBI(CvMat * im, float col, float row)<br /> {<br /> int irow, icol;<br /> float rfrac, cfrac;<br /> float row1 = 0, row2 = 0;<br /> int width=im->cols;<br /> int height=im->rows;<br /> #define ImMat(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]</p> <p> irow = (int) row;<br /> icol = (int) col;</p> <p> if (irow < 0 || irow >= height<br /> || icol < 0 || icol >= width)<br /> return 0;<br /> if (row > height - 1)<br /> row = height - 1;<br /> if (col > width - 1)<br /> col = width - 1;<br /> rfrac = 1.0 - (row - (float) irow);<br /> cfrac = 1.0 - (col - (float) icol);<br /> if (cfrac < 1)<br /> {<br /> row1 = cfrac * ImMat(irow,icol) + (1.0 - cfrac) * ImMat(irow,icol+1);<br /> }<br /> else<br /> {<br /> row1 = ImMat(irow,icol);<br /> }<br /> if (rfrac < 1)<br /> {<br /> if (cfrac < 1)<br /> {<br /> row2 = cfrac * ImMat(irow+1,icol) + (1.0 - cfrac) * ImMat(irow+1,icol+1);<br /> } else<br /> {<br /> row2 = ImMat(irow+1,icol);<br /> }<br /> }<br /> return rfrac * row1 + (1.0 - rfrac) * row2;<br /> }</p> <p>//矩阵归一化<br /> void normalizeMat(CvMat* mat)<br /> {<br /> #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]<br /> float sum = 0;</p> <p> for (unsigned int j = 0; j < mat->rows; j++)<br /> for (unsigned int i = 0; i < mat->cols; i++)<br /> sum += Mat(j,i);<br /> for ( j = 0; j < mat->rows; j++)<br /> for (unsigned int i = 0; i < mat->rows; i++)<br /> Mat(j,i) /= sum;<br /> }</p> <p>//向量归一化<br /> void normalizeVec(float* vec, int dim)<br /> {<br /> unsigned int i;<br /> float sum = 0;<br /> for ( i = 0; i < dim; i++)<br /> sum += vec[i];<br /> for ( i = 0; i < dim; i++)<br /> vec[i] /= sum;<br /> }</p> <p>//得到向量的欧式长度，2-范数<br /> float GetVecNorm( float* vec, int dim )<br /> {<br /> float sum=0.0;<br /> for (unsigned int i=0;i<dim;i++)<br /> sum+=vec[i]*vec[i];<br /> return sqrt(sum);<br /> }</p> <p>//产生1D高斯核<br /> float* GaussianKernel1D(float sigma, int dim)<br /> {</p> <p> unsigned int i;<br /> //printf("GaussianKernel1D(): Creating 1x%d vector for sigma=%.3f gaussian kernel/n", dim, sigma);</p> <p> float *kern=(float*)malloc( dim*sizeof(float) );<br /> float s2 = sigma * sigma;<br /> int c = dim / 2;<br /> float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);<br /> double v;<br /> for ( i = 0; i < (dim + 1) / 2; i++)<br /> {<br /> v = m * exp(-(1.0*i*i)/(2.0 * s2)) ;<br /> kern[c+i] = v;<br /> kern[c-i] = v;<br /> }<br /> // normalizeVec(kern, dim);<br /> // for ( i = 0; i < dim; i++)<br /> // printf("%f ", kern[i]);<br /> // printf("/n");<br /> return kern;<br /> }</p> <p>//产生2D高斯核矩阵<br /> CvMat* GaussianKernel2D(float sigma)<br /> {<br /> // int dim = (int) max(3.0f, GAUSSKERN * sigma);<br /> int dim = (int) max(3.0f, 2.0 * GAUSSKERN *sigma + 1.0f);<br /> // make dim odd<br /> if (dim % 2 == 0)<br /> dim++;<br /> //printf("GaussianKernel(): Creating %dx%d matrix for sigma=%.3f gaussian/n", dim, dim, sigma);<br /> CvMat* mat=cvCreateMat(dim, dim, CV_32FC1);<br /> #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]<br /> float s2 = sigma * sigma;<br /> int c = dim / 2;<br /> //printf("%d %d/n", mat.size(), mat[0].size());<br /> float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);<br /> for (int i = 0; i < (dim + 1) / 2; i++)<br /> {<br /> for (int j = 0; j < (dim + 1) / 2; j++)<br /> {<br /> //printf("%d %d %d/n", c, i, j);<br /> float v = m * exp(-(1.0*i*i + 1.0*j*j) / (2.0 * s2));<br /> Mat(c+i,c+j) =v;<br /> Mat(c-i,c+j) =v;<br /> Mat(c+i,c-j) =v;<br /> Mat(c-i,c-j) =v;<br /> }<br /> }<br /> // normalizeMat(mat);<br /> return mat;<br /> }</p> <p>//x方向像素处作卷积<br /> float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y)<br /> {<br /> #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]<br /> unsigned int i;<br /> float pixel = 0;<br /> int col;<br /> int cen = dim / 2;<br /> //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);<br /> for ( i = 0; i < dim; i++)<br /> {<br /> col = x + (i - cen);<br /> if (col < 0)<br /> col = 0;<br /> if (col >= src->cols)<br /> col = src->cols - 1;<br /> pixel += kernel[i] * Src(y,col);<br /> }<br /> if (pixel > 1)<br /> pixel = 1;<br /> return pixel;<br /> }</p> <p>//x方向作卷积<br /> void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst)<br /> {<br /> #define DST(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]<br /> unsigned int i,j;</p> <p> for ( j = 0; j < src->rows; j++)<br /> {<br /> for ( i = 0; i < src->cols; i++)<br /> {<br /> //printf("%d, %d/n", i, j);<br /> DST(j,i) = ConvolveLocWidth(kern, dim, src, i, j);<br /> }<br /> }<br /> }</p> <p>//y方向像素处作卷积<br /> float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y)<br /> {<br /> #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]<br /> unsigned int j;<br /> float pixel = 0;<br /> int cen = dim / 2;<br /> //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);<br /> for ( j = 0; j < dim; j++)<br /> {<br /> int row = y + (j - cen);<br /> if (row < 0)<br /> row = 0;<br /> if (row >= src->rows)<br /> row = src->rows - 1;<br /> pixel += kernel[j] * Src(row,x);<br /> }<br /> if (pixel > 1)<br /> pixel = 1;<br /> return pixel;<br /> }</p> <p>//y方向作卷积<br /> void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst)<br /> {<br /> #define Dst(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]<br /> unsigned int i,j;<br /> for ( j = 0; j < src->rows; j++)<br /> {<br /> for ( i = 0; i < src->cols; i++)<br /> {<br /> //printf("%d, %d/n", i, j);<br /> Dst(j,i) = ConvolveLocHeight(kern, dim, src, i, j);<br /> }<br /> }<br /> }</p> <p>//卷积模糊图像<br /> int BlurImage(CvMat * src, CvMat * dst, float sigma)<br /> {<br /> float* convkernel;<br /> int dim = (int) max(3.0f, 2.0 * GAUSSKERN * sigma + 1.0f);<br /> CvMat *tempMat;<br /> // make dim odd<br /> if (dim % 2 == 0)<br /> dim++;<br /> tempMat = cvCreateMat(src->rows, src->cols, CV_32FC1);<br /> convkernel = GaussianKernel1D(sigma, dim);</p> <p> Convolve1DWidth(convkernel, dim, src, tempMat);<br /> Convolve1DHeight(convkernel, dim, tempMat, dst);<br /> cvReleaseMat(&tempMat);<br /> return dim;<br /> }

五个步骤

    ok，接下来，进入重点部分，咱们依据上文介绍的sift算法的几个步骤，来一一实现这些函数。
    为了版述清晰，再贴一下，主函数，顺便再加强下对sift 算法的五个步骤的认识：

1、SIFT算法第一步：图像预处理
CvMat *ScaleInitImage(CvMat * im) ;                  //金字塔初始化
2、SIFT算法第二步：建立高斯金字塔函数
ImageOctaves* BuildGaussianOctaves(CvMat * image) ;  //建立高斯金字塔
3、SIFT算法第三步：特征点位置检测，最后确定特征点的位置
int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);
4、SIFT算法第四步：计算高斯图像的梯度方向和幅值，计算各个特征点的主方向
void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);
5、SIFT算法第五步：抽取各个特征点处的特征描述字
void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);

ok，接下来一一具体实现这几个函数：
SIFT算法第一步
    SIFT算法第一步：扩大图像，预滤波剔除噪声，得到金字塔的最底层-第一阶的第一层：
CvMat *ScaleInitImage(CvMat * im)<br /> {<br /> double sigma,preblur_sigma;<br /> CvMat *imMat;<br /> CvMat * dst;<br /> CvMat *tempMat;<br /> //首先对图像进行平滑滤波，抑制噪声<br /> imMat = cvCreateMat(im->rows, im->cols, CV_32FC1);<br /> BlurImage(im, imMat, INITSIGMA);<br /> //针对两种情况分别进行处理：初始化放大原始图像或者在原图像基础上进行后续操作<br /> //建立金字塔的最底层<br /> if (DOUBLE_BASE_IMAGE_SIZE)<br /> {<br /> tempMat = doubleSizeImage2(imMat);//对扩大两倍的图像进行二次采样，采样率为0.5，采用线性插值<br /> #define TEMPMAT(ROW,COL) ((float *)(tempMat->data.fl + tempMat->step/sizeof(float) * (ROW)))[(COL)]</p> <p> dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);<br /> preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);<br /> BlurImage(tempMat, dst, preblur_sigma); </p> <p> // The initial blurring for the first image of the first octave of the pyramid.<br /> sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );<br /> // sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);<br /> //printf("Init Sigma: %f/n", sigma);<br /> BlurImage(dst, tempMat, sigma); //得到金字塔的最底层-放大2倍的图像<br /> cvReleaseMat( &dst );<br /> return tempMat;<br /> }<br /> else<br /> {<br /> dst = cvCreateMat(im->rows, im->cols, CV_32FC1);<br /> //sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA);<br /> preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);<br /> sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );<br /> //printf("Init Sigma: %f/n", sigma);<br /> BlurImage(imMat, dst, sigma); //得到金字塔的最底层：原始图像大小<br /> return dst;<br /> }<br /> }  

SIFT算法第二步
    SIFT第二步，建立Gaussian金字塔，给定金字塔第一阶第一层图像后，计算高斯金字塔其他尺度图像，
每一阶的数目由变量SCALESPEROCTAVE决定，给定一个基本图像，计算它的高斯金字塔图像，返回外部向量是阶梯指针，内部向量是每一个阶梯内部的不同尺度图像。
//SIFT算法第二步<br /> ImageOctaves* BuildGaussianOctaves(CvMat * image)<br /> {<br /> ImageOctaves *octaves;<br /> CvMat *tempMat;<br /> CvMat *dst;<br /> CvMat *temp;</p> <p> int i,j;<br /> double k = pow(2, 1.0/((float)SCALESPEROCTAVE)); //方差倍数<br /> float preblur_sigma, initial_sigma , sigma1,sigma2,sigma,absolute_sigma,sigma_f;<br /> //计算金字塔的阶梯数目<br /> int dim = min(image->rows, image->cols);<br /> int numoctaves = (int) (log((double) dim) / log(2.0)) - 2; //金字塔阶数<br /> //限定金字塔的阶梯数<br /> numoctaves = min(numoctaves, MAXOCTAVES);<br /> //为高斯金塔和DOG金字塔分配内存<br /> octaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );<br /> DOGoctaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );</p> <p> printf("BuildGaussianOctaves(): Base image dimension is %dx%d/n", (int)(0.5*(image->cols)), (int)(0.5*(image->rows)) );<br /> printf("BuildGaussianOctaves(): Building %d octaves/n", numoctaves);</p> <p> // start with initial source image<br /> tempMat=cvCloneMat( image );<br /> // preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);<br /> initial_sigma = sqrt(2);//sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );<br /> // initial_sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);</p> <p> //在每一阶金字塔图像中建立不同的尺度图像<br /> for ( i = 0; i < numoctaves; i++)<br /> {<br /> //首先建立金字塔每一阶梯的最底层，其中0阶梯的最底层已经建立好<br /> printf("Building octave %d of dimesion (%d, %d)/n", i, tempMat->cols,tempMat->rows);<br /> //为各个阶梯分配内存<br /> octaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 3) * sizeof(ImageLevels) );<br /> DOGoctaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 2) * sizeof(ImageLevels) );<br /> //存储各个阶梯的最底层<br /> (octaves[i].Octave)[0].Level=tempMat;</p> <p> octaves[i].col=tempMat->cols;<br /> octaves[i].row=tempMat->rows;<br /> DOGoctaves[i].col=tempMat->cols;<br /> DOGoctaves[i].row=tempMat->rows;<br /> if (DOUBLE_BASE_IMAGE_SIZE)<br /> octaves[i].subsample=pow(2,i)*0.5;<br /> else<br /> octaves[i].subsample=pow(2,i);</p> <p> if(i==0)<br /> {<br /> (octaves[0].Octave)[0].levelsigma = initial_sigma;<br /> (octaves[0].Octave)[0].absolute_sigma = initial_sigma;<br /> printf("0 scale and blur sigma : %f /n", (octaves[0].subsample) * ((octaves[0].Octave)[0].absolute_sigma));<br /> }<br /> else<br /> {<br /> (octaves[i].Octave)[0].levelsigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].levelsigma;<br /> (octaves[i].Octave)[0].absolute_sigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].absolute_sigma;<br /> printf( "0 scale and blur sigma : %f /n", ((octaves[i].Octave)[0].absolute_sigma) );<br /> }<br /> sigma = initial_sigma;<br /> //建立本阶梯其他层的图像<br /> for ( j = 1; j < SCALESPEROCTAVE + 3; j++)<br /> {<br /> dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用于存储高斯层<br /> temp = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用于存储DOG层<br /> // 2 passes of 1D on original<br /> // if(i!=0)<br /> // {<br /> // sigma1 = pow(k, j - 1) * ((octaves[i-1].Octave)[j-1].levelsigma);<br /> // sigma2 = pow(k, j) * ((octaves[i].Octave)[j-1].levelsigma);<br /> // sigma = sqrt(sigma2*sigma2 - sigma1*sigma1);<br /> sigma_f= sqrt(k*k-1)*sigma;<br /> // }<br /> // else<br /> // {<br /> // sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4)*pow(k,j);<br /> // }<br /> sigma = k*sigma;<br /> absolute_sigma = sigma * (octaves[i].subsample);<br /> printf("%d scale and Blur sigma: %f /n", j, absolute_sigma);</p> <p> (octaves[i].Octave)[j].levelsigma = sigma;<br /> (octaves[i].Octave)[j].absolute_sigma = absolute_sigma;<br /> //产生高斯层<br /> int length=BlurImage((octaves[i].Octave)[j-1].Level, dst, sigma_f);//相应尺度<br /> (octaves[i].Octave)[j].levelsigmalength = length;<br /> (octaves[i].Octave)[j].Level=dst;<br /> //产生DOG层<br /> cvSub( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp, 0 );<br /> // cvAbsDiff( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp );<br /> ((DOGoctaves[i].Octave)[j-1]).Level=temp;<br /> }<br /> // halve the image size for next iteration<br /> tempMat = halfSizeImage( ( (octaves[i].Octave)[SCALESPEROCTAVE].Level ) );<br /> }<br /> return octaves;<br /> }  

SIFT算法第三步
    SIFT算法第三步，特征点位置检测，最后确定特征点的位置检测DOG金字塔中的局部最大值，找到之后，还要经过两个检验才能确认为特征点：一是它必须有明显的差异，二是他不应该是边缘点，（也就是说，在极值点处的主曲率比应该小于某一个阈值）。
//SIFT算法第三步，特征点位置检测，<br /> int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr)<br /> {<br /> //计算用于DOG极值点检测的主曲率比的阈值<br /> double curvature_threshold;<br /> curvature_threshold= ((CURVATURE_THRESHOLD + 1)*(CURVATURE_THRESHOLD + 1))/CURVATURE_THRESHOLD;<br /> #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]</p> <p> int keypoint_count = 0;<br /> for (int i=0; i<numoctaves; i++)<br /> {<br /> for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层<br /> {<br /> //在图像的有效区域内寻找具有显著性特征的局部最大值<br /> //float sigma=(GaussianPyr[i].Octave)[j].levelsigma;<br /> //int dim = (int) (max(3.0f, 2.0*GAUSSKERN *sigma + 1.0f)*0.5);<br /> int dim = (int)(0.5*((GaussianPyr[i].Octave)[j].levelsigmalength)+0.5);<br /> for (int m=dim;m<((DOGoctaves[i].row)-dim);m++)<br /> for(int n=dim;n<((DOGoctaves[i].col)-dim);n++)<br /> {<br /> if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )<br /> {</p> <p> if ( ImLevels(i,j,m,n)!=0.0 ) //1、首先是非零<br /> {<br /> float inf_val=ImLevels(i,j,m,n);<br /> if(( (inf_val <= ImLevels(i,j-1,m-1,n-1))&&<br /> (inf_val <= ImLevels(i,j-1,m ,n-1))&&<br /> (inf_val <= ImLevels(i,j-1,m+1,n-1))&&<br /> (inf_val <= ImLevels(i,j-1,m-1,n ))&&<br /> (inf_val <= ImLevels(i,j-1,m ,n ))&&<br /> (inf_val <= ImLevels(i,j-1,m+1,n ))&&<br /> (inf_val <= ImLevels(i,j-1,m-1,n+1))&&<br /> (inf_val <= ImLevels(i,j-1,m ,n+1))&&<br /> (inf_val <= ImLevels(i,j-1,m+1,n+1))&& //底层的小尺度9</p> <p> (inf_val <= ImLevels(i,j,m-1,n-1))&&<br /> (inf_val <= ImLevels(i,j,m ,n-1))&&<br /> (inf_val <= ImLevels(i,j,m+1,n-1))&&<br /> (inf_val <= ImLevels(i,j,m-1,n ))&&<br /> (inf_val <= ImLevels(i,j,m+1,n ))&&<br /> (inf_val <= ImLevels(i,j,m-1,n+1))&&<br /> (inf_val <= ImLevels(i,j,m ,n+1))&&<br /> (inf_val <= ImLevels(i,j,m+1,n+1))&& //当前层8</p> <p> (inf_val <= ImLevels(i,j+1,m-1,n-1))&&<br /> (inf_val <= ImLevels(i,j+1,m ,n-1))&&<br /> (inf_val <= ImLevels(i,j+1,m+1,n-1))&&<br /> (inf_val <= ImLevels(i,j+1,m-1,n ))&&<br /> (inf_val <= ImLevels(i,j+1,m ,n ))&&<br /> (inf_val <= ImLevels(i,j+1,m+1,n ))&&<br /> (inf_val <= ImLevels(i,j+1,m-1,n+1))&&<br /> (inf_val <= ImLevels(i,j+1,m ,n+1))&&<br /> (inf_val <= ImLevels(i,j+1,m+1,n+1)) //下一层大尺度9<br /> ) ||<br /> ( (inf_val >= ImLevels(i,j-1,m-1,n-1))&&<br /> (inf_val >= ImLevels(i,j-1,m ,n-1))&&<br /> (inf_val >= ImLevels(i,j-1,m+1,n-1))&&<br /> (inf_val >= ImLevels(i,j-1,m-1,n ))&&<br /> (inf_val >= ImLevels(i,j-1,m ,n ))&&<br /> (inf_val >= ImLevels(i,j-1,m+1,n ))&&<br /> (inf_val >= ImLevels(i,j-1,m-1,n+1))&&<br /> (inf_val >= ImLevels(i,j-1,m ,n+1))&&<br /> (inf_val >= ImLevels(i,j-1,m+1,n+1))&&</p> <p> (inf_val >= ImLevels(i,j,m-1,n-1))&&<br /> (inf_val >= ImLevels(i,j,m ,n-1))&&<br /> (inf_val >= ImLevels(i,j,m+1,n-1))&&<br /> (inf_val >= ImLevels(i,j,m-1,n ))&&<br /> (inf_val >= ImLevels(i,j,m+1,n ))&&<br /> (inf_val >= ImLevels(i,j,m-1,n+1))&&<br /> (inf_val >= ImLevels(i,j,m ,n+1))&&<br /> (inf_val >= ImLevels(i,j,m+1,n+1))&& </p> <p> (inf_val >= ImLevels(i,j+1,m-1,n-1))&&<br /> (inf_val >= ImLevels(i,j+1,m ,n-1))&&<br /> (inf_val >= ImLevels(i,j+1,m+1,n-1))&&<br /> (inf_val >= ImLevels(i,j+1,m-1,n ))&&<br /> (inf_val >= ImLevels(i,j+1,m ,n ))&&<br /> (inf_val >= ImLevels(i,j+1,m+1,n ))&&<br /> (inf_val >= ImLevels(i,j+1,m-1,n+1))&&<br /> (inf_val >= ImLevels(i,j+1,m ,n+1))&&<br /> (inf_val >= ImLevels(i,j+1,m+1,n+1))<br /> ) ) //2、满足26个中极值点<br /> {<br /> //此处可存储<br /> //然后必须具有明显的显著性，即必须大于CONTRAST_THRESHOLD=0.02<br /> if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )<br /> {<br /> //最后显著处的特征点必须具有足够的曲率比，CURVATURE_THRESHOLD=10.0，首先计算Hessian矩阵<br /> // Compute the entries of the Hessian matrix at the extrema location.<br /> /*<br /> 1 0 -1<br /> 0 0 0<br /> -1 0 1 *0.25<br /> */<br /> // Compute the trace and the determinant of the Hessian.<br /> //Tr_H = Dxx + Dyy;<br /> //Det_H = Dxx*Dyy - Dxy^2;<br /> float Dxx,Dyy,Dxy,Tr_H,Det_H,curvature_ratio;<br /> Dxx = ImLevels(i,j,m,n-1) + ImLevels(i,j,m,n+1)-2.0*ImLevels(i,j,m,n);<br /> Dyy = ImLevels(i,j,m-1,n) + ImLevels(i,j,m+1,n)-2.0*ImLevels(i,j,m,n);<br /> Dxy = ImLevels(i,j,m-1,n-1) + ImLevels(i,j,m+1,n+1) - ImLevels(i,j,m+1,n-1) - ImLevels(i,j,m-1,n+1);<br /> Tr_H = Dxx + Dyy;<br /> Det_H = Dxx*Dyy - Dxy*Dxy;<br /> // Compute the ratio of the principal curvatures.<br /> curvature_ratio = (1.0*Tr_H*Tr_H)/Det_H;<br /> if ( (Det_H>=0.0) && (curvature_ratio <= curvature_threshold) ) //最后得到最具有显著性特征的特征点<br /> {<br /> //将其存储起来，以计算后面的特征描述字<br /> keypoint_count++;<br /> Keypoint k;<br /> /* Allocate memory for the keypoint. */<br /> k = (Keypoint) malloc(sizeof(struct KeypointSt));<br /> k->next = keypoints;<br /> keypoints = k;<br /> k->row = m*(GaussianPyr[i].subsample);<br /> k->col =n*(GaussianPyr[i].subsample);<br /> k->sy = m; //行<br /> k->sx = n; //列<br /> k->octave=i;<br /> k->level=j;<br /> k->scale = (GaussianPyr[i].Octave)[j].absolute_sigma;<br /> }//if >curvature_thresh<br /> }//if >contrast<br /> }//if inf value<br /> }//if non zero<br /> }//if >contrast<br /> } //for concrete image level col<br /> }//for levels<br /> }//for octaves<br /> return keypoint_count;<br /> }</p> <p>//在图像中，显示SIFT特征点的位置<br /> void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr)<br /> {</p> <p> Keypoint p = keypoints; // p指向第一个结点<br /> while(p) // 没到表尾<br /> {<br /> cvLine( image, cvPoint((int)((p->col)-3),(int)(p->row)),<br /> cvPoint((int)((p->col)+3),(int)(p->row)), CV_RGB(255,255,0),<br /> 1, 8, 0 );<br /> cvLine( image, cvPoint((int)(p->col),(int)((p->row)-3)),<br /> cvPoint((int)(p->col),(int)((p->row)+3)), CV_RGB(255,255,0),<br /> 1, 8, 0 );<br /> // cvCircle(image,cvPoint((uchar)(p->col),(uchar)(p->row)),<br /> // (int)((GaussianPyr[p->octave].Octave)[p->level].absolute_sigma),<br /> // CV_RGB(255,0,0),1,8,0);<br /> p=p->next;<br /> }<br /> }</p> <p>// Compute the gradient direction and magnitude of the gaussian pyramid images<br /> void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr)<br /> {<br /> // ImageOctaves *mag_thresh ;<br /> mag_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );<br /> grad_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );<br /> // float sigma=( (GaussianPyr[0].Octave)[SCALESPEROCTAVE+2].absolute_sigma ) / GaussianPyr[0].subsample;<br /> // int dim = (int) (max(3.0f, 2 * GAUSSKERN *sigma + 1.0f)*0.5+0.5);<br /> #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]<br /> for (int i=0; i<numoctaves; i++)<br /> {<br /> mag_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );<br /> grad_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );<br /> for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层<br /> {<br /> CvMat *Mag = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);<br /> CvMat *Ori = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);<br /> CvMat *tempMat1 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);<br /> CvMat *tempMat2 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);<br /> cvZero(Mag);<br /> cvZero(Ori);<br /> cvZero(tempMat1);<br /> cvZero(tempMat2);<br /> #define MAG(ROW,COL) ((float *)(Mag->data.fl + Mag->step/sizeof(float) *(ROW)))[(COL)]<br /> #define ORI(ROW,COL) ((float *)(Ori->data.fl + Ori->step/sizeof(float) *(ROW)))[(COL)]<br /> #define TEMPMAT1(ROW,COL) ((float *)(tempMat1->data.fl + tempMat1->step/sizeof(float) *(ROW)))[(COL)]<br /> #define TEMPMAT2(ROW,COL) ((float *)(tempMat2->data.fl + tempMat2->step/sizeof(float) *(ROW)))[(COL)]<br /> for (int m=1;m<(GaussianPyr[i].row-1);m++)<br /> for(int n=1;n<(GaussianPyr[i].col-1);n++)<br /> {<br /> //计算幅值<br /> TEMPMAT1(m,n) = 0.5*( ImLevels(i,j,m,n+1)-ImLevels(i,j,m,n-1) ); //dx<br /> TEMPMAT2(m,n) = 0.5*( ImLevels(i,j,m+1,n)-ImLevels(i,j,m-1,n) ); //dy<br /> MAG(m,n) = sqrt(TEMPMAT1(m,n)*TEMPMAT1(m,n)+TEMPMAT2(m,n)*TEMPMAT2(m,n)); //mag<br /> //计算方向<br /> ORI(m,n) =atan( TEMPMAT2(m,n)/TEMPMAT1(m,n) );<br /> if (ORI(m,n)==CV_PI)<br /> ORI(m,n)=-CV_PI;<br /> }<br /> ((mag_pyr[i].Octave)[j-1]).Level=Mag;<br /> ((grad_pyr[i].Octave)[j-1]).Level=Ori;<br /> cvReleaseMat(&tempMat1);<br /> cvReleaseMat(&tempMat2);<br /> }//for levels<br /> }//for octaves<br /> }
SIFT算法第四步
//SIFT算法第四步：计算各个特征点的主方向，确定主方向<br /> void AssignTheMainOrientation(int numoctaves, ImageOctaves *GaussianPyr,ImageOctaves *mag_pyr,ImageOctaves *grad_pyr)<br /> {<br /> // Set up the histogram bin centers for a 36 bin histogram.<br /> int num_bins = 36;<br /> float hist_step = 2.0*PI/num_bins;<br /> float hist_orient[36];<br /> for (int i=0;i<36;i++)<br /> hist_orient[i]=-PI+i*hist_step;<br /> float sigma1=( ((GaussianPyr[0].Octave)[SCALESPEROCTAVE].absolute_sigma) ) / (GaussianPyr[0].subsample);//SCALESPEROCTAVE+2<br /> int zero_pad = (int) (max(3.0f, 2 * GAUSSKERN *sigma1 + 1.0f)*0.5+0.5);<br /> //Assign orientations to the keypoints.<br /> #define ImLevels(OCTAVES,LEVELS,ROW,COL) ((float *)((GaussianPyr[(OCTAVES)].Octave[(LEVELS)].Level)->data.fl + (GaussianPyr[(OCTAVES)].Octave[(LEVELS)].Level)->step/sizeof(float) *(ROW)))[(COL)]</p> <p> int keypoint_count = 0;<br /> Keypoint p = keypoints; // p指向第一个结点</p> <p> while(p) // 没到表尾<br /> {<br /> int i=p->octave;<br /> int j=p->level;<br /> int m=p->sy; //行<br /> int n=p->sx; //列<br /> if ((m>=zero_pad)&&(m<GaussianPyr[i].row-zero_pad)&&<br /> (n>=zero_pad)&&(n<GaussianPyr[i].col-zero_pad) )<br /> {<br /> float sigma=( ((GaussianPyr[i].Octave)[j].absolute_sigma) ) / (GaussianPyr[i].subsample);<br /> //产生二维高斯模板<br /> CvMat* mat = GaussianKernel2D( sigma );<br /> int dim=(int)(0.5 * (mat->rows));<br /> //分配用于存储Patch幅值和方向的空间<br /> #define MAT(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]</p> <p> //声明方向直方图变量<br /> double* orienthist = (double *) malloc(36 * sizeof(double));<br /> for ( int sw = 0 ; sw < 36 ; ++sw)<br /> {<br /> orienthist[sw]=0.0;<br /> }<br /> //在特征点的周围统计梯度方向<br /> for (int x=m-dim,mm=0;x<=(m+dim);x++,mm++)<br /> for(int y=n-dim,nn=0;y<=(n+dim);y++,nn++)<br /> {<br /> //计算特征点处的幅值<br /> double dx = 0.5*(ImLevels(i,j,x,y+1)-ImLevels(i,j,x,y-1)); //dx<br /> double dy = 0.5*(ImLevels(i,j,x+1,y)-ImLevels(i,j,x-1,y)); //dy<br /> double mag = sqrt(dx*dx+dy*dy); //mag<br /> //计算方向<br /> double Ori =atan( 1.0*dy/dx );<br /> int binIdx = FindClosestRotationBin(36, Ori); //得到离现有方向最近的直方块<br /> orienthist[binIdx] = orienthist[binIdx] + 1.0* mag * MAT(mm,nn);//利用高斯加权累加进直方图相应的块<br /> }<br /> // Find peaks in the orientation histogram using nonmax suppression.<br /> AverageWeakBins (orienthist, 36);<br /> // find the maximum peak in gradient orientation<br /> double maxGrad = 0.0;<br /> int maxBin = 0;<br /> for (int b = 0 ; b < 36 ; ++b)<br /> {<br /> if (orienthist[b] > maxGrad)<br /> {<br /> maxGrad = orienthist[b];<br /> maxBin = b;<br /> }<br /> }<br /> // First determine the real interpolated peak high at the maximum bin<br /> // position, which is guaranteed to be an absolute peak.<br /> double maxPeakValue=0.0;<br /> double maxDegreeCorrection=0.0;<br /> if ( (InterpolateOrientation ( orienthist[maxBin == 0 ? (36 - 1) : (maxBin - 1)],<br /> orienthist[maxBin], orienthist[(maxBin + 1) % 36],<br /> &maxDegreeCorrection, &maxPeakValue)) == false)<br /> printf("BUG: Parabola fitting broken");</p> <p> // Now that we know the maximum peak value, we can find other keypoint<br /> // orientations, which have to fulfill two criterias:<br /> //<br /> // 1. They must be a local peak themselves. Else we might add a very<br /> // similar keypoint orientation twice (imagine for example the<br /> // values: 0.4 1.0 0.8, if 1.0 is maximum peak, 0.8 is still added<br /> // with the default threshhold, but the maximum peak orientation<br /> // was already added).<br /> // 2. They must have at least peakRelThresh times the maximum peak<br /> // value.<br /> bool binIsKeypoint[36];<br /> for ( b = 0 ; b < 36 ; ++b)<br /> {<br /> binIsKeypoint[b] = false;<br /> // The maximum peak of course is<br /> if (b == maxBin)<br /> {<br /> binIsKeypoint[b] = true;<br /> continue;<br /> }<br /> // Local peaks are, too, in case they fulfill the threshhold<br /> if (orienthist[b] < (peakRelThresh * maxPeakValue))<br /> continue;<br /> int leftI = (b == 0) ? (36 - 1) : (b - 1);<br /> int rightI = (b + 1) % 36;<br /> if (orienthist[b] <= orienthist[leftI] || orienthist[b] <= orienthist[rightI])<br /> continue; // no local peak<br /> binIsKeypoint[b] = true;<br /> }<br /> // find other possible locations<br /> double oneBinRad = (2.0 * PI) / 36;<br /> for ( b = 0 ; b < 36 ; ++b)<br /> {<br /> if (binIsKeypoint[b] == false)<br /> continue;<br /> int bLeft = (b == 0) ? (36 - 1) : (b - 1);<br /> int bRight = (b + 1) % 36;<br /> // Get an interpolated peak direction and value guess.<br /> double peakValue;<br /> double degreeCorrection;</p> <p> double maxPeakValue, maxDegreeCorrection;<br /> if (InterpolateOrientation ( orienthist[maxBin == 0 ? (36 - 1) : (maxBin - 1)],<br /> orienthist[maxBin], orienthist[(maxBin + 1) % 36],<br /> °reeCorrection, &peakValue) == false)<br /> {<br /> printf("BUG: Parabola fitting broken");<br /> }</p> <p> double degree = (b + degreeCorrection) * oneBinRad - PI;<br /> if (degree < -PI)<br /> degree += 2.0 * PI;<br /> else if (degree > PI)<br /> degree -= 2.0 * PI;<br /> //存储方向，可以直接利用检测到的链表进行该步主方向的指定;<br /> //分配内存重新存储特征点<br /> Keypoint k;<br /> /* Allocate memory for the keypoint Descriptor. */<br /> k = (Keypoint) malloc(sizeof(struct KeypointSt));<br /> k->next = keyDescriptors;<br /> keyDescriptors = k;<br /> k->descrip = (float*)malloc(LEN * sizeof(float));<br /> k->row = p->row;<br /> k->col = p->col;<br /> k->sy = p->sy; //行<br /> k->sx = p->sx; //列<br /> k->octave = p->octave;<br /> k->level = p->level;<br /> k->scale = p->scale;<br /> k->ori = degree;<br /> k->mag = peakValue;<br /> }//for<br /> free(orienthist);<br /> }<br /> p=p->next;<br /> }<br /> }</p> <p>//寻找与方向直方图最近的柱，确定其index<br /> int FindClosestRotationBin (int binCount, float angle)<br /> {<br /> angle += CV_PI;<br /> angle /= 2.0 * CV_PI;<br /> // calculate the aligned bin<br /> angle *= binCount;<br /> int idx = (int) angle;<br /> if (idx == binCount)<br /> idx = 0;<br /> return (idx);<br /> }</p> <p>// Average the content of the direction bins.<br /> void AverageWeakBins (double* hist, int binCount)<br /> {<br /> // TODO: make some tests what number of passes is the best. (its clear<br /> // one is not enough, as we may have something like<br /> // ( 0.4, 0.4, 0.3, 0.4, 0.4 ))<br /> for (int sn = 0 ; sn < 2 ; ++sn)<br /> {<br /> double firstE = hist[0];<br /> double last = hist[binCount-1];<br /> for (int sw = 0 ; sw < binCount ; ++sw)<br /> {<br /> double cur = hist[sw];<br /> double next = (sw == (binCount - 1)) ? firstE : hist[(sw + 1) % binCount];<br /> hist[sw] = (last + cur + next) / 3.0;<br /> last = cur;<br /> }<br /> }<br /> }</p> <p>// Fit a parabol to the three points (-1.0 ; left), (0.0 ; middle) and<br /> // (1.0 ; right).<br /> // Formulas:<br /> // f(x) = a (x - c)^2 + b<br /> // c is the peak offset (where f'(x) is zero), b is the peak value.<br /> // In case there is an error false is returned, otherwise a correction<br /> // value between [-1 ; 1] is returned in 'degreeCorrection', where -1<br /> // means the peak is located completely at the left vector, and -0.5 just<br /> // in the middle between left and middle and > 0 to the right side. In<br /> // 'peakValue' the maximum estimated peak value is stored.<br /> bool InterpolateOrientation (double left, double middle,double right, double *degreeCorrection, double *peakValue)<br /> {<br /> double a = ((left + right) - 2.0 * middle) / 2.0; //抛物线捏合系数a<br /> // degreeCorrection = peakValue = Double.NaN;</p> <p> // Not a parabol<br /> if (a == 0.0)<br /> return false;<br /> double c = (((left - middle) / a) - 1.0) / 2.0;<br /> double b = middle - c * c * a;<br /> if (c < -0.5 || c > 0.5)<br /> return false;<br /> *degreeCorrection = c;<br /> *peakValue = b;<br /> return true;<br /> }</p> <p>//显示特征点处的主方向<br /> void DisplayOrientation (IplImage* image, ImageOctaves *GaussianPyr)<br /> {<br /> Keypoint p = keyDescriptors; // p指向第一个结点<br /> while(p) // 没到表尾<br /> {<br /> float scale=(GaussianPyr[p->octave].Octave)[p->level].absolute_sigma;<br /> float autoscale = 3.0;<br /> float uu=autoscale*scale*cos(p->ori);<br /> float vv=autoscale*scale*sin(p->ori);<br /> float x=(p->col)+uu;<br /> float y=(p->row)+vv;<br /> cvLine( image, cvPoint((int)(p->col),(int)(p->row)),<br /> cvPoint((int)x,(int)y), CV_RGB(255,255,0),<br /> 1, 8, 0 );<br /> // Arrow head parameters<br /> float alpha = 0.33; // Size of arrow head relative to the length of the vector<br /> float beta = 0.33; // Width of the base of the arrow head relative to the length</p> <p> float xx0= (p->col)+uu-alpha*(uu+beta*vv);<br /> float yy0= (p->row)+vv-alpha*(vv-beta*uu);<br /> float xx1= (p->col)+uu-alpha*(uu-beta*vv);<br /> float yy1= (p->row)+vv-alpha*(vv+beta*uu);<br /> cvLine( image, cvPoint((int)xx0,(int)yy0),<br /> cvPoint((int)x,(int)y), CV_RGB(255,255,0),<br /> 1, 8, 0 );<br /> cvLine( image, cvPoint((int)xx1,(int)yy1),<br /> cvPoint((int)x,(int)y), CV_RGB(255,255,0),<br /> 1, 8, 0 );<br /> p=p->next;<br /> }<br /> }

SIFT算法第五步
    SIFT算法第五步：抽取各个特征点处的特征描述字，确定特征点的描述字。描述字是Patch网格内梯度方向的描述，旋转网格到主方向，插值得到网格处梯度值。
一个特征点可以用2*2*8=32维的向量，也可以用4*4*8=128维的向量更精确的进行描述。
void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr)<br /> {<br /> // The orientation histograms have 8 bins<br /> float orient_bin_spacing = PI/4;<br /> float orient_angles[8]={-PI,-PI+orient_bin_spacing,-PI*0.5, -orient_bin_spacing,<br /> 0.0, orient_bin_spacing, PI*0.5, PI+orient_bin_spacing};<br /> //产生描述字中心各点坐标<br /> float *feat_grid=(float *) malloc( 2*16 * sizeof(float));<br /> for (int i=0;i<GridSpacing;i++)<br /> {<br /> for (int j=0;j<2*GridSpacing;++j,++j)<br /> {<br /> feat_grid[i*2*GridSpacing+j]=-6.0+i*GridSpacing;<br /> feat_grid[i*2*GridSpacing+j+1]=-6.0+0.5*j*GridSpacing;<br /> }<br /> }<br /> //产生网格<br /> float *feat_samples=(float *) malloc( 2*256 * sizeof(float));<br /> for ( i=0;i<4*GridSpacing;i++)<br /> {<br /> for (int j=0;j<8*GridSpacing;j+=2)<br /> {<br /> feat_samples[i*8*GridSpacing+j]=-(2*GridSpacing-0.5)+i;<br /> feat_samples[i*8*GridSpacing+j+1]=-(2*GridSpacing-0.5)+0.5*j;<br /> }<br /> }<br /> float feat_window = 2*GridSpacing;<br /> Keypoint p = keyDescriptors; // p指向第一个结点<br /> while(p) // 没到表尾<br /> {<br /> float scale=(GaussianPyr[p->octave].Octave)[p->level].absolute_sigma;</p> <p> float sine = sin(p->ori);<br /> float cosine = cos(p->ori);<br /> //计算中心点坐标旋转之后的位置<br /> float *featcenter=(float *) malloc( 2*16 * sizeof(float));<br /> for (int i=0;i<GridSpacing;i++)<br /> {<br /> for (int j=0;j<2*GridSpacing;j+=2)<br /> {<br /> float x=feat_grid[i*2*GridSpacing+j];<br /> float y=feat_grid[i*2*GridSpacing+j+1];<br /> featcenter[i*2*GridSpacing+j]=((cosine * x + sine * y) + p->sx);<br /> featcenter[i*2*GridSpacing+j+1]=((-sine * x + cosine * y) + p->sy);<br /> }<br /> }<br /> // calculate sample window coordinates (rotated along keypoint)<br /> float *feat=(float *) malloc( 2*256 * sizeof(float));<br /> for ( i=0;i<64*GridSpacing;i++,i++)<br /> {<br /> float x=feat_samples[i];<br /> float y=feat_samples[i+1];<br /> feat[i]=((cosine * x + sine * y) + p->sx);<br /> feat[i+1]=((-sine * x + cosine * y) + p->sy);<br /> }<br /> //Initialize the feature descriptor.<br /> float *feat_desc = (float *) malloc( 128 * sizeof(float));<br /> for (i=0;i<128;i++)<br /> {<br /> feat_desc[i]=0.0;<br /> // printf("%f ",feat_desc[i]);<br /> }<br /> //printf("/n");<br /> for ( i=0;i<512;++i,++i)<br /> {<br /> float x_sample = feat[i];<br /> float y_sample = feat[i+1];<br /> // Interpolate the gradient at the sample position<br /> /*<br /> 0 1 0<br /> 1 * 1<br /> 0 1 0 具体插值策略如图示<br /> */<br /> float sample12=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample-1);<br /> float sample21=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample-1, y_sample);<br /> float sample22=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample);<br /> float sample23=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample+1, y_sample);<br /> float sample32=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample+1);<br /> //float diff_x = 0.5*(sample23 - sample21);<br /> //float diff_y = 0.5*(sample32 - sample12);<br /> float diff_x = sample23 - sample21;<br /> float diff_y = sample32 - sample12;<br /> float mag_sample = sqrt( diff_x*diff_x + diff_y*diff_y );<br /> float grad_sample = atan( diff_y / diff_x );<br /> if(grad_sample == CV_PI)<br /> grad_sample = -CV_PI;<br /> // Compute the weighting for the x and y dimensions.<br /> float *x_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));<br /> float *y_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));<br /> float *pos_wght=(float *) malloc( 8*GridSpacing * GridSpacing * sizeof(float));;<br /> for (int m=0;m<32;++m,++m)<br /> {<br /> float x=featcenter[m];<br /> float y=featcenter[m+1];<br /> x_wght[m/2] = max(1 - (fabs(x - x_sample)*1.0/GridSpacing), 0);<br /> y_wght[m/2] = max(1 - (fabs(y - y_sample)*1.0/GridSpacing), 0); </p> <p> }<br /> for ( m=0;m<16;++m)<br /> for (int n=0;n<8;++n)<br /> pos_wght[m*8+n]=x_wght[m]*y_wght[m];<br /> free(x_wght);<br /> free(y_wght);<br /> //计算方向的加权，首先旋转梯度场到主方向，然后计算差异<br /> float diff[8],orient_wght[128];<br /> for ( m=0;m<8;++m)<br /> {<br /> float angle = grad_sample-(p->ori)-orient_angles[m]+CV_PI;<br /> float temp = angle / (2.0 * CV_PI);<br /> angle -= (int)(temp) * (2.0 * CV_PI);<br /> diff[m]= angle - CV_PI;<br /> }<br /> // Compute the gaussian weighting.<br /> float x=p->sx;<br /> float y=p->sy;<br /> float g = exp(-((x_sample-x)*(x_sample-x)+(y_sample-y)*(y_sample-y))/(2*feat_window*feat_window))/(2*CV_PI*feat_window*feat_window);</p> <p> for ( m=0;m<128;++m)<br /> {<br /> orient_wght[m] = max((1.0 - 1.0*fabs(diff[m%8])/orient_bin_spacing),0);<br /> feat_desc[m] = feat_desc[m] + orient_wght[m]*pos_wght[m]*g*mag_sample;<br /> }<br /> free(pos_wght);<br /> }<br /> free(feat);<br /> free(featcenter);<br /> float norm=GetVecNorm( feat_desc, 128);<br /> for (int m=0;m<128;m++)<br /> {<br /> feat_desc[m]/=norm;<br /> if (feat_desc[m]>0.2)<br /> feat_desc[m]=0.2;<br /> }<br /> norm=GetVecNorm( feat_desc, 128);<br /> for ( m=0;m<128;m++)<br /> {<br /> feat_desc[m]/=norm;<br /> printf("%f ",feat_desc[m]);<br /> }<br /> printf("/n");<br /> p->descrip = feat_desc;<br /> p=p->next;<br /> }<br /> free(feat_grid);<br /> free(feat_samples);<br /> }</p> <p>//为了显示图象金字塔，而作的图像水平拼接<br /> CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 )<br /> {<br /> int row,col;<br /> CvMat *mosaic = cvCreateMat( max(im1->rows,im2->rows),(im1->cols+im2->cols),CV_32FC1);<br /> #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]<br /> #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]<br /> #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]<br /> cvZero(mosaic);<br /> /* Copy images into mosaic1. */<br /> for ( row = 0; row < im1->rows; row++)<br /> for ( col = 0; col < im1->cols; col++)<br /> Mosaic(row,col)=Im11Mat(row,col) ;<br /> for ( row = 0; row < im2->rows; row++)<br /> for ( col = 0; col < im2->cols; col++)<br /> Mosaic(row, (col+im1->cols) )= Im22Mat(row,col) ;<br /> return mosaic;<br /> }</p> <p>//为了显示图象金字塔，而作的图像垂直拼接<br /> CvMat* MosaicVertical( CvMat* im1, CvMat* im2 )<br /> {<br /> int row,col;<br /> CvMat *mosaic = cvCreateMat(im1->rows+im2->rows,max(im1->cols,im2->cols), CV_32FC1);<br /> #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]<br /> #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]<br /> #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]<br /> cvZero(mosaic);</p> <p> /* Copy images into mosaic1. */<br /> for ( row = 0; row < im1->rows; row++)<br /> for ( col = 0; col < im1->cols; col++)<br /> Mosaic(row,col)= Im11Mat(row,col) ;<br /> for ( row = 0; row < im2->rows; row++)<br /> for ( col = 0; col < im2->cols; col++)<br /> Mosaic((row+im1->rows),col)=Im22Mat(row,col) ;</p> <p> return mosaic;<br /> }

ok，为了版述清晰，再贴一下上文所述的主函数（注，上文已贴出，此是为了版述清晰，重复造轮）：

int main( void )<br /> {<br /> //声明当前帧IplImage指针<br /> IplImage* src = NULL;<br /> IplImage* image1 = NULL;<br /> IplImage* grey_im1 = NULL;<br /> IplImage* DoubleSizeImage = NULL;</p> <p> IplImage* mosaic1 = NULL;<br /> IplImage* mosaic2 = NULL; </p> <p> CvMat* mosaicHorizen1 = NULL;<br /> CvMat* mosaicHorizen2 = NULL;<br /> CvMat* mosaicVertical1 = NULL;</p> <p> CvMat* image1Mat = NULL;<br /> CvMat* tempMat=NULL;</p> <p> ImageOctaves *Gaussianpyr;<br /> int rows,cols;</p> <p>#define Im1Mat(ROW,COL) ((float *)(image1Mat->data.fl + image1Mat->step/sizeof(float) *(ROW)))[(COL)]</p> <p> //灰度图象像素的数据结构<br /> #define Im1B(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3]<br /> #define Im1G(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+1]<br /> #define Im1R(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+2]</p> <p> storage = cvCreateMemStorage(0); </p> <p> //读取图片<br /> if( (src = cvLoadImage( "street1.jpg", 1)) == 0 ) // test1.jpg einstein.pgm back1.bmp<br /> return -1;</p> <p> //为图像分配内存<br /> image1 = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U,3);<br /> grey_im1 = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U,1);<br /> DoubleSizeImage = cvCreateImage(cvSize(2*(src->width), 2*(src->height)), IPL_DEPTH_8U,3);</p> <p> //为图像阵列分配内存，假设两幅图像的大小相同，tempMat跟随image1的大小<br /> image1Mat = cvCreateMat(src->height, src->width, CV_32FC1);<br /> //转化成单通道图像再处理<br /> cvCvtColor(src, grey_im1, CV_BGR2GRAY);<br /> //转换进入Mat数据结构,图像操作使用的是浮点型操作<br /> cvConvert(grey_im1, image1Mat);</p> <p> double t = (double)cvGetTickCount();<br /> //图像归一化<br /> cvConvertScale( image1Mat, image1Mat, 1.0/255, 0 );</p> <p> int dim = min(image1Mat->rows, image1Mat->cols);<br /> numoctaves = (int) (log((double) dim) / log(2.0)) - 2; //金字塔阶数<br /> numoctaves = min(numoctaves, MAXOCTAVES);</p> <p> //SIFT算法第一步，预滤波除噪声，建立金字塔底层<br /> tempMat = ScaleInitImage(image1Mat) ;<br /> //SIFT算法第二步，建立Guassian金字塔和DOG金字塔<br /> Gaussianpyr = BuildGaussianOctaves(tempMat) ;</p> <p> t = (double)cvGetTickCount() - t;<br /> printf( "the time of build Gaussian pyramid and DOG pyramid is %.1f/n", t/(cvGetTickFrequency()*1000.) );</p> <p>#define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]<br /> //显示高斯金字塔<br /> for (int i=0; i<numoctaves;i++)<br /> {<br /> if (i==0)<br /> {<br /> mosaicHorizen1=MosaicHorizen( (Gaussianpyr[0].Octave)[0].Level, (Gaussianpyr[0].Octave)[1].Level );<br /> for (int j=2;j<SCALESPEROCTAVE+3;j++)<br /> mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[0].Octave)[j].Level );<br /> for ( j=0;j<NUMSIZE;j++)<br /> mosaicHorizen1=halfSizeImage(mosaicHorizen1);<br /> }<br /> else if (i==1)<br /> {<br /> mosaicHorizen2=MosaicHorizen( (Gaussianpyr[1].Octave)[0].Level, (Gaussianpyr[1].Octave)[1].Level );<br /> for (int j=2;j<SCALESPEROCTAVE+3;j++)<br /> mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (Gaussianpyr[1].Octave)[j].Level );<br /> for ( j=0;j<NUMSIZE;j++)<br /> mosaicHorizen2=halfSizeImage(mosaicHorizen2);<br /> mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );<br /> }<br /> else<br /> {<br /> mosaicHorizen1=MosaicHorizen( (Gaussianpyr[i].Octave)[0].Level, (Gaussianpyr[i].Octave)[1].Level );<br /> for (int j=2;j<SCALESPEROCTAVE+3;j++)<br /> mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[i].Octave)[j].Level );<br /> for ( j=0;j<NUMSIZE;j++)<br /> mosaicHorizen1=halfSizeImage(mosaicHorizen1);<br /> mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );<br /> }<br /> }<br /> mosaic1 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height), IPL_DEPTH_8U,1);<br /> cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0, 0 );<br /> cvConvertScaleAbs( mosaicVertical1, mosaic1, 1, 0 );</p> <p> // cvSaveImage("GaussianPyramid of me.jpg",mosaic1);<br /> cvNamedWindow("mosaic1",1);<br /> cvShowImage("mosaic1", mosaic1);<br /> cvWaitKey(0);<br /> cvDestroyWindow("mosaic1");<br /> //显示DOG金字塔<br /> for ( i=0; i<numoctaves;i++)<br /> {<br /> if (i==0)<br /> {<br /> mosaicHorizen1=MosaicHorizen( (DOGoctaves[0].Octave)[0].Level, (DOGoctaves[0].Octave)[1].Level );<br /> for (int j=2;j<SCALESPEROCTAVE+2;j++)<br /> mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[0].Octave)[j].Level );<br /> for ( j=0;j<NUMSIZE;j++)<br /> mosaicHorizen1=halfSizeImage(mosaicHorizen1);<br /> }<br /> else if (i==1)<br /> {<br /> mosaicHorizen2=MosaicHorizen( (DOGoctaves[1].Octave)[0].Level, (DOGoctaves[1].Octave)[1].Level );<br /> for (int j=2;j<SCALESPEROCTAVE+2;j++)<br /> mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (DOGoctaves[1].Octave)[j].Level );<br /> for ( j=0;j<NUMSIZE;j++)<br /> mosaicHorizen2=halfSizeImage(mosaicHorizen2);<br /> mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );<br /> }<br /> else<br /> {<br /> mosaicHorizen1=MosaicHorizen( (DOGoctaves[i].Octave)[0].Level, (DOGoctaves[i].Octave)[1].Level );<br /> for (int j=2;j<SCALESPEROCTAVE+2;j++)<br /> mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[i].Octave)[j].Level );<br /> for ( j=0;j<NUMSIZE;j++)<br /> mosaicHorizen1=halfSizeImage(mosaicHorizen1);<br /> mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );<br /> }<br /> }<br /> //考虑到DOG金字塔各层图像都会有正负，所以，必须寻找最负的，以将所有图像抬高一个台阶去显示<br /> double min_val=0;<br /> double max_val=0;<br /> cvMinMaxLoc( mosaicVertical1, &min_val, &max_val,NULL, NULL, NULL );<br /> if ( min_val<0.0 )<br /> cvAddS( mosaicVertical1, cvScalarAll( (-1.0)*min_val ), mosaicVertical1, NULL );<br /> mosaic2 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height), IPL_DEPTH_8U,1);<br /> cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0/(max_val-min_val), 0 );<br /> cvConvertScaleAbs( mosaicVertical1, mosaic2, 1, 0 );</p> <p> // cvSaveImage("DOGPyramid of me.jpg",mosaic2);<br /> cvNamedWindow("mosaic1",1);<br /> cvShowImage("mosaic1", mosaic2);<br /> cvWaitKey(0);</p> <p> //SIFT算法第三步：特征点位置检测，最后确定特征点的位置<br /> int keycount=DetectKeypoint(numoctaves, Gaussianpyr);<br /> printf("the keypoints number are %d ;/n", keycount);<br /> cvCopy(src,image1,NULL);<br /> DisplayKeypointLocation( image1 ,Gaussianpyr);</p> <p> cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );<br /> cvNamedWindow("image1",1);<br /> cvShowImage("image1", DoubleSizeImage);<br /> cvWaitKey(0);<br /> cvDestroyWindow("image1");</p> <p> //SIFT算法第四步：计算高斯图像的梯度方向和幅值，计算各个特征点的主方向<br /> ComputeGrad_DirecandMag(numoctaves, Gaussianpyr);<br /> AssignTheMainOrientation( numoctaves, Gaussianpyr,mag_pyr,grad_pyr);<br /> cvCopy(src,image1,NULL);<br /> DisplayOrientation ( image1, Gaussianpyr);</p> <p> // cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );<br /> cvNamedWindow("image1",1);<br /> // cvResizeWindow("image1", 2*(image1->width), 2*(image1->height) );<br /> cvShowImage("image1", image1);<br /> cvWaitKey(0);</p> <p> //SIFT算法第五步：抽取各个特征点处的特征描述字<br /> ExtractFeatureDescriptors( numoctaves, Gaussianpyr);<br /> cvWaitKey(0);</p> <p> //销毁窗口<br /> cvDestroyWindow("image1");<br /> cvDestroyWindow("mosaic1");<br /> //释放图像<br /> cvReleaseImage(&image1);<br /> cvReleaseImage(&grey_im1);<br /> cvReleaseImage(&mosaic1);<br /> cvReleaseImage(&mosaic2);<br /> return 0;<br /> }

最后，再看一下，运行效果（图中美女为老乡+朋友，何姐08年照）：

完。 

updated

    有很多朋友都在本文评论下要求要本程序的完整源码包（注：本文代码未贴全，复制粘贴编译肯定诸多错误），但由于时隔太久，这份代码我自己也找不到了，不过，我可以提供一份sift + KD + BBF，且可以编译正确的代码供大家参考学习，有pudn帐号的朋友可以前去下载：http://www.pudn.com/downloads340/sourcecode/graph/texture_mapping/detail1486667.html（没有pudn账号的同学请加群：169056165，验证信息：sift，至群共享下载），然后用两幅不同的图片做了下匹配（当然，运行结果显示是不匹配的），效果还不错：http://weibo.com/1580904460/yDmzAEwcV#1348475194313! July、二零一二年十月十一日。