Python spatial.cKDTree方法代码示例

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def mi(x, y, k=3, base=2):
  """Mutual information of x and y.
  x,y should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
  if x is a one-dimensional scalar and we have four samples.
  """
  assert len(x)==len(y), 'Lists should have same length.'
  assert k <= len(x) - 1, 'Set k smaller than num samples - 1.'
  intens = 1e-10 # Small noise to break degeneracy, see doc.
  x = [list(p + intens*nr.rand(len(x[0]))) for p in x]
  y = [list(p + intens*nr.rand(len(y[0]))) for p in y]
  points = zip2(x,y)
  # Find nearest neighbors in joint space, p=inf means max-norm.
  tree = ss.cKDTree(points)
  dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points]
  a = avgdigamma(x,dvec)
  b = avgdigamma(y,dvec)
  c = digamma(k)
  d = digamma(len(x))
  return (-a-b+c+d) / log(base) 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def cmi(x, y, z, k=3, base=2):
  """Mutual information of x and y, conditioned on z
  x,y,z should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]]
  if x is a one-dimensional scalar and we have four samples
  """
  assert len(x)==len(y), 'Lists should have same length.'
  assert k <= len(x) - 1, 'Set k smaller than num samples - 1.'
  intens = 1e-10 # Small noise to break degeneracy, see doc.
  x = [list(p + intens*nr.rand(len(x[0]))) for p in x]
  y = [list(p + intens*nr.rand(len(y[0]))) for p in y]
  z = [list(p + intens*nr.rand(len(z[0]))) for p in z]
  points = zip2(x,y,z)
  # Find nearest neighbors in joint space, p=inf means max-norm.
  tree = ss.cKDTree(points)
  dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points]
  a = avgdigamma(zip2(x,z), dvec)
  b = avgdigamma(zip2(y,z), dvec)
  c = avgdigamma(z,dvec)
  d = digamma(k)
  return (-a-b+c+d) / log(base) 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def kldiv(x, xp, k=3, base=2):
  """KL Divergence between p and q for x~p(x),xp~q(x).
  x, xp should be a list of vectors, e.g. x = [[1.3],[3.7],[5.1],[2.4]]
  if x is a one-dimensional scalar and we have four samples
  """
  assert k <= len(x) - 1, 'Set k smaller than num samples - 1.'
  assert k <= len(xp) - 1, 'Set k smaller than num samples - 1.'
  assert len(x[0]) == len(xp[0]), 'Two distributions must have same dim.'
  d = len(x[0])
  n = len(x)
  m = len(xp)
  const = log(m) - log(n-1)
  tree = ss.cKDTree(x)
  treep = ss.cKDTree(xp)
  nn = [tree.query(point, k+1, p=float('inf'))[0][k] for point in x]
  nnp = [treep.query(point, k, p=float('inf'))[0][k-1] for point in x]
  return (const + d * np.mean(map(log, nnp)) \
    - d * np.mean(map(log, nn))) / log(base)

