1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495 ``````from scipy.spatial import ConvexHull from scipy.spatial.distance import cdist from .swcFuncs import readSWC_numpy import numpy as np from pyemd import emd def cdist_1d_centripetal(list1, list2, center): """ Return a matrix of pairwise displacements of entries of list1 from corresponding entries in list2, with displacements towards center having positive value and those away having negative values :param list1: iterable of floats or ints :param list2: iterable of floats or ints :param center: float or int :return: numpy.ndarray of shape (, ) """ assert all([isinstance(x, (float, int)) for x in list1]) assert all([isinstance(x, (float, int)) for x in list2]) assert isinstance(center, (float, int)) mesh2 = np.array([list2] * len(list1)) mesh1 = np.array([list1] * len(list2)).T mesh1CenteredAbs = np.abs(mesh1 - center) mesh2CenteredAbs = np.abs(mesh2 - center) return mesh2CenteredAbs - mesh1CenteredAbs def calcMaxDistances(swcList): """ Compute the convex hull of the union of points from all swcs in swcList. For each vertex of this convex hull, compute the distance of the farthest point among the vertices of the hull and return them. :param swcList: list of valid SWC files on the file system, list of strings. :return: list of maximum distances """ swcPointSets = [] for swc in swcList: headr, swcData = readSWC_numpy(swc) swcPointSets.append(swcData[:, 2:5]) unionWithDuplicates = np.concatenate(swcPointSets, axis=0) if any(np.abs(unionWithDuplicates).max(axis=0) == 0): raise(ValueError("The list of SWCs all lie on a plane or on a line and hence do not " "for a 3D point cloud. Such SWCs are not supported.")) hull = ConvexHull(unionWithDuplicates) vertices = unionWithDuplicates[hull.vertices, :] distMatrix = cdist(unionWithDuplicates, vertices) maxDistances = distMatrix.max(axis=1).tolist() return maxDistances def maxDistEMD(swcList): """ Calculate the maxDistance based metric. It is the size normalized Earth mover distance between the distribution of maxDistances (see calcMaxDistances above) of the pooled collection of points of all swcs in swcList and the distribution of pooled diagonal maxDistances of individual swcs in swcList :param swcList: list of valid swc files on the system, list of strings :return: float """ individualMaxDistances = [calcMaxDistances([swc]) for swc in swcList] pooledIndividualMaxDistances = np.concatenate(individualMaxDistances, axis=0) meanPIMD = pooledIndividualMaxDistances.mean() PIMDNorm = (pooledIndividualMaxDistances - meanPIMD) / meanPIMD maxDistancesAllPts = np.array(calcMaxDistances(swcList)) MDAPNorm = (maxDistancesAllPts - meanPIMD) / meanPIMD binWidth = 1 / meanPIMD bins = np.arange(MDAPNorm.min() - 0.5 * binWidth, MDAPNorm.max() + 0.5 * binWidth, binWidth) hist1, bins1 = np.histogram(MDAPNorm, bins) hist2, bins2 = np.histogram(PIMDNorm, bins1) hist1Normed = hist1 / float(hist1.sum()) hist2Normed = hist2 / float(hist2.sum()) dist_metric = cdist_1d_centripetal(bins1, bins2, center=0) emd_val = emd(np.asarray(hist2Normed, dtype=np.float64), np.asarray(hist1Normed, dtype=np.float64), np.asarray(dist_metric, dtype=np.float64)) return emd_val ``````