Python numpy.meshgrid() Examples

The following are 30 code examples of numpy.meshgrid(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_inside_volume(self, sphere3_msh):
        X, Y, Z = np.meshgrid(np.linspace(-100, 100, 100),
                              np.linspace(-40, 40, 10), [0])
        np.random.seed(0)
        points = np.vstack([X.reshape(-1), Y.reshape(-1), Z.reshape(-1)]).T
        points += np.random.random_sample(points.shape) - .5
        dist = np.linalg.norm(points, axis=1)

        msh = sphere3_msh.crop_mesh([4, 5])
        points_outside = points[(dist > 95) + (dist < 84)]
        inside = msh.test_inside_volume(points_outside)
        assert np.all(~inside)

        points_inside = points[(dist <= 94) * (dist >= 86)]
        inside = msh.test_inside_volume(points_inside)
        assert np.all(inside) 
Example #2
Source File: test_electrode_placement.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_draw_polygon_2D(self):
        X, Y = np.meshgrid(np.linspace(-8, 8, 25, dtype=float), np.linspace(-8, 8, 25, dtype=float))
        nodes = np.vstack([X.reshape(-1), Y.reshape(-1)]).T
        tri = scipy.spatial.Delaunay(nodes)
        poly = np.array([[5.5, 5.5], [5.5, -5.5], [-5.5, -5.5], [-5.5, 5.5]])
        #angles = np.linspace(0, 2*np.pi, endpoint=False, num=20)
        #poly = np.vstack([5*np.sin(angles), 5*np.cos(angles)]).T
        trs = tri.simplices[:, [1,0,2]]
        new_points, _ = electrode_placement._draw_polygon_2D(
            poly, tri.points, trs, ends=True)
        bar = np.mean(new_points[trs], axis=1)
        m = new_points[trs[:, 1:]] -\
            new_points[trs[:, 0]][:, None, :]
        area = .5 * -np.linalg.det(m)
        inside = electrode_placement._point_inside_polygon(poly, bar, tol=1e-3)
        #plt.triplot(new_points[:, 0], new_points[:, 1], tri.simplices.copy())
        #plt.show()
        assert np.isclose(np.sum(area[inside]), 11**2) 
Example #3
Source File: depth_utils.py    From DOTA_models with Apache License 2.0 6 votes vote down vote up
def get_point_cloud_from_z(Y, camera_matrix):
  """Projects the depth image Y into a 3D point cloud.
  Inputs:
    Y is ...xHxW
    camera_matrix
  Outputs:
    X is positive going right
    Y is positive into the image
    Z is positive up in the image
    XYZ is ...xHxWx3
  """
  x, z = np.meshgrid(np.arange(Y.shape[-1]),
                     np.arange(Y.shape[-2]-1, -1, -1))
  for i in range(Y.ndim-2):
    x = np.expand_dims(x, axis=0)
    z = np.expand_dims(z, axis=0)
  X = (x-camera_matrix.xc) * Y / camera_matrix.f
  Z = (z-camera_matrix.zc) * Y / camera_matrix.f
  XYZ = np.concatenate((X[...,np.newaxis], Y[...,np.newaxis],
                        Z[...,np.newaxis]), axis=X.ndim)
  return XYZ 
Example #4
Source File: test_electrode_placement.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_inside_complex_polygon(self):
        X, Y = np.meshgrid(np.linspace(-8, 8, 17, dtype=float), np.linspace(-8, 8, 17, dtype=float))
        nodes = np.vstack([X.reshape(-1), Y.reshape(-1)]).T
        tri = scipy.spatial.Delaunay(nodes)
        poly = np.array([[5, 5], [5, -5], [-5, -5], [-5, 5]])
        hole1 = np.array([[2, 2], [2, -2], [-2, -2], [-2, 2]])
        hole2 = np.array([[4, 4], [4, 3], [3, 3], [3, 4]])
        trs = tri.simplices[:, [1,0,2]]
        inside = electrode_placement._inside_complex_polygon(poly, tri.points,
                                                             trs,
                                                             holes=[hole1, hole2],
                                                             tol=1e-3)
        m = tri.points[trs[inside, 1:]] -\
            tri.points[trs[inside, 0]][:, None, :]
        area = .5 * -np.linalg.det(m)
        assert np.isclose(np.sum(area), 100-16-1) 
Example #5
Source File: test_PlanetPopulation.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_sma_radius(self):
        """
        Test that sma and radius values outside of the range have zero probability
        """
        
        for mod in self.allmods:
            if 'dist_sma_radius' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)
                
                a = np.logspace(np.log10(pp.arange[0].value/10.),np.log10(pp.arange[1].value*100),100)
                Rp = np.logspace(np.log10(pp.Rprange[0].value/10.),np.log10(pp.Rprange[1].value*100),100)
                
                aa, RR = np.meshgrid(a,Rp)
                
                fr = pp.dist_sma_radius(aa,RR)
                self.assertTrue(np.all(fr[aa < pp.arange[0].value] == 0),'dist_sma_radius low bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[aa > pp.arange[1].value] == 0),'dist_sma_radius high bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR < pp.Rprange[0].value] == 0),'dist_sma_radius low bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR > pp.Rprange[1].value] == 0),'dist_sma_radius high bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(aa > pp.arange[0].value) & (aa < pp.arange[1].value) & (RR > pp.Rprange[0].value) & (RR < pp.Rprange[1].value)] > 0),'dist_sma_radius is improper pdf for %s'%mod.__name__) 
