def terminal_path_lengths_per_neurite(neurites, neurite_type=NeuriteType.all):
    """Get the path lengths to each terminal point per neurite in a collection"""
    return list(
        sectionfunc.section_path_length(s)
        for n in iter_neurites(neurites, filt=is_type(neurite_type))
        for s in iter_sections(n, iterator_type=Tree.ileaf)
    )

def trunk_origin_elevations(nrn, neurite_type=NeuriteType.all):
    '''Get a list of all the trunk origin elevations of a neuron or population

    The elevation is defined as the angle between x-axis and the
    vector defined by (initial tree point - soma center)
    on the x-y half-plane.

    The range of the elevation angle [-pi/2, pi/2] radians
    '''
    neurite_filter = is_type(neurite_type)
    nrns = neuron_population(nrn)

    def _elevation(section, soma):
        '''Elevation of a section'''
        vector = morphmath.vector(section[0], soma.center)
        norm_vector = np.linalg.norm(vector)

        if norm_vector >= np.finfo(type(norm_vector)).eps:
            return np.arcsin(vector[COLS.Y] / norm_vector)
        else:
            raise ValueError("Norm of vector between soma center and section is almost zero.")

    return [_elevation(s.root_node.points, n.soma)
            for n in nrns
            for s in n.neurites if neurite_filter(s)]

def partition_pairs(neurites, neurite_type=NeuriteType.all):
    '''Partition pairs at bifurcation points of a collection of neurites.
    Partition pait is defined as the number of bifurcations at the two
    daughters of the bifurcating section'''
    return map(_bifurcationfunc.partition_pair,
               iter_sections(neurites,
                             iterator_type=Tree.ibifurcation_point,
                             neurite_filter=is_type(neurite_type)))

def section_radial_distances(neurites, neurite_type=NeuriteType.all, origin=None):
    """Remote bifurcation angles in a collection of neurites"""
    dist = []
    for n in iter_neurites(neurites, filt=is_type(neurite_type)):
        pos = n.root_node.points[0] if origin is None else origin
        dist.extend([sectionfunc.section_radial_distance(s, pos) for s in n.iter_sections()])

    return dist

def segment_lengths(neurites, neurite_type=NeuriteType.all):
    """Lengths of the segments in a collection of neurites"""

    def _seg_len(sec):
        """list of segment lengths of a section"""
        return np.linalg.norm(np.diff(sec.points[:, : COLS.R], axis=0), axis=1)

    neurite_filter = is_type(neurite_type)
    return [s for ss in iter_sections(neurites, neurite_filter=neurite_filter) for s in _seg_len(ss)]

def map_segments(func, neurites, neurite_type):
    ''' Map `func` to all the segments in a collection of neurites

        `func` accepts a section and returns list of values corresponding to each segment.
    '''
    neurite_filter = is_type(neurite_type)
    return [
        s for ss in iter_sections(neurites, neurite_filter=neurite_filter) for s in func(ss)
    ]

def principal_direction_extents(neurites, neurite_type=NeuriteType.all, direction=0):
    '''Principal direction extent of neurites in neurons'''
    def _pde(neurite):
        '''Get the PDE of a single neurite'''
        # Get the X, Y,Z coordinates of the points in each section
        points = neurite.points[:, :3]
        return mm.principal_direction_extent(points)[direction]

    return map(_pde, iter_neurites(neurites, filt=is_type(neurite_type)))

def segment_lengths(neurites, neurite_type=NeuriteType.all):
    '''Lengths of the segments in a collection of neurites'''
    def _seg_len(sec):
        '''list of segment lengths of a section'''
        vecs = np.diff(sec.points, axis=0)[:, :3]
        return np.sqrt([np.dot(p, p) for p in vecs])

    neurite_filter = is_type(neurite_type)
    return [s for ss in iter_sections(neurites, neurite_filter=neurite_filter)
            for s in _seg_len(ss)]

def segment_radii(neurites, neurite_type=NeuriteType.all):
    """arithmetic mean of the radii of the points in segments in a collection of neurites"""

    def _seg_radii(sec):
        """vectorized mean radii"""
        pts = sec.points[:, COLS.R]
        return np.divide(np.add(pts[:-1], pts[1:]), 2.0)

    neurite_filter = is_type(neurite_type)
    return [s for ss in iter_sections(neurites, neurite_filter=neurite_filter) for s in _seg_radii(ss)]

def segment_midpoints(neurites, neurite_type=NeuriteType.all):
    """Return a list of segment mid-points in a collection of neurites"""

    def _seg_midpoint(sec):
        """Return the mid-points of segments in a section"""
        pts = sec.points
        return np.divide(np.add(pts[:-1], pts[1:])[:, :3], 2.0)

    neurite_filter = is_type(neurite_type)
    return [s for ss in iter_sections(neurites, neurite_filter=neurite_filter) for s in _seg_midpoint(ss)]

