Python scipy.spatial() Examples

def distance_from_goal(self, goal: np.ndarray, state: dict) -> float:
        """
        Given a state, check its distance from the goal

        :param goal: a numpy array representing the goal
        :param state: a dict representing the state
        :return: the distance from the goal
        """
        state_value = self.goal_from_state(state)

        # calculate distance
        if self.distance_metric == self.DistanceMetric.Cosine:
            dist = scipy.spatial.distance.cosine(goal, state_value)
        elif self.distance_metric == self.DistanceMetric.Euclidean:
            dist = scipy.spatial.distance.euclidean(goal, state_value)
        elif self.distance_metric == self.DistanceMetric.Manhattan:
            dist = scipy.spatial.distance.cityblock(goal, state_value)
        elif callable(self.distance_metric):
            dist = self.distance_metric(goal, state_value)
        else:
            raise ValueError("The given distance metric for the goal is not valid.")

        return dist 

def _in_hull(p, hull):
  ''' 
  Tests if points in `p` are in the convex hull made up by `hull`
  '''
  dim = p.shape[1]
  # if there are not enough points in `hull` to form a simplex then 
  # return False for each point in `p`.
  if hull.shape[0] <= dim:
    return np.zeros(p.shape[0], dtype=bool)
  
  if dim >= 2:
    hull = scipy.spatial.Delaunay(hull)
    return hull.find_simplex(p)>=0
  else:
    # one dimensional points
    min = np.min(hull)
    max = np.max(hull)
    return (p[:, 0] >= min) & (p[:, 0] <= max) 

def project_to_sphere(points, center, radius):
    """
    Projects the elements of points onto the sphere defined
    by center and radius.

    Parameters
    ----------
    points : array of floats of shape (npoints, ndim)
             consisting of the points in a space of dimension ndim
    center : array of floats of shape (ndim,)
            the center of the sphere to project on
    radius : float
            the radius of the sphere to project on

    returns: array of floats of shape (npoints, ndim)
            the points projected onto the sphere
    """

    lengths = scipy.spatial.distance.cdist(points, np.array([center]))
    return (points - center) / lengths * radius + center 

def project_to_sphere(points, center, radius):
    """
    Projects the elements of points onto the sphere defined
    by center and radius.

    Parameters
    ----------
    points : array of floats of shape (npoints, ndim)
             consisting of the points in a space of dimension ndim
    center : array of floats of shape (ndim,)
            the center of the sphere to project on
    radius : float
            the radius of the sphere to project on

    returns: array of floats of shape (npoints, ndim)
            the points projected onto the sphere
    """

    lengths = scipy.spatial.distance.cdist(points, np.array([center]))
    return (points - center) / lengths * radius + center 

def _compute_gram_matrix(X, kernel_type, params):
    """
    """
    if kernel_type == 'rbf':
        if 'gamma' in params:
            gamma = params['gamma']
        else:
            gamma = 1.0 / X.shape[1]

        pairwise_dist = scipy.spatial.distance.pdist(X, metric='sqeuclidean')
        pairwise_dist = scipy.spatial.distance.squareform(pairwise_dist)
        gram_matrix = np.exp( - gamma * pairwise_dist )

        np.fill_diagonal(gram_matrix, 1)
    else:
        pass

    return(gram_matrix) 

def distance_compare_unit(data):
        host = data[0]
        host_words = host['gt'].split(" ")
        host_mean = host['mu']
        host_mean = host_mean / LA.norm(host_mean)

        guest = data[1:]
        num_guest = len(guest)
        bag = []
        for idx, g in enumerate(guest):
            g_mean = g['mu']
            g_words = g['gt'].split(" ")
            # jaccard = comp_jaccard_distance(g_words, host_words)
            cos_distance = -1 * (scipy.spatial.distance.cosine(g_mean, host_mean) - 1)
            bag.append([0, cos_distance])
        return bag 

def distance_compare_unit(data):
        host = data[0]
        host_words = host['gt'].split(" ")
        host_mean = host['mean']
        host_mean = host_mean / LA.norm(host_mean)

        guest = data[1:]
        num_guest = len(guest)
        bag = []
        for idx, g in enumerate(guest):
            g_mean = g['mean']
            g_words = g['gt'].split(" ")
            # jaccard = comp_jaccard_distance(g_words, host_words)
            cos_distance = -1 * (scipy.spatial.distance.cosine(g_mean, host_mean) - 1)
            bag.append([0, cos_distance])
        return bag 