Example #6
Source File: test_common.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 
Example #7
Source File: test_common.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 
Example #8
Source File: Stark.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def calcfbetaInput(self):
        # table 17 in Leinert et al. (1998)
        # Zodiacal Light brightness function of solar LON (rows) and LAT (columns)
        # values given in W m−2 sr−1 μm−1 for a wavelength of 500 nm
        path = os.path.split(inspect.getfile(self.__class__))[0]
        Izod = np.loadtxt(os.path.join(path, 'Leinert98_table17.txt'))*1e-8 # W/m2/sr/um
        # create data point coordinates
        lon_pts = np.array([0., 5, 10, 15, 20, 25, 30, 35, 40, 45, 60, 75, 90,
                105, 120, 135, 150, 165, 180]) # deg
        lat_pts = np.array([0., 5, 10, 15, 20, 25, 30, 45, 60, 75, 90]) # deg
        y_pts, x_pts = np.meshgrid(lat_pts, lon_pts)
        points = np.array(list(zip(np.concatenate(x_pts), np.concatenate(y_pts))))
        # create data values, normalized by (90,0) value
        z = Izod/Izod[12,0]
        values = z.reshape(z.size)
        return  points, values 
Example #9
Source File: images.py    From neuropythy with GNU Affero General Public License v3.0 6 votes vote down vote up
def image_reslice(image, spec, method=None, fill=0, dtype=None, weights=None, image_type=None):
    '''
    image_reslice(image, spec) yields a duplicate of the given image resliced to have the voxels
      indicated by the given image spec. Note that spec may be an image itself.

    Optional arguments that can be passed to image_interpolate() (asside from affine) are allowed
    here and are passed through.
    '''
    if image_type is None and is_image(image): image_type = to_image_type(image)
    spec = to_image_spec(spec)
    image = to_image(image)
    # we make a big mesh and interpolate at these points...
    imsh = spec['image_shape']
    (args, kw) = ([np.arange(n) for n in imsh[:3]], {'indexing': 'ij'})
    ijk = np.asarray([u.flatten() for u in np.meshgrid(*args, **kw)])
    ijk = np.dot(spec['affine'], np.vstack([ijk, np.ones([1,ijk.shape[1]])]))[:3]
    # interpolate here...
    u = image_interpolate(image, ijk, method=method, fill=fill, dtype=dtype, weights=weights)
    return to_image((np.reshape(u, imsh), spec), image_type=image_type) 
Example #10
Source File: test_common.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 
Example #11
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_find_tetrahedron_with_points_non_convex(self, sphere3_msh):
        X, Y, Z = np.meshgrid(np.linspace(-100, 100, 100),
                              np.linspace(-40, 40, 100), [0])
        points = np.vstack([X.reshape(-1), Y.reshape(-1), Z.reshape(-1)]).T
        dist = np.linalg.norm(points, axis=1)

        msh = sphere3_msh.crop_mesh(5)
        points_outside = points[(dist > 95) + (dist < 89)]
        th_with_points, bar = msh.find_tetrahedron_with_points(points_outside)
        assert np.all(th_with_points == -1)
        assert np.allclose(bar, 0)

        points_inside = points[(dist <= 94) * (dist >= 91)]
        th_with_points, bar = msh.find_tetrahedron_with_points(points_inside)
        eps = 1e-3
        assert np.all(th_with_points != -1)
        assert np.all(bar >= 0 - eps)
        assert np.all(bar <= 1. + eps)
        th_coords = \
            msh.nodes[msh.elm[th_with_points]]
        assert np.allclose(np.einsum('ikj, ik -> ij', th_coords, bar), points_inside) 
