Version 1.6.0

New example.
Bugfix for plotting.
Code cleaning.

Author: hbldh

CHANGELOG.md

@@ -4,6 +4,17 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [1.6.0] (2021-01-22)
+
+### Added
+
+- Added a demo for 3D surfaces with cylindrical symmetries. (`examples/example1.py`)
+
+### Fixes
+
+- Fixes incorrectly plotted curves when no `imshow` has been called.
+- Fixes ugly coefficient calculation code.
+
 ## [1.5.1] (2021-01-22)
 
 ### Added

examples/example1.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+
+An example showing how to use pyefd for fitting points along a closed curve.
+
+The demo curves are taken from:
+"Optimized Fourier representations for three-dimensional magnetic surfaces"
+by S. P. Hirshman and H. K. Meier,
+Phys. Fluids 28 (5) 1985 (https://doi.org/10.1063/1.864972)
+
+
+Created by Jonathan Schilling <jonathan.schilling@mail.de> on 2021-12-09.
+
+"""
+
+import numpy as np
+
+# hack to get the import from within module directory to work
+# see also: https://stackoverflow.com/a/16985066
+import sys
+import os
+SCRIPT_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..')
+if not SCRIPT_DIR in sys.path:
+    sys.path.append(SCRIPT_DIR)
+
+import pyefd
+
+# bean
+rCos = np.array([3.0, 1.042, 0.502, -0.0389])
+zSin = np.array([0.0, 1.339, 0.296, -0.0176])
+
+# diamond
+# rCos = np.array([0.5, 0.427, 0.0,  0.0732])
+# zSin = np.array([0.0, 0.427, 0.0, -0.0732])
+
+# D
+# rCos = np.array([3.0, 0.991,  0.136])
+# zSin = np.array([0.0, 1.409, -0.118])
+
+# belt
+# rCos = np.array([3.0, 0.453, 0.0, 0.0  ])
+# zSin = np.array([0.0, 0.6  , 0.0, 0.196])
+
+# ellipse
+# rCos = np.array([3.0, 1.0])
+# zSin = np.array([0.0, 3.0])
+
+# evaluate curve geometry given as Fourier coefficients above
+n = 300
+contour = np.zeros([n,2])
+theta = np.linspace(0.0, 2.0*np.pi, n)
+mMax = len(rCos)
+for m in range(mMax):
+    contour[:,0] += rCos[m] * np.cos(m*theta)
+    contour[:,1] += zSin[m] * np.sin(m*theta)
+
+# apply pyefd to get elliptic Fourier descriptors
+coeffs = pyefd.elliptic_fourier_descriptors(contour)
+a0, c0 = pyefd.calculate_dc_coefficients(contour)
+pyefd.plot_efd(coeffs, locus=(a0,c0), contour=contour)

pyefd.py

@@ -60,10 +60,11 @@ def elliptic_fourier_descriptors(
     orders = np.arange(1, order + 1)
     consts = T / (2 * orders * orders * np.pi * np.pi)
     phi = phi * orders.reshape((order, -1))
+
     d_cos_phi = np.cos(phi[:, 1:]) - np.cos(phi[:, :-1])
     d_sin_phi = np.sin(phi[:, 1:]) - np.sin(phi[:, :-1])
-    cos_phi = (dxy[:, 0] / dt) * d_cos_phi
-    a = consts * np.sum(cos_phi, axis=1)
+
+    a = consts * np.sum((dxy[:, 0] / dt) * d_cos_phi, axis=1)
     b = consts * np.sum((dxy[:, 0] / dt) * d_sin_phi, axis=1)
     c = consts * np.sum((dxy[:, 1] / dt) * d_cos_phi, axis=1)
     d = consts * np.sum((dxy[:, 1] / dt) * d_sin_phi, axis=1)
@@ -247,10 +248,23 @@ def plot_efd(coeffs, locus=(0.0, 0.0), image=None, contour=None, n=300):
         )
         ax = plt.subplot2grid((n_rows, N_half), (n // N_half, n % N_half))
         ax.set_title(str(n + 1))
-        if contour is not None:
-            ax.plot(contour[:, 1], contour[:, 0], "c--", linewidth=2)
-        ax.plot(yt, xt, "r", linewidth=2)
+
         if image is not None:
+            # A background image of shape [rows, cols] gets transposed
+            # by imshow so that the first dimension is vertical
+            # and the second dimension is horizontal.
+            # This implies swapping the x and y axes when plotting a curve.
+            if contour is not None:
+                ax.plot(contour[:, 1], contour[:, 0], "c--", linewidth=2)
+            ax.plot(yt, xt, "r", linewidth=2)
             ax.imshow(image, plt.cm.gray)
+        else:
+            # Without a background image, no transpose is implied.
+            # This case is useful when (x,y) point clouds
+            # without relation to an image are to be handled.
+            if contour is not None:
+                ax.plot(contour[:, 0], contour[:, 1], "c--", linewidth=2)
+            ax.plot(xt, yt, "r", linewidth=2)
+            ax.axis("equal")
 
     plt.show()

setup.py

@@ -22,7 +22,7 @@
 EMAIL = "henrik.blidh@nedomkull.com"
 AUTHOR = "Henrik Blidh"
 REQUIRES_PYTHON = ">=2.7.10"
-VERSION = "1.5.1"
+VERSION = "1.6.0"
 
 REQUIRED = ["numpy>=1.7.0"]