r"""
        Work out the voxels inside the convex hull of the voronoi vertices of
        each pore
501989
def in_hull_volume(self):
        
        """
        i = self.network['pore.internal']
        s = self.network['pore.surface']
        d = self.network['pore.delaunay']
        Ps = self.network.pores()[np.logical_and(d, np.logical_or(i, s))]
        inds = self.network._map(ids=self['pore._id'][Ps], element='pore',
                                 filtered=True)
        # Get the fiber image
        self._get_fiber_image(inds)
        hull_image = np.ones_like(self._fiber_image, dtype=np.uint16)*-1
        self._hull_image = hull_image
        for pore in Ps:
            logger.info("Processing Pore: "+str(pore+1)+" of "+str(len(Ps)))
            verts = self['pore.vertices'][pore]
            verts = np.asarray(unique_list(np.around(verts, 6)))
            verts /= self.network.resolution
            self.inhull(verts, pore)
        self._process_pore_voxels()
527802
def computeDelaunayTriangulation(points):
    
    siteList = SiteList(points)
    context  = Context()
    context.triangulate = True
    voronoi(siteList,context)
    return context.triangles
1240381
def q12d(grid, spacing=None, E=1e5, nu=0.3, dirichlet_boundary=True,
         format=None):
    
    X, Y = tuple(grid)

    if X < 1 or Y < 1:
        raise ValueError('invalid grid shape')

    if dirichlet_boundary:
        X += 1
        Y += 1

    pts = np.mgrid[0:X+1, 0:Y+1]
    pts = np.hstack((pts[0].T.reshape(-1, 1) - X / 2.0,
                     pts[1].T.reshape(-1, 1) - Y / 2.0))

    if spacing is None:
        DX, DY = 1, 1
    else:
        DX, DY = tuple(spacing)
        pts *= [DX, DY]

    # compute local stiffness matrix
    lame = E * nu / ((1 + nu) * (1 - 2*nu))  # Lame's first parameter
    mu = E / (2 + 2*nu)                   # shear modulus

    vertices = np.array([[0, 0], [DX, 0], [DX, DY], [0, DY]])
    K = q12d_local(vertices, lame, mu)

    nodes = np.arange((X+1)*(Y+1)).reshape(X+1, Y+1)
    LL = nodes[:-1, :-1]
    Id = (2*LL).repeat(K.size).reshape(-1, 8, 8)
    J = Id.copy()
    Id += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1))
    J += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1)).T
    V = np.tile(K, (X*Y, 1))

    Id = np.ravel(Id)
    J = np.ravel(J)
    V = np.ravel(V)

    # sum duplicates
    A = coo_matrix((V, (Id, J)), shape=(pts.size, pts.size)).tocsr()
    A = A.tobsr(blocksize=(2, 2))

    del Id, J, V, LL, nodes

    B = np.zeros((2 * (X+1)*(Y+1), 3))
    B[0::2, 0] = 1
    B[1::2, 1] = 1
    B[0::2, 2] = -pts[:, 1]
    B[1::2, 2] = pts[:, 0]

    if dirichlet_boundary:
        mask = np.zeros((X+1, Y+1), dtype='bool')
        mask[1:-1, 1:-1] = True
        mask = np.ravel(mask)
        data = np.zeros(((X-1)*(Y-1), 2, 2))
        data[:, 0, 0] = 1
        data[:, 1, 1] = 1
        indices = np.arange((X-1)*(Y-1))
        indptr = np.concatenate((np.array([0]), np.cumsum(mask)))
        P = bsr_matrix((data, indices, indptr),
                       shape=(2*(X+1)*(Y+1), 2*(X-1)*(Y-1)))
        Pt = P.T
        A = P.T * A * P

        B = Pt * B

    return A.asformat(format), B
520439
def euclidean3d(v1, v2):
    """"""
    if not len(v1) == 3 and len(v2) == 3:
        print("Vectors are not in 3D space. Returning None.")
        return None
    return np.sqrt((v1[0] - v2[0]) ** 2 + (v1[1] - v2[1]) ** 2 + (v1[2] - v2[2]) ** 2)
1174939
def convex_conj(self):
        
        conj_exp = conj_exponent(self.pointwise_norm.exponent)
        return GroupL1Norm(self.domain, exponent=conj_exp)
633895
def voxel_loop(self):
        ''''''
        # Prob not the most efficient, but the best I can do for now:
        for x in xrange(len(self.data)):
            for y in xrange(len(self.data[x])):
                for z in xrange(len(self.data[x][y])):
                    yield ((x,y,z),self.data[x][y][z])