def project_to_sphere(points, center, radius):
    """
    Projects the elements of points onto the sphere defined
    by center and radius.

    Parameters
    ----------
    points : array of floats of shape (npoints, ndim)
             consisting of the points in a space of dimension ndim
    center : array of floats of shape (ndim,)
            the center of the sphere to project on
    radius : float
            the radius of the sphere to project on

    returns: array of floats of shape (npoints, ndim)
            the points projected onto the sphere
    """

    lengths = scipy.spatial.distance.cdist(points, np.array([center]))
    return (points - center) / lengths * radius + center 

def get_extend_hull(self):
        ext_points = []
        # 点集转换为numpy数组
        points = np.array([self.lons, self.lats]).T
        if len(points) < 3:
            return ext_points
        # 获取点集的凸多边形轮廓
        hull = scipy.spatial.ConvexHull(points)

        for simplex in hull.simplices:
            # 设置初值 以获得两个解
            if simplex[1] == 0:
                continue
            pairs = [True, False]
            for pair in pairs:
                # print(pair)
                extend_point = self.equations(points[simplex], pair)
                # 在边界内的点排除
                if extend_point and not self.point_in_path(extend_point):
                    ext_points.append([extend_point[0], extend_point[1], self.zvalues[simplex[0]]])
        return ext_points 

def dHdxy(self):
        """
        Calculate the spatial derivative of water depth in each direction
        (xi and eta).

        Parameters
        ----------
        None

        Returns
        -------
        dHdxi : ndarray,
          Slope in x-direction
        dHdeta : ndarray,
          Slope in eta-direction
        """
        dHdxi = np.zeros(self.h.shape)
        dHdeta = np.zeros(self.h.shape)
        dHdxi[:, :-1] = -np.diff(self.h, axis=1) * self.pm[:, 1:]
        dHdxi[:, -1] = dHdxi[:, -2]
        dHdeta[:-1, :] = -np.diff(self.h, axis=0) * self.pn[1:, :]
        dHdeta[-1, :] = dHdeta[-2, :]

        return dHdxi, dHdeta 

def calc_dist_mat(subject, bipolar=False, snap=False):
    from scipy.spatial import distance

    pos_fname = 'electrodes{}_{}positions.npz'.format('_bipolar' if bipolar else '', 'snap_' if snap else '')
    pos_fname = op.join(MMVT_DIR, subject, 'electrodes', pos_fname)
    output_fname = 'electrodes{}_{}dists.npy'.format('_bipolar' if bipolar else '', 'snap_' if snap else '')
    output_fname = op.join(MMVT_DIR, subject, 'electrodes', output_fname)

    if not op.isfile(pos_fname):
        return False
    d = np.load(pos_fname)
    pos = d['pos']
    x = np.zeros((len(pos), len(pos)))
    for i in range(len(pos)):
        for j in range(len(pos)):
            x[i,j] = np.linalg.norm(pos[i]- pos[j]) # np.sqrt((pos[i]- pos[j])**2)
            # assert(x[i,j]==np.linalg.norm(pos[i]- pos[j]))
    # x = distance.cdist(pos, pos, 'euclidean')
    np.save(output_fname, x)
    return op.isfile(output_fname) 

def solve_dist_mat(residue_position_args):
    residue_a, all_residues = residue_position_args
    ret = list()
    for residue_b in all_residues:
        dists = scipy.spatial.distance.cdist(residue_a, residue_b)
        ret.append(np.min(dists))
    return ret 

def nearest(self, lon, lat, grid="rho"):
        """
        Find the indices nearest to each point in the given list of
        longitudes and latitudes.

        Parameters
        ----------
        lon : ndarray,
            longitude of points to find
        lat : ndarray
            latitude of points to find
        grid : string, optional,
            "rho", "u", or "v" grid to search

        Returns
        -------
        indices : tuple of ndarray
            The indices for each dimension of the grid that are closest
            to the lon/lat points specified
        """

        glat = getattr(self, "lat_" + grid)
        glon = getattr(self, "lon_" + grid)
        xy = np.dstack([glat.ravel(), glon.ravel()])[0]
        pts = np.dstack([np.atleast_1d(lat), np.atleast_1d(lon)])[0]
        grid_tree = scipy.spatial.cKDTree(xy)
        dist, idx = grid_tree.query(pts)
        return np.unravel_index(idx, glat.shape) 