Example #12
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_to_nonlinear_grid(self, sphere3_msh):
        import nibabel
        data = sphere3_msh.elm.tag1
        f = mesh_io.ElementData(data, mesh=sphere3_msh)
        affine = np.array([[1, 0, 0, -100.5],
                           [0, 1, 0, -5],
                           [0, 0, 1, 0],
                           [0, 0, 0, 1]], dtype=float)
        x, y, z = np.meshgrid(np.arange(-100, 100),
                              np.arange(-5, 5),
                              np.arange(0, 1),
                              indexing='ij')
        nonl_transform = np.concatenate(
            (x[..., None], y[..., None], z[..., None]), axis=3).astype(float)
        img = nibabel.Nifti1Pair(nonl_transform, affine)
        tempf = tempfile.NamedTemporaryFile(suffix='.nii', delete=False)
        fn = tempf.name
        tempf.close()
        nibabel.save(img, fn)
        interp = f.to_deformed_grid(fn, fn, method='assign')
        interp = interp.get_data()
        assert np.isclose(interp[100, 5, 0], 3)
        assert np.isclose(interp[187, 5, 0], 4)
        assert np.isclose(interp[193, 5, 0], 5)
        assert np.isclose(interp[198, 5, 0], 0) 
Example #13
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_to_nonlinear_grid_crop(self, sphere3_msh):
        import nibabel
        data = sphere3_msh.elm.tag1
        f = mesh_io.ElementData(data, mesh=sphere3_msh)
        affine = np.array([[1, 0, 0, -100.5],
                           [0, 1, 0, -5],
                           [0, 0, 1, 0],
                           [0, 0, 0, 1]], dtype=float)
        x, y, z = np.meshgrid(np.arange(-100, 100),
                              np.arange(-5, 5),
                              np.arange(0, 1),
                              indexing='ij')
        nonl_transform = np.concatenate(
            (x[..., None], y[..., None], z[..., None]), axis=3).astype(float)
        img = nibabel.Nifti1Pair(nonl_transform, affine)
        tempf = tempfile.NamedTemporaryFile(suffix='.nii', delete=False)
        fn = tempf.name
        tempf.close()
        nibabel.save(img, fn)
        interp = f.to_deformed_grid(fn, fn, tags=3, method='assign')
        interp = interp.get_data()
        assert np.isclose(interp[100, 5, 0], 3)
        assert np.isclose(interp[187, 5, 0], 0)
        assert np.isclose(interp[193, 5, 0], 0)
        assert np.isclose(interp[198, 5, 0], 0) 
Example #14
Source File: test_topfarm.py    From TOPFARM with GNU Affero General Public License v3.0 6 votes vote down vote up
def test_2x3(self):
        # Loading the water depth map
        dat = loadtxt('data/WaterDepth1.dat')
        X, Y = meshgrid(linspace(0., 1000., 50), linspace(0., 1000., 50))
        depth = array(zip(X.flatten(), Y.flatten(), dat.flatten()))
        borders = array([[200, 200], [150, 500], [200, 800], [600, 900], [700, 700], [900, 500], [800, 200], [500, 100], [200, 200]])
        baseline = array([[587.5, 223.07692308], [525., 346.15384615], [837.5, 530.76923077], [525., 530.76923077], [525., 838.46153846], [837.5, 469.23076923]])

        wt_desc = WTDescFromWTG('data/V80-2MW-offshore.wtg').wt_desc
        wt_layout = GenericWindFarmTurbineLayout([WTPC(wt_desc=wt_desc, position=pos) for pos in baseline])

        t = Topfarm(
            baseline_layout = wt_layout,
            borders = borders,
            depth_map = depth,
            dist_WT_D = 5.0,
            distribution='spiral',
            wind_speeds=[4., 8., 20.],
            wind_directions=linspace(0., 360., 36)[:-1]
        )

        t.run()

        self.fail('make save function')
        t.save() 
