diff --git a/CHANGELOG.md b/CHANGELOG.md
index 985ffc0..e6757f3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Documentation correction. Merged #12.
 
+### Removed
+
+- Deleted broken example script `scikit_image.py`.
+
 ## [1.4.1] (2020-09-28)
 
 ### Added
diff --git a/docs/index.rst b/docs/index.rst
index 51d7974..b467470 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -54,13 +54,13 @@ This section describes the general usage patterns of :py:mod:`pyefd`.
 The coefficients returned are the :math:`a_n`, :math:`b_n`, :math:`c_n` and :math:`d_n` of
 the following Fourier series representation of the shape.
 
-The coefficients returned are by default normalized so that they are
-rotation and size-invariant. This can be overridden by calling:
+The coefficients returned can be normalized so that they are
+rotation and size-invariant. This can be achieved by calling:
 
 .. code:: python
 
     from pyefd import elliptic_fourier_descriptors
-    coeffs = elliptic_fourier_descriptors(contour, order=10, normalize=False)
+    coeffs = elliptic_fourier_descriptors(contour, order=10, normalize=True)
 
 Normalization can also be done afterwards:
 
diff --git a/pyefd.py b/pyefd.py
index f328adf..3c64e43 100644
--- a/pyefd.py
+++ b/pyefd.py
@@ -34,20 +34,25 @@
     _range = range
 
 
-def elliptic_fourier_descriptors(contour, order=10, normalize=False):
+def elliptic_fourier_descriptors(
+    contour, order=10, normalize=False, return_transformation=False
+):
     """Calculate elliptical Fourier descriptors for a contour.
 
     :param numpy.ndarray contour: A contour array of size ``[M x 2]``.
     :param int order: The order of Fourier coefficients to calculate.
     :param bool normalize: If the coefficients should be normalized;
         see references for details.
-    :return: A ``[order x 4]`` array of Fourier coefficients.
-    :rtype: :py:class:`numpy.ndarray`
+    :param bool return_transformation: If the normalization parametres should be returned.
+        Default is ``False``.
+    :return: A ``[order x 4]`` array of Fourier coefficients and optionally the
+        transformation parametres ``scale``, ``psi_1`` (rotation) and ``theta_1`` (phase)
+    :rtype: ::py:class:`numpy.ndarray` or (:py:class:`numpy.ndarray`, (float, float, float))
 
     """
     dxy = np.diff(contour, axis=0)
     dt = np.sqrt((dxy ** 2).sum(axis=1))
-    t = np.concatenate([([0.]), np.cumsum(dt)])
+    t = np.concatenate([([0.0]), np.cumsum(dt)])
     T = t[-1]
 
     phi = (2 * np.pi * t) / T
@@ -74,12 +79,12 @@ def elliptic_fourier_descriptors(contour, order=10, normalize=False):
     )
 
     if normalize:
-        coeffs = normalize_efd(coeffs)
+        coeffs = normalize_efd(coeffs, return_transformation=return_transformation)
 
     return coeffs
 
 
-def normalize_efd(coeffs, size_invariant=True):
+def normalize_efd(coeffs, size_invariant=True, return_transformation=False):
     """Normalizes an array of Fourier coefficients.
 
     See [#a]_ and [#b]_ for details.
@@ -87,8 +92,11 @@ def normalize_efd(coeffs, size_invariant=True):
     :param numpy.ndarray coeffs: A ``[n x 4]`` Fourier coefficient array.
     :param bool size_invariant: If size invariance normalizing should be done as well.
         Default is ``True``.
-    :return: The normalized ``[n x 4]`` Fourier coefficient array.
-    :rtype: :py:class:`numpy.ndarray`
+    :param bool return_transformation: If the normalization parametres should be returned.
+        Default is ``False``.
+    :return: The normalized ``[n x 4]`` Fourier coefficient array and optionally the
+        transformation parametres ``scale``, :math:`psi_1` (rotation) and :math:`theta_1` (phase)
+    :rtype: :py:class:`numpy.ndarray` or (:py:class:`numpy.ndarray`, (float, float, float))
 
     """
     # Make the coefficients have a zero phase shift from
@@ -136,11 +144,15 @@ def normalize_efd(coeffs, size_invariant=True):
             )
         ).flatten()
 
+    size = coeffs[0, 0]
     if size_invariant:
         # Obtain size-invariance by normalizing.