def build_truth_for_one(self, id):
        index_reference, sequence, structure = self.get_index_ref_seq_and_struct(id)

        residue_positions = GroundTruthBuilder.get_residue_positions(structure)
        pdb_dist_mat = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(residue_positions, 'euclidean'))

        ret_mat = GroundTruthBuilder.get_fixed_mat(pdb_dist_mat, index_reference, len(sequence))
        ret_mat = self.fix_astral_206_ground_truth(id, ret_mat)
        return ret_mat 

def in_hull(p, hull):
    """
    :param p: (N, K) test points
    :param hull: (M, K) M corners of a box
    :return (N) bool
    """
    try:
        if not isinstance(hull, Delaunay):
            hull = Delaunay(hull)
        flag = hull.find_simplex(p) >= 0
    except scipy.spatial.qhull.QhullError:
        print('Warning: not a hull %s' % str(hull))
        flag = np.zeros(p.shape[0], dtype=np.bool)

    return flag 

def _cal_silhouette_score(X, y, sample_size, metric):
    """
    """
    n_samples = X.shape[0]
    if sample_size is not None:
       assert(1 < sample_size < n_samples)
       rand_inx = np.random.choice(X.shape[0], size=sample_size)
       X_samp, y_samp = X[rand_inx], y[rand_inx]

    else:
       X_samp, y_samp = X, y

    n_clust_samp = len(np.unique(y_samp))
    pair_dist = scipy.spatial.distance.pdist(X_samp, metric=metric)

    pair_dist_matrix = scipy.spatial.distance.squareform(pair_dist)

    ## Compute average intra-cluster distances 
    ## for each of the selected samples
    a_arr = _intra_cluster_distances(pair_dist_matrix, y_samp, np.arange(y_samp.shape[0]))

    ## Compute Avg. Distabces to the Neighboring Clusters
    b_arr =  _neighboring_cluster_distances(pair_dist_matrix, y_samp, np.arange(y_samp.shape[0]))

    comb_arr = np.vstack((a_arr, b_arr))
    return((b_arr-a_arr)/np.max(comb_arr, axis=0), a_arr, b_arr) 

def query(self, centroids):
        if self.entity_neighbors is not None:
            distances, indices = self.entity_neighbors.kneighbors(centroids)

            return distances, indices
        else:
            pairwise_distances = scipy.spatial.distance.cdist(
                centroids, self.entity_representations,
                metric=self.entity_representation_distance)

            distances = np.sort(pairwise_distances, axis=1)
            indices = pairwise_distances.argsort(axis=1)\
                .argsort(axis=1).argsort(axis=1)

            return distances, indices 

def getHausdorff(testImage, resultImage):
    """Compute the Hausdorff distance."""
        
    # Edge detection is done by ORIGINAL - ERODED, keeping the outer boundaries of lesions. Erosion is performed in 2D
    eTestImage   = sitk.BinaryErode(testImage, (1,1,0) )
    eResultImage = sitk.BinaryErode(resultImage, (1,1,0) )
    
    hTestImage   = sitk.Subtract(testImage, eTestImage)
    hResultImage = sitk.Subtract(resultImage, eResultImage)    
    
    hTestArray   = sitk.GetArrayFromImage(hTestImage)
    hResultArray = sitk.GetArrayFromImage(hResultImage)   
        
    # Convert voxel location to world coordinates. Use the coordinate system of the test image
    # np.nonzero   = elements of the boundary in numpy order (zyx)
    # np.flipud    = elements in xyz order
    # np.transpose = create tuples (x,y,z)
    # testImage.TransformIndexToPhysicalPoint converts (xyz) to world coordinates (in mm)
    testCoordinates   = np.apply_along_axis(testImage.TransformIndexToPhysicalPoint, 1, np.transpose( np.flipud( np.nonzero(hTestArray) )).astype(int) )
    resultCoordinates = np.apply_along_axis(testImage.TransformIndexToPhysicalPoint, 1, np.transpose( np.flipud( np.nonzero(hResultArray) )).astype(int) )
            
    # Use a kd-tree for fast spatial search
    def getDistancesFromAtoB(a, b):    
        kdTree = scipy.spatial.KDTree(a, leafsize=100)
        return kdTree.query(b, k=1, eps=0, p=2)[0]
    
    # Compute distances from test to result; and result to test
    dTestToResult = getDistancesFromAtoB(testCoordinates, resultCoordinates)
    dResultToTest = getDistancesFromAtoB(resultCoordinates, testCoordinates)    
    
    return max(np.percentile(dTestToResult, 95), np.percentile(dResultToTest, 95)) 

