Geometry & Recognition

먼 곳에서 $v$의 속력으로  달에 접근하는 우주선이 달(속력: $U$)의 중력에 의해서 그림과 같은 궤도를 그리면서 움직인다(궤도는 이심률이 매우 큰 타원으로 간주한다). 다시 우주선이 달에서 충분히 멀어졌을 때 속력은 얼마쯤 될까? 단, 달의 질량은 우주선에 비해 매우 크다고 생각할 수 있다.

1. $\approx v$

2. $\approx v+U$

3. $\approx v+2U$

우주선의 속력이 증가했다면 누가 에너지를 제공한 것인가?

https://www.youtube.com/watch?v=E01J_uAPM64 

마찰이 없는 레일 위에 설치된 대포에서 포탄을 발사한다. 포신이 수평에서 30도 윗방향을 향할 때 포탄의 날아가는 실제 각은?(지상에서 보는 기준)

1. 30도

2. 대포가 밀리므로 30도 보다 작다.

3. 대포가 밀리므로 30도 보다 크다.

모래가 섞인 얼음을 물이 담긴 컵에 넣었더니 바닥에 가라 앉았다. 얼음이 다 녹으면 수면의 높이는?

1. 올라간다.
2. 변함없다.
3. 내려간다.
4. 모래의 양에 따라 다르다.

참고: https://kipl.tistory.com/120

동물이 달릴 때 다리와 바닥 사이의 마찰력으로 가속을 하거나 감속할 수 있다. 개가 얼마나 가속/감속을 할 수 있는가는 사냥개로써의 특성에 중요한 요소 중 하나이다. 그럼 어떤 신체적인 요건이 이를 결정하는 알아보기 위해 우선 개를 그림과 같이 간단히 강체로 근사를 하여 가능한 가속도를 구해보자.

개의 가속도가 $a_x$일 때, 뉴턴의 운동 방정식을 세우면

$$ \sum F_x = f_f + f_h = m a_x ,$$

$$ \sum F_y = N_f + N_h -m g = 0$$

그리고 넘어지지 않고 안정적으로 가속하기 위해서는 회전하지(넘어지지) 않아야 한다. 질량중심에 대한 회전운동방정식을 적으면

$$ \sum \tau = N_f L_f - N_h L_h + (f_f + f_h) L_\text{leg} = 0$$

가속하는 동안 개의 발이 땅에 접촉을 하기 있어야 하므로 앞 뒤 다리에 걸리는 수직항력은 $N_f \ge  0$, $N_h \ge 0$이어야 한다. 위 운동방정식에서 앞 다리와 뒷다리의 수직항력이 각각

$$ N_f = \frac{m L_\text{leg}}{L_f + L_h} \Big( g \frac{L_h}{L_\text{leg}} - a_x \Big),$$

$$ N_h = \frac{m L_\text{leg}}{L_f + L_h} \Big( g \frac{L_f}{L_\text{leg}} + a_x \Big)$$

이다. 가속하는 동안에는 항상 $N_h>0$이므로 $N_f \ge 0$ 조건에서 가능한 가속도는 

$$ a_x \le  g\frac{L_h}{L_\text{leg}}$$

이고, 또 감속($a_x<0$)하는 동안에는 항상 $N_f >0$이므로  $N_h\ge 0$인 조건에서 가능한 가속도는

$$ |a_x| \le  g\frac{L_f}{L_\text{leg}}$$

을 만족해야 한다. 즉, 가속/감속을 잘하는 사냥개의 조건은 허리가 길고 다리가 짧아야 한다.

참고: https://kipl.tistory.com/265

주어진 점집합을 원으로 피팅하기 위해 이차식

$$A(x^2 + y^2) + Bx + Cy +D =0$$

을 이용하자. 원의 경우는 $A=0$인 경우는 직선을 나타내고, $A\ne0$인 경우가 원을 표현한다. 물론 $A=1$로 설정을 할 수 있으나 이 경우는 주어진 데이터가 원의 일부 부분만 나타내는 경우 문제가 생긴다. 원은 중심과 반지름을 알면 결정되므로 3개의 변수가 필요한데 위의 이차 다항식은 4개의 변수를 가진다. 원을 제대로 표현하기 위해서는 제한 조건이 들어오는데 

$$ \Big(x + \frac {B}{2A} \Big)^2 + \Big(y + \frac {C}{2A} \Big)^2 = \frac {B^2 +C^2 -4A  D}{4A^2}$$

으로 쓸 수 있으므로 좌변이 양수이므로 

