FFT를 이용한 영상의 미분

주어진 함수 $y(x)$를 구간 $0\le x < L$에서 일정한 간격으로 샘플링해서 얻은 $N$개의 data $\{y_n = y(nL/N)|n=0, 1,...,N-1\} $의 discrete Fourier transform(DFT) $Y_k$는

$$ Y_ k = \frac {1}{N} \sum _{n = 0}^{N-1} y_n e^{ -\frac {2\pi i}{N} n k}, $$

로 정의된다. 역변환(IDFT)은 $Y_k$에서 $y_n$을 얻는 과정으로 다음 식으로 계산된다:

$$ y_n = \sum_{k=0}^{N-1} Y_k e^{+\frac {2\pi i}{N}  k n}. $$

이산적인 data $y_n$만을 가지고 $x_n = nL/N$에서 미분계수를 구하려면 어떻게 해야 할까? 샘플링 점 사이 구간에서 값을 추정할 수 있도록 보간 함수를 만들어 사용해야 한다. 간단하게는 IDFT 식에서 $n$번째 샘플링 위치인 $x_n = n(L/N)$을 실수 $x$로 대체하여 연속함수를 얻는 방법을 생각할 수 있다. 그리고 IDFT 식에서 $ Y_k e^{+\frac {2\pi i}{N} k n}$ 대신에 추가적으로 위상 요소를 곱해서 $Y_k e^{+\frac {2\pi i}{N} (k + m_k N) n}, ~m_k =\text {integer}$ 로 바꾸더라도 동일한 원 데이터 $y_n$을 얻는다는 사실을 이용하도록 하자. 이 두 가지 관찰을 이용해서 만든 보간 함수는

$$ y(x) = \sum _{k=0}^{N-1} Y_k e^{+\frac{2\pi i}{L} (k +m_k N ) x}  $$

이다. 추가 위상이 $y_n$은 바꾸지 않지만 어떤 선택을 하는지에 따라 샘플링 사이 구간에서 다른 진동을 보이므로 미분계수가 영향을 받게 된다. 그러면 어떤 정수 $m_k$을 선택해야 하는가? 여러 가지 선택이 있을 수 있지만 전구간에서 도함수의 제곱 평균을 최소로 만드는 $m_k$을 선택하도록 하자: 

$$\underset {m_k} {\text {argmin} } ~ \frac {1}{L} \int_0^L | y'(x)|^2 dx  $$

그런데 $$\frac{1}{L} \int_0^L | y'(x)|^2 dx  = \Big(\frac {2\pi}{L} \Big)^2 ~\sum_{k=0}^{N-1}  (k+m_k N)^2|Y_k|^2  $$ 이므로 각 $k$에 대해서 $|k +m_k N|^2$이 최소가 되도록 선택해야 한다. $k<N/2$일 때는 $m_k=0$을 선택하면 되고, $k>N/2$일 때는 $m_k=-1$을 선택하면 된다. $N=\text{even}$일 때는 $m_{N/2}=0 \text{ or }1$이 같은 결과를 주는데, 두 기여를 균등하게 되도록 선택한다: 

$$ Y_{N/2} e^{\frac{2\pi i}{L} (N/2 + m_{N/2} N)x} \rightarrow \frac{1}{2} Y_{N/2}\left( e^{\frac{\pi i}{L} Nx} + e^{-\frac{\pi i}{L} Nx} \right) = Y_{N/2} \cos \Big( \frac{\pi}{L} Nx\Big) $$ 따라서 $m_k$ 선택은 정리하면 다음과 같다:

$$  m_k= \left\{\begin{array}{ll}   0, & 0\le k < N/2   \\ -1, &  N/2 <k <N  \end{array} \right.$$

따라서 주어진 샘플링 데이터 $\{y_n\}$의 DFT 계수 $\{ Y_n\}$을 구하면 샘플링 데이터를 보간하는 함수(Fourier interpolation)는

$$ y(x) = Y_0 + \sum_{0 <k <N/2}  \left( Y_k e^{+\frac {2\pi i}{L}  k x} + Y_{N-k} e^{ -\frac {2\pi i}{L}  k x} \right)  +  Y_{N/2} \cos\Big( \frac{\pi}{L}Nx\Big) $$

로 얻어진다(마지막 항은 $N$이 짝수일 때만 더해짐). 이 보간함수의 도함수는 

$$y'(x)= \sum_{0 <k <N/2} \frac {2\pi i }{L} k \Big(Y_k e^{+\frac {2\pi i}{L}  k x} - Y_{N-k} e^{ -\frac{2\pi i}{L}  k x}  \Big) - \frac{\pi}{L} N Y_{N/2} \sin \Big( \frac{\pi}{L} Nx\Big)$$

로 계산된다(마지막 항은 $N$이 짝수일 때만 더해짐). 샘플링 위치에서 미분계수는

$$ y'_n = y'(nL/N) = \sum_{0<k<N/2} \frac{2\pi i}{L} k \Big( Y_k e^{+\frac{2\pi i}{N} kn} - Y_{N-k} e^{ -\frac{2\pi i}{N} kn}  \Big) = \sum_{k=0}^{N-1} Y'_k e^{+\frac{2\pi i}{N} kn  }$$

처럼 쓸 수 있는데 이 식은 $Y'_k$의 IDFT 형태이다. 

따라서 DFT를 이용해서 미분계수를 구하는 전체적인 알고리즘은

1. DFT로 $Y_k$을 구하고,

2. $Y'_k$는 $$  Y_k' = \left\{  \begin {array}{ll}    \frac {2\pi i}{L} k  \times Y_k & 0 \le k <N/2\\ \frac{2\pi i}{L} (k-N) \times Y_k &  k > N/2 \\ 0 & k=N/2 ~\text {when } N=\text {even} \end {array} \right.  $$

3. $Y'_k$에 IDFT를 해서 $y'_n$을 구한다.

 

위에서 구한 보간 함수를 이용하면 이차 미분 값도 DFT를 이용해서 구할 수 있다.

 

아래 그림은 6개의 주파수 $f_k = 2,5,9,11,21,29\text{Hz}$가 섞인 신호 함수의 미분을 FFT를 이용해서 구한 결과를 정확한 도함수(green)와 함께 보여준다.

$$f(x) = \sum_{k=0}^{5} (k+1) \cos (2\pi f_k x)$$

$0\le x < 1=L$ 구간에서 일정한 간격으로 샘플링한 $N=64$개의 데이터를 사용하였다. 

void fftgrad1d() {
    const double TWOPI = 8.0 * atan(1.0);
    const int samples = 64;
    double re[samples], im[samples];
    const int nfreqs = 6;
    const double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
    for ( int i = 0; i < samples; i++ ) {
        double ti = double(i) / samples;
        double signal = 0;
        for ( int k = 0; k < nfreqs; k++ )
            signal += (k + 1) * cos ( TWOPI * freq[k] * ti);
        re[i] = signal;
        im[i] = 0.;
    }
    // int fft1d(int dir, int nn, double x[], double y[]);
    fft1d(1, samples, re, im);
    const double norm = TWOPI / 1.0;
    for (int i = 0; i < samples; i++) {
        double rv = re[i];
        double iv = im[i];
        double fac = norm * double (i > samples / 2 ? i - samples: i);  
        if (i == samples / 2) fac = 0;
        re[i] = - fac * iv;
        im[i] = fac * rv;
    }
    //invers FFT: re[] stores derivatives.
    fft1d(-1, samples, re, im);
}

이차원에서 DFT를 이용해서 gradient을 구하는 경우는 복소함수 $f = Re(x, y) + i Im(x, y)$의 gradient가

$$\Big(\frac{\partial}{ \partial x} +i \frac{\partial}{\partial y} \Big) (Re + i Im)= i\Big( \frac{\partial Re}{\partial y} + \frac{\partial Im}{\partial x} \Big) - \Big(- \frac{\partial Re}{\partial x} + \frac{\partial Im}{\partial y}\Big)  $$

임을 이용하면 된다.

영상(실수 함수)의 IDFT 식에 적용하면 미분 연산자가 우변에 전체적으로 $i$를 곱하는 효과를 주게 되고, 우변식의 실수부와 허수가 각각 영상의 $x$-방향/ $y$-방향 미분에 해당한다.

spectrum

Canny edge

 

void FFTGrad(CRaster& raster, CRaster& out) {
    const double TWOPI = 8.0 * atan(1.0);
    CSize sz = raster.GetSize();
    int nx = sz.cx, ny = sz.cy;
    int stride = nx << 1;
    std::vector<double> data(stride * sz.cy);
    for (int y = 0; y < ny; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        double *re = &data[y * stride];
        double *im = re + nx;
        for (int x = 0; x < nx; x++) {
            re[x] = *p++;
            im[x] = 0;
        }
    }
    // forward FFT(1): https://kipl.tistory.com/22
    fft2d(&data[0], nx, ny, 1);
    double normx = TWOPI / double(nx);
    double normy = TWOPI / double(ny);
    int hnx = nx >> 1;
    int hny = ny >> 1;
    for (int y = 0; y < ny; y++) {
        double *re = &data[y * stride];
        double *im = re + nx;
        for (int x = 0; x < nx; x++) {
            double rv = re[x];
            double iv = im[x];
            double xx = normx * double(x > hnx ? x - nx: x);
            double yy = normy * double(y > hny ? y - ny: y);
            if (x == hnx) xx = 0;  // N = even;
            if (y == hny) yy = 0;
            re[x] =  yy * rv + xx * iv;
            im[x] = -xx * rv + yy * iv;
        }
    }
    /* inverse FFT of {Y'_k} */
    fft2d(&data[0], nx, ny, -1);
    std::vector<double> mag(nx * ny);
    double magmax = 0, magmin = 1.e+20;
    for (int y = 0, k = 0; y < ny; y++) {
        double *re = &data[y * stride];
        double *im = re + nx;
        for (int x = 0; x < nx; x++) {
            double m = hypot(re[x], im[x]);
            if (m > magmax) magmax = m;
            if (m < magmin) magmin = m;
            mag[k++] = m;
        }
    }
    // stretch the magnitude;
    for (int y = 0, k = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            double a = 255 * (mag[k++] - magmin) / (magmax - magmin);
            out.SetPixel(x, y, 255-int(a));
        }
}