def __init__(self, points: np.array, epsilon: float,
                 cut_off: Optional[float] = None,
                 num_eigenpairs: Optional[int] = default.num_eigenpairs,
                 normalize_kernel: Optional[bool] = True,
                 renormalization: Optional[float] = default.renormalization,
                 kdtree_options: Optional[Dict] = None,
                 use_cuda: Optional[bool] = default.use_cuda) -> None:
        self.points = points
        self.epsilon = epsilon

        if cut_off is not None:
            import warnings
            warnings.warn('A cut off was specified for dense diffusion maps.')

        distance_matrix_squared = scipy.spatial.distance.pdist(points, metric='sqeuclidean')  # noqa
        distance_matrix_squared = scipy.spatial.distance.squareform(distance_matrix_squared)
        kernel_matrix = np.exp(-distance_matrix_squared / (2.0 * epsilon))
        if normalize_kernel is True:
            kernel_matrix = self.normalize_kernel_matrix(kernel_matrix,
                                                         renormalization)
        self.kernel_matrix = kernel_matrix
        self.renormalization = renormalization if normalize_kernel else None

        if use_cuda is True:
            import warnings
            warnings.warn('Dense diffusion maps are not implemented on the '
                          'GPU. Using the CPU instead.')

        ew, ev = self.solve_eigenproblem(self.kernel_matrix, num_eigenpairs,
                                         use_cuda=False)
        if np.linalg.norm(ew.imag > 1e2 * sys.float_info.epsilon, np.inf):
            raise ValueError('Eigenvalues have non-negligible imaginary part')
        self.eigenvalues = ew.real
        self.eigenvectors = ev.real 

def compute_kdtree(points: np.array, kdtree_options: Optional[Dict]) \
            -> None:
        """Compute kd-tree from points.

        """
        if kdtree_options is None:
            kdtree_options = dict()

        return scipy.spatial.cKDTree(points, **kdtree_options) 

def get_irsa_catalog(ra=165.86, dec=34.829694, tab=None, radius=3, catalog='allwise_p3as_psd', wise=False, twomass=False, ROW_LIMIT=500000, TIMEOUT=3600):
    """Query for objects in the `AllWISE <http://wise2.ipac.caltech.edu/docs/release/allwise/>`_ source catalog

    Parameters
    ----------
    ra, dec : float
        Center of the query region, decimal degrees

    radius : float
        Radius of the query, in arcmin

    Returns
    -------
    table : `~astropy.table.Table`
        Result of the query

    """
    from astroquery.irsa import Irsa
    Irsa.ROW_LIMIT = ROW_LIMIT
    Irsa.TIMEOUT = TIMEOUT

    #all_wise = 'wise_allwise_p3as_psd'
    #all_wise = 'allwise_p3as_psd'
    if wise:
        catalog = 'allwise_p3as_psd'
    elif twomass:
        catalog = 'fp_psc'

    coo = coord.SkyCoord(ra*u.deg, dec*u.deg)

    table = Irsa.query_region(coo, catalog=catalog, spatial="Cone",
                              radius=radius*u.arcmin, get_query_payload=False)

    return table 

def spatial_hash(coordinates):
        import scipy, scipy.spatial
        try:              from scipy.spatial import cKDTree as shash
        except Exception: from scipy.spatial import KDTree  as shash
        return shash(coordinates)
    # Static helper function 

def cdist(A, B, **kwargs):
    return scipy.spatial.distance.cdist(A, B, **kwargs) 

