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

동일한 디스크 두 개의 중심이 가벼운 막대로 연결되어 있고, 앞 디스크는 처음 각속도 ${\omega }_0$로 시계방향 회전을 한다. 이 두 디스크를 바닥에 놓았을 때 앞바퀴와 달리 뒷바퀴에 작용하는 마찰이 충분히 커서 미끄러지지 않고 구르는 운동을 한다. 두 디스크는 처음 질량중심이 오른쪽 가속도를 가지지만 A의 각속도는 감소를 한다. 

1. A에 작용하는 운동마찰력의 방향은 오른쪽: $f_k=\mu mg$
2. B에 작용하는 정지마찰력($f$) 방향은 왼쪽
3. B의 cm에 대한 회전운동(막대의 장력은 토크 기여가 없음): 정지마찰력만 토크에 기여  $$fR=\frac{1}{2}m{R}^2\frac{a}{R}\to f=\frac{m}{2}a$$
4. 계의 수평방향 cm 병진운동(수평방향 외력=마찰력)  $$f_k-f=(m+m)a\to a=\frac{2}{5}\mu g$$
5. 따라서 B에 작용하는 정지마찰력: $f=\frac{1}{5}\mu mg$
6. 바닥의 임의의 한 지점에 대한 각운동량 보존됨을 이용하면(바닥에서 작용하는 수평방향 마찰력은 토크를 만들지 못함) 최종적으로 두 바퀴가 공통으로 회전하는 각속도를 구할 수 있다. $$L_i=L_{A,i}=\frac{m{R}^2}{2}{\omega }_0$$ 같은 각속도로 구르기 시작할 때, 두 바퀴의 질량중심운동과 질량중심축에 대한 회전운동이 각운동량에 기여하므로 $$L_f=\frac{3}{2}m{R}^2\omega \times 2$$ $$\to \omega =\frac{1}{6}{\omega }_0$$

무거운 줄이 도르래에 걸쳐있고, 한쪽 끝에는 사람이 매달려 있다. 사람이 갑자기 줄에 대한 상대속도 $v_{rel}$로 위쪽으로 올라간다. 이때 사람의 지상에 대한 속도는? 줄의 길이는 $L$, 단위길이당 밀도는 $\lambda$, 그리고 사람의 질량은 $M$이다. 

사람이 올라가기 위해서는 줄에 힘(충격량=$J$)을 주어야 하고, 이 힘의 반작용으로 위로 올라간다. 사람이 준 힘 때문에 줄은 아래로 움직이게 된다. 

사람의 운동량 변화(위쪽 =+) $Mv-0=J$

줄의 운동량 변화(아래쪽=+) $(\lambda L)u-0=J$

따라서 $Mv=\lambda Lu$이고, $v_{rel}=v-(-u)=v+u$이므로 

$$v=\frac{\lambda L}{M+\lambda L}{v}_{rel}$$

Q2: 사람이 올라가기 시작하면 사람쪽 줄의 무게가 더 크므로 더 빨리 내려가려고 할 것이다. 따라서 일정한 시간이 지나면 지상에서 볼 때 사람은 더 이상 위로 올라가지 못하게 된다. 그때가 언제인가?

이를 해결하기 위해 사람과 줄의 운동방정식을 만들자. 바닥에서 잰 사람의 높이를 $y$(위쪽+), 도르래에서 잰 왼쪽 줄의 끝을 $y_1$(아래+), 오른쪽 줄의 끝을 $y_2$(아래+)라면, $y_1+y_2=L=const$이다. 그리고 사람이 줄로 부터 받는 힘을 $f(t)$라면 사람의 운동방정식은 

$$M{y}^{″}(t)=f(t)-Mg$$

그리고 줄의 운동은(아래쪽+, 줄의 총질량 $m=\lambda L$)

$$m{y}_1^{″}(t)=\lambda (y_1(t)-y_2(t))g+f(t)=2\lambda g{y}_1(t)-mg+f(t)$$

상대속도가 일정하므로 $y(t)$와 $y_1(t)$사이의 관계를 만들 수 있다. $$v_{rel}=y^{\prime }(t)+y_1^{\prime }(t)=const$$$$\to y_1^{″}(t)=-y^{″}(t)$$$$\text{and}y(t)+y_1(t)=v_{rel}t+C$$

상수 $C$는 줄과 사람이 처음 평형상태였음을 이용하면 $\lambda y_2(0)g=\lambda y_1(0)g+Mg$ 에서 

$$y_1(0)=\frac{m-M}{2\lambda }$$ 또 사람의 출발높이가 $y(0)=0$이므로 

$$C=y_1(0)=\frac{m-M}{2\lambda }$$

이제 앞에서 얻은 방정식을 이용해서 $y(t)$의 운동방정식에서 $f(t)$을 소거하면

$$\frac{m+M}{2\lambda g}{y}^{″}(t)-y(t)=-v_{rel}t$$

이 방정식의 일반해는 

$$y(t)=A\cosh (\alpha t)+B\sinh (\alpha t)+v_{rel}t,{\alpha }^2=\frac{2\lambda g}{m+M}$$

인데, $t=0$일 때 $y(0)=0$이므로 $A=0$이고, $y^{\prime }(0)=\frac{m}{m+M}{v}_{rel}$ 였으므로 $B=-\frac{1}{\alpha }\frac{M}{m+M}{v}_{rel}$. 따라서 사람의 높이는 

$$y(t)=v_{rel}(t-\frac{1}{\alpha }\frac{M}{m+M}\sinh (\alpha t))$$

사람이 더 이상 높이 올라갈 수 없는 상태가 되면 속도 $y^{\prime }(t)=0$이어야 한다. 출발에서 그때까지 걸린 시간은

 $$y^{\prime }(t)=v_{rel}(1-\frac{M}{m+M}\cosh (\alpha t))=0$$

$$\to t=\sqrt{\frac{m+M}{2\lambda g}}{\cosh }^{-1}\frac{m+M}{m}$$