$$ B^2+C^2 - 4A  D>0$$이어야 한다. 따라서 $$ B^2 + C^2 - 4A  D=1$$인 제한 조건을 둘 수 있다.

주어진 각 점들이 정확히 추정하는 원 상에 있지 않으므로 이차식의 좌변의 값이 반드시 0이 되지는 않는다. 잘 피팅하기 위해서는 제한조건을 만족하면 각 점들에서 차이를 최소로 만드는 파라미터를 찾으면 된다:

$$ \epsilon(A,B, C, D)= \sum | A(x_i ^2 + y_i^2) + Bx_i + Cy_i +D |^2 - \lambda (B^2 +C^2 - 4A D-1) \longrightarrow \text {min}$$

여기서 $\lambda$는 Lagrange multiplier이다. 식을 행렬로 표현하기 위해

$$z_i= x_i^2 + y_i^2 , ~M_{zz}= \frac {1}{n} \sum z_i^2 , ~M_{xx} = \frac{1}{n} \sum x_i^2, ~M_{yy}=\frac{1}{n}\sum y_i^2,  $$

$$M_{zx}= \frac {1}{n} \sum z_i x_i  , ~M_{zy} = \frac{1}{n} \sum z_i y_i , ~M_{xy}=\frac{1}{n}\sum x_i y_i, $$

$$M_{z}= M_{xx}+M_{yy} , ~M_{x} = \frac{1}{n} \sum x_i, ~M_{y}=\frac {1}{n}\sum y_i,  $$

$${\tt x}= \left(\begin {array}{c} A \\ B \\ C \\D\end {array}\right) , ~ \tt M = \left( \begin{array}{cccc} M_{zz}& M_{zx}& M_{zy} & M_z \\ M_{zx} & M_{xx}& M_{xy} & M_{x}\\ M_{zy} & M_{xy} & M_{yy} &M_y \\M_z &M_x &M_y& 1 \end{array} \right), ~\tt  N=\left( \begin{array}{cccc} 0& 0& 0& -2 \\ 0& 1& 0 & 0\\ 0& 0& 1& 0\\ -2&0&0& 0 \end{array} \right)$$

로 놓으면 ($\tt M=\text {scattering matrix}$) $$ n^{-1} \epsilon= {\tt x}^T. {\tt M}. {\tt x} - \lambda ( {\tt x}^T. {\tt N}. {\tt x} - 1)$$

로 표현된다. 질량중심 좌표계에서 계산을 하면 $M_x= M_y =0$이므로 행렬이 더 간단해진다.

$\epsilon$을 $\tt x^T$에 대해서 미분하면

$$ \tt M.x = \lambda \tt N.x$$

인 lagrange multiplier가 고윳값이 되는 일반화된 고윳값 방정식을 얻는다. 이 고윳값 식이 해를 가지려면

$$\text {det}( \tt M -\lambda \tt N)=0$$

이어야 한다. 이 특성 방정식은 4차 다항식이므로 closed form 형태의 해를 쓸 수 있으나 실질적으로 수치해석적으로 구하는 것이 더 편하다(Newton-Raphson). $\text {min}(\epsilon)$은 제한조건 때문에

$$ n^{-1}\epsilon = \lambda \tt x^T. N. x = \lambda > 0$$

이므로 $0$보다 큰 최소 고윳값이 원하는 해이다. ($\lambda = 0$은 $\text {det}(\tt M)=0$, 즉, $\tt M$이 singular 한 경우)

그런데 주어진 고윳값 $\lambda$ 해당하는 고유 벡터가 $\tt x$일 때, 임의의 0이 아닌 상수 $\mu$에 대해 $\mu \tt x$도 역시 고유벡터이다. 고유 벡터의 크기를 고정하기 위해서 제한조건을 적용하면 $1=\mu^2 \tt x^T .N .x = \mu^2 \tt x^T. M. x/\lambda$이므로 $\mu = \sqrt{\lambda/ \tt x^T. M. x}$로 선택하면 된다. 그렇지만 해가 원을 나타내는 경우 원의 중심 $(-B/2A, -C/2A)$, 반지름 $\sqrt{ (B/2A)^2 + (C/2A)^2 - D/A}$은 계수의 비에만 의존하므로 굳이 normalized 된 고유벡터를 구할 필요가 없다. 따라서, 해 중에 직선인 경우는 무시하고($A=0$) 원만을 선택할 때는 $A=1$로 놓고 계산을 해도 영향을 받지 않는다. 또한 이 경우 고유벡터 성분을 구체적으로 적을 수 있어 수치해석에 의존할 필요도 없어진다.

