Harris Corner Detector

Harris response:

$$ R=\text{Det}(M)-k\text{Trace}(M)^{2}~\quad k \in [0.04,0.06]$$

참고: http://www.ipol.im/pub/art/2018/229/article_lr.pdf

number of corners=86

nonmaximal suppression를 하기 전 conrner의 response map

double PI = 4.0 * atan(1.0);
double Gaussian(double x, double y, double sigma) {
    double sigma2 = sigma * sigma;
    double t = (x * x + y * y) / (2 * sigma2);
    double n = 1.0 / (2.0 * PI * sigma2);
    return n * exp( -t );
};

struct Corner {
    double x, y;
    double strength;
    int    flag;
};
int HarrisCornerDetector(BYTE *img, int w, int h,
                    double sigma/*=1.2*/, double kvalue/*=0.04~0.06*/,
                    double minMeasure/*=1e-1*/,
                    std::vector<Corner>& corners)
{
    // Sobel gradient 3x3
    const int nn[] = { -w - 1, -w, -w + 1, -1, 0, 1, w - 1, w, w + 1}; // 3x3 window
    const int sobelX[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1};
    const int xmax = w - 1, ymax = h - 1;
    std::vector<double> Gx(w * h, 0);
    std::vector<double> Gy(w * h, 0);
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                int v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            Gx[pos] = sx / (4 * 255);
            Gy[pos] = sy / (4 * 255);
        }
        pos++; // skip x=xmax;
    }
    // precompute the coefficients of the Gaussian filter
    const int radius = int(2 * sigma);
    const int window = 1 + 2 * radius;
    std::vector<double> gFilter(window * window);
    for (int j = -radius, pos = 0; j <= radius; j++)
        for (int i = -radius; i <= radius; i++, pos++)
            gFilter[pos] = Gaussian(i, j, sigma);

    // compute the Harris measure for each pixel of the image
    std::vector<double> harrismap(w * h, 0);
    double hmax = 0;
    for (int y = 0, pos = 0; y < h; y++) {
        for (int x = 0; x < w; x++, pos++) {
            // Convolve gradient with Gaussian filter:
            double filtersum = 0, sxx = 0, syy = 0, sxy = 0;
            for (int yk = y - radius, wpos = 0; yk <= y + radius; yk++) {
                for (int xk = x - radius; xk <= x + radius; xk++, wpos++) {
                    if (yk < 0 || yk >= h) continue;
                    if (xk < 0 || xk >= w) continue;
                    // Gaussian weight
                    double gauss = gFilter[wpos];
                    filtersum += gauss;
                    // convolution;
                    int loc = yk * w + xk; //image location;
                    sxx += gauss * Gx[loc] * Gx[loc];
                    syy += gauss * Gy[loc] * Gy[loc];
                    sxy += gauss * Gx[loc] * Gy[loc];
                }
            }
            sxx /= filtersum;
            syy /= filtersum;
            sxy /= filtersum;
            // Harris corner measure = Det(M) - kvalue.(Trace(M))^2
            double hmeasure = sxx * syy - sxy * sxy - kvalue * (sxx + syy) * (sxx + syy);
            if (hmeasure <= 0) continue;
            if (hmeasure > hmax) hmax = hmeasure;
            harrismap[pos] = hmeasure;
        }
    }
    // log scale
    double scale = 255 / log(1 + hmax);
    for (int pos = w * h; pos-- > 0;)
        harrismap[pos] = scale * log(1 + harrismap[pos]);
    // Nonmaximal suppression: keep only local maxima;
    minMeasure *= 255; //255 = the maximum value of Harris measure(in log scale);
    const int nn8[] = { -w - 1, -w, -w + 1, 1, w + 1, w, w - 1, -1}; // 8-neighbors;
    for (int y = 1, pos = w; y < ymax; y++) {
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            // thresholding: Harris measure > epsilon
            double hmeasure = harrismap[pos];
            if (hmeasure <= minMeasure) continue;
            bool isLocalMax = true;
            for (int k = 0; k < 8; k++) {
                double hk = harrismap[pos + nn8[k]];
                if (hk >= hmeasure) {
                    isLocalMax = false;
                    break;
                }
            }
            if (isLocalMax) {
                Corner c;
                c.x = x, c.y = y, c.strength = hmeasure, c.flag = 1;
                corners.push_back(c);
            }
        }
        pos++; //skip x=w-1;
    }

    // remove corners too close to each other: keep the highest measure
    double minDistance = 8;
    std::vector<Corner>::iterator iter = corners.begin();
    while (iter != corners.end()) {
        for (int i = 0; i < corners.size(); i++) {
            Corner& c1 = corners[i];
            if (c1.flag == 0) continue;
            if ((c1.x == iter->x && c1.y == iter->y) ) continue;
            if (sqrt((c1.x-iter->x)*(c1.x-iter->x) + (c1.y-iter->y)*(c1.y-iter->y)) > minDistance) 
                continue;
            if (c1.strength < iter->strength) continue;
            iter->flag = 0;    // iter = corners.erase(iter); berased = true;
            break;
        }
        iter++;
    }
    return 1;
};