\begin{gather}BF[I]_{\bf p} = \frac{1}{W_\bf{p}}  \sum_{ {\bf q} \in S} G_{\sigma_s} (||{\bf p} - {\bf q}||) G_{\sigma_r}(|I_{\bf p} - I_{\bf q} |) I_{\bf q} \\  W_{\bf p} = \sum_{{\bf q}\in S} G_{\sigma_s} (||{\bf p}-{\bf q}||) G_{\sigma_r}(|I_{\bf p} - I_{\bf q} |) \\ G_\sigma ({\bf r}) = e^{ - ||\bf{r}||^2/2\sigma^2 }\end{gather}

 

Bilateral Filter
Gaussian Filter

 

smoothing based on the nonlinear heat eq

// sigmar controls the intensity range that is smoothed out. 
// Higher values will lead to larger regions being smoothed out. 
// The sigmar value should be selected with the dynamic range of the image pixel values in mind.
// sigmas controls smoothing factor. Higher values will lead to more smoothing.
// convolution through using lookup tables.
int BilateralFilter(BYTE *image, int width, int height, 
    double sigmas, double sigmar, int ksize, BYTE* out) {
    //const double sigmas = 1.7;
    //const double sigmar = 50.;
    double sigmas_sq = sigmas * sigmas;
    double sigmar_sq = sigmar * sigmar;
    //const int ksize = 7;
    const int hksz = ksize / 2;
    ksize = hksz * 2 + 1;
    std::vector<double> smooth(width * height, 0);
    // LUT for spatial gaussian;
    std::vector<double> spaceKer(ksize * ksize, 0);
    for (int j = -hksz, pos = 0; j <= hksz; j++) 
        for (int i = -hksz; i <= hksz; i++) 
            spaceKer[pos++] = exp(- 0.5 * double(i * i + j * j)/ sigmas_sq); 
    // LUT for image similarity gaussian;
    double pixelKer[256];
    for (int i = 0; i < 256; i++)
        pixelKer[i] = exp(- 0.5 * double(i * i) / sigmar_sq);

    for (int y = 0, imgpos = 0; y < height; y++) {
        int top = y - hksz;
        int bot = y + hksz;
        for (int x = 0; x < width; x++) {
            int left = x - hksz;
            int right = x + hksz;
            // convolution;
            double wsum = 0;
            double fsum = 0; 	
            int refVal = image[imgpos];
            for (int yy = top, kpos = 0; yy <= bot; yy++) {
                for (int xx = left; xx <= right; xx++) {
                    // check whether the kernel rect is inside the image;
                    if ((yy >= 0) && (yy < height) && (xx >= 0) && (xx < width)) {
                        int curVal = image[yy * width + xx];
                        int idiff = curVal - refVal;
                        double weight = spaceKer[kpos] * pixelKer[abs(idiff)];
                        wsum += weight;
                        fsum += weight * curVal;
                    }
                    kpos++;
                }
            }
            smooth[imgpos++] = fsum / wsum;
        }
    }

    for (int k = smooth.size(); k-- > 0;) {
        int a = int(smooth[k]);
        out[k] = a < 0 ? 0: a > 255 ? 255: a;
    }
    return 1;
}
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

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

Cubic Spline Kernel  (1) 2024.03.12
Ellipse Fitting  (0) 2024.03.02
영상에서 Impulse Noise 넣기  (2) 2023.02.09
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Posted by helloktk
,

열방정식은 매질에서 열이 전달되는 과정을 기술한다. 열은 온도차이가 있을 때 전달되는 에너지로 온도 분포 $u(\vec{r},t)$가 주어진 경우 주변으로 나가는 열에너지는 온도의 gradient $\nabla u$에 비례한다. 

$$ \text{heat flux: }\quad  \vec{j}= \kappa (u)\nabla u$$

