dipy/dipy

class mam(Metric):
""" Computes the mam distance between two streamlines. """
def init(self):

    # For simplicity, features will be the vector between endpoints of a streamline.
    super(mam, self).__init__(feature=VectorOfEndpointsFeature())

def are_compatible(self, shape1, shape2):
    """ Checks if two features are vectors of same dimension.

    Basically this method exists so we don't have to do this check
    inside the `dist` method (speedup).
    """
    return shape1 == shape2 and shape1[0] == 1

def dist(self, v1, v2):

    track1 = np.ascontiguousarray(v1, dtype=np.float32)
    t1_len = track1.shape[0]
    track2 = np.ascontiguousarray(v2, dtype=np.float32)
    t2_len = track2.shape[0]
    # preallocate buffer array for track distance calculations

    #distances_buffer = np.zeros((t1_len ,), dtype=np.float32)
    min_t2t1 = np.zeros((t2_len ,), dtype=np.float32)
    min_t1t2 = np.zeros((t1_len ,), dtype=np.float32) 

    for t2_pi in range(0,t2_len):
        min_t2t1[t2_pi] = np.inf
    for t1_pi in range(0,t1_len):
        min_t1t2[t1_pi] = np.inf
    # pointer to current point in track 1
    t1_pt = track1
    t2_pt = track2
    # calculate min squared distance between each point in the two
    # lines.  Squared distance to delay doing the sqrt until after this
    # speed-critical loop
    for t1_pi in range(0,t1_len):
        # pointer to current point in track 2
        for t2_pi in range(0,t2_len):
            d0 = t1_pt[t1_pi][0] - t2_pt[t2_pi][0]
            d1 = t1_pt[t1_pi][1] - t2_pt[t2_pi][1]
            delta2 = d0*d0 + d1*d1 #+ d2*d2
            if delta2 < min_t1t2[t1_pi]:
                min_t1t2[t1_pi]=delta2

    for t2_pi in range(0,t2_len):
        # pointer to current point in track 2
        for t1_pi in range(0,t1_len):
            d0 = t1_pt[t1_pi][0] - t2_pt[t2_pi][0]
            d1 = t1_pt[t1_pi][1] - t2_pt[t2_pi][1]
            delta2 = d0*d0 + d1*d1 #+ d2*d2
            if delta2 < min_t2t1[t2_pi]:
                min_t2t1[t2_pi]=delta2

    # sqrt to get Euclidean distance from squared distance
    for t1_pi in range(0,t1_len):
        min_t1t2[t1_pi]=math.sqrt(min_t1t2[t1_pi])
    for t2_pi in range(0,t2_len):
        min_t2t1[t2_pi]=math.sqrt(min_t2t1[t2_pi])

    mean_t2t1 = 0
    mean_t1t2 = 0

    for t1_pi in range(0, t1_len):
        mean_t1t2+=min_t1t2[t1_pi]
    mean_t1t2=mean_t1t2 / t1_len
    for t2_pi in range(0, t2_len):
        mean_t2t1+=min_t2t1[t2_pi]
    mean_t2t1=mean_t2t1 / t2_len
    return np.min((mean_t2t1,mean_t1t2))

metric = mam()
qb2 = QuickBundles(threshold=0.15, metric=metric)
clus = qb2.cluster(streamlines)