Geometry & Recognition

평면을 한 점 $(x,y)$이 원점을 축으로 하는 회전에 의해서 $(x^{\prime },y^{\prime })$ 위치로 옮겼다면 두 점 사이의 관계는 

$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)={\mathbf{R}}_{\theta }\left(\begin{array}{c}x\\ y\end{array}\right)$$

로 표현된다. 여기서 $\theta$는 시계방향으로 회전한 각이다. 이 회전 변환은 3번의 shear 변환의 곱으로 다음처럼 쓸 수 있다.

$${\mathbf{R}}_{\theta }={\mathbf{S}}_x.{\mathbf{S}}_y.{\mathbf{S}}_x$$

$${\mathbf{S}}_x=\left(\begin{array}{cc}1& \tan (\theta /2)\\ 0& 1\end{array}\right),{\mathbf{S}}_y=\left(\begin{array}{cc}1& 0\\ -\sin \theta & 1\end{array}\right)$$

이제 Fourier변환과 shear 변환 사이의 관계를 알아보자. 일반적인 2차원 Fourier 변환은

$$\begin{array}{c}F(u,v)=\iint f(x,y)e^{-2\pi i(ux+vy)}dxdy=\mathcal{F}[f(x,y)]\\ ={\mathcal{F}}_{\mathcal{x}}[{\mathcal{F}}_{\mathcal{y}}[f(x,y)]]={\mathcal{F}}_{\mathcal{y}}[{\mathcal{F}}_{\mathcal{x}}[f(x,y)]]=F_1(u)F_2(v)\end{array}$$

여기서 ${\mathcal{F}}_{\mathcal{x}}\mathcal{,}{\mathcal{F}}_{\mathcal{y}}$는 각각 $x$, $y$-방향 Fourier 변환이고 순서에 상관이 없다.

이제 $f(x,y)$에 $x$-방향으로 $a$ 만큼 shear 변환을 한 결과가 $g(x,y)$일 때,

$$g(x,y)=f(x+ay,y)$$

로 쓸 수 있고, 여기에 $x$ 방향으로 Fourier 변환을 하면

$$G_1(u,y)={\mathcal{F}}_{\mathcal{x}}[g(x,y)]={\mathcal{F}}_{\mathcal{x}}[f(x+ay,y)]=F_1(u,y)e^{-2\pi iuay}=e^{-2\pi iuay}{\mathcal{F}}_{\mathcal{x}}[f(x,y)]$$ 

임을 알 수 있다. 즉, shear 변환된 함수의 Fourier 변환은 원함수의 Fourier 변환에 shear 파라미터에 의존하는 위상 요소 $e^{-2\pi iuay}$만 곱해주면 쉽게 구할 수 있음을 알 수 있다. 또한 shear 변환에 의해서 스펙트럼이 위상만 변화고 크기에는 변화가 생기지 않음도 알 수 있다.

참고: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.190.4473&rep=rep1&type=pdf

아래의 코드는 이미지의 중앙을 기준으로 회전을 한 경우를 보여준다: $\theta ={60}^{\circ }$. 중간 단계에서 artifact를 없애기 위해 원본 이미지의 가장자리를 $0$으로 채워 두 배 크기로 만든 후 수행한다. 임의의 지점에 대한 회전도 shear 변환 코드를 조금 손보면 쉽게 구현할 수 있다.

중간단계: 