$\kappa$는 매질에서 열이 얼마나 잘 전도는지를 표현하는 일반적으로 $u$에 의존할 수 있다. 균일하고 등방적인 매질인 경우는 $\kappa$가 일정한 상수에 해당한다. 국소적으로 단위부피당 단위시간당 빠져나가는 열에너지는 thermal flux의 divergence에 비례하고 이 값이 0이 아니면 그 지점의 온도 변화를 만든다. 따라서 온도는 다음과 같은 방정식을 만족한다.

$$ u_t = \nabla \cdot \vec{j}  = \nabla\cdot(\kappa \nabla u)$$

이 방정식은 물질의 확산 현상에도 적용이 가능한데 이 경우 $u$는 물질의 밀도분포에 해당할 것이다.

$\kappa$가 상수인 경우 이 방정식은 선형방정식으로 초기 온도분포 $u(x,t=0)=f(x)$가 주어지면 해는 사이즈가 $\sigma = \sqrt{2\kappa t}$인 가우시안 커널(heat kernel 또는 Greeen function)과 초기 온도분포의 convolution으로 표현할 수 있다. 1차원 선형 열방정식의 해는 구체적으로 다음과 같이 표현된다.

\begin{align} u(x,t) =  \int \frac{1}{\sqrt{4\pi \kappa t}} e^{- (x-y)^2 / 4\kappa t } f(y) dy = ( G_{\sigma} * f)(x,t)\\G_\sigma (x) = \frac{1}{\sqrt{2\pi\sigma} }  e^{- x^2/2\sigma^2} \end{align}

열방정식은 뜨거운 커피가 담긴 잔을 방에 놓았을 때 방 안의 온도가 어떻게 변할지도 알 수 있게 해 준다. 물론 경험적으로 커피의 온도는 내려가고 방안은 잔의 주변부터 점차로 데워져서 나중에서는 모든 부분이 일정한 온도를 가지는 평형상태에 도달한다. 이는  가우시안 커널의 크기가 $\sqrt{t}$에 비례하므로 시간이 흐를수록 모든 지점에서 온도는 일정한 값으로 수렴할 것임을 쉽게 예측할 수 있다.

이 열방정식을 이미지에 대해서도 적용할 수 있다. 이 경우 초기 온도분포는 주어진 이미지의 명암값에 해당한다. 열방정식의 우변을 이미지에 적용하면 Gaussian 블러링에 해당하고, 시간은 블러링의 순차적인 적용 횟수를 의미한다. 이미지가 블러링 되는 경우 영상이 담고 있는 정보는 점점 잃게 되는데 이를 방지하기 위해서는 영상의 에지 정보를 보존하는 방법으로 열방정식을 변형시켜 보자. 영상에서 에지영역에서는 gradient 값이 커지므로 이 지점에서는 커널사이즈를 작게 (열방정식 관점에서는 열전도도를 줄임) 만들면 블러링 효과가 줄어들게 된다. 이를 위해서 다음과 같은 열전도도 함수를 도입하자.

$$ \kappa (u) = \frac {1}{1+ |\nabla u|^2/\lambda^2 }$$

여기서 $\lambda$은 contrast 파라미터로 에지의 세기가 $|\nabla| \gg \lambda$ 영역에서는 블러링의 영향을 거의 받지 않게 되어 열방정식을 적용하더라도 에지에 대한 정보 손실은 거의 일어나지 않는다. 그러나 $\kappa $가 상수가 아닌 경우는 열방정식은 비선형이므로 해를 명시적으로 쓸 수 없고 오직 수치해석적으로만 구할 수 있다.

먼저 미분연산자(Sobel operator: $\lambda$의 설정은 미분연산자의 normalization을 고려해야 함)을 이용해서 이미지의 gradient을 구하면 flux $\vec{j}$을 구할 수 있고, 다시 $j_x$에 대해서 $x-$방향 미분, $j_y$에 대해서 $y-$방향 미분을 해서 $\vec{j}$의 divergence, 즉 이미지의 일반화된 Laplace 연산 결과를 얻을 수 있다. 이 결과에 시간간격을 곱한 양을 이전 이미지에 더하여 새로운 이미지를 얻는다. 이 과정을 반복적으로 수행하면 된다.

