728x90

공기 저항이 없을 때 물체를 $v_0$ 속력으로 위로 던지면 최고점에 올라가는데 걸리는 시간과 다시 내려오는데 걸리는 시간은 동일하게 $t_{ff} = v_0/g$로 주어진다. 공기 저항이 있는 경우는 어떻게 될지 구체적으로 계산해보자.

반지름이 $R$인 공 모양의 물체가 속력의 제곱에 비례하는 공기 저항(끌림힘) $D= \frac{1}{2} C\rho_{air} A v^2 \approx  0.2 \rho_{air} \pi R^2 v^2$을 받을 때, 올라가는 동안 운동 방정식은 (물체가 받는 공기의 부력도 고려해야 하지만 여기서는 무시한다. 부력은 $g$을 약간 줄이는 효과를 만든다)

$$  m \frac{dv}{dt} = -mg -\ 0.2   \rho_{air} \pi R^2  v^2, $$

$$ \rightarrow  \quad \frac{dv}{dt} = -g (1 + \gamma^2 v^2) , \quad\quad  \gamma^2 =  0.2 \frac{\pi R^2 \rho_{air} }{mg} .$$

$\gamma$는 내려오는 과정에서 terminal speed의 역수를 의미한다. 시간에 대해 적분을 해서 속도를 구하면,

$$ v(t) = \frac{1}{\gamma} \tan \Big( -\gamma gt + \arctan(\gamma v_0) \Big)$$

최고점에 도달하는데 걸리는 시간은 $v(t_{up})=0$ 에서

$$ t_{up} = \frac{1}{\gamma g} \arctan(\gamma v_0) = \frac{v_0}{g} \Big(1 - \frac{(\gamma v_0)^2}{3}+....\Big) $$

로 주어지므로 공기저항이 없을 때보다 더 짧다. 최고점의 높이는

$$ h_{max} = \int_0^{t_{up}} v(t) dt = \frac{1}{\gamma^2 g} \ln \sqrt{ 1 + (\gamma v_0)^2 } =\frac{v_0^2}{2g} \Big( 1- \frac{(\gamma v_0)^2}{2 }+...\Big)$$

로 주어지므로 역시 공기 저항이 없을 때보다 낮다.

 

다시 내려오는 과정에서 중력과 끌림힘이 반대방향이므로 운동 방정식은

$$ \frac{dv}{dt}= -g (1 -\gamma^2 v^2)$$

로 주어지고, 이를 적분하면 (출발 시간을 $t=0$으로)

$$ v(t) = -\frac{1}{\gamma} \tanh (\gamma g t)$$

이를 다시 적분하면 낙하시간에 따른 높이를 얻을 수 있다: $h(0)=h_{max}$

$$ h(t) = \frac{1}{\gamma^2 g} \ln \frac{ \sqrt{ 1+ (\gamma v_0)^2 } }{ \cosh( \gamma gt)}$$

따라서 바닥에 떨어지는데 걸리는 시간 $t_{dn}$은: $ h(t_{dn}) = 0$

$$t_{dn} = \frac{1}{\gamma g}\text{arccosh}\sqrt{1+(\gamma v_0)^2} =\frac{1}{\gamma g} \ln \left(\gamma v_0 + \sqrt{1+ (\gamma v_0)^2 } \right) = \frac{v_0}{g}\Big( 1  -  \frac{ (\gamma v_0)^2 }{6}+...\Big)$$

 

두 시간을 비교해보면 위로 던져진 물체가 최고점에 올라가는데 걸리는 시간보다 다시 내려오는데 걸리는 시간이 더 길다는 것을 볼 수 있다:

$$ {t_{dn} - t_{up}} \approx  \frac{(\gamma v_0)^2}{6} t_{ff}$$

 

$v_0 =10\text{m/s}$로 ($\rightarrow t_{ff} = 1.02\text{s}$, $v_{terminal} \approx 36.6\text{m/s}$) 야구공을 공중으로 던지는 경우를 예로 들면, 차이는 대략 0.013초 정도로 계산된다. 던지는 속력이 종단속력에 가깝거나 더 크면 근사식을 사용할 수 없고 정확한 계산식을 이용해야 한다.

Posted by helloktk

댓글을 달아 주세요

728x90

이미지 처리 과정에서 미분에 해당하는 그래디언트 필드(gradient field: $g_x$, $g_y$ )를 이용하면 이미지 상의 특징인 corner, edge, ridge 등의 정보를 쉽게 얻을 수 있다. 이미지의 한 지점이 이러한 특징을 가지는 특징점이 되기 위해서는 그래디언트 필드의 값이 그 점 주변에서 (3x3나 5x5정도 크기의 window) 일정한 패턴을 유지해야 한다. 이 패턴을 찾기 위해서 그래디언트 필드에 PCA를 적용해보자. 수평과 수직방향의 그래디언트 field인 $g_x$와 $g_y$ 사이의 covariance 행렬은 다음 식으로 정의된다:

$$ \Sigma = \left [ \begin {array}{cc} < g_x^2 > & < g_x g_y > \\    <g_x g_y> & <g_y^2 > \end {array}\right] =\left [ \begin {array}{cc} s_{xx} & s_{xy} \\ s_{xy} & s_{yy}\end {array}\right];$$

$<...> = \int_{W}(...) dxdy$는 픽셀 윈도에 대한 적분을 의미한다. $\Sigma$의 eigenvalue는 항상 음이 아닌 값을 갖게 되는데 (matrix 자체가 positive semi-definitive), 두 eigenvalue이 $λ_1$, $λ_2$면

$$λ_1 + λ_2 = s_{xx} + s_{yy} \ge 0, \quad \quad    λ_1  λ_2 = s_{xx} s_{yy} - s_{xy}^2 \ge0 $$

을 만족한다 (완전히 상수 이미지를 배제하면 0인 경우는 없다). eigenvalue $λ_1$, $λ_2$는 principal axis 방향으로 그래디언트 필드의 변동(분산)의 크기를 의미한다. edge나 ridge의 경우는 그 점 주변에서 잘 정의된 방향성을 가져야 하고, corner의 경우는 방향성이 없어야 한다. edge나 ridge처럼 일방향성의 그래디언트 특성을 갖거나 corner처럼 방향성이 없는 특성을 서로 구별할 수 있는 measure가 필요한데, $λ_1$과 $λ_2$를 이용하면 차원이 없는 measure을 만들 수 있다. 가장 간단한 차원이 없는 측도(dimensionless measure)는  eigenvalue의 기하평균과 산술평균의 비를 비교하는 것이다.

$$ Q = \frac { {λ_{1} λ_{2}} }{ \left( \frac {λ_{1}+λ_{2}}{2} \right)^2} = 4\frac { s_{xx} s_{yy} - s_{xy}^2}{(s_{xx} + s_{yy})^2};$$

기하평균은 산술평균보다도 항상 작으므로

$$ 0 \le Q \le 1 $$

의 범위를 갖는다. 그리고 $Q$의 complement로

$$P = 1-Q = \frac{(s_{xx}-s_{yy})^2 + 4 s_{xy}^2}{(s_{xx}+s_{yy})^2};$$를 정의할 수 있는 데 $0 \le P \le 1$이다. $Q$와 $P$의 의미는 무엇인가? 자세히 증명을 할 수 있지만 간단히 살펴보면 한 지점에서 $Q \rightarrow 1$이려면 $λ_{1} \approx λ_{2}$이어야 하고, 이는 두 주축이 동등하다는 의미이므로 그 점에서는 방향성이 없는 코너의 특성을 갖게 된다. 반대로 $Q \rightarrow 0$이면 강한 방향성을 갖게 되어서 edge(ridge) 특성을 갖게 된다.

 

실제적인 응용으로는 지문 인식에서 지문 영역을 알아내거나 (이 경우는 상당이 큰 윈도를 사용해야 한다) 또는 이미지 텍스쳐 특성을 파악하기 위해서는 이미지를 작은 블록으로 나누고 그 블록 내의 미분 연산자의 균일성을 파악할 필요가 있는데 이 차원이 없는 측도는 이미지의 상태에 상관없이 좋은 기준을 주게 된다.

 

참고 논문:

Image field categorization and edge/corner detection from gradient covariance
Ando, S.

Pattern Analysis and Machine Intelligence, IEEE Transactions on
Volume 22, Issue 2, Feb 2000 Page(s):179 - 190

 

** 네이버 블로그 이전;

'Image Recognition > Fundamental' 카테고리의 다른 글

Is Power of 2  (0) 2021.02.12
Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Gradient Covariance  (0) 2021.01.27
점증적인 cosine/sine 값 계산  (0) 2020.12.28
Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
Posted by helloktk

댓글을 달아 주세요

728x90

기하 알고리즘을 특히 폴리곤 관련 알고리즘을 테스트하기 위해서는 폴리곤 데이터를 만들어야 한다. 여기서 2차원의 점들을 이용하여 간단히 simple polygon을 만드는 방법을 소개한다.

 

단순 폴리곤을 만들기 위해서는 점들을 정렬을 해야 한다. 각각의 점들을 x 축에 프로젝션을 시키면 x 축에 대해서 정렬이 된다(같은 x 값인 경우에는 y의 크기대로 정렬). 따라서 x 값으로 정렬된 점들을 쭉 이어나가면 하나의 poly-line을 만들 수 있다. poly-line의 끝을 잇는다고 해서 단순 폴리곤이 항상 만들어지지는 않는다. 이를 해결하기 위해서는 점들을  한 직선을 기준으로 위-아래로 분리하고, 직선의 윗부분은 x 값이 큰 순서대로 정렬을 하고, 아래 부분은 x 값이 작은 순서대로 정열을 하면 폴리 라인이 아닌 단순 폴리곤을 얻을 수 있다. 직선은 어떻게 잡을까? 이것은 가장 작은 x 값을 갖는 점과 가장 큰 x 값을 같은 점을 잇는 직선을 생각하면 편리하다.