# DISCRETE ESTIMATORS 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def evaluate_pbc_fast(self, points):
        grid = points
        pos = self.pos
        box = self.box
        d = self.sigma * 2.5
        results = np.zeros(grid.shape[1], dtype=float)
        gridT = grid[::-1].T[:]
        tree = cKDTree(gridT, boxsize=box)
        # the indices of grid elements within distane d from each of the pos
        scale = 2. * self.sigma**2
        indlist = tree.query_ball_point(pos, d)
        for n, ind in enumerate(indlist):
            dr = gridT[ind, :] - pos[n]
            cond = np.where(dr > box / 2.)
            dr[cond] -= box[cond[1]]
            cond = np.where(dr < -box / 2.)
            dr[cond] += box[cond[1]]
            dens = np.exp(-np.sum(dr * dr, axis=1) / scale)
            results[ind] += dens

        return results 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def _create_mesh(self):
        """ Mesh assignment method

            Based on a target value, determine a mesh size for the testlines
            that is compatible with the simulation box.
            Create the grid and initialize a cKDTree object with it to
            facilitate fast searching of the gridpoints touched by molecules.
        """
        box = utilities.get_box(self.universe, self.normal)
        n, d = utilities.compute_compatible_mesh_params(self.target_mesh, box)
        self.mesh_nx = n[0]
        self.mesh_ny = n[1]
        self.mesh_dx = d[0]
        self.mesh_dy = d[1]

        _x = np.linspace(0, box[0], num=int(self.mesh_nx), endpoint=False)
        _y = np.linspace(0, box[1], num=int(self.mesh_ny), endpoint=False)
        _X, _Y = np.meshgrid(_x, _y)
        self.meshpoints = np.array([_X.ravel(), _Y.ravel()]).T
        self.meshtree = cKDTree(self.meshpoints, boxsize=box[:2]) 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def adi(R_est, t_est, R_gt, t_gt, pts):
  """Average Distance of Model Points for objects with indistinguishable views
  - by Hinterstoisser et al. (ACCV'12).

  :param R_est: 3x3 ndarray with the estimated rotation matrix.
  :param t_est: 3x1 ndarray with the estimated translation vector.
  :param R_gt: 3x3 ndarray with the ground-truth rotation matrix.
  :param t_gt: 3x1 ndarray with the ground-truth translation vector.
  :param pts: nx3 ndarray with 3D model points.
  :return: The calculated error.
  """
  pts_est = misc.transform_pts_Rt(pts, R_est, t_est)
  pts_gt = misc.transform_pts_Rt(pts, R_gt, t_gt)

  # Calculate distances to the nearest neighbors from vertices in the
  # ground-truth pose to vertices in the estimated pose.
  nn_index = spatial.cKDTree(pts_est)
  nn_dists, _ = nn_index.query(pts_gt, k=1)

  e = nn_dists.mean()
  return e 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def adi(R_est, t_est, R_gt, t_gt, model):
    """
    Average Distance of Model Points for objects with indistinguishable views
    - by Hinterstoisser et al. (ACCV 2012).

    :param R_est, t_est: Estimated pose (3x3 rot. matrix and 3x1 trans. vector).
    :param R_gt, t_gt: GT pose (3x3 rot. matrix and 3x1 trans. vector).
    :param model: Object model given by a dictionary where item 'pts'
    is nx3 ndarray with 3D model points.
    :return: Error of pose_est w.r.t. pose_gt.
    """
    pts_est = misc.transform_pts_Rt(model['pts'], R_est, t_est)
    pts_gt = misc.transform_pts_Rt(model['pts'], R_gt, t_gt)

    # Calculate distances to the nearest neighbors from pts_gt to pts_est
    nn_index = spatial.cKDTree(pts_est)
    nn_dists, _ = nn_index.query(pts_gt, k=1)

    e = nn_dists.mean()
    return e 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def bench_build(self):
        print()
        print('        Constructing kd-tree')
        print('=====================================')
        print(' dim | # points |  KDTree  | cKDTree ')

        for (m, n, repeat) in [(3,10000,3), (8,10000,3), (16,10000,3)]:
            print('%4s | %7s ' % (m, n), end=' ')
            sys.stdout.flush()

            data = np.concatenate((np.random.randn(n//2,m),
                                   np.random.randn(n-n//2,m)+np.ones(m)))

            print('| %6.3fs ' % (measure('T1 = KDTree(data)', repeat) / repeat), end=' ')
            sys.stdout.flush()
            print('| %6.3fs' % (measure('T2 = cKDTree(data)', repeat) / repeat), end=' ')
            sys.stdout.flush()
            print('') 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def _find_correspondence(surf, ref_surf, eps=0, n_jobs=1, use_cell=False):

    if use_cell:
        points = compute_cell_center(surf)
        ref_points = compute_cell_center(ref_surf)
    else:
        points = get_points(surf)
        ref_points = get_points(ref_surf)

    tree = cKDTree(ref_points, leafsize=20, compact_nodes=False,
                   copy_data=False, balanced_tree=False)
    d, idx = tree.query(points, k=1, eps=0, n_jobs=n_jobs,
                        distance_upper_bound=eps+np.finfo(np.float).eps)

    if np.isinf(d).any():
        raise ValueError('Cannot find correspondences. Try increasing '
                         'tolerance.')

    return idx 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def _view_kd(cls, data, data_view, target_view, x_tolerance=1e-3, y_tolerance=1e-3):
        """
        Appropriate if:
            - Grid is regular
            - Data is regular
            - Grid and data have same cellsize OR no interpolation is needed
        """
        nodata = target_view.nodata
        view = np.full(target_view.shape, nodata)
        viewrows, viewcols = target_view.grid_indices()
        rows, cols = data_view.grid_indices()
        ytree = spatial.cKDTree(rows[:, None])
        xtree = spatial.cKDTree(cols[:, None])
        ydist, y_ix = ytree.query(viewrows[:, None])
        xdist, x_ix = xtree.query(viewcols[:, None])
        y_passed = ydist < y_tolerance
        x_passed = xdist < x_tolerance
        view[np.ix_(y_passed, x_passed)] = data[y_ix[y_passed]][:, x_ix[x_passed]]
        return view 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def _view_kd_2d(cls, data, data_view, target_view, x_tolerance=1e-3, y_tolerance=1e-3):
        t_xmin, t_ymin, t_xmax, t_ymax = target_view.bbox
        d_xmin, d_ymin, d_xmax, d_ymax = data_view.bbox
        nodata = target_view.nodata
        view = np.full(target_view.shape, nodata)
        yx_tolerance = np.sqrt(x_tolerance**2 + y_tolerance**2)
        viewrows, viewcols = target_view.grid_indices()
        rows, cols = data_view.grid_indices()
        row_bool = (rows <= t_ymax + y_tolerance) & (rows >= t_ymin - y_tolerance)
        col_bool = (cols <= t_xmax + x_tolerance) & (cols >= t_xmin - x_tolerance)
        yx_tree = np.vstack(np.dstack(np.meshgrid(rows[row_bool], cols[col_bool], indexing='ij')))
        yx_query = np.vstack(np.dstack(np.meshgrid(viewrows, viewcols, indexing='ij')))
        tree = spatial.cKDTree(yx_tree)
        yx_dist, yx_ix = tree.query(yx_query)
        yx_passed = yx_dist < yx_tolerance
        view.flat[yx_passed] = data[np.ix_(row_bool, col_bool)].flat[yx_ix[yx_passed]]
        return view 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def adi(pose_est, pose_gt, model):
    """
    Average Distance of Model Points for objects with indistinguishable views
    - by Hinterstoisser et al. (ACCV 2012).

    :param pose_est: Estimated pose given by a dictionary:
    {'R': 3x3 rotation matrix, 't': 3x1 translation vector}.
    :param pose_gt: The ground truth pose given by a dictionary (as pose_est).
    :param model: Object model given by a dictionary where item 'pts'
    is nx3 ndarray with 3D model points.
    :return: Error of pose_est w.r.t. pose_gt.
    """
    pts_gt = misc.transform_pts_Rt(model['pts'], pose_gt['R'], pose_gt['t'])
    pts_est = misc.transform_pts_Rt(model['pts'], pose_est['R'], pose_est['t'])

    # Calculate distances to the nearest neighbors from pts_gt to pts_est
    nn_index = spatial.cKDTree(pts_est)
    nn_dists, _ = nn_index.query(pts_gt, k=1)

    e = nn_dists.mean()
    return e 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def test_ridges(self, name):
        # Check that the ridges computed by Voronoi indeed separate
        # the regions of nearest neighborhood, by comparing the result
        # to KDTree.

        points = DATASETS[name]

        tree = KDTree(points)
        vor = qhull.Voronoi(points)

        for p, v in vor.ridge_dict.items():
            # consider only finite ridges
            if not np.all(np.asarray(v) >= 0):
                continue

            ridge_midpoint = vor.vertices[v].mean(axis=0)
            d = 1e-6 * (points[p[0]] - ridge_midpoint)

            dist, k = tree.query(ridge_midpoint + d, k=1)
            assert_equal(k, p[0])

            dist, k = tree.query(ridge_midpoint - d, k=1)
            assert_equal(k, p[1]) 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def test_query_ball_point_multithreading():
    np.random.seed(0)
    n = 5000
    k = 2
    points = np.random.randn(n,k)
    T = cKDTree(points)
    l1 = T.query_ball_point(points,0.003,n_jobs=1)
    l2 = T.query_ball_point(points,0.003,n_jobs=64)
    l3 = T.query_ball_point(points,0.003,n_jobs=-1)
    
    for i in range(n):
        if l1[i] or l2[i]:
            assert_array_equal(l1[i],l2[i])
        
    for i in range(n):
        if l1[i] or l3[i]:
            assert_array_equal(l1[i],l3[i]) 

# 需要导入模块: from scipy import spatial [as 别名]
# 或者: from scipy.spatial import cKDTree [as 别名]
def setup_method(self):
        n = 50
        m = 4
        np.random.seed(0)
        data1 = np.random.randn(n,m)
        data2 = np.random.randn(n,m)
        self.T1 = cKDTree(data1,leafsize=2)
        self.T2 = cKDTree(data2,leafsize=2)
        self.ref_T1 = KDTree(data1, leafsize=2)
        self.ref_T2 = KDTree(data2, leafsize=2)
        self.r = 0.5
        self.n = n
        self.m = m
        self.data1 = data1
        self.data2 = data2
        self.p = 2