$$ u(x,y, t+\Delta t) = u(x, y, t) + \Delta t \times \nabla \cdot (\kappa (|\nabla u| ) \vec{j})$$

$\kappa$가 상수일 때 $t \to \kappa t$로 시간변수를 바꾸면 열방정식은 차원이 없는 형태로 되고, $\kappa <1$이므로 시간간격은 $\Delta t \ll 1$를 만족하도록 선택해야 한다.

 
 
 

 

$\lambda=50,~\Delta t=0.01$일 때 5 스텝 간격으로의 변화

 
double conductivity(double gmagsq, double lambdasq) {
    return 1.0/ (1.0 + gmagsq/lambdasq);
}
int LaplaceImage(double *img, int w, int h, double lambda, double *laplace) {
    const double lambdasq = lambda * lambda;
    // Sobel gradient 3x3
    const int nn[] = { -w - 1, -w, -w + 1, -1, 0, 1, w - 1, w, w + 1};
    const int sobelX[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1};
    std::vector<double> Gx(w * h, 0);
    std::vector<double> Gy(w * h, 0);
    std::vector<double> Gmag2(w * h, 0);
    const int xmax = w - 1, ymax = h - 1;
    // gradient-x, gradient-y, and mag of gradient.
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                double v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            Gx[pos] = sx;
            Gy[pos] = sy;
            // gx^2 + gy^2;
            Gmag2[pos] = sx * sx + sy * sy;
        }
        pos++; // skip x=xmax;
    }
    // flux;
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double kappa = conductivity(Gmag2[pos], lambdasq);
            // multiply conductivity;
            Gx[pos] *= kappa;
            Gy[pos] *= kappa;
        }
        pos++; // skip x=xmax;
    }
    // divergence -> laplace;
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                sx += sobelX[k] * Gx[pos + nn[k]];
                sy += sobelY[k] * Gy[pos + nn[k]];
            }
            laplace[pos] = sx + sy;
        }
        pos++; // skip x=xmax;
    }
    return 1;
}
void EpSmooth(BYTE *image, int w, int h, BYTE *out) {
    const double lambda = 50.0;
    const double dt = 0.01; //time step;
    const int maxiter = 100;
    std::vector<double> fimage(w * h, 0);
    std::vector<double> laplace(w * h, 0);
    for (int k = fimage.size(); k-- > 0;) fimage[k] = image[k];
    for (int iter = maxiter; iter-- > 0;) {
        LaplaceImage(&fimage[0], w, h, lambda, &laplace[0]);
        // update image;
        for (int k = fimage.size(); k-- > 0; ) 
            fimage[k] += dt * laplace[k];   
    }
    // preparing output;
    for (int k = fimage.size(); k-- > 0; ) {
        int a = int(fimage[k]);
        out[k] = a > 255 ? 255: a < 0 ? 0: a;
    }
}

 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

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

Anisotropic Diffusion Filter  (0) 2024.02.23
Watershed Segmentation  (0) 2021.02.27
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (3) 2021.02.15
Fixed-Point Bicubic Interpolation  (1) 2021.01.19
Posted by helloktk
,

1차원 열방정식은 

$$ u_t (x,t)= \kappa u_{x x}(x,t), \qquad -\infty < x < \infty,~t>0 $$

으로 초기조건이 주어지거나 경계조건이 주어야 한다. $u(x,t)$가 해라면 $u(\lambda  x , \lambda^2 t)$ 또한 해가 됨음 쉽게 확인해 볼 수 있다. 즉, 열방정식은 scale 변경에 불변이 특성을 가진다. 이 특성을 이용해서 열방정식의 특별한 해를 구해보자. scale 불변성을 보면 해의 $x,t$에 대한 의존성은 $x/\sqrt{t}\to \lambda x/\sqrt{\lambda^2 t} = x/\sqrt{t}$의 조합형태로 나타나야 될 것으로 보인다.