Example #15
Source File: spectral_graph_partition.py    From LanczosNetwork with MIT License 6 votes vote down vote up
def get_L_cluster_cut(L, node_label):
  adj = L - np.diag(np.diag(L))
  adj[adj != 0] = 1.0
  num_nodes = adj.shape[0]
  idx_row, idx_col = np.meshgrid(range(num_nodes), range(num_nodes))
  idx_row, idx_col = idx_row.flatten().astype(
      np.int64), idx_col.flatten().astype(np.int64)
  mask = (node_label[idx_row] == node_label[idx_col]).reshape(
      num_nodes, num_nodes).astype(np.float)

  adj_cluster = adj * mask
  adj_cut = adj - adj_cluster
  L_cut = get_laplacian(adj_cut, graph_laplacian_type='L4')
  L_cluster = get_laplacian(adj_cluster, graph_laplacian_type='L4')

  return L_cluster, L_cut 
Example #16
Source File: region_proposal_network.py    From easy-faster-rcnn.pytorch with MIT License 6 votes vote down vote up
def generate_anchors(self, image_width: int, image_height: int, num_x_anchors: int, num_y_anchors: int) -> Tensor:
        center_ys = np.linspace(start=0, stop=image_height, num=num_y_anchors + 2)[1:-1]
        center_xs = np.linspace(start=0, stop=image_width, num=num_x_anchors + 2)[1:-1]
        ratios = np.array(self._anchor_ratios)
        ratios = ratios[:, 0] / ratios[:, 1]
        sizes = np.array(self._anchor_sizes)

        # NOTE: it's important to let `center_ys` be the major index (i.e., move horizontally and then vertically) for consistency with 2D convolution
        # giving the string 'ij' returns a meshgrid with matrix indexing, i.e., with shape (#center_ys, #center_xs, #ratios)
        center_ys, center_xs, ratios, sizes = np.meshgrid(center_ys, center_xs, ratios, sizes, indexing='ij')

        center_ys = center_ys.reshape(-1)
        center_xs = center_xs.reshape(-1)
        ratios = ratios.reshape(-1)
        sizes = sizes.reshape(-1)

        widths = sizes * np.sqrt(1 / ratios)
        heights = sizes * np.sqrt(ratios)

        center_based_anchor_bboxes = np.stack((center_xs, center_ys, widths, heights), axis=1)
        center_based_anchor_bboxes = torch.from_numpy(center_based_anchor_bboxes).float()
        anchor_bboxes = BBox.from_center_base(center_based_anchor_bboxes)

        return anchor_bboxes 
Example #17
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_from_data_grid_vec_extra_args(self, sphere3_msh):
        X, Y, Z = np.meshgrid(np.linspace(-100, 100, 201),
                              np.linspace(-100, 100, 201),
                              np.linspace(-1, 1, 3),
                              indexing='ij')
        V = np.stack([X, Y, Z], axis=3)
        affine = np.array([[1, 0, 0, -100],
                           [0, 1, 0, -100],
                           [0, 0, 1, -1],
                           [0, 0, 0, 1]], dtype=float)
        ed = mesh_io.ElementData.from_data_grid(sphere3_msh, V, affine,
                                             cval=0.0, order=1)
        bar = sphere3_msh.elements_baricenters().value
        in_area = (bar[:, 2] >= -1) * (bar[:, 2] <= 1)
        assert np.allclose(bar[in_area], ed.value[in_area])
        assert np.allclose(ed.value[~in_area], 0) 
Example #18
Source File: test_common.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 
Example #19
Source File: test_sq.py    From pyqmc with MIT License 6 votes vote down vote up
def test_big_cell():
    import time

    a = 1
    ncell = (2, 2, 2)
    Lvecs = np.diag(ncell) * a
    unit_cell = np.zeros((4, 3))
    unit_cell[1:] = (np.ones((3, 3)) - np.eye(3)) * a / 2

    grid = np.meshgrid(*map(np.arange, ncell), indexing="ij")
    shifts = np.stack(list(map(np.ravel, grid)), axis=1)
    supercell = (shifts[:, np.newaxis] + unit_cell[np.newaxis]).reshape(1, -1, 3)

    configs = supercell.repeat(1000, axis=0)
    configs += np.random.randn(*configs.shape) * 0.1

    df = run(Lvecs, configs, 8)
    df = df.groupby("qmag").mean().reset_index()

    large_q = df[-35:-10]["Sq"]
    mean = np.mean(large_q - 1)
    rms = np.sqrt(np.mean((large_q - 1) ** 2))
    assert np.abs(mean) < 0.01, mean
    assert rms < 0.1, rms 
