【影象處理筆記】SIFT演演算法原理與原始碼分析

class buildDoGPyramidComputer : public ParallelLoopBody
{
public:
    buildDoGPyramidComputer(
        int _nOctaveLayers,
        const std::vector<Mat>& _gpyr,
        std::vector<Mat>& _dogpyr)
        : nOctaveLayers(_nOctaveLayers),
        gpyr(_gpyr),
        dogpyr(_dogpyr) { }

    void operator()(const cv::Range& range) const CV_OVERRIDE
    {
        CV_TRACE_FUNCTION();

        const int begin = range.start;
        const int end = range.end;

        for (int a = begin; a < end; a++)
        {
            const int o = a / (nOctaveLayers + 2);//確定該影象屬於哪個Octave
            const int i = a % (nOctaveLayers + 2);//確定該影象屬於Octave的第幾層

            const Mat& src1 = gpyr[o * (nOctaveLayers + 3) + i];
            const Mat& src2 = gpyr[o * (nOctaveLayers + 3) + i + 1];
            Mat& dst = dogpyr[o * (nOctaveLayers + 2) + i];
            subtract(src2, src1, dst, noArray(), DataType<sift_wt>::type);//第i+1層減去第i層
        }
    }

private:
    int nOctaveLayers;
    const std::vector<Mat>& gpyr;
    std::vector<Mat>& dogpyr;
};

void SIFT_Impl::buildDoGPyramid(const std::vector<Mat>& gpyr, std::vector<Mat>& dogpyr) const
{
    CV_TRACE_FUNCTION();

    int nOctaves = (int)gpyr.size() / (nOctaveLayers + 3);
    dogpyr.resize(nOctaves * (nOctaveLayers + 2));//儲存所有的DoG影象
    // 平行計算，共nOctaves個Octave，每個Octave有nOctaveLayers+3層，共nOctaves * (nOctaveLayers + 2)個差分計算
    parallel_for_(Range(0, nOctaves * (nOctaveLayers + 2)), buildDoGPyramidComputer(nOctaveLayers, gpyr, dogpyr));
}