$$ u(x,t) = h(x/\sqrt{t})$$

그런데 에너지 보존에 의해서 계의 전체 열에너지는 시간에 무관하게 일정한 값을 가져야 한다. 열에너지가 온도에 비례하므로 총 열에너지는 다음 적분값에 비례한다.

$$ \int_{-\infty}^\infty u(x,t) dx = \text{const}$$

따라서 $u = h(x/\sqrt{t})$와 같은 형태의 해는 에너지 보존을 만족시킬 수 없다.

$$ \int _{-\infty}^\infty h(x/\sqrt{t}) dx = \sqrt{t} \int_{-\infty}^\infty h(y) dy \propto \sqrt{t}$$

그렇지만 $u(x,t) dx $가 스케일 불변꼴이면 가능하므로 

$$u(x,t) = \frac{1}{\sqrt{t}} h( x / \sqrt{t})$$

와 같은 형태로 시도하자. 그러면 열방정식의 좌우변의 항이

$$ u_{xx} = t^{-3/2} h'' (x/\sqrt{t}), \qquad u_t = -\frac{1}{2} t^{-3/2}\Big[ h(x/\sqrt{t}) + h'(x/\sqrt{t}) x/\sqrt{t} \Big] $$

이므로 이를 열방정식에 대입하면 ($y= x/\sqrt{t}$)

\begin{align} \kappa h''(y) + \frac{1}{2} h(y) + \frac{1}{2} y h'(y)=0 \\ \to ~ h''(y) + \frac{1}{2\kappa} (y h(y) )' = 0 \\ \to~ h'(y) + \frac{y}{2\kappa } h(y) = \text{const} \to 0 \end{align}

여기서 $\text{const}=0$ 일 때 해만 고려한다. 이 경우

$$ h(y) = A e^{- y^2/4\kappa} \qquad \to \qquad u(x,t) = \frac{A}{\sqrt{t}} e^{- x^2 / 4\kappa t } $$

인 해를 얻는다. $\int _{-\infty}^\infty u(x,t) dx = 1$로 정규화시키면 적분상수 $A$는

$$ A= \frac{1}{\sqrt{4\pi \kappa } }$$

임을 알 수 있다. 이렇게 정규화된 해를 열방정식의 기본해, Green 함수 또는 heat kernel이라고 부른다.

$$K(x,t) = \frac{1}{\sqrt{4\pi\kappa t} }e^{- x^2/4\kappa t}$$

$K(x,t)$는 $t\to 0^+$일 때 $\delta(x)$와 같은 특성을 가지는데, 물리적으로는 $t=0$ 시점에 원점에 $\delta-$함수로 표현되는 열원을 놓았을 때 열의 확산을 기술하는 해이다. 열방정식이 선형이므로  $x-$축을 따라 분포된 각각의 $\{ \delta(x-x_i)| i=1,2,\dots\}$ 열원의 선형결합도 해가 된다. 

$$ u(x,t) = \sum_i K(x-x_i, t) f_i$$

또한 연속적으로 분포하는 열원의 경우에는 합을 적분형태로 변경하면 되므로 $t=0$일 때 온도분포가 $f(x)$로 주어진 경우 다음 표현이

$$ u(x,t) = \int_{-\infty}^\infty K(x-y, t)f(y)dy$$

열방정식을 만족시키고, 초기조건도 만족시킴을 확인할 수 있다.

$$ \lim_{t\to 0} \int_{-\infty}^\infty K(x-y, t) f(y)dy = \int_{-\infty}^\infty \delta(x-y) f(y) dy = f(y)$$

 

해의 몇가지 재미있는 특징을 살펴보자. $t=0$일 때 열원이 유한한 영역에서만 0이 아니더라도 $K(x,t>0)$가 0이 아니기 때문에 매우 짧은 시간에 모든 영역으로 열이 퍼져나간다. 즉 열이 퍼지는 속도가 무한대에 해당한다. 열방정식은 비상대론적인 방정식이다. 두 번째로 $t>0$일 때 $K(x-y,t)$가 $x=y$인 지점에서 벗어나는 경우에 매우 빨리 0으로 변하므로 해는 무한번 미분가능하다. 즉, 처음 주어진 온도 분포가 불연속적이더라도 바로 매끄러운 온도분포로 변하게 만드는 특성이 있다. 이 열방정식을 잡음이 들어있는 신호나 영상에 적용하면 smoothing 효과를 줄 수 있다. 그리고

\begin{align}u(\lambda x, \lambda^2 t) &= \int K(\lambda x - y, \lambda ^2 t) f(y)dy \\ &= \int K(\lambda x - \lambda y' , \lambda^2 t )  {\lambda } f(\lambda y') dy' \\ &= \int K(x-y',t)  {\lambda }f(\lambda y') dy' \end{align}

이므로 스케일링된 해는 초기 온도분포가 $\lambda f(\lambda y)$인 경우와 같음을 알 수 있다. 초기에 국소적으로 가열된 경우라면 $\lambda f(\lambda x)$의 분포는 거의 $\delta$ 함수에 접근할 것이므로 충분히 시간이 지나고 먼 거리에서는 온도분포는 초기조건에 무관하게 됨을 알 수 있고, 이는 smoothing 효과로 이해할 수 있다.

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

Generalized eigenvalues problem  (0) 2024.03.02
수치적으로 보다 정밀한 이차방정식의 해  (0) 2024.02.23
지구의 나이는?  (0) 2024.02.11
n 차원 공의 부피  (2) 2024.02.07
타원의 둘레길이  (1) 2024.02.06
Posted by helloktk
,

지구의 나이는?

Mathematics 2024. 2. 11. 22:16

초기의 지구는 매우 뜨거워 모든 암석들이 녹아있는 상태였을 것이다. 시간이 흐르면서 지구는 표면에서 차가운 우주(거의 0K이다)로 열이 배출되어 점차 식게 된다. 처음에서 표면과 내부의 온도가 같았지만 열의 배출과정에서 표면의 온도는 내부의 온도보다 낮게 되어 지구 내부에서 온도의 gradient가 만들어진다. 실험적으로 지표면 근처에서 이 기울기를 측정하면 지구 나이를 추정할 수 있다.

 

지구 내부에서 온도의 분포는 열방정식에 의해서 구할 수 있다. 우선 지구가 매우 크므로 평평한 판(그리고 매우 두꺼운)으로 근사를 하자. 그러면 열이 식는 과정은 1차원 열방정식으로 설명할 수 있다. 지구 내부 깊이 $x>0$인 지점에서 온도 분포 $u(x, t)$는 다음의 열방정식에 의해서 결정된다.

\begin{align}  u_t(x, t) = \kappa u_{xx} (x,t), \quad u(x, t=0)=f(x), \qquad x>0\end{align}

여기서 $f(x)$는 초기 지구 내부 온도로 일정한 값으로 생각할 수 있다. 지표면 밖(우주)은 무한히 큰 차가운 열원으로 생각할 수 있으므로 온도가 $0\rm K$으로 고정되어 있다고 볼 수 있다. 즉, $u(x,t)$은 모든 시간에서 다음의 경계조건을 만족해야 한다.

$$ u(x=0, t)=0$$

이 $x\ge 0$에서 열방정식은 전체 실수구간에서의 열방정식을 푸는 문제로 변환하면 보다 쉽게 해결할 수 있다. 단, $x=0$에서 경곗값을 유지해야 하므로 $t=0$일 때 초기값은 원 문제의 기함수로 확장을 쓰면 된다. 물리적으로 $t=0$ 일 때 $x=0$을 기준으로 $T$와 $-T$(물론 음의 절대온도는 없다!)를 가지는 영역을 열접촉을 시키고 난 후의 온도 분포를 찾는 문제로 환원시킨 것과 같다. 따라서 다음의 열방정식의 해 $v(x,t)$를 우선 고려하자.

\begin{align} v_t (x,t) = \kappa v_{xx}(x,t),     \quad v(x,0)= {f}_o(x)=\left\{ \begin{matrix} f(x) & x>0 \\ -f(-x) & x <0\end{matrix}\right. \end{align}

이 방정식의 해는 다음처럼 써진다.(Fourier 변환과 convolution을 사용하면 쉽게 구할 수 있다). $\hat{v}(s,t)$와 $\hat{f_o}(s)$를 각각 $v(x,t)$와 ${f_o}(x)$의 공간에 대한 Fourier 변환이라면 열방정식은 

$$\hat{v}_t (s,t)= -\kappa s^2 \hat{v}(s,t),\quad \hat{v}(s, 0) = \hat{f_o}(s)$$

$t$에 대해 적분하면 다음처럼 두 함수의 곱으로 표현된다.

\begin{align} \hat{v}(s,t) =  e^{- \kappa s^2 t /2} \hat{f_o}(s) = \hat{K}(s,t) \hat{f_o}(s)\end{align}

따라서 $v(x,t)$는 heat kernel $K(x,t)$와 ${f}_o(s)$의 convolution으로 주어지는데, $e^{-\kappa s^2 t}$의 역변환이 $K(x,t)=\frac{1}{\sqrt{4\pi\kappa t}} e^{- x^2/4\kappa t}$이므로 

\begin{align} v(x,t) &= (K * {f}_o)(x,t)\\ &=\frac{1}{\sqrt{ 4\pi\kappa t}} \int_{-\infty}^\infty e^{- (x-y)^2/4\kappa t } {f}_o(y) dy \\ &= \frac{1}{\sqrt{4\pi\kappa t}} \int_0^\infty \Big( e^{-(x-y)^2/4\kappa t} - e^{-(x+y)^2 / 4\kappa t} \Big) f(y)dy \end{align}

항상 $v(0,t) = 0$을 만족함을 쉽게 확인할 수 있다. 따라서 $x\ge 0$에서 $u(x,t) = v(x,t)$이다. 초기온도 분포가 일정한 ($f(x)=T$) 경우는 적분을 정리하면 error function으로 써짐을 알 수 있다.

$$u(x,t) = \frac{T}{\sqrt{\pi \kappa t}} \int_0^x e^{-s^2 /4\kappa t } ds = T\times  \text{Erf}(x/2\sqrt{\kappa t}) \quad (x\ge0)$$

현시점($t=t_0$)에서 $x=0$에서 온도 기울기 $\beta$를 구하면 

$$\beta= u_x (0, t_0) = \frac{T}{\sqrt{\pi \kappa t_0}}  \quad \rightarrow \quad t_0 = \frac{T^2}{\pi \kappa \beta^2} $$

지구 초기 온도로 $T = 3900 \rm K$,  지각의 평균 열확산도로 $\kappa = 1.2 \times 10^{-2} \rm cm^2/s$, 그리고 현시점의 지표면 근처에서 실험적으로 측정된 온도 기울기가 $\beta = 3.6\times 10^{-4} \rm K/cm$ 정도이므로

$$ t _0 =3.1131 \times 10^{15} \rm s =  9.7 \times 10^7 yr$$

이 추정에 의하면 지구나이가 1억 년이 안된다. 방사선 동위원소법에 의해서 지구 나이가 45억 년쯤 되는 것으로 나오는데 이와 큰 차이가 있다. 구형 지구를 고려하면 3차원 열방정식 문제가 되는데 이 방법을 쓰더라도 평평한 지구 모형을 쓴 결과와 큰 차이가 나지 않는다. 이는 지구가 생성된 이후 추가적인 열의 공급 없이 식는 과정만 고려한 모델인데 실제로는 지구내부에 있는 방사성 원소가 붕괴할 때 발생하는 열이 지속적으로 공급되고 있으므로 이에 대한 고려를 해야 좀 더 정확한 나이를 추정할 수 있을 것이다. 열방정식을 이용한 지구 나이 추정은 영국의 수리물리학자이자 공학자인 Kelvin 남작에 의해 시도되었다.

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90
Posted by helloktk
,

반지름이 $r$인 $n$차원 공의 부피 $V_n(r)$은 $r^n$에 비례하므로 

$$ V_n(r) = K_n r^n$$

처럼 쓸 수 있다. 계수 $K_n$는 반지름이 1인 공의 부피를 나타낸다. $K_n$을 구하기 위해서 다음의 $n$ 차원 gaussian 적분을 고려하자.

$$ I_n = \int_{-\infty}^\infty dx_1 \int_{-\infty}^\infty dx_2 \cdots \int_{-\infty}^\infty dx_n e^{-x_1^2 - x_2^2 - \cdots - x_n^2}$$

이 적분은 

$$I_n= \Big( \int_{-\infty}^\infty  e^{-x^2} dx\Big)^n = (\sqrt{\pi})^{n}$$

임을 쉽게 알 수 있다.  $n$ 공간에서 보면 적분인자가 중심에서 $(x_1, x_2,\cdots, x_n)$까지의 거리의 제곱 $x_1^2 + x_2^2 +\cdots + x_n^2$에만 의존하는 구대칭을 가진다. 따라서 $I_n$은 공간을 미소구각으로 나눈 후 $r$로 적분만 하면 된다. 이 경우 적분요소는 

$$ dx_1 dx_2 \cdots dx_n =( \text{구각면적}) \times dr= \frac{dV_n}{dr} dr = nK_n r^{n-1} dr$$ 

$$ I_n = \int_0^\infty e^{-r^2 }  nK _n r^{n-1} dr = \frac{n K_n }{2} \int_0^\infty e^{-t} t^{n/2-1} dt =K_n \frac{n}{2} \left(n/2-1\right)! = { K_ n } \times  (n/2)!$$

따라서 

$$ K_n = \frac{\pi^{n/2}}{(n/2)!}=\frac{\pi^{n/2}}{\Gamma(n/2+1)}$$

$(-1/2)!=\sqrt{\pi}$, $(1/2)!= \sqrt{\pi}/2$임을 고려하면

\begin{align} K_0 &= 1 \\ K_1 &= 2 \\  K_2 &= \pi=3.1415 \\ K_3 &= \frac{4\pi }{3}= 4.18879 \\ K_4 &= \frac{\pi^2}{2} = 4.9348\\ K_5 &= \frac{8\pi^2}{15} =  5.26379 \\ K_6 &= \frac{\pi^3}{6}  = 5.16771 \\ K_7 &= \frac{16\pi^3}{105}=4.72477 \\ K_ 8 &= \frac{\pi^4}{24}=4.05871 \\ K_{2n}& = \frac{\pi^n }{ n!},\qquad K_{2n+1} = \frac{2(2\pi)^n}{(2n+1)!!} \end{align} 

이 결과를 보면 차원이 증가하면 단위공의 부피값이 증가하다가 5차원에서 최대가 된 후 다시 감소한다.

 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

열방정식의 Green function  (0) 2024.02.12
지구의 나이는?  (0) 2024.02.11
타원의 둘레길이  (1) 2024.02.06
왜 땅속은 여름에 시원하고 겨울에는 따뜻한가?  (1) 2024.02.03
Rejection Sampling  (0) 2023.02.20
Posted by helloktk
,