int fft1d(int dir, int nn, double x[], double y[]); /* https://kipl.tistory.com/22 */
#define HORIZONTAL 0
#define VERTICAL   1
/* perform shear transform along "axis" direction; */
void shearImage(int axis, double shear,
                int width, int height, double *input, double *output) 
{
    int size, halfsz, othersize, halfothersz;
    const double TWOPI = 8.0 * atan(1.);
    /* size is the size of the image in the direction of the axis of shear. 
    ** othersize the size of the image in the orthogonal direction 
    */
    if (axis == HORIZONTAL) {
        halfsz = (size = width) >> 1;
        halfothersz = (othersize = height) >> 1;
    }
    else {
        halfsz = (size = height) >> 1;
        halfothersz = (othersize = width) >> 1;
    }

    std::vector<double> re_sig(size, 0);
    std::vector<double> im_sig(size, 0);
    for (int i = 0; i < othersize; i++){
        if (axis == HORIZONTAL)
            memcpy(&re_sig[0],  &input[i * size], size * sizeof(double));
        else
            for (int j = 0; j < size; j++) re_sig[j] = output[j * othersize + i];

        std::fill(im_sig.begin(), im_sig.end(), 0);  // im: fill 0;
        fft1d(1, size, &re_sig[0], &im_sig[0]);
        double shift = - (i - halfothersz - .5) * shear  * TWOPI;
        // No need to touch j = 0;
        for (int j = 1; j < halfsz; j++) {
            double cc = cos(j * shift / size);
            double ss = sin(j * shift / size);
            double reval = re_sig[j];
            double imval = im_sig[j];
            re_sig[j] = cc * reval - ss * imval;
            im_sig[j] = ss * reval + cc * imval;
            re_sig[size - j] =  re_sig[j];
            im_sig[size - j] = -im_sig[j];
        }
        re_sig[halfsz] = cos(shift/2) * re_sig[halfsz] - sin(shift/2) * im_sig[halfsz];
        im_sig[halfsz] = 0.;

        fft1d(-1, size, &re_sig[0], &im_sig[0]);
        if (axis == HORIZONTAL)
            memcpy(&output[i * size], &re_sig[0], size * sizeof(double));
        else
            for (int j = 0; j < size; j++) output[j * othersize + i] = re_sig[j];
    }
};

주어진 함수 $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}nk},$$

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

$$y_n=\sum _{k=0}^{N-1}{Y}_k{e}^{+\frac{2\pi i}{N}kn}.$$

이산적인 data $y_n$만을 가지고 $x_n=nL/N$에서 미분계수를 구하려면 어떻게 해야 할까? 샘플링 점 사이 구간에서 값을 추정할 수 있도록 보간 함수를 만들어 사용해야 한다. 간단하게는 IDFT 식에서 $n$번째 샘플링 위치인 $x_n=n(L/N)$을 실수 $x$로 대체하여 연속함수를 얻는 방법을 생각할 수 있다. 그리고 IDFT 식에서 $Y_k{e}^{+\frac{2\pi i}{N}kn}$ 대신에 추가적으로 위상 요소를 곱해서 $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^{\prime }(x){|}^{2}dx$$

그런데 $$\frac{1}{L}{\int }_0^L|y^{\prime }(x){|}^{2}dx=(\frac{2\pi }{L}{)}^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}\to \frac{1}{2}{Y}_{N/2}(e^{\frac{\pi i}{L}Nx}+e^{-\frac{\pi i}{L}Nx})=Y_{N/2}\cos (\frac{\pi }{L}Nx)$$ 따라서 $m_k$ 선택은 정리하면 다음과 같다:

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

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

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

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

$$y^{\prime }(x)=\sum _{0<k<N/2}\frac{2\pi i}{L}k(Y_k{e}^{+\frac{2\pi i}{L}kx}-Y_{N-k}{e}^{-\frac{2\pi i}{L}kx})-\frac{\pi }{L}N{Y}_{N/2}\sin (\frac{\pi }{L}Nx)$$

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

$$y_n^{\prime }=y^{\prime }(nL/N)=\sum _{0<k<N/2}\frac{2\pi i}{L}k(Y_k{e}^{+\frac{2\pi i}{N}kn}-Y_{N-k}{e}^{-\frac{2\pi i}{N}kn})=\sum _{k=0}^{N-1}{Y}_k^{\prime }{e}^{+\frac{2\pi i}{N}kn}$$

처럼 쓸 수 있는데 이 식은 $Y_k^{\prime }$의 IDFT 형태이다. 

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

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

2. $Y_k^{\prime }$는 $$Y_k^{\prime }=\{\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}$$

3. $Y_k^{\prime }$에 IDFT를 해서 $y_n^{\prime }$을 구한다.

위에서 구한 보간 함수를 이용하면 이차 미분 값도 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)+iIm(x,y)$의 gradient가

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

임을 이용하면 된다.

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

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));
        }
}

최소 자승법 문제는 추정치와 실제의 측정치와의 차이를 최소로 만드는 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,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}(\underset{\text{nonzero}}{\underbrace{w_0,w_1,...}},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{array}{c}{\mathbf{v}}_k^T.\mathbf{x}={\mathbf{u}}_k^T.\mathbf{b}/w_k,w_k\ne 0,\\ {\mathbf{v}}_k^T.\mathbf{x}=0,w_k=0,\\ \to \mathbf{x}=\sum _{w_k\ne 0}({\mathbf{u}}_k^T.\mathbf{b}/w_k){\mathbf{v}}_k,\end{array}$$

이어서 위의 해를 구할 수 있다. 이 해는 $|\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;
}