def triangle_area(a,b,c):
    '''
    triangle_area(a, b, c) yields the area of the triangle whose vertices are given by the points a,
    b, and c.
    '''
    (a,b,c) = [np.asarray(x) for x in (a,b,c)]
    sides = np.sqrt(np.sum([(p1 - p2)**2 for (p1,p2) in zip([b,c,a],[c,a,b])], axis=1))
    s = 0.5 * np.sum(sides, axis=0)
    s = np.clip(s * np.prod(s - sides, axis=0), 0.0, None)
    return np.sqrt(s) def triangle_address(fx, pt):
    '''
    triangle_address(FX, P) yields an address coordinate (t,r) for the point P in the triangle
    defined by the (3 x d)-sized coordinate matrix FX, in which each row of the matrix is the
    d-dimensional vector representing the respective triangle vertx for triangle [A,B,C]. The
    resulting coordinates (t,r) (0 <= t <= 1, 0 <= r <= 1) address the point P such that, if t gives
    the fraction of the angle from vector AB to vector AC that is made by the angle between vectors
    AB and AP, and r gives the fraction ||AP||/||AR|| where R is the point of intersection between
    lines AP and BC. If P is a (d x n)-sized matrix of points, then a (2 x n) matrix of addresses
    is returned.
    '''
    fx = np.asarray(fx)
    pt = np.asarray(pt)
    # The triangle vectors...
    ab = fx[1] - fx[0]
    ac = fx[2] - fx[0]
    bc = fx[2] - fx[1]
    ap = np.asarray([pt_i - a_i for (pt_i, a_i) in zip(pt, fx[0])])
    # get the unnormalized distance...
    r = np.sqrt((ap ** 2).sum(0))
    # now we can find the angle...
    unit = 1 - r.astype(bool)
    t0 = vector_angle(ab, ac)
    t = vector_angle(ap + [ab_i * unit for ab_i in ab], ab)
    sint = np.sin(t)
    sindt = np.sin(t0 - t)
    # finding r0 is tricker--we use this fancy formula based on the law of sines
    q0 = np.sqrt((bc ** 2).sum(0))          # B->C distance
    beta = vector_angle(-ab, bc)            # Angle at B
    sinGamma = np.sin(math.pi - beta - t0)
    sinBeta  = np.sin(beta)
    r0 = q0 * sinBeta * sinGamma / (sinBeta * sindt + sinGamma * sint)
    return np.asarray([t/t0, r/r0]) def alignment_matrix_3D(u, v):
    '''
    alignment_matrix_3D(u, v) yields a 3x3 numpy array that rotates the vector u to the vector v
    around the origin.
    '''
    # normalize both vectors:
    u = normalize(u)
    v = normalize(v)
    # get the cross product of u cross v
    w = np.cross(u, v)
    # the angle we need to rotate
    th = vector_angle(u, v)
    # we rotate around this vector by the angle between them
    return rotation_matrix_3D(w, th) def rotation_matrix_3D(u, th):
    """
    rotation_matrix_3D(u, t) yields a 3D numpy matrix that rotates any vector about the axis u
    t radians counter-clockwise.
    """
    # normalize the axis:
    u = normalize(u)
    # We use the Euler-Rodrigues formula;
    # see https://en.wikipedia.org/wiki/Euler-Rodrigues_formula
    a = math.cos(0.5 * th)
    s = math.sin(0.5 * th)
    (b, c, d) = -s * u
    (a2, b2, c2, d2) = (a*a, b*b, c*c, d*d)
    (bc, ad, ac, ab, bd, cd) = (b*c, a*d, a*c, a*b, b*d, c*d)
    return np.array([[a2 + b2 - c2 - d2, 2*(bc + ad),         2*(bd - ac)],
                     [2*(bc - ad),       a2 + c2 - b2 - d2,   2*(cd + ab)],
                     [2*(bd + ac),       2*(cd - ab),         a2 + d2 - b2 - c2]]) def cartesian_to_barycentric_3D(tri, xy):
    '''
    cartesian_to_barycentric_3D(tri,xy) is identical to cartesian_to_barycentric_2D(tri,xy) except
    it works on 3D data. Note that if tri is a 3 x 3 x n, a 3 x n x 3 or an n x 3 x 3 matrix, the
    first dimension must always be the triangle vertices and the second 3-sized dimension must be
    the (x,y,z) coordinates.
    '''
    xy = np.asarray(xy)
    tri = np.asarray(tri)
    if len(xy.shape) == 1:
        return cartesian_to_barycentric_3D(np.transpose(np.asarray([tri]), (1,2,0)),
                                           np.asarray([xy]).T)[:,0]
    xy = xy if xy.shape[0] == 3 else xy.T
    if tri.shape[0] == 3:
        tri = tri if tri.shape[1] == 3 else np.transpose(tri, (0,2,1))
    elif tri.shape[1] == 3:
        tri = tri.T if tri.shape[0] == 3 else np.transpose(tri, (1,2,0))
    elif tri.shape[2] == 3:
        tri = np.transpose(tri, (2,1,0) if tri.shape[1] == 3 else (2,0,1))
    if tri.shape[0] != 3 or tri.shape[1] != 3:
        raise ValueError('Triangle array did not have dimensions of sizes 3 and 3')
    if xy.shape[0] != 3:
        raise ValueError('coordinate matrix did not have a dimension of size 3')
    if tri.shape[2] != xy.shape[1]:
        raise ValueError('number of triangles and coordinates must match')
    # The algorithm here is borrowed from this stack-exchange post:
    # http://gamedev.stackexchange.com/questions/23743
    # in which it is attributed to Christer Ericson's book Real-Time Collision Detection.
    v0 = tri[1] - tri[0]
    v1 = tri[2] - tri[0]
    v2 = xy - tri[0]
    d00 = np.sum(v0 * v0, axis=0)
    d01 = np.sum(v0 * v1, axis=0)
    d11 = np.sum(v1 * v1, axis=0)
    d20 = np.sum(v2 * v0, axis=0)
    d21 = np.sum(v2 * v1, axis=0)
    den = d00*d11 - d01*d01
    zero = np.isclose(den, 0)
    unit = 1 - zero
    den += zero
    l2 = unit * (d11 * d20 - d01 * d21) / den
    l3 = unit * (d00 * d21 - d01 * d20) / den
    return np.asarray([1.0 - l2 - l3, l2]) def cartesian_to_barycentric_2D(tri, xy, atol=1e-8):
    '''
    cartesian_to_barycentric_2D(tri, xy) yields a (2 x n) barycentric coordinate matrix 
    (or just a tuple if xy is a single (x, y) coordinate) of the first two barycentric coordinates
    for the triangle coordinate array tri. The array tri should be 3 (vertices) x 2 (coordinates) x
    n (triangles) unless xy is a tuple, in which case it should be a (3 x 2) matrix.
    '''
    xy = np.asarray(xy)
    tri = np.asarray(tri)
    if len(xy.shape) == 1:
        return cartesian_to_barycentric_2D(np.transpose(np.asarray([tri]), (1,2,0)),
                                           np.asarray([xy]).T)[:,0]
    xy = xy if xy.shape[0] == 2 else xy.T
    if tri.shape[0] == 3:
        tri = tri if tri.shape[1] == 2 else np.transpose(tri, (0,2,1))
    elif tri.shape[1] == 3:
        tri = tri.T if tri.shape[0] == 2 else np.transpose(tri, (1,2,0))
    elif tri.shape[2] == 3:
        tri = np.transpose(tri, (2,1,0) if tri.shape[1] == 2 else (2,0,1))
    if tri.shape[0] != 3 or tri.shape[1] != 2:
        raise ValueError('Triangle array did not have dimensions of sizes 3 and 2')
    if xy.shape[0] != 2:
        raise ValueError('coordinate matrix did not have a dimension of size 2')
    if tri.shape[2] != xy.shape[1]:
        raise ValueError('number of triangles and coordinates must match')
    # Okay, everything's the right shape...
    (x,y) = xy
    ((x1,y1), (x2,y2), (x3,y3)) = tri
    x_x3  = x  - x3
    x1_x3 = x1 - x3
    x3_x2 = x3 - x2
    y_y3  = y  - y3
    y1_y3 = y1 - y3
    y2_y3 = y2 - y3
    num1 = (y2_y3*x_x3  + x3_x2*y_y3)
    num2 = (-y1_y3*x_x3 + x1_x3*y_y3)
    den  = (y2_y3*x1_x3 + x3_x2*y1_y3)
    zero = np.isclose(den, 0, atol=atol)
    den += zero
    unit = 1 - zero
    l1 = unit * num1 / den
    l2 = unit * num2 / den
    return np.asarray((l1, l2)) def vector_angle_cos(u, v):
    '''
    vector_angle_cos(u, v) yields the cosine of the angle between the two vectors u and v. If u
    or v (or both) is a (d x n) matrix of n vectors, the result will be a length n vector of the
    cosines.
    '''
    u = np.asarray(u)
    v = np.asarray(v)
    return (u * v).sum(0) / np.sqrt((u ** 2).sum(0) * (v ** 2).sum(0)) def segment_intersection_2D(p12arg, p34arg, atol=1e-8):
    '''
    segment_intersection((a, b), (c, d)) yields the intersection point between the line segments
    that pass from point a to point b and from point c to point d. If there is no intersection
    point, then (numpy.nan, numpy.nan) is returned.
    '''
    (p1,p2) = p12arg
    (p3,p4) = p34arg
    pi = np.asarray(line_intersection_2D(p12arg, p34arg, atol=atol))
    p1 = np.asarray(p1)
    p2 = np.asarray(p2)
    p3 = np.asarray(p3)
    p4 = np.asarray(p4)
    u12 = p2 - p1
    u34 = p4 - p3
    cfn = lambda px,iis: (px if iis is None or len(px.shape) == 1 or px.shape[1] == len(iis) else
                          px[:,iis])
    dfn = lambda a,b:     a[0]*b[0] + a[1]*b[1]
    sfn = lambda a,b:     ((a-b)                 if len(a.shape) == len(b.shape) else
                           (np.transpose([a])-b) if len(a.shape) <  len(b.shape) else
                           (a - np.transpose([b])))
    fn  = lambda px,iis:  (1 - ((dfn(cfn(u12,iis), sfn(         px, cfn(p1,iis))) > 0) *
                                (dfn(cfn(u34,iis), sfn(         px, cfn(p3,iis))) > 0) *
                                (dfn(cfn(u12,iis), sfn(cfn(p2,iis),          px)) > 0) *
                                (dfn(cfn(u34,iis), sfn(cfn(p4,iis),          px)) > 0)))
    if len(pi.shape) == 1:
        if not np.isfinite(pi[0]): return (np.nan, np.nan)
        bad = fn(pi, None)
        return (np.nan, np.nan) if bad else pi
    else:
        nonpar = np.where(np.isfinite(pi[0]))[0]
        bad = fn(cfn(pi, nonpar), nonpar)
        (xi,yi) = pi
        bad = nonpar[np.where(bad)[0]]
        xi[bad] = np.nan
        yi[bad] = np.nan
        return (xi,yi) def segments_overlapping(ab, cd, atol=1e-8):
    '''
    segments_overlapping((a, b), (c, d)) yields True if the line segments (a,b) and (c,d) are both
    colinear and have a non-finite overlap. If (a,b) and (c,d) touch at a single point, they are not
    considered overlapping.
    '''
    (a,b) = ab
    (c,d) = cd
    ss = [point_in_segment(ab, c, atol=atol), point_in_segment(ab, d, atol=atol),
          point_in_segment(cd, a, atol=atol), point_in_segment(cd, b, atol=atol)]
    return (~(points_close(a,b) | points_close(c,d, atol=atol)) & 
            ((np.sum(ss, axis=0) > 1) |
             (points_close(a,c, atol=atol) & points_close(b,d, atol=atol)) |
             (points_close(a,d, atol=atol) & points_close(b,c, atol=atol)))) def points_close(a,b, atol=1e-8):
    '''
    points_close(a,b) yields True if points a and b are close to each other and False otherwise.
    '''
    (a,b) = [np.asarray(u) for u in (a,b)]
    if len(a.shape) == 2 or len(b.shape) == 2: (a,b) = [np.reshape(u,(len(u),-1)) for u in (a,b)]
    return np.isclose(np.sqrt(np.sum((a - b)**2, axis=0)), 0, atol=atol) def point_in_segment(ac, b, atol=1e-8):
    '''
    point_in_segment((a,b), c) yields True if point x is in segment (a,b) and False otherwise. Note
    that this differs from point_on_segment in that a point that if c is equal to a or b it is
    considered 'on' but not 'in' the segment.
    The option atol can be given and is used only to test for difference from 0; by default it is
    1e-8.
    '''
    (a,c) = ac
    abc = [np.asarray(u) for u in (a,b,c)]
    if any(len(u.shape) > 1 for u in abc): (a,b,c) = [np.reshape(u,(len(u),-1)) for u in abc]
    else:                                  (a,b,c) = abc
    vab = b - a
    vbc = c - b
    vac = c - a
    dab = np.sqrt(np.sum(vab**2, axis=0))
    dbc = np.sqrt(np.sum(vbc**2, axis=0))
    dac = np.sqrt(np.sum(vac**2, axis=0))
    return (np.isclose(dab + dbc - dac, 0, atol=atol) &
            ~np.isclose(dac - dab, 0, atol=atol) &
            ~np.isclose(dac - dbc, 0, atol=atol)) def point_on_segment(ac, b, atol=1e-8):
    '''
    point_on_segment((a,b), c) yields True if point x is on segment (a,b) and False otherwise. Note
    that this differs from point_in_segment in that a point that if c is equal to a or b it is
    considered 'on' but not 'in' the segment.
    The option atol can be given and is used only to test for difference from 0; by default it is
    1e-8.
    '''
    (a,c) = ac
    abc = [np.asarray(u) for u in (a,b,c)]
    if any(len(u.shape) > 1 for u in abc): (a,b,c) = [np.reshape(u,(len(u),-1)) for u in abc]
    else:                                  (a,b,c) = abc
    vab = b - a
    vbc = c - b
    vac = c - a
    dab = np.sqrt(np.sum(vab**2, axis=0))
    dbc = np.sqrt(np.sum(vbc**2, axis=0))
    dac = np.sqrt(np.sum(vac**2, axis=0))
    return np.isclose(dab + dbc - dac, 0, atol=atol) def barycentric_to_cartesian(tri, bc):
    '''
    barycentric_to_cartesian(tri, bc) yields the d x n coordinate matrix of the given barycentric
    coordinate matrix (also d x n) bc interpolated in the n triangles given in the array tri. See
    also cartesian_to_barycentric. If tri and bc represent one triangle and coordinate, then just
    the coordinate and not a matrix is returned. The value d, dimensions, must be 2 or 3.
    '''
    bc = np.asarray(bc)
    tri = np.asarray(tri)
    if len(bc.shape) == 1:
        return barycentric_to_cartesian(np.transpose(np.asarray([tri]), (1,2,0)),
                                        np.asarray([bc]).T)[:,0]
    bc = bc if bc.shape[0] == 2 else bc.T
    if bc.shape[0] != 2: raise ValueError('barycentric matrix did not have a dimension of size 2')
    n = bc.shape[1]
    # we know how many bc's there are now; lets reorient tri to match with the last dimension as n
    if len(tri.shape) == 2:
        tri = np.transpose([tri for _ in range(n)], (1,2,0))
    # the possible orientations of tri:
    if tri.shape[0] == 3:
        if tri.shape[1] in [2,3] and tri.shape[2] == n:
            pass # default orientation
        elif tri.shape[1] == n and tri.shape[2] in [2,3]:
            tri = np.transpose(tri, (0,2,1))
        else: raise ValueError('could not deduce triangle dimensions')
    elif tri.shape[1] == 3:
        if tri.shape[0] in [2,3] and tri.shape[2] == n:
            tri = np.transpose(tri, (1,0,2))
        elif tri.shape[0] == n and tri.shape[2] in [2,3]:
            tri = np.transpose(tri, (1,2,0))
        else: raise ValueError('could not deduce triangle dimensions')
    elif tri.shape[2] == 3:
        if tri.shape[0] in [2,3] and tri.shape[1] == n:
            tri = np.transpose(tri, (2,0,1))
        elif tri.shape[0] == n and tri.shape[1] in [2,3]:
            tri = np.transpose(tri, (2,1,0))
        else: raise ValueError('could not deduce triangle dimensions')
    else: raise ValueError('At least one dimension of triangles must be 3')
    if tri.shape[0] != 3 or (tri.shape[1] not in [2,3]):
        raise ValueError('Triangle array did not have dimensions of sizes 3 and (2 or 3)')
    if tri.shape[2] != n:
        raise ValueError('number of triangles and coordinates must match')
    (l1,l2) = bc
    (p1, p2, p3) = tri
    l3 = (1 - l1 - l2)
    return np.asarray([x1*l1 + x2*l2 + x3*l3 for (x1,x2,x3) in zip(p1, p2, p3)]) def point_in_triangle(tri, pt, atol=1e-13):
    tri = np.asarray(tri)
    pt  = np.asarray(pt)
    if len(tri.shape) == 2 and len(pt.shape) == 1:
        if len(pt) == 2:
            tol = atol
            v0 = tri[2] - tri[0]
            v1 = tri[1] - tri[0]
            v2 = pt - tri[0]
            d00 = np.dot(v0, v0)
            d01 = np.dot(v0, v1)
            d02 = np.dot(v0, v2)
            d11 = np.dot(v1, v1)
            d12 = np.dot(v1, v2)
            invDenom = (d00*d11 - d01*d01)
            if np.isclose(invDenom, 0, atol=atol): return False
            s = (d11*d02 - d01*d12) / invDenom
            if (s + tol) < 0 or (s - tol) >= 1: return False
            t = (d00*d12 - d01*d02) / invDenom
            #return False if (t + tol) < 0 or (s + t - tol) > 1 else True
            s = s + t
            tp = np.isclose(t, 0, atol=tol) or t >= 0
            tn = np.isclose(s - 1, 0, atol=tol) or s <= 1
            return tn and tp
        else:
            dp1 = np.dot(pt - tri[0], np.cross(tri[0], tri[1] - tri[0]))
            dp2 = np.dot(pt - tri[1], np.cross(tri[1], tri[2] - tri[1]))
            db3 = np.dot(pt - tri[2], np.cross(tri[2], tri[0] - tri[2]))
            return ((dp1 > 0 or np.isclose(dp1, 0, atol=atol)) and
                    (dp2 > 0 or np.isclose(dp2, 0, atol=atol)) and
                    (dp3 > 0 or np.isclose(dp3, 0, atol=atol)))
    elif len(tri.shape) == 3 and len(pt.shape) == 2:
        if len(pt) != len(tri):
            raise ValueError('the number of triangles and points must be equal')
        if pt.shape[1] == 2:
            tol = atol
            v0 = tri[:,2] - tri[:,0]
            v1 = tri[:,1] - tri[:,0]
            v2 = pt - tri[:,0]
            d00 = np.sum(v0*v0, axis=1)
            d01 = np.sum(v0*v1, axis=1)
            d02 = np.sum(v0*v2, axis=1)
            d11 = np.sum(v1*v1, axis=1)
            d12 = np.sum(v1*v2, axis=1)
            invDenom = (d00*d11 - d01*d01)
            zeros = np.isclose(invDenom, 0)
            invDenom[zeros] = 1.0
            s = (d11*d02 - d01*d12) / invDenom
            t = (d00*d12 - d01*d02) / invDenom
            #return ~((s + tol < 0) | (s - tol >= 1) | (t + tol < 0) | (s + t - tol > 1) | zeros)
            #return (~(s + tol < 0) & ~(s - tol >= 1) & ~(t + tol < 0) & ~(s + t - tol > 1) & ~zeros)
            #return ((s + tol >= 0) & (s - tol < 1) & (t + tol >= 0) & (s + t - tol <= 1) & ~zeros)
            st = s + t
            q1 = np.isclose(s, 0, atol=atol) | (s > 0)
            q2 = ~np.isclose(s - 1, 0, atol=atol) & (s < 1)
            q3 = np.isclose(t, 0, atol=atol) | (t > 0)
            q4 = np.isclose(st - 1, 0, atol=atol) | (st < 1)
            return q1 * q2 & q3 & q4 & ~zeros
        else:
            x0 = np.sum((pt - tri[:,0]) * np.cross(tri[:,0], tri[:,1] - tri[:,0], axis=1), axis=1)
            x1 = np.sum((pt - tri[:,1]) * np.cross(tri[:,1], tri[:,2] - tri[:,1], axis=1), axis=1)
            x2 = np.sum((pt - tri[:,2]) * np.cross(tri[:,2], tri[:,0] - tri[:,2], axis=1), axis=1)
            return (((x0 > 0) | np.isclose(x0, 0, atol=atol)) &
                    ((x1 > 0) | np.isclose(x1, 0, atol=atol)) &
                    ((x2 > 0) | np.isclose(x2, 0, atol=atol)))
    elif len(tri.shape) == 3 and len(pt.shape) == 1:
        return point_in_triangle(tri, np.asarray([pt for _ in tri]))
    elif len(tri.shape) == 2 and len(pt.shape) == 1:
        return point_in_triangle(np.asarray([tri for _ in pt]), pt)
    else:
        raise ValueError('triangles and pts do not have parallel shapes') def line_segment_intersection_2D(p12arg, p34arg, atol=1e-8):
    '''
    line_segment_intersection((a, b), (c, d)) yields the intersection point between the line
    passing through points a and b and the line segment that passes from point c to point d. If
    there is no intersection point, then (numpy.nan, numpy.nan) is returned.
    '''
    (p1,p2) = p12arg
    (p3,p4) = p34arg
    pi = np.asarray(line_intersection_2D(p12arg, p34arg, atol=atol))
    p3 = np.asarray(p3)
    u34 = p4 - p3
    cfn = lambda px,iis: (px if iis is None or len(px.shape) == 1 or px.shape[1] == len(iis) else
                          px[:,iis])
    dfn = lambda a,b:     a[0]*b[0] + a[1]*b[1]
    sfn = lambda a,b:     ((a-b)                 if len(a.shape) == len(b.shape) else
                           (np.transpose([a])-b) if len(a.shape) <  len(b.shape) else
                           (a - np.transpose([b])))
    fn  = lambda px,iis:  (1 - ((dfn(cfn(u34,iis), sfn(         px, cfn(p3,iis))) > 0) *
                                (dfn(cfn(u34,iis), sfn(cfn(p4,iis),          px)) > 0)))
    if len(pi.shape) == 1:
        if not np.isfinite(pi[0]): return (np.nan, np.nan)
        bad = fn(pi, None)
        return (np.nan, np.nan) if bad else pi
    else:
        nonpar = np.where(np.isfinite(pi[0]))[0]
        bad = fn(cfn(pi, nonpar), nonpar)
        (xi,yi) = pi
        bad = nonpar[np.where(bad)[0]]
        xi[bad] = np.nan
        yi[bad] = np.nan
        return (xi,yi) def segments_touch_2D(ab, cd, atol=1e-8):
    '''
    segmentss_touch_2D((a,b), (c,d)) is equivalent to segments_colinear((a,b), (c,d)) |
    numpy.isfinite(segment_intersection_2D((a,b), (c,d))[0])
    '''
    return (segments_colinear(ab, cd, atol=atol) |
            np.isfinite(segment_intersection_2D(ab, cd, atol=atol)[0])) def to_hemi_str(s):
    def __init__(self, meta_data=None):
    def meta_data(md):
    def meta(self, name, missing=None):
    def with_meta(self, *args, **kwargs):
    def wout_meta(self, *args, **kwargs):
    def normalize(self):
    def denormalize(self, params):
