ENH Add ellipse_axes function

This was part of eccentricity computation, but can be useful on its own.

Also, added the shape functions to the ``mahotas.features`` namespace.

ChangeLog

@@ -1,5 +1,7 @@
 Version 1.1.0+
 	* Export locmax|locmin at the mahotas namespace level
+	* Break away ellipse_axes from eccentricity code as it can be useful on
+	its own
 
 Version 1.1.0 2014-02-12 by luispedro
 	* Better error checking

mahotas/features/__init__.py

@@ -1,12 +1,21 @@
+# -*- coding: utf-8 -*-
+# Copyright (C) 2012-2014, Luis Pedro Coelho <luis@luispedro.org>
+# vim: set ts=4 sts=4 sw=4 expandtab smartindent:
+#
+
 from .texture import haralick
 from .tas import tas, pftas
 from .zernike import zernike, zernike_moments
 from .lbp import lbp
+from .shape import roundness, eccentricity, ellipse_axes
 
 __all__ = [
+    'eccentricity',
+    'ellipse_axes',
     'haralick',
     'lbp',
     'pftas',
+    'roundness',
     'tas',
     'zernike',
     'zernike_moments',

mahotas/features/shape.py

@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# Copyright (C) 2012-2013, Luis Pedro Coelho <luis@luispedro.org>
+# Copyright (C) 2012-2014, Luis Pedro Coelho <luis@luispedro.org>
 # vim: set ts=4 sts=4 sw=4 expandtab smartindent:
 #
 from __future__ import division
@@ -10,7 +10,8 @@
 
 __all__ = [
     'roundness',
-    'eccentricity'
+    'eccentricity',
+    'ellipse_axes',
     ]
 
 def roundness(bw):
@@ -36,11 +37,13 @@ def roundness(bw):
     return float(perim)*perim/4./np.pi/area
 
 
-def eccentricity(bwimage):
-    """
-    ecc = eccentricity(bwimage)
+def ellipse_axes(bwimage):
+    ''' Parameters of the 'image ellipse'
 
-    Compute eccentricity
+    semimajor,semiminor = ellipse_axes(bwimage)
+
+    Returns the parameters of the constant intensity ellipse with the same mass
+    and second order moments as the original image.
 
     Parameters
     ----------
@@ -49,34 +52,53 @@ def eccentricity(bwimage):
 
     Returns
     -------
-    r : float
-        Eccentricity measure
-    """
+    semimajor : float
+    semiminor : float
+
+    References
+    ----------
+    Prokop, RJ, and Reeves, AP.  1992. CVGIP: Graphical Models and Image
+    Processing 54(5):438-460
+
+    '''
     from .moments import moments
     bwimage = _make_binary(bwimage)
 
     if not np.any(bwimage):
-        return 0
+        return 0.,0.
 
     cof = mh.center_of_mass(bwimage)
     hull_mu00 = moments(bwimage, 0, 0, cof)
     hull_mu11 = moments(bwimage, 1, 1, cof)
     hull_mu02 = moments(bwimage, 0, 2, cof)
     hull_mu20 = moments(bwimage, 2, 0, cof)
 
-    # Parameters of the 'image ellipse'
-    #   (the constant intensity ellipse with the same mass and
-    #   second order moments as the original image.)
-    #   From Prokop, RJ, and Reeves, AP.  1992. CVGIP: Graphical
-    #   Models and Image Processing 54(5):438-460
     semimajor = np.sqrt((2 * (hull_mu20 + hull_mu02 + \
                     np.sqrt((hull_mu20 - hull_mu02)**2 + \
                     4 * hull_mu11**2)))/hull_mu00)
 
     semiminor = np.sqrt((2 * (hull_mu20 + hull_mu02 - \
                     np.sqrt((hull_mu20 - hull_mu02)**2 + \
                     4 * hull_mu11**2)))/hull_mu00)
+    return semimajor, semiminor
+
+def eccentricity(bwimage):
+    """
+    ecc = eccentricity(bwimage)
+
+    Compute eccentricity
+
+    Parameters
+    ----------
+    bwimage : ndarray
+        Interpreted as a boolean image
 
+    Returns
+    -------
+    r : float
+        Eccentricity measure
+    """
+    semimajor, semiminor = ellipse_axes(bwimage)
     if semimajor == 0.:
         return 0.
     return  np.sqrt(semimajor**2 - semiminor**2) / semimajor

mahotas/tests/test_features_shape.py

@@ -28,3 +28,20 @@ def test_zeros():
     I = np.zeros((16,16))
     I[8:4:12] = 1
     assert eccentricity(I) == 0
+
+def test_ellipse_axes():
+    Y,X = np.mgrid[:1024,:1024]
+    Y = Y/1024.
+    X = X/1024.
+    im = ((2.*(Y - .5)**2 + (X - .5)**2) < .2)
+    maj,min = mh.features.ellipse_axes(im)
+    assert np.abs(2 - (maj/min)**2) < .01
+
+    maj2,min2 = mh.features.ellipse_axes(im.T)
+
+    assert np.abs(maj - maj2) < .001
+    assert np.abs(min - min2) < .001
+
+    im = (((Y - .5)**2 + (X - .5)**2) < .2)
+    maj,min = mh.features.ellipse_axes(im)
+    assert np.abs(1 - (maj/min)**2) < .01