学步园

 作者：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金字塔中的局部最大值，找到之后，还要经过两个检验才能确认为特征点：一是它必须有明显的差异，二是他不应该是边缘点，（也就是说，在极值点处的主曲率比应该小于某一个阈值）。