void process(const cv::Range& range)
{
    CV_TRACE_FUNCTION();

    const int begin = range.start;
    const int end = range.end;

    static const int n = SIFT_ORI_HIST_BINS;
    float CV_DECL_ALIGNED(CV_SIMD_WIDTH) hist[n];

    const Mat& img = dog_pyr[idx];
    const Mat& prev = dog_pyr[idx - 1];
    const Mat& next = dog_pyr[idx + 1];

    for (int r = begin; r < end; r++)
    {
        const sift_wt* currptr = img.ptr<sift_wt>(r);
        const sift_wt* prevptr = prev.ptr<sift_wt>(r);
        const sift_wt* nextptr = next.ptr<sift_wt>(r);
        int c = SIFT_IMG_BORDER;

        // vector loop reminder, better predictibility and less branch density
        for (; c < cols - SIFT_IMG_BORDER; c++)
        {
            sift_wt val = currptr[c];
            if (std::abs(val) <= threshold)
                continue;

            sift_wt _00, _01, _02;
            sift_wt _10, _12;
            sift_wt _20, _21, _22;
            _00 = currptr[c - step - 1]; _01 = currptr[c - step]; _02 = currptr[c - step + 1];
            _10 = currptr[c - 1];                        _12 = currptr[c + 1];
            _20 = currptr[c + step - 1]; _21 = currptr[c + step]; _22 = currptr[c + step + 1];

            bool calculate = false;
            if (val > 0)//當val大於0，找極大值
            {
                sift_wt vmax = std::max(std::max(std::max(_00, _01), std::max(_02, _10)), std::max(std::max(_12, _20), std::max(_21, _22)));
                if (val >= vmax)//當前點和8鄰域相比是否為極大值
                {
                    _00 = prevptr[c - step - 1]; _01 = prevptr[c - step]; _02 = prevptr[c - step + 1];
                    _10 = prevptr[c - 1];                        _12 = prevptr[c + 1];
                    _20 = prevptr[c + step - 1]; _21 = prevptr[c + step]; _22 = prevptr[c + step + 1];
                    vmax = std::max(std::max(std::max(_00, _01), std::max(_02, _10)), std::max(std::max(_12, _20), std::max(_21, _22)));
                    if (val >= vmax)//當前點和上一層的8鄰域相比是否為極大值
                    {
                        _00 = nextptr[c - step - 1]; _01 = nextptr[c - step]; _02 = nextptr[c - step + 1];
                        _10 = nextptr[c - 1];                        _12 = nextptr[c + 1];
                        _20 = nextptr[c + step - 1]; _21 = nextptr[c + step]; _22 = nextptr[c + step + 1];
                        vmax = std::max(std::max(std::max(_00, _01), std::max(_02, _10)), std::max(std::max(_12, _20), std::max(_21, _22)));
                        if (val >= vmax)//當前點和下一層的8鄰域相比是否為極大值
                        {
                            sift_wt _11p = prevptr[c], _11n = nextptr[c];
                            calculate = (val >= std::max(_11p, _11n));//當前點和上下兩層對應點相比是否為極大值
                        }
                    }
                }

            }
            else { // val cant be zero here (first abs took care of zero), must be negative// val要不是正的，要不是負的，不可能是0
                sift_wt vmin = std::min(std::min(std::min(_00, _01), std::min(_02, _10)), std::min(std::min(_12, _20), std::min(_21, _22)));
                if (val <= vmin)
                {
                    _00 = prevptr[c - step - 1]; _01 = prevptr[c - step]; _02 = prevptr[c - step + 1];
                    _10 = prevptr[c - 1];                        _12 = prevptr[c + 1];
                    _20 = prevptr[c + step - 1]; _21 = prevptr[c + step]; _22 = prevptr[c + step + 1];
                    vmin = std::min(std::min(std::min(_00, _01), std::min(_02, _10)), std::min(std::min(_12, _20), std::min(_21, _22)));
                    if (val <= vmin)
                    {
                        _00 = nextptr[c - step - 1]; _01 = nextptr[c - step]; _02 = nextptr[c - step + 1];
                        _10 = nextptr[c - 1];                        _12 = nextptr[c + 1];
                        _20 = nextptr[c + step - 1]; _21 = nextptr[c + step]; _22 = nextptr[c + step + 1];
                        vmin = std::min(std::min(std::min(_00, _01), std::min(_02, _10)), std::min(std::min(_12, _20), std::min(_21, _22)));
                        if (val <= vmin)
                        {
                            sift_wt _11p = prevptr[c], _11n = nextptr[c];
                            calculate = (val <= std::min(_11p, _11n));
                        }
                    }
                }
            }
    ...//後面還有求亞畫素精度極值，梯度方向直方圖
}

static bool adjustLocalExtrema(
        const std::vector<Mat>& dog_pyr, KeyPoint& kpt, int octv,
        int& layer, int& r, int& c, int nOctaveLayers,
        float contrastThreshold, float edgeThreshold, float sigma
    )
{
    CV_TRACE_FUNCTION();

    const float img_scale = 1.f / (255 * SIFT_FIXPT_SCALE);
    const float deriv_scale = img_scale * 0.5f;
    const float second_deriv_scale = img_scale;
    const float cross_deriv_scale = img_scale * 0.25f;

    float xi = 0, xr = 0, xc = 0, contr = 0;
    int i = 0;

    for (; i < SIFT_MAX_INTERP_STEPS; i++)
    {
        int idx = octv * (nOctaveLayers + 2) + layer;
        const Mat& img = dog_pyr[idx];
        const Mat& prev = dog_pyr[idx - 1];
        const Mat& next = dog_pyr[idx + 1];
        //D(x,y,σ)在樣本點點的導數，用中心差分計算
        Vec3f dD((img.at<sift_wt>(r, c + 1) - img.at<sift_wt>(r, c - 1)) * deriv_scale,
            (img.at<sift_wt>(r + 1, c) - img.at<sift_wt>(r - 1, c)) * deriv_scale,
            (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c)) * deriv_scale);
        // 基於中心差分的二階導數
        float v2 = (float)img.at<sift_wt>(r, c) * 2;
        float dxx = (img.at<sift_wt>(r, c + 1) + img.at<sift_wt>(r, c - 1) - v2) * second_deriv_scale;
        float dyy = (img.at<sift_wt>(r + 1, c) + img.at<sift_wt>(r - 1, c) - v2) * second_deriv_scale;
        float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2) * second_deriv_scale;
        float dxy = (img.at<sift_wt>(r + 1, c + 1) - img.at<sift_wt>(r + 1, c - 1) -
            img.at<sift_wt>(r - 1, c + 1) + img.at<sift_wt>(r - 1, c - 1)) * cross_deriv_scale;
        float dxs = (next.at<sift_wt>(r, c + 1) - next.at<sift_wt>(r, c - 1) -
            prev.at<sift_wt>(r, c + 1) + prev.at<sift_wt>(r, c - 1)) * cross_deriv_scale;
        float dys = (next.at<sift_wt>(r + 1, c) - next.at<sift_wt>(r - 1, c) -
            prev.at<sift_wt>(r + 1, c) + prev.at<sift_wt>(r - 1, c)) * cross_deriv_scale;
        //海森矩陣
        Matx33f H(dxx, dxy, dxs,
            dxy, dyy, dys,
            dxs, dys, dss);
        //取D(X)關於X的導數並令其為零，可求得極值的位置X, X = (x, y, σ)T
        //即，HX=-dD
        Vec3f X = H.solve(dD, DECOMP_LU);

        xi = -X[2];
        xr = -X[1];
        xc = -X[0];
        //如果由泰勒級數插值得到的三個座標的偏移量都小於0.5，說明已經找到特徵點，則退出迭代
        if (std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f)
            break;
        //如果三個座標偏移量中任意一個大於一個很大的數，則說明該極值點不是特徵點，函數返回
        if (std::abs(xi) > (float)(INT_MAX / 3) ||
            std::abs(xr) > (float)(INT_MAX / 3) ||
            std::abs(xc) > (float)(INT_MAX / 3))
            return false;
        //由上面得到的偏移量重新定義插值中心的座標位置,用於下次迭代，上面偏移量小於0.5的話，四捨五入還在layer層的(r，c)點，就不用下面的幾句
        c += cvRound(xc);
        r += cvRound(xr);
        layer += cvRound(xi);
        //如果新的座標超出了金字塔的座標範圍，則說明該極值點不是特徵點，函數返回
        if (layer < 1 || layer > nOctaveLayers ||
            c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER ||
            r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER)
            return false;
    }

    // ensure convergence of interpolation //如果迭代次數超過最大迭代次數，這個極值點不是特徵點
    if (i >= SIFT_MAX_INTERP_STEPS)
        return false;
    // 剔除具有低對比度的不穩定極值 D(X)=D+0.5*dD*X 
    {
        int idx = octv * (nOctaveLayers + 2) + layer;
        const Mat& img = dog_pyr[idx];
        const Mat& prev = dog_pyr[idx - 1];
        const Mat& next = dog_pyr[idx + 1];
        Matx31f dD((img.at<sift_wt>(r, c + 1) - img.at<sift_wt>(r, c - 1)) * deriv_scale,
            (img.at<sift_wt>(r + 1, c) - img.at<sift_wt>(r - 1, c)) * deriv_scale,
            (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c)) * deriv_scale);
        float t = dD.dot(Matx31f(xc, xr, xi));

        contr = img.at<sift_wt>(r, c) * img_scale + t * 0.5f;
        if (std::abs(contr) * nOctaveLayers < contrastThreshold)//當響應值D(X)小於contrastThreshold時剔除
            return false;
        // 剔除曲折度較大的點，利用Hessian的跡和行列式計算曲折度，和上面不同，下面的Hessian是2×2的
        // principal curvatures are computed using the trace and det of Hessian 
        float v2 = img.at<sift_wt>(r, c) * 2.f;
        float dxx = (img.at<sift_wt>(r, c + 1) + img.at<sift_wt>(r, c - 1) - v2) * second_deriv_scale;
        float dyy = (img.at<sift_wt>(r + 1, c) + img.at<sift_wt>(r - 1, c) - v2) * second_deriv_scale;
        float dxy = (img.at<sift_wt>(r + 1, c + 1) - img.at<sift_wt>(r + 1, c - 1) -
            img.at<sift_wt>(r - 1, c + 1) + img.at<sift_wt>(r - 1, c - 1)) * cross_deriv_scale;
        float tr = dxx + dyy;
        float det = dxx * dyy - dxy * dxy;
        //若行列式是負的，則不同曲折度具有不同的符號，且討論的關鍵點不可能是一個極值，因此要丟棄它。曲折度大就說明兩個特徵值相差大，就是邊緣，丟棄
        if (det <= 0 || tr * tr * edgeThreshold >= (edgeThreshold + 1) * (edgeThreshold + 1) * det)
            return false;
    }
    //儲存特徵點資訊
    kpt.pt.x = (c + xc) * (1 << octv);
    kpt.pt.y = (r + xr) * (1 << octv);
    kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5) * 255) << 16);
    kpt.size = sigma * powf(2.f, (layer + xi) / nOctaveLayers) * (1 << octv) * 2;
    kpt.response = std::abs(contr);

    return true;
}

