등고선을 그리는 subroutine

CONREC (Contouring for Rectangularly Spaced Data)는 Paul Bourke 가 소개한 등고선(contour) 생성 알고리즘으로, 직사각형 격자 위에 주어진 데이터 값으로부터 등고선(level)을 계산하는 방법이다.
(ref: https://paulbourke.net/papers/conrec/)

1. 격자를 작은 사각형 셀(cell)로 분할: 데이터가 \(M\times N\) 격자에 주어져 있으면, 인접한 네 개의 격자점을 하나의 셀로 간주하여 차례대로 처리한다.
2. 각 셀을 네 개의 삼각형으로 분할: 셀의 중심점을 추가하고, 중심점과 네 꼭짓점을 연결하여 네 개의 삼각형을 만든다. 중심점의 함수값은 네 꼭짓점 값의 평균을 사용한다.
3. 각 등고선 레벨과 삼각형 교차 검사: 삼각형의 세 꼭짓점 함수값이 등고선 레벨보다 큰지 작은지를 판별하여 가능한 경우를 분류한다.
4. 선형보간으로 교점 찾기: 한 변의 양 끝 \(\mathbf{p}_1\), \(\mathbf{p}_2\)를 잇는 선분과 등고선이 교차하는 경우, 교차점 \(\mathbf{p}\)는 선형보간으로 구한다. 변의 양 끝에서 함수값이 \(f_1, f_2\)이고, 등고선 레벨이 \(c\)이면 선형보간 매개변수는 \[ t= \frac{c-f_1}{f_2- f_1}\]로 찾을 수 있고, 교점은 \[\mathbf{p} = (1-t) \mathbf{p}_1 + t \mathbf{p}_2\]로 계산된다.
5. 교점들을 연결하여 등고선을 생성: 각 삼각형에서 얻어진 교점들을 연결하면 하나의 짧은 선분이 만들어지고, 모든 셀을 처리하면 전체 등고선이 완성된다.

CONREC 알고리즘의 장점은 구현이 간단하고 계산이 빠르며, 모호한 경우가 거의 발생하지 않는다는 점이다. 이러한 이유로 과학 계산과 영상 처리 분야에서 직사각형 격자 데이터의 등고선을 생성하는 데 널리 사용되어 왔다.
예제는 그레이 영상의 픽셀값의 등고선을 pseudo coloring(jet(64) of matlab)을 이용해서 그린 것이다.

//=============================================================================
//     CONREC is a contouring subroutine for rectangularily spaced data.
//     It emits calls to a line drawing subroutine supplied by the user
//     which draws a contour map corresponding to real*4data on a randomly
//     spaced rectangular grid. The coordinates emitted are in the same
//     units given in the x() and y() arrays.
//     Any number of contour levels may be specified but they must be
//     in order of increasing value.
//=============================================================================
// function for rounding off the pixels 
int round(double n) { 
    return n > 0 ? int(n+.5): -int(-n+0.5);
} 
struct Ridge {
    double xs, ys;
    double xe, ye;
    double level;
    Ridge() {}
    Ridge(double xs, double ys, double xe, double ye, double level) :
                xs(xs), ys(ys), xe(xe), ye(ye), level(level) {}

};
// Function for line generation 
void DDALine(double x1, double y1, double x2, double y2, 
             int val, CRaster& raster) 
{ 
    double dx = (x2 - x1), dy = (y2 - y1);
    double step;
    if (abs(dx) >= abs(dy))	step = abs(dx);
    else		            step = abs(dy);
    if (step == 0) step = 1;
    dx /= step; dy /= step;
    double x = x1, y = y1;
    int i = 0;

    while (i++ <= step) {
        raster.SetPixel(round(x), round(y), val);
        x += dx; y += dy;
    }
} 
#define xsect(p1, p2) (h[p2] * xh[p1] - h[p1] * xh[p2]) / (h[p2] - h[p1])
#define ysect(p1, p2) (h[p2] * yh[p1] - h[p1] * yh[p2]) / (h[p2] - h[p1])
int conrec(std::vector<int*> &d,  // matrix of data to contour
           CRect rect,            // left <= x < right; top <= y < bottom;
           std::vector<int> &x,   // data matrix column coordinates
           std::vector<int> &y,   // data matrix row coordinates
           std::vector<double> &z,// contour levels in increasing order
           std::vector<Ridge> &ridges)
{
    int		m1, m2, m3, sh[5];
    double  x1, x2, y1, y2;
    double	h[5], xh[5], yh[5];
    const int im[4] = {0, 1, 1, 0}; //curr; right; right-bottom;bottom;
    const int jm[4] = {0, 0, 1, 1};
    const int castab[3][3][3] = {{{ 0, 0, 8}, { 0, 2, 5}, { 7, 6, 9}}, 
                                 {{ 0, 3, 4}, { 1, 3, 1}, { 4, 3, 0}},
                                 {{ 9, 6, 7}, { 5, 2, 0}, { 8, 0, 0}}};

    for (int j = (rect.bottom-2); j >= rect.top; j--) {
        for (int i = rect.left; i <= rect.right-2; i++) {
            double dmin = min(min(d[j][i], d[j+1][i]), min(d[j][i+1], d[j+1][i+1]));
            double dmax = max(max(d[j][i], d[j+1][i]), max(d[j][i+1], d[j+1][i+1]));
            if (dmax < z.front() || dmin > z.back()) continue;
            for (int k = z.size(); k-->0;) {
                if (z[k] < dmin || z[k] > dmax) continue;
                for (int m = 4; m >= 0; m--) {
                    if (m > 0) {
                        h[m]  = d[j+jm[m-1]][i+im[m-1]] - z[k];
                        xh[m] = x[i+im[m-1]];
                        yh[m] = y[j+jm[m-1]];
                    } else {
                        h[0]  = 0.25*(h[1] + h[2] + h[3] + h[4]);
                        xh[0] = 0.5*(x[i] + x[i+1]);
                        yh[0] = 0.5*(y[j] + y[j+1]);
                    }

                    if (h[m] > 0)		 sh[m] =  1; 
                    else if (h[m] < 0)   sh[m] = -1; 
                    else				 sh[m] =  0;
                }
                for (int m = 1; m <= 4; m++) {
                    m1 = m; m2 = 0;
                    if (m != 4) m3 = m + 1;
                    else		m3 = 1;	 
                    int case_value = castab[sh[m1]+1][sh[m2]+1][sh[m3]+1];
                    if (case_value == 0) continue; 
                    switch (case_value) {
                    case 1:
                        x1 = xh[m1]; y1 = yh[m1];
                        x2 = xh[m2]; y2 = yh[m2];
                        break;
                    case 2:
                        x1 = xh[m2]; y1 = yh[m2];
                        x2 = xh[m3]; y2 = yh[m3];
                        break;
                    case 3:
                        x1 = xh[m3]; y1 = yh[m3];
                        x2 = xh[m1]; y2 = yh[m1];
                        break;
                    case 4:
                        x1 = xh[m1]; y1 = yh[m1];
                        x2 = xsect(m2, m3); y2 = ysect(m2, m3);
                        break;
                    case 5:
                        x1 = xh[m2]; y1 = yh[m2];
                        x2 = xsect(m3, m1); y2 = ysect(m3, m1);
                        break;
                    case 6:
                        x1 = xh[m3]; y1 = yh[m3];
                        x2 = xsect(m1, m2); y2 = ysect(m1, m2);
                        break;
                    case 7:
                        x1 = xsect(m1, m2); y1 = ysect(m1, m2);
                        x2 = xsect(m2, m3); y2 = ysect(m2, m3);
                        break;
                    case 8:
                        x1 = xsect(m2, m3); y1 = ysect(m2, m3);
                        x2 = xsect(m3, m1); y2 = ysect(m3, m1);
                        break;
                    case 9:
                        x1 = xsect(m3, m1); y1 = ysect(m3, m1);
                        x2 = xsect(m1, m2); y2 = ysect(m1, m2);
                        break;
                    default:
                        break;
                    }
                    ridges.push_back(Ridge(x1, y1, x2, y2, z[k]));
                }
            }
        }
    }
    return 0;
};
void contourMap(CRaster &raster, int level) {
    CSize sz = raster.GetSize();
    int W = sz.cx, H = sz.cy;
    level = max(2, min(255, level)) ;

    std::vector<int*> data(H);
    std::vector<int> buff(W * H);
    for (int y = 0; y < H; y++) 
        data[y] = &buff[y * W];

    for (int y = 0; y < H; y++) 
        for (int x = 0; x < W; x++) 
            data[y][x] = raster.GetPixel0(x,y);

    // 주어진 데이터의 sampling position;
    std::vector<int> X(W), Y(H);
    std::vector<double> Z(level);
    for (int x = 0; x < W; x++) X[x] = x;
    for (int y = 0; y < H; y++) Y[y] = y;
    // height of levels;
    for (int k = 0; k < level; k++) Z[k] = double(k) * 255.0 / (level-1) ;
    
    //rc: bottom과 right는 포함안됨;
    CRaster levelImg;
    levelImg.SetDimensions(sz, 8);
    std::vector<Ridge> ridges;
    conrec(data, raster.GetRect(), X, Y, Z, ridges);
    //drawing level-line;
    for (int i = ridges.size(); i-->0;) {
        Ridge &r = ridges.at(i);
        int a = int(r.level); 
        a = a < 0 ? 0: a > 255 ? 255: a;
        DDALine(r.xs, r.ys, r.xe, r.ye, a, levelImg);
    }

    // pseudo coloring;
    RGBQUAD *pPal = levelImg.GetPalettePtr();
    for (int i = 0; i < 64; i++)
        for (int k = 0; k < 4; k++) {
            int j = i * 4 + k;
            pPal[j].rgbRed   = colortable[3*i+0];
            pPal[j].rgbGreen = colortable[3*i+1];
            pPal[j].rgbBlue  = colortable[3*i+2];
        }
    //
    raster.SwapBuffer(levelImg);
}