# Source code for dclab.features.inert_ratio

#!/usr/bin/python
# -*- coding: utf-8 -*-
"""Computation of inertia ratio from contour data"""
from __future__ import division, print_function, unicode_literals

import numpy as np
import scipy.spatial as ssp

[docs]def get_inert_ratio_cvx(cont):
"""Compute the inertia ratio of the convex hull of a contour

The inertia ratio is computed from the central second order of moments
along x (mu20) and y (mu02) via sqrt(mu20/mu02).

Parameters
----------
cont: ndarray or list of ndarrays of shape (N,2)
A 2D array that holds the contour of an event (in pixels)
e.g. obtained using mm.contour where  mm is an instance
of RTDCBase. The first and second columns of cont
correspond to the x- and y-coordinates of the contour.

Returns
-------
inert_ratio: float or ndarray of size N
The inertia ratio of the contour

Notes
-----
The contour moments mu20 and mu02 are computed the same way they
are computed in OpenCV's moments.cpp.

--------
get_inert_ratio_raw: Compute inertia ratio of a raw contour
"""
if isinstance(cont, np.ndarray):
# If cont is an array, it is not a list of contours,
# because contours can have different lengths.
cont = [cont]
ret_list = False
else:
ret_list = True

length = len(cont)

inert_ratio = np.zeros(length, dtype=float) * np.nan

for ii in range(length):
try:
chull = ssp.ConvexHull(cont[ii])
except ssp.qhull.QhullError:
pass
else:
hull = cont[ii][chull.vertices,:]
inert_ratio[ii] = get_inert_ratio_raw(hull)

if not ret_list:
inert_ratio = inert_ratio[0]

return inert_ratio

[docs]def get_inert_ratio_raw(cont):
"""Compute the inertia ratio of a contour

The inertia ratio is computed from the central second order of moments
along x (mu20) and y (mu02) via sqrt(mu20/mu02).

Parameters
----------
cont: ndarray or list of ndarrays of shape (N,2)
A 2D array that holds the contour of an event (in pixels)
e.g. obtained using mm.contour where  mm is an instance
of RTDCBase. The first and second columns of cont
correspond to the x- and y-coordinates of the contour.

Returns
-------
inert_ratio: float or ndarray of size N
The inertia ratio of the contour

Notes
-----
The contour moments mu20 and mu02 are computed the same way they
are computed in OpenCV's moments.cpp.

--------
get_inert_ratio_cvx: Compute inertia ratio of the convex hull of
a contour
"""
if isinstance(cont, np.ndarray):
# If cont is an array, it is not a list of contours,
# because contours can have different lengths.
cont = [cont]
ret_list = False
else:
ret_list = True

length = len(cont)

inert_ratio = np.zeros(length, dtype=float) * np.nan

for ii in range(length):
moments = cont_moments_cv(cont[ii])
if moments is not None:
inert_ratio[ii] = np.sqrt(moments["mu20"]/moments["mu02"])

if not ret_list:
inert_ratio = inert_ratio[0]

return inert_ratio

def cont_moments_cv(cont,
flt_epsilon=1.19209e-07,
dbl_epsilon=2.2204460492503131e-16):
"""Compute the moments of a contour

The moments are computed in the same way as they are computed
in OpenCV's contourMoments in moments.cpp.

Parameters
----------
cont: array of shape (N,2)
The contour for which to compute the moments.
flt_epsilon: float
The value of FLT_EPSILON in OpenCV/gcc.
dbl_epsilon: float
The value of DBL_EPSILON in OpenCV/gcc.

Returns
-------
moments: dict
A dictionary of moments. If the moment m00 is smaller
than half of flt_epsilon, None is returned.
"""
# Make sure we have no unsigned integers
if np.issubdtype(cont.dtype, np.unsignedinteger):
cont = cont.astype(np.int)

xi = cont[:,0]
yi = cont[:,1]

xi_1 = np.roll(xi, -1)
yi_1 = np.roll(yi, -1)

xi_12 = xi_1**2
yi_12 = yi_1**2

xi2 = xi**2
yi2 = yi**2

dxy = xi_1 * yi - xi * yi_1

xii_1 = xi_1 + xi
yii_1 = yi_1 + yi

a00 = np.sum(dxy)
a10 = np.sum(dxy * xii_1)
a01 = np.sum(dxy * yii_1)
a20 = np.sum(dxy * (xi_1 * xii_1 + xi2))
a11 = np.sum(dxy * (xi_1 * (yii_1 + yi_1) + xi * (yii_1 + yi)))
a02 = np.sum(dxy * (yi_1 * yii_1 + yi2))
a30 = np.sum(dxy * xii_1 * (xi_12 + xi2))
a03 = np.sum(dxy * yii_1 * (yi_12 + yi2))
a21 = np.sum(dxy * (xi_12 * (3 * yi_1 + yi) + 2 * xi * xi_1 * yii_1 + xi2 * (yi_1 + 3 * yi)))
a12 = np.sum(dxy * (yi_12 * (3 * xi_1 + xi) + 2 * yi * yi_1 * xii_1 + yi2 * (xi_1 + 3 * xi)))

if abs(a00) > flt_epsilon:
db1_2 = 0.5
db1_6 = 0.16666666666666666666666666666667
db1_12 = 0.083333333333333333333333333333333
db1_24 = 0.041666666666666666666666666666667
db1_20 = 0.05
db1_60 = 0.016666666666666666666666666666667

if a00 < 0:
db1_2 *= -1
db1_6 *= -1
db1_12 *= -1
db1_24 *= -1
db1_20 *= -1
db1_60 *= -1

m = dict(m00 = a00 * db1_2,
m10 = a10 * db1_6,
m01 = a01 * db1_6,
m20 = a20 * db1_12,
m11 = a11 * db1_24,
m02 = a02 * db1_12,
m30 = a30 * db1_20,
m21 = a21 * db1_60,
m12 = a12 * db1_60,
m03 = a03 * db1_20,
)

if m["m00"] > dbl_epsilon:
# Center of gravity
cx = m["m10"]/m["m00"]
cy = m["m01"]/m["m00"]
else:
cx = 0
cy = 0

# central second order moments
m["mu20"] = m["m20"] - m["m10"]*cx
m["mu11"] = m["m11"] - m["m10"]*cy
m["mu02"] = m["m02"] - m["m01"]*cy

m["mu30"] = m["m30"] - cx*(3*m["mu20"] + cx*m["m10"])
m["mu21"] = m["m21"] - cx*(2*m["mu11"] + cx*m["m01"]) - cy*m["mu20"]
m["mu12"] = m["m12"] - cy*(2*m["mu11"] + cy*m["m10"]) - cx*m["mu02"]
m["mu03"] = m["m03"] - cy*(3*m["mu02"] + cy*m["m01"])
return m
else:
return None