// cross(BC, BA)
// A->B-C: 반시계(C가 AB 왼편: > 0), 시계(C가 AB오른편: < 0), 일직선(0);
int CCW(CPoint A, CPoint B, CPoint C) {
    C.x -= B.x; C.y -= B.y; 
    A.x -= B.x; A.y -= B.y;
    return C.x * A.y - C.y * A.x;
}
int comp(const void *A, const void *B) { //x->y의 크기 순서대로 정렬;
    int ax = ((CPoint *)A)->x - ((POINT *)B)->x ;
    if (ax > 0) return  1;
    if (ax < 0) return -1;
    int ay = ((CPoint *)A)->y - ((CPoint *)B)->y ;
    if (ay > 0) return  1;
    if (ay < 0) return -1;
    return 0;
}
void MakeSimplePolygon(std::vector<CPoint>& pts, std::vector<CPoint>& spoly) {
    if (pts.size() < 1) return;
    int maxid = 0, minid = 0;
    for (int i = 1; i < pts.size(); i++) {
        if (pts[i].x > pts[maxid].x) maxid = i;
        if (pts[i].x < pts[minid].x) minid = i;
    }
    std::vector<CPoint> LP, RP;
    for (int i = 0; i < pts.size(); i++) {
        if (i == maxid || i == minid) continue;
        // 기준선의 왼편(LP)/오른편(RP)에 놓인 점 분리;
        if (CCW(pts[minid], pts[maxid], pts[i]) >= 0) LP.push_back(pts[i]); 
        else RP.push_back(pts[i]);
    }
    if (LP.size() > 0) 
        qsort(&LP[0], LP.size(), sizeof(CPoint), comp);
    if (RP.size() > 0)    
        qsort(&RP[0], RP.size(), sizeof(CPoint), comp);
	
    spoly.clear();
    spoly.push_back(pts[minid]);
    for (int i = 0; i < LP.size(); i++) spoly.push_back(LP[i]);
    spoly.push_back(pts[maxid]);
    for (int i = RP.size() - 1; i >= 0; i--) spoly.push_back(RP[i]);
    // spoly = clockwise simple polygon;
};

다른 방법으로는 위와 같은 직선을 만들고, 아랫부분의 점들의 convex hull을 구성하고, 나머지 점들을 마찬가지의 방법으로 정렬을 하여 폴리곤을 만들거나, 아니면 한 점을 잡아서(맨 아래 오른쪽) 그 점을 기준으로 해서 나머지 점들을 각도에 대해서 정렬을 한 다음 그 정렬된 순서대로 폴리곤을 구성하면 된다.

기준점에 대한 각도를 정렬하여서 만든 예(동일한 각도가 생기는 경우에는 기준점에서의 거리로 비교)

**네이버 블로그에서 이전;

'Computational Geometry' 카테고리의 다른 글

Minimum Bounding Rectangle  (0) 2021.03.15
Minimum Enclosing Circle  (0) 2021.03.01
단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
Posted by helloktk

댓글을 달아 주세요

728x90

평면 상에 주어진 점들의 convex hull을 구하는 알고리즘 중의 하나인 Graham scan에서는 먼저 주어진 점들을 한 점을 기준으로 각도로 정렬하는 과정이 필요했다. 그러면 점들이 순차적으로 연결된 단순 다각형에서는 sorting과정이 없이 Graham scan 알고리즘을 바로 적용이 가능한 것일까? 이에 대한 counter example은 쉽게  찾을 수 있다. 단순 폴리곤에 대해서도 항상 각도에 대한 정렬 과정이 필요한 것인가? 답은 아니다. 정렬 과정이 없이도 단순 폴리곤에 대해서는 쉽게 convex hull을 구할 수 있는 알고리즘이 존재한다. 정렬 과정이 없이 단순 폴리곤의 convex hull을 찾는 알고리즘에 대해서 알아보자.

Melkman Algorithm;

우선 폴리곤의 이웃하는 세 꼭짓점을 잡아서, 반시계 방향의 삼각형을 구성한다. 이 삼각형을 deque에 넣는다 (bottom = top). 폴리곤을 순환하면서 새 점이 들어올 때 이미 만들어진 convex hull의 내부점이 아니면 이 점이 포함되도록 convex hull을 업데이트한다: Graham scan 알고리즘의 scanning 과정을 bottom을 기준으로 반시계 방향으로 convexity를 만족시킬 때까지 bottom을 제거하고, top을 기준으로는 시계방향으로 convexity를 만족시킬 때까지 top을 제거한다. 이 과정이 끝나면 새 점을 deque에 추가한다. 이 과정을 나머지 모든 점들에 대해서도 진행한다.

int LEFTSIDE(CPoint A, CPoint B, CPoint C){ // cross(AB, AC) > 0: C가 AB 대해서 왼쪽
	return ((B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x)) > 0;
}
int convexhull_spolygon(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    int N = pts.size();
    if (N < 3) return 0;
    std::vector<CPoint> Q(2 * N + 1);  // deque;
    int bot = N - 2, top = bot + 3;
    // |0|1|......|N-2|N-1|N+0|N+1|N+2|..............|2N|;
    //              bot    top    
    Q[bot] = Q[top] = pts[2];
    if (LEFTSIDE(pts[0], pts[1], pts[2])) { //2->0->1->2; Q에 ccw로 입력;
        Q[bot + 1] = pts[0]; Q[bot + 2] = pts[1];
    } else {
        Q[bot + 1] = pts[1]; Q[bot + 2] = pts[0];
    }
    int i = 2; // last touched index;
    while (++i < N) {
        // 기존의 convex_hull에 들어있으면 제외.
        if (LEFTSIDE(Q[bot + 0], Q[bot + 1], pts[i]) && 
            LEFTSIDE(Q[top - 1], Q[top + 0], pts[i]))
            continue;
        // bottom에 대해서 ccw 방향으로 체크(bot 증가 방향)
        // pts[i]가 (bot)->(bot+1)라인의 오른쪽에 있으면, 기존의 bottom을 제외;
        while (!LEFTSIDE(Q[bot], Q[bot + 1], pts[i]) && (bot < top)) bot++;
        Q[--bot] = pts[i];
        // 추가점에서 top->top-1의 ccw 방향으로 convexity 체크하여 만족하지 않은 경우
        // top을 감소
        while (!LEFTSIDE(Q[top - 1], Q[top], pts[i]) && (top > bot)) top-- ;
        Q[++top] = pts[i];
    }
    hull.resize(top - bot);
    for (int i = hull.size() - 1; i >= 0; --i) hull[i] = Q[i + bot];
    return hull.size();
};

**네이버 블로그 이전;

 

 

'Computational Geometry' 카테고리의 다른 글

Minimum Enclosing Circle  (0) 2021.03.01
단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Posted by helloktk

댓글을 달아 주세요

728x90

단순 다각형의 무게중심(center of gravity or center of mass)은 다각형을 균일한 밀도의 판으로 생각했을 때 판의 무게중심과 같다. 가장 단순한 다각형인 삼각형의 무게중심은 세 꼭짓점의 산술평균으로 표현된다.

$$ \text{CoG} = \frac{1}{3} ({{\bf P} + {\bf Q} + {\bf R}}).$$

다각형은 삼각형으로 분할되므로 이 분할된 삼각형의 무게중심을 이용하면 쉽게 계산할 수 있다. 분할된  삼각형의 무게중심을 면적으로 가중치를 준 평균값이 다각형의 무게중심이 된다. 

 

실제 계산에서는 다각형을 삼각 분할하지 않고도 간단한 방법에 의해서 무게중심을 구할 수 있다. 원점과 다각형의 각 변의 꼭짓점을 이용해서 삼각형들을 구성하면 원래의 다각형을 겹치게(원점이 내부에 있으면 겹침이 없다) 분할할 수 있다. 분할된 삼각형으로 무게중심을 구할 때 겹치는 영역의 기여를 제거해야 한다. 그런데 다각형 밖의 영역을 분할하는 삼각형은 다각형 내부를 분할하는 삼각형과는 다른 orientation을 가지게 된다. 삼각형의 면적은 한 꼭짓점을 공유하는 두 변의 외적에 비례하므로, 반대의 orientation을 갖는 삼각형은 자동으로 반대 부호의 면적을 가지게 된다. 따라서 분할된 삼각형 면적 가중치를 외적으로 주어서 무게중심을 구하면 겹치는 영역이 자동으로 상쇄되는 효과를 얻을 수 있다.

노란색 영역에 있는 분할 삼각형은 음의 면적

$$\begin{align} \text{CoG} &= \frac{1}{\text{다각형 면적}} \sum (\text{삼각형 CoG})  (\text{면적})  \\ &= \frac{1}{\text{다각형 면적}} \sum \frac{1}{3} \left( {\bf V}_i + {\bf V}_{i+1} + {\bf O}\right ) \frac{ {\bf V}_{i} \times {\bf V}_{i+1} }{2} \\ &= \frac{1}{3}\frac{1}{   \text{다각형 면적} }\sum ( {\bf V}_{i} + {\bf V}_{i+1}) \frac{{\bf V}_{i} \times {\bf V}_{i+1}}{2} \end{align}$$ 

다각형의 면적($=\sum \frac{1}{2}({\bf V}_i \times {\bf V}_{i+1})$)을 구할 때 삼각형과 동일하게 orientation에 따라 부호를 포함하도록 설정하면 다각형의 면적 부호가 삼각형의 면적 부호로 상쇄되므로 다각형의 orientation에 무관하게 성립하는 공식이 된다.

POINT polygon_centroid(POINT V[], int N) {
    double cx = 0, cy = 0, area2 = 0;
    for(int i = 0, j = N - 1; i < N; j = i++) {
        double tri_area2 = V[i].x * V[j].y - V[i].y * V[j].x;
        cx += (V[i].x + V[j].x) * tri_area2;
        cy += (V[i].x + V[j].x) * tri_area2;
        area2 += tri_area2;
    }
    cx /= 3 * area2;
    cy /= 3 * area2;
    return POINT(cx, cy)
};

** 네이버 블로그에서 수정 이전;

'Computational Geometry' 카테고리의 다른 글

단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Posted by helloktk

댓글을 달아 주세요

728x90

평면상의 다각형(모서리의 교차가 없는 단순 다각형)의 면적을 구하는 것은 단순하지 않을 것처럼 보이지만 계산식은 무척이나 간단하게 주어진다. 기본적인 아이디어는 다각형에 임의의 점을 찍으면 이 점과 이웃한 두 개의 꼭짓점으로 형성이 되는 삼각형의 합으로 다각형을 분할할 수 있다. 분할된 삼각형의 면적을 구하여 합산하면 다각형의 면적을 구할 수 있다.

 

