作者: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;
- }
#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;
- }
- }
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;
- }
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))&&