Geometry & Recognition

표면이 $+q$ 전하로 균일하게 대전된 반지름 $r$인 반구(hemispherical shell)가 있다.  반구의 대척점 P에서 전기장은 점전하 $+q$에서 $r$만큼 떨어진 지점에서 전기장의 몇 배일까?

1. $1-1/\sqrt{2}$

2. $1-1/\sqrt{2}$보다 작다.

3. $1-1/\sqrt{2}$보다 크다.

힌트: 복잡한 적분이 필요없다. solid angle 계산이 가능하면 된다.

매끄러운 벽에 그림처럼 막대를 $60^\circ$ 기울인 상태로 잡고 있다 가만히 놓는다. 바닥의 마찰도 없다. 손을 놓은 직후 미끄러지기 시작하는 막대의 (반시계방향) 회전 각가속도는?

1. $\frac{3g}{4L}$

2. $\frac{3g}{4L}$ 보다 크다.

3. $\frac{3g}{4L}$ 보다 작다.

4. 막대의 회전관성을 결정하는 질량/길이에 따라 다르다.

힌트: 질량중심+질량중심에 대한 회전운동 방정식을 사용해도 되고, 미끄러지는 동안 막대가 순간적인 회전을 한다고 보고 해결할 수도 있다(물론 회전축의 위치는 연속적으로 변한다).

막대는 언제 벽에서 떨어질까?(이전에 비슷한 질문이 있었다). 구체적으로 운동방정식을 풀어서 계산할 수 있는가?

최소 자승법 문제는 추정치와 실제의 측정치와의 차이를 최소로 만드는 parameter vector를 구하는 것이다.
$$\underset{\mathbf{x}} {\text{argmin}} ~|\mathbf{A}. \mathbf{x} - \mathbf{b}|^2.$$

여기서 design matrix $\mathbf{A}$는 추정에 사용이 된 basis 함수를 각각의 독립변수에서 계산한 값이고, $\mathbf{x}$는 구하고자 하는 parameter vector이며, $\mathbf{b}$는 측정값 vector이다. 예를 들면 주어진 측정값을 $n-1$차의 다항식을 이용하여서 피팅하려고 하는 경우에는 parameter는 다항식의 계수가 되며 $n$-차원의 vector space을 형성하게 되며, $\mathbf{A}$는

$$   A_{ij} = (x_i)^j ,  \quad j=0,..., n-1$$

가 일 것이다. 일반적으로 $A_{ij}$는 $x_i$에서 계산된 basis-함수의 값이 된다. 위의 식을 $\mathbf{x}$에 대해서 미분을 하면 극값을 취하면 최소자승 문제는 아래의 행렬식을 푸는 문제가 된다

$$    (\mathbf{A}^T. \mathbf{A}) .\mathbf{x} =  \mathbf{A}^T.  \mathbf{b}.$$

$\mathbf{A}^T. \mathbf{A}$ 은 $n\times n$ matrix다. 이 행렬이 역행렬을 가지게 되면

 $$ \mathbf{x} = (\mathbf{A}^T. \mathbf{A})^{-1} . (\mathbf{A}. \mathbf{b}),$$

를 하여서 원하는 parameter vector를 얻을 수 있다. 그러나 피팅 문제에 따라 행렬 $\mathbf{A}^T. \mathbf{A}$가 매우 singular 해져 역행렬을 구할 수 없게 되는 경우에 종종 생긴다. 예를 들면, 저주파의 신호를 고주파 기저 함수를 이용하여서 최소자승법을 사용하고자 하는 경우 등에 이러한 문제에 부딪히게 된다. 이런 경우에는 직접적으로 $\mathbf{A}^T. \mathbf{A}$의 역행렬을 구하는 방법을 이용할 수 없고

$$   \mathbf{A} .\mathbf{x} =  \mathbf{b}$$

