|
| 1 | +#!/usr/bin/env python |
| 2 | +# -*- coding: utf-8 -*- |
| 3 | +""" |
| 4 | +:mod:`efd` |
| 5 | +================== |
| 6 | +
|
| 7 | +Created by hbldh <henrik.blidh@nedomkull.com> |
| 8 | +Created on 2016-01-30 |
| 9 | +
|
| 10 | +A Python implementation of the method described in [1]_ for |
| 11 | +calculating Fourier coefficients for characterizing |
| 12 | +closed contours. |
| 13 | +
|
| 14 | +References: |
| 15 | +----------- |
| 16 | +
|
| 17 | +.. [1] F. P. Kuhl and C. R. Giardina, “Elliptic Fourier Features of a |
| 18 | + Closed Contour," Computer Vision, Graphics and Image Processing, |
| 19 | + Vol. 18, pp. 236-258, 1982. |
| 20 | +
|
| 21 | +.. [2] Oivind Due Trier, Anil K. Jain and Torfinn Taxt, “Feature Extraction |
| 22 | + Methods for Character Recognition - A Survey”, Pattern Recognition |
| 23 | + Vol. 29, No.4, pp. 641-662 (1996) |
| 24 | +
|
| 25 | +""" |
| 26 | + |
| 27 | +from __future__ import division |
| 28 | +from __future__ import print_function |
| 29 | +from __future__ import unicode_literals |
| 30 | +from __future__ import absolute_import |
| 31 | + |
| 32 | +import numpy as np |
| 33 | + |
| 34 | +try: |
| 35 | + _range = xrange |
| 36 | +except NameError: |
| 37 | + _range = range |
| 38 | + |
| 39 | + |
| 40 | +def elliptical_fourier_descriptors(contour, order=10, normalize=False): |
| 41 | + """Calculate elliptical Fourier descriptors for a contour. |
| 42 | +
|
| 43 | + :param contour: A contour array of size [M x 2]. |
| 44 | + :type contour: :py:class:`numpy.ndarray` |
| 45 | + :param order: The order of Fourier coefficients to calculate. |
| 46 | + :type order: int |
| 47 | + :param normalize: If the coefficients should be normalized |
| 48 | + as is described in [1]_ and [2]_. |
| 49 | + :type normalize: bool |
| 50 | + :return: A [n x 4] array of Fourier coefficients. |
| 51 | + :rtype: :py:class:`numpy.ndarray` |
| 52 | +
|
| 53 | + """ |
| 54 | + dxy = np.diff(contour, axis=0) |
| 55 | + dt = np.sqrt((dxy ** 2).sum(axis=1)) |
| 56 | + t = np.cumsum(dt) |
| 57 | + T = t[-1] |
| 58 | + |
| 59 | + phi = np.concatenate([([0., ]), (2 * np.pi * t) / T]) |
| 60 | + |
| 61 | + coeffs = np.zeros((order, 4)) |
| 62 | + for n in _range(1, order + 1): |
| 63 | + const = T / (2 * n * n * np.pi * np.pi) |
| 64 | + phi_n = phi * n |
| 65 | + d_cos_phi_n = np.cos(phi_n[1:]) - np.cos(phi_n[:-1]) |
| 66 | + d_sin_phi_n = np.sin(phi_n[1:]) - np.sin(phi_n[:-1]) |
| 67 | + a_n = const * np.sum((dxy[:, 0] / dt) * d_cos_phi_n) |
| 68 | + b_n = const * np.sum((dxy[:, 0] / dt) * d_sin_phi_n) |
| 69 | + c_n = const * np.sum((dxy[:, 1] / dt) * d_cos_phi_n) |
| 70 | + d_n = const * np.sum((dxy[:, 1] / dt) * d_sin_phi_n) |
| 71 | + coeffs[n - 1, :] = a_n, b_n, c_n, d_n |
| 72 | + |
| 73 | + if normalize: |
| 74 | + coeffs = normalize_efd(coeffs) |
| 75 | + |
| 76 | + return coeffs |
| 77 | + |
| 78 | + |
| 79 | +def normalize_efd(coeffs, size_invariant=True): |
| 80 | + """Normalizes an array of Fourier coefficients. |
| 81 | +
|
| 82 | + See details in [1]_ or [2]_. |
| 83 | +
|
| 84 | + :param coeffs: A [n x 4] Fourier coefficient array. |
| 85 | + :type coeffs: :py:class:`numpy.ndarray` |
| 86 | + :return: The normalized [n x 4] Fourier coefficient array. |
| 87 | + :rtype: :py:class:`numpy.ndarray` |
| 88 | +
|
| 89 | + """ |
| 90 | + # Make the coefficients have a zero phase shift from |
| 91 | + # the first major axis. Theta_1 is that shift angle. |
| 92 | + theta_1 = 0.5 * np.arctan2( |
| 93 | + 2 * ((coeffs[0, 0] * coeffs[0, 1]) + (coeffs[0, 2] * coeffs[0, 3])), |
| 94 | + ((coeffs[0, 0] ** 2) - (coeffs[0, 1] ** 2) + (coeffs[0, 2] ** 2) - (coeffs[0, 3] ** 2))) |
| 95 | + # Rotate all coefficients by theta_1. |
| 96 | + for n in _range(1, coeffs.shape[0] + 1): |
| 97 | + coeffs[n - 1, :] = np.dot( |
| 98 | + np.array([[coeffs[n - 1, 0], coeffs[n - 1, 1]], |
| 99 | + [coeffs[n - 1, 2], coeffs[n - 1, 3]]]), |
| 100 | + np.array([[np.cos(n * theta_1), -np.sin(n * theta_1)], |
| 101 | + [np.sin(n * theta_1), np.cos(n * theta_1)]])).flatten() |
| 102 | + |
| 103 | + # Make the coefficients rotation invariant by rotating so that |
| 104 | + # the semi-major axis is parallel to the x-axis. |
| 105 | + psi_1 = np.arctan2(coeffs[0, 2], coeffs[0, 0]) |
| 106 | + psi_rotation_matrix = np.array([[np.cos(psi_1), np.sin(psi_1)], |
| 107 | + [-np.sin(psi_1), np.cos(psi_1)]]) |
| 108 | + # Rotate all coefficients by -psi_1. |
| 109 | + for n in _range(1, coeffs.shape[0] + 1): |
| 110 | + coeffs[n - 1, :] = psi_rotation_matrix.dot( |
| 111 | + np.array([[coeffs[n - 1, 0], coeffs[n - 1, 1]], |
| 112 | + [coeffs[n - 1, 2], coeffs[n - 1, 3]]])).flatten() |
| 113 | + |
| 114 | + if size_invariant: |
| 115 | + # Obtain size-invariance by normalizing. |
| 116 | + coeffs /= np.abs(coeffs[0, 0]) |
| 117 | + |
| 118 | + return coeffs |
| 119 | + |
| 120 | + |
| 121 | +def plot_efd(contour, coeffs, locus=(0., 0.), n=300): |
| 122 | + """Plot a [2 x (n/2)] grid of successive truncations of the series. |
| 123 | +
|
| 124 | + :param coeffs: [n x 4] Fourier coefficient array. |
| 125 | + :type coeffs: :py:class:`numpy.ndarray` |
| 126 | + :param locus: The A_0 and C_0 elliptic locus in [1]_ and [2]_. |
| 127 | + :type locus: list, tuple or :py:class:`numpy.ndarray` |
| 128 | + :param n: Number of points to use for plotting of Fourier series. |
| 129 | + :type n: int |
| 130 | +
|
| 131 | + """ |
| 132 | + try: |
| 133 | + import matplotlib.pyplot as plt |
| 134 | + except ImportError: |
| 135 | + print("Cannot plot: matplotlib was not installed.") |
| 136 | + return |
| 137 | + |
| 138 | + N = coeffs.shape[0] |
| 139 | + t = np.linspace(0, 1.0, n) |
| 140 | + xt = np.ones((n,)) * locus[0] |
| 141 | + yt = np.ones((n,)) * locus[1] |
| 142 | + |
| 143 | + ax = plt.subplot2grid((3, N // 2), (0, 0), colspan=N//2) |
| 144 | + ax.imshow(contour, plt.cm.gray) |
| 145 | + for n in _range(coeffs.shape[0]): |
| 146 | + xt += (coeffs[n, 0] * np.cos(2 * (n + 1) * np.pi * t)) + \ |
| 147 | + (coeffs[n, 1] * np.sin(2 * (n + 1) * np.pi * t)) |
| 148 | + yt += (coeffs[n, 2] * np.cos(2 * (n + 1) * np.pi * t)) + \ |
| 149 | + (coeffs[n, 3] * np.sin(2 * (n + 1) * np.pi * t)) |
| 150 | + ax = plt.subplot2grid((3, N // 2), (n // (N // 2) + 1, n % (N // 2))) |
| 151 | + ax.set_title(str(n + 1)) |
| 152 | + ax.plot(yt, -xt, 'r') |
| 153 | + |
| 154 | + plt.show() |
| 155 | + |
| 156 | + |
0 commit comments