1
"""
2
Image registration tools
3

4
Copyright (C) 2016-2019 Jiri Borovec <jiri.borovec@fel.cvut.cz>
5
"""
6

7 30
import numpy as np
8 30
from skimage.transform import AffineTransform
9

10

11 30
def transform_points(points, matrix):
12
    """ transform points according to given transformation matrix
13

14
    :param ndarray points: set of points of shape (N, 2)
15
    :param ndarray matrix: transformation matrix of shape (3, 3)
16
    :return ndarray: warped points  of shape (N, 2)
17
    """
18 30
    points = np.array(points)
19
    # Pad the data with ones, so that our transformation can do translations
20 30
    pts_pad = np.hstack([points, np.ones((points.shape[0], 1))])
21 30
    points_warp = np.dot(pts_pad, matrix.T)
22 30
    return points_warp[:, :-1]
23

24

25 30
def estimate_affine_transform(points_0, points_1):
26
    """ estimate Affine transformations and warp particular points sets
27
    to the other coordinate frame
28

29
    :param ndarray points_0: set of points of shape (N, 2)
30
    :param ndarray points_1: set of points of shape (N, 2)
31
    :return tuple(ndarray,ndarray,ndarray,ndarray): transform. matrix & inverse
32
        and warped point sets
33

34
    >>> pts0 = np.array([[4., 116.], [4., 4.], [26., 4.], [26., 116.]], dtype=int)
35
    >>> pts1 = np.array([[61., 56.], [61., -56.], [39., -56.], [39., 56.]])
36
    >>> mx, mx_inv, pts0_w, pts1_w = estimate_affine_transform(pts0, pts1)
37
    >>> np.round(mx, 2) + 0.  # eliminate negative zeros
38
    array([[ -1.,   0.,  65.],
39
           [  0.,   1., -60.],
40
           [  0.,   0.,   1.]])
41
    >>> pts0_w
42
    array([[ 61.,  56.],
43
           [ 61., -56.],
44
           [ 39., -56.],
45
           [ 39.,  56.]])
46
    >>> pts1_w
47
    array([[   4.,  116.],
48
           [   4.,    4.],
49
           [  26.,    4.],
50
           [  26.,  116.]])
51
    """
52
    # SEE: https://stackoverflow.com/questions/20546182
53 30
    nb = min(len(points_0), len(points_1))
54
    # Pad the data with ones, so that our transformation can do translations
55 30
    x = np.hstack([points_0[:nb], np.ones((nb, 1))])
56 30
    y = np.hstack([points_1[:nb], np.ones((nb, 1))])
57

58
    # Solve the least squares problem X * A = Y to find our transform. matrix A
59 30
    matrix = np.linalg.lstsq(x, y, rcond=-1)[0].T
60 30
    matrix[-1, :] = [0, 0, 1]
61
    # invert the transformation matrix
62 30
    matrix_inv = np.linalg.pinv(matrix.T).T
63 30
    matrix_inv[-1, :] = [0, 0, 1]
64

65 30
    points_0_warp = transform_points(points_0, matrix)
66 30
    points_1_warp = transform_points(points_1, matrix_inv)
67

68 30
    return matrix, matrix_inv, points_0_warp, points_1_warp
69

70

71 30
def get_affine_components(matrix):
72
    """ get the main components of Affine transform
73

74
    :param ndarray matrix: affine transformation matrix for 2d
75
    :return dict: dictionary of float values
76

77
    >>> mtx = np.array([[ -0.95,   0.1,  65.], [  0.1,   0.95, -60.], [  0.,   0.,   1.]])
78
    >>> import  pandas as pd
79
    >>> aff = pd.Series(get_affine_components(mtx)).sort_index()
80
    >>> aff  # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
81
    rotation                           173.9...
82
    scale                    (0.95..., 0.95...)
83
    shear                              -3.14...
84
    translation                   (65.0, -60.0)
85
    dtype: object
86
    """
87 30
    aff = AffineTransform(matrix)
88 30
    norm_rotation = norm_angle(np.rad2deg(aff.rotation), deg=True)
89 30
    comp = {
90
        'rotation': float(norm_rotation),
91
        'translation': tuple(aff.translation.tolist()),
92
        'scale': aff.scale,
93
        'shear': aff.shear,
94
    }
95 30
    return comp
96

97

98 30
def norm_angle(angle, deg=True):
99
    """ normalize angles to be in range -half -> +half
100

101
    :param float angle: input angle
102
    :param bool deg: using degree or radian
103
    :return float:
104

105
    >>> norm_angle(60)
106
    60
107
    >>> norm_angle(200)
108
    -160
109
    >>> norm_angle(-540)
110
    180
111
    """
112 30
    unit = 180 if deg else np.pi
113 30
    while angle <= -unit:
114 30
        angle += 2 * unit
115 30
    while angle > unit:
116 30
        angle -= 2 * unit
117 30
    return angle

Read our documentation on viewing source code .

Loading