Geometry & Recognition

질량 $M_*$이고, 반지름이 $R_*$인 행성의 표면에 그림과 같이 밑변이 서로 연결된 $N$개의 물기둥(communicating vessels)을 고려하자. 각 물기둥의 단면적은 $A_i$이고, 처음 물기둥의 (별의 중심에서 잰) 높이는 $R_i  \ge R_*$로 되어 있다. 이웃하는 물기둥을 연결하는 밸브를 열기 전 물기둥에 담긴 물의 중력위치에너지는 (물의 밀도 $\rho$)

$$ U_\text{initial} = -\sum \int_{R_*}^{R_i} \frac{GM_*  \rho A_i dr}{r} = -GM_* \rho \sum A_i \ln \left( \frac{R_i}{R_*}  \right) $$

이제 관을 연결하는 밸브를 모두 열면 수압의 차이에 의해서 각 물기둥 사이에 물의 교환이 생기게 되고, 결국에는 모든 물기둥의 높이는 공통의 값 $R_i \to \overline{R}$을 가지게 되어 더 이상 물의 교환이 생기지 않는 평형상태에 도달한다. 이 과정에서 물의 총량(총 부피)은 보존이 되므로 

$$ \sum A_i (R_i - R_*) = \sum A_i ( \overline{R} - R_*) $$

를 만족하므로 평형상태에서 물기둥의 높이는
$$ \overline{R}  = \frac{\sum A_i R_i }{\sum A_i} = \sum p_i  R_i,\qquad p_i \equiv \frac{A_i}{\sum A_k }$$

물기둥이 평형에 이르기 위해서는 일부 에너지를 마찰등에 의해서 일어버리는 과정이 있어야 한다. 그렇지 않으면 각각의 물기둥은 영원히 출렁거리게 된다. 평형이 이르는 물리적인 프로세스가 다음의 부등식이 항상 성립함을 보증함을 알 수 있다.

$$ U_\text{initial} \ge U_\text{final}\\ -GM_*\rho\sum A_i  \ln \left( \frac{R_i }{R_*} \right) \ge  -GM_* \rho \sum A_i \ln \left( \frac{\bar{R}}{R_*} \right) $$

등호는 처음 물기둥의 높이가 모두 같은 경우에 성립한다.  이 식을 다시 정리하면 각각의 물기둥의 단면적으로 가중치를 준 물기둥 높이의 산술평균(arithmatic mean)은 기하평균(geometric mean)보다도 크거나 같다는 것을 보여준다.

$$ \overline{R}  \ge     \prod R_i^{p_i}$$

$R_i$가 일반적인 양의 실수이므로 위 부등식은 일반적인 양의 실수 사이의 가중치 산술평균과 기하평균 사이에 다음의 부등식이 성립함을 의미한다.

$$ \sum p_i R_i \ge \prod R_i^{p_i},\qquad \left( \sum p_i =1, ~p_i \ge0 \right)$$

이 부등식은 수학적인 증명 대신 물리계인 communicating vessels 시스템이 평형상태에 도달하는 과정에서 물의 양은 보존되지만, 평형상태에 이르는 과정에서 계가 항상 더 낮은 에너지 상태로 옳겨간다는 사실만을 이용하여 보였다.

이제 물기둥 시스템의 크기가 행성의 크기보다 매우 작을 때를 고려하자. 즉, 행성 표면에서 잰 물기둥의 높이가 $h_i = R_i - R_*  \ll R_*$인 경우 앞의 결과를 $h_i^2$의 order까지 전개하면 다음과 같은 식을 얻을 수 있다.

$$-\frac{GM_*\rho}{R_*^2 }  \sum A_i \left(  R_*\bar{h}  -\frac{1}{2}\bar{h}^2  \right) \ge 
-\frac{GM_* \rho}{R_*^2 }\sum A_i   \left(R_* {h_i} -\frac{1}{2} {h_i^2}\right)$$

