n 차 Bezier 곡선은 두 개의 (n - 1) 차 Bezier 곡선의 선형보간으로 표현할 수 있다. Bezier 곡선은 Bernstein 다항식을 이용해서도 표현할 수도 있지만, 높은 찻수의 곡선일 때는 De Casteljau's Algorithm을 이용하는 것이 수치적으로 보다 안정적인 결과를 준다.
// De Casteljau's algorithm; recursive version; slow for larger deg;
double Bezier(int deg, double Q[], double t) {
if (deg == 0) return Q[0];
else if (deg == 1) return (1-t)*Q[0] + t*Q[1];
else if (deg == 2) return (1-t)*((1-t)*Q[0] + t*Q[1]) + t*((1-t)*Q[1] + t*Q[2]);
else return (1 - t) * Bezier(deg-1, &Q[0], t) + t * Bezier(deg-1, &Q[1], t);
}
// De Casteljau's algorithm(degree=n-1);
// non-recursive. Bezier() modifies Q's;
double Bezier(int deg, double Q[], double t) {
if (deg==0) return Q[0];
else if (deg==1) return (1-t)*Q[0] + t*Q[1];
else if (deg==2) return (1-t)*((1-t)*Q[0] + t*Q[1]) + t*((1-t)*Q[1] + t*Q[2]);
for (int k = 0; k < deg; k++)
for (int j = 0; j < (deg - k); j++)
Q[j] = (1 - t) * Q[j] + t * Q[j + 1];
return Q[0];
}
void BezierCurve(std::vector<CfPt> &cntls,
int segments, std::vector<CfPt> &curves) {
std::vector<double> xp(cntls.size()), yp(cntls.size());
curves.resize(segments + 1);
for (int i = 0; i <= segments; ++i) {
double t = double(i) / segments;
// clone control points; non-rec version modifies inputs;
for (int k = cntls.size(); k-->0;) {
xp[k] = cntls[k].x; yp[k] = cntls[k].y;
}
curves[i] = CfPt(Bezier(xp.size()-1, &xp[0], t), Bezier(yp.size()-1, &yp[0], t));
}
}
평면 위에서 주어진 점집합을 포함하는 가장 작은 원을 찾는 문제는 가장 단순하게는 임의의 두 점이 지름이 되는 원과 임의의 세 점이 만드는 외접원을 모두 조사하여 찾으면 된다. 따라서 O(n^4)의 복잡도를 가질 것이다. 그러나 이 문제는 점집합의 갯수에 비례하는 시간 복잡도(O(n))를 가지는 알고리즘이 알려져 있고, 여기서는 재귀적 방법을 이용한 Welzl's algorithm을 구현한다.
int MinEncCir ( CfPt *P, int n, CfPt* boundary, int b, CfPt& center, double& rad ) {
// exiting cases
if ( b == 3 ) CalcCircle3 ( boundary, center, rad ); // a circle passing 3 points;
else if ( ( n == 1 ) && ( b == 0 ) ) {
rad = 0; b = 1;
center = boundary[0] = P[0];
}
else if ( ( n == 2 ) && ( b == 0 ) ) {
boundary[0] = P[0]; boundary[1] = P[1]; b = 2;
CalcCircle2 (boundary, center, rad ); // a circle with diagonal consisting of 2 points;
}
else if ( ( n == 0 ) && ( b == 2 ) ) {
CalcCircle2 ( boundary, center, rad );
}
else if ( ( n == 1 ) && ( b == 1 ) ) {
boundary[1] = P[0]; b = 2;
CalcCircle2 ( boundary, center, rad );
}
else {// general case; ( b < 3 ) && ( n + b > 2 )
// choose a random pivot;
int k = rand() % n;
if ( k != 0 ) SWAP( P[0], P[k] );
int b1 = MinEncCir ( &P[1], n - 1, boundary, b, center, rad );
if (!InCircle(P[0], center, rad) ) {
// Now, P[0] belongs to the boundary.
boundary[b++] = P[0];
return MinEncCir ( &P[1], n - 1, boundary, b, center, rad );
} else return b1;
}
return b;
}
FFT는 입력 데이터의 개수$(N)$가 2의 지수승으로 주어질 때 $O(N\log N)$의 연산만으로 빠르게 DFT을 수행하는 알고리즘이다. 주어진 $N$개의 data $\{F_0, F_1, ..., F_{N-1}\}$의 DFT $\{a_0, a_1,..., a_{N-1}\}$는
$$ a_k = \sum_{n = 0}^{N-1} F_n \exp\Big( -i \frac{2\pi nk }{N} \Big)=\sum_{n = 0}^{N-1} F_n (W^n )^k , \quad W = \exp\Big( -i \frac{2\pi}{N}\Big)$$ 로 표현된다. $N$이 2의 지수승으로 주어지면 이 합은 $n$이 짝수인 항과 홀수인 항으로 아래처럼 분리할 수 있다:
\begin{align} a_k &= \sum _{n=0}^{N/2-1} F_{2n} (W^{2n})^k+W^{k} \sum_{n=0}^{N/2-1} F_{2n+1}(W^{2n})^k\\ &=\text{(N/2개 even 항에 대한 DFT)} + W^k \times \text{(N/2개 odd 항에 대한 DFT)}\\ &= e_k + W^k o_k \end{align}
이는 원래의 문제가 처음 입력 데이터의 절반만 가지고 처리하는 두 개의 동일한 문제로 줄일 수 있음을 보여준다. 전형적인 divide and conquer 기법을 적용할 수 있는 구조이므로 매우 빠르게 연산을 수행할 수 있다(Cooley-Tukey algorithm). 역변환의 경우에는 twiddle factor $W$를 $W^*$로 바꾸면 되므로(전체 normalization factor $\frac{1}{N}$가 덧붙여진다) 코드에서 쉽게 수정할 수 있다: $\tt W_{im} \rightarrow - { W}_{im}$. 아래는 FFT를 재귀적으로 구현한 코드이다. 비재귀적 구현은 kipl.tistory.com/22에서 찾을 수 있다.
$F_n$ 인덱스의 이진비트가 twiddlefactor를 어떻게 곱해야 하는지 알려주고, 비트를 역으로 읽으면 계수가 나타나는 순서를 알 수 있게 한다.
#define TWOPI 6.28318530717958647693
void split(int N, double* data) {
if (N == 8) {
double t1 = data[4]; data[4] = data[1];
double t2 = data[2]; data[2] = t1; data[1] = t2;
t1 = data[5]; data[5] = data[3];
t2 = data[6]; data[6] = t1; data[3] = t2;
return;
}
double *tmp = new double [N / 2];
for (int i = 0; i < N / 2; i++) //copy odd elements to tmp buffer;
tmp[i] = data[2 * i + 1];
for (int i = 0; i < N / 2; i++) // move even elements to lower half;
data[i] = data[2 * i];
for (int i = 0; i < N / 2; i++) // move odd elements to upper half;
data[i + N / 2] = tmp[i];
delete [] tmp;
}
void fft4(double *re, double *im) {
double tr1 = re[0]-re[2]+im[1]-im[3];
double ti1 = im[0]-im[2]-re[1]+re[3];
double tr2 = re[0]+re[2]-re[1]-re[3];
double ti2 = im[0]+im[2]-im[1]-im[3];
double ti3 = im[0]-im[2]+re[1]-re[3];
double tr3 = re[0]-re[2]-im[1]+im[3];
re[0] = re[0]+re[2]+re[1]+re[3];
im[0] = im[0]+im[2]+im[1]+im[3];
re[1] = tr1; im[1] = ti1;
re[2] = tr2; im[2] = ti2;
re[3] = tr3; im[3] = ti3;
}
int fft ( int N, double* re, double* im ) {
if (N <= 0 || N & (N-1)) return 0; // N is not a power of 2;
if ( N < 2 ) return 1;
else if (N == 2) {
double tr1 = re[0]-re[1], ti1 = im[0]-im[1];
re[0] = re[0]+re[1]; im[0] = im[0]+im[1];
re[1] = tr1; im[1] = ti1;
return 1;
} else if (N == 4) {
fft4(re, im);
return 1;
}
split ( N, re);
split ( N, im);
fft ( N / 2, &re[ 0], &im[ 0] );
fft ( N / 2, &re[N / 2], &im[N / 2] );
for ( int k = 0; k < N / 2; k++ ) {
double Ere = re[k];
double Eim = im[k];
double Ore = re[k + N / 2];
double Oim = im[k + N / 2];
double Wre = cos ( TWOPI * k / N ); //real part of twiddle factor:W
double Wim = -sin ( TWOPI * k / N ); //imag part of twiddlw factor:W
double WOre = Wre * Ore - Wim * Oim;
double WOim = Wre * Oim + Wim * Ore;
re[k] = Ere + WOre;
im[k] = Eim + WOim;
re[k + N / 2] = Ere - WOre;
im[k + N / 2] = Eim - WOim;
}
return 1;
};
6개의 주파수 $\{f_i = 2, 5, 9, 11, 21, 29\text{Hz}\}$가 섞인 신호 $s(t)$를 $[0,1]$ 초 구간에서 일정한 간격으로 sampling 해서 64개 data를 얻고(Nyquist frequency = 32Hz), 이에 대해 FFT를 수행한다.
$$s(t)=\sum_{i=0}^{5} (i+1) \cos(2\pi f_i t) $$
int fft_test_main() {
const int samples = 64;
double re[samples], im[samples];
const int nfreqs = 6;
double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
for (int i = 0; i < samples; i++) {
double signal = 0;
for (int k = 0; k < nfreqs; k++)
signal += (k + 1) * cos(TWOPI * freq[k] * i / samples);
re[i] = signal;
im[i] = 0.;
}
fft(samples, &re[0], &im[0]);
for (int i = 0; i < samples; i++) TRACE("%d, %f, %f\n", i, re[i], im[i]);
return 0;
}
Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.
이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram($\tt ch$)과 평균을 계산할 수 있는 1차 culmuative histogram($\tt cxh$)을 입력으로 사용한다.