1236975
def find_windows(elements,
                 coordinates,
                 processes=None,
                 mol_size=None,
                 adjust=1,
                 pore_opt=True,
                 increment=1.0,
                 **kwargs):
    """"""
    # Copy the coordinates as will perform many opertaions on them
    coordinates = deepcopy(coordinates)
    # Center of our cartesian system is always at origin
    origin = np.array([0, 0, 0])
    # Initial center of mass to reverse translation at the end
    initial_com = center_of_mass(elements, coordinates)
    # Shift the cage to the origin using either the standard center of mass
    # or if pore_opt flag is True, the optimised pore center as center of mass
    if pore_opt is True:
        # Normally the pore is calculated from the COM of a molecule.
        # So, esentially the molecule's COM is the pore center.
        # To shift the molecule so that the center of the optimised pore
        # is at the origin of the system and not the center of the not
        # optimised one, we need to adjust the shift. We also have to update
        # the initial com.
        com_adjust = initial_com - opt_pore_diameter(elements, coordinates, **
                                                     kwargs)[2]
        initial_com = initial_com - com_adjust
        coordinates = shift_com(elements, coordinates, com_adjust=com_adjust)
    else:
        # Otherwise, we just shift the cage to the origin.
        coordinates = shift_com(elements, coordinates)
    # We create an array of vdw radii of elements.
    elements_vdw = np.array([[atomic_vdw_radius[x.upper()]] for x in elements])
    # We calculate maximum diameter of a molecule to determine the radius
    # of a sampling sphere neccessary to enclose the whole molecule.
    shpere_radius = max_dim(elements, coordinates)[2] / 2
    sphere_surface_area = 4 * np.pi * shpere_radius**2
    # Here we determine the number of sampling points necessary for a fine
    # sampling. Smaller molecules require more finner density of sampling
    # points on the sampling sphere's surface, whereas largen require less.
    # This formula was created so that larger molecule do not take much longer
    # to analyse, as number_sampling_points*length_of_sampling_vectors
    # results in quadratic increase of sampling time. The 250 factor was
    # specificly determined to produce close to 1 sampling point /Angstrom^2
    # for a sphere of radius ~ 24 Angstrom. We can adjust how fine is the
    # sampling by changing the adjust factor.
    number_of_points = int(np.log10(sphere_surface_area) * 250 * adjust)
    # Here I use code by Alexandre Devert for spreading points on a sphere:
    # http://blog.marmakoide.org/?p=1
    golden_angle = np.pi * (3 - np.sqrt(5))
    theta = golden_angle * np.arange(number_of_points)
    z = np.linspace(1 - 1.0 / number_of_points, 1.0 / number_of_points - 1.0,
                    number_of_points)
    radius = np.sqrt(1 - z * z)
    points = np.zeros((number_of_points, 3))
    points[:, 0] = radius * np.cos(theta) * shpere_radius
    points[:, 1] = radius * np.sin(theta) * shpere_radius
    points[:, 2] = z * shpere_radius
    # Here we will compute the eps parameter for the sklearn.cluster.DBSCAN
    # (3-dimensional spatial clustering algorithm) which is the mean distance
    # to the closest point of all points.
    values = []
    tree = KDTree(points)
    for i in points:
        dist, ind = tree.query(i.reshape(1, -1), k=10)
        values.extend(dist)
    mean_distance = np.mean(values)
    # The best eps is parametrized when adding the mean distance and it's root.
    eps = mean_distance + mean_distance**0.5
    # Here we either run the sampling points vectors analysis in serial
    # or parallel. The vectors that go through molecular pores return
    # as analysed list with the increment at vector's path with largest
    # included sphere, coordinates for this narrow channel point. vectors
    # that find molecule on theirs path are return as NoneType object.
    # Parralel analysis on user's defined number of CPUs.
    if processes:
        pool = Pool(processes=processes)
        parallel = [
            pool.apply_async(
                vector_preanalysis,
                args=(
                    point,
                    coordinates,
                    elements_vdw, ),
                kwds={'increment': increment}) for point in points
        ]
        results = [p.get() for p in parallel if p.get() is not None]
        pool.terminate()
        # Dataset is an array of sampling points coordinates.
        dataset = np.array([x[5:8] for x in results])
    else:
        results = [
            vector_preanalysis(
                point, coordinates, elements_vdw, increment=increment)
            for point in points
        ]
        results = [x for x in results if x is not None]
        dataset = np.array([x[5:8] for x in results])
    # If not a single vector was returned from the analysis it mean that
    # no molecular channels (what we call windows here) connects the
    # molecule's interior with the surroungsings (exterior space).
    # The number of windows in that case equals zero and zero is returned.
    # Otherwise we continue our search for windows.
    if len(results) == 0:
        return None
    else:
        # Perfomr DBSCAN to cluster the sampling points vectors.
        # the n_jobs will be developed later.
        # db = DBSCAN(eps=eps, n_jobs=_ncpus).fit(dataset)
        db = DBSCAN(eps=eps).fit(dataset)
        core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
        core_samples_mask[db.core_sample_indices_] = True
        labels = set(db.labels_)
        # Assing cluster label to a sampling point.
        clusters = [[i, j] for i, j in zip(results, db.labels_)]
        clustered_results = {label: [] for label in labels}
        # Create a dictionary of clusters with points listed.
        [clustered_results[i[1]].append(i[0]) for i in clusters]
        # No for the sampling point vector in each cluster that had
        # the widest channel's 'neck' is assumed to pass the closest
        # to the window's center and therefore will be passed to
        # window analysis function.
        # We also pass user defined settings for window analysis.
        # Again either in serlia or in parallel.
        # Noisy points get a cluster label -1, therefore we have to exclude it.
        if processes:
            pool = Pool(processes=processes)
            parallel = [
                pool.apply_async(
                    window_analysis,
                    args=(np.array(clustered_results[cluster]), elements,
                          coordinates, elements_vdw),
                    kwds=kwargs) for cluster in clustered_results
                if cluster != -1
            ]
            window_results = [p.get() for p in parallel if p.get() is not None]
            pool.terminate()
        else:
            window_results = [
                window_analysis(
                    np.array(clustered_results[cluster]), elements,
                    coordinates, elements_vdw, **kwargs)
                for cluster in clustered_results if cluster != -1
            ]
    # The function returns two numpy arrays, one with windows diameters
    # in Angstrom, second with corresponding windows center's coordinates
    windows = np.array([result[0] for result in window_results
                        if result is not None])
    windows_coms = np.array(
        [np.add(result[1], initial_com) for result in window_results
         if result is not None])
    # Safety measures, if one of the windows is None or negative a warning
    # should be raised.
    for result in window_results:
        if result is None:
            msg_ = " ".join(
                ['Warning. One of the analysed windows has',
                 'returned as None. See manual.']
            )
            # print(msg_)
        elif result[0] < 0:
            msg_ = " ".join(
                ['Warning. One of the analysed windows has a vdW',
                 'corrected diameter smaller than 0. See manual.']
            )
            # print(msg_)
    return (windows, windows_coms)
