등고선을 그리는 subroutine

Corec는 2차원 곡면의 등고선을 그려주는 procedure이다. 알고리즘에 대한 구체적인 정보는 https://paulbourke.net/papers/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);
}