Image rotation by FFT

평면을 한 점 $(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 변환 코드를 조금 손보면 쉽게 구현할 수 있다.

중간단계: 

$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];
    }
};