Marching Squares

2차원 데이터 격자에서 등고선을 그리는 Marching Squares 알고리즘:

•  2차원 배열 데이터를 2x2 크기의 작은 사각형(셀)으로 분할한다.
• 각 셀의 꼭짓점에서 기준값보다 높은지 또는 낮은지에 따라 총 16가지 교차 패턴으로 분류한다. 
• 셀의 변과 레벨선의 교차점은 변의 두 꼭짓점에서의 값을 선형 보간해서 찾은 후 16가지 교차 패턴에 해당하는 등고선 조각을 만든다.  셀의 한 변의 꼭짓점이 ${\bf p}_A$와 ${\bf p}_B$이고, 꼭짓점에서 데이터 값이 각각 $v_A$, $v_B$인 경우, 등고선 $L$과 만나는 지점 ${\bf p}$는 다음과 같이 정의된다:
  $$ {\bf p} = {\bf p}_A + \mu ( {\bf p}_B - {\bf p}_A), \quad \mu = \frac{L - v_A}{v_B - v_A}$$ $v_A=v_B$인 경우 그 변이 등고선의 일부이므로 교차점이 정의가 되지 않는 상황이지만, 여기서는 두 꼭짓점의 중점을 교차점으로 잡는 convention을 취한다.
• 교차패턴:  셀의 꼭지점을 (left, bottom)=0, (right, bottom)=1, (right, top)=2, (left,top)=3으로 할 때
  

pattern	기준값 이상인 꼭짓점	교차 변
0, 15	없음 / 전부	없음
1, 14	0 / 1,2,3	bottom–left
2, 13	1 / 0,2,3	bottom–right
3, 12	0,1 / 2,3	left–right
4, 11	2 / 0,1,3	right–top
5	0,2	ambiguous
6, 9	1,2 / 0,3	bottom–top
7, 8	0,1,2 / 3	left–top
10	1,3	ambiguous

패턴 5와 10에서는 셀 내부에 saddle point가 존재하며, saddle point에서의 값에 따라 등고선의 연결 방식이 달라진다. 간단하게는 네 꼭짓점 값의 평균을 셀 중심에서의 값으로 보고 연결 방향을 결정할 수도 있다. 그러나 이러한 판별이 항상 올바른 등고선 topology를 주는 것은 아니다.

등고선과 셀의 변이 만나는 위치를 구할 때 선형보간을 사용하므로, 셀 내부의 값도 이와 일관되게 양선형보간(bilinear interpolation)으로 표현하는 것이 자연스럽다. 부호 판별을 간단히 하기 위해 각 꼭짓점의 값에서 등고선의 높이 $L$을 뺀 값을 사용하자. 셀의 좌하단 꼭짓점에서 반시계 방향으로 각 꼭짓점의 값을 $v_0,v_1,v_2,v_3$라 하고,$$q_i=v_i-L,\qquad i=0,1,2,3$$

로 정의한다. 그러면 각 $q_i$의 부호는 해당 꼭짓점이 등고선보다 높은지 낮은지를 나타낸다. 셀 내부의 양선형보간 함수를$$g(x,y)=a+bx+cy+dxy$$로 놓으면,$$g(0,0)=q_0,\qquad g(1,0)=q_1,\qquad g(1,1)=q_2,\qquad g(0,1)=q_3$$이므로 계수 $a,b,c,d$가 결정되어$$g(x,y) =q_0+(q_1-q_0)x+(q_3-q_0)y +(q_0-q_1-q_3+q_2)xy$$를 얻는다. Saddle point는$$\frac{\partial g}{\partial x} = \frac{\partial g}{\partial y} =0$$을 만족하는 점이며, 이 점에서 $g$의 값은$$g_s= \frac{q_0q_2-q_1q_3} {q_0-q_1-q_3+q_2}$$임을 보일 수 있다.

패턴 5에서는

$$q_0>0,\qquad q_2>0,\qquad q_1<0,\qquad q_3<0$$이므로 $g_s$의 분모 $q_0-q_1-q_3+q_2$는 항상 양수이다. 따라서 $g_s$의 부호는 분자 $D=q_0q_2-q_1q_3$의 부호만으로 결정할 수 있다. $g_s>0\,$이면 등고선은 (left, top)과 (bottom, right)을 각각 연결하고, $g_s<0\,$이면 (left, bottom)과 (top, right)을 각각 연결한다.

패턴 10에서는 $$q_0<0,\qquad q_2<0,\qquad q_1>0,\qquad q_3>0$$이므로 위와 반대로 $g_s$의 분모가 항상 음수가 된다. 따라서 분자 $D=q_0q_2-q_1q_3$의 부호와 $g_s$의 부호가 반대가 되며, 이에 따라 패턴 5와 반대의 방식으로 등고선의 연결 방향을 결정한다.

아래 코드에서는  saddle point의 값을 이용하는 asymptotic decider를 적용하였다. $D=0\,$이면 등고선이 saddle point를 정확히 지나므로 connectivity가 퇴화한다. 구현에서는 일관된 tie-breaking convention을 사용한다.

// 등고선 선분(시작점과 끝점)을 나타내는 구조체
struct Segment {
    CfPt p1, p2;
    double level;
    Segment() {}
    Segment(const CfPt& a, const CfPt& b, double z)
        : p1(a), p2(b), level(z) {}
};

// 두 점 사이에서 목표 높이(isoValue)가 위치하는 지점을 선형 보간하는 함수
CfPt interpolate(CfPt pA, CfPt pB, double valA, double valB, double isoValue) {
    if (fabs(valA - valB) < 1e-6) {
        return CfPt((pA.x + pB.x) / 2.0, (pA.y + pB.y) / 2.0);
    } else {
        double mu = (isoValue - valA) / (valB - valA);
        return CfPt(pA.x + mu * (pB.x - pA.x), pA.y + mu * (pB.y - pA.y));
    }
}