def normalize(data):
def denormalize(data):
def to_affine(aff, dims=None):
def apply_affine(aff, coords):
def is_dataframe(d):
def to_dataframe(d, **kw):
def dataframe_select(df, *cols, **filters):
def dataframe_except(df, *cols, **filters):
    def __init__(self, *args, **kwargs):
    def __missing__(self, key):
def auto_dict(ival=None, miss=None):
def simplex_summation_matrix(simplices, weight=None, inverse=False):
def simplex_averaging_matrix(simplices, weight=None, inverse=False):
def is_image(image):
def is_image_header(x):
def is_address(data):
def address_data(data, dims=None, surface=0.5, strict=True):
def address_interpolate(addr, prop, method=None, surface='midgray', strict=False, null=np.nan,
                        tess=None):
def numel(x):
def rows(x):
def check_sparsity(x, fraction=0.6):
def unbroadcast(a, b):
        def fix(sp,nm):
def cplus(*args):
def plus(*args):
def cminus(a, b):
def minus(a, b):
def ctimes(*args):
def times(*args):
def inv(x):
def zinv(x, null=0):
def cdivide(a, b):
def divide(a, b):
def czdivide(a, b, null=0):
def zdivide(a, b, null=0):
def cpower(a,b):
def power(a,b):
def inner(a,b):
def sine(x):
def cosine(x):
def tangent(x, null=(-np.inf, np.inf), rtol=default_rtol, atol=default_atol):
def cotangent(x, null=(-np.inf, np.inf), rtol=default_rtol, atol=default_atol):
def secant(x, null=(-np.inf, np.inf), rtol=default_rtol, atol=default_atol):
def cosecant(x, null=(-np.inf, np.inf), rtol=default_rtol, atol=default_atol):
def arcsine(x, null=(-np.inf, np.inf)):
def arccosine(x, null=(-np.inf, np.inf)):
def arctangent(y, x=None, null=0, broadcast=False, rtol=default_rtol, atol=default_atol):
def flattest(x):
def flatter(x, k=1):
def part(x, *args):
def hstack(tup):
def vstack(tup):
def repmat(x, r, c):
def replace_close(x, xhat, rtol=default_rtol, atol=default_atol, copy=True):
def chop(x, rtol=default_rtol, atol=default_atol, copy=True):
def nan_compare(f, x, y, nan_nan=False, nan_val=False, val_nan=False):
def naneq(x, y, nan_nan=False, nan_val=False, val_nan=False):
def nanne(x, y, nan_nan=False, nan_val=False, val_nan=False):
def nanlt(x, y, nan_nan=False, nan_val=False, val_nan=False):
def nanle(x, y, nan_nan=False, nan_val=False, val_nan=False):
def nangt(x, y, nan_nan=False, nan_val=False, val_nan=False):
def nange(x, y, nan_nan=False, nan_val=False, val_nan=False):
def nanlog(x, null=np.nan):
def library_path():
    def __init__(self, x, y=None,
                 order=1, weights=None, smoothing=None, periodic=False,
                 distances=None,
                 meta_data=None):
    def coordinates(x):
    def order(o):
    def smoothing(s):
    def weights(w):
    def periodic(p):
    def distances(ds):
    def check_distances(distances, coordinates, periodic):
    def t(distances,coordinates):
    def splrep(coordinates, t, order, weights, smoothing, periodic):
    def __repr__(self):
    def __call__(self, t, derivative=0):
    def curve_length(self, start=None, end=None, precision=0.01):
    def linspace(self, n=100, derivative=0):
    def even_out(self, precision=0.001):
    def reverse(self):
    def subcurve(self, t0, t1):
