Geometry & Recognition

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);
}

Forward fft: 

$$ {\tt fft}(x+iy)=X + i Y = \sum(x +i y) (W_r + i W_i)$$

$$=\sum  (x W_r - y W_i ) + i (x W_i + y W_r)$$

Inverse fft:

$$ {\tt ifft}(X+iY)= \frac{1}{N}\sum(X +i Y) (W_r - i W_i)$$

$$= \frac{1}{N} \sum  (X W_r + Y W_i ) + i (-X W_i + Y W_r)$$

$$ =\frac{1}{N}\sum \left(XW_r - (-Y) W_i\right) +(-i) \left(XW_i + (-Y)W_r\right) $$

$$ = \frac{1}{N}\left[{\tt fft}(X-iY) \right]^* = \frac{1}{N}\left[{\tt fft}((X+iY)^*)\right]^*$$

또는 

$$ \frac{1}{N} {\tt fft}(Y+iX)= \frac{1}{N} \sum(Y +i X) (W_r + i W_i)$$

$$=\frac{1}{N} \sum  (Y W_r - X W_i ) + i (Y W_i + X W_r)\\= y + i x$$

$$X_m = \sum_{n=0}^{N/2-1} x_{2n} W_N^{mn} +W_N^m \sum _{n=0}^{N/2-1} x_{2n+1} W_N^{mn}$$

$$X_{m+N/2} = \sum_{n=0}^{N/2-1} x_{2n} W_N^{mn} -W_N^m \sum _{n=0}^{N/2-1} x_{2n+1} W_N^{mn}$$

$$m=0,1,...,N/2-1$$

Butterfly Structure:

FFT for N=8:

/* sin(), cos() 함수를 $\log_2(N/4)$번만 호출하면 된다. 이마저도 2배각 공식을 쓰면 sqrt() 호출로 대체할 수 있다 */

#define PI	3.1415926535897932384
#define SWAP(a, b) {double temp = (a); (a) = (b); (b)=temp;}
int fft(int n, double x[], double y[]) {
    if ((n !=0) && ((n & (n-1)) != 0) ) return 0;
    // bit-reversal;
    for (int i = 0, j = 0; i < n-1; i++) {
        if (i < j) {
            SWAP(x[i], x[j]); SWAP(y[i], y[j]);
        }
        int k = n >> 1;
        while (k <= j) {
            j -= k; k >>= 1;
        }
        j += k;
    }
    // Danileson-Lanczos section;
    double rotwr = -1, rotwi = 0;
    for (int loop = 1; loop < n; loop <<= 1) {	
        int block = loop << 1;				// block-size;2->4->8->...
        if (loop > 2) {
            // rotation factor of twiddle factor;
            double theta = PI / loop;
            rotwr = cos(theta);
            rotwi = -sin(theta);
        } else if (loop == 2) {
            rotwr = 0;
            rotwi = -1;
        }	
        // starting twiddle factor;
        double wr = 1;
        double wi = 0;
        for (int j = 0; j < loop; ++j) { 
            // 각 block의 같은 위치에서 동일한 twiddle factor가 곱해진다는 
            // 사실을 이용함;
            for (int i = j; i < n; i += block) {
                int ip = i + loop;
                // step-1; X[ip]*W
                double xwre = x[ip] * wr - y[ip] * wi;
                double xwim = x[ip] * wi + y[ip] * wr;
                // step-2; X[ip] = X[i] - X[ip]*W;
                x[ip] = x[i] - xwre;
                y[ip] = y[i] - xwim;
                // step-3; X[i] = X[i] + X[ip]*W;
                x[i] += xwre;
                y[i] += xwim;
            };
            // 각 block의 다음 차례에 곱해지는 twiddle factor 계산; T = W * rotW;
            double tre = wr * rotwr - wi * rotwi;
            double tim = wi * rotwr + wr * rotwi;
            wr = tre;
            wi = tim;
        }
    }
    return 1;
    // 역변환: (1/N)*conjugate[FFT[conjugate(X)]] 
}