세 점 ${\bf P, Q, R}$(이 순서대로 반시계방향으로 배열)이 만드는 삼각형의 면적은

$$\text {삼각형의 면적}=  \frac{1}{2} ({\bf R} - {\bf  Q})  \times ( {\bf P} - {\bf Q}); \quad  (\Rightarrow \text {시계 방향이면 면적이 음수})$$

로 주어지므로, 꼭짓점이 ${\bf P}_0(x_0, y_0), {\bf P}_1(x_1, y_1),....$(반시계 방향)으로 주어지는 $N$각형의 면적은 아래와 같이 주어진다. 

$$\begin{align} \text{다각형 면적} &= \sum \text{각 삼각형의 면적} \\ &= \frac{1}{2}\sum ({\bf P}_{i+1}-{\bf Q})\times ({\bf P}_{i}-{\bf Q})\quad \quad ({\bf Q}\text{는  임의의 점})\\ &= \frac{1}{2} \sum\left(  {\bf P}_{i+1} \times {\bf P}_{i} - {\bf P}_{i+1}\times {\bf Q} + {\bf P}_{i} \times {\bf Q}\right) \\ &=\frac{1}{2} \sum {\bf P}_{i+1} \times {\bf P}_{i} \\ &= \frac{1}{2}\sum \left( x_{i+1} y_{i} -x_{i}y_{i+1} \right) \end{align}$$

이 결과는 $\bf Q$에 무관하다. 다각형의 꼭짓점이 시계 방향으로 정렬이 된 경우는 면적이 음수로 나온다(윈도우 DC는 위-아래가 역전되어 있으므로 orientation이 반대로 보인다). 그리고 이 공식은 단순 다각형에만 적용이 되고 모서리의 교차가 있는 경우에는 적용이 되지 않는다.

double simple_polygon_area2D(POINT point[], int N) {
    double area = 0;
    for (int i = 0, j = N - 1; i < N; j = i++) 
        area += (point[i].x * point[j].y) - (point[i].y * point[j].x);
    area /= 2;
    // return area;  // signed area;
    return area < 0 ? -area: area;
}

**네이버 블로그 수정 이전;

'Computational Geometry' 카테고리의 다른 글

단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Incremental Delaunay Triangulation  (1) 2020.12.01
Posted by helloktk

댓글을 달아 주세요

728x90

마찰 없이 바닥에서 움직일 수 있는 경사면$(M)$에 올려진 물체$(m)$가 정지상태에서 운동을 시작한다. 바닥에 내려왔을 때의  속력을 구하기 위해서 다음 과정을 따라 계산을 한다.

1. 다 내려와서는 물체나 경사면 모두 바닥에서 서로 반대 방향으로 움직인다(물체: $v_b$, 경사면:$V$)

2. 운동량 보존을 고려하면,

$$ m v_b = MV$$

3. 역학적에너지가 보존되므로

$$\frac {1}{2} m v_b^2 + \frac {1}{2} MV^2 = mgh$$

두 식을 이용하면 경사면이 움직이는 속력은

$$V = \sqrt{  \frac {2gh}{(1+\frac {M}{m})\frac {M}{m}} }$$

로 주어진다.

그런데 kipl.tistory.com/118 에서 얻은 결과를 이용하면 (좀 복잡하지만: $\theta \rightarrow \alpha$)

$$V =\sqrt{ \frac {2gh}{(1+\frac {M}{m})( \frac {M}{m}  + \sin^2 \alpha) } }\cos\alpha$$

로 주어진다.

두 결과를 비교하면 다르다. 어디서 잘못이 들어온 것일까? 결과만 놓고 보았을 때 어느 계산이 맞는지 직관적으로 설명할 수 있는가?

 

Posted by helloktk

댓글을 달아 주세요

728x90

막대의 끝과 쇠공은 같은 높이에서 떨어진다. 바닥에 먼저 닿는 물체는?

1. 막대

2. 쇠공

3. 동시에

 

Posted by helloktk

댓글을 달아 주세요

728x90

(매우 무거운) 롤 화장지를 줄을 이용해서 벽에 걸어두었다. 벽과 화장지 사이에 마찰은 줄에 걸리는 장력을 줄이는데 도움이 될까? 

1. 도움이 된다.

2. 오히려 장력이 더 커진다.

3. 마찰이 있던 없든 장력은 영향을 안 받는다.

화장실을 사용할 때마다 고민을 해보면 변비 해소에 도움이 될 수도 있다

Posted by helloktk

댓글을 달아 주세요

728x90

아날로그 신호로부터 디지털 데이터로 바꿀 때 보통 시간당(sound) 또는 공간적 거리당(image) 일정한 횟수로  데이터를 수집한다. 이 비율을 sampling rate이라 한다. 그런데 sampling rate을 바꾸어야 할 경우가 있다. 1초마다 한 번씩 데이터를 모았는데 실제로 0.5초마다 수집된 데이터가 필요가 한 경우에는 기존의 수집된 데이터와 데이터 사이에서 원 아날로그 신호 값을 알아내야 한다. sampling rate을 바꿀 때 기존의 데이터를 이용해서 원 아날로그 신호에서 수집한 것처럼 새로운 데이터를 만들어 내는 과정이 resampling이다. 이미지에서는 원본을  확대하는 up sampling이나 축소하는 down sampling 등이 resampling의 일종이다. resampling을 하기 위해서는 기존 데이터를 이용해서 데이터와 데이터 사이에서의 값을 추정하는 interpolation이 필요하게 된다. 많이 사용하는 interpolation 방법으로는 nearest neighbor interpolation, linear interpolation, cubic interpolation, spline interpolation 등이 있다.  여기서는 cubic interpolation을 이용한 이미지 resampling을 살펴본다.

 

$x$축 위의 4 점 $x=-1$, $x=0$, $x=1$, $x=2$에서 값이 $p_0$, $p_1$, $p_2$, $p_3$으로 주어진 경우 $0 \le x \le1$ 구간에서 cubic interpolation 함수는 

$$ f(x)=ax^3 +b x^2 + cx +d$$

라고 쓰면, 우선 $x=0$에서 $p_1$, $x=1$에서 $p_2$를 지나므로

$$f(0) = d = p_1, \quad\quad f(1) = a + b+ c+d = p_2$$

를 만족해야 하고, 양끝에서 도함수 값은 주변 입력점의 평균변화율로 주면, 다음식을 만족한다.

$$  f'(0) = c = \frac{p_2 - p_0}{2}, \quad \quad f'(1) = 3a + 2b + c =\frac{p_3 - p_1}{2}$$

이를 이용해서 $a,b,c,d$를 구하면,

$$a= \frac{1}{2}( -p_0 +3 p_1 -3p_2 +p_3), ~b = \frac{1}{2}(2p_0 -5 p_1 + 4p_2 -p_3),~ c=\frac{1}{2} (-p_0 +p_2), ~d=p_1. $$

따라서 cubic interpolation 함수는

\begin{align} f(x) &=\frac {1}{2}\big(- x^3 + 2x^2 -x\big) p_0+\frac {1}{2}\big(3x^3 - 5x^2 +2 \big) p_1\\&+\frac {1}{2}\big( -3x^3 + 4 x^2 + x\big) p_2 +\frac {1}{2}\big( x^3 - x^2 \big) p_3 \end{align}

로 쓰인다. 이는 다음처럼 정의된 kernel 함수를(Catmull-Rom filter) convolution한 형태로도 표현할 수 있다.

