Image rotation by FFT

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

$$ \begin{pmatrix} x' \\ y'\end{pmatrix} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\end{pmatrix}\begin{pmatrix} x\\ y\end{pmatrix} = \mathbf{R}_\theta \begin{pmatrix} x\\ y\end{pmatrix}$$

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

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

$$ \mathbf{S}_x =\begin{pmatrix} 1 & \tan (\theta/2) \\0 &1\end{pmatrix}, ~~\mathbf{S}_y =\begin{pmatrix}1 & 0 \\ -\sin \theta & 1\end{pmatrix}$$

 

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

\begin{gather} F(u, v)  =\iint   f(x, y) e^{-2\pi i (ux + vy)}dxdy={\cal F}[ f(x, y)] \\= {\cal F_x}\left[{\cal F_y} [f(x, y)]\right]={\cal F_y}\left[  {\cal  F_x}  [f(x, y)]\right] = F_1 (u) F_2 (v)\end{gather}

여기서 $\cal F_{x}, F_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) = {\cal F_x}[g(x, y)]= {\cal F_x}[f(x+ay, y)] = F_1(u, y) e^{-2\pi i ua y} =e^{-2\pi i  ua y}  {\cal F_x}[ f(x, y)]$$ 

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

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

 

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

중간단계: 

$x$ 방향 shear 변환 결과

추가적인 $y$ 방향 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];
    }
};