Python - Contingency tables using sparse coordinate matrices(Applications of sparse matrices: image transformations)

Elegant Scipy

開発環境

• OS X El Capitan - Apple (OS)
• Emacs (Text Editor)
• Python 3.5 (プログラミング言語)

Elegant Scipy (Juan Nunez-iglesias (著)、Stéfan Van Der Walt (著)、Harriet Dashnow (著)、 O'Reilly Media)のChapter 5.(Contingency tables using sparse coordinate matrices)、Applications of sparse matrices: image transformations の Exercise(No. 3282) を取り組んでみる。

Exercise(No. 3282)

コード(Emacs)

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import scipy
from skimage import data
import matplotlib.pyplot as plt
from itertools import product


def transform(image, angle):
    c = np.cos(np.deg2rad(angle))
    s = np.sin(np.deg2rad(angle))

    H = np.array([[c, -s, 0],
                  [s, c, 0],
                  [0, 0, 1]])

    x, y = np.array(image.shape) / 2

    Htr = np.array([[1, 0, -x],
                    [0, 1, -y],
                    [0, 0, 1]])
    Htr1 = np.array([[1, 0, x],
                     [0, 1, y],
                     [0, 0, 1]])

    return Htr1 @ H @ Htr


def homography(tf, image_shape):
    H = np.linalg.inv(tf)
    m, n = image_shape
    row, col, values = [], [], []
    for sparse_op_row, (out_row, out_col) in (
            enumerate(product(range(m), range(n)))):
        in_row, in_col, in_abs = H @ [out_row, out_col, 1]
        in_row /= in_abs
        in_col /= in_abs
        if not 0 <= in_row < m - 1 or not 0 <= in_col < n - 1:
            continue
        top = int(np.floor(in_row))
        left = int(np.floor(in_col))
        t = in_row - top
        u = in_col - left
        row.extend([sparse_op_row] * 4)
        sparse_op_col = np.ravel_multi_index(([top, top, top + 1, top + 1],
                                              [left, left + 1, left, left + 1]),
                                             dims=(m, n))
        col.extend(sparse_op_col)
        values.extend([(1 - t) * (1 - u), (1 - t) * u, t * (1 - u), t * u])
    operator = scipy.sparse.coo_matrix(
        (values, (row, col)), shape=(m * n, m * n)).tocsr()
    return operator


def apply_transform(image, tf):
    return (tf @ image.flat).reshape(image.shape)

if __name__ == '__main__':
    plt.figure(figsize=(4,4))
    image = data.camera()
    angle = 30
    H = transform(image, angle)
    tf = homography(H, image.shape)
    out = apply_transform(image, tf)
    plt.imshow(out, cmap='gray')
    plt.savefig('rotated_camera.png')
    plt.show()

入出力結果(Terminal, IPython)

$ ./sample3.py
$