Geometry & Recognition

영상처리에서 Color Quantization는 연속적이거나 매우 넓은 색 영역을 표현하는 이미지를, 시각적 손실을 최소화하면서 제한된 개수의 대표 색상으로 매핑하는 이산화 과정이다. 이를 위해 다양한 알고리즘이 개발되어 왔으면, 여기서는 군집화 품질과 디테일 보존 능력이 뛰어나 학술 및 산업계에서 널리 활용되는 PNN(Pairwise Nearest Neighbor) 컬러 양자화 알고리즘을 구현한다. PNN 알고리즘은 통계학의 상향식 계층적 군집화(Agglomerative Hierarchical Clustering)를 컬러 양자화에 적용한 방법이다. 처음부터 목표 클러스터 수 $K$를 정하고 각 데이터를 $K$개의 클러스터 중 하나에 할당하는 K-means와 달리, PNN은 많은 수의 작은 클러스터에서 출발하여 가장 적절한 두 클러스터를 반복적으로 병합한다. 따라서 데이터 공간의 국소적 구조를 단계적으로 통합하면서 최종적으로 원하는 수의 대표 색상을 얻는다.

알고리즘의 기본 과정은 다음과 같다.

PNN 알고리즘이 임의의 두 클러스터를 병합할 때, 군집 내 오차 제곱합(SSE, Sum of Squared Errors)의 증가량을 최소화하는 방향으로 결정하는 Ward's Criterion을 적용한다. 두 클러스터 $C₁$과 $C₂$를 하나로 병합할 때 발생하는 오차 증가량(병합 비용) $ΔE$는$$\Delta E= \frac{w_1  w_2}{w_1 + w_2 } || \mathbf{\mu}_1 - \mathbf{\mu}_2||^2 $$인데, $\mu_i=(r_i, g_i, b_i)$로 $i$-번째 클러스터의 평균 컬러이고, $w_i$는 그 클러스터에 소속된 픽셀의 수이다. 이 식에서 중요한 점은 단순한 RGB  거리뿐 아니라 각 클러스터에 포함된 픽셀의 수까지 병합 비용에 반영된다는 것이다. 예를 들어 $w_1 \ll w_2$인 경우

$$ \frac{ w_1 w_2}{w_1 + w_2} \approx w_1$$이므로 작은 클러스터와 매우 큰 클러스터의 병합 비용은 대략 작은 클러스터의 크기에 비례한다. 반면 크기가 거의 같은 두 클러스터  $w_1 \approx w_2 \approx w$의 경우에는 $$\frac{w_1 w_2}{w_1+ w_2} \approx \frac{{w}}{2}$$가 된다. 즉, Ward 기준은 색상 중심 사이의 거리와 클러스터의 크기를 동시에 고려하여 병합 여부를 결정한다.

특히 픽셀 수가 적은 색상은 큰 영역을 찾지 하는 색상보다 병합에 따른 전체 SSE 증가에 미치는 영향이 작으므로 비교적 일찍 병합될 수 있다. 그러나 단순히 빈도가 낮다는 이유만으로 가장 큰 배경 클러스터에 흡수되는 것은 아니면, RGB 공간에서 가까운 다른 소규묘 클러스터가 존재한다면 그쪽과 먼저 병합될 수 있다. 이러한 특성 때문에 PNN/Ward 방식은 색상 빈도와 색상 거리를 함께 고려하면서 영상의 색 분포를 점진적으로 압축할 수 있다.

// 실수연산용 RGB color
struct ColorRGB {
    double r, g, b;
    ColorRGB() : r(0), g(0), b(0) {}
    ColorRGB(double rr, double gg, double bb) : r(rr), g(gg), b(bb) {}
};
struct FastCluster {
    ColorRGB	centroid;
    int			weight;
    bool		bActive;
    int			nnIndex;   //nearest neighbor index;
    double		minMergeCost;   
    FastCluster() 
     : weight(0), bActive(false), nnIndex(-1), minMergeCost(DBL_MAX) { }
};
// Disjoint-Set
struct DisjointSet {
    std::vector<int> parent;
    DisjointSet(int n) {
        parent.resize(n);
        for (int i = 0; i < n; ++i)	parent[i] = i; 
    }
    int Find(int i) {
        if (parent[i] == i) return i;
        return parent[i] = Find(parent[i]); // 재귀적으로 루트를 찾으며 경로 압축
    }
    void Union(int dst, int src) {
        dst = Find(dst); src = Find(src);
        if (dst != src) parent[src] = dst;
    }
};

