Geometry & Recognition

이전 포스트에서 타원의 한 초점에서 나오는 빛은 타원에서 반사가 되면 다른 초점으로 모임을 보였다. 타원의 한 지점에서 접선벡터와 입사광선이 이루는 각이 접선벡터와 반사광의 이루는 각과 같음을 의미한다. 이를 물리적으로 증명하자. 그림과 같이 평면에 높인 타원이 있고, 타원의 한 지점에 길이 $L$인 줄을 걸친 후 두 초점을 통과하도록 만든다. 이제 초점을 통과한 줄 끝에 같은 무게의 추를 매단다. 두 추의 무게가 같으므로 추는 움직이지 않는 평형상태이고, 초점에서 매듭까지 거리를 각각 $d_1,d_2$라면 평면 아래로 늘어진 길이가 $L-d_1-d_2$이므로 두 추의 총 중력위치에너지는 $U=(d_1+d_2-L)mg$로 써진다.

그런데 타원상의 임의의 지점에서 $d_1+d_2=\text{const}$이므로 중력위치에너지는 일정하게 된다. 따라서 추의 중력 때문에 줄에 생기는 장력은 타원의 접선방향 성분은 없고 오직 타원의 접선에 수직한 성분만을 만든다(힘은 위치에너지의 그래디언트임). 줄이 매듭에 작용하는 장력(매듭에서 각 초점을 향하는 방향이다)을 각각 ${\overrightarrow{T}}_1$, ${\overrightarrow{T}}_2$라면, ${\overrightarrow{T}}_1+{\overrightarrow{T}}_2$는 매듭 위치에서 타원의 법선방향이어야 한다. 이는 ${\overrightarrow{T}}_1$과 ${\overrightarrow{T}}_2$의 접선성분이 같음을 의미하므로 두 줄이 접선과 이루는 각이 같게 되어 반사법칙을 보일 수 있게 된다.

초점의 위치가 각각 $\overrightarrow{a}$, $\overrightarrow{b}$인 타원이 있을 때 두 초점에서 타원상의 한 위치 $\overrightarrow{r}$까지 거리의 합은 항상 일정한 값을 가진다.

$$|\overrightarrow{r}-\overrightarrow{a}|+|\overrightarrow{r}-\overrightarrow{b}|=\text{const}$$

타원은 적당한 매개변수를 이용해서 매개화시킬 수 있으므로 이 매개변수에 대해서 윗식을 미분하자.

$$\dot{\overrightarrow{r}}\cdot \frac{\overrightarrow{r}-\overrightarrow{a}}{|\overrightarrow{r}-\overrightarrow{a}|}+\dot{\overrightarrow{r}}\cdot \frac{\overrightarrow{r}-\overrightarrow{b}}{|\overrightarrow{r}-\overrightarrow{b}|}=0$$

$\frac{\overrightarrow{r}-\overrightarrow{a}}{|\overrightarrow{r}-\overrightarrow{a}|}$는 초점 $\overrightarrow{a}$에서 $\overrightarrow{r}$을 향하는 단위벡터이고,  $\frac{\overrightarrow{r}-\overrightarrow{b}}{|\overrightarrow{r}-\overrightarrow{b}|}$는 초점 $\overrightarrow{b}$에서 $\overrightarrow{r}$을 향하는 단위벡터다. 이 두 단위벡터가 $\overrightarrow{r}$ 위치에서 접선벡터 $\dot{\overrightarrow{r}}$과 이루는 각을 각각 ${\theta }_a$, ${\theta }_b$라면

$$\cos {\theta }_a+\cos {\theta }_b=0$$임을 의미한다. $\cos {\theta }_a+\cos {\theta }_b=2\cos \frac{{\theta }_a+{\theta }_b}{2}\cos \frac{{\theta }_a-{\theta }_b}{2}$이므로 

$${\theta }_a+{\theta }_b=\pi$$

임을 알 수 있다.