1324338
def grid_coords_from_corners(upper_left_corner, lower_right_corner, size):
    '''  '''
    assert upper_left_corner.wkt == lower_right_corner.wkt
    geotransform = np.array([upper_left_corner.lon, -(upper_left_corner.lon - lower_right_corner.lon) / float(size[1]), 0,
                            upper_left_corner.lat, 0, -(upper_left_corner.lat - lower_right_corner.lat) / float(size[0])])
    return GridCoordinates(geotransform=geotransform,
                               wkt=upper_left_corner.wkt,
                               y_size=size[0],
                               x_size=size[1])
1158381
def pointsToVoronoiGridArray(lat, lon, extent=None):
    
    voronoi_centroids = _get_voronoi_centroid_array(lat, lon, extent)
    # find nodes surrounding polygon centroid
    # sort nodes in counterclockwise order
    # create polygon perimeter through nodes
    log("Building Voronoi polygons...")
    # compute voronoi
    voronoi_manager = Voronoi(voronoi_centroids)
    voronoi_vertices = voronoi_manager.vertices
    voronoi_regions = voronoi_manager.regions
    feature_list = []
    for point_id, region_index in enumerate(voronoi_manager.point_region):
        vert_index_list = np.array(voronoi_regions[region_index])
        voronoi_centroid = voronoi_centroids[point_id]
        voronoi_poly_points = _get_voronoi_poly_points(vert_index_list,
                                                       voronoi_vertices,
                                                       voronoi_centroid)

        if len(voronoi_poly_points) == 4:
            feature_list.append({'polygon': Polygon(voronoi_poly_points),
                                 'lon': voronoi_centroid[0],
                                 'lat': voronoi_centroid[1]})

    return feature_list
