diff --git a/Licenses/scikit-image.txt b/Licenses/scikit-image.txt new file mode 100644 index 0000000..108806d --- /dev/null +++ b/Licenses/scikit-image.txt @@ -0,0 +1,81 @@ +Copyright (C) 2019, the scikit-image team +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the + distribution. + 3. Neither the name of skimage nor the names of its contributors may be + used to endorse or promote products derived from this software without + specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR +IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, +INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. + +skimage/_shared/version_requirements.py:_check_version + + Copyright (c) 2013 The IPython Development Team + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + * Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +skimage/_shared/version_requirements.py:is_installed: + + Original Copyright (C) 2009-2011 Pierre Raybaut + + Permission is hereby granted, free of charge, to any person obtaining + a copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to + permit persons to whom the Software is furnished to do so, subject to + the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/facetracker.spec b/facetracker.spec new file mode 100644 index 0000000..a62d869 --- /dev/null +++ b/facetracker.spec @@ -0,0 +1,53 @@ +# -*- mode: python ; coding: utf-8 -*- + +block_cipher = None + + +a = Analysis(['facetracker.py'], + pathex=['C:\\OpenSeeFaceBuild'], + binaries=[('dshowcapture/dshowcapture_x86.dll', '.'), ('dshowcapture/dshowcapture_x64.dll', '.'), ('escapi/escapi_x86.dll', '.'), ('escapi/escapi_x64.dll', '.'), ('run.bat', '.')], + datas=[], + hiddenimports=[], + hookspath=[], + runtime_hooks=[], + excludes=['matplotlib', 'mpl-data', 'PyInstaller', 'pywt', 'skimage', 'scipy', 'pyinstaller'], + win_no_prefer_redirects=False, + win_private_assemblies=False, + cipher=block_cipher, + noarchive=False) + +remove_bin = [] +for bin in a.binaries: + if bin[0].startswith("opencv_video") or bin[0].startswith("PyInstaller"): + remove_bin.append(bin[0]) + print(bin) +remove_dat = [] +for bin in a.datas: + if bin[0].startswith("opencv_video") or bin[0].startswith("PyInstaller"): + remove_dat.append(bin[0]) + print("data ", bin) + +a.binaries = [x for x in a.binaries if not x[0] in remove_bin] +a.datas = [x for x in a.datas if not x[0] in remove_dat] + +pyz = PYZ(a.pure, a.zipped_data, + cipher=block_cipher) +exe = EXE(pyz, + a.scripts, + [], + exclude_binaries=True, + name='facetracker', + debug=False, + bootloader_ignore_signals=False, + strip=False, + upx=True, + console=True ) + +coll = COLLECT(exe, + a.binaries, + a.zipfiles, + a.datas, + strip=False, + upx=True, + upx_exclude=[], + name='facetracker') diff --git a/make_exe.bat b/make_exe.bat index cc5118b..e3a8b4d 100644 --- a/make_exe.bat +++ b/make_exe.bat @@ -1,15 +1,17 @@ -REM pyinstaller --onefile ^ +REM pyinstaller facetracker.spec --onefile ^ REM --add-binary dshowcapture/dshowcapture_x86.dll;. ^ REM --add-binary dshowcapture/dshowcapture_x64.dll;. ^ REM --add-binary escapi/escapi_x86.dll;. ^ REM --add-binary escapi/escapi_x64.dll;. ^ -REM --add-binary run.bat;. ^ -REM facetracker.py +REM --add-binary run.bat;. -pyinstaller --onedir ^ +pyinstaller facetracker.spec --onedir ^ --add-binary dshowcapture/dshowcapture_x86.dll;. ^ --add-binary dshowcapture/dshowcapture_x64.dll;. ^ --add-binary escapi/escapi_x86.dll;. ^ --add-binary escapi/escapi_x64.dll;. ^ - --add-binary run.bat;. ^ - facetracker.py \ No newline at end of file + --add-binary run.bat;. + +del dist\facetracker\cv2\opencv_videoio_ffmpeg420_64.dll + +pause \ No newline at end of file diff --git a/similaritytransform.py b/similaritytransform.py new file mode 100644 index 0000000..4174787 --- /dev/null +++ b/similaritytransform.py @@ -0,0 +1,1370 @@ +# From scikit-image +import math +import numpy as np +import sys +import textwrap + +def get_bound_method_class(m): + """Return the class for a bound method. + + """ + return m.im_class if sys.version < '3' else m.__self__.__class__ + +def safe_as_int(val, atol=1e-3): + """ + Attempt to safely cast values to integer format. + + Parameters + ---------- + val : scalar or iterable of scalars + Number or container of numbers which are intended to be interpreted as + integers, e.g., for indexing purposes, but which may not carry integer + type. + atol : float + Absolute tolerance away from nearest integer to consider values in + ``val`` functionally integers. + + Returns + ------- + val_int : NumPy scalar or ndarray of dtype `np.int64` + Returns the input value(s) coerced to dtype `np.int64` assuming all + were within ``atol`` of the nearest integer. + + Notes + ----- + This operation calculates ``val`` modulo 1, which returns the mantissa of + all values. Then all mantissas greater than 0.5 are subtracted from one. + Finally, the absolute tolerance from zero is calculated. If it is less + than ``atol`` for all value(s) in ``val``, they are rounded and returned + in an integer array. Or, if ``val`` was a scalar, a NumPy scalar type is + returned. + + If any value(s) are outside the specified tolerance, an informative error + is raised. + + Examples + -------- + >>> safe_as_int(7.0) + 7 + + >>> safe_as_int([9, 4, 2.9999999999]) + array([9, 4, 3]) + + >>> safe_as_int(53.1) + Traceback (most recent call last): + ... + ValueError: Integer argument required but received 53.1, check inputs. + + >>> safe_as_int(53.01, atol=0.01) + 53 + + """ + mod = np.asarray(val) % 1 # Extract mantissa + + # Check for and subtract any mod values > 0.5 from 1 + if mod.ndim == 0: # Scalar input, cannot be indexed + if mod > 0.5: + mod = 1 - mod + else: # Iterable input, now ndarray + mod[mod > 0.5] = 1 - mod[mod > 0.5] # Test on each side of nearest int + + try: + np.testing.assert_allclose(mod, 0, atol=atol) + except AssertionError: + raise ValueError("Integer argument required but received " + "{0}, check inputs.".format(val)) + + return np.round(val).astype(np.int64) + +def _to_ndimage_mode(mode): + """Convert from `numpy.pad` mode name to the corresponding ndimage mode.""" + mode_translation_dict = dict(edge='nearest', symmetric='reflect', + reflect='mirror') + if mode in mode_translation_dict: + mode = mode_translation_dict[mode] + return mode + + +def _center_and_normalize_points(points): + """Center and normalize image points. + + The points are transformed in a two-step procedure that is expressed + as a transformation matrix. The matrix of the resulting points is usually + better conditioned than the matrix of the original points. + + Center the image points, such that the new coordinate system has its + origin at the centroid of the image points. + + Normalize the image points, such that the mean distance from the points + to the origin of the coordinate system is sqrt(2). + + Parameters + ---------- + points : (N, 2) array + The coordinates of the image points. + + Returns + ------- + matrix : (3, 3) array + The transformation matrix to obtain the new points. + new_points : (N, 2) array + The transformed image points. + + References + ---------- + .. [1] Hartley, Richard I. "In defense of the eight-point algorithm." + Pattern Analysis and Machine Intelligence, IEEE Transactions on 19.6 + (1997): 580-593. + + """ + + centroid = np.mean(points, axis=0) + + rms = math.sqrt(np.sum((points - centroid) ** 2) / points.shape[0]) + + norm_factor = math.sqrt(2) / rms + + matrix = np.array([[norm_factor, 0, -norm_factor * centroid[0]], + [0, norm_factor, -norm_factor * centroid[1]], + [0, 0, 1]]) + + pointsh = np.row_stack([points.T, np.ones((points.shape[0]),)]) + + new_pointsh = (matrix @ pointsh).T + + new_points = new_pointsh[:, :2] + new_points[:, 0] /= new_pointsh[:, 2] + new_points[:, 1] /= new_pointsh[:, 2] + + return matrix, new_points + + +def _umeyama(src, dst, estimate_scale): + """Estimate N-D similarity transformation with or without scaling. + + Parameters + ---------- + src : (M, N) array + Source coordinates. + dst : (M, N) array + Destination coordinates. + estimate_scale : bool + Whether to estimate scaling factor. + + Returns + ------- + T : (N + 1, N + 1) + The homogeneous similarity transformation matrix. The matrix contains + NaN values only if the problem is not well-conditioned. + + References + ---------- + .. [1] "Least-squares estimation of transformation parameters between two + point patterns", Shinji Umeyama, PAMI 1991, :DOI:`10.1109/34.88573` + + """ + + num = src.shape[0] + dim = src.shape[1] + + # Compute mean of src and dst. + src_mean = src.mean(axis=0) + dst_mean = dst.mean(axis=0) + + # Subtract mean from src and dst. + src_demean = src - src_mean + dst_demean = dst - dst_mean + + # Eq. (38). + A = dst_demean.T @ src_demean / num + + # Eq. (39). + d = np.ones((dim,), dtype=np.double) + if np.linalg.det(A) < 0: + d[dim - 1] = -1 + + T = np.eye(dim + 1, dtype=np.double) + + U, S, V = np.linalg.svd(A) + + # Eq. (40) and (43). + rank = np.linalg.matrix_rank(A) + if rank == 0: + return np.nan * T + elif rank == dim - 1: + if np.linalg.det(U) * np.linalg.det(V) > 0: + T[:dim, :dim] = U @ V + else: + s = d[dim - 1] + d[dim - 1] = -1 + T[:dim, :dim] = U @ np.diag(d) @ V + d[dim - 1] = s + else: + T[:dim, :dim] = U @ np.diag(d) @ V + + if estimate_scale: + # Eq. (41) and (42). + scale = 1.0 / src_demean.var(axis=0).sum() * (S @ d) + else: + scale = 1.0 + + T[:dim, dim] = dst_mean - scale * (T[:dim, :dim] @ src_mean.T) + T[:dim, :dim] *= scale + + return T + + +class GeometricTransform(object): + """Base class for geometric transformations. + + """ + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : (N, 2) array + Source coordinates. + + Returns + ------- + coords : (N, 2) array + Destination coordinates. + + """ + raise NotImplementedError() + + def inverse(self, coords): + """Apply inverse transformation. + + Parameters + ---------- + coords : (N, 2) array + Destination coordinates. + + Returns + ------- + coords : (N, 2) array + Source coordinates. + + """ + raise NotImplementedError() + + def residuals(self, src, dst): + """Determine residuals of transformed destination coordinates. + + For each transformed source coordinate the euclidean distance to the + respective destination coordinate is determined. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + residuals : (N, ) array + Residual for coordinate. + + """ + return np.sqrt(np.sum((self(src) - dst)**2, axis=1)) + + def __add__(self, other): + """Combine this transformation with another. + + """ + raise NotImplementedError() + + +class FundamentalMatrixTransform(GeometricTransform): + """Fundamental matrix transformation. + + The fundamental matrix relates corresponding points between a pair of + uncalibrated images. The matrix transforms homogeneous image points in one + image to epipolar lines in the other image. + + The fundamental matrix is only defined for a pair of moving images. In the + case of pure rotation or planar scenes, the homography describes the + geometric relation between two images (`ProjectiveTransform`). If the + intrinsic calibration of the images is known, the essential matrix describes + the metric relation between the two images (`EssentialMatrixTransform`). + + References + ---------- + .. [1] Hartley, Richard, and Andrew Zisserman. Multiple view geometry in + computer vision. Cambridge university press, 2003. + + Parameters + ---------- + matrix : (3, 3) array, optional + Fundamental matrix. + + Attributes + ---------- + params : (3, 3) array + Fundamental matrix. + + """ + + def __init__(self, matrix=None): + if matrix is None: + # default to an identity transform + matrix = np.eye(3) + if matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix") + self.params = matrix + + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : (N, 2) array + Source coordinates. + + Returns + ------- + coords : (N, 3) array + Epipolar lines in the destination image. + + """ + coords_homogeneous = np.column_stack([coords, np.ones(coords.shape[0])]) + return coords_homogeneous @ self.params.T + + def inverse(self, coords): + """Apply inverse transformation. + + Parameters + ---------- + coords : (N, 2) array + Destination coordinates. + + Returns + ------- + coords : (N, 3) array + Epipolar lines in the source image. + + """ + coords_homogeneous = np.column_stack([coords, np.ones(coords.shape[0])]) + return coords_homogeneous @ self.params + + def _setup_constraint_matrix(self, src, dst): + """Setup and solve the homogeneous epipolar constraint matrix:: + + dst' * F * src = 0. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + F_normalized : (3, 3) array + The normalized solution to the homogeneous system. If the system + is not well-conditioned, this matrix contains NaNs. + src_matrix : (3, 3) array + The transformation matrix to obtain the normalized source + coordinates. + dst_matrix : (3, 3) array + The transformation matrix to obtain the normalized destination + coordinates. + + """ + if src.shape != dst.shape: + raise ValueError('src and dst shapes must be identical.') + if src.shape[0] < 8: + raise ValueError('src.shape[0] must be equal or larger than 8.') + + # Center and normalize image points for better numerical stability. + try: + src_matrix, src = _center_and_normalize_points(src) + dst_matrix, dst = _center_and_normalize_points(dst) + except ZeroDivisionError: + self.params = np.full((3, 3), np.nan) + return 3 * [np.full((3, 3), np.nan)] + + # Setup homogeneous linear equation as dst' * F * src = 0. + A = np.ones((src.shape[0], 9)) + A[:, :2] = src + A[:, :3] *= dst[:, 0, np.newaxis] + A[:, 3:5] = src + A[:, 3:6] *= dst[:, 1, np.newaxis] + A[:, 6:8] = src + + # Solve for the nullspace of the constraint matrix. + _, _, V = np.linalg.svd(A) + F_normalized = V[-1, :].reshape(3, 3) + + return F_normalized, src_matrix, dst_matrix + + def estimate(self, src, dst): + """Estimate fundamental matrix using 8-point algorithm. + + The 8-point algorithm requires at least 8 corresponding point pairs for + a well-conditioned solution, otherwise the over-determined solution is + estimated. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + + F_normalized, src_matrix, dst_matrix = \ + self._setup_constraint_matrix(src, dst) + + # Enforcing the internal constraint that two singular values must be + # non-zero and one must be zero. + U, S, V = np.linalg.svd(F_normalized) + S[2] = 0 + F = U @ np.diag(S) @ V + + self.params = dst_matrix.T @ F @ src_matrix + + return True + + def residuals(self, src, dst): + """Compute the Sampson distance. + + The Sampson distance is the first approximation to the geometric error. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + residuals : (N, ) array + Sampson distance. + + """ + src_homogeneous = np.column_stack([src, np.ones(src.shape[0])]) + dst_homogeneous = np.column_stack([dst, np.ones(dst.shape[0])]) + + F_src = self.params @ src_homogeneous.T + Ft_dst = self.params.T @ dst_homogeneous.T + + dst_F_src = np.sum(dst_homogeneous * F_src.T, axis=1) + + return np.abs(dst_F_src) / np.sqrt(F_src[0] ** 2 + F_src[1] ** 2 + + Ft_dst[0] ** 2 + Ft_dst[1] ** 2) + + +class EssentialMatrixTransform(FundamentalMatrixTransform): + """Essential matrix transformation. + + The essential matrix relates corresponding points between a pair of + calibrated images. The matrix transforms normalized, homogeneous image + points in one image to epipolar lines in the other image. + + The essential matrix is only defined for a pair of moving images capturing a + non-planar scene. In the case of pure rotation or planar scenes, the + homography describes the geometric relation between two images + (`ProjectiveTransform`). If the intrinsic calibration of the images is + unknown, the fundamental matrix describes the projective relation between + the two images (`FundamentalMatrixTransform`). + + References + ---------- + .. [1] Hartley, Richard, and Andrew Zisserman. Multiple view geometry in + computer vision. Cambridge university press, 2003. + + Parameters + ---------- + rotation : (3, 3) array, optional + Rotation matrix of the relative camera motion. + translation : (3, 1) array, optional + Translation vector of the relative camera motion. The vector must + have unit length. + matrix : (3, 3) array, optional + Essential matrix. + + Attributes + ---------- + params : (3, 3) array + Essential matrix. + + """ + + def __init__(self, rotation=None, translation=None, matrix=None): + if rotation is not None: + if translation is None: + raise ValueError("Both rotation and translation required") + if rotation.shape != (3, 3): + raise ValueError("Invalid shape of rotation matrix") + if abs(np.linalg.det(rotation) - 1) > 1e-6: + raise ValueError("Rotation matrix must have unit determinant") + if translation.size != 3: + raise ValueError("Invalid shape of translation vector") + if abs(np.linalg.norm(translation) - 1) > 1e-6: + raise ValueError("Translation vector must have unit length") + # Matrix representation of the cross product for t. + t_x = np.array([0, -translation[2], translation[1], + translation[2], 0, -translation[0], + -translation[1], translation[0], 0]).reshape(3, 3) + self.params = t_x @ rotation + elif matrix is not None: + if matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix") + self.params = matrix + else: + # default to an identity transform + self.params = np.eye(3) + + def estimate(self, src, dst): + """Estimate essential matrix using 8-point algorithm. + + The 8-point algorithm requires at least 8 corresponding point pairs for + a well-conditioned solution, otherwise the over-determined solution is + estimated. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + + E_normalized, src_matrix, dst_matrix = \ + self._setup_constraint_matrix(src, dst) + + # Enforcing the internal constraint that two singular values must be + # equal and one must be zero. + U, S, V = np.linalg.svd(E_normalized) + S[0] = (S[0] + S[1]) / 2.0 + S[1] = S[0] + S[2] = 0 + E = U @ np.diag(S) @ V + + self.params = dst_matrix.T @ E @ src_matrix + + return True + + +class ProjectiveTransform(GeometricTransform): + r"""Projective transformation. + + Apply a projective transformation (homography) on coordinates. + + For each homogeneous coordinate :math:`\mathbf{x} = [x, y, 1]^T`, its + target position is calculated by multiplying with the given matrix, + :math:`H`, to give :math:`H \mathbf{x}`:: + + [[a0 a1 a2] + [b0 b1 b2] + [c0 c1 1 ]]. + + E.g., to rotate by theta degrees clockwise, the matrix should be:: + + [[cos(theta) -sin(theta) 0] + [sin(theta) cos(theta) 0] + [0 0 1]] + + or, to translate x by 10 and y by 20:: + + [[1 0 10] + [0 1 20] + [0 0 1 ]]. + + Parameters + ---------- + matrix : (3, 3) array, optional + Homogeneous transformation matrix. + + Attributes + ---------- + params : (3, 3) array + Homogeneous transformation matrix. + + """ + + _coeffs = range(8) + + def __init__(self, matrix=None): + if matrix is None: + # default to an identity transform + matrix = np.eye(3) + if matrix.shape != (3, 3): + raise ValueError("invalid shape of transformation matrix") + self.params = matrix + + @property + def _inv_matrix(self): + return np.linalg.inv(self.params) + + def _apply_mat(self, coords, matrix): + coords = np.array(coords, copy=False, ndmin=2) + + x, y = np.transpose(coords) + src = np.vstack((x, y, np.ones_like(x))) + dst = src.T @ matrix.T + + # below, we will divide by the last dimension of the homogeneous + # coordinate matrix. In order to avoid division by zero, + # we replace exact zeros in this column with a very small number. + dst[dst[:, 2] == 0, 2] = np.finfo(float).eps + # rescale to homogeneous coordinates + dst[:, :2] /= dst[:, 2:3] + + return dst[:, :2] + + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : (N, 2) array + Source coordinates. + + Returns + ------- + coords : (N, 2) array + Destination coordinates. + + """ + return self._apply_mat(coords, self.params) + + def inverse(self, coords): + """Apply inverse transformation. + + Parameters + ---------- + coords : (N, 2) array + Destination coordinates. + + Returns + ------- + coords : (N, 2) array + Source coordinates. + + """ + return self._apply_mat(coords, self._inv_matrix) + + def estimate(self, src, dst): + """Estimate the transformation from a set of corresponding points. + + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source and destination coordinates must match. + + The transformation is defined as:: + + X = (a0*x + a1*y + a2) / (c0*x + c1*y + 1) + Y = (b0*x + b1*y + b2) / (c0*x + c1*y + 1) + + These equations can be transformed to the following form:: + + 0 = a0*x + a1*y + a2 - c0*x*X - c1*y*X - X + 0 = b0*x + b1*y + b2 - c0*x*Y - c1*y*Y - Y + + which exist for each set of corresponding points, so we have a set of + N * 2 equations. The coefficients appear linearly so we can write + A x = 0, where:: + + A = [[x y 1 0 0 0 -x*X -y*X -X] + [0 0 0 x y 1 -x*Y -y*Y -Y] + ... + ... + ] + x.T = [a0 a1 a2 b0 b1 b2 c0 c1 c3] + + In case of total least-squares the solution of this homogeneous system + of equations is the right singular vector of A which corresponds to the + smallest singular value normed by the coefficient c3. + + In case of the affine transformation the coefficients c0 and c1 are 0. + Thus the system of equations is:: + + A = [[x y 1 0 0 0 -X] + [0 0 0 x y 1 -Y] + ... + ... + ] + x.T = [a0 a1 a2 b0 b1 b2 c3] + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + + try: + src_matrix, src = _center_and_normalize_points(src) + dst_matrix, dst = _center_and_normalize_points(dst) + except ZeroDivisionError: + self.params = np.nan * np.empty((3, 3)) + return False + + xs = src[:, 0] + ys = src[:, 1] + xd = dst[:, 0] + yd = dst[:, 1] + rows = src.shape[0] + + # params: a0, a1, a2, b0, b1, b2, c0, c1 + A = np.zeros((rows * 2, 9)) + A[:rows, 0] = xs + A[:rows, 1] = ys + A[:rows, 2] = 1 + A[:rows, 6] = - xd * xs + A[:rows, 7] = - xd * ys + A[rows:, 3] = xs + A[rows:, 4] = ys + A[rows:, 5] = 1 + A[rows:, 6] = - yd * xs + A[rows:, 7] = - yd * ys + A[:rows, 8] = xd + A[rows:, 8] = yd + + # Select relevant columns, depending on params + A = A[:, list(self._coeffs) + [8]] + + _, _, V = np.linalg.svd(A) + # if the last element of the vector corresponding to the smallest + # singular value is close to zero, this implies a degenerate case + # because it is a rank-defective transform, which would map points + # to a line rather than a plane. + if np.isclose(V[-1, -1], 0): + return False + + H = np.zeros((3, 3)) + # solution is right singular vector that corresponds to smallest + # singular value + H.flat[list(self._coeffs) + [8]] = - V[-1, :-1] / V[-1, -1] + H[2, 2] = 1 + + # De-center and de-normalize + H = np.linalg.inv(dst_matrix) @ H @ src_matrix + + self.params = H + + return True + + def __add__(self, other): + """Combine this transformation with another. + + """ + if isinstance(other, ProjectiveTransform): + # combination of the same types result in a transformation of this + # type again, otherwise use general projective transformation + if type(self) == type(other): + tform = self.__class__ + else: + tform = ProjectiveTransform + return tform(other.params @ self.params) + elif (hasattr(other, '__name__') + and other.__name__ == 'inverse' + and hasattr(get_bound_method_class(other), '_inv_matrix')): + return ProjectiveTransform(other.__self__._inv_matrix @ self.params) + else: + raise TypeError("Cannot combine transformations of differing " + "types.") + + def __nice__(self): + """common 'paramstr' used by __str__ and __repr__""" + npstring = np.array2string(self.params, separator=', ') + paramstr = 'matrix=\n' + textwrap.indent(npstring, ' ') + return paramstr + + def __repr__(self): + """Add standard repr formatting around a __nice__ string""" + paramstr = self.__nice__() + classname = self.__class__.__name__ + classstr = classname + return '<{}({}) at {}>'.format(classstr, paramstr, hex(id(self))) + + def __str__(self): + """Add standard str formatting around a __nice__ string""" + paramstr = self.__nice__() + classname = self.__class__.__name__ + classstr = classname + return '<{}({})>'.format(classstr, paramstr) + + +class AffineTransform(ProjectiveTransform): + """2D affine transformation. + + Has the following form:: + + X = a0*x + a1*y + a2 = + = sx*x*cos(rotation) - sy*y*sin(rotation + shear) + a2 + + Y = b0*x + b1*y + b2 = + = sx*x*sin(rotation) + sy*y*cos(rotation + shear) + b2 + + where ``sx`` and ``sy`` are scale factors in the x and y directions, + and the homogeneous transformation matrix is:: + + [[a0 a1 a2] + [b0 b1 b2] + [0 0 1]] + + Parameters + ---------- + matrix : (3, 3) array, optional + Homogeneous transformation matrix. + scale : {s as float or (sx, sy) as array, list or tuple}, optional + Scale factor(s). If a single value, it will be assigned to both + sx and sy. + + .. versionadded:: 0.17 + Added support for supplying a single scalar value. + rotation : float, optional + Rotation angle in counter-clockwise direction as radians. + shear : float, optional + Shear angle in counter-clockwise direction as radians. + translation : (tx, ty) as array, list or tuple, optional + Translation parameters. + + Attributes + ---------- + params : (3, 3) array + Homogeneous transformation matrix. + + """ + + _coeffs = range(6) + + def __init__(self, matrix=None, scale=None, rotation=None, shear=None, + translation=None): + params = any(param is not None + for param in (scale, rotation, shear, translation)) + + if params and matrix is not None: + raise ValueError("You cannot specify the transformation matrix and" + " the implicit parameters at the same time.") + elif matrix is not None: + if matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix.") + self.params = matrix + elif params: + if scale is None: + scale = (1, 1) + if rotation is None: + rotation = 0 + if shear is None: + shear = 0 + if translation is None: + translation = (0, 0) + + if np.isscalar(scale): + sx = sy = scale + else: + sx, sy = scale + + self.params = np.array([ + [sx * math.cos(rotation), -sy * math.sin(rotation + shear), 0], + [sx * math.sin(rotation), sy * math.cos(rotation + shear), 0], + [ 0, 0, 1] + ]) + self.params[0:2, 2] = translation + else: + # default to an identity transform + self.params = np.eye(3) + + @property + def scale(self): + sx = math.sqrt(self.params[0, 0] ** 2 + self.params[1, 0] ** 2) + sy = math.sqrt(self.params[0, 1] ** 2 + self.params[1, 1] ** 2) + return sx, sy + + @property + def rotation(self): + return math.atan2(self.params[1, 0], self.params[0, 0]) + + @property + def shear(self): + beta = math.atan2(- self.params[0, 1], self.params[1, 1]) + return beta - self.rotation + + @property + def translation(self): + return self.params[0:2, 2] + + + + + +class EuclideanTransform(ProjectiveTransform): + """2D Euclidean transformation. + + Has the following form:: + + X = a0 * x - b0 * y + a1 = + = x * cos(rotation) - y * sin(rotation) + a1 + + Y = b0 * x + a0 * y + b1 = + = x * sin(rotation) + y * cos(rotation) + b1 + + where the homogeneous transformation matrix is:: + + [[a0 b0 a1] + [b0 a0 b1] + [0 0 1]] + + The Euclidean transformation is a rigid transformation with rotation and + translation parameters. The similarity transformation extends the Euclidean + transformation with a single scaling factor. + + Parameters + ---------- + matrix : (3, 3) array, optional + Homogeneous transformation matrix. + rotation : float, optional + Rotation angle in counter-clockwise direction as radians. + translation : (tx, ty) as array, list or tuple, optional + x, y translation parameters. + + Attributes + ---------- + params : (3, 3) array + Homogeneous transformation matrix. + + """ + + def __init__(self, matrix=None, rotation=None, translation=None): + params = any(param is not None + for param in (rotation, translation)) + + if params and matrix is not None: + raise ValueError("You cannot specify the transformation matrix and" + " the implicit parameters at the same time.") + elif matrix is not None: + if matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix.") + self.params = matrix + elif params: + if rotation is None: + rotation = 0 + if translation is None: + translation = (0, 0) + + self.params = np.array([ + [math.cos(rotation), - math.sin(rotation), 0], + [math.sin(rotation), math.cos(rotation), 0], + [ 0, 0, 1] + ]) + self.params[0:2, 2] = translation + else: + # default to an identity transform + self.params = np.eye(3) + + def estimate(self, src, dst): + """Estimate the transformation from a set of corresponding points. + + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source and destination coordinates must match. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + + self.params = _umeyama(src, dst, False) + + return True + + @property + def rotation(self): + return math.atan2(self.params[1, 0], self.params[1, 1]) + + @property + def translation(self): + return self.params[0:2, 2] + + +class SimilarityTransform(EuclideanTransform): + """2D similarity transformation. + + Has the following form:: + + X = a0 * x - b0 * y + a1 = + = s * x * cos(rotation) - s * y * sin(rotation) + a1 + + Y = b0 * x + a0 * y + b1 = + = s * x * sin(rotation) + s * y * cos(rotation) + b1 + + where ``s`` is a scale factor and the homogeneous transformation matrix is:: + + [[a0 b0 a1] + [b0 a0 b1] + [0 0 1]] + + The similarity transformation extends the Euclidean transformation with a + single scaling factor in addition to the rotation and translation + parameters. + + Parameters + ---------- + matrix : (3, 3) array, optional + Homogeneous transformation matrix. + scale : float, optional + Scale factor. + rotation : float, optional + Rotation angle in counter-clockwise direction as radians. + translation : (tx, ty) as array, list or tuple, optional + x, y translation parameters. + + Attributes + ---------- + params : (3, 3) array + Homogeneous transformation matrix. + + """ + + def __init__(self, matrix=None, scale=None, rotation=None, + translation=None): + params = any(param is not None + for param in (scale, rotation, translation)) + + if params and matrix is not None: + raise ValueError("You cannot specify the transformation matrix and" + " the implicit parameters at the same time.") + elif matrix is not None: + if matrix.shape != (3, 3): + raise ValueError("Invalid shape of transformation matrix.") + self.params = matrix + elif params: + if scale is None: + scale = 1 + if rotation is None: + rotation = 0 + if translation is None: + translation = (0, 0) + + self.params = np.array([ + [math.cos(rotation), - math.sin(rotation), 0], + [math.sin(rotation), math.cos(rotation), 0], + [ 0, 0, 1] + ]) + self.params[0:2, 0:2] *= scale + self.params[0:2, 2] = translation + else: + # default to an identity transform + self.params = np.eye(3) + + def estimate(self, src, dst): + """Estimate the transformation from a set of corresponding points. + + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source and destination coordinates must match. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + + self.params = _umeyama(src, dst, True) + + return True + + @property + def scale(self): + # det = scale**(# of dimensions), therefore scale = det**(1/2) + return np.sqrt(np.linalg.det(self.params)) + + +class PolynomialTransform(GeometricTransform): + """2D polynomial transformation. + + Has the following form:: + + X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) + Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) + + Parameters + ---------- + params : (2, N) array, optional + Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So, + a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`. + + Attributes + ---------- + params : (2, N) array + Polynomial coefficients where `N * 2 = (order + 1) * (order + 2)`. So, + a_ji is defined in `params[0, :]` and b_ji in `params[1, :]`. + + """ + + def __init__(self, params=None): + if params is None: + # default to transformation which preserves original coordinates + params = np.array([[0, 1, 0], [0, 0, 1]]) + if params.shape[0] != 2: + raise ValueError("invalid shape of transformation parameters") + self.params = params + + def estimate(self, src, dst, order=2): + """Estimate the transformation from a set of corresponding points. + + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source and destination coordinates must match. + + The transformation is defined as:: + + X = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) + Y = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) + + These equations can be transformed to the following form:: + + 0 = sum[j=0:order]( sum[i=0:j]( a_ji * x**(j - i) * y**i )) - X + 0 = sum[j=0:order]( sum[i=0:j]( b_ji * x**(j - i) * y**i )) - Y + + which exist for each set of corresponding points, so we have a set of + N * 2 equations. The coefficients appear linearly so we can write + A x = 0, where:: + + A = [[1 x y x**2 x*y y**2 ... 0 ... 0 -X] + [0 ... 0 1 x y x**2 x*y y**2 -Y] + ... + ... + ] + x.T = [a00 a10 a11 a20 a21 a22 ... ann + b00 b10 b11 b20 b21 b22 ... bnn c3] + + In case of total least-squares the solution of this homogeneous system + of equations is the right singular vector of A which corresponds to the + smallest singular value normed by the coefficient c3. + + Parameters + ---------- + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. + order : int, optional + Polynomial order (number of coefficients is order + 1). + + Returns + ------- + success : bool + True, if model estimation succeeds. + + """ + xs = src[:, 0] + ys = src[:, 1] + xd = dst[:, 0] + yd = dst[:, 1] + rows = src.shape[0] + + # number of unknown polynomial coefficients + order = safe_as_int(order) + u = (order + 1) * (order + 2) + + A = np.zeros((rows * 2, u + 1)) + pidx = 0 + for j in range(order + 1): + for i in range(j + 1): + A[:rows, pidx] = xs ** (j - i) * ys ** i + A[rows:, pidx + u // 2] = xs ** (j - i) * ys ** i + pidx += 1 + + A[:rows, -1] = xd + A[rows:, -1] = yd + + _, _, V = np.linalg.svd(A) + + # solution is right singular vector that corresponds to smallest + # singular value + params = - V[-1, :-1] / V[-1, -1] + + self.params = params.reshape((2, u // 2)) + + return True + + def __call__(self, coords): + """Apply forward transformation. + + Parameters + ---------- + coords : (N, 2) array + source coordinates + + Returns + ------- + coords : (N, 2) array + Transformed coordinates. + + """ + x = coords[:, 0] + y = coords[:, 1] + u = len(self.params.ravel()) + # number of coefficients -> u = (order + 1) * (order + 2) + order = int((- 3 + math.sqrt(9 - 4 * (2 - u))) / 2) + dst = np.zeros(coords.shape) + + pidx = 0 + for j in range(order + 1): + for i in range(j + 1): + dst[:, 0] += self.params[0, pidx] * x ** (j - i) * y ** i + dst[:, 1] += self.params[1, pidx] * x ** (j - i) * y ** i + pidx += 1 + + return dst + + def inverse(self, coords): + raise Exception( + 'There is no explicit way to do the inverse polynomial ' + 'transformation. Instead, estimate the inverse transformation ' + 'parameters by exchanging source and destination coordinates,' + 'then apply the forward transformation.') + + +TRANSFORMS = { + 'euclidean': EuclideanTransform, + 'similarity': SimilarityTransform, + 'affine': AffineTransform, + 'projective': ProjectiveTransform, + 'fundamental': FundamentalMatrixTransform, + 'essential': EssentialMatrixTransform, + 'polynomial': PolynomialTransform, +} + + +def estimate_transform(ttype, src, dst, **kwargs): + """Estimate 2D geometric transformation parameters. + + You can determine the over-, well- and under-determined parameters + with the total least-squares method. + + Number of source and destination coordinates must match. + + Parameters + ---------- + ttype : {'euclidean', similarity', 'affine', 'piecewise-affine', \ + 'projective', 'polynomial'} + Type of transform. + kwargs : array or int + Function parameters (src, dst, n, angle):: + + NAME / TTYPE FUNCTION PARAMETERS + 'euclidean' `src, `dst` + 'similarity' `src, `dst` + 'affine' `src, `dst` + 'piecewise-affine' `src, `dst` + 'projective' `src, `dst` + 'polynomial' `src, `dst`, `order` (polynomial order, + default order is 2) + + Also see examples below. + + Returns + ------- + tform : :class:`GeometricTransform` + Transform object containing the transformation parameters and providing + access to forward and inverse transformation functions. + + Examples + -------- + >>> import numpy as np + >>> from skimage import transform + + >>> # estimate transformation parameters + >>> src = np.array([0, 0, 10, 10]).reshape((2, 2)) + >>> dst = np.array([12, 14, 1, -20]).reshape((2, 2)) + + >>> tform = transform.estimate_transform('similarity', src, dst) + + >>> np.allclose(tform.inverse(tform(src)), src) + True + + >>> # warp image using the estimated transformation + >>> from skimage import data + >>> image = data.camera() + + >>> warp(image, inverse_map=tform.inverse) # doctest: +SKIP + + >>> # create transformation with explicit parameters + >>> tform2 = transform.SimilarityTransform(scale=1.1, rotation=1, + ... translation=(10, 20)) + + >>> # unite transformations, applied in order from left to right + >>> tform3 = tform + tform2 + >>> np.allclose(tform3(src), tform2(tform(src))) + True + + """ + ttype = ttype.lower() + if ttype not in TRANSFORMS: + raise ValueError('the transformation type \'%s\' is not' + 'implemented' % ttype) + + tform = TRANSFORMS[ttype]() + tform.estimate(src, dst, **kwargs) + + return tform + + +def matrix_transform(coords, matrix): + """Apply 2D matrix transform. + + Parameters + ---------- + coords : (N, 2) array + x, y coordinates to transform + matrix : (3, 3) array + Homogeneous transformation matrix. + + Returns + ------- + coords : (N, 2) array + Transformed coordinates. + + """ + return ProjectiveTransform(matrix)(coords) diff --git a/tracker.py b/tracker.py index 2d22374..64fa892 100644 --- a/tracker.py +++ b/tracker.py @@ -7,7 +7,7 @@ import time import queue import threading import copy -from skimage.transform import SimilarityTransform +from similaritytransform import SimilarityTransform from retinaface import RetinaFaceDetector from remedian import remedian