經典算法研究：手把手教你用C語言從零實現SIFT算法（下）

編輯：關於C語言

作者：July | 日期：2011年3月12日 | 接上篇繼續核心函數實現

核心工具函數實現

接續上文，逐一實現主函數調用的所有關鍵函數，這是本文核心部分。

1. 圖像縮放操作

下采樣（縮小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)]
    
    // 先填充所有像素
    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
       A,C,H,J為原始像素；B,E,G,I,F為插值像素 */
    
    // 插值 B 和 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));
    
    // 插值 E 和 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));
    
    // 插值 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;
}

2. 雙線性插值獲取亞像素灰度值

float getPixelBI(CvMat* im, float col, float row) {
    int irow = (int)row, icol = (int)col;
    float rfrac = 1.0 - (row - irow);
    float cfrac = 1.0 - (col - icol);
    int width = im->cols, height = im->rows;
    
    #define ImMat(ROW,COL) ((float*)(im->data.fl + im->step/sizeof(float)*(ROW)))[(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;
    
    float row1 = 0, row2 = 0;
    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;
}

3. 歸一化操作

矩陣歸一化

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 (unsigned int j = 0; j < mat->rows; j++)
        for (unsigned int i = 0; i < mat->cols; i++)
            Mat(j,i) /= sum;
}

向量歸一化

void normalizeVec(float* vec, int dim) {
    float sum = 0;
    for (unsigned int i = 0; i < dim; i++) sum += vec[i];
    for (unsigned int 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);
}

4. 高斯核生成

1D高斯核

float* GaussianKernel1D(float sigma, int dim) {
    unsigned int i;
    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); // 可選歸一化
    return kern;
}

2D高斯核矩陣

CvMat* GaussianKernel2D(float sigma) {
    int dim = (int)max(3.0f, 2.0 * GAUSSKERN * sigma + 1.0f);
    if (dim % 2 == 0) dim++; // 確保為奇數
    
    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;
    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++) {
            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;
}

5. 卷積運算（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, cen = dim / 2;
    
    for (i = 0; i < dim; i++) {
        col = x + (i - cen);
        if (col >= 0 && col < src->cols)
            pixel += kernel[i] * Src(y, col);
    }
    return pixel;
}