최소의 고윳값에 대해서 고유 방정식을 풀면 고유벡터의 성분은

$$ C_{xy} \equiv  M_{xx} M_{yy}- M_{xy}^2 , ~~\Delta \equiv \lambda^2 - M_z   \lambda +C_{xy}$$

$$A=1$$

$$B = \frac {M_{zx}(M_{yy} - \lambda)-M_{zy} M_{xy}}{ \Delta }$$

$$C = \frac {M_{zy}(M_{xx} - \lambda)-M_{zx} M_{xy}}{ \Delta }$$

$$D=-M_z -2\lambda$$로 주어진다. 따라서 원의 중심은

$$ c_x = -\frac {B}{2A}= \frac{    M_{zx}(M_{yy} -\lambda) - M_{zy}M_{xy}}{2  \Delta}, \quad c_y= -\frac{C}{2A}= \frac{ M_{zy}(M_{xx} -\lambda )- M_{zx} M_{xy}}{2  \Delta} $$

원의 반지름은 

$$ r =\sqrt { (B/2A)^2+(C/2A)^2 -D/A} = \sqrt{ c_x^2 + c_y^2 +M_{z} +2\lambda}$$ 

로 주어진다.

** 추가: Sylvester's law of inertia에 의하면 위에서 구하는 고유값 $\lambda$의 부호 개수는  $\tt N$의 고유값 부호별 개수와 같아야 하는데, $\tt N$의 고유값이 $\{-2,1,1,2\}$이므로 양수인 고유값 $\lambda$는 3개가 존재한다.

Ref: V. Pratt, Direct least-squares fitting of algebraic surfaces, Computer Graphics 21, 1987, 145–152.

double circleFit(std::vector<CPoint>& data, double& centerx, double& centery, double& radius) {
    const int maxIter = 99;
    double mx = 0, my = 0;
    for (int i = data.size(); i-- > 0;)
        mx += data[i].x, my += data[i].y;

    // center of mass;
    mx /= data.size(); my /= data.size();

    // moment calculation;
    double Mxy, Mxx, Myy, Mzx, Mzy, Mzz;
    Mxx = Myy = Mxy = Mzx = Mzy = Mzz = 0.;
    for (int i = data.size(); i-- > 0;) {
        double xi = data[i].x - mx;   //  center of mass coordinate
        double yi = data[i].y - my;   //  center of mass coordinate
        double zi = xi * xi + yi * yi;
        Mxy += xi * yi; Mxx += xi * xi;
        Myy += yi * yi; Mzx += xi * zi;
        Mzy += yi * zi; Mzz += zi * zi;
    }
    Mxx /= data.size(); Myy /= data.size();
    Mxy /= data.size(); Mzx /= data.size();
    Mzy /= data.size(); Mzz /= data.size();
    double Mz = Mxx + Myy;
    double Cxy = Mxx * Myy - Mxy * Mxy;
    double Vz = Mzz - Mz * Mz;
    // coefficients of characteristic polynomial;
    double C2 = 4 * Cxy - 3 * Mz * Mz - Mzz; 
    double C1 = Vz * Mz + 4 * Cxy * Mz - Mzx * Mzx - Mzy * Mzy;
    double C0 = Mzx * (Mzx * Myy - Mzy * Mxy) + Mzy * (Mzy * Mxx - Mzx * Mxy) - Vz * Cxy;
    // Newton's method starting at lambda = 0
    double lambda = 0;
    double y = C0;  // det(lambda = 0)
    for (int iter = 0; iter < maxIter; iter++) {
        double Dy = C1 + lambda * (2. * C2 + 16.*lambda * lambda);
        double lambdaNew = lambda - y / Dy;
        if ((lambdaNew == lambda)) break;
        double ynew = C0 + lambdaNew * (C1 + lambdaNew * (C2 + 4 * lambdaNew * lambdaNew));
        if (fabs(ynew) >= fabs(y))  break;
        lambda = lambdaNew;
        y = ynew;
    }
    double DEL = lambda * lambda - lambda * Mz + Cxy;
    double cx = (Mzx * (Myy - lambda) - Mzy * Mxy) / DEL / 2;
    double cy = (Mzy * (Mxx - lambda) - Mzx * Mxy) / DEL / 2;
    centerx = cx + mx;   // recover origianl coordinate;
    centery = cy + my;
    radius = sqrt(cx * cx + cy * cy + Mz + 2 * lambda);
    return fitError(data, centerx, centery, radius);
}