Classical Algorithm Research Series: 9: Continue: Teach you to implement the sift algorithm step by step using C language

Author: July, March 12, 2011
Source: http://blog.csdn.net/v_july_v.
Refer to the sift library maintained by Rob Hess.
Environment: windows xp + vc6.0
Condition: C language.
Note: This BLOG will implement all the classic algorithms one by one.
------------------------

This article teaches you how to use the C language to implement the sift algorithm step by step:
Function writing
Okay. Next, let's write all the functions involved in the main function one by one, which is also the key part of this article:
// Subsample the original image and return a 2x smaller image
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 For (I = 0; I <w; I ++)
Imnew (j, I) = Im (j * 2, I * 2 );
Return imnew;
}

// Sample the original image and return an image of 2x larger size
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 For (I = 0; I <w; I ++)
Imnew (j, I) = Im (j/2, I/2 );

Return imnew;
}

// Sample the original image and return a linear interpolation image of 2x larger size
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 {
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 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 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 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;
}

// Returns the gray value between pixels through bilinear interpolation.
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;
}

// Matrix Normalization
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;
}

// Vector Normalization
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;
}

// Obtain the Euclidean length of the vector, 2-norm.
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 );
}

// Generates 1D Gaussian Kernel
Float * GaussianKernel1D (float sigma, int dim)
{
 
Unsigned int I;
// Printf ("GaussianKernel1D (): Creating 1x % d vector for sigma = %. 3f gaussian kernel", 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)/(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 ("");
Return kern;
}

// Generate a 2D Gaussian Kernel Matrix
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", 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", 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", c, I, j );
Float v = m * exp (-(1.0 * I + 1.0 * 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;
}

// Convolution at pixels in the x direction
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)", x, y );
For (I = 0; I <dim; I ++)
{
Col = x + (I-cen );
If (col &