$$ K(x) = \left\{ \begin{array}{ll} \frac{1}{2} (3|x|^3 - 5|x|^2 +2) &  |x| <1 \\  \frac{1}{2}(-|x|^3+5|x|^2 -8|x|+4 ) & 1 \le |x|<2 \\ 0 & \text{otherwise} \end{array} \right. $$

$$  f(x) = p_0 K(-1-x) + p_1 K(-x) + p_2 K(1-x) + p_3 K(2-x),\quad (0\le x <1)$$

여기서는 좀 더 간단히 표현하기 위해서 $$ f(x) = m_0 p_0 + m_1 p_1 + m_2 p_2 + m_3 p_3 $$ 로 표현한다. 물론 $m_i$는 $0\le x <1$의 함수이다. $f(x)$는 그림에서 보는 것처럼 $p_1$과 $p_2$를 통과한다.

구간 끝에서 도함수를 $$f'(0) = \frac{p_2 - p_0}{2}, ~f'(1) = \frac{p_3 - p_1}{2}$$로 설정했기 때문에 인접 구간(eg. $[p_2,p_3]$)에서의 interpolation과 경계에서 같은  도함수를 가지게 되어 곡선이 꺽임없이 부드럽게 연결이 되는 특성을 가진다(예를 들면, $p_2$와 $p_3$을 보간하는 cubic interpolation함수는 $p_2$에서 미분값이 $(p_3 - p_1)/2$인데, 이는 $p_1,p_2$구간의 interpolation  함수가 $p_2$에서 가지는 미분값과 같다).

 

이를 2차원 공간으로 확장하면, 평면 위의 격자점 $\{(x, y)| x, y=-1,0,1, 2\}$에서 16개의 값 $\{ p_{ij} |i, j=0,1,2, 3\}$가 주어진 경우 정사각형 영역 $ 0\le x, y \le 1$에서 bicubic interpolation은 먼저 4개의 행 $y=-1,0,1,2$에서 각각 $x$ 방향의 cubic interpolation 된 값을 구한다:

\begin{align} q_0 =f(x, y=-1)= m_0 p_{00} + m_1 p_{10} + m_2 p_{20} + m_3 p_{30} \\ q_1 =f(x, y=0)= m_0 p_{01} + m_1 p_{11} + m_2 p_{21} + m_3 p_{31} \\ q_2 =f(x, y=1)= m_0 p_{02} + m_1 p_{12} + m_2 p_{22} + m_3 p_{32} \\ q_3 =f(x, y=2)= m_0 p_{03} + m_1 p_{13} + m_2 p_{23} + m_3 p_{33} \end{align}

4개 행에서 구해서 값을 이용해서 다시 $y$ 방향으로 cubic interpolation 한 값을 구하면 bicubic interpolation이 완성된다: $$ q =f(x,y)= q_0 n_0 + q_1 n_1 + q_2 n_2 + q_3 n_3$$ 여기서 $n_i$는 $m_i$에서 $x$를 $y$로 치환한 값이다.

 

원본 이미지를 확대/축소 또는 변형 등의 변환 과정을 거쳐서 출력 이미지를 만들 때 원본 이미지의 픽셀 값을 resampling 하는 과정이 들어온다. 이때 원본 이미지의 픽셀과 픽셀 사이의 값이 필요한 경우가 발생하여 적절한 interpolation이 필요한데 많이 쓰이는 interpolation 중의 하나가 bicubic interpolation이다. 아래는 bicubic interpolation을 이용한 이미지 resampling을 수행하는 코드다.

 

interpolation에서 실수(float) 연산을 하지 않고 정수 연산만 사용하면서도 일정한 정밀도로 소수점 아래 부분의 정보를 유지하기 위해 나누기 전에 미리 256을 곱한 후에 나눈다(256 대신에 적당히 큰 수를 선택해도 되는데 2의 지수승을 잡으면 곱셈/나눗셈을 shift 연산으로 대체할 수 있는 장점이 있다). 이렇게 하면 나눈 몫에 $\tt 0xFF$ 비트 마스크를 적용할 때 남는 부분이 소수점 아랫부분을 표현한다. 정밀도는 $1/256$ 까지다. 중간 과정에서 소수점 이하 부분끼리 곱셈을 할 때는 항상 $256$으로 나누어서 $255$ 이하가 되도록 만들어야 한다. 최종 결과에서도 다시 $256$으로 나누기를 해주어야 된다. bicubic인 경우는 $x/y$ 양방향으로 적용하므로 $256\times 256$으로 나누어야 하고 cubic interpolation 계수에서 $1/2$이 들어오므로 추가로 $4$만큼 더 나누어 주어야 한다(코드의 마지막 결과에서 shift 연산 "$\tt >> 18$"이 들어온 이유다).

 

bicubic interpolation을 적용할 때 $\tt y=0$이나 $\tt x=0$에서는 이전 행이나 열이 없으므로 자신을 반복하는 방식으로 처리해 주어야 하고, 또 $\tt y=rows-2, rows-1$이나 $\tt x=cols-2, cols-1$일 때도 비슷한 처리가 있어야 한다.

int resample_bicubic ( BYTE *src, int cols, int rows,
                       BYTE *des, int newcols, int newrows ) {
    if (cols < 4 || rows < 4)
        return resample_bilinear(src, cols, rows, des, newcols, newrows);

    int ixn = cols - 4;       
    BYTE *pa, *pb, *pc, *pd;
    for (int j = 0; j < newrows; j++) {
        int yy = ( ( j * rows ) << 8 ) / newrows;
        int yp = yy >> 8;                        // src pixel y-position;
        int dy = yy & 0xFF;
        int dy2 = (dy * dy) >> 8;
        int dy3 = (dy2 * dy) >> 8;
        int n0 = -dy3 + 2 * dy2 - dy;
        int n1 = 3 * dy3 - 5 * dy2 + 2 * 256;
        int n2 = -3 * dy3 + 4 * dy2 + dy;
        int n3 = dy3 - dy2;
        //
        if (yp == 0) pa = src + yp * cols;
        else         pa = src + (yp - 1) * cols;
        pb = src + yp * cols;
        if (yy < rows - 2) {
            pc = src + (yp + 1) * cols;
            pd = src + (yp + 2) * cols;
        } else if (yp < rows - 1) {
            pc = src + (yp + 1) * cols;
            pd = pc;
        } else {
            pc = src + yp * cols;
            pd = pc;
        }
        for (int i = 0; i < newcols; i++) {
            int xx = ( ( i * cols ) << 8 ) / newcols;
            int xp = xx >> 8;                    // src pixel x-position;
            int dx = xx & 0xFF;
            int dx2 = ( dx * dx ) >> 8;
            int dx3 = ( dx2 * dx ) >> 8;
            int m0 = -dx3 + 2 * dx2 - dx;
            int m1 = 3 * dx3 - 5 * dx2 + 2 * 256;
            int m2 = -3 * dx3 + 4 * dx2 + dx;
            int m3 = dx3 - dx2;
            int p = (xp == 0) ? 0 : (xp < ixn) ? xp - 1: ixn;    // p+3 <= ixn+3=cols-1;
            int a = ((m0 * pa[p] + m1 * pa[p + 1] + m2 * pa[p + 2] + m3 * pa[p + 3]) * n0 +
                     (m0 * pb[p] + m1 * pb[p + 1] + m2 * pb[p + 2] + m3 * pb[p + 3]) * n1 +
                     (m0 * pc[p] + m1 * pc[p + 1] + m2 * pc[p + 2] + m3 * pc[p + 3]) * n2 +
                     (m0 * pd[p] + m1 * pd[p + 1] + m2 * pd[p + 2] + m3 * pd[p + 3]) * n3) >> 18;
            *des++ = (a > 0xFF) ? 0xFF: (a < 0) ? 0: a;
        }
    }
	return 1;
}

bicubic interpolation을 하기 위해서는 4점이 필요하므로 폭이나 높이가 이보다 작은 경우는 bilinear interpolation을 사용한다. 다음 코드는 fixed-point bilinear interpolation을 구현한 코드다.

더보기
int resample_bilinear(BYTE *src, int cols, int rows,
                      BYTE *des, int newcols, int newrows ) {
    for (int j = 0; j < newrows; j++ ) {
        int yy = ((j * rows ) << 8) / newrows;
        int yp = yy >> 8;   // integer part;
        int dy = yy & 0xFF; // fractional part;
        for (int i = 0; i < newcols; i++) {
            int xx = ((i * cols ) << 8) / newcols;
            int xp = xx >> 8;
            int dx = xx & 0xFF;
            int pos = yp * cols + xp; //src position;
            int p00 = *(src + pos);
            int p10 = *(src + pos + 1);
            int p01 = *(src + pos + cols);
            int p11 = *(src + pos + cols + 1);
            int val = ((p11 * dx + p01 * (256 - dx)) * dy
                    + (p10 * dx + p00 * (256 - dx)) * (256 - dy)) >> 16;
            *des++ = val > 255 ? 0xFF: val < 0 ? 0 : val;
        }
    }
    return 1;
}

원본 컬러 이미지
2.5배 확대 영상

kipl.tistory.com/55

 

Bicubic Interpolation

이미지 처리에서 픽셀 좌표는 간격이 1인 2차원 그리드의 교차점으로 볼 수 있다. 이미지를 확대하거나 축소할 때, 픽셀과 픽셀 사이 구간에서 값이 필요한 경우가 생긴다. 간단히는 주변의 가장

kipl.tistory.com

 

'Image Recognition' 카테고리의 다른 글

Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

진자의 주기와 진폭

Physics 2021. 1. 18. 11:49
728x90

진자의 주기를 구할 때 보통 작은 진동 근사를 사용한다. 진자의 진폭이 크지 않는 경우 주기는 진폭에 무관하게 일정한 값 $T=2\pi \sqrt {\frac {\ell}{g}}$를 갖는다. 그럼 진폭이 커지는 경우는 어떻게 될까?

운동 방정식을 써도 되지만 역학적 에너지가 보존되므로 이를 이용하면(회전 관성: $I=m\ell^2$, 진폭=$\theta_0$)

$$ \frac {1}{2} I \Big(\frac {d\theta}{dt}\Big)^2 + mg \ell (1 - \cos\theta)=const= mg\ell (1- \cos \theta_0) \\  \rightarrow \quad \Big(\frac {d\theta}{dt} \Big)^2  =\frac {2g}{\ell} (\cos \theta- \cos \theta_0).$$

우변을 $\theta_0, ~\theta$에 대해서 전개하면

$$ \Big( \frac { d\theta}{dt } \Big)^2   = \frac {g}{\ell}\Big(\theta_0^2 -\frac {1}{12} \theta_0^4 - \theta^2 + \frac {1}{12} \theta^4+...\Big) =\frac{g}{\ell}(\theta_0^2 -\theta^2) \Big( 1 - \frac{1}{12} (\theta_0^2 + \theta^2)+...\Big)$$로 써지는데 작은 각 근사를 벗어났을 때 가장 큰 기여를 하는 $-(\theta_0^2 + \theta^2 ) /12$항이  음의 기여를 한다. 이는 같은 위치에서 작은 각 근사를 할 때보다 각속도가 더 작아짐을 의미한다. 따라서 진자가 더 느리게 움직여서 주기가 길어질 것이라는 예측을 구체적인 계산 없이도 할 수 있게 된다.

 

이제 주기를 구해보자. 에너지 보존식에서 변수 분리를 해서 적분하면

$$T = \int dt = 4 \sqrt {\frac {\ell}{2g}} \int_0^{\theta_0} {\frac {d\theta}{\sqrt {\cos \theta - \cos \theta_0}}}$$을 얻는다. 여기서 $\sin(\theta/2) = \sin (\theta_0/2) \sin(\varphi )$로 치환을 하면

$$T = 4\sqrt { \frac { \ell }{g}} \int_0^{\pi/2} {\frac {d \varphi}{\sqrt {1 - k^2 \sin^2 \varphi}}}, \quad k^2 = \sin^2(\theta_0/2).$$

진폭이 작은 경우($\theta_0  \ll 1 ~\rightarrow k=0)$는 적분 값이 $\frac {\pi}{2}$이므로 $T = 2\pi \sqrt {\frac {\ell}{g} }$가 됨을 확인할 수 있다.  위 적분은 타원 적분이라고 부르고 $k$가 주어지면 수치 연산을 통해서 그 값을 얻을 수 있다. 

 

좀 더 직관적으로 진폭에 따른 주기의 변화를 보기 위해서 (진자의 경우 $k^2 \le \frac {1}{2}$이므로) 급수 전개를 하면, 

$$\frac {1}{\sqrt {1-k^2 \sin^2\varphi}}   = 1 +\frac {1}{2} k^2\sin^2 \varphi +\frac {1}{2}\frac {3}{2} k^4 \sin^4 \varphi +\dots $$

이므로 주기는 

$$T = 2\pi \sqrt { \frac {\ell}{g} } \left [ 1 + \Big( \frac {1}{2} \Big)^2 k^2 + \Big( \frac {1}{2} \frac {3}{4} \Big)^2 k^4 + \dots \right]$$

로 표현된다. 이 식은 진자의 진폭이 커지면 주기도 길어진다는 것을 명확히 보여준다.

 

강의동영상을 볼 수 있는 곳:

youtu.be/34zcw_nNFGU

 

'Physics' 카테고리의 다른 글

바닥에 먼저 닿는 물체는?  (0) 2021.01.21
마찰력은 도움이 될까?  (0) 2021.01.21
진자의 주기와 진폭  (0) 2021.01.18
물이 새는 두레박이 달린 진자의 주기는  (0) 2021.01.17
왜 공은 움직이지 않을까?  (0) 2021.01.17
두레박이 달린 진자  (0) 2021.01.17
Posted by helloktk

댓글을 달아 주세요

728x90

물이 담긴 두레박 진자가 있다. 두레박에 구멍을 뚫어 물이 새게 하면 진자의 주기는

1. 변함없다.

2. 줄어든다.

3. 늘어난다.

참고로 진자의 주기는 질량에는 무관하다. 그렇지만...

두레박에 페인트를 담아 그린 패턴

설명 동영상:

youtu.be/wGS97z0WFmU

'Physics' 카테고리의 다른 글

마찰력은 도움이 될까?  (0) 2021.01.21
진자의 주기와 진폭  (0) 2021.01.18
물이 새는 두레박이 달린 진자의 주기는  (0) 2021.01.17
왜 공은 움직이지 않을까?  (0) 2021.01.17
두레박이 달린 진자  (0) 2021.01.17
왜 이런 일은 안 일어나는가?  (0) 2021.01.17
Posted by helloktk

댓글을 달아 주세요

728x90

경사진 벽에 그림처럼 공을 바싹 붙여서 놓으면 바닥과 벽은 수직방향으로 힘을 작용할 것이다. 물론 공은 중력도 받는다. 경사진 벽이 주는 힘은 수평 성분과 수직 성분이 있을 것이므로 공은 가만히 있지 못하고 오른쪽으로 가야 할 것 같은데 그렇지 않다. 바닥을 아주 미끄럽게 한 후 관찰해도 마찬가지이므로 마찰 때문은 아닐 것이다. 왜 그런가?

'Physics' 카테고리의 다른 글

진자의 주기와 진폭  (0) 2021.01.18
물이 새는 두레박이 달린 진자의 주기는  (0) 2021.01.17
왜 공은 움직이지 않을까?  (0) 2021.01.17
두레박이 달린 진자  (0) 2021.01.17
왜 이런 일은 안 일어나는가?  (0) 2021.01.17
불타는 양초 시소  (3) 2021.01.17
Posted by helloktk

댓글을 달아 주세요

두레박이 달린 진자

Physics 2021. 1. 17. 17:57
728x90

줄에 매달린 두레박에 물을 넣고 조금 흔들면 일정한 주기로 진동을 한다. 이 두레박의 물을 얼리면 진동의 주기는 어떻게 변할까?

1. 변함없다

2. 증가한다.

3. 감소한다.

'Physics' 카테고리의 다른 글

물이 새는 두레박이 달린 진자의 주기는  (0) 2021.01.17
왜 공은 움직이지 않을까?  (0) 2021.01.17
두레박이 달린 진자  (0) 2021.01.17
왜 이런 일은 안 일어나는가?  (0) 2021.01.17
불타는 양초 시소  (3) 2021.01.17
용수철 저울의 눈금은?  (0) 2021.01.17
Posted by helloktk

댓글을 달아 주세요

728x90

지구가 물체에 작용하는 중력의 세기는 지구의 중심(무게중심)에서 거리의 제곱에 반비례하게 작용한다. 컵 속에 각설탕을 넣는다고 하자. 각설탕은 지구 중력과 컵 중력을 받고 아래로 내려가지만 컵의 무게중심에 가까워질수록 거리가 작아지므로 컵이 작용하는 중력이 매우 커지게 된다. 각설탕이 어찌어찌해서 컵의 무게중심에 매우 가까운 아래쪽 적당한 위치에 도달하는 경우(B:경우) 위쪽으로 작용하는 컵의 중력이 아래쪽으로 작용하는 지구의 중력을 상쇄할 수 있다. 이 경우 각설탕은 무게중심 약간 아래쪽에서 떠있어야 하는데 이런 현상은 누구도 본 적이 없다. 왜 그럴까?

'Physics' 카테고리의 다른 글

왜 공은 움직이지 않을까?  (0) 2021.01.17
두레박이 달린 진자  (0) 2021.01.17
왜 이런 일은 안 일어나는가?  (0) 2021.01.17
불타는 양초 시소  (3) 2021.01.17
용수철 저울의 눈금은?  (0) 2021.01.17
태양계 행성의 hodograph  (0) 2021.01.16
Posted by helloktk

댓글을 달아 주세요

불타는 양초 시소

Physics 2021. 1. 17. 16:12
728x90

양초로 시소를 만들고 살짝 흔들었을 때 진동을 한다(시소 축을 양초의 중심축에서 약간 비켜서 통과하게 만들어야 한다. 왜?). 양초의 양쪽에 불을 붙여서 태우면 진동의 주기는 어떻게 변할까?

1. 그대로

2. 증가한다.

3. 짧아진다.

'Physics' 카테고리의 다른 글

두레박이 달린 진자  (0) 2021.01.17
왜 이런 일은 안 일어나는가?  (0) 2021.01.17
불타는 양초 시소  (3) 2021.01.17
용수철 저울의 눈금은?  (0) 2021.01.17
태양계 행성의 hodograph  (0) 2021.01.16
태양계의 행성 운동에서 각운동량 보존법칙  (0) 2021.01.15
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2021.01.17 19:47 신고  댓글주소  수정/삭제  댓글쓰기

    주기는 짧아질 것 같은데,
    왜 중심축을 비켜야 할 지는 모르겠네요

728x90

두 바스켓이 움직이지 않으면(예를 들면, 양쪽 줄을 서로 묶어서) 용수철 저울의 눈금은 3kg으로 나올 것이다(도르래, 줄의 무게 무시). 두 바스켓이 움직임을 시작하면 용수철 저울의 눈금은 

1. 3kg 

2. 3kg 보다 크다.

3. 3kg 보다 작다.

 

Posted by helloktk

댓글을 달아 주세요

728x90

Distance Transform은 이미지 상의 각 픽셀에서 가장 가까이 놓인 물체(전경)의 경계까지 거리를 구하는 연산이다. 보통의 경우에 이진 영상에만 적용이 된다. 구해진 거리의 정보를 이용해서 영상에 놓인 전경 물체의 skeleton을 구하거나 (OCR), 로봇의 팔 등을 자동 제어하는 데 있어서 움직이는 경로에 놓인 방해물의 위치를 예측할 수 있으며, 또한 물체까지의 경로 찾기(path finding) 등에 이용이 된다. distance transform은 사용되는 거리 측도(measure)에 따라 그 특성이 달라진다. 우선 가장 단순한 Euclidean 측도를 이용하도록 하자.

 

distance transform은 전경의 경계선을 모두 찾아서 배경의 각 픽셀에서 거리의 최솟값을 할당하면 된다. 그러나 경계가 많은 픽셀로 구성이 되는 이미지에서는 연산 비용이 매우 커지게 된다. 여기서는 이미지의 크기에 선형으로 비례하는 알고리즘에 대해서 살펴보도록 하겠다. 픽셀 수에 선형이기 위해서는 알고리즘이 raster 스캔 방식을 따라야 한다.

 

Euclidean 측도를 사용하는 distance transform은 한 점 $(x, y)$에서 전경 픽셀 $(i, j)$까지 거리의 제곱이 $z = (x-i)^2 + (y-j)^2$로 표현된다. 이미지 상의 각 전경 픽셀 꼭짓점으로 하는 paraboloid를 생각하면 distance transform은 한 픽셀에서 가장 밑에 위치한 paraboloid를 찾는 문제로 환원된다. 그러나 이차원 영역에서 이를 찾기는 쉽지 않은 작업이므로 이를 수평과 수직 방향으로 분리가 가능한지를 찾아보아야 한다. 우선 $x$-방향으로 스캔하면서 $x$-축으로 사영된 거리를 저장하면서 다시 $y$-축 방향으로 스캔하면서 저장된 $x$-방향의 거리의 제곱의 정보를 이용하면 된다. 이것은  $x$-방향으로 스캔할 때 거리의 제곱 값을 이미지의 그레이 값으로 저장하면 보다 쉽게 해결이 된다. 그러면 $y$-방향으로 스캔할 때는 꼭짓점이 이전에 $x$-방향 스캔 시 저장된 거리의 제곱만큼 이동이 된 이차곡선들 중에서 맨 아래에 놓인 것을 찾으면 된다.  따라서 2차원 distance transform은 1차원의 distance transform을 행방향 그리고 열방향으로 2번 시행을 하면 얻을 수 있다.

 

1차원의 distance transform은 일직선에 놓인 각 픽셀 위치를 꼭짓점으로 하는 이차 곡선들 중에서 가장 아래에 놓인 곡선을 선택하여 그 값을 취하면 된다. 거리의 최솟값을 구할 때 배경에 해당되는 픽셀은 선택되지 않기 위해서 배경 위치에서 이차곡선에 매우 큰 (수학적으로는 무한대) offset을 준다:

$$\text {1-dim distance transform}(p)=\underset {q\in\text {pixels}}{\text {min}}  \left( (q-p)^2+ F(q) \right)$$

초기 offset $F(q)$ 값은 전경 픽셀에서는 0, 배경 픽셀에서는 +INFINITY로 설정한다.  전경의 경계에서 형성이 된 2차 곡선 가장 낮게 놓이므로 그것에 의해서 거리가 결정이 된다.

일반적인 F(q)의 값이 주어지는 경우에는 문제는 각 픽셀점에서 보는 가장 낮게 놓인 이차곡선이 어느 것이고, 그것이 성립이 되는 영역을 찾는 문제로 귀결이 된다. 위치 $p$, $q$에 꼭짓점이 있는 이차곡선의 교차점이

$$s=\frac { (F(p) + p^2)-(F(q) + q^2) }{ 2(p-q)}$$ 표현되므로 쉽게 가장 낮은 2차 곡선의 영역을 찾을 수 있다.

 

아래 그림은 일차원 이미지의 distance transform의 결과이다: 붉은 선은 이미지(0:전경, 255:배경)이고, 녹색 선은 거리 값을 나타내는데 직관적으로 맞음을 확인할 수 있다.

네이버 블로그 이전

이 일차원 distance transform을 각각의 행에 대해서 한 후에, 그 결과를 다시 각각의 열에 또다시 적용하면 원하는 결과를 얻는다. 아래의 그림은 글씨에 대해서 skeleton을 구하기 위해서 distance transform을 적용하였다(skeleton을 얻기 위해서 반대로 함: 흰색=배경, 검은색=전경)

** 참고: en.wikipedia.org/wiki/Distance_transform

** cs.brown.edu/people/pfelzens/dt/에서 관련 논문과 원 코드를 찾을 수 있다;

#define INF     (1E+20)
#define SQR(a)	((a) * (a))
/* distanceTransform of 1d function using the Euclidean distance */
struct para_envelop {
    int x;		   // x-coord of the vertex of parabolar;
    double start;  // 현재 parabola가 가장 아래 높인 구간의 시작;
};
double para_intersection(int q, int p, double *dblimg) {
    return ((dblimg[q] + SQR(q)) - (dblimg[p] + SQR(p))) / (2. * q - 2. * p);
}
void DistanceTransform_1D(double *dblimg, int n, double* dist) {
    para_envelop *env = new para_envelop [n + 1] ;
    int k = 0;				// index of right_most parabola lower envelop;
    env[k].x = 0; env[k].start = -INF;
    env[k + 1].start = +INF;
    for (int q = 1; q < n; ++q) {
        double s = para_intersection(q, env[k].x, dblimg);
        // q에 꼭지점이 있는 parabola가 가장 밑에 있는 영역의 시작 위치(start)를 찾음;
        while (s <= env[k].start) s = para_intersection(q, env[--k].x, dblimg);
        env[++k].x = q ;    // 가장 아래에 있는 parabola 갱신;
        env[k].start = s ;
        env[k + 1].start = +INF;
    }
    k = 0;
    for (int q = 0; q < n; ++q) {
        while (env[k + 1].start < q) k++;//현재 위치에서 가장 낮은 parabola를 찾음;
        // Euclidean 거리 정보를 기록;
        dist[q] = SQR(q - env[k].x) + dblimg[env[k].x];
    }
    delete[] env;
}

'Image Recognition' 카테고리의 다른 글

Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

지상에서 포물선 운동을 하는 물체의 속도 벡터를  모두 모아 시작을 같게 만들면 속도 벡터의 끝이 그리는 자취(hodograph라 함)는 직선이 된다. 이는 속도의 차이가 가속도이고 지상에서 중력가속도는 크기와 방향이 일정하기 때문이다.

$$\vec {v}(t) =\vec {v}_0 - g\hat {j} t$$

태양계에서 행성은 태양을 한 초점으로 하는 타원 궤도상에서 움직인다. 그러면 행성의 각 지점에서 속도 벡터의 시작점을 한 지점에 모을 때 벡터의 끝점이 그리는 궤적은 무엇일까? 원 궤도를 그리면 속력이 일정하므로 당연히 원이 될 것으로 예상할 수 있지만,  타원 궤도에서는 속력은 에너지 보존을 고려하면 일정할 수 없다. 그런데 타원궤도를 그리는 경우에도 속도의 hodograph는 원궤도에서와 마찬가지로 벡터 공간에서 원으로 표현된다. 왜 그럴까?

행성의 운동 방정식에서 가속도는 

$$ \frac {d\vec {v}}{dt}= -\frac {GM}{r^2 }\hat {r} \quad \rightarrow\quad |\Delta \vec {v}| =\frac {GM}{r^2} \Delta t$$

이고, 각운동량 보존법칙에서 (각 $\varphi$는 태양을 원점으로 잼)

$$ L = |\vec{r}\times \vec {v}| = m r^2 \frac {d\varphi}{dt}  \quad\rightarrow\quad \Delta\varphi =\frac {L}{mr^2}\Delta t$$

두 식에서 $r^2$을 소거하면,

$$ |\Delta \vec{v}| = \frac {GMm}{L} \Delta \varphi$$ 

각운동량의 크기가 일정하므로 속도의 증가분이 회전각의 증가분과 비례함을 의미하고, 기하학적으로는 속도 벡터를 모아둔 벡터 공간에서 각이 증가하면서 움직이는 벡터의 끝이 반지름 $GMm/L$인 원주 위에 있을 때만 가능한 관계이다. 벡터 공간에서 벡터의 시작점이 각을 재는 원점이 아니다. hodograph의 반지름은 근일점에서 속력 $v_P$와 원일점에서 속력 $v_A$의 절반임으로 표현되는데 $v_R = (v_A + v_P)/2 = \frac {GMm}{L}$임을 근일점, 원일점에서 역학적 에너지 보존과 각운동량 보존식을  이용해서 확인할 수 있다.

설명 동영상:

youtu.be/xdIjYBtnvZU

Posted by helloktk
TAG hodograph

댓글을 달아 주세요

728x90

태양계의 행성 운동에서 각운동량은 보존이 된다(Kepler의 제2법칙). 이는 태양이 행성에 작용하는 힘인 만유인력이 중심력의 형태를 띠고 있기 때문이다. 식으로는 태양에서 행성까지 위치 벡터를 $\vec {r}$이라면 태양이 행성에 작용하는 만유인력은 $\vec {F} = \frac {GMm}{r^3}\vec {r}$로 쓸 수 있으므로 만유인력이 만드는 토크가 $\vec {\tau} = \vec {r}\times \vec {F}=0$임을 쉽게 알 수 있다. 따라서 각운동량 보존은 자명해진다.

만약 행성의 위치를 재는 원점을 태양이 아니라 다른 지점으로 잡으면 어떻게 될까? 이 경우 만유인력의 방향과 위치 벡터의 방향이 나란하지 않으므로 토크가 0이 안된다. 그럼 각운동량은 원점을 어디로 잡는가에 따라 보존되기도 하고 안되기도 하는 물리량일까? 무엇을 놓치고 있는 것일까?

Posted by helloktk

댓글을 달아 주세요

728x90

사람이 걷는 동작은 복잡하지만 몇 가지 가정을 하면 단순한 강체의 운동으로 근사를 할 수 있다. 우선 사람이 걷는 동안 항상 한 발은 땅을 딛고 있다. 그리고 땅을 딛고 있는 발을 기준으로 몸의 질량중심은 원호를 그리면서 거의 일정한 속도로 움직이게 된다. 사람을 질량중심에 모든 질량이 뭉친 점으로 근사를 하면 걷는 동작은 거의 질량이 없는 막대(다리: 길이=$L$)에 매달린 역진자(inverted pendulum) 운동과 비슷하다고 볼 수 있다. 또한 원호를 따라 중심이 이동하는 속력은 거의 일정하다고 근사할 수 있다.

이 경우에 질량중심은 등속 원운동을 한다고 볼 수 있고, 뉴턴의 제2법칙에 의해서 질량중심에서 발 쪽을 향하는 힘 성분이 구심력 역할을 한다. 구심력 역학을 할 수 있는 힘은 중력, 수직항력, 그리고 마찰력이 있다.

수직에서 $\theta$만큼 기울어졌을 때 등속원운동의 가속도 벡터는

$$ \vec{a} =\frac{V^2}{L} (-\sin \theta \hat{x} -\cos \theta \hat{y})$$

이고,  운동 방정식의 수식 성분을 보면 중력과 수직항력이 기여하므로

$$ \sum F_y = N - mg = m a_y = -\frac{m V^2}{L} \cos \theta $$

움직이는 속력이 너무 빠르면 달리는 동작이 되고, 이 경우 양발이 땅에서 떨어지게 된다. 발이 땅에서 떨어지지 않으려면  수직항력 $N\ge0$인 조건을 만족하도록 걷는 속력을 조절해야 한다. 

$$N \ge  0 \quad \rightarrow  \quad V \le  \sqrt{\frac{gL }{\cos \theta}}$$

원호를 그리며 움직이는 동안 이 관계가 성립해야 하므로 걷기가 가능한 최대속력은:

$$V_\text{max} = \sqrt{gL}$$

사람을 너무 단순화시킨 감은 있지만 왜 다리가 긴 사람이 더 빨리 걸을 수 있는지를 물리 법칙으로 추정할 수 있다는 것을 보여준다. 사람을 좀 더 복잡한 강체로 근사하더라도 $\sqrt{gL}$ 앞의 factor만 바뀔 것이다.  사람이 걷는 동안 땅을 딛지 않는 발은 약간 구부러진 상태로 엉덩이를 축으로 일종의 물리진자처럼 행동하는데 이를 이용해도 역시 같은 추정치를 얻을 수 있다.

 

다리 길이가 $L=1\text{m}$이면,  $V_\text{max}= \sqrt{(1 \text{m}) \times (9.8\text{m/s}^2)} = 3.13\text{m/s}=11.3\text{km/h}$. 경보 세계기록이 대략 $15\text{km/h}$이므로 합리적인 추정이 된다.

 

youtu.be/HypJY1XWkWY

Posted by helloktk
TAG 역학

댓글을 달아 주세요

728x90

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). 더불어 주기성 때문에 $a_k$의 계산도 절반만 하면 된다: 