Example #20
Source File: snippets.py    From Collaborative-Learning-for-Weakly-Supervised-Object-Detection with MIT License 6 votes vote down vote up
def generate_anchors_pre(height, width, feat_stride, anchor_scales=(8,16,32), anchor_ratios=(0.5,1,2)):
  """ A wrapper function to generate anchors given different scales
    Also return the number of anchors in variable 'length'
  """
  anchors = generate_anchors(ratios=np.array(anchor_ratios), scales=np.array(anchor_scales))
  A = anchors.shape[0]
  shift_x = np.arange(0, width) * feat_stride
  shift_y = np.arange(0, height) * feat_stride
  shift_x, shift_y = np.meshgrid(shift_x, shift_y)
  shifts = np.vstack((shift_x.ravel(), shift_y.ravel(), shift_x.ravel(), shift_y.ravel())).transpose()
  K = shifts.shape[0]
  # width changes faster, so here it is H, W, C
  anchors = anchors.reshape((1, A, 4)) + shifts.reshape((1, K, 4)).transpose((1, 0, 2))
  anchors = anchors.reshape((K * A, 4)).astype(np.float32, copy=False)
  length = np.int32(anchors.shape[0])

  return anchors, length 
Example #21
Source File: pc_util.py    From H3DNet with MIT License 5 votes vote down vote up
def gaussian_3d(x_mean, y_mean, z_mean, vxy, vz, dev=2.0):
  x, y, z = np.meshgrid(np.arange(vxy), np.arange(vxy), np.arange(vz))
  #z=(1.0/(2.0*np.pi*dev*dev))*np.exp(-((x-x_mean)**2+ (y-y_mean)**2)/(2.0*dev**2))
  m=np.exp(-((x-x_mean)**2 + (y-y_mean)**2+(z-z_mean)**2)/(2.0*dev**2))
  return m 
Example #22
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def test_to_nonlinear_grid_nodedata(self, sphere3_msh):
        import nibabel
        data = sphere3_msh.nodes.node_coord[:, 0]
        f = mesh_io.NodeData(data, mesh=sphere3_msh)
        affine = np.array([[1, 0, 0, -100.5],
                           [0, 1, 0, -5],
                           [0, 0, 1, 0],
                           [0, 0, 0, 1]], dtype=float)
        x, y, z = np.meshgrid(np.arange(-100, 100),
                              np.arange(-5, 5),
                              np.arange(0, 1),
                              indexing='ij')
        nonl_transform = np.concatenate(
            (x[..., None], y[..., None], z[..., None]), axis=3).astype(float)
        img = nibabel.Nifti1Pair(nonl_transform, affine)
        tempf = tempfile.NamedTemporaryFile(suffix='.nii', delete=False)
        fn = tempf.name
        tempf.close()
        nibabel.save(img, fn)
        interp = f.to_deformed_grid(fn, fn, tags=3, order=1)
        interp = interp.get_data()
        '''
        import matplotlib.pyplot as plt
        plt.figure()
        plt.imshow(np.squeeze(interp))
        plt.colorbar()
        plt.show()
        '''
        assert np.isclose(interp[150, 5, 0], 50)
        assert np.isclose(interp[190, 5, 0], 0)
        assert np.isclose(interp[193, 5, 0], 0)
        assert np.isclose(interp[198, 5, 0], 0) 
Example #23
Source File: test_transformations.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def nonl_large_rotate():
    x, y, z = np.meshgrid(np.linspace(-120, 120, 10),
                          np.linspace(-120, 120, 10),
                          np.linspace(-120, 120, 10),
                          indexing='ij')
    nonl_transform = np.concatenate(
        (-y[..., None], x[..., None], z[..., None]), axis=3).astype(float)
    s = 240 / 9.
    affine = np.array([[s, 0, 0, -120],
                       [0, s, 0, -120],
                       [0, 0, s, -120],
                       [0, 0, 0, 1]], dtype=float)
    return nonl_transform, affine 