두 초점에서 각각 타원의 한 지점까지 연결하는 선분이 그 점에서 접선과 이루는 각이 같으므로, ${\theta }_a=\pi -{\theta }_b$, 거울이 타원면 형태의 곡면으로 만들어졌다면 한 초점에서 발사된 빛은 타원에서 반사한 후 다시 다른 초점으로 모일 수 있음을 의미한다. 빛이 움직이는 경로의 길이가 모두 같으므로 한 초점에서 동시에 발사된 빛은 어느 방향을 향하는지에 상관없이 같은 시간에 다른 초점에 모이게 된다. 지붕과 벽면이 타원형으로 설계된 공연장에서 가수가 한 초점에서 노래를 부르면 반대편 초점에 앉은 관객에게는 마이크가 없더라도 다른 위치보다 상대적으로 더 잘 들리게 된다.

타원은 베지어 곡선을 이용해서 근사적으로 표현할 수 있음을 잘 알고 있다(https://kipl.tistory.com/313). 표준형 타원의 경우 1사분면에서 표현만 알면 대칭성에 의해서 나머지 사분면에서는 바로 구할 수 있으므로 거리를 구하는 점이 1사분면에 있는 경우만 알아보면 충분하다. 장축 반지름이 $a$, 단축 반지름이 $b$인 타원은 다음과 같이 베이지 곡선으로 표현된다.

$$\mathbf{B}(t)=(1-t^3){\mathbf{P}}_0+3t(1-t{)}^2{\mathbf{P}}_1+3t^2(1-t){\mathbf{P}}_2+t^3{\mathbf{P}}_3=\left(\begin{array}{c}ax(t)\\ by(t)\end{array}\right)$$

$$x(t)=3k(1-t{)}^{2}t+3(1-t)t^2+t^3$$

$$y(t)=x(1-t),0\le t\le 1$$

$k$ 값은 $t=1/2$일 때 $(a/\sqrt{2},b/\sqrt{2})$을 통과하는 조건을 부여하여 정하면

$$k=\frac{4}{3}(\sqrt{2}-1)=0.5522847498...$$임을 알 수 있다. 

한 점 $\mathbf{P}=(p,q)$에서 베지어 곡선 위의 한 점 $\mathbf{B}(t)$까지의 거리가 최단거리가 되기 위해서는 다음 직교 조건

$$f(t)=(\mathbf{B}(t)-\mathbf{P})\cdot {\mathbf{B}}^{\prime }(t)=0$$

$$f(t)=a^{2}x(t)x^{\prime }(t)-b^{2}y(t)y^{\prime }(t)-ap{x}^{\prime }(t)-bq{y}^{\prime }(t)=0$$

을 만족하는 근을 찾으면 된다. 풀어야 할 방정식이 5차이므로 직접적으로 근을 찾는 것은 불가능하므로 Newton-Raphson 방법을 사용하도록하자. 앞선 정확한 거리 계산에서와는 달리 초기 $t$ 값 설정에 민감하게 의존하지 않는다.

//cubic bezier_x: P0(0, 1), P1(k, 1),P2(1, k), P3(1, 0);
double B(double t) { //t in [0:1]
    const double k = 4.*(sqrt(2.)-1)/3;
    return t*(3*k + t*(3 - 6*k + (-2 + 3*k)*t));
}
// derivative of B(t);
double DB(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 3*k + t*(6 - 12*k + (-6 + 9*k)*t);
}
// derivative of DB(t);
double D2B(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 6 - 12*k + (-12 + 18*k)*t;
}
// ellipse radii=(a, b);
double dist2EllipseBezier3(double p, double q, double a, double b,
                           double& xt, double& yt) {
    if (a == b) return dist2Circle(p, q, a, xt, yt);
    double x = fabs(p), y = fabs(q);
    const double eps = 0.001 / max(a, b);
    double t = 0.5;  // mid
    while (1) {
        // Newton-Raphson;
        double f = a*a*B(t)*DB(t)-b*b*B(1-t)*DB(1-t)-a*x*DB(t)+b*y*DB(1-t);
        if (f == 0) break;
        double df = a*a*(DB(t)*DB(t)+B(t)*D2B(t))+b*b*(DB(1-t)*DB(1-t)+B(1-t)*D2B(1-t))
                    -a*x*D2B(t)-b*y*D2B(1-t);
        double dt = f / df;
        t -= dt;
        t = max(0, min(1, t));
        if (abs(dt) < eps ) break;
    }
 
    xt = a * B(t); yt = b * B(1-t);
    xt = p >= 0? xt: -xt;
    yt = q >= 0? yt: -yt;
    return hypot(p - xt, q - yt);
}