$$a_{k+ N/2} = e_k - W^k o_k, \quad k = 0,1,..., N/2-1.$$

역변환의 경우에는 twiddle factor $W$를 $W^*$로 바꾸면 되므로(전체 normalization factor $\frac{1}{N}$가 덧붙여진다) 코드에서 쉽게 수정할 수 있다: $\tt W_{im}  \rightarrow  - { W}_{im}$. 아래는 FFT를 재귀적으로 구현한 코드이다. 비재귀적 구현은 kipl.tistory.com/22에서 찾을 수 있다.

void split(int N, double* data) {
    double *tmp = new double [N / 2];
    for (int i = 0; i < N / 2; i++)   //copy odd elements to tmp buffers;
        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 fft(int N, double* re, double* im) {
    if (N < 2) return;
    else {
        split(N, re);                             // divide;
        split(N, im);        
        fft(N / 2, &re[    0], &im[    0]);       // conquer;
        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(2. * PI * k / N);   // real part of twiddle factor:W
            double Wim =  -sin(2. * PI * k / N);  // imag part of twiddle factor:W
            double WOre = Wre * Ore - Wim * Oim;  // WO = W * O
            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;
}

 6개의 주파수$\{f=2,5,9,11,21,29\text{Hz}\}$가 섞인 신호에서 [0,1]초 사이에 일정한 간격으로 sampling한  64개의 data을 FFT한 결과 (Nyquist frequency=32Hz);

$$s(t)=\sum_{i=0}^{5} (i+1) \cos(2\pi f_i t) $$

int fft_test_main() {
    const double PI = 4. * atan(1.);
    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 ( 2.0 * PI * 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;
}

FFT결과:

real part (imaginary part = 0)

'Image Recognition' 카테고리의 다른 글

Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

Edge를 보존하는 low-pass filter로 edge의 방향으로 픽셀 평균값으로 대체한다. 방향은 8방향 (45도) 기준이다. 에지의 방향은 그림처럼 중심 픽셀을 기준으로 주변의 8개를 픽셀을 이용해서 45도씩 회전하면서  차이(평균 변화율)를 계산하여 가장 작은 값을 주는 방향이 된다. 중심 픽셀이 line[1][1]일 때

     0도 에지: abs(line[1][0] - line[1][2])가 최소

   45도 에지: abs(line[0][0] - line[2][2])가 최소

   90도 에지: abs(|ine[0][1] - line[2][1])가 최소

 135도 에지: abs(line[0][2] - line[2][0])가 최소 

edge 방향이 정해졌으면 에지 방향의 3 픽셀의 평균값으로 중심 픽셀 값을 대체하면 에지는 보존하면서 평균 필터를 적용한 효과를 얻을 수 있다.

 

void smooth_ep_3x3(BYTE *src, int width, int height, BYTE* dst) {
    const int xend = width - 2, yend = height - 2;
    BYTE *line[3];
    line[0] = src; line[1] = line[0] + width; line[2] = line[1] + width;
    BYTE *dst_line = dst + width;        // starting dst row;
    for (int y = 0; y < yend; ++y) {
        for (int x = 0; x < xend; ++x) {
            int diff1 = line[1][x] - line[1][x + 2];
            if (diff1 < 0) diff1 = - diff1;
            int diff2 = line[0][x] - line[2][x + 2];
            if (diff2 < 0) diff2 = - diff2;
            int diff3 = line[0][x + 1] - line[2][x + 1];
            if (diff3 < 0) diff3 = - diff3;
            int diff4 = line[0][x + 2] - line[2][x];
            if (diff4 < 0) diff4 = - diff4;
            
            if (diff1 <= diff2 && diff1 <= diff3 && diff1 <= diff4)         //0-도
                dst_line[x + 1] = (line[1][x] + line[1][x + 1] + line[1][x + 2]) / 3;
            else if (diff2 <= diff3 && diff2 <= diff4)                      //45-도
                dst_line[x + 1] = (line[0][x] + line[1][x + 1] + line[2][x + 2]) / 3;
            else if (diff3 <= diff4)                                        //90-도;
                dst_line[x + 1] = (line[0][x + 1] + line[1][x + 1] + line[2][x + 1]) / 3;
            else                                                            //135-도;
                dst_line[x + 1] = (line[0][x + 2] + line[1][x + 1] + line[2][x]) / 3;
                
            dst_line += width;                              //move to next dst line;
        }
        // increases src line ptr;
        BYTE *tptr = line[2] + width;
        line[0] = line[1]; line[1] = line[2]; line[2] = tptr;
    }
};

9번 반복 적용 결과: 지문처럼 규칙적인 패턴이 나오는 경우는 Gabor 필터를 이용하면 보다 좋은 결과를 얻을 수 있다.

네이버 블로그 이전

iteration에 따른 correlation의 계산:

iter=0, corr w.r.t original=0.925144, corr w.r.t previous=0.925144
iter=1, corr w.r.t original=0.903661, corr w.r.t previous=0.985120
iter=2, corr w.r.t original=0.888224, corr w.r.t previous=0.994821
iter=3, corr w.r.t original=0.876582, corr w.r.t previous=0.997403
iter=4, corr w.r.t original=0.866254, corr w.r.t previous=0.998183
iter=5, corr w.r.t original=0.856620, corr w.r.t previous=0.998596
iter=6, corr w.r.t original=0.847365, corr w.r.t previous=0.998841
iter=7, corr w.r.t original=0.838400, corr w.r.t previous=0.998981
iter=8, corr w.r.t original=0.829703, corr w.r.t previous=0.999088

'Image Recognition' 카테고리의 다른 글

Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Posted by helloktk

댓글을 달아 주세요

728x90

컬러를 양자화하는 방법으로 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 노드의 수를 줄이는 과정을 시행한다. 이 병합 과정은 원래의 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 값과 가장 가까운 유클리디안 거리를 주는 팔레트의 인덱스를 할당하는 방법도 있다. 이 방법을 이용하면 아래의 결과와 약간 차이가 생긴다.

양자화된 컬러 이미지
L2 error;

** 네이버 블로그에서 이전;

source code(C++): web.archive.org/web/20050306011057/www.drmay.net/octree/

'Image Recognition' 카테고리의 다른 글

FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Posted by helloktk

댓글을 달아 주세요

728x90

컬러 영상을 처리할 때 가장 흔히 사용하는 컬러 표현은 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공간에서의 컬러 분포가 몇 개의 대표적인 컬러(예:박스의 중앙값)로 줄어들면(양자화 과정:: 공간에 smooth 하게 분포한 것이 몇 개의 점으로 대체됨), 원본 영상의 각각의 픽셀에서의 대체 컬러 값은 원래의 컬러와 가장 유사한, 죽 RGB 공간에서 유클리디안 거리가 가장 작은 박스의 컬러로 대체하면 된다. 

 

그러나 너무 적은 수의 컬러로 줄이게 되면 인접 픽셀 간의 컬러 값의 차이가 눈에 띄게 나타나는 현상이 생기게 된다. 이러한 것을 줄이기 위해서는 디더링(dithering) 과정과 같은 후처리가 필요하게 된다.

 

**네이버 블로그에서 이전;

 

source code: dev.w3.org/Amaya/libjpeg/jquant2.c

'Image Recognition' 카테고리의 다른 글

Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Posted by helloktk

댓글을 달아 주세요

728x90

Union-Find 알고리즘을 이용하여 영역 분할 과정이다. 각각의 픽셀에서 4방향으로 연결된 픽셀이 속한 영역 merge 할 수 있는지를 시도한다. merge 조건은 현재 픽셀의 그레이 값과 인접 영역의 평균 그레이 값의 차이가 주어진 임계값보다 작은 경우다.

$$\text{merge 조건: }\quad | 그레이 값- \text{인접 영역 평균값}| < threshold$$

컬러 영상의 경우에는 RGB 채널에 조건을 부여하거나 gray level만 이용해서 판단할 수 있다. root node 수만큼의 분할 영역이 만들어지고, 임계값이 클수록 분할된 각 영역의 사이즈가 커진다.  

gray = 0.2989 * r + 0.5870 * g + 0.1140 * b ;

같은 평균 픽셀값을 가지고 있더라도 4-방향으로 서로 연결된 영역이 아니면 합병되지 못하고 서로 다른 영역으로 남는다. 이를 하나의 영역으로 만들려면 분할된 영역을 다시 검사를 하는 과정이 필요하다. 자세한 과정은 Statistical Region Merge(SRM) 알고리즘(kipl.tistory.com/98)에 잘 구현이 되어 있다.

int FindCompress ( int *p, int node ) { // find root node;
    if ( p[node] < 0 ) return node;     // root-node;
    else return p[node] = FindCompress ( p, p[node] ); // 찾음과 동시에 depth을 줄인다;      
}
// root가 r1인 집합(단일픽셀)을 root가 r2인 집합으로 합병을 시도한다;
// 합병이 되면, 전체 픽셀값의 합과 픽셀수를 update한다;
BOOL MergePixel(int r1, int r2, int thresh, int *sums, int *sizes, int *p) {
    double diff = sums[r2] / (double) sizes[r2] - sums[r1] / (double) sizes[r1];
    if ( fabs(diff) < thresh ) { // merge two sets;
        sums[r2]  += sums[r1];
        sizes[r2] += sizes[r1];
        p[r1]     = r2;
        return TRUE;
    }
    return FALSE;
}
// root 노드는 분할영역의 픽셀 갯수와 픽셀 값의 합을 저장한다.
// 처음 각 픽셀이 root 이고, 픽셀 갯수는 1, 픽셀값은 자신의 픽셀값을 저장;
BOOL UnionFindAverage ( BYTE *image, int width, int height, int thresh ){
    int npixs = width * height;
    int* p      = new int [npixs];    // parent table;
    int* sums   = new int [npixs];  
    int* sizes  = new int [npixs];
    // initial setting;
    for ( int k = 0; k < npixs; k++ ) {
        p[k]     = -1;                // 각 픽셀이 root임;
        sums[k]  = image[k];
        sizes[k] = 1;
    }
    // 4-connectivity:
    // 첫행; LEFT 픽셀이 속한 영역의 평균값과 차이가 thresh 이내이면 left로 합병함;
    for ( int x = 1; x < width; x++ ) {
        int left = FindCompress ( p, x - 1 );
        int curr = FindCompress ( p, x );
        MergePixel(curr, left, thresh, sums, sizes, p);
    }
    for ( int y = 1, pos = y * width; y < height; y++ ) {
        // x = 0; up 픽셀이 속학 영역과 합병 체크;
        int up = FindCompress ( p, pos - width );
        int curr = FindCompress ( p, pos );
        MergePixel(curr, up, thresh, sums, sizes, p);
        pos++;
        // UP-LEFT 픽셀 영역과 합병 체크;
        for ( int x = 1; x < width; x++, pos++ ) {
            int left = FindCompress ( p, pos  - 1 );
            int curr = FindCompress ( p, pos );
            // left와 합병이 되면 left가 up과 합병이 되는지 확인; 안되면 up과 합병시도;
            if (MergePixel(curr, left, thresh, sums, sizes, p)) curr = left;
            int up = FindCompress ( p, pos - width );
            if (curr != up) MergePixel(curr, up, thresh, sums, sizes, p);
        }
    }
    // 평균 이미지를 만듬(input is overwritten);
    for ( int k = 0; k < npixs; k++ ) {
        int pos = k;
        // find root node: p[root] = -1;
        while ( p[pos] >= 0 ) pos = p[pos];
        int a = ( int ) ( sums[pos] / ( double ) sizes[pos] + 0.5 );
        image[k] = a > 255 ? 255 : a;
    };
    delete [] p;
    delete [] sums;
    delete [] sizes;
    return TRUE;
};

네이버 블로그 이전
statistical region merge 알고리즘을 적용한 결과

'Image Recognition' 카테고리의 다른 글

Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Kuwahara Filter  (2) 2020.12.28
Posted by helloktk

댓글을 달아 주세요

728x90

Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.

이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram(ch)과 평균을 계산할 수 있는 1차 culmuative histogram(cxh)을 입력으로 사용한다. 

$$ th= \text {argmax} \left( \sigma^2_B = \sum_{j=0}^{nclass-1} \omega_j m_j^2 \right)$$

 

* Otsu 알고리즘을 이용한 이미지 이진화 코드: kipl.tistory.com/17

* 좋은 결과를 얻으려면 히스토그램에 적당한 필터를 적용해서 smooth하게 만드는 과정이 필요하다.

// 0 <= start < n;
double histo_partition(int nclass, double cxh[], int ch[], int n, int start, int th[]) {
    if (nclass < 1) return 0;
    if (nclass == 1) {
        int ws; double ms;
        if (start == 0) {
            ws = ch[n - 1];
            ms = cxh[n - 1] / ws;
        } else {
            ws = ch[n - 1] - ch[start - 1];             // start부터 끝까지 돗수;
            ms = (cxh[n - 1] - cxh[start - 1]) / ws;    // start부터 끝까지 평균값;
        }
        th[0] = n - 1;
        return ws * ms * ms;                            // weighted mean;
    }

    double gain_max = -1;
    int *tt = new int [nclass - 1];
    for (int j = start; j < n; j++) {
        int wj; double mj;
        if (start == 0) {
            wj = ch[j]; 
            mj = cxh[j];                    //mj = cxh[j] / wj;
        }
        else {
            wj = ch[j] - ch[start - 1];     //start부터 j까지 돗수;
            mj = (cxh[j] - cxh[start - 1]); //mj = (cxh[j] - cxh[start - 1]) / wj;
        }
        if (wj == 0) continue;
        mj /= wj;                           //start부터 j까지 평균;
        double gain = wj * mj * mj + histo_partition(nclass - 1, cxh, ch, n, j + 1, tt);
        if (gain > gain_max) {
            th[0] = j;
            for (int k = nclass - 1; k > 0; k--) th[k] = tt[k - 1];
            gain_max = gain;
        }
    }
    delete [] tt;
    return gain_max;
};

trimodal 분포의 분리;

class0: 0~th[0]

class1: (th[0]+1)~th[1],

class2: (th[1]+1)~th[2]=255

th[0]=103, th[1]=172  (th[2]=255)
th[0]=88, th[1]=176, th[2]=255

 

'Image Recognition' 카테고리의 다른 글

Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Kuwahara Filter  (2) 2020.12.28
Moving Average을 이용한 Thresholding  (0) 2020.11.26
Posted by helloktk

댓글을 달아 주세요

728x90

책상에 세워둔 연필이 넘어져서 바닥에 닿는데 걸리는 시간은 어떻게 될까? 대략 경험적으로 보면 연필 길이가 길수록 넘어지는데 시간이 오래 걸린다. 왜 그럴까? 넘어지는데 걸리는 시간은 연필의 운동을 결정하는 물리량과 관련 있는데, 우선 크기와 관련 있는 길이$(L)$가 있고, 넘어지려면 토크를 작용해야 하는데 이는 중력(가속도: $g$)과 관련이 있다.(마찰/수직 항력은 회전축에 걸리므로 무관하다). 물론 질량에도 의존할 수 있는데, 중력 가속도, 길이, 질량의 물리량으로 시간의 단위를 만들 수 있는 조합은 $\text{const} \times\sqrt { L /g}$ 밖에 없다. 즉, 넘어지는데 걸리는 시간은 길이에 제곱근에 비례한다. 

구체적으로 얼마나 걸리는지 계산을 시도해 보자. 연필이 넘어지면 수직과 이루는 각 $\theta$가 증가한다: $0\rightarrow \pi/2$. 뉴턴 방정식을 이용하면 $\theta$가 만족해야 할 미분방정식을 얻을 수 있지만, 실질적으로 일하는 힘이 중력밖에 없으므로 역학적 에너지가 보존된다는 사실을 이용하는 것이 좀 더 수월하다. 연필을 균일한 막대로 근사하면 넘어지는 과정은 연필심에 대한 회전운동이다. 역학적 에너지 보존식을 쓰면

$$\frac {1}{2} I_{tip} \Big(\frac {d\theta}{dt}\Big)^2 + \frac {1}{2} MgL \cos \theta = \text {const} = \frac {1}{2} MgL.$$ 여기서 각속도를 구하면

$$ \\ \frac{d\theta}{dt} = \sqrt { \frac {3g}{L}  (1- \cos \theta ) } = \sqrt {\frac {6g}{L}} \sin \frac {\theta}{2},$$

이고, 적분하여 수직 상태에서 바닥에 도달하는데 걸리는 시간을 구하면,

$$T= \sqrt { \frac {L}{6g}} \int_0^{\pi/2} \frac {d \theta }{\sin \frac {\theta}{2} } \rightarrow \infty.$$

기대(?)와는 다르지만 이 결과는 구체적으로 계산하지 않더라도 예상할 수 있다. 왜냐면 완전히 똑바로 서 있으면 회전의 시작에 필요한 토크를 생성하는 힘이 없기 때문이다. 연필에는 중력이 작용하고 끝에 마찰력이나 수직 항력이 있긴 하지만 토크를 만들지는 못한다. 따라서 연필을 넘어뜨리기 위해서는 처음에 약간의 충격을 주던지(속도 제공) 아니면 약간 기울인 상태에서 시작해야 한다. 처음 $\theta_0$의 각도에서 시작하였다면 걸리는 시간은 위 적분식에서

$$ T = \sqrt{\frac {2L}{3g} }\log \frac { \tan \frac {\pi}{8}}{\tan \frac {\theta_0}{4} }.$$

$\theta_0\rightarrow 0$이면 시간은 $T\sim -\sqrt { \frac {2L}{3g}} \log \theta_0 \rightarrow \infty$이고, 유한할 때는 길이의 제곱근에 비례함도 확인할 수 있다.

youtu.be/oowAdPjDY5M

 

Posted by helloktk

댓글을 달아 주세요

Center of Percussion

Physics 2021. 1. 8. 16:07
728x90

야구 배트로 공을 칠 때 공을 맞추는 지점을 잘 선택하면 배트를 잡는 위치(회전축)의 움직임이 거의 없게(즉, 손이 충격을 안 받게) 만들 수 있다. 이는 배트가 힘을 받았을 때 운동이 질량중심과 같은 속도로 움직이는 병진 운동과 질량중심에 대한 회전운동의 합으로 표현되는 점을 고려하면 쉽게 이해할 수 있다. 배트의 각 지점의 실제 속도는 질량중심의 속도와 회전운동에 의한 속도의 벡터 합이므로 손잡이 부분에서 0인 조건을 만족되면 가능하게 된다. 이 경우 배트의 운동은 손잡이 위치에 대한 순수 회전운동으로 표현할 수 있다. 그리고 공을 맞추는 지점을 center of percussion(COP)이라고 한다(손잡이 위치에 대한 상대적 위치임)

공이 배트에 준 힘을 $F$라면 (손이 준 힘이 없는 조건에서)

$$\text {cm-병진}:  ~F = M a,$$

$$\text {cm-회전}: ~Fb = I_{cm} \alpha. $$

손잡이 지점의 가속도가 없는 조건을 쓰면

$$ a_{손잡이}= a- \alpha p = 0\quad \rightarrow ~ \frac {F}{M} = \frac {Fb} {I_{cm}} p \quad \therefore b =  \frac {I_{cm}}{pM}.$$

$b$와 $p$는 정확히 대칭적 역할을 한다. COP 지점을 손으로 잡고 원래 손잡이 위치에 공을 맞추어도 손에 충격이 오지 않는다.

회전관성을 직접 구하지 않고 어떻게 COP을 알아낼 수 있을까? 손잡이나 그것의 COP 지점을 회전축으로 해서 배트를 진동시키면 일정한 주기를 갖고 진동을 한다(물리 진자). 이때 진동의 주기는

$$ T ={2\pi} \sqrt{ \frac{I_{cm} + M p^2 }{ Mgp}} = 2\pi \sqrt{\frac{b+p}{g}} .$$  

이 값은 길이가 $\ell=b+p$인 단순진자의 주기와 같다. 따라서 단순 진자의 길이를 바꾸면서 배트의 주기와 비교해서 같아지는 경우를 찾으면 된다(배트의 질량중심은 균형을 이용하면 쉽게 찾을 수 있다)

https://youtu.be/Dw3UpKQVhVY

 

Posted by helloktk

댓글을 달아 주세요

728x90

동일한 비커에 같은 양의 물을 채운 후 한쪽 비커에는 바닥에 연결된 줄로 탁구공을 완전히 잠기도록 만들었고, 다른 비커에는 스탠드와 줄을 이용해서 쇠공이 바닥에 닿지 않고 물속에 완전히 잠겨 있게 만들었다. 이 두 비커를 그대로 양팔 저울 위에 올리면 어느 쪽으로 기울까? 단, 탁구공과 쇠공의 부피는 같다. (이 문제는 googling을 하면 찾아볼 수 있을 정도로 유명하다)

1. 수평을 유지: 탁구공과 쇠공이 밀어낸 물의 부피가 같아 비커 바닥이 받는 압력도 동일하게 증가하므로

2. 탁구공 쪽: 쇠공은 스탠드가 지탱하여 영향이 없지만 탁구공의 무게는 결국 바닥이 받기 때문에

3. 쇠공 쪽: 부력은 같지만 탁구공에 연결된 줄은 위로 당겨 바닥을 누르는 힘이 감소시키므로

Posted by helloktk

댓글을 달아 주세요