Example #24
Source File: test_mesh_io.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def test_to_nonlinear_grid_linear_interp(self, sphere3_msh):
        import nibabel
        data = sphere3_msh.elements_baricenters().value[:, 0]
        f = mesh_io.ElementData(data, mesh=sphere3_msh)
        affine = np.array([[1, 0, 0, -100.5],
                           [0, 1, 0, -5],
                           [0, 0, 1, 0],
                           [0, 0, 0, 1]], dtype=float)
        x, y, z = np.meshgrid(np.arange(-100, 100),
                              np.arange(-5, 5),
                              np.arange(0, 1),
                              indexing='ij')
        nonl_transform = np.concatenate(
            (x[..., None], y[..., None], z[..., None]), axis=3).astype(float)
        img = nibabel.Nifti1Pair(nonl_transform, affine)
        tempf = tempfile.NamedTemporaryFile(suffix='.nii',delete=False)
        fn = tempf.name
        tempf.close()
        nibabel.save(img, fn)
        interp = f.to_deformed_grid(fn, fn, order=1,
                                    method='linear', continuous=True)
        interp = interp.get_data()
        '''
        import matplotlib.pyplot as plt
        plt.figure()
        plt.imshow(np.squeeze(interp))
        plt.colorbar()
        plt.show()
        '''
        assert np.isclose(interp[150, 5, 0], 50, atol=1e-2)
        assert np.isclose(interp[190, 5, 0], 90, atol=1e-1)
        assert np.isclose(interp[191, 5, 0], 91, atol=5e-1)
        assert np.isclose(interp[198, 5, 0], 0) 
Example #25
Source File: omnitest.py    From pymoo with Apache License 2.0 5 votes vote down vote up
def _calc_pareto_set(self, n_pareto_points=500):
        # The Omni-test problem has 3^D Pareto subsets
        num_ps = int(3 ** self.n_var)
        h = int(n_pareto_points / num_ps)
        PS = np.zeros((num_ps * h, self.n_var))

        candidates = np.array([np.linspace(2 * m + 1, 2 * m + 3 / 2, h) for m in range(3)])
        # generate combination indices
        candidates_indices = [[0, 1, 2] for _ in range(self.n_var)]
        a = np.meshgrid(*candidates_indices)
        combination_indices = np.array(a).T.reshape(-1, self.n_var)
        # generate 3^D combinations
        for i in range(num_ps):
            PS[i * h:i * h + h, :] = candidates[combination_indices[i]].T
        return PS 
Example #26
Source File: post_proc.py    From HorizonNet with MIT License 5 votes vote down vote up
def fuv2img(fuv, coorW=1024, floorW=1024, floorH=512):
    '''
    Project 1d signal in uv space to 2d floor plane image
    '''
    floor_plane_x, floor_plane_y = np.meshgrid(range(floorW), range(floorH))
    floor_plane_x, floor_plane_y = -(floor_plane_y - floorH / 2), floor_plane_x - floorW / 2
    floor_plane_coridx = (np.arctan2(floor_plane_y, floor_plane_x) / (2 * PI) + 0.5) * coorW - 0.5
    floor_plane = map_coordinates(fuv, floor_plane_coridx.reshape(1, -1), order=1, mode='wrap')
    floor_plane = floor_plane.reshape(floorH, floorW)
    return floor_plane 