$\sum A_i \overline{h}= \sum A_i h_i$이므로 

$$ \frac{ \sum A_i h_i^2}{ \sum A_j} \ge \bar{h}^2$$

따라서, $x_i^2= A_i h_i^2$, $y_i^2 = A_i$로 놓으면, $\bar{h}^2 = (\sum A_i  h_i  )^2/ (\sum A_k )^2 = (\sum x_i y_i )^2/(\sum y_i^2 )^2$이므로 위의 결과에서 Cauchy-Schwartz 부등식을 얻을 수 있음을 알려준다. $$\left( \sum x_i^2  \right)\left(\sum y_i^2\right) \ge \left(\sum x_i y_i\right)^2 $$

Ref: https://grail.cs.washington.edu/projects/digital-matting/image-matting/  

// knockout method;
// trimap을 만들 때 전경/배경/중간 영역(경계) 구분;
#define BACKGROUND 0    // 배경 마스크 색상;
#define FOREGROUND 255  // 전경 마스크 색상;
#define UNKNOWN    128  // 경계 영역의 생상;
#define EDGEMARK   3    // 전경과 배경의 경계선을 표시하는 색상
// 주어진 픽셀에서 경계점까지 떨어진 거리(제곱)에 따라 경계점을 
// 정렬하기 위해 필요함;
struct PtDist {
    int x, y, d;
    PtDist() {}
    PtDist(int x, int y): x(x), y(y), d(0) {}
};
struct RGBTRIO {
    double rgb[3];  // b-g-r
};
static 
bool compare(const PtDist& a, const PtDist& b) {
    // 거리 제곱의 ascending order로 정렬;
    return a.d < b.d;
}
// 최소거리의 DISTANCE_THRESHOLD배 이내에 있는 경계 픽셀만을 이용해서 
// alpha값을 추정한다:
#define DISTANCE_THRESHOLD  (5.0)
static
RGBTRIO estimateColor(int x0, int y0, CRaster& raster, 
                   std::vector<PtDist>& pts) {
    for (int i = pts.size(); i-->0;) {
        int x = pts[i].x - x0; 
        int y = pts[i].y - y0;
        pts[i].d = x * x + y * y; 
    }
    std::sort(pts.begin(), pts.end(), compare);
    const double minDist = sqrt(double(pts[0].d));
    ASSERT(minDist > 0);
    const double maxDist = DISTANCE_THRESHOLD * minDist;
    const double inv_diff = 1.0/ (maxDist - minDist);
    double wsum = 0, rsum = 0, gsum = 0, bsum = 0, d;
    for (int i = 0; i < pts.size(); i++) { // sorted data;
        if ((d = sqrt(double(pts[i].d))) < maxDist) {
            double w = (maxDist - d) * inv_diff;
            DWORD q = raster.GetPixel0(pts[i].x, pts[i].y);
            bsum += w * ((q      ) & 0xFF); // blue;
            gsum += w * ((q >>  8) & 0xFF); // green;
            rsum += w * ((q >> 16) & 0xFF); // red;
            wsum += w;
        }
        else break;
    }
    if (wsum == 0) wsum = 1;
    RGBTRIO Q;
    Q.rgb[0] = (bsum / wsum);
    Q.rgb[1] = (gsum / wsum);
    Q.rgb[2] = (rsum / wsum);
    return Q;
}
#define AddVec(A,B,C)	{C[0]=A[0]+B[0];C[1]=A[1]+B[1];C[2]=A[2]+B[2];}
#define SubVec(A,B,C)	{C[0]=A[0]-B[0];C[1]=A[1]-B[1];C[2]=A[2]-B[2];}
#define NormSq(A)       (A[0]*A[0]+A[1]*A[1]+A[2]*A[2])
#define DotProd(A,B)	(A[0]*B[0]+A[1]*B[1]+A[2]*B[2])
#define ScaleVec(A,s)   {A[0]*=s;A[1]*=s;A[2]*=s;}
static
double estimateAlpha(double F[3], double B[3], double C[3]) {
    double N[3], V[3], Bp[3];
    SubVec(F, B, N);						//N=F-B;
    double nn = NormSq(N);					//nn=N.N;
    if (nn == 0) return 1;	
    SubVec(C, B, V);						//V=C-B;
    double proj = DotProd (N, V) / nn;		//proj=(V.N)/(N.N);
    ScaleVec(N, proj);						//N=N(V.N)/(N.N);
    SubVec(C, N, Bp);						//Bp=C-N(V.N)/(N.N);
    //
    SubVec(F, Bp, N);						//N=F-Bp;
    nn = NormSq(N);							//nn=N.N;
    if (nn == 0) return 1;
    SubVec(C, Bp, V);						//V=C-Bp;
    proj = DotProd(N, V) / nn;				//proj=(N.V)/(N.N);

    return proj < 0 ? 0: proj > 1 ? 1: proj;
}
void FindMatt(CRaster& raster, CRaster& trimap, CRaster& alphaMap) {
    ASSERT((trimap.GetBPP()==8) && (raster.GetBPP()==24));
    const CSize sz = raster.GetSize();
    alphaMap.SetDimensions(sz, 32);// Bitmap with alpha-channel;
    std::vector<PtDist> ptFore = binaryTracing(trimap, FOREGROUND);
    std::vector<PtDist> ptBack = binaryTracing(trimap, BACKGROUND);
    
    double C[3];
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)trimap.GetLinePtr(y);
        BYTE *q = (BYTE *)raster.GetLinePtr(y);
        BYTE *r = (BYTE *)alphaMap.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            if (UNKNOWN == *p) {// unknown pixel(0=background,255=foreground);
                RGBTRIO Qf = estimateColor(x, y, raster, ptFore);
                RGBTRIO Qb = estimateColor(x, y, raster, ptBack);
                C[0] = q[0]; C[1] = q[1]; C[2] = q[2];
                double alpha = estimateAlpha(Qf.rgb, Qb.rgb, C);
                // if ( alpha >= 0.5)
                r[0] = (int)max(0, min(255, Qf.rgb[0]));
                r[1] = (int)max(0, min(255, Qf.rgb[1]));
                r[2] = (int)max(0, min(255, Qf.rgb[2]));
                r[3] = int(alpha * 255);       //[0:1]->[0:255]
            } else if (FOREGROUND == *p) {
                r[0] = q[0];
                r[1] = q[1];
                r[2] = q[2];
                r[3] = 0xFF;
            }
            p++; q += 3; r += 4;
        }
    }
}