1062529
def _ncores(meshes, twod):
    """"""
    nnpb = len(meshes)  # number of nodes per block
    nns = [1, 1, 1]  # number of nodes in x, y, z directions
    if twod is None or 'X' in twod:
        while (nnpb > 1 and
               meshes[nns[0]]['X'][0, 0, 0] ==
               meshes[nns[0] - 1]['X'][-1, 0, 0]):
            nns[0] += 1
            nnpb -= 1
    cpu = lambda icy: icy * nns[0]
    if twod is None or 'Y' in twod:
        while (nnpb > 1 and
               meshes[cpu(nns[1])]['Y'][0, 0, 0] ==
               meshes[cpu(nns[1] - 1)]['Y'][0, -1, 0]):
            nns[1] += 1
            nnpb -= nns[0]
    cpu = lambda icz: icz * nns[0] * nns[1]
    while (nnpb > 1 and
           meshes[cpu(nns[2])]['Z'][0, 0, 0] ==
           meshes[cpu(nns[2] - 1)]['Z'][0, 0, -1]):
        nns[2] += 1
        nnpb -= nns[0] * nns[1]
    return np.array(nns)
710856
def isosurface_from_data(data, isolevel, origin, spacing):
    
    spacing = np.array(extent/resolution)
    if isolevel >= 0:
        triangles = marching_cubes(data, isolevel)
    else: # Wrong traingle unwinding roder -- god only knows why
        triangles = marching_cubes(-data, -isolevel)
    faces = []
    verts = []
    for i, t in enumerate(triangles):
       faces.append([i * 3, i * 3 +1, i * 3 + 2])
       verts.extend(t)
       
    faces = np.array(faces)
    verts = origin + spacing/2 + np.array(verts)*spacing
    return verts, faces
975442
def convex_hull(features):
    """"""
    points = sorted([s.point() for s in features])
    l = reduce(_keep_left, points, [])
    u = reduce(_keep_left, reversed(points), [])
    return l.extend(u[i] for i in xrange(1, len(u) - 1)) or l
544303
def vt2esofspy(vesseltree, outputfilename="tracer.txt", axisorder=[0, 1, 2]):
    

    if (type(vesseltree) == str) and os.path.isfile(vesseltree):
        import io3d
        vt = io3d.misc.obj_from_file(vesseltree)
    else:
        vt = vesseltree
    print(vt['general'])
    print(vt.keys())
    vtgm = vt['graph']['microstructure']
    lines = []
    vs = vt['general']['voxel_size_mm']
    sh = vt['general']['shape_px']

    # switch axis
    ax = axisorder

    lines.append("#Tracer+\n")
    lines.append("#voxelsize mm %f %f %f\n" % (vs[ax[0]], vs[ax[1]], vs[ax[2]]))
    lines.append("#shape %i %i %i\n" % (sh[ax[0]], sh[ax[1]], sh[ax[2]]))
    lines.append(str(len(vtgm) * 2)+"\n")

    i = 1
    for id in vtgm:
        # edge['']
        try:
            nda = vtgm[id]['nodeA_ZYX']
            ndb = vtgm[id]['nodeB_ZYX']
            lines.append("%i\t%i\t%i\t%i\n" % (nda[ax[0]], nda[ax[1]], nda[ax[2]], i))
            lines.append("%i\t%i\t%i\t%i\n" % (ndb[ax[0]], ndb[ax[1]], ndb[ax[2]], i))
            i += 1
        except:
            pass


    lines.append("%i\t%i\t%i\t%i" % (0, 0, 0, 0))
    lines[3] = str(i - 1) + "\n"
    with open(outputfilename, 'wt') as f:
        f.writelines(lines)
