作者：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 dont 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 sum+=vec[i]*vec[i];
return sqrt(sum);
}