-        coeffs /= np.abs(coeffs[0, 0])
+        coeffs /= np.abs(size)
 
-    return coeffs
+    if return_transformation:
+        return coeffs, (size, psi_1, theta_1)
+    else:
+        return coeffs
 
 
 def calculate_dc_coefficients(contour):
@@ -153,7 +165,7 @@ def calculate_dc_coefficients(contour):
     """
     dxy = np.diff(contour, axis=0)
     dt = np.sqrt((dxy ** 2).sum(axis=1))
-    t = np.concatenate([([0.]), np.cumsum(dt)])
+    t = np.concatenate([([0.0]), np.cumsum(dt)])
     T = t[-1]
 
     xi = np.cumsum(dxy[:, 0]) - (dxy[:, 0] / dt) * t[1:]
@@ -199,7 +211,7 @@ def reconstruct_contour(coeffs, locus=(0, 0), num_points=300):
     return reconstruction
 
 
-def plot_efd(coeffs, locus=(0., 0.), image=None, contour=None, n=300):
+def plot_efd(coeffs, locus=(0.0, 0.0), image=None, contour=None, n=300):
     """Plot a ``[2 x (N / 2)]`` grid of successive truncations of the series.
 
     .. note::
diff --git a/setup.py b/setup.py
index f2a8864..ae89ee6 100644
--- a/setup.py
+++ b/setup.py
@@ -22,7 +22,7 @@
 EMAIL = "henrik.blidh@nedomkull.com"
 AUTHOR = "Henrik Blidh"
 REQUIRES_PYTHON = ">=2.7.10"
-VERSION = "1.4.1"
+VERSION = "1.5.0"
 
 REQUIRED = ["numpy>=1.7.0"]
 
diff --git a/tests.py b/tests.py
index b1a41b3..9d73d62 100644
--- a/tests.py
+++ b/tests.py
@@ -18,6 +18,7 @@
 
 import numpy as np
 from scipy.spatial.distance import directed_hausdorff
+from math import pi
 
 import pyefd
 
@@ -1008,13 +1009,38 @@ def test_normalizing_2():
     np.testing.assert_almost_equal(c[0, 2], 0.0, decimal=14)
 
 
+def test_normalizing_3():
+    # rotate and scale contour_1 by a known amount
+    theta = np.radians(30)
+    c, s = np.cos(theta), np.sin(theta)
+    R = np.array(((c, -s), (s, c))) * 1.5
+    contour_2 = np.transpose(np.dot(R, np.transpose(contour_1)))
+
+    c1, t1 = pyefd.elliptic_fourier_descriptors(
+        contour_1, normalize=True, return_transformation=True
+    )
+    c2, t2 = pyefd.elliptic_fourier_descriptors(
+        contour_2, normalize=True, return_transformation=True
+    )
+
+    # check if coefficients are equal (invariance)
+    np.testing.assert_almost_equal(c1, c2, decimal=12)
+    # check if normalization parametres match the initial transform
+    np.testing.assert_almost_equal(t1[0] * 1.5, t2[0], decimal=12)
+    np.testing.assert_almost_equal(
+        (t1[1] + np.radians(30)) % (2 * pi), t2[1], decimal=12
+    )
+
+
 def test_locus():
     locus = pyefd.calculate_dc_coefficients(contour_1)
     np.testing.assert_array_almost_equal(locus, np.mean(contour_1, axis=0), decimal=0)
 
 
 def test_reconstruct_simple_contour():
-    simple_polygon = np.array([[1., 1.], [0., 1.], [0., 0.], [1., 0.], [1., 1.]])
+    simple_polygon = np.array(
+        [[1.0, 1.0], [0.0, 1.0], [0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]
+    )
     number_of_points = simple_polygon.shape[0]
     locus = pyefd.calculate_dc_coefficients(simple_polygon)
     coeffs = pyefd.elliptic_fourier_descriptors(simple_polygon, order=30)
@@ -1049,7 +1075,7 @@ def for_loop_efd(contour, order=10, normalize=False):
         """
         dxy = np.diff(contour, axis=0)
         dt = np.sqrt((dxy ** 2).sum(axis=1))
-        t = np.concatenate([([0.]), np.cumsum(dt)])
+        t = np.concatenate([([0.0]), np.cumsum(dt)])
         T = t[-1]
 
         phi = (2 * np.pi * t) / T