컬러를 양자화하는 방법으로 median-cut 알고리즘 이외에도 Octree을 이용한 알고리즘도 많이 쓰인다. (페인트샵 프로를 보면 이 두 가지 방법과 더불어 간단한 websafe 양자화 방법을 이용한다). Octree는 8개의 서브 노드를 가지는 데이터 구조이다. 이는 R, G, B의 비트 평면에서의 비트 값(각각의 레벨에서 8가지가 나옴)을 서브 노드에 할당하기 편리한 구조이다. 풀컬러를 이 Octree를 이용해서 표현을 하면 깊이가 9인 트리가 생성이 된다 (root + 각각의 비트 평면 레벨 = 8^8 = num of true colors). Octree quantization은 이 트리의 깊이를 조절하여서 leaf의 수가 원하는 컬러 수가 되도록 조절하는 알고리즘이다. 전체 이미지를 스캔하면서 트리를 동적으로 구성하고 난 후 트리의 leaf의 수가 원하는 컬러수보다 작으면 각각의 leaf가 표현하는 컬러를 팔레트로 만들면 된다. 그러나 컬러를 삽입하면서 만들어지는 leaf의 수가 원하는 컬러 개수보다 많아지는 경우는 가장 깊은 leaf 노드를 부모 노드에 병합시켜서 leaf 노드의 수를 줄이는 과정을 시행한다. 이 병합 과정은 원래 leaf 노드들의 컬러 평균을 취하여 부모 노드에 전달한다. 따라서 Octree 알고리즘은 RGB 컬러에서 상위 비트가 컬러에 대표성을 가지도록 하위 비트 값을 병합하여서 컬러 수를 줄이는 방법을 이용한 것이다. 특정 RGB 비트 평면에서 R, G, B비트 값을 이용하여서 child-노드의 인덱스를 만드는 방법은
child-node index = (R-비트<<2)|(G-비트<<1)|(B-비트) ;
(R, G, B)=(109,204,170)의 경우;
node | level | 0 1 2 3 4 5 6 7 ------------------------- R | 0 1 1 0 1 1 0 1 G | 1 1 0 0 1 1 0 0 B | 1 0 1 0 1 0 1 0 ------|------------------ child | 3 6 5 0 7 6 1 4 index | Octree 알고리즘은 컬러 이미지를 스캔하면서 트리를 동적으로 구성할 수 있으므로 median-cut 알고리즘보다도 메모리 사용에 더 효과적인 방법을 제공한다. 그리고 트리 깊이를 제한하여 일정 깊이 이상이 안되도록 만들 수 있으므로 빠른 알고리즘을 구현할 수 있다. 최종적으로 팔레트를 구성하는 컬러 값은 leaf 노드에 축적된 컬러의 평균값이 된다 (일반적으로 leaf-노드가 나타내는 비트 값과는 불일치가 생긴다). 이 알고리즘은 재귀적인 방법을 이용하면 쉽게 구현이 가능하고, 또한 웹상에서 구현된 소스도 찾을 수 있다.
참고 논문: "A Simple Method for Color Quantization: Octree Quantization." by M. Gervautz and W. Purgathofer 1988.
256 컬러 이미지(트리 깊이 = 5): RGB 값이 주어지면 전체 트리에서 해당 leaf을 찾아서 팔레트 인덱스를 얻는다. 이와는 다른 방법으로는 팔레트가 주어졌으므로 주어진 RGB 값과 가장 가까운 유클리디안 거리를 주는 팔레트의 인덱스를 할당하는 방법도 있다. 이 방법을 이용하면 아래의 결과와 약간 차이가 생긴다.
컬러 영상을 처리할 때 가장 흔히 사용하는 컬러 표현은 RGB 컬러이다. 이것은 R, G, B에 각각 8-비트를 할당하여 256-단계를 표현할 수 있게 하여, 전체적으로 256x256x256=16777216가지의 컬러를 표현할 수 있게 하는 것이다. 그러나 인간의 눈은 이렇게 많은 컬러를 다 구별할 수 없으므로 24-비트 RGB 컬러를 사용하는 경우는 대부분의 경우에 메모리의 낭비와 연산에 오버헤드를 가져오는 경우가 많이 생긴다. RGB 컬러 영상을 R, G, B를 각각 한축으로 하는 3차원의 컬러 공간에서의 벡터(점)로 표현이 가능하다. 컬러 영상의 픽셀들이 RGB삼차원의 공간에 골고루 펴져 있는 경우는 매우 드물고, 많은 경우에는 이 컬러 공간에서 군집(groups)을 이루면서 분포하게 된다. 하나의 군(group)은 유사한 RGB 값을 갖는 픽셀들로 구성이 되므로, 이 군에 포함이 되는 픽셀들에게 (군의 중앙에 해당하는) 대표적인 컬러 값을 대체하면 그 군에 포함이 된 픽셀은 이젠 RGB공간에서 한 점으로 표현이 되고, RGB공간상에는 픽셀 수만큼의 점이 있는 것이 아니라, 대표 RGB 값에 해당하는 점만이 존재하게 된다. 따라서 적당한 Lookup테이블(colormap)을 이용하면, 적은 메모리 공간만을 가지고도 원본의 컬러와 유사한 컬러를 구현할 수 있다.
이제 문제는 원래의 컬러 공간을 어떻게 군집화 하는가에 대한 것으로 바뀌었다. 간단한 방법으로는 원래의 컬러 영상이 차지는 하는 RGB공간에서의 영역을 감싸는 최소의 박스를 찾고, 이것을 원하는 최종적으로 원하는 컬러수만큼의 박스로 분할을 하는 것이다. 그러나 박스를 어떨게 분할을 해야만 제대로 컬러를 나누고, 또한 효율적으로 할 수 있는가를 고려해야 한다. 분할된 박스의 크기가 너무 크면 제대로 된 대푯값을 부여하기가 힘들어지고, 너무 작게 만들면 원하는 수에 맞추기가 어렵다.
Median-Cut 양자화(quantization)에서 사용하는 방법은 박스의 가장 긴축을 기준으로 그 축으로 projection 된 컬러 히스토그램의 메디안 값을 기준으로 분할을 하여서 근사적으로 픽셀들을 절반 정도 되게 분리를 한다 (한축에 대한 메디안이므로 정확히 반으로 분리되지 않는다). 두 개의 박스가 이렇게 해서 새로 생기는데, 다시 가장 많은 픽셀을 포함하는 박스를 다시 위의 과정을 통해서 분할을 하게 된다. 이렇게 원하는 수의 박스를 얻을 때까지 위의 과정을 반복적으로 시행을 하게 된다.
여기서 원래의 컬러 값을 모두 이용하게 되면 계산에 필요한 히스토그램을 만들기 위해서 너무 많은 메모리를 사용하게 되고 이것이 또한 연산의 오버헤드로 작용하게 되므로 RGB 컬러 비트에서 적당히 하위 비트를 버리면, 초기의 RGB공간에서의 histogram의 크기를 줄일 수 있게 된다.(보통 하위 3-비트를 버려서, 각각 5-비트만 이용하여, 전체 컬러의 개수를 32x32x32= 32768로 줄인다)
이렇게 RGB공간에서의 컬러 분포가 몇 개의 대표적인 컬러(예:박스의 중앙값)로 줄어들면(양자화 과정:: 공간에 연속적인 분포가 몇 개의 점으로 대체됨), 원본 영상의 각각의 픽셀에서의 대체 컬러 값은 원래의 컬러와 가장 유사한, 즉 RGB 공간에서 Euclidean 거리가 가장 작은 박스의 컬러로 대체하면 된다.
그러나 너무 적은 수의 컬러로 줄이게 되면 인접 픽셀 간의 컬러 값의 차이가 눈에 띄게 나타나는 현상이 생기게 된다. 이러한 것을 줄이기 위해서는 디더링(dithering) 과정과 같은 후처리가 필요하게 된다.
int MedCutQuantizer(CRaster& raster, CRaster& out8) {
const int MAX_CUBES = 256;
const int HSIZE =32 * 32 * 32; // # of 15-bit colors;
int hist[HSIZE] = {0};
CSize sz = raster.GetSize();
for (int y = 0; y < sz.cy; y++) {
BYTE *p = (BYTE *)raster.GetLinePtr(y);
for (int x = 0; x < sz.cx; x++) {
int b = *p++; int g = *p++; int r = *p++;
hist[bgr15(b, g, r)]++;
}
}
std::vector<int> hist_ptr;
for (int color15 = 0; color15 < HSIZE; color15++)
if (hist[color15] != 0) {
// 0이 아닌 히스토그램의 bin에 해당하는 15비트 컬러를 줌;
hist_ptr.push_back(color15);
}
std::vector<Cube> cubeList;
Cube cube(0, hist_ptr.size()-1, sz.cx * sz.cy, 0);
shrinkCube(cube, &hist_ptr[0]);
cubeList.push_back(cube);
while (cubeList.size() < MAX_CUBES) {
int level = 255, splitpos = -1;
for (int k = cubeList.size(); k-->0;) {
if (cubeList[k].lower == cubeList[k].upper) ; // single color; cannot be split
else if (cubeList[k].level < level) {
level = cubeList[k].level;
splitpos = k;
}
}
if (splitpos == -1) // no more cubes to split
break;
// Find longest dimension of this cube
Cube cube = cubeList[splitpos];
int dr = cube.rmax - cube.rmin;
int dg = cube.gmax - cube.gmin;
int db = cube.bmax - cube.bmin;
int longdim = 0;
if (dr >= dg && dr >= db) longdim = 0;
if (dg >= dr && dg >= db) longdim = 1;
if (db >= dr && db >= dg) longdim = 2;
// 가장 넓게 분포한 color를 colorTab의 15비트 컬러의 상위비트 이동;
reorderColors(&hist_ptr[0], cube.lower, cube.upper, longdim);
// hist_ptr값을 증가 순서로 바꿈->가장 폭이 넓은 컬러 방향으로 정렬됨;
quickSort(&hist_ptr[0], cube.lower, cube.upper);
// hist_ptr을 다시 원컬러로 복구;
restoreColorOrder(&hist_ptr[0], cube.lower, cube.upper, longdim);
// Find median
int count = 0;
int med = 0;
for (med = cube.lower; med < cube.upper; med++) {
if (count >= cube.count/2) break;
count += hist[hist_ptr[med]];
}
// cube에 들어있는 color의 median을 기준으로 두개 cube로 쪼김;
// 낮은 쪽은 원래 cube를 대체하고 높은 쪽은 새로 list에 추가;
Cube cubeLow(cube.lower, med - 1, count, cube.level + 1);
shrinkCube(cubeLow, &hist_ptr[0]);
cubeList[splitpos] = cubeLow; // add in old slot
//
Cube cubeHigh(med, cube.upper, cube.count - count, cube.level + 1);
shrinkCube(cubeHigh, &hist_ptr[0]);
cubeList.push_back(cubeHigh);
}
// 각 cube에 포함된 컬러를 하나의 대표값으로 설정;
// 대표값은 cube내 컬러의 평균이다;
BYTE rLUT[256], gLUT[256], bLUT[256];
for (int k = cubeList.size(); k-->0;) {
int rsum = 0, gsum = 0, bsum = 0;
Cube &cube = cubeList[k];
for (int i = cube.lower; i <= cube.upper; i++) {
int color15 = hist_ptr[i];
rsum += getRed(color15) * hist[color15];
gsum += getGreen(color15) * hist[color15];
bsum += getBlue(color15) * hist[color15];
}
rLUT[k] = int(rsum / (double)cube.count + .5);
gLUT[k] = int(gsum / (double)cube.count + .5);
bLUT[k] = int(bsum / (double)cube.count + .5);
}
// 각 cube에 포함된 color를 대표하는 컬러값을 쉽게 찾도록 lookup table
// 설정함;
BYTE inverseMap[HSIZE];
for (int k = cubeList.size(); k-->0;) {
// 각 cube 내의 컬러는 하나의 대표컬러로...;
Cube& cube = cubeList[k];
for (int i = cube.lower; i <= cube.upper; i++)
inverseMap[hist_ptr[i]] = (BYTE)k;
}
// 8-bit indexed color 이미지 생성;
out8.SetDimensions(sz, 8);
for (int y = 0; y < sz.cy; y++) {
BYTE *p = (BYTE *)raster.GetLinePtr(y);
BYTE *q = (BYTE *)out8.GetLinePtr(y);
for (int x = 0; x < sz.cx; x++) {
int b = *p++; int g = *p++; int r = *p++;
*q++ = inverseMap[bgr15(b, g, r)];
}
}
// 8비트 이미지 팔레트 설정;
for (int i = 0; i < 256; i++) {
out8.m_palette[i].rgbBlue = bLUT[i];
out8.m_palette[i].rgbGreen = gLUT[i];
out8.m_palette[i].rgbRed = rLUT[i];
}
// SaveRaster(out8, "res_medcut.bmp");
return 1;
}
같은 평균 픽셀값을 가지고 있더라도 4-방향으로 서로 연결된 영역이 아니면 합병되지 못하고 서로 다른 영역으로 남는다. 이를 하나의 영역으로 만들려면 분할된 영역을 다시 검사를 하는 과정이 필요하다. 자세한 과정은 Statistical Region Merge(SRM) 알고리즘(kipl.tistory.com/98)에 잘 구현이 되어 있다.
struct Universe {
int n;
std::vector<int> sizes;
std::vector<int> sums;
std::vector<int> parent;
Universe(int n, BYTE *image) {
this->n = n;
parent.resize(n, -1); // all nodes are root;
sizes.resize(n, 1);
sums.resize(n);
for (int i = 0; i < n; i++)
sums[i] = image[i];
}
int FindCompress (int node) { // find root node;
if ( parent[node] < 0 ) return node; // root-node;
else return parent[node] = FindCompress (parent[node]);
// 찾음과 동시에 depth을 줄인다;
}
int Find(int node) { // find root;
while (parent[node] >= 0) node = parent[node];
return node;
}
bool Predicate(int root1, int root2, int thresh) {
double diff = double(sums[root2]) / sizes[root2]
- double(sums[root1]) / sizes[root1];
return fabs(diff) < thresh;
}
void Merge(int root1, int root2) {
sums[root2] += sums[root1];
sizes[root2] += sizes[root1];
parent[root1] = root2; // r2 트리로 합병함;
}
};
// root 노드는 분할영역의 픽셀 갯수와 픽셀 값의 합을 저장한다.
// 처음 각 픽셀이 root 이고, 픽셀 갯수는 1, 픽셀값은 자신의 픽셀값을 저장;
BOOL UnionFindAverage(BYTE *image, int width, int height, int thresh) {
Universe *UF = new Universe(width * height, image);
// 4-connectivity:
// 첫행; LEFT 픽셀이 속한 영역의 평균값과 차이가 thresh 내이면 left로 합병;
for ( int x = 1; x < width; x++ ) {
int left = UF->FindCompress (x - 1);
int curr = UF->FindCompress (x);
if (UF->Predicate(curr, left, thresh))
UF->Merge(curr, left);
}
//Flatten(y=0);
for (int x = 0; x < width; x++) UF->FindCompress(x);
for ( int y = 1, pos = y * width; y < height; y++ ) {
// x = 0; TOP 픽셀이 속학 영역과 합병 체크;
int up = UF->FindCompress (pos - width);
int curr = UF->FindCompress (pos);
if (UF->Predicate(curr, up, thresh))
UF->Merge(curr, up);
pos++;
// TOP-LEFT 픽셀 영역과 합병 체크;
for ( int x = 1; x < width; x++, pos++ ) {
int left = UF->FindCompress(pos - 1);
int curr = UF->FindCompress(pos);
// left와 합병이 되면 left가 up과 합병이 되는지 확인;
if (UF->Predicate(curr, left, thresh)) {
UF->Merge(curr, left);
curr = left;
}
int up = UF->FindCompress (pos - width);
if ((curr != up) && (UF->Predicate(curr, up, thresh)))
UF->Merge(curr, up);
}
// Flatten(y>0)
for (int x = 0; x < width; x++) UF->FindCompress(x + y * width);
}
// 평균 이미지 생성;
for ( int k = width*height; k-->0;) {
int root = UF->Find(k);
int avg = int(double(UF->sums[root]) / UF->sizes[root] + 0.5);
image[k] = avg > 255 ? 255 : avg;
}
delete UF;
return TRUE;
};
Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.
이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram($\tt ch$)과 평균을 계산할 수 있는 1차 culmuative histogram($\tt cxh$)을 입력으로 사용한다.
책상에 세워둔 연필이 넘어져서 바닥에 닿는데 걸리는 시간은 어떻게 될까? 대략 경험적으로 보면 연필 길이가 길수록 넘어지는데 시간이 오래 걸린다. 왜 그럴까? 넘어지는데 걸리는 시간은 연필의 운동을 결정하는 물리량과 관련 있는데, 우선 크기와 관련 있는 길이$(L)$가 있고, 넘어지려면 토크를 작용해야 하는데 이는 중력(가속도: $g$)과 관련이 있다.(마찰/수직 항력은 회전축에 걸리므로 무관하다). 물론 질량에도 의존할 수 있는데, 중력 가속도, 길이, 질량의 물리량으로 시간의 단위를 만들 수 있는 조합은 $\text{const} \times\sqrt { L /g}$ 밖에 없다. 즉, 넘어지는데 걸리는 시간은 길이에 제곱근에 비례한다.
구체적으로 얼마나 걸리는지 계산을 시도해 보자. 연필이 넘어지면 수직과 이루는 각 $\theta$가 증가한다: $0\rightarrow \pi/2$. 뉴턴 방정식을 이용하면 $\theta$가 만족해야 할 미분방정식을 얻을 수 있지만, 실질적으로 일하는 힘이 중력밖에 없으므로 역학적 에너지가 보존된다는 사실을 이용하는 것이 좀 더 수월하다. 연필을 균일한 막대로 근사하면 넘어지는 과정은 연필심에 대한 회전운동이다. 역학적 에너지 보존식을 쓰면
기대(?)와는 다르지만 이 결과는 구체적으로 계산하지 않더라도 예상할 수 있다. 왜냐면 완전히 똑바로 서 있으면 회전의 시작에 필요한 토크를 생성하는 힘이 없기 때문이다. 연필에는 중력이 작용하고 끝에 마찰력이나 수직 항력이 있긴 하지만 토크를 만들지는 못한다. 따라서 연필을 넘어뜨리기 위해서는 처음에 약간의 충격을 주던지(속도 제공) 아니면 약간 기울인 상태에서 시작해야 한다. 처음 $\theta_0$의 각도에서 시작하였다면 걸리는 시간은 위 적분식에서