$$\text{Cauchy-Schwartz Inequality:}~~|\mathbf{a}|^2 |\mathbf{b}|^2  \ge \left( \mathbf{a}\cdot \mathbf{b} \right)^2$$

이 부등식을 간단한 물리법칙을 이용해서 증명하자. 우선 그림과 같이 $N$개의 축전기가 연결된 회로가 있다. 처음 축전기에는 각각 일정한 전하가 충전되어 있고, 회로를 연결하는 스위치는 모두 열린 상태이다.

각 축전기의 전기용량을 $C_i$, 충전된 전하를 $Q_i$라면 축전기에 양극판에 걸리는 전위차는 $V_i = Q_i /C_i$로 주어진다. 전하의 상대적인 부호에 따라 전위차도 +/- 부호를 가질 수 있다. 충전과정에서 축전기에는 에너지가 제공되므로 충전된 축전기는 에너지를 가진다. 각 축전기에 저장된 에너지가 $E_i = \frac{1}{2} Q_i V_i \ge0$이므로 축전기 회로에 저장된 총에너지는

$$ E_\text{initial} = \frac{1}{2} \sum Q_i V_i$$

이제 회로의 모든 스위치를 닫으면 각 축전기의 전위차 때문에 전하의 재분배가 발생하는 데 ($Q_i \to Q'_i$), 이 과정은 각 축전기의 전위차가 모두 같을 때까지( $V_i \to V_\text{eq}$) 일어난다. 물론 재분배 과정에서 회로의 저항때문에 에너지의 일부를 잃어버리지만 전하보존에 의해서 총전하량은 변함이 없어야 한다. 스위치를 닫았을 때 모든 축전기는 병렬연결이므로 등가전기용량은 $C_\text{eq} = \sum C_i = \sum Q_i/V_i $이므로 평형상태에서의 공통 전위차는 

$$ C_\text{eq} V_\text{eq} = Q_\text{total} = \sum Q_i$$

을 만족한다. 그리고 평형상태에 도달했을 때 축전기에 저장된 총 에너지를 스위치를 닫기 전 전하와 전위차로 표현하면

$$ E_\text{final} = \frac{1}{2} \sum Q'_i    V_\text{eq} = \frac{1}{2}\left(\sum Q_i \right) V_\text{eq}  = \frac{1}{2} \frac{\left( \sum Q_i \right)^2}{  C_\text{eq}} =\frac{1}{2} \frac{\left(\sum   Q_i \right)^2}{\sum Q_i/ V_i } $$

전하의 재분배 과정에서 에너지의 일부를 잃어버리므로 평형상태에서의 에너지는 항상 처음보다 클 수는 없다.

$$ E_\text{initial} \ge E_\text{final}$$

따라서

$$  \left( \sum Q_i V_i \right) \left( \sum \frac{Q_i}{V_i} \right) \ge \left( \sum Q_i \right)^2$$

축전기의 전기용량이 물리적으로 항상 양수이므로 전하와 전위차는 같은 부호를 가져야 한다. 따라서  $a_i^2 = Q_i V_i$, $b_i^2 = Q_i /V_i$로 놓으면 다음과 같은 Cauchy-Schwartz 부등식을 얻을 수 있다.

$$ \left(\sum a_i^2 \right) \left(\sum b_i^2 \right) \ge \left( \sum   a_i b_i \right)^2$$

여기서 등호는 축전기의 처음 전위차가 모두 같은 경우로 이때는 스위치를 닫아도 전류의 흐름이 생기지 않아 에너지의 손실이 발생할 수 없는 경우에 해당한다. 이 부등식은 전하보존법칙과 물리계가 평형에 이르는 과정에서 에너지적으로 가장 낮은 상태로 옳겨간다는 사실을 이용하여 증명하였다.

In computational geometry, an alpha shape, or α-shape, is a family of piecewise linear simple curves in the Euclidean plane associated with the shape of a finite set of points. https://en.wikipedia.org/wiki/Alpha_shape

class AlphaShape {
    double alpha;
    std::set<std::pair<int,int>> Border; //pair;
public:
    AlphaShape() { }
    AlphaShape(CDC *pDC, std::vector<CfPt>& Q, const double radius = 50) {
        if (Q.size() < 2) return;
        alpha = radius;
        const double alphaSq = alpha * alpha;
        const double diameter = 2 * alpha;
        for (int i = 0; i < Q.size(); ++i) {
            const CfPt &p = Q[i];
            for (int j = 0; j < i; ++j) {
                const CfPt &q = Q[j];
                if (q == p) continue;
                const double edist = p.Dist(q); 
                // 두 점의 사잇거리가 지름보다 크면 edge가 안됨;
                if (edist > diameter) continue;                     

                // 두 점을 원주에 포함하는 반지름 alpha인 두 원을 찾음;
                // note that c1 == c2 if edist == diameter;
                double scale = sqrt(alphaSq - (edist/2)*(edist/2));
                CfPt n = (q - p) * scale / edist;
                CfPt mid = (q + p) / 2;
                CfPt c1 = mid + n.Rot90(); // CfPt(mid.x - ny, mid.y + nx);
                CfPt c2 = mid - n.Rot90(); // CfPt(mid.x + ny, mid.y - nx);

                // alpha shape의 변이 되기 위해서는 두 원 중 하나는
                // 다른 점들을 포함하지 않아야 함;
                bool c1_empty = true, c2_empty = true;
                for (int k = Q.size(); (k-->0) && (c1_empty || c2_empty);) {
                    if (Q[k] == p|| Q[k]== q) continue; //i==k || j==k?
                    if (c1.DistSq(Q[k]) < alphaSq) c1_empty = false;
                    if (c2.DistSq(Q[k]) < alphaSq) c2_empty = false;            
                }

                if (c1_empty || c2_empty) {                  
                    if (c1_empty) drawCircle(pDC, c1);
                    if (c2_empty) drawCircle(pDC, c2);
                    // draw edge;
                    drawEdge(pDC, p, q);
                    //Border.insert(std::make_pair(i,j));
                }
            }
        }
    } 
    void drawEdge(CDC *pDC, const CfPt &p, const CfPt &q) {
        CPen pen(PS_SOLID,2,RGB(255,0,255)); //magenta;
        CPen *pOld = pDC->SelectObject(&pen);
        pDC->MoveTo(p); pDC->LineTo(q);
        pDC->SelectObject(pOld);
    }
    void drawCircle(CDC *pDC, const CfPt &c) {
        CPen pen(PS_SOLID, 1, RGB(0,200,200));
        CPen *pOld = pDC->SelectObject(&pen);
        pDC->SelectObject(GetStockObject(NULL_BRUSH));
        pDC->Ellipse(c.x-alpha, c.y-alpha, c.x+alpha, c.y+alpha);
        pDC->SelectObject(pOld);		
    }
    std::set<std::pair<int,int>> getBorder() {
        return Border;
    }
};

https://en.wikipedia.org/wiki/Viviani%27s_theorem

정삼각형 내부 한 지점에서 세 변에 내린 수선의 길이 합은 항상 일정하다(역도 성립한다).  기본적인 물리법칙을 이용하면 이 사실을 쉽게 보일 수 있다. 같은 무게의 추가 달린 세 줄을 묶아 매듭을 만든 후 정삼각형의 각 변에 줄이 걸치도록 설치한다. 각 줄에 걸리는 추의 중력에 의한 장력이 같으므로 매듭을 기준으로 세 줄은 서로 120도의 각을 이룰 때 평형이 되어 더 이상 움직이지 않게 된다(Lami의 법칙). 120도 조건만 만족하면 힘의 평형상태이므로 매듭이 정삼각형의 어느 지점에 있더라도 움직이지 않게 된다. 그런데 평형상태는 중력 위치에너지가 가장 낮아지는 상태에 해당하고, 그 값은 삼각형의 아래로 늘어진 줄 길이의 합에 비례한다. 매듭의 위치가 바뀌면 각 줄의 늘어진 길이는 달라질 수 있지만 중력위치에너지는 변하지 않아야 하므로 그 합은 같아야 한다(추가: 수선 조건을 유지하면 매듭의 위치를 미소변위만큼 이동시킬 때 위치에너지 차이가 있다면 매듭에 힘이 작용한다는 의미이므로 평형조건에 위배된다. 힘은 위치에너지함수의 -gradient이기 때문이다). 그런데 각 줄의 길이가 고정되어 있으므로 삼각형 위쪽에 있는 줄의 길이의 합도 고정되어야  함을 알 수 있다. 이 줄의 합은 실제로 한 꼭짓점에서 내린 수선의 길이와 같음은 매듭이 한 꼭짓점에 접근할 때를 생각해 보면 쉽게 알 수 있다. 힘의 평형조건을 고려하면 역도 성립함을 쉽게 알 수 있다. 그리고 일반적인 등각다각형에 대해서도 성립한다.