def curve_spline(x, y=None, weights=None, order=1, even_out=True,
                 smoothing=None, periodic=False, meta_data=None):
def is_curve_spline(obj):
def to_curve_spline(obj):
def curve_intersection(c1, c2, grid=16):
            def f(t): return np.sum((c1(t[0]) - c2(t[1]))**2)
    def f(t): return np.sum((c1(t[0]) - c2(t[1]))**2)
def close_curves(*crvs, **kw):
    def __init__(self, **kw):    self.__dict__.update(kw)
    def __setattr__(self, k, v): raise ValueError('DataStruct objects are immutable')
    def __delattr__(self, k):    raise ValueError('DataStruct objects are immutable')
    def set(**kw):
    def delete(*args):
def data_struct(*args, **kw):
def tmpdir(prefix='npythy_tempdir_', delete=True):
def dirpath_to_list(p):
def try_until(*args, **kw):
class ObjectWithMetaData(object):
class AutoDict(dict):
class CurveSpline(ObjectWithMetaData):
class DataStruct(object): class config(object):
    '''
    neuropythy.util.conf.config is a class that manages configuration items for neuropythy and their
    interface. This class reads in the user's neuropythy-rc file at startup, which by default is in
    the user's home directory named "~/.npythyrc" though it may be altered by setting the
    environment variable NPYTHYRC. Environment variables that are associated with configurable
    variables always override the values in the RC file, and any direct set

    To declare a configurable variable in Neuropythy, you must use config.declare().
    '''
    _rc = None
    @staticmethod
    def rc():
        '''
        config.rc() yields the data imported from the Neuropythy rc file, if any.
        '''
        if config._rc is None:
            # First: We check to see if we have been given a custom nptyhrc file:
            npythyrc_path = os.path.expanduser('~/.npythyrc')
            if 'NPYTHYRC' in os.environ:
                npythyrc_path = os.path.expanduser(os.path.expandvars(os.environ['NPYTHYRC']))
            # the default config:
            if os.path.isfile(npythyrc_path):
                try:
                    config._rc = loadrc(npythyrc_path)
                    config._rc['npythyrc_loaded'] = True
                except Exception as err:
                    warnings.warn('Could not load neuropythy RC file: %s' % npythyrc_path)
                    config._rc = {'npythyrc_loaded':False,
                                  'npythyrc_error': err}
            else:
                config._rc = {'npythyrc_loaded':False}
            config._rc['npythyrc'] = npythyrc_path
        return config._rc
    _vars = {}
    @staticmethod
    def declare(name, rc_name=None, environ_name=None, filter=None, merge=None, default_value=None):
        '''
        config.declare(name) registers a configurable variable with the given name to the neuropythy
          configuration system. This allows the variable to be looked up in the neuropythy RC-file
          and the neuropythy environment variables.

        By default, the variable will be assumed to have the identical name in the neuropythy
        RC-file and the environment variable ('NPYTHY_' + name.upper()) is searched for in the
        environment. The environment variable will always overwrite the RC-file value if both are
        provided. Note that all inputs from the environment or the RC-file are parsed as JSON inputs
        prior to being put in the config object itself.

        The following optional arguments may be provided:
          * rc_name (default: None) specifies the name of the variable in the neuropythy RC-file; if
            None, then uses the name.lower().
          * environ_name (default: None) specifies the name of the variable to look for in the
            os.environ dictionary; if None, then uses ('NPYTHY_' + name.upper()).
          * filter (default: None) specifies a function f that is passed the provided value u; the 
            resulting config value is f(u).
          * merge (default: None) specifies how the environment variable should be merged with the
            value found in the rc file if both are defined. If this value is None or False, then the
            environment variable value always overrides the rc-file value; otherwise, this must be
            a function f(rc_value, environ_value) and it must return the value that should be used.
          * default_value (default: None) specifies the default value the configuration item should
            take if not provided.
        '''
        if rc_name is None: rc_name = name.lower()
        if environ_name is None: environ_name = 'NPYTHY_' + name.upper()
        # okay, see if the name exists
        if name in config._vars: raise ValueError('Multiple config items declared for %s' % name)
        if merge is False: merge = None
        config._vars[name] = (rc_name, environ_name, filter, default_value, merge)
        return True
    @staticmethod
    def declare_credentials(name, rc_name=None, environ_name=None, default_value=None, merge=None,
                            extra_environ=None, filenames=None, aws_profile_name=None):
        '''
        config.declare_credentials(name) registers a configurable variable with the given name to
          the neuropythy configuration system (see also config.declare()). The variable is
          specifically created as a credentials-type variable such that if the config variable is 
          not set, a call to detect_credentials is used to attempt to find the credentials; this
          effectively allows for additional filenames and environment variables to be searched while
          ensuring that strings, files, and json representations all get coerced to (key, secret)
          if possible.

        Note that if you set default_value to something that is neither None nor a valid credentials
        representation, then the config item will still end up as None; this is because the
        default value must fall through the to_credentials() function.

        The following optional arguments from detect_credentials are accepted and passed along to it
        if necessary:
          * extra_environ
          * filenames
          * aws_profile_name
          * default_value (Ellipsis is not allowed via this interface).
        The following optional arguments from config.declare() are accepted and passed along to it:
          * rc_name
          * environ_name
          * merge

        See also help(config.declare) and help(detect_credentials).
        '''
        if default_value is Ellipsis:
            raise ValueError('Ellipsis is not a valid credentials default_value when calling the'
                             ' declare_credentials() method.')
        def _check_creds(creds):
            if creds is None:
                return detect_credentials(None,
                                          extra_environ=extra_environ,
                                          filenames=filenames,
                                          aws_profile_name=aws_profile_name,
                                          default=default_value)
            else: return to_credentials(creds)
        return config.declare(name,
                              rc_name=rc_name,
                              environ_name=environ_name,
                              filter=_check_creds,
                              merge=merge,
                              default_value=None)
    @staticmethod
    def declare_path(name, 
                     rc_name=None, environ_name=None, filter=None, default_value=None, merge=None,
                     fail_value=None):
        '''
        config.declare_path(...) is equivalent to config.declare(...) except that it requires in
          addition that the value found for the config item be an existing file or dir on the local
          filesystem. Only if this is true is the value passed to any provided filter.

        This function (similar to config.declare_file and config.declare_dir) is generally expected
        to return some data--the contents of the path--not the path-name itself. For this reason
        the handling of the default_value argument is expanded relative to that of config.declare().
        If no value is provided or if all provided values are either not paths or fail the filter,
        then the default_value is checked (as a path); if this also fails then the fail_value is
        used as the config value. By default this is None, but note that by using a filter function
        that loads the provided path (which will be guaranteed to exist) and yields the data along
        with a fail_value that resembles 'default' data, one can effectively guarantee that data is
        either loaded from a source on the filesystem or is in a coherent default state.

        See also config.declare_dir and config.declare_file.
        '''
        def _check_path(path):
            if path is None:
                # check the default_value...
                if default_value is not None:
                    try: return _check_path(default_value)
                    except Exception: pass
                return None
            path = os.path.expanduser(os.path.expandvars(path))
            if not os.path.exists(path): raise ValueError('Provided file not found: %s' % path)
            path = os.path.abspath(path)
            try:              return path if filter is None else filter(path)
            except Exception: return fail_value
        return config.declare(name, rc_name=rc_name, environ_name=environ_name, filter=_check_path,
                              default_value=None, merge=merge)
    @staticmethod
    def declare_dir(name, 
                    rc_name=None, environ_name=None, filter=None, merge=None,
                    default_value=None, fail_value=None, use_temp=False):
        '''
        config.declare_dir(...) is equivalent to config.declare_path(...) except that it
          additionally requires that the path be a directory.
        '''
        def _check_dir(path):
            if not os.path.isdir(path):
                raise ValueError('Path exists but is not a directory: %s' % path)
            return path if filter is None else filter(path)
        return config.declare_path(name, rc_name=rc_name, environ_name=environ_name,
                                   filter=_check_dir, default_value=default_value, merge=merge,
                                   fail_value=fail_value)
    @staticmethod
    def declare_file(name, 
                     rc_name=None, environ_name=None, filter=None, merg=None,
                     default_value=None, fail_value=None):
        '''
        config.declare_file(...) is equivalent to config.declare_path(...) except that it
          additionally requires that the path be a file (i.e., not a directory).
        '''
        def _check_file(path):
            if not os.path.isfile(path):
                raise ValueError('Path exists but is not a file: %s' % path)
            return path if filter is None else filter(path)
        return config.declare_path(name, rc_name=rc_name, environ_name=environ_name,
                                   filter=_check_file, default_value=default_value, merge=merge,
                                   fail_value=fail_value)
    @staticmethod
    def declare_json(name,
                     rc_name=None, environ_name=None, filter=None, merge=None,
                     default_value=None, fail_value=None):
        '''
        config.declare_json(...) is equivalent to config.declare_file(...) except that it
          additionally requires that the file be a json file and it loads the file then passes the
          json object through the given filter function, if any.
        '''
        def _check_json(path):
            with open(path, 'r') as fl: dat = json.load(fl)
            return dat if filter is None else filter(dat)
        return config.declare_file(name, rc_name=rc_name, environ_name=environ_name,
                                   filter=_check_hson, default_value=default_value, merge=merge,
                                   fail_value=fail_value)
    _vals = {}
    @staticmethod
    def _getitem(self, name):
        if name not in config._vars: raise KeyError(name)
        if name not in config._vals:
            (rcname, envname, fltfn, dval, merge) = config._vars[name]
            val = dval
            rcdat = config.rc()
            # see where it's defined:
            if envname in os.environ:
                val = os.environ[envname]
                try: val = json.loads(val)
                except Exception: pass # it's a string if it can't be json'ed
                # it could be in both env and rc: use merge if needed
                if merge is not None and rcname in rcdat: val = merge(rcdat[rcname], val)
            elif rcname in rcdat: val = rcdat[rcname]
            # if there's a filter, run it
            if fltfn is not None:
                try: val = fltfn(val)
                except Exception: val = dval # failure--reset to default
            config._vals[name] = val
        return config._vals[name]
    @staticmethod
    def _setitem(self, name, val):
        if name not in config._vars:
            raise ValueError('Configurable neuropythy key "%s" not declared' % name)
        (rcname, envname, fltfn, dval, merge) = config._vars[name]
        self._vals[name] = val if fltfn is None else fltfn(val)
    @staticmethod
    def _iter(self): return six.iterkeys(self._vars)
    @staticmethod
    def _len(self): return len(self._vars)
    @staticmethod
    def keys(): return config._vars.keys()
    @staticmethod
    def values(): return map(lambda k:config[k], config.keys())
    @staticmethod
    def items(): return map(lambda k:(k,config[k]), config.keys())
    @staticmethod
    def todict(): return {k:config[k] for k in config.keys()} def load(filename, format=None, **kwargs):
    '''
    load(filename) yields the data contained in the file referenced by the given filename in a
      neuropythy or neuropythy-friendly format.
    load(filename, format) specifies that the given format should be used; this should be the name
      of the importer (though a file extension that is recognized also will work).

    Additionally, functions located in load.<format> may be used; so, for example, the following
    are equivalent calls:
      load(filename, 'nifti')
      load.nifti(filename)
    In fact, the load.nifti function is just the nifti importer, so help(load.nifti) will also
    yield documentation for the nifti importer.

    Keyword options may be passed to load; these must match those accepted by the given import
    function.
    '''
    from neuropythy.util import (ObjectWithMetaData, to_pseudo_path)
    # use the pseudo-path interface for the filename:
    pdir = to_pseudo_path(filename)
    filename = pdir.local_path()
    obj = None
    if format is None:
        # Try to guess first; if this fails, we try others...
        format = guess_import_format(filename, **kwargs)
        if format is not None:
            try:
                (f,_,_,md) = importers[format]
                obj = f(filename, **kwargs)
            except Exception: format = None
        # no good guess, try them all:
        if format is None:
            for (k,(f,_,_,md)) in six.iteritems(importers):
                if obj is not None: break
                if not bool(md.get('auto', True)): continue
                try:              obj = f(filename, **kwargs)
                except Exception: pass
        if obj is None:
            raise ValueError('Could not deduce format of file %s' % filename)
    else:
        format = format.lower()
        if format not in importers:
            raise ValueError('Format \'%s\' not recognized by neuropythy' % format)
        (f,_,_,_) = importers[format]
        obj = f(filename, **kwargs)
    if isinstance(obj, ObjectWithMetaData): obj = obj.with_meta(source_filename=filename)
    return obj def importer(name, extensions=None, sniff=None, auto=True):
    '''
    @importer(name) is a decorator that declares that the following function is an file loading
      function that should be registered with the neuropythy load function. See also the
      forget_importer function.

    Any importer function must take, as its first argument, a filename; after that it may take any
    number of keyword arguments, but no other non-keyword arguments. These keyword arguments can be
    passed to the neuropythy load function.
    
    The following options are accepted:
      * extensions (default: None) may be a string or a collection of strings that indicate possible
        file extensions for files of this type.
      * sniff (default: None) may optionally be a function f(s) that yields True when the given
        string s is a filename for a file of this type. If no sniff is given, this type can still
        be detected by running the importer and catching any raised exception.
      * auto (default: True) may optionally be set to false to specify that the importer should not
        be auto-guessed--i.e., either the sniff must yield true or the format must be explicitly
        specified for the importer to be invoked; it will never be invoked by trial and error.
    '''
    name = name.lower()
    if name in importers:
        raise ValueError('An importer for type %s already exists; see forget_importer' % name)
    if extensions is None:         extensions = ()
    elif pimms.is_str(extensions): (extensions,)
    else:                          extensions = tuple(extensions)
    def _importer(f):
        global importers
        importers = importers.set(name, (f, extensions, sniff, pyr.m(auto=auto)))
        setattr(load, name, f)
        return f
    return _importer def exporter(name, extensions=None, sniff=None):
    '''
    @exporter(name) is a decorator that declares that the following function is an file saveing
      function that should be registered with the neuropythy save function. See also the
      forget_exporter function.

    Any exporter function must take, as its first argument, a filename and, as its second argument,
    the object to be exported; after that it may take any number of keyword arguments, but no other
    non-keyword arguments. These keyword arguments can be passed to the neuropythy save function.
    
    The following options are accepted:
      * extensions (default: None) may be a string or a collection of strings that indicate possible
        file extensions for files of this type.
      * sniff (default: None) may optionally be a function f(s, d) that yields True when the given
        string s is a filename for a file of this type and/or the given object d is an object that
        can be exported as this type. If no sniff is given, this type can still be detected by
        running all exporters exhaustively and catching any raised exceptions; though this may
        result in partial files written to disk, so is not used by save.
    '''
    name = name.lower()
    if name in exporters:
        raise ValueError('An exporter for type %s already exists; use forget_exporter' % name)
    extensions = (extensions,) if pimms.is_str(extensions) else \
                 ()            if extensions is None       else \
                 tuple(extensions)
    def _exporter(f):
        global exporters
        exporters = exporters.set(name, (f, extensions, sniff))
        setattr(save, name, f)
        return f
    return _exporter