def _distance_mahalanobis(X=None):
    cov = X.cov().values
    cov = scipy.linalg.inv(cov)

    col_means = X.mean().values

    dist = np.full(len(X), np.nan)
    for i in range(len(X)):
        dist[i] = scipy.spatial.distance.mahalanobis(X.iloc[i, :].values, col_means, cov) ** 2
    return dist 

def _find_closest_surface_pos(pos, mesh, surface_tags):
    nodes_in_surface = _get_nodes_in_surface(mesh, surface_tags)
    if len(nodes_in_surface) == 0:
        raise ValueError('Surface with tags: {0} not found'.format(surface_tags))
    kd_tree = scipy.spatial.cKDTree(mesh.nodes[nodes_in_surface])
    _, center_idx = kd_tree.query(pos)
    #center_idx = np.argmin(np.linalg.norm(mesh.nodes[nodes_in_surface] - pos))
    pos = mesh.nodes.node_coord[nodes_in_surface - 1][center_idx]
    return pos 

def _calc_kdtree(triangles, nodes):
    tr_baricenters = _triangle_baricenters(triangles, nodes)
    return scipy.spatial.cKDTree(tr_baricenters) 

def in_hull(p, hull):
    """
    :param p: (N, K) test points
    :param hull: (M, K) M corners of a box
    :return (N) bool
    """
    try:
        if not isinstance(hull, Delaunay):
            hull = Delaunay(hull)
        flag = hull.find_simplex(p) >= 0
    except scipy.spatial.qhull.QhullError:
        print('Warning: not a hull %s' % str(hull))
        flag = np.zeros(p.shape[0], dtype=np.bool)

    return flag 

def clip_lists(input, output, clip=20):
    """TBD

    Clip [x,y] arrays of objects that don't have a match within `clip` pixels
    in either direction
    """
    import scipy.spatial

    tree = scipy.spatial.cKDTree(input, 10)

    # Forward
    N = output.shape[0]
    dist, ix = np.zeros(N), np.zeros(N, dtype=int)
    for j in range(N):
        dist[j], ix[j] = tree.query(output[j, :], k=1,
                                    distance_upper_bound=np.inf)

    ok = dist < clip
    out_arr = output[ok]
    if ok.sum() == 0:
        print('No matches within `clip={0:f}`'.format(clip))
        return False

    # Backward
    tree = scipy.spatial.cKDTree(out_arr, 10)

    N = input.shape[0]
    dist, ix = np.zeros(N), np.zeros(N, dtype=int)
    for j in range(N):
        dist[j], ix[j] = tree.query(input[j, :], k=1,
                                    distance_upper_bound=np.inf)

    ok = dist < clip
    in_arr = input[ok]

    return in_arr, out_arr 

def update_coverage_of_endpoints(self, map, points):
        self.points = points
        self.__vor = scipy.spatial.Voronoi(self.points)
        self.regions, self.vertices = self.voronoi_finite_polygons_2d(self.__vor)
        self.cells = [map.to_pixels(self.vertices[region]) for region in self.regions]
        cmap = plt.cm.Set3
        self.colors_cells = cmap(np.linspace(0., 1., len(self.points)))[:, :3] 

def drizzle_footprint(weight_image, shrink=10, ext=0, outfile=None, label=None):
    """
    Footprint of image pixels where values > 0.  Works best with drizzled
    weight images.
    """
    from scipy.spatial import ConvexHull

    im = pyfits.open(weight_image)
    wcs = pywcs.WCS(im[ext].header, fobj=im)
    sh = np.array(im[ext].data.shape)//shrink

    yp, xp = np.indices(tuple(sh))*shrink
    nonzero = im[ext].data[yp, xp] > 0

    h = ConvexHull(np.array([xp[nonzero], yp[nonzero]]).T)
    hx = xp[nonzero][h.vertices]
    hy = yp[nonzero][h.vertices]

    hrd = wcs.all_pix2world(np.stack([hx, hy]).T, 0)

    pstr = 'polygon('+','.join(['{0:.6f}'.format(i) for i in hrd.flatten()])+')'
    if label is not None:
        pstr += ' # text={{{0}}}'.format(label)

    if outfile is None:
        return pstr

    fp = open(outfile, 'w')
    fp.write('fk5\n')

    fp.write(pstr+'\n')
    fp.close()