def isolines(obj, prop, val,<EOL>outliers=None,  data_range=None,  clipped=np.inf,<EOL>weights=None,   weight_min=<NUM_LIT:0>,     weight_transform=Ellipsis,<EOL>mask=None,      valid_range=None, transform=None,<EOL>smooth=False,   yield_addresses=False):
import scipy.sparse as sps<EOL>from neuropythy import (is_subject, is_list, is_tuple, is_topo, is_mesh, is_tess)<EOL>from neuropythy.util import flattest<EOL>if is_subject(obj):<EOL><INDENT>kw = dict(outliers=outliers,     data_range=data_range,<EOL>clipped=clipped,       weights=weights,<EOL>weight_min=weight_min, weight_transform=weight_transform,<EOL>mask=mask,             valid_range=valid_range,<EOL>transform=transform,   smooth=smooth,<EOL>yield_addresses=yield_addresses)<EOL>return pimms.lazy_map({h:curry(lambda h: isolines(mesh.hemis[h],prop,val,**kw),h)<EOL>for h in six.iterkeys(mesh.hemis)})<EOL><DEDENT>elif not (is_topo(obj) or is_mesh(obj)):<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>if   smooth is True:          smooth = ()<EOL>elif smooth in [False,None]:  smooth = None<EOL>elif pimms.is_vector(smooth): smooth = tuple(smooth)<EOL>else: raise ValueError('<STR_LIT>')<EOL>fs = obj.tess.indexed_faces if not is_tess(obj) else obj.indexed_faces<EOL>N  = obj.vertex_count<EOL>p  = obj.property(prop,<EOL>outliers=outliers,     data_range=data_range,<EOL>clipped=clipped,       weights=weights,<EOL>weight_min=weight_min, weight_transform=weight_transform,<EOL>mask=mask,             valid_range=valid_range,<EOL>transform=transform)<EOL>ii  = np.isfinite(p)<EOL>fii = np.where(np.sum([ii[f] for f in fs], axis=<NUM_LIT:0>) == <NUM_LIT:3>)[<NUM_LIT:0>]<EOL>fs  = fs[:,fii]<EOL>fp  = np.asarray([p[f] for f in fs])<EOL>lt  = (fp <= val)<EOL>slt = np.sum(lt, <NUM_LIT:0>)<EOL>(ii1,ii2) = [(slt == k) for k in (<NUM_LIT:1>,<NUM_LIT:2>)]<EOL>wmtx = sps.lil_matrix((N, N), dtype='<STR_LIT:float>')<EOL>emtx = sps.dok_matrix((N*N, N*N), dtype='<STR_LIT:bool>')<EOL>for (ii,fl) in zip([ii1,ii2], [True,False]):<EOL><INDENT>tmp = np.array(lt.T)<EOL>tmp[~ii,:] = False<EOL>if fl:<EOL><INDENT>w = fs.T[tmp]<EOL>tmp[ii] = ~tmp[ii]<EOL>(u,v) = np.reshape(fs.T[tmp], (-<NUM_LIT:1>,<NUM_LIT:2>)).T<EOL><DEDENT>else:<EOL><INDENT>(u,v) = np.reshape(fs.T[tmp], (-<NUM_LIT:1>,<NUM_LIT:2>)).T<EOL>tmp[ii] = ~tmp[ii]<EOL>w = fs.T[tmp]<EOL><DEDENT>(pu,pv,pw) = [p[q] for q in (u,v,w)]<EOL>(wu, wv) = [(val - pw) / (px - pw) for px in (pu,pv)]<EOL>wmtx[u,w] = wu<EOL>wmtx[v,w] = wv<EOL>wmtx[w,u] = (<NUM_LIT:1> - wu)<EOL>wmtx[w,v] = (<NUM_LIT:1> - wv)<EOL>e1 = np.sort([u,w],axis=<NUM_LIT:0>).T.dot([N,<NUM_LIT:1>])<EOL>e2 = np.sort([v,w],axis=<NUM_LIT:0>).T.dot([N,<NUM_LIT:1>])<EOL>emtx[e1, e2] = True<EOL>emtx[e2, e1] = True<EOL><DEDENT>wmtx = wmtx.tocsr()<EOL>(rs,cs,xs) = sps.find(emtx)<EOL>((rs,ridx),(cs,cidx)) = [np.unique(x, return_inverse=True) for x in (rs,cs)]<EOL>assert(len(ridx) == len(cidx) and len(cidx) == len(xs) and np.array_equal(rs,cs))<EOL>em = sps.csr_matrix((xs, (ridx,cidx)), dtype='<STR_LIT:bool>')<EOL>ecounts = flattest(em.sum(axis=<NUM_LIT:1>))<EOL>ccs = sps.csgraph.connected_components(em, directed=False)[<NUM_LIT:1>]<EOL>lines = {}<EOL>loops = {}<EOL>for (lbl,st) in zip(*np.unique(ccs, return_index=True)):<EOL><INDENT>if lbl < <NUM_LIT:0>: continue<EOL>o0 = sps.csgraph.depth_first_order(em, st, directed=False)[<NUM_LIT:0>]<EOL>o = rs[o0]<EOL>(u,v) = (np.floor(o/N).astype('<STR_LIT:int>'), np.mod(o, N).astype('<STR_LIT:int>'))<EOL>w = flattest(wmtx[u,v])<EOL>if ecounts[o0[-<NUM_LIT:1>]] == <NUM_LIT:2>:<EOL><INDENT>loops[lbl] = (u,v,w)<EOL><DEDENT>else:<EOL><INDENT>if ecounts[st] != <NUM_LIT:1>:<EOL><INDENT>o0 = sps.csgraph.depth_first_order(em, o0[-<NUM_LIT:1>], directed=False)[<NUM_LIT:0>]<EOL>o = rs[o0]<EOL>(u,v) = (np.floor(o/N).astype('<STR_LIT:int>'), np.mod(o, N).astype('<STR_LIT:int>'))<EOL>w = flattest(wmtx[u,v])<EOL><DEDENT>lines[lbl] = (u,v,w)<EOL><DEDENT><DEDENT>addrs = []<EOL>for (lbl,(us,vs,ws)) in six.iteritems(lines):<EOL><INDENT>(u0s,v0s,u1s,v1s) = (us[:-<NUM_LIT:1>],vs[:-<NUM_LIT:1>],us[<NUM_LIT:1>:],vs[<NUM_LIT:1>:])<EOL>qs = [np.setdiff1d([u1,v1], [u0,v0])[<NUM_LIT:0>] for (u0,v0,u1,v1) in zip(u0s,v0s,u1s,v1s)]<EOL>fs = np.asarray([u0s, qs, v0s], dtype=np.int)<EOL>fend = (u1s[-<NUM_LIT:1>], np.setdiff1d(fs[:,-<NUM_LIT:1>], (u1s[-<NUM_LIT:1>],v1s[-<NUM_LIT:1>])), v1s[-<NUM_LIT:1>])<EOL>fs = np.hstack([fs, np.reshape(fend, (<NUM_LIT:3>,<NUM_LIT:1>))])<EOL>fs = np.asarray([obj.labels[f] for f in fs])<EOL>ws = np.asarray([ws, <NUM_LIT:0>*ws])<EOL>addrs.append({'<STR_LIT>': fs, '<STR_LIT>': ws})<EOL><DEDENT>addrs = list(sorted(addrs, key=lambda a:a['<STR_LIT>'].shape[<NUM_LIT:0>]))<EOL>if yield_addresses or not is_mesh(obj): return addrs<EOL>xs = [obj.unaddress(addr) for addr in addrs]<EOL>if smooth is not None: xs = smooth_lines(xs, *smooth)<EOL>return xs<EOL>
isolines(mesh, prop, val) yields a 2 x D x N array of points that represent the lines of
  equality of the given property to the given value, according to linear interpolation
  of the property values over the mesh; D represents the dimensions of the mesh (2D or
  3D) and N is the number of faces whose prop values surround val.
isolines(cortex, prop, val) yields the lines as addresses instead of coordinates.
isolines(subject, prop, val) yields a lazy map whose keys are the keys of subject.hemis
  and whose values are equivalent to isolines(subject.hemis[key], prop, val).

The following optional arguments may be given:
  * yield_addresses (default: False) may be set to True to instruct isolines to return the
    addresses instead of the coordinates, even for meshes; if the first argument is not a
    mesh then this argument is ignored.
  * smooth (default: None) may optionally be True or a value n or a tuple (n,m) that are
    used as arguments to the smooth_lines() function.
  * Almost all options that can be passed to the to_property() function can be passed to
    this function and are in turn passed along to to_property():
      * outliers       * data_range          * clipped        * weights
      * weight_min     * weight_transform    * mask           * valid_range
      * transform
    Weights, if specified, are not used in the isoline calculation; they are only used for
    thresholding.