1 ```""" ``` 2 ```Image registration tools ``` 3 4 ```Copyright (C) 2016-2019 Jiri Borovec ``` 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 .