1202191
function ( inPolygonData ) {
		// positional slots for neighboring trapezoid-diagonals
		var DIAG_UL = 0, DIAG_UM = 1, DIAG_ULR = 2, DIAG_UR = 3;
		var DIAG_DR = 4, DIAG_DM = 5, DIAG_DLR = 6, DIAG_DL = 7;

		var i, j;
		var nbVertices = inPolygonData.nbVertices();

		// initialize arrays for neighboring trapezoid-diagonals and vertices
		var myVisibleDiagonals	= new Array(nbVertices);
		for ( i = 0; i < nbVertices; i++ ) {
			myVisibleDiagonals[i] = new Array(DIAG_DL+1);
		}
		// create the list of neighboring trapezoid-diagonals
		//	put into their positional slots
		var myExternalNeighbors = new Array(nbVertices);
		for ( i = 0, j = this.trapArray.length; i < j; i++ ) {
			var curTrap = this.trapArray[i];
			var highPos = curTrap.uL ?
						( curTrap.uR ? DIAG_DM : DIAG_DL ) :
						( curTrap.uR ? DIAG_DR : DIAG_DLR );
			var lowPos = curTrap.dL ?
						( curTrap.dR ? DIAG_UM : DIAG_UL ) :
						( curTrap.dR ? DIAG_UR : DIAG_ULR );

			if ( ( curTrap.depth % 2 ) == 1 ) {		// inside ?
				if ( ( highPos == DIAG_DM ) || ( lowPos == DIAG_UM ) ||
					 ( ( highPos == DIAG_DL ) && ( lowPos == DIAG_UR ) ) ||
					 ( ( highPos == DIAG_DR ) && ( lowPos == DIAG_UL ) ) ) {
					var lhDiag = inPolygonData.appendDiagonalsEntry( {
									vFrom: curTrap.vLow, vTo: curTrap.vHigh,
									mprev: null, mnext: null, marked: false } );
					var hlDiag = inPolygonData.appendDiagonalsEntry( {
									vFrom: curTrap.vHigh, vTo: curTrap.vLow, revDiag: lhDiag,
									mprev: null, mnext: null, marked: false } );
					lhDiag.revDiag = hlDiag;
					myVisibleDiagonals[ curTrap.vLow.id][ lowPos] = lhDiag;
					myVisibleDiagonals[curTrap.vHigh.id][highPos] = hlDiag;
				}
			} else {		// outside, hole
				if ( curTrap.vHigh.id !== null )	myExternalNeighbors[curTrap.vHigh.id] = highPos;
				if ( curTrap.vLow.id  !== null )	myExternalNeighbors[ curTrap.vLow.id] = lowPos;
			}
		}
		// create the list of outgoing diagonals in the right order (CW)
		//	from the ordered list of neighboring trapezoid-diagonals
		//	- starting from an external one
		// and connect those incoming to
		var curDiag, curDiags, firstElem, fromVertex, lastIncoming;
		for ( i = 0; i < nbVertices; i++ ) {
			curDiags  = myVisibleDiagonals[i];
			firstElem = myExternalNeighbors[i];
			if ( firstElem == null )	continue;		// eg. skipped vertices (zero length, co-linear		// NOT: === !
			j = firstElem;
			lastIncoming = null;
			do {
				if ( j++ > DIAG_DL )			j = DIAG_UL;	// circular positional list
				if ( curDiag = curDiags[j] ) {
					if ( lastIncoming ) {
						curDiag.mprev = lastIncoming;
						lastIncoming.mnext = curDiag;
					} else {
						fromVertex = curDiag.vFrom;
						fromVertex.firstOutDiag = curDiag;
					}
					lastIncoming = curDiag.revDiag;
				}
			} while ( j != firstElem );
			if ( lastIncoming )		fromVertex.lastInDiag = lastIncoming;
		}
	}
50359
def convex_hull_image(image):
    ''''''
    labels = image.astype(int)
    points, counts = convex_hull(labels, np.array([1]))
    output = np.zeros(image.shape, int)
    for i in range(counts[0]):
        inext = (i+1) % counts[0]
        draw_line(output, points[i,1:], points[inext,1:],1)
    output = fill_labeled_holes(output)
    return output == 1
368544
def direction_cossim(layer, vec, batch=None):
  """"""
  def inner(T):
    act_mags = tf.sqrt(tf.reduce_sum(T(layer)**2, -1, keepdims=True))
    vec_mag = tf.sqrt(tf.reduce_sum(vec**2))
    mags = act_mags * vec_mag
    if batch is None:
      return tf.reduce_mean(T(layer) * vec.reshape([1, 1, 1, -1]) / mags)
    else:
      return tf.reduce_mean(T(layer)[batch] * vec.reshape([1, 1, -1]) / mags)
  return inner
1237012
def _interpolate_cube(self, lon, lat, egy=None, interp_log=True):
        

        shape = np.broadcast(lon, lat, egy).shape
        lon = lon * np.ones(shape)
        lat = lat * np.ones(shape)
        theta = np.pi / 2. - np.radians(lat)
        phi = np.radians(lon)
        vals = []
        for i, _ in enumerate(self.hpx.evals):
            v = hp.pixelfunc.get_interp_val(self.counts[i], theta,
                                            phi, nest=self.hpx.nest)
            vals += [np.expand_dims(np.array(v, ndmin=1), -1)]

        vals = np.concatenate(vals, axis=-1)

        if egy is None:
            return vals.T

        egy = egy * np.ones(shape)

        if interp_log:
            xvals = utils.val_to_pix(np.log(self.hpx.evals), np.log(egy))
        else:
            xvals = utils.val_to_pix(self.hpx.evals, egy)

        vals = vals.reshape((-1, vals.shape[-1]))
        xvals = np.ravel(xvals)
        v = map_coordinates(vals, [np.arange(vals.shape[0]), xvals],
                            order=1)
        return v.reshape(shape)