눈은 작지만 유한한 크기를 가지고 있으므로 물체의 한 점에서는 나오는 여러 광선을 받아들인다. 눈에 들어오는 광선들의 연장선이 만나는 지점에 물체의 겉보기 위치가 된다. 그리고 같이 광선 1과 2의 경로를 분석해보자. 눈이 작으므로 두 광선은 매우 가까이 있다( $\delta\ll 1$: 그림은 과장되게 그려진 것이다). 광선 1에 Snell의 법칙을 적용하면
$$ n \sin\theta_i = \sin \theta_r$$
이다. 이보다 조금 다른 각도 $(\theta_r + d\theta_r)$ 로 들어오는 광선 2에 대해서도 Snell의 법칙을 적용한 후 $((n \sin (\theta_i + d\theta_i) = \sin(\theta_r +d\theta_r) )$ 광선 1과의 차이를 구하면
$$ n \cos\theta_i d\theta_i =\cos \theta_r d \theta_r \quad \rightarrow \quad \frac{d\theta_i}{d\theta_r} = \frac{\cos\theta_r}{n\cos \theta_i}$$
그림에서 광선 1의 실제 경로에 대해서 물체의 수평 거리$(x)$와 깊이$(y)$의 관계는 $x=y\tan \theta_i$, 광선 2에 대해서는 $x+\delta = y \tan(\theta_i + d\theta_i)$이므로 둘의 차이를 계산하면
$$ \delta = y \sec^2 \theta_i d \theta_i$$
마찬가지로 광선1 과 2의 겉보기 경로에 대해 겉보기 수평 거리$(x_a)$와 깊이$(y_a)$의 관계를 정리하면
여기서 $\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 변환에 의해서 스펙트럼이 위상만 변화고 크기에는 변화가 생기지 않음도 알 수 있다.
아래의 코드는 이미지의 중앙을 기준으로 회전을 한 경우를 보여준다: $\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;
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: padding 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];
}
};