// PNN 병합 비용 계산 함수 (Ward's Method)
inline double calcMergeCost(const FastCluster& c1, const FastCluster& c2) {
    double w1 = (double)c1.weight;
    double w2 = (double)c2.weight;
    double dr = c1.centroid.r - c2.centroid.r;
    double dg = c1.centroid.g - c2.centroid.g;
    double db = c1.centroid.b - c2.centroid.b;
    return (w1 * w2 / (w1 + w2)) * (dr * dr + dg * dg + db * db);
}
// 지정 클러스터의 nearest neighbor를 처음부터 다시 검색
static void updateNearestNeighbor(std::vector<FastCluster>& clusters,
    int index) {
    FastCluster& c = clusters[index];
    c.minMergeCost = DBL_MAX;
    c.nnIndex = -1;

    const int n = (int)clusters.size();
    for (int j = 0; j < n; ++j) {
        if (j == index || !clusters[j].bActive)	continue;
        const double cost = calcMergeCost(c, clusters[j]);
        if (cost < c.minMergeCost) {
            c.minMergeCost = cost;
            c.nnIndex = j;
        }
    }
}
inline int index444(int r, int g, int b) {
    return ((r >> 4) << 8)|((g >> 4) << 4)|((b >> 4));
}
// 알고리즘 main: disjoint-set 사용;
void FastPnnQuantization(const CRaster& raster, int nMaxColors, CRaster& out) {
    if (raster.IsEmpty() || raster.GetBPP() != 24) return;
    CSize sz = raster.GetSize();
    nMaxColors = max(1, nMaxColors);

    // 24비트 색상을 RGB 각각 4비트(총 12비트, 4096가지 색상)로 축소하여
    // 히스토그램 생성(연산속도 증가). 초기 클러스터 개수 압축.
    const int HIST_SIZE = 16 * 16 * 16;
    std::vector<int>		histogram(HIST_SIZE, 0);
    std::vector<ColorRGB>	colorSum(HIST_SIZE, ColorRGB(0,0,0));

    for (int y = 0; y < sz.cy; ++y) {
        BYTE* p = (BYTE*)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; ++x, p += 3) {
            int bin = index444(p[2], p[1], p[0]);
            ++histogram[bin];
            colorSum[bin].r += p[2];
            colorSum[bin].g += p[1];
            colorSum[bin].b += p[0];
        }
    }

    std::vector<FastCluster> clusters;
    clusters.reserve(HIST_SIZE);
    std::vector<int> colorMap(HIST_SIZE, -1);   

    for (int i = histogram.size(); i-->0;) {
        const int count = histogram[i];
        if (count == 0) continue;
        FastCluster c;
        const double inv_count = 1.0 / (double)count;
        c.centroid = ColorRGB(colorSum[i].r * inv_count, 
                              colorSum[i].g * inv_count,
                              colorSum[i].b * inv_count);
        c.weight    = count;
        c.bActive   = true;
        colorMap[i] = clusters.size();
        clusters.push_back(c);
    }
    int nActiveClusters = clusters.size();
    if (nActiveClusters == 0) return;

    // 클러스터가 목표이하면 최종 컬러수 조정;
    if (nMaxColors > nActiveClusters)
        nMaxColors = nActiveClusters;

    // Disjoint-Set 초기화 (초기 고유 색상 수만큼)
    const int nClusters = clusters.size();
    DisjointSet ds(nClusters);
    for (int i = 0; i < nClusters; ++i) {
        clusters[i].minMergeCost = DBL_MAX;
        clusters[i].nnIndex = -1;
    }
    // min(merge cost)cluster 찾기;
    for (int i = 0; i < nClusters; ++i) {
        for (int j = 0; j < i; ++j) {
            double cost = calcMergeCost(clusters[i], clusters[j]);
            if (cost < clusters[i].minMergeCost) {
                clusters[i].minMergeCost = cost;
                clusters[i].nnIndex = j;
            }
            if (cost < clusters[j].minMergeCost) {
                clusters[j].minMergeCost = cost;
                clusters[j].nnIndex = i;
            }
        }
    }

    // Ward merging
    while (nActiveClusters > nMaxColors) {
        double best_cost = DBL_MAX;
        int merge_i = -1;
        // 각 cluster가 가지고 있는 NN 중 global minimum 검색
        for (int i = 0; i < nClusters; ++i) {
            if (!clusters[i].bActive) 	 continue;
            if (clusters[i].nnIndex < 0) continue;

            if (clusters[i].minMergeCost < best_cost) {
                best_cost = clusters[i].minMergeCost;
                merge_i = i;
            }
        }
        if (merge_i < 0) break;
        int merge_j = clusters[merge_i].nnIndex;

        if (merge_j < 0 || merge_j >= nClusters || !clusters[merge_j].bActive) {
            updateNearestNeighbor(clusters, merge_i);
            if (clusters[merge_i].nnIndex < 0) break;
            continue;
        }
        // merge_j -> merge_i
        FastCluster& c1 = clusters[merge_i];
        FastCluster& c2 = clusters[merge_j];
        int w1 = c1.weight;
        int w2 = c2.weight;

        double inv_weight = 1.0 / double(w1 + w2);
        c1.centroid.r = (c1.centroid.r * w1 + c2.centroid.r * w2) * inv_weight;
        c1.centroid.g =	(c1.centroid.g * w1 + c2.centroid.g * w2) * inv_weight;
        c1.centroid.b = (c1.centroid.b * w1 + c2.centroid.b * w2) * inv_weight;
        c1.weight     = w1 + w2;
        c2.bActive    = false;
        ds.Union(merge_i, merge_j);
        --nActiveClusters;

        // nearest-neighbor 정보 갱신
        // merge_i의 centroid가 변경됨: 자신의 NN은 처음부터 다시 계산
        c1.minMergeCost = DBL_MAX;
        c1.nnIndex = -1;
        for (int k = 0; k < nClusters; ++k) {
            if (k == merge_i || !clusters[k].bActive) continue;
            // 기존 NN이 변경/삭제된 cluster였다면 전체 NN 재검색 필요
            if (clusters[k].nnIndex == merge_i || clusters[k].nnIndex == merge_j)
                updateNearestNeighbor(clusters, k);
            // 아니면 기존 NN은 여전히 유효 새로 위치가 변한 merge_i만 추가 비교
            else {
                double cost = calcMergeCost(clusters[k], c1);
                if (cost < clusters[k].minMergeCost) {
                    clusters[k].minMergeCost = cost;
                    clusters[k].nnIndex = merge_i;
                }
            }
            // 동시에 merge_i 자신의 NN 검색
            double cost = calcMergeCost(c1, clusters[k]);
            if (cost < c1.minMergeCost) {
                c1.minMergeCost = cost;
                c1.nnIndex = k;
            }
        }
    }

    // 각 클러스의 대표컬러를 담은 palette 생성;
    std::vector<RGBQUAD> palette;
    palette.reserve(nMaxColors);
    std::vector<int> clusterMap(clusters.size(), -1);

    for (int cid = clusters.size(), pid = 0; cid-->0;) {
        if (!clusters[cid].bActive) continue;
        const ColorRGB& color = clusters[cid].centroid;
        clusterMap[cid] = pid++;
        RGBQUAD q = {0,0,0,0};
        q.rgbRed   = max(0, min(255, int(color.r + 0.5)));
        q.rgbGreen = max(0, min(255, int(color.g + 0.5)));
        q.rgbBlue  = max(0, min(255, int(color.b + 0.5)));
        palette.push_back(q);
    }

    // Disjoint-Set의 Find를 사용하여 컬러 단위 최종 인덱스 맵 생성
    // color444 -> colorMap -> cluster_root -> palette_index;
    std::vector<int> inverseMap(colorMap.size());
    for (int color = colorMap.size(); color-->0; ) {
        int clusterId  = colorMap[color];
        if (clusterId < 0) continue;    //음수 = 한번도 안나타난 컬러;
        inverseMap[color] = clusterMap[ds.Find(clusterId)];
    }
    // 출력용 24-bit 영상 생성;
    out.SetDimensions(sz, 24);
    for (int y = 0, i = 0; y < sz.cy; ++y) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)out.GetLinePtr(y);
        for (int x = 0; x < sz.cx; ++x, p += 3) {
            int idx = inverseMap[index444(p[2], p[1], p[0])];
            *q++ = palette[idx].rgbBlue;
            *q++ = palette[idx].rgbGreen;
            *q++ = palette[idx].rgbRed;;
        }
    }
}