// Marching Squares 알고리즘을 이용해 등고선 세그먼트를 생성하는 함수
std::vector<Segment>
extractContourSegments(const std::vector<int*>& grid, int cols, int rows,
            const std::vector<double>& isoValues) {
    std::vector<Segment> segPool;
    if (rows < 2 || cols < 2) return segPool;

    // 격자의 각 셀(Square)을 순회
    for (int r = 0; r < rows - 1; ++r) {
        for (int c = 0; c < cols - 1; ++c) {
            // 네 꼭짓점의 값
            double v0 = grid[r][c];         // 좌하
            double v1 = grid[r][c + 1];     // 우하
            double v2 = grid[r + 1][c + 1]; // 우상
            double v3 = grid[r + 1][c];     // 좌상
            
            double vmin = min(min(v0, v1), min(v2, v3));
            double vmax = max(max(v0, v1), max(v2, v3));
            if (vmin == vmax) continue;
            
            // check vmin < v <= vmax;
            int k2 = isoValues.size() - 1;
            while (k2 > 0 && isoValues[k2] > vmax) k2--;
            int k1 = 0;
            while (k1 < k2 && isoValues[k1] <= vmin) k1++;  
            
            // 각 꼭짓점의 좌표 (단위 격자 크기를 1)
            CfPt p0 = CfPt(c,     r);
            CfPt p1 = CfPt(c + 1, r);
            CfPt p2 = CfPt(c + 1, r + 1);
            CfPt p3 = CfPt(c,     r + 1);
            for (int k = k1; k <= k2; k++) {
                const double isoValue = isoValues[k];
                double q0 = v0 - isoValue;
                double q1 = v1 - isoValue;
                double q2 = v2 - isoValue;
                double q3 = v3 - isoValue;
                double decider = q0 * q2 - q1 * q3;
                
                // 각 꼭짓점에서 값과 기준선 값 비교 --> 코드 생성
                int pattern = 0;
                if (q0 >= 0.0) pattern |= 1;
                if (q1 >= 0.0) pattern |= 2;
                if (q2 >= 0.0) pattern |= 4;
                if (q3 >= 0.0) pattern |= 8;
                
                // 변 위의 보간점들 정의
                CfPt bottom = interpolate(p0, p1, v0, v1, isoValue);
                CfPt right  = interpolate(p1, p2, v1, v2, isoValue);
                CfPt top    = interpolate(p3, p2, v3, v2, isoValue);
                CfPt left   = interpolate(p0, p3, v0, v3, isoValue);

                // 셀의 변과의 교차 패턴에 따른 등고선 조각 선택
                switch (pattern) {
                case 0: case 15:
                    // 모든 점이 위이거나 아래이므로 선 없음
                    break;
                case 1: case 14:
                    segPool.push_back(Segment(bottom, left, isoValue));
                    break;
                case 2: case 13:
                    segPool.push_back(Segment(bottom, right, isoValue));
                    break;
                case 3: case 12:
                    segPool.push_back(Segment(left, right, isoValue));
                    break;
                case 4: case 11:
                    segPool.push_back(Segment(top, right, isoValue));
                    break;
                case 6: case 9:
                    segPool.push_back(Segment(bottom, top, isoValue));
                    break;
                case 7: case 8:
                    segPool.push_back(Segment(left, top, isoValue));
                    break;

                case 5:
                    // saddle point(모호) -> saddle point -> asymptotic decider;
                    if (decider >= 0.0) {
                        // saddle pt가 높음: 0번과 2번의 high 영역이 서로 연결
                        segPool.push_back(Segment(bottom, right, isoValue));
                        segPool.push_back(Segment(top, left, isoValue));
                    } else {
                        // saddle pt가 낮음: 0번과 2번은 서로 분리
                        segPool.push_back(Segment(bottom, left, isoValue));
                        segPool.push_back(Segment(top, right, isoValue));
                    }
                    break;
                case 10:
                    // saddle point -> asymptotic decider;
                    if (decider <= 0.0) {
                        // saddle pt 높음: 1번과 3번의 high 영역이 서로 연결
                        segPool.push_back(Segment(bottom, left, isoValue));
                        segPool.push_back(Segment(top, right, isoValue));
                    } else {
                        // saddle pt가 낮음: 1번과 3번은 서로 분리
                        segPool.push_back(Segment(bottom, right, isoValue));
                        segPool.push_back(Segment(top, left, isoValue));
                    }
                    break;
                }
            }
        }
    }

    return segPool;
}
void drawContourMap(const CRaster& raster, int nLevels, CRaster& levelImg) {
    if (raster.IsEmpty() || raster.GetBPP() != 8) return;
    CSize sz = raster.GetSize();
    nLevels = max(2, min(255, nLevels));
    
    std::vector<int*> data(sz.cy);
    std::vector<int> buff(sz.cx * sz.cy);
    for (int y = 0; y < sz.cy; y++) {
        data[y] = &buff[y * sz.cx];
        BYTE *q = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++)
            data[y][x] = *q++;
    }
    
    // height of levels
    std::vector<double> zLevel(nLevels);

    for (int k = 0; k < nLevels; k++) 
        zLevel[k] = 255.0 * k / (nLevels - 1);

    std::vector<Segment> segPool = extractContourSegments(data, sz.cx, sz.cy, zLevel);

    // drawing level-line
    levelImg.SetDimensions(sz, 8);

    for (int i = segPool.size(); i-- > 0;)
        DDALine(segPool[i], levelImg);

    // pseudo coloring(64 colors);
    falseColoring(levelImg, 64);
}