def section_radial_distances(neurites, neurite_type=NeuriteType.all, origin=None,
                             iterator_type=Tree.ipreorder):
    '''Section radial distances in a collection of neurites.
    The iterator_type can be used to select only terminal sections (ileaf)
    or only bifurcations (ibifurcation_point).'''
    dist = []
    for n in iter_neurites(neurites, filt=is_type(neurite_type)):
        pos = n.root_node.points[0] if origin is None else origin
        dist.extend(sectionfunc.section_radial_distance(s, pos)
                    for s in iter_sections(n,
                                           iterator_type=iterator_type))
    return dist

def neurite_volume_density(neurites, neurite_type=NeuriteType.all):
    """Get the volume density per neurite

    The volume density is defined as the ratio of the neurite volume and
    the volume of the neurite's enclosung convex hull
    """

    def vol_density(neurite):
        """volume density of a single neurite"""
        return neurite.volume / convex_hull(neurite).volume

    return list(vol_density(n) for n in iter_neurites(neurites, filt=is_type(neurite_type)))

def segment_radial_distances(neurites, neurite_type=NeuriteType.all, origin=None):
    '''Lengths of the segments in a collection of neurites'''
    def _seg_rd(sec, pos):
        '''list of radial distances of all segments of a section'''
        mid_pts = np.divide(np.add(sec.points[:-1], sec.points[1:])[:, :3], 2.0)
        return np.sqrt([mm.point_dist2(p, pos) for p in mid_pts])

    dist = []
    for n in iter_neurites(neurites, filt=is_type(neurite_type)):
        pos = n.root_node.points[0] if origin is None else origin
        dist.extend([s for ss in n.iter_sections() for s in _seg_rd(ss, pos)])

    return dist

def segment_taper_rates(neurites, neurite_type=NeuriteType.all):
    """taper rates of the segments in a collection of neurites

    The taper rate is defined as the absolute radii differences divided by length of the section
    """

    def _seg_taper_rates(sec):
        """vectorized taper rates"""
        pts = sec.points[:, : COLS.TYPE]
        diff = np.diff(pts, axis=0)
        distance = np.linalg.norm(diff[:, : COLS.R], axis=1)
        return np.divide(2 * np.abs(diff[:, COLS.R]), distance)

    neurite_filter = is_type(neurite_type)
    return [s for ss in iter_sections(neurites, neurite_filter=neurite_filter) for s in _seg_taper_rates(ss)]

def section_path_lengths(neurites, neurite_type=NeuriteType.all):
    """Path lengths of a collection of neurites """
    # Calculates and stores the section lengths in one pass,
    # then queries the lengths in the path length iterations.
    # This avoids repeatedly calculating the lengths of the
    # same sections.
    dist = {}
    neurite_filter = is_type(neurite_type)

    for s in iter_sections(neurites, neurite_filter=neurite_filter):
        dist[s] = s.length

    def pl2(node):
        """Calculate the path length using cached section lengths"""
        return sum(dist[n] for n in node.iupstream())

    return map_sections(pl2, neurites, neurite_type=neurite_type)

def trunk_origin_azimuths(nrn, neurite_type=NeuriteType.all):
    '''Get a list of all the trunk origin azimuths of a neuron or population

    The azimuth is defined as Angle between x-axis and the vector
    defined by (initial tree point - soma center) on the x-z plane.

    The range of the azimuth angle [-pi, pi] radians
    '''
    neurite_filter = is_type(neurite_type)
    nrns = nrn.neurons if hasattr(nrn, 'neurons') else [nrn]

    def _azimuth(section, soma):
        '''Azimuth of a section'''
        vector = mm.vector(section[0], soma.center)
        return np.arctan2(vector[COLS.Z], vector[COLS.X])

    return [_azimuth(s.root_node.points, n.soma)
            for n in nrns for s in n.neurites if neurite_filter(s)]

def n_neurites(neurites, neurite_type=NeuriteType.all):
    """Number of neurites in a collection of neurites"""
    return sum(1 for _ in iter_neurites(neurites, filt=is_type(neurite_type)))

def n_segments(neurites, neurite_type=NeuriteType.all):
    """Number of segments in a collection of neurites"""
    return sum(len(s.points) - 1 for s in iter_sections(neurites, neurite_filter=is_type(neurite_type)))

def trunk_origin_radii(nrn, neurite_type=NeuriteType.all):
    '''radii of the trunk sections of neurites in a neuron'''
    neurite_filter = is_type(neurite_type)
    return [s.root_node.points[0][COLS.R] for s in nrn.neurites if neurite_filter(s)]

def total_volume_per_neurite(neurites, neurite_type=NeuriteType.all):
    """Get the volume per neurite in a collection"""
    return list(sum(s.volume for s in n.iter_sections()) for n in iter_neurites(neurites, filt=is_type(neurite_type)))