import numpy                        as np
import numpy.matlib                 as npml
import numpy.linalg                 as npla
import scipy                        as sp
import scipy.spatial                as space
import scipy.sparse                 as sps
import scipy.optimize               as spopt
import nibabel                      as nib
import nibabel.freesurfer.mghformat as fsmgh
import pyrsistent                   as pyr
import os, sys, six, types, logging, warnings, gzip, json, pimms
    import neuropythy.mri as mri
    import scipy.sparse as sps
from .util  import (triangle_area, triangle_address, alignment_matrix_3D, rotation_matrix_3D,
                    cartesian_to_barycentric_3D, cartesian_to_barycentric_2D, vector_angle_cos,
                    segment_intersection_2D, segments_overlapping, points_close, point_in_segment,
                    point_on_segment, barycentric_to_cartesian, point_in_triangle,
                    line_segment_intersection_2D, segments_touch_2D)
from ..util import (ObjectWithMetaData, to_affine, zinv, is_image, is_address, address_data, curry,
                    curve_spline, CurveSpline, chop, zdivide, flattest, inner, config, library_path,
                    dirpath_to_list, to_hemi_str, is_tuple, is_list, is_set, close_curves,
                    normalize, denormalize, AutoDict, auto_dict, times, apply_affine, nanlt)
