Harris Corner Detector

Harris Corner Detector는 영상에서 코너(corner)를 검출하는 대표적인 특징점(feature point) 검출 알고리즘이다. 코너는 작은 윈도를 어느 방향으로 이동하더라도 영상의 밝기 변화가 크게 나타나는 점으로 정의된다.

이를 위해 Harris 알고리즘은 각 픽셀 주변의 gradient를 이용하여 구조 텐서(structure tensor 또는 gradient covariance matrix)를 계산한다. 구조 텐서의 두 고유값이 모두 크면 모든 방향에서 밝기 변화가 크므로 코너로 판단하고, 한 고유값만 크면 에지(edge), 두 고유값이 모두 작으면 평탄한 영역으로 판단한다.

실제 구현에서는 고유값을 직접 계산하지 않고 Harris 응답 함수(response)

$$ R=\text{Det}(M)-k\text{Trace}(M)^{2}$$

$$M = \begin{pmatrix}\sum w(x) I_x^2  & \sum w(x) I_x I_y\\ \sum w(x) I_x I_y & \sum w(x) I_y ^2 \end{pmatrix}$$

를 이용하여 코너를 검출한다. 여기서 \(M\)은 구조 텐서이고, \(k\)는 보통 0.04~0.06 정도의 값을 사용한다. \(R\)이 큰 양수이면 코너, 음수이면 에지, 0에 가까우면 평탄한 영역으로 분류한다.

Harris Corner Detector는 계산량이 적고 회전에 대해 안정적이며, 영상 정합(image matching), 카메라 보정, 객체 추적 등 다양한 컴퓨터 비전 분야에서 널리 사용되는 코너 검출 알고리즘이다.

참고: 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;
};