// Computes a gradient orientation histogram at a specified pixel
static float calcOrientationHist(
        const Mat& img, Point pt, int radius,
        float sigma, float* hist, int n
)
{
    CV_TRACE_FUNCTION();

    int i, j, k, len = (radius*2+1)*(radius*2+1);

    float expf_scale = -1.f/(2.f * sigma * sigma);

    cv::utils::BufferArea area;
    float *X = 0, *Y = 0, *Mag, *Ori = 0, *W = 0, *temphist = 0;
    area.allocate(X, len, CV_SIMD_WIDTH);
    area.allocate(Y, len, CV_SIMD_WIDTH);
    area.allocate(Ori, len, CV_SIMD_WIDTH);
    area.allocate(W, len, CV_SIMD_WIDTH);
    area.allocate(temphist, n+4, CV_SIMD_WIDTH);
    area.commit();
    temphist += 2;
    Mag = X;

    for( i = 0; i < n; i++ ) //初始化直方圖
        temphist[i] = 0.f;

    for( i = -radius, k = 0; i <= radius; i++ )//兩個for迴圈，計算(radius*2+1)*(radius*2+1)個dx，dy和w
    {
        int y = pt.y + i;
        if( y <= 0 || y >= img.rows - 1 )
            continue;
        for( j = -radius; j <= radius; j++ )
        {
            int x = pt.x + j;
            if( x <= 0 || x >= img.cols - 1 )
                continue;

            float dx = (float)(img.at<sift_wt>(y, x+1) - img.at<sift_wt>(y, x-1));//計算相對值，不用除以2了
            float dy = (float)(img.at<sift_wt>(y-1, x) - img.at<sift_wt>(y+1, x));

            X[k] = dx; Y[k] = dy; W[k] = (i*i + j*j)*expf_scale;//這裡的W為高斯函數的e的指數
            k++;
        }
    }

    len = k;
    // 計算鄰域中所有元素的高斯加權值W，梯度幅角Ori 和梯度幅值Mag
    // compute gradient values, orientations and the weights over the pixel neighborhood
    cv::hal::exp32f(W, W, len);
    cv::hal::fastAtan2(Y, X, Ori, len, true);
    cv::hal::magnitude32f(X, Y, Mag, len);

    k = 0;
    for( ; k < len; k++ )//漲知識了，我每次都要寫老長的判斷，判斷是在第幾象限
    {
        int bin = cvRound((n/360.f)*Ori[k]);//判斷鄰域畫素的梯度幅角屬於n個柱體的哪一個，和Lowe推薦的方法不同，它選了最近的柱體，然後再平滑
        //如果超出範圍，則利用圓周迴圈確定其真正屬於的那個柱體，比如dy=-1，dx=1，Ori=-45°，假設n=8，bin=-1，bin+=n，就是第7個柱體
        if( bin >= n )
            bin -= n;
        if( bin < 0 )
            bin += n;
        temphist[bin] += W[k]*Mag[k];//直方圖放入的值是帶高斯加權的幅值，增強特徵點近的鄰域點對關鍵點方向的作用，並減少突變的影響
    }
    
    // smooth the histogram
    temphist[-1] = temphist[n-1];//為了圓周迴圈，提前填充好直方圖前後各兩個變數
    temphist[-2] = temphist[n-2];//公式不是要取前面兩個和後面兩個一起平滑嘛
    temphist[n] = temphist[0];//也就是說，第一個應該取最後兩個一起計算，所以把最後兩個拷到前面去
    temphist[n+1] = temphist[1];
    // 平滑直方圖
    i = 0;
    for( ; i < n; i++ )
    {
        hist[i] = (temphist[i-2] + temphist[i+2])*(1.f/16.f) +
            (temphist[i-1] + temphist[i+1])*(4.f/16.f) +
            temphist[i]*(6.f/16.f);
    }
    //計算直方圖的主峰值
    float maxval = hist[0];
    for( i = 1; i < n; i++ )
        maxval = std::max(maxval, hist[i]);

    return maxval;
}

#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

int main() {
    Mat src1 = imread("./1.jpg", 0);
    Mat src2 = imread("./2.jpg", 0);

    Mat img1 = imread("./2.jpg");
    Mat img2 = imread("./2.jpg");

    SIFT S;
    Ptr<SIFT> pSIFT = SIFT::create(55);
    vector<KeyPoint> points1, points2;
    pSIFT->detect(src1, points1);
    pSIFT->detect(src2, points2);
    drawKeypoints(src1, points1, img1, Scalar(0,255,0));
    drawKeypoints(src2, points2, img2, Scalar(0, 255, 0));
    Mat discriptions1, discriptions2;
    pSIFT->compute(src1, points1, discriptions1);
    pSIFT->compute(src2, points2, discriptions2);
    Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce");
    vector<DMatch> matches;
    matcher->match(discriptions1, discriptions2, matches);
    Mat match_img;
    drawMatches(src1, points1, src2, points2, matches, match_img);
    imshow("markImg", match_img);
    waitKey(0);
    return 0; 
}