from ..io   import (load, importer, exporter)
from functools import reduce
try:              from StringIO import StringIO
except Exception: from io import StringIO
    from neuropythy.util import (naneq, nanle, nanlt)
        from scipy.sparse import csr_matrix
        from scipy.sparse import csr_matrix
        from scipy.sparse import csr_matrix
        from scipy.sparse import csr_matrix, find
        from neuropythy.mri import image_interpolate
        from scipy.sparse import csr_matrix, find, csgraph
        from neuropythy.mri import is_subject
        from neuropythy.mri import is_subject
    from neuropythy.mri import is_subject
    from neuropythy import (is_subject, is_list, is_tuple, is_topo, is_mesh, is_tess)
    from neuropythy.util import flattest
            from neuropythy import is_cortex
    @staticmethod
    def estimate_distances(addresses, mesh):
        '''
        Path.estimate_distances(addresses, mesh) estimates all the distances between the vertices in
          the mesh and the path implied by the given set of addresses using a minimum-graph-distance
          algorithm over the mesh edges.
        '''
        # we're going to estimate an unsigned distance for every vertex, one can add the sign in for
        # a closed path using the contained_faces or similar if desired
        d = np.full(mesh.vertex_count, np.inf) # initial distances
        # for starters, we need to get distances from all the addresses
        (faces, coords) = address_data(addresses, 2)
        coords = np.vstack([coords, [1 - np.sum(coords, axis=0)]])
        faces  = mesh.tess.index(faces)
        # each face implies a minimum distance to each of its vertices
        for (f,x) in zip(faces.T, coords.T):
            fx = mesh.coordinates[:,f]
            pt = np.sum(fx * x, axis=1)
            # distances between...
            dd = np.sqrt(np.sum((fx.T - pt)**2, axis=1))
            ii = (d[f] > dd)
            d[f[ii]] = dd[ii]
        # okay, starting distances obtained; now let's do a search!
        ii0 = np.where(np.isfinite(d))[0]
        if len(ii0) == 0: raise ValueError('No distances obtained from addresses!')
        # we iterate until we cannot decrease the path-length of anything
        (max_iter,it) = (mesh.vertex_count + 1, 0)
        (u,v) = np.hstack([mesh.tess.indexed_edges, np.flipud(mesh.tess.indexed_edges)])
        elens = np.concatenate([mesh.edge_lengths, mesh.edge_lengths])
        updated = ii0
        while len(updated) > 0:
            if it > max_iter: raise ValueError('Max iterations exceeded in min-distance search loop')
            else: it += 1
            # what edges might be affected by this update
            ii = np.where(np.isin(u, updated))[0]
            (uu,vv) = (u[ii],v[ii])
            dnew = d[uu] + elens[ii]
            # we want to compare into vv, but vv will have repeats, so we must min across these;
            # easy way to do this is to put them all in a sparse matrix and min/argmin across
            (vv,ridx) = np.unique(vv, return_inverse=True)
            # we want to take the smallest value that is greater than 0, so we have to use a trick
            mx = np.max(dnew) + 1.0
            mtx = sps.csr_matrix((mx - dnew, (ridx, np.arange(len(dnew)))))
            dnew = mx - flattest(mtx.max(axis=1))
            # okay, find out where this distance is less than the previous dist
            jj = np.where(d[vv] > dnew)[0]
            updated = vv[jj]
            d[updated] = dnew[jj]
        return d
    @pimms.value
    def estimated_distances(addresses, surface):
        '''
        path.estimated_distances is a vector of estimated distance values from the given path. The
        distance estimates are upper bounds on the actual distance but are not exact.
        '''
        if is_topo(surface):
            return pimms.lazy_map(
                {k:curry(lambda a,k: Path.estimate_distances(a, surface.surfaces[k]), addresses, k)
                 for k in six.iterkeys(surface.surfaces)})
        else: return Path.estimate_distances(addresses, surface)
    @pimms.value
    def estimated_signed_distances(estimated_distances, contained_faces, surface):
        '''
        path.estimated_distances is a vector of estimated distance values from the given path. The
        distance estimates are upper bounds on the actual distance but are not exact.
        '''
        fs = surface.tess.index(np.unique(contained_faces))
        mlt = np.ones(surface.vertex_count, dtype=np.int)
        mlt[fs] = -1
        if is_topo(surface):
            return pimms.lazy_map(
                {k:curry(lambda k: estimated_distances[k]*mlt, k)
                 for k in six.iterkeys(estimated_distances)})
        else: return estimated_distances * mlt
    def reverse(self, meta_data=None):
        '''
        path.reverse() yields a path that is equivalent to the given path but reversed (thus it is
          considered to contain all the vertices not contained by path if path is closed).
        '''
        addrs = {k:np.fliplr(v) for (k,v) in six.iteritems(self.addresses)}
        return Path(self.surface, addrs, meta_data=meta_data, atol=atol)
    def boundary_vertices(self, distance, outer=False, indices=False):
        '''
        path.boundary_vertices(d) yields an array of vertex labels representing all vertices that
          are both on the on the inner/left-hand side of the given path and within a distance d of
          the boundary of the path.
        
        The following options may be given:
          * outer (default: False) may be set to True in order to use the right-hand/outer side of
            the path instead of the inner/left-hand side.
          * indices (default: False) may be set to True in order to yield indices instead of labels.
        '''
        #TODO
        raise NotImplementedError('Path.boundary_vertices is not yet implemented')
def is_path(p):
    '''
    is_path(p) yields True if p is a Path object and False otherwise.
    '''
    return isinstance(p, Path)

# Utility function for below
def bcfull(bc):
    '''
    bcfull(bc) converts the 2D or 2xN barycentric coordinate data into 3D or 3xN where the third
      dimension is simply 1 - bc[0] - bc[1].
    '''
    if len(bc) == 3: return np.asarray(bc)
    bc = np.asarray(bc)
    bc = [bc, [1 - bc[0] - bc[1]]]
    if len(bc[0].shape) == 1: bc = np.concatenate(bc)
    else: bc = np.vstack(bc)
    return bc
def bcfix(bc, atol=1e-8):
    '''
    bcfix(bc) rounds near-zero values of the given barycentric coordinates to 0 and returns the
      rescaled coordinates.
    '''
    bc = bcfull(bc)
    bc = bc * zinv(np.sum(np.abs(bc), axis=0))
    bc[np.isclose(bc, 0, atol=atol)] = 0

    return bc * zinv(np.sum(np.abs(bc), axis=0))