Numpyで画像の幾何学的変換を実装してみます。

まず、使用する画像を読み込みます。

import numpy as np
import matplotlib.pyplot as plt

original_image = plt.imread(image_name)
if np.issubdtype(original_image.dtype, np.integer):
  original_image = original_image / np.iinfo(original_image.dtype).max
plt.imshow(original_image)

画像の変換は同次座標を用いて3x3の行列で表現することができるので、画像に行列を適用する関数を定義します。
この関数内では、変換後の画像座標に対して逆行列をかけて原画像での位置を求め、二アレストネイバー法で補間しています。

def apply_matrix(image, matrix, background):
  if image.ndim == 2 and background.ndim == 3:
    image = np.tile(image.reshape(image.shape + (1,)), (1, 1, background.shape[2]))
  if image.ndim == 3 and background.ndim == 2:
    background = np.tile(background.reshape(background.shape + (1,)), (1, 1, image.shape[2]))

  coord = np.fromfunction(lambda y, x: np.dstack((x + 0.5, y + 0.5, np.ones(background.shape[:2]))), background.shape[:2])
  sample_coord = np.einsum('ijk,lk->ijl', coord, np.linalg.inv(matrix))

  sample_coord_x = np.round(sample_coord[:,:,0]).astype(int)
  sample_coord_y = np.round(sample_coord[:,:,1]).astype(int)
  clip_sample_coord_x = np.clip(sample_coord_x, 0, image.shape[1] - 1)
  clip_sample_coord_y = np.clip(sample_coord_y, 0, image.shape[0] - 1)

  samplable = (sample_coord_x >= 0) & (sample_coord_x < image.shape[1]) & (sample_coord_y >=0) & (sample_coord_y < image.shape[0])
  if image.ndim == 3:
    samplable = np.tile(samplable.reshape(samplable.shape + (1,)), (1, 1, image.shape[2]))

  return np.where(samplable, image[clip_sample_coord_y, clip_sample_coord_x], background)

この記事では画像の原点は左上になっています。

平行移動

x方向に$t_x$、y方向に$t_y$だけ平行移動させる変換を表す行列は次のようになります。

\begin{pmatrix}
   1 & 0 & t_x \\
   0 & 1 & t_y \\
   0 & 0 & 1 \\
\end{pmatrix}

def translate(x, y):
  return np.array([[1, 0, x],
                   [0, 1, y],
                   [0, 0, 1]])

平行移動を適用します。

background = np.full((512, 512, 3), 0.5)
translate_matrix = translate(100, 200)
translate_image = apply_matrix(original_image, translate_matrix, background)
plt.figure(figsize=(8, 8))
plt.imshow(translate_image)

拡大・縮小

x方向に$s_x$倍、y方向に$s_y$倍の拡大・縮小を表す行列は次のようになります。

\begin{pmatrix}
   s_x & 0 & 0 \\
   0 & s_y & 0 \\
   0 & 0 & 1 \\
\end{pmatrix}

def scale(x, y):
  return np.array([[x, 0, 0],
                   [0, y, 0],
                   [0, 0, 1]])

拡大・縮小を適用します。

background = np.full((512, 512, 3), 0.5)
scale_matrix = scale(1.3, 1.6)
scale_image = apply_matrix(original_image, scale_matrix, background)
plt.figure(figsize=(8, 8))
plt.imshow(scale_image)

回転

原点を中心に$\theta$だけ回転させる変換を表す行列は次のようになります。

\begin{pmatrix}
   \cos\theta & -\sin\theta & 0 \\
   \sin\theta & \cos\theta & 0 \\
   0 & 0 & 1 \\
\end{pmatrix}

def rotate(angle):
  c = np.cos(angle)
  s = np.sin(angle)
  return np.array([[c, -s, 0],
                   [s, c, 0],
                   [0, 0, 1]])

回転を適用します。

background = np.full((512, 512, 3), 0.5)
rotate_matrix = rotate(np.pi / 6.0)
rotate_image = apply_matrix(original_image, rotate_matrix, background)
plt.figure(figsize=(8, 8))
plt.imshow(rotate_image)

変換の合成

変換を合成するには行列の積をとります。

background = np.full((512, 512, 3), 0.5)

matrix1 = translate(-original_image.shape[0] / 2, -original_image.shape[1] / 2)
matrix2 = rotate(np.pi / 2)
matrix3 = translate(original_image.shape[0] / 2, original_image.shape[1] / 2)
matrix4 = scale(2.0, 2.0)

composite_matrix = matrix4 @ matrix3 @ matrix2 @ matrix1

composite_image = apply_matrix(original_image, composite_matrix, background)
plt.figure(figsize=(8, 8))
plt.imshow(composite_image)

実装したコードはGoogle Colaboratoryに置いてあります。