|
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
|
Borda