의 식을 $\mathbf{A}$의 SVD(Singular Value Decomposition)를 이용하여서 풀 수가 있다. $\mathbf{A}$를 SVD 하면 $\mathbf{A}_{m\times n}=\mathbf{U}_{m\times n} . \mathbf{w}_{n\times n}. \mathbf{V}_{n\times n}^T $의 형태로 분해할 수 있다. 여기서 $\mathbf{w}=\text{diag}(\underbrace{w_0, w_1,...}_{\text{nonzero}},0,..,0)$로 쓰여지는 대각행렬이다. matrix $\mathbf{U}$와 $\mathbf{V}$의 column vector를 사용하면
$$ \mathbf{A}  =\sum_{w_k \ne 0} w_k \mathbf{u}_k \otimes \mathbf{v}_k^T$$

의 형태로 쓰인다. $\mathbf{u}_k$는 $\mathbf{U}$의 $k$-번째 열벡터이고, $\mathbf{v}_k$는 $\mathbf{V}$의 $k$-번째 열벡터로 각각 orthonormal basis를 형성한다. parameter 벡터를 $\{ \mathbf{v}_k \}$ basis로 전개를 하면 영이 아닌 singularvalue에 해당하는 성분만 가지게 된다. 구체적으로 위의 $\mathbf{A}$ 분해와 $\mathbf{u}_j^T.\mathbf{u}_k=\delta_{jk}$, 그리고 $\sum_k \mathbf{v}_k \otimes \mathbf{v}_k^T= \mathbf{I}_{n\times n}$임을 이용하면,

\begin{gather}  \mathbf{v}_k^T . \mathbf{x} = \mathbf{u}_k^T . \mathbf{b} / w_k, \quad w_k \ne 0, \\                    \mathbf{v}_k^T . \mathbf{x} = 0, \quad w_k = 0, \\  \rightarrow ~~\mathbf{x} = \sum _{w_k \ne 0 } ( \mathbf{u}_k^T . \mathbf{b} / w_k)  \mathbf{ v} _k , \end{gather}

이어서 위의 해를 구할 수 있다. 이 해는 $|\mathbf{A} . \mathbf{x} -  \mathbf{b}|^2$를 최소화한다.

int svd(double *A, int m, int n, double* w, double *V); // from cmath libary.
void fit_func(double x, double val[], int n) {          // polynomial fitting sample;
    val[0] = 1;
    for(int i = 1; i < n; ++i)
        val[i] = x * val[i - 1];
}
#define EPSILON 1.E-8
int svd_fit(const double x[], const double y[], const int m, const int n,
            void (*fit_func)(double , double [], int ),
            double params[],
            double *error)
{
    double *A = new double [m * n];
    double *w = new double [n];
    double *V = new double [n * n];
    // evaluate design matrix;
    for (int i = 0; i < m; ++i)
        fit_func(x[i], &A[i * n + 0], n) ;

    svd(A, m, n, w, V);
    // now A becomes U matrix;
    // truncate small singular values;
    double wmax = 0;
    for (int i = 0; i < n; ++i)
        if (w[i] > wmax) wmax = w[i];
    double thresh = wmax * EPSILON;
    for (int i = 0; i < n; ++i)
        if (w[i] < thresh) w[i] = 0;
    
    // back substitution;
    double *tmp = new double [n];
    for (int j = 0; j < n; ++j) {
        double s = 0;
        if (w[j]) {
            for (int i = 0; i < m; ++i)
                s += A[i * n + j] * y[i];
            s /= w[j];
        }
        tmp[j] = s;
    }
    for (int j = 0; j < n; ++j) {
        double s = 0;
        for (int jj = 0; jj < n; ++jj)
            s += V[j * n + jj] * tmp[jj];
        params[j] = s;
    };

    //estimate error;
    *error = 0;
    for (int i = 0; i < m; ++i) {
        fit_func(x[i], &A[i * n + 0], n); //use A as a tmp buffer;
        double sum = 0;
        for (int j = 0; j < n; ++j) sum += params[j] * A[i * n + j] ;
        double err = (y[i] - sum);
        *error += err * err ;
    }
    delete[] A; delete[] w; delete[] V;
    delete[] tmp;
    return 1;
}

RGB 컬러 이미지의 gray level의 cdf를 이용해서 histogram equalization을 수행한다. 컬러 이미지의 gray level은 

$$gray = \frac{r + g+ b}{3}$$

으로 정의한다.

void test_color_equalize_new(CRaster& raster, CRaster& out) {
    int hist[256] = {0};
    int chist[256], lut[256], partition[256 + 1];
    CSize sz = raster.GetSize();
    out = raster; // clone; 
    // pdf of gray level: g = (r + g + b) / 3
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++){
            int a = *p++; a += *p++; a += *p++;
            hist[a / 3]++;
        }
    };
    // cdf;
    chist[0] = hist[0];
    for (int i = 1; i < 256; i++)
        chist[i] = hist[i] + chist[i - 1];
    
    /* assign equal number of pixels in each gray levels;*/
    int num_pixels = sz.cx * sz.cy;
    double pixels_per_level = double(num_pixels) / 256;

    /* first and last in partition */
    partition[0]   = 0;
    partition[256] = 256;
    /* intermediate; */
    for (int j = 0, i = 1; i < 256; i++) {
        double desired = i * pixels_per_level;			
        while (chist[j + 1] <= desired) j++;
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j + 1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* fill equalization LUT */
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i + 1] <= j) i++;
        lut[j] = i;
    } 
    // needs hue preserving processes;
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)out.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            *q++ = lut[*p++]; *q++ = lut[*p++]; *q++ = lut[*p++];
        }
    }
}

void test_color_equalize_HSV(CRaster& raster, CRaster& out) {
    int hist[256] = {0};
    int chist[256] = {0}, lut[256] = {0}, partition[256 + 1];
    CSize sz = raster.GetSize();
    out = raster; // cloning;
    std::vector<double> fH(sz.cx * sz.cy);
    std::vector<double> fV(sz.cx * sz.cy);
    std::vector<double> fS(sz.cx * sz.cy);

    int n = sz.cx * sz.cy;
    double r, g, b;
    for (int y = 0, k = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++, k++){
            b = *p++, g = *p++, r = *p++;
            RGBtoHSV(r / 255., g / 255., b / 255., fH[k], fS[k], fV[k]);
        }
    };	
    // make histogram of V-color;
    for (int k = 0; k < n; k++)
        hist[int(fV[k] * 255)]++;

    // cdf;
    chist[0] = hist[0];
    for (int i = 1; i < 256; i++) {
        chist[i] = hist[i] + chist[i - 1];
    }
    /* assign equal number of pixels in each V-color;*/
    int num_pixels = sz.cx * sz.cy;
    double pixels_per_level = double(num_pixels) / 256;

    /* first and last in partition */
    partition[0]   = 0;
    partition[256] = 256;
    /* intermediate; */
    for (int j = 0, i = 1; i < 256; i++) {
        double desired = i * pixels_per_level;			
        while (chist[j + 1] <= desired) j++;
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j + 1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* fill equalization LUT */
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i + 1] <= j) i++;
        lut[j] = i;
    } 
    for (int k = 0; k < n; k++)
        fV[k]= lut[int(fV[k] * 255)] / 255.;

    for (int y = 0, k = 0; y < sz.cy; y++) {
        BYTE *q = (BYTE *)out.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++, k++) {
            HSVtoRGB(fH[k], fS[k], fV[k], r, g, b);
            *q++ = int(b * 255); 
            *q++ = int(g * 255);
            *q++ = int(r * 255);
        }
    }
}

/* fR Red component, range: [0, 1]
** fG Green component, range: [0, 1]
** fB Blue component, range: [0, 1]
** fH Hue component, range: [0, 360]
** fS Hue component, range: [0, 1]
** fV Hue component, range: [0, 1] */
void RGBtoHSV(double fR, double fG, double fB, double& fH, double& fS, double& fV);
void HSVtoRGB(double fH, double fS, double fV, double& fR, double& fG, double& fB);