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 (<len of list1>, <len of list2>)
- """
- 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
|