分享一下我老师大神的人工智能教程！零基础，通俗易懂！http://blog.csdn.net/jiangjunshow
也欢迎大家转载本篇文章。分享知识，造福人民，实现我们中华民族伟大复兴！
教你一步一步用c语言实现sift算法、下
作者：July、二零一一年三月十二日
出处：http://blog.csdn.net/v_JULY_v。
参考：Rob Hess维护的sift 库
环境：windows xp+vc6.0
条件：c语言实现。
说明：本BLOG内会陆续一一实现所有经典算法。
------------------------
本文接上，教你一步一步用c语言实现sift算法、上，而来：
函数编写
ok，接上文，咱们一个一个的来编写main函数中所涉及到所有函数，这也是本文的关键部分：
- //下采样原来的图像，返回缩小2倍尺寸的图像
- CvMat * halfSizeImage(CvMat * im)
- {
- unsigned int i,j;
- int w = im->cols/2;
- int h = im->rows/2;
- CvMat *imnew = cvCreateMat(h, w, CV_32FC1);
- #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]
- #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]
- for ( j = 0; j < h; j++)
- for ( i = 0; i < w; i++)
- Imnew(j,i)=Im(j*2, i*2);
- return imnew;
- }
- //上采样原来的图像，返回放大2倍尺寸的图像
- CvMat * doubleSizeImage(CvMat * im)
- {
- unsigned int i,j;
- int w = im->cols*2;
- int h = im->rows*2;
- CvMat *imnew = cvCreateMat(h, w, CV_32FC1);
- #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]
- #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]
- for ( j = 0; j < h; j++)
- for ( i = 0; i < w; i++)
- Imnew(j,i)=Im(j/2, i/2);
- return imnew;
- }
- //上采样原来的图像，返回放大2倍尺寸的线性插值图像
- CvMat * doubleSizeImage2(CvMat * im)
- {
- unsigned int i,j;
- int w = im->cols*2;
- int h = im->rows*2;
- CvMat *imnew = cvCreateMat(h, w, CV_32FC1);
- #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]
- #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]
- // fill every pixel so we don't have to worry about skipping pixels later
- for ( j = 0; j < h; j++)
- {
- for ( i = 0; i < w; i++)
- {
- Imnew(j,i)=Im(j/2, i/2);
- }
- }
- /*
- A B C
- E F G
- H I J
- pixels A C H J are pixels from original image
- pixels B E G I F are interpolated pixels
- */
- // interpolate pixels B and I
- for ( j = 0; j < h; j += 2)
- for ( i = 1; i < w - 1; i += 2)
- Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2, i/2+1));
- // interpolate pixels E and G
- for ( j = 1; j < h - 1; j += 2)
- for ( i = 0; i < w; i += 2)
- Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2+1, i/2));
- // interpolate pixel F
- for ( j = 1; j < h - 1; j += 2)
- for ( i = 1; i < w - 1; i += 2)
- 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));
- return imnew;
- }
- //双线性插值，返回像素间的灰度值
- float getPixelBI(CvMat * im, float col, float row)
- {
- int irow, icol;
- float rfrac, cfrac;
- float row1 = 0, row2 = 0;
- int width=im->cols;
- int height=im->rows;
- #define ImMat(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]
- irow = (int) row;
- icol = (int) col;
- if (irow < 0 || irow >= height
- || icol < 0 || icol >= width)
- return 0;
- if (row > height - 1)
- row = height - 1;
- if (col > width - 1)
- col = width - 1;
- rfrac = 1.0 - (row - (float) irow);
- cfrac = 1.0 - (col - (float) icol);
- if (cfrac < 1)
- {
- row1 = cfrac * ImMat(irow,icol) + (1.0 - cfrac) * ImMat(irow,icol+1);
- }
- else
- {
- row1 = ImMat(irow,icol);
- }
- if (rfrac < 1)
- {
- if (cfrac < 1)
- {
- row2 = cfrac * ImMat(irow+1,icol) + (1.0 - cfrac) * ImMat(irow+1,icol+1);
- } else
- {
- row2 = ImMat(irow+1,icol);
- }
- }
- return rfrac * row1 + (1.0 - rfrac) * row2;
- }
- //矩阵归一化
- void normalizeMat(CvMat* mat)
- {
- #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]
- float sum = 0;
- for (unsigned int j = 0; j < mat->rows; j++)
- for (unsigned int i = 0; i < mat->cols; i++)
- sum += Mat(j,i);
- for ( j = 0; j < mat->rows; j++)
- for (unsigned int i = 0; i < mat->rows; i++)
- Mat(j,i) /= sum;
- }
- //向量归一化
- void normalizeVec(float* vec, int dim)
- {
- unsigned int i;
- float sum = 0;
- for ( i = 0; i < dim; i++)
- sum += vec[i];
- for ( i = 0; i < dim; i++)
- vec[i] /= sum;
- }
- //得到向量的欧式长度，2-范数
- float GetVecNorm( float* vec, int dim )
- {
- float sum=0.0;
- for (unsigned int i=0;i<dim;i++)
- sum+=vec[i]*vec[i];
- return sqrt(sum);
- }
- //产生1D高斯核
- float* GaussianKernel1D(float sigma, int dim)
- {
- unsigned int i;
- //printf("GaussianKernel1D(): Creating 1x%d vector for sigma=%.3f gaussian kernel/n", dim, sigma);
- float *kern=(float*)malloc( dim*sizeof(float) );
- float s2 = sigma * sigma;
- int c = dim / 2;
- float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);
- double v;
- for ( i = 0; i < (dim + 1) / 2; i++)
- {
- v = m * exp(-(1.0*i*i)/(2.0 * s2)) ;
- kern[c+i] = v;
- kern[c-i] = v;
- }
- // normalizeVec(kern, dim);
- // for ( i = 0; i < dim; i++)
- // printf("%f ", kern[i]);
- // printf("/n");
- return kern;
- }
- //产生2D高斯核矩阵
- CvMat* GaussianKernel2D(float sigma)
- {
- // int dim = (int) max(3.0f, GAUSSKERN * sigma);
- int dim = (int) max(3.0f, 2.0 * GAUSSKERN *sigma + 1.0f);
- // make dim odd
- if (dim % 2 == 0)
- dim++;
- //printf("GaussianKernel(): Creating %dx%d matrix for sigma=%.3f gaussian/n", dim, dim, sigma);
- CvMat* mat=cvCreateMat(dim, dim, CV_32FC1);
- #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]
- float s2 = sigma * sigma;
- int c = dim / 2;
- //printf("%d %d/n", mat.size(), mat[0].size());
- float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);
- for (int i = 0; i < (dim + 1) / 2; i++)
- {
- for (int j = 0; j < (dim + 1) / 2; j++)
- {
- //printf("%d %d %d/n", c, i, j);
- float v = m * exp(-(1.0*i*i + 1.0*j*j) / (2.0 * s2));
- Mat(c+i,c+j) =v;
- Mat(c-i,c+j) =v;
- Mat(c+i,c-j) =v;
- Mat(c-i,c-j) =v;
- }
- }
- // normalizeMat(mat);
- return mat;
- }
- //x方向像素处作卷积
- float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y)
- {
- #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]
- unsigned int i;
- float pixel = 0;
- int col;
- int cen = dim / 2;
- //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);
- for ( i = 0; i < dim; i++)
- {
- col = x + (i - cen);
- if (col < 0)
- col = 0;
- if (col >= src->cols)
- col = src->cols - 1;
- pixel += kernel[i] * Src(y,col);
- }
- if (pixel > 1)
- pixel = 1;
- return pixel;
- }
- //x方向作卷积
- void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst)
- {
- #define DST(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]
- unsigned int i,j;
- for ( j = 0; j < src->rows; j++)
- {
- for ( i = 0; i < src->cols; i++)
- {
- //printf("%d, %d/n", i, j);
- DST(j,i) = ConvolveLocWidth(kern, dim, src, i, j);
- }
- }
- }
- //y方向像素处作卷积
- float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y)
- {
- #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]
- unsigned int j;
- float pixel = 0;
- int cen = dim / 2;
- //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);
- for ( j = 0; j < dim; j++)
- {
- int row = y + (j - cen);
- if (row < 0)
- row = 0;
- if (row >= src->rows)
- row = src->rows - 1;
- pixel += kernel[j] * Src(row,x);
- }
- if (pixel > 1)
- pixel = 1;
- return pixel;
- }
- //y方向作卷积
- void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst)
- {
- #define Dst(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]
- unsigned int i,j;
- for ( j = 0; j < src->rows; j++)
- {
- for ( i = 0; i < src->cols; i++)
- {
- //printf("%d, %d/n", i, j);
- Dst(j,i) = ConvolveLocHeight(kern, dim, src, i, j);
- }
- }
- }
- //卷积模糊图像
- int BlurImage(CvMat * src, CvMat * dst, float sigma)
- {
- float* convkernel;
- int dim = (int) max(3.0f, 2.0 * GAUSSKERN * sigma + 1.0f);
- CvMat *tempMat;
- // make dim odd
- if (dim % 2 == 0)
- dim++;
- tempMat = cvCreateMat(src->rows, src->cols, CV_32FC1);
- convkernel = GaussianKernel1D(sigma, dim);
- Convolve1DWidth(convkernel, dim, src, tempMat);
- Convolve1DHeight(convkernel, dim, tempMat, dst);
- cvReleaseMat(&tempMat);
- return dim;
- }
五个步骤
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)
- {
- double sigma,preblur_sigma;
- CvMat *imMat;
- CvMat * dst;
- CvMat *tempMat;
- //首先对图像进行平滑滤波，抑制噪声
- imMat = cvCreateMat(im->rows, im->cols, CV_32FC1);
- BlurImage(im, imMat, INITSIGMA);
- //针对两种情况分别进行处理：初始化放大原始图像或者在原图像基础上进行后续操作
- //建立金字塔的最底层
- if (DOUBLE_BASE_IMAGE_SIZE)
- {
- tempMat = doubleSizeImage2(imMat);//对扩大两倍的图像进行二次采样，采样率为0.5，采用线性插值
- #define TEMPMAT(ROW,COL) ((float *)(tempMat->data.fl + tempMat->step/sizeof(float) * (ROW)))[(COL)]
- dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);
- preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);
- BlurImage(tempMat, dst, preblur_sigma);
- // The initial blurring for the first image of the first octave of the pyramid.
- sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );
- // sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);
- //printf("Init Sigma: %f/n", sigma);
- BlurImage(dst, tempMat, sigma); //得到金字塔的最底层-放大2倍的图像
- cvReleaseMat( &dst );
- return tempMat;
- }
- else
- {
- dst = cvCreateMat(im->rows, im->cols, CV_32FC1);
- //sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA);
- preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);
- sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );
- //printf("Init Sigma: %f/n", sigma);
- BlurImage(imMat, dst, sigma); //得到金字塔的最底层：原始图像大小
- return dst;
- }
- }
SIFT算法第二步
SIFT第二步，建立Gaussian金字塔，给定金字塔第一阶第一层图像后，计算高斯金字塔其他尺度图像，
每一阶的数目由变量SCALESPEROCTAVE决定，给定一个基本图像，计算它的高斯金字塔图像，返回外部向量是阶梯指针，内部向量是每一个阶梯内部的不同尺度图像。
- //SIFT算法第二步
- ImageOctaves* BuildGaussianOctaves(CvMat * image)
- {
- ImageOctaves *octaves;
- CvMat *tempMat;
- CvMat *dst;
- CvMat *temp;
- int i,j;
- double k = pow(2, 1.0/((float)SCALESPEROCTAVE)); //方差倍数
- float preblur_sigma, initial_sigma , sigma1,sigma2,sigma,absolute_sigma,sigma_f;
- //计算金字塔的阶梯数目
- int dim = min(image->rows, image->cols);
- int numoctaves = (int) (log((double) dim) / log(2.0)) - 2; //金字塔阶数
- //限定金字塔的阶梯数
- numoctaves = min(numoctaves, MAXOCTAVES);
- //为高斯金塔和DOG金字塔分配内存
- octaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );
- DOGoctaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );
- printf("BuildGaussianOctaves(): Base image dimension is %dx%d/n", (int)(0.5*(image->cols)), (int)(0.5*(image->rows)) );
- printf("BuildGaussianOctaves(): Building %d octaves/n", numoctaves);
- // start with initial source image
- tempMat=cvCloneMat( image );
- // preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);
- initial_sigma = sqrt(2);//sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );
- // initial_sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);
- //在每一阶金字塔图像中建立不同的尺度图像
- for ( i = 0; i < numoctaves; i++)
- {
- //首先建立金字塔每一阶梯的最底层，其中0阶梯的最底层已经建立好
- printf("Building octave %d of dimesion (%d, %d)/n", i, tempMat->cols,tempMat->rows);
- //为各个阶梯分配内存
- octaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 3) * sizeof(ImageLevels) );
- DOGoctaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 2) * sizeof(ImageLevels) );
- //存储各个阶梯的最底层
- (octaves[i].Octave)[0].Level=tempMat;
- octaves[i].col=tempMat->cols;
- octaves[i].row=tempMat->rows;
- DOGoctaves[i].col=tempMat->cols;
- DOGoctaves[i].row=tempMat->rows;
- if (DOUBLE_BASE_IMAGE_SIZE)
- octaves[i].subsample=pow(2,i)*0.5;
- else
- octaves[i].subsample=pow(2,i);
- if(i==0)
- {
- (octaves[0].Octave)[0].levelsigma = initial_sigma;
- (octaves[0].Octave)[0].absolute_sigma = initial_sigma;
- printf("0 scale and blur sigma : %f /n", (octaves[0].subsample) * ((octaves[0].Octave)[0].absolute_sigma));
- }
- else
- {
- (octaves[i].Octave)[0].levelsigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].levelsigma;
- (octaves[i].Octave)[0].absolute_sigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].absolute_sigma;
- printf( "0 scale and blur sigma : %f /n", ((octaves[i].Octave)[0].absolute_sigma) );
- }
- sigma = initial_sigma;
- //建立本阶梯其他层的图像
- for ( j = 1; j < SCALESPEROCTAVE + 3; j++)
- {
- dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用于存储高斯层
- temp = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用于存储DOG层
- // 2 passes of 1D on original
- // if(i!=0)
- // {
- // sigma1 = pow(k, j - 1) * ((octaves[i-1].Octave)[j-1].levelsigma);
- // sigma2 = pow(k, j) * ((octaves[i].Octave)[j-1].levelsigma);
- // sigma = sqrt(sigma2*sigma2 - sigma1*sigma1);
- sigma_f= sqrt(k*k-1)*sigma;
- // }
- // else
- // {
- // sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4)*pow(k,j);
- // }
- sigma = k*sigma;
- absolute_sigma = sigma * (octaves[i].subsample);
- printf("%d scale and Blur sigma: %f /n", j, absolute_sigma);
- (octaves[i].Octave)[j].levelsigma = sigma;
- (octaves[i].Octave)[j].absolute_sigma = absolute_sigma;
- //产生高斯层
- int length=BlurImage((octaves[i].Octave)[j-1].Level, dst, sigma_f);//相应尺度
- (octaves[i].Octave)[j].levelsigmalength = length;
- (octaves[i].Octave)[j].Level=dst;
- //产生DOG层
- cvSub( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp, 0 );
- // cvAbsDiff( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp );
- ((DOGoctaves[i].Octave)[j-1]).Level=temp;
- }
- // halve the image size for next iteration
- tempMat = halfSizeImage( ( (octaves[i].Octave)[SCALESPEROCTAVE].Level ) );
- }
- return octaves;
- }
SIFT算法第三步
SIFT算法第三步，特征点位置检测，最后确定特征点的位置检测DOG金字塔中的局部最大值，找到之后，还要经过两个检验才能确认为特征点：一是它必须有明显的差异，二是他不应该是边缘点，（也就是说，在极值点处的主曲率比应该小于某一个阈值）。
- //SIFT算法第三步，特征点位置检测，
- int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr)
- {
- //计算用于DOG极值点检测的主曲率比的阈值
- double curvature_threshold;
- curvature_threshold= ((CURVATURE_THRESHOLD + 1)*(CURVATURE_THRESHOLD + 1))/CURVATURE_THRESHOLD;
- #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]
- int keypoint_count = 0;
- for (int i=0; i<numoctaves; i++)
- {
- for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层
- {
- //在图像的有效区域内寻找具有显著性特征的局部最大值
- //float sigma=(GaussianPyr[i].Octave)[j].levelsigma;
- //int dim = (int) (max(3.0f, 2.0*GAUSSKERN *sigma + 1.0f)*0.5);
- int dim = (int)(0.5*((GaussianPyr[i].Octave)[j].levelsigmalength)+0.5);
- for (int m=dim;m<((DOGoctaves[i].row)-dim);m++)
- for(int n=dim;n<((DOGoctaves[i].col)-dim);n++)
- {
- if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )
- {
- if ( ImLevels(i,j,m,n)!=0.0 ) //1、首先是非零
- {
- float inf_val=ImLevels(i,j,m,n);
- if(( (inf_val <= ImLevels(i,j-1,m-1,n-1))&&
- (inf_val <= ImLevels(i,j-1,m ,n-1))&&
- (inf_val <= ImLevels(i,j-1,m+1,n-1))&&
- (inf_val <= ImLevels(i,j-1,m-1,n ))&&
- (inf_val <= ImLevels(i,j-1,m ,n ))&&
- (inf_val <= ImLevels(i,j-1,m+1,n ))&&
- (inf_val <= ImLevels(i,j-1,m-1,n+1))&&
- (inf_val <= ImLevels(i,j-1,m ,n+1))&&
- (inf_val <= ImLevels(i,j-1,m+1,n+1))&& //底层的小尺度9
- (inf_val <= ImLevels(i,j,m-1,n-1))&&
- (inf_val <= ImLevels(i,j,m ,n-1))&&
- (inf_val <= ImLevels(i,j,m+1,n-1))&&
- (inf_val <= ImLevels(i,j,m-1,n ))&&
- (inf_val <= ImLevels(i,j,m+1,n ))&&
- (inf_val <= ImLevels(i,j,m-1,n+1))&&
- (inf_val <= ImLevels(i,j,m ,n+1))&&
- (inf_val <= ImLevels(i,j,m+1,n+1))&& //当前层8
- (inf_val <= ImLevels(i,j+1,m-1,n-1))&&
- (inf_val <= ImLevels(i,j+1,m ,n-1))&&
- (inf_val <= ImLevels(i,j+1,m+1,n-1))&&
- (inf_val <= ImLevels(i,j+1,m-1,n ))&&
- (inf_val <= ImLevels(i,j+1,m ,n ))&&
- (inf_val <= ImLevels(i,j+1,m+1,n ))&&
- (inf_val <= ImLevels(i,j+1,m-1,n+1))&&
- (inf_val <= ImLevels(i,j+1,m ,n+1))&&
- (inf_val <= ImLevels(i,j+1,m+1,n+1)) //下一层大尺度9
- ) ||
- ( (inf_val >= ImLevels(i,j-1,m-1,n-1))&&
- (inf_val >= ImLevels(i,j-1,m ,n-1))&&
- (inf_val >= ImLevels(i,j-1,m+1,n-1))&&
- (inf_val >= ImLevels(i,j-1,m-1,n ))&&
- (inf_val >= ImLevels(i,j-1,m ,n ))&&
- (inf_val >= ImLevels(i,j-1,m+1,n ))&&
- (inf_val >= ImLevels(i,j-1,m-1,n+1))&&
- (inf_val >= ImLevels(i,j-1,m ,n+1))&&
- (inf_val >= ImLevels(i,j-1,m+1,n+1))&&
- (inf_val >= ImLevels(i,j,m-1,n-1))&&
- (inf_val >= ImLevels(i,j,m ,n-1))&&
- (inf_val >= ImLevels(i,j,m+1,n-1))&&
- (inf_val >= ImLevels(i,j,m-1,n ))&&
- (inf_val >= ImLevels(i,j,m+1,n ))&&
- (inf_val >= ImLevels(i,j,m-1,n+1))&&
- (inf_val >= ImLevels(i,j,m ,n+1))&&
- (inf_val >= ImLevels(i,j,m+1,n+1))&&
- (inf_val >= ImLevels(i,j+1,m-1,n-1))&&
- (inf_val >= ImLevels(i,j+1,m ,n-1))&&
- (inf_val >= ImLevels(i,j+1,m+1,n-1))&&
- (inf_val >= ImLevels(i,j+1,m-1,n ))&&
- (inf_val >= ImLevels(i,j+1,m ,n ))&&
- (inf_val >= ImLevels(i,j+1,m+1,n ))&&
- (inf_val >= ImLevels(i,j+1,m-1,n+1))&&
- (inf_val >= ImLevels(i,j+1,m ,n+1))&&
- (inf_val >= ImLevels(i,j+1,m+1,n+1))
- ) ) //2、满足26个中极值点
- {
- //此处可存储
- //然后必须具有明显的显著性，即必须大于CONTRAST_THRESHOLD=0.02
- if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )
- {
- //最后显著处的特征点必须具有足够的曲率比，CURVATURE_THRESHOLD=10.0，首先计算Hessian矩阵
- // Compute the entries of the Hessian matrix at the extrema location.
- /*
- 1 0 -1
- 0 0 0
- -1 0 1 *0.25
- */
- // Compute the trace and the determinant of the Hessian.
- //Tr_H = Dxx + Dyy;
- //Det_H = Dxx*Dyy - Dxy^2;
- float Dxx,Dyy,Dxy,Tr_H,Det_H,curvature_ratio;
- Dxx = ImLevels(i,j,m,n-1) + ImLevels(i,j,m,n+1)-2.0*ImLevels(i,j,m,n);
- Dyy = ImLevels(i,j,m-1,n) + ImLevels(i,j,m+1,n)-2.0*ImLevels(i,j,m,n);
- 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);
- Tr_H = Dxx + Dyy;
- Det_H = Dxx*Dyy - Dxy*Dxy;
- // Compute the ratio of the principal curvatures.
- curvature_ratio = (1.0*Tr_H*Tr_H)/Det_H;
- if ( (Det_H>=0.0) && (curvature_ratio <= curvature_threshold) ) //最后得到最具有显著性特征的特征点
- {
- //将其存储起来，以计算后面的特征描述字
- keypoint_count++;
- Keypoint k;
- /* Allocate memory for the keypoint. */
- k = (Keypoint) malloc(sizeof(struct KeypointSt));
- k->next = keypoints;
- keypoints = k;
- k->row = m*(GaussianPyr[i].subsample);
- k->col =n*(GaussianPyr[i].subsample);
- k->sy = m; //行
- k->sx = n; //列
- k->octave=i;
- k->level=j;
- k->scale = (GaussianPyr[i].Octave)[j].absolute_sigma;
- }//if >curvature_thresh
- }//if >contrast
- }//if inf value
- }//if non zero
- }//if >contrast
- } //for concrete image level col
- }//for levels
- }//for octaves
- return keypoint_count;
- }
- //在图像中，显示SIFT特征点的位置
- void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr)
- {
- Keypoint p = keypoints; // p指向第一个结点
- while(p) // 没到表尾
- {
- cvLine( image, cvPoint((int)((p->col)-3),(int)(p->row)),
- cvPoint((int)((p->col)+3),(int)(p->row)), CV_RGB(255,255,0),
- 1, 8, 0 );
- cvLine( image, cvPoint((int)(p->col),(int)((p->row)-3)),
- cvPoint((int)(p->col),(int)((p->row)+3)), CV_RGB(255,255,0),
- 1, 8, 0 );
- // cvCircle(image,cvPoint((uchar)(p->col),(uchar)(p->row)),
- // (int)((GaussianPyr[p->octave].Octave)[p->level].absolute_sigma),
- // CV_RGB(255,0,0),1,8,0);
- p=p->next;
- }
- }
- // Compute the gradient direction and magnitude of the gaussian pyramid images
- void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr)
- {
- // ImageOctaves *mag_thresh ;
- mag_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );
- grad_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );
- // float sigma=( (GaussianPyr[0].Octave)[SCALESPEROCTAVE+2].absolute_sigma ) / GaussianPyr[0].subsample;
- // int dim = (int) (max(3.0f, 2 * GAUSSKERN *sigma + 1.0f)*0.5+0.5);
- #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]
- for (int i=0; i<numoctaves; i++)
- {
- mag_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );
- grad_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );
- for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中间的scaleperoctave个层
- {
- CvMat *Mag = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);
- CvMat *Ori = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);
- CvMat *tempMat1 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);
- CvMat *tempMat2 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);
- cvZero(Mag);
- cvZero(Ori);
- cvZero(tempMat1);
- cvZero(tempMat2);
- #define MAG(ROW,COL) ((float *)(Mag->data.fl + Mag->step/sizeof(float) *(ROW)))[(COL)]
- #define ORI(ROW,COL) ((float *)(Ori->data.fl + Ori->step/sizeof(float) *(ROW)))[(COL)]
- #define TEMPMAT1(ROW,COL) ((float *)(tempMat1->data.fl + tempMat1->step/sizeof(float) *(ROW)))[(COL)]
- #define TEMPMAT2(ROW,COL) ((float *)(tempMat2->data.fl + tempMat2->step/sizeof(float) *(ROW)))[(COL)]
- for (int m=1;m<(GaussianPyr[i].row-1);m++)
- for(int n=1;n<(GaussianPyr[i].col-1);n++)
- {
- //计算幅值
- TEMPMAT1(m,n) = 0.5*( ImLevels(i,j,m,n+1)-ImLevels(i,j,m,n-1) ); //dx
- &nb
给我老师的人工智能教程打call！http://blog.csdn.net/jiangjunshow
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