Example #27
Source File: ada_lanczos_net.py    From LanczosNetwork with MIT License 5 votes vote down vote up
def _get_graph_laplacian(self, node_feat, adj_mask):
    """ Compute graph Laplacian

      Args:
        node_feat: float tensor, shape B X N X D
        adj_mask: float tensor, shape B X N X N, binary mask, should contain self-loop

      Returns:
        L: float tensor, shape B X N X N
    """
    batch_size = node_feat.shape[0]
    num_node = node_feat.shape[1]
    dim_feat = node_feat.shape[2]

    # compute pairwise distance
    idx_row, idx_col = np.meshgrid(range(num_node), range(num_node))
    idx_row, idx_col = torch.Tensor(idx_row.reshape(-1)).long().to(node_feat.device), torch.Tensor(
        idx_col.reshape(-1)).long().to(node_feat.device)

    diff = node_feat[:, idx_row, :] - node_feat[:, idx_col, :]  # shape B X N^2 X D
    dist2 = (diff * diff).sum(dim=2)  # shape B X N^2
    
    # sigma2, _ = torch.median(dist2, dim=1, keepdim=True) # median is sometimes 0
    # sigma2 = sigma2 + 1.0e-7

    sigma2 = torch.mean(dist2, dim=1, keepdim=True)

    A = torch.exp(-dist2 / sigma2)  # shape B X N^2
    A = A.reshape(batch_size, num_node, num_node) * adj_mask  # shape B X N X N
    row_sum = torch.sum(A, dim=2, keepdim=True)
    pad_row_sum = torch.zeros_like(row_sum)
    pad_row_sum[row_sum == 0.0] = 1.0    
    alpha = 0.5
    D = 1.0 / (row_sum + pad_row_sum).pow(alpha)  # shape B X N X 1
    L = D * A * D.transpose(1, 2)  # shape B X N X N

    return L 
Example #28
Source File: interpolation.py    From scarlet with MIT License 5 votes vote down vote up
def get_filter_coords(filter_values, center=None):
    """Create filter coordinate grid needed for the apply filter function

    Parameters
    ----------
    filter_values: array
        The 2D array of the filter to apply.
    center: tuple
        The center (y,x) of the filter. If `center` is `None` then
        `filter_values` must have an odd number of rows and columns
        and the center will be set to the center of `filter_values`.

    Returns
    -------
    coords: array
        The coordinates of the pixels in `filter_values`,
        where the coordinates of the `center` pixel are `(0,0)`.
    """
    if len(filter_values.shape) != 2:
        raise ValueError("`filter_values` must be 2D")
    if center is None:
        if filter_values.shape[0] % 2 == 0 or filter_values.shape[1] % 2 == 0:
            msg = """Ambiguous center of the `filter_values` array,
                     you must use a `filter_values` array
                     with an odd number of rows and columns or
                     calculate `coords` on your own."""
            raise ValueError(msg)
        center = [filter_values.shape[0]//2, filter_values.shape[1]//2]
    x = np.arange(filter_values.shape[1])
    y = np.arange(filter_values.shape[0])
    x, y = np.meshgrid(x, y)
    x -= center[1]
    y -= center[0]
    coords = np.dstack([y, x])
    return coords 
Example #29
Source File: operator.py    From scarlet with MIT License 5 votes vote down vote up
def sort_by_radius(shape, center=None):
    """Sort indices distance from the center

    Given a shape, calculate the distance of each
    pixel from the center and return the indices
    of each pixel, sorted by radial distance from
    the center, which need not be in the center
    of the image.

    Parameters
    ----------
    shape: `tuple`
        Shape (y,x) of the source frame.

    center: array-like
        Location of the center pixel.

    Returns
    -------
    didx: `~numpy.array`
        Indices of elements in an image with shape `shape`,
        sorted by distance from the center.
    """
    # Get the center pixels
    if center is None:
        cx = (shape[1] - 1) >> 1
        cy = (shape[0] - 1) >> 1
    else:
        cy, cx = int(center[0]), int(center[1])
    # Calculate the distance between each pixel and the peak
    x = np.arange(shape[1])
    y = np.arange(shape[0])
    X, Y = np.meshgrid(x, y)
    X = X - cx
    Y = Y - cy
    distance = np.sqrt(X ** 2 + Y ** 2)
    # Get the indices of the pixels sorted by distance from the peak
    didx = np.argsort(distance.flatten())
    return didx 
Example #30
Source File: ops_test.py    From object_detector_app with MIT License 5 votes vote down vote up
def test_meshgrid_numpy_comparison(self):
    """Tests meshgrid op with vectors, for which it should match numpy."""
    x = np.arange(4)
    y = np.arange(6)
    exp_xgrid, exp_ygrid = np.meshgrid(x, y)
    xgrid, ygrid = ops.meshgrid(x, y)
    with self.test_session() as sess:
      xgrid_output, ygrid_output = sess.run([xgrid, ygrid])
      self.assertAllEqual(xgrid_output, exp_xgrid)
      self.assertAllEqual(ygrid_output, exp_ygrid)