Python numpy.mgrid() Examples

def _tf_fspecial_gauss(size, sigma):
    """Function to mimic the 'fspecial' gaussian MATLAB function
    """
    x_data, y_data = np.mgrid[-size//2 + 1:size//2 + 1, -size//2 + 1:size//2 + 1]

    x_data = np.expand_dims(x_data, axis=-1)
    x_data = np.expand_dims(x_data, axis=-1)

    y_data = np.expand_dims(y_data, axis=-1)
    y_data = np.expand_dims(y_data, axis=-1)

    x = tf.constant(x_data, dtype=tf.float32)
    y = tf.constant(y_data, dtype=tf.float32)

    g = tf.exp(-((x**2 + y**2)/(2.0*sigma**2)))
    return g / tf.reduce_sum(g) 

def test_order_zero_2d_fails_on_extrapolation():
    a = np.mgrid[0:5, 0:5][0].reshape(25)
    b = np.mgrid[0:5, 0:5][1].reshape(25)
    df = pd.DataFrame({'a': a + 0.5, 'b': b + 0.5, 'c': b*3, 'garbage': ['test']*len(a)})
    df = make_bin_edges(df, 'a')
    df = make_bin_edges(df, 'b')
    df = df.sample(frac=1)  # Shuffle table to assure interpolation works given unsorted input

    i = Interpolation(df, ('garbage',), [('a', 'a_left', 'a_right'), ('b', 'b_left', 'b_right')],
                      order=0, extrapolate=False)

    column = np.arange(4, step=0.011)
    query = pd.DataFrame({'a': column, 'b': column, 'garbage': ['test']*(len(column))})

    with pytest.raises(ValueError) as error:
        i(query)

    message = error.value.args[0]

    assert 'Extrapolation' in message and 'a' in message 

def __init__(self, anchors, shape, sigma=0.5, randrange=None):
        """
        Args:
            anchors (list): list of center coordinates in range [0,1].
            shape(list or tuple): image shape in [h, w].
            sigma (float): sigma for Gaussian weight
            randrange (int): offset range. Defaults to shape[0] / 8
        """
        logger.warn("GaussianDeform is slow. Consider using it with 4 or more prefetching processes.")
        super(GaussianDeform, self).__init__()
        self.anchors = anchors
        self.K = len(self.anchors)
        self.shape = shape
        self.grid = np.mgrid[0:self.shape[0], 0:self.shape[1]].transpose(1, 2, 0)
        self.grid = self.grid.astype('float32')  # HxWx2

        gm = GaussianMap(self.shape, sigma=sigma)
        self.gws = np.array([gm.get_gaussian_weight(ank)
                             for ank in self.anchors], dtype='float32')  # KxHxW
        self.gws = self.gws.transpose(1, 2, 0)  # HxWxK
        if randrange is None:
            self.randrange = self.shape[0] / 8
        else:
            self.randrange = randrange
        self.sigma = sigma 

def same_transform(taff, ornt, shape):
    # Applying transformations implied by `ornt` to a made-up array
    # ``arr`` of shape `shape`, results in ``t_arr``. When the point
    # indices from ``arr`` are transformed by (the inverse of) `taff`,
    # and we index into ``t_arr`` with these transformed points, then we
    # should get the same values as we would from indexing into arr with
    # the untransformed points. 
    shape = np.array(shape)
    size = np.prod(shape)
    arr = np.arange(size).reshape(shape)
    # apply ornt transformations
    t_arr = apply_orientation(arr, ornt)
    # get all point indices in arr
    i,j,k = shape
    arr_pts = np.mgrid[:i,:j,:k].reshape((3,-1))
    # inverse of taff takes us from point index in arr to point index in
    # t_arr
    itaff = np.linalg.inv(taff)
    # apply itaff so that points indexed in t_arr should correspond 
    o2t_pts = np.dot(itaff[:3,:3], arr_pts) + itaff[:3,3][:,None]
    assert np.allclose(np.round(o2t_pts), o2t_pts)
    # fancy index out the t_arr values
    vals = t_arr[list(o2t_pts.astype('i'))]
    return np.all(vals == arr.ravel()) 

def _compute_q_vectors(self, box):
        """ Compute the q-vectors compatible with the current box dimensions.

            Calculated quantities:
            q_vectors : two 2D arrays forming the grid of q-values, similar
                        to a meshgrid
            Qxy       : array of the different q-vectors
            Q         : squared module of Qxy with the first element missing
                        (no Q = 0.0)
        """
        self.box = np.roll(box, 2 - self.normal)
        nmax = list(map(int, np.ceil(self.box[0:2] / self.alpha)))
        self.q_vectors = np.mgrid[0:nmax[0], 0:nmax[1]] * 1.0
        self.q_vectors[0] *= 2. * np.pi / box[0]
        self.q_vectors[1] *= 2. * np.pi / box[1]
        self.modes_shape = self.q_vectors[0].shape
        qx = self.q_vectors[0][:, 0]
        qy = self.q_vectors[1][0]
        Qx = np.repeat(qx, len(qy))
        Qy = np.tile(qy, len(qx))
        self.Qxy = np.vstack((Qx, Qy)).T
        self.Q = np.sqrt(np.sum(self.Qxy * self.Qxy, axis=1)[1:]) 

def plot(self, smearing=0.1):
        renderer = Renderer()
        N = len(self.structure)
        V = self.structure.get_volume()
        xs = np.mgrid[0.5:self.limit:1000j]
        dr = xs[1] - xs[0]
        norms = [ ( (x+dr/2)**3 - (x-dr/2)**3) for x in xs ]
        for pair in self.pairs:
            e1, e2 = pair
            vals = np.zeros(xs.shape)
            nanb = self.structure.comp[e1]*self.structure.comp[e2]
            prefactor = 1.0/(smearing*np.sqrt(2*np.pi))
            #prefactor = V**2 * (N-1)/N
            for w, d in zip(self.weights[pair], self.distances[pair]):
                if not d:
                    continue
                vals += np.exp(-(d-xs)**2/(2*smearing**2))*w
            vals = prefactor*vals/norms
            vals = [ v if v > 1e-4 else 0.0 for v in vals ]
            line = Line(zip(xs, vals), label='%s-%s' % (e1, e2))
            renderer.add(line)

        renderer.xaxis.label = 'interatomic distance [Å]'
        return renderer 

def quality(func, mesh, interpolator='nn', n=33):
    """Compute a quality factor (the quantity r**2 from TOMS792).

    interpolator must be in ('linear', 'nn').
    """
    fz = func(mesh.x, mesh.y)
    tri = Triangulation(mesh.x, mesh.y)
    intp = getattr(tri,
                   interpolator + '_extrapolator')(fz, bbox=(0., 1., 0., 1.))
    Y, X = np.mgrid[0:1:complex(0, n), 0:1:complex(0, n)]
    Z = func(X, Y)
    iz = intp[0:1:complex(0, n), 0:1:complex(0, n)]
    #nans = np.isnan(iz)
    #numgood = n*n - np.sum(np.array(nans.flat, np.int32))
    numgood = n * n

    SE = (Z - iz) ** 2
    SSE = np.sum(SE.flat)
    meanZ = np.sum(Z.flat) / numgood
    SM = (Z - meanZ) ** 2
    SSM = np.sum(SM.flat)

    r2 = 1.0 - SSE / SSM
    print(func.func_name, r2, SSE, SSM, numgood)
    return r2 

def random_warp( image ):
    assert image.shape == (256,256,3)
    range_ = numpy.linspace( 128-80, 128+80, 5 )
    mapx = numpy.broadcast_to( range_, (5,5) )
    mapy = mapx.T

    mapx = mapx + numpy.random.normal( size=(5,5), scale=5 )
    mapy = mapy + numpy.random.normal( size=(5,5), scale=5 )

    interp_mapx = cv2.resize( mapx, (80,80) )[8:72,8:72].astype('float32')
    interp_mapy = cv2.resize( mapy, (80,80) )[8:72,8:72].astype('float32')

    warped_image = cv2.remap( image, interp_mapx, interp_mapy, cv2.INTER_LINEAR )

    src_points = numpy.stack( [ mapx.ravel(), mapy.ravel() ], axis=-1 )
    dst_points = numpy.mgrid[0:65:16,0:65:16].T.reshape(-1,2)
    mat = umeyama( src_points, dst_points, True )[0:2]

    target_image = cv2.warpAffine( image, mat, (64,64) )

    return warped_image, target_image 

def __init__(self, anchors, shape, sigma=0.5, randrange=None):
        """
        :param anchors: in [0,1] coordinate
        :param shape: image shape in [h, w]
        :param sigma: sigma for Gaussian weight
        :param randrange: default to shape[0] / 8
        """
        logger.warn("GaussianDeform is slow. Consider using it with 4 or more prefetching processes.")
        super(GaussianDeform, self).__init__()
        self.anchors = anchors
        self.K = len(self.anchors)
        self.shape = shape
        self.grid = np.mgrid[0:self.shape[0], 0:self.shape[1]].transpose(1,2,0)
        self.grid = self.grid.astype('float32') # HxWx2

        gm = GaussianMap(self.shape, sigma=sigma)
        self.gws = np.array([gm.get_gaussian_weight(ank)
                             for ank in self.anchors], dtype='float32') # KxHxW
        self.gws = self.gws.transpose(1, 2, 0)  #HxWxK
        if randrange is None:
            self.randrange = self.shape[0] / 8
        else:
            self.randrange = randrange 

def _():
        """ additional tests

        here we generate a paraboloid (x^2+y^2) and a hyperbolic paraboloid
        (x^2-y^2) to check that the curvature code gives the right answers for
        the Gaussian (4, -4) and mean (2, 0) curvatures

        >>> import pytim
        >>> x,y=np.mgrid[-5:5,-5:5.]/2.
        >>> p = np.asarray(list(zip(x.flatten(),y.flatten())))
        >>> z1 = p[:,0]**2+p[:,1]**2
        >>> z2 = p[:,0]**2-p[:,1]**2
        >>>
        >>> for z in [z1, z2]:
        ...     pp = np.asarray(list(zip(x.flatten()+5,y.flatten()+5,z)))
        ...     curv = pytim.observables.Curvature(cutoff=1.,warning=False).compute(pp)
        ...     val =  (curv[np.logical_and(p[:,0]==0,p[:,1]==0)])
        ...     # add and subtract 1e3 to be sure to have -0 -> 0
        ...     print(np.array_str((val+1e3)-1e3, precision=2, suppress_small=True))
        [[4. 2.]]
        [[-4.  0.]]


        """
# 

def hyperball(ndim, radius):
    """Return a binary morphological filter containing pixels within `radius`.

    Parameters
    ----------
    ndim : int
        The number of dimensions of the filter.
    radius : int
        The radius of the filter.

    Returns
    -------
    ball : array of bool, shape [2 * radius + 1,] * ndim
        The required structural element
    """
    size = 2 * radius + 1
    center = [(radius,) * ndim]

    coords = np.mgrid[[slice(None, size),] * ndim].reshape(ndim, -1).T
    distances = np.ravel(spatial.distance_matrix(coords, center))
    selector = distances <= radius

    ball = np.zeros((size,) * ndim, dtype=bool)
    ball.ravel()[selector] = True
    return ball 

def __init__(self, pos, size):
        """
        pos is (...,3) array of the bar positions (the corner of each bar)
        size is (...,3) array of the sizes of each bar
        """
        nCubes = reduce(lambda a,b: a*b, pos.shape[:-1])
        cubeVerts = np.mgrid[0:2,0:2,0:2].reshape(3,8).transpose().reshape(1,8,3)
        cubeFaces = np.array([
            [0,1,2], [3,2,1],
            [4,5,6], [7,6,5],
            [0,1,4], [5,4,1],
            [2,3,6], [7,6,3],
            [0,2,4], [6,4,2],
            [1,3,5], [7,5,3]]).reshape(1,12,3)
        size = size.reshape((nCubes, 1, 3))
        pos = pos.reshape((nCubes, 1, 3))
        verts = cubeVerts * size + pos
        faces = cubeFaces + (np.arange(nCubes) * 8).reshape(nCubes,1,1)
        md = MeshData(verts.reshape(nCubes*8,3), faces.reshape(nCubes*12,3))
        
        GLMeshItem.__init__(self, meshdata=md, shader='shaded', smooth=False) 

def quality(func, mesh, interpolator='nn', n=33):
    """Compute a quality factor (the quantity r**2 from TOMS792).

    interpolator must be in ('linear', 'nn').
    """
    fz = func(mesh.x, mesh.y)
    tri = Triangulation(mesh.x, mesh.y)
    intp = getattr(tri,
                   interpolator + '_extrapolator')(fz, bbox=(0., 1., 0., 1.))
    Y, X = np.mgrid[0:1:complex(0, n), 0:1:complex(0, n)]
    Z = func(X, Y)
    iz = intp[0:1:complex(0, n), 0:1:complex(0, n)]
    #nans = np.isnan(iz)
    #numgood = n*n - np.sum(np.array(nans.flat, np.int32))
    numgood = n * n

    SE = (Z - iz) ** 2
    SSE = np.sum(SE.flat)
    meanZ = np.sum(Z.flat) / numgood
    SM = (Z - meanZ) ** 2
    SSM = np.sum(SM.flat)

    r2 = 1.0 - SSE / SSM
    print(func.func_name, r2, SSE, SSM, numgood)
    return r2 

def _generate_anchors_one_layer(h_I, w_I, h_l, w_l):
    """
    generate anchors on on layer
    return a ndarray with shape (h_l, w_l, 4), and the last dimmension in the order:[cx, cy, w, h]
    """
    y, x = np.mgrid[0: h_l, 0:w_l]
    cy = (y + config.anchor_offset) / h_l * h_I
    cx = (x + config.anchor_offset) / w_l * w_I
    
    anchor_scale = _get_scale(w_I, w_l)
    anchor_w = np.ones_like(cx) * anchor_scale
    anchor_h = np.ones_like(cx) * anchor_scale # cx.shape == cy.shape
    
    anchors = np.asarray([cx, cy, anchor_w, anchor_h])
    anchors = np.transpose(anchors, (1, 2, 0))
    
    return anchors 

def test_2d_interpolation():
    a = np.mgrid[0:5, 0:5][0].reshape(25)
    b = np.mgrid[0:5, 0:5][1].reshape(25)
    df = pd.DataFrame({'a': a, 'b': b, 'c': b, 'd': a})
    df = df.sample(frac=1)  # Shuffle table to assure interpolation works given unsorted input

    i = Interpolation(df, (), ('a', 'b'), 1, True)

    query = pd.DataFrame({'a': np.arange(4, step=0.01), 'b': np.arange(4, step=0.01)})

    assert np.allclose(query.b, i(query).c)
    assert np.allclose(query.a, i(query).d) 

def __getitem__(self, idx):
        intrinsics, _, _, _ = util.parse_intrinsics(os.path.join(self.instance_dir, "intrinsics.txt"),
                                                                  trgt_sidelength=self.img_sidelength)
        intrinsics = torch.Tensor(intrinsics).float()

        rgb = data_util.load_rgb(self.color_paths[idx], sidelength=self.img_sidelength)
        rgb = rgb.reshape(3, -1).transpose(1, 0)

        pose = data_util.load_pose(self.pose_paths[idx])

        if self.has_params:
            params = data_util.load_params(self.param_paths[idx])
        else:
            params = np.array([0])

        uv = np.mgrid[0:self.img_sidelength, 0:self.img_sidelength].astype(np.int32)
        uv = torch.from_numpy(np.flip(uv, axis=0).copy()).long()
        uv = uv.reshape(2, -1).transpose(1, 0)

        sample = {
            "instance_idx": torch.Tensor([self.instance_idx]).squeeze(),
            "rgb": torch.from_numpy(rgb).float(),
            "pose": torch.from_numpy(pose).float(),
            "uv": uv,
            "param": torch.from_numpy(params).float(),
            "intrinsics": intrinsics
        }
        return sample 

def test_ball_point_ints():
    """Regression test for #1373."""
    x, y = np.mgrid[0:4, 0:4]
    points = list(zip(x.ravel(), y.ravel()))
    tree = KDTree(points)
    assert_equal(sorted([4, 8, 9, 12]),
                 sorted(tree.query_ball_point((2, 0), 1)))
    points = np.asarray(points, dtype=np.float)
    tree = KDTree(points)
    assert_equal(sorted([4, 8, 9, 12]),
                 sorted(tree.query_ball_point((2, 0), 1)))

# cKDTree is specialized to type double points, so no need to make
# a unit test corresponding to test_ball_point_ints() 

def work(self, fig=None, ax=None):
        """Draw a two dimensional kernel density plot.
        You can specify either a figure or an axis to draw on.

        Parameters:
        -----------
        fig: matplotlib figure object
        ax: matplotlib axis object to draw on

        Returns:
        --------
        fig, ax: matplotlib figure and axis objects
        """
        if ax is None:
            if fig is None:
                return fig, ax
            else:
                ax = fig.gca()
        x = self.data[self.aes['x']]
        y = self.data[self.aes['y']]
        rvs = np.array([x, y])
        x_min = x.min()
        x_max = x.max()
        y_min = y.min()
        y_max = y.max()
        X, Y = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]
        positions = np.vstack([X.ravel(), Y.ravel()])
        values = np.vstack([x, y])
        import scipy.stats as stats
        kernel = stats.gaussian_kde(values)
        Z = np.reshape(kernel(positions).T, X.shape)
        ax.contour(Z, extent=[x_min, x_max, y_min, y_max])
        return fig, ax 

def genGTmap(h, w, pos_x, pos_y, sigma_h=1, sigma_w=1, init=None):
    """
    Compute the heat-map of size (w x h) with a gaussian distribution fit in
    position (pos_x, pos_y) and a convariance matix defined by the related
    sigma values.
    The resulting heat-map can be summed to a given heat-map init.
    """
    if pos_x>0 and pos_y>0:
        init = init if init is not None else []

        cov_matrix = np.eye(2) * ([sigma_h**2, sigma_w**2])

        x, y = np.mgrid[0:h, 0:w]
        pos = np.dstack((x, y))
        rv = multivariate_normal([pos_x, pos_y], cov_matrix)

        tmp = rv.pdf(pos)
        hmap = np.multiply(
            tmp, np.sqrt(np.power(2 * np.pi, 2) * np.linalg.det(cov_matrix))
        )
        idx = np.where(hmap.flatten() <= np.exp(-4.6052))
        hmap.flatten()[idx] = 0

        if np.size(init) == 0:
            return hmap

        assert (np.shape(init) == hmap.shape)
        hmap += init
        idx = np.where(hmap.flatten() > 1)
        hmap.flatten()[idx] = 1
        return hmap
    else:
        return np.zeros((h, w)) 

def generate_deformation_field_gauss(shape, points, max_deform=DEFORMATION_MAX,
                                     deform_smooth=DEFORMATION_SMOOTH):
    """ generate deformation field as combination of positive and
    negative Galatians densities scaled in range +/- max_deform

    :param tuple(int,int) shape: tuple of size 2
    :param points: <nb_points, 2> list of landmarks
    :param float max_deform: maximal deformation distance in any direction
    :param float deform_smooth: smoothing the deformation by Gaussian filter
    :return: np.array<shape>
    """
    ndim = len(shape)
    x, y = np.mgrid[0:shape[0], 0:shape[1]]
    pos_grid = np.rollaxis(np.array([x, y]), 0, 3)
    # initialise the deformation
    deform = np.zeros(shape)
    for point in points:
        sign = np.random.choice([-1, 1])
        cov = np.random.random((ndim, ndim))
        cov[np.eye(ndim, dtype=bool)] = 100 * np.random.random(ndim)
        # obtain a positive semi-definite matrix
        cov = np.dot(cov, cov.T) * (0.1 * np.mean(shape))
        gauss = stats.multivariate_normal(point, cov)
        deform += sign * gauss.pdf(pos_grid)
    # normalise the deformation and multiply by the amplitude
    deform *= max_deform / np.abs(deform).max()
    # set boundary region to zeros
    fix_deform_bounds = DEFORMATION_BOUNDARY_COEF * deform_smooth
    deform[:fix_deform_bounds, :] = 0
    deform[-fix_deform_bounds:, :] = 0
    deform[:, :fix_deform_bounds] = 0
    deform[:, -fix_deform_bounds:] = 0
    # smooth the deformation field
    deform = ndimage.gaussian_filter(deform, sigma=deform_smooth, order=0)
    return deform 

def test_mgrid_single_element(self, level=rlevel):
        """Ticket #339"""
        assert_array_equal(np.mgrid[0:0:1j], [0])
        assert_array_equal(np.mgrid[0:0], []) 

def cluster_kmeans(sim):
    """Run k-means on similarity matrix and segment"""
    sim_dim = sim.shape[0]
    sim = sim.reshape(-1, 1)

    # Augment with spatial coordinates
    sim_aug = np.concatenate(
        [sim,
         np.mgrid[:sim_dim, :sim_dim].reshape(-1, sim_dim ** 2).T],
        axis=1
    )

    # Empirical metric for number of loop-closures given number of images
    # in sequence (assumption: equally-spaced samples):
    n_clusters = int(np.sqrt(sim_dim))
    print('Performing clustering via KMeans(n={}).'.format(n_clusters))

    km = KMeans(n_clusters=n_clusters, n_jobs=2,
                max_iter=300)
    labels = km.fit_predict(sim_aug)
    print('Got cluster labels')

    for i in range(n_clusters):
        lab_idx = (labels == i)
        if lab_idx.size:
            cc = sim[lab_idx].mean()
            # cc = sim[lab_idx].max()
            sim[lab_idx] = cc

    # Re-normalize and reshape
    sim = sim.reshape(sim_dim, sim_dim) / sim.max()
    return sim 

def generate_deformation_field_rbf(shape, points, max_deform=DEFORMATION_MAX,
                                   nb_bound_points=25):
    """ generate deformation field as thin plate spline  deformation
    in range +/- max_deform

    :param tuple(int,int) shape: tuple of size 2
    :param points: np.array<nb_points, 2> list of landmarks
    :param float max_deform: maximal deformation distance in any direction
    :param int nb_bound_points: number of fix boundary points
    :return: np.array<shape>
    """
    # x_point = points[:, 0]
    # y_point = points[:, 1]
    # generate random shifting
    move = (np.random.random(points.shape[0]) - 0.5) * max_deform

    # fix boundary points
    # set the boundary points
    bound = np.ones(nb_bound_points - 1)
    x_bound = np.linspace(0, shape[0] - 1, nb_bound_points)
    y_bound = np.linspace(0, shape[1] - 1, nb_bound_points)
    x_point = np.hstack((points[:, 0], 0 * bound, x_bound[:-1],
                         (shape[0] - 1) * bound, x_bound[::-1][:-1]))
    y_point = np.hstack((points[:, 1], y_bound[:-1], (shape[1] - 1) * bound,
                         y_bound[::-1][:-1], 0 * bound))
    # the boundary points sex as 0 shift
    move = np.hstack((move, np.zeros(4 * nb_bound_points - 4)))
    # create the interpolation function
    smooth = 0.2 * max_deform
    rbf = interpolate.Rbf(x_point, y_point, move, function='thin-plate',
                          epsilon=1, smooth=smooth)
    # interpolate in regular grid
    x_grid, y_grid = np.mgrid[0:shape[0], 0:shape[1]].astype(np.int32)
    # FIXME: it takes to much of RAM memory, for sample image more that 8GM !
    deform = rbf(x_grid, y_grid)
    return deform 

def sample_normal_expectations(self, gaussian_state):
        """Returns the expectation value E(f) and the variance var(f)
        for some normal distribution X~N(mu, cov).

        Args:
            mu (array): means vector
            cov (array): covariance matrix
            func (function): function acting on the random variables X, P, XP,
                returning a second order polynomial

        Returns:
            tuple: tuple of expectation value and variance.
        """

        def _sample(func, correction=0, mu=None, cov=None):
            """wrapped function"""
            if mu is None:
                mu = gaussian_state[0]

            if cov is None:
                cov = gaussian_state[1]

            X, P = np.mgrid[-7:7:0.01, -7:7:0.01]
            grid = np.dstack((X, P))
            XP = np.prod(grid, axis=2)

            poly = func(X, P, XP)
            PDF = multivariate_normal.pdf(grid, mu, cov)

            Ex = simps(simps(poly * PDF, P[0]), X.T[0])
            ExSq = simps(simps(poly ** 2 * PDF, P[0]), X.T[0])

            var = ExSq - Ex ** 2 + correction

            return Ex, var

        return _sample 

def drawSphere(xCenter, yCenter, zCenter, r):
    # draw sphere
    u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
    x = np.cos(u) * np.sin(v)
    y = np.sin(u) * np.sin(v)
    z = np.cos(v)
    # shift and scale sphere
    x = r * x + xCenter
    y = r * y + yCenter
    z = r * z + zCenter
    return (x, y, z) 

def drawSphere(xCenter, yCenter, zCenter, r):
    # draw sphere
    u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
    x = np.cos(u) * np.sin(v)
    y = np.sin(u) * np.sin(v)
    z = np.cos(v)
    # shift and scale sphere
    x = r * x + xCenter
    y = r * y + yCenter
    z = r * z + zCenter
    return (x, y, z) 

def draw_flow(img, flow, step=16):
    h, w = img.shape[:2]
    y, x = np.mgrid[step / 2:h:step, step / 2:w:step].reshape(2, -1).astype(int)
    fx, fy = flow[y, x].T
    lines = np.vstack([x, y, x + fx, y + fy]).T.reshape(-1, 2, 2)
    lines = np.int32(lines + 0.5)
    vis = cv2.cvtColor(img, cv2.COLOR_GRAY2BGR)
    cv2.polylines(vis, lines, 0, (0, 255, 0))
    for (x1, y1), (x2, y2) in lines:
        cv2.circle(vis, (x1, y1), 1, (0, 255, 0), -1)
    return vis 

def logistic2d(size, radius, roll=True, logres=None):
    """Create a circle with blurred edges

    Set roll=False to have the circle centered in the middle of the
    matrix. Default is to center at the extremes (best for convolution).

    The transition width of the blur scales with the size of the grid.
    I'm not actually sure of the math behind it, but it's what was presented
    in the code from:
    https://0fps.net/2012/11/19/conways-game-of-life-for-curved-surfaces-part-1/
    """
    y, x = size
    # Get coordinate values of each point
    yy, xx = np.mgrid[:y, :x]
    # Distance between each point and the center
    radiuses = np.sqrt((xx - x/2)**2 + (yy - y/2)**2)
    # Scale factor for the transition width
    if logres is None:
        logres = math.log(min(*size), 2)
    with np.errstate(over="ignore"):
        # With big radiuses, the exp overflows,
        # but 1 / (1 + inf) == 0, so it's fine
        logistic = 1 / (1 + np.exp(logres * (radiuses - radius)))
    if roll:
        logistic = np.roll(logistic, y//2, axis=0)
        logistic = np.roll(logistic, x//2, axis=1)
    return logistic 

def path(self, origin, end):
        """
        Gets a 1D array of values for a line connecting `origin` and `end`.
        """
        path_dens = []
        origin = np.array([ float(i) for i in origin ])
        end = np.array([ float(i) for i in end ])
        result = []
        for i in np.mgrid[0:1:50j]:
            point = (1-i)*origin + i*end
            result.append((i,self.interpolate(point)))
        return result 

def deform_image_landmarks(image, points, max_deform=DEFORMATION_MAX):
    """ deform the image by randomly generated deformation field
    and compute new positions for all landmarks

    :param image: np.array<height, width, 3>
    :param points: np.array<nb_points, 2>
    :param float max_deform: maximal deformation distance in any direction
    :return: np.array<height, width, 3>, np.array<nb_points, 2>
    """
    x, y = np.mgrid[0:image.shape[0], 0:image.shape[1]]
    # generate the deformation field
    nb_fix_points = int(np.max(image.shape) / max_deform * 2.)
    x_deform = generate_deformation_field_rbf(image.shape[:2], points,
                                              max_deform, nb_fix_points)
    # TODO: look for another elastic deformation which is friendly to Memory usage
    # -> generate random elastic deformation and using this field get new landmarks
    y_deform = generate_deformation_field_rbf(image.shape[:2], points,
                                              max_deform, nb_fix_points)
    # interpolate the image
    img_warped = interpolate.griddata(zip(x.ravel(), y.ravel()),
                                      image.reshape(-1, 3),
                                      (x + x_deform, y + y_deform),
                                      method='linear', fill_value=1.)
    # compute new positions of landmarks
    x_new = x - x_deform
    y_new = y - y_deform
    pts_warped = np.array([[x_new[pt[0], pt[1]], y_new[pt[0], pt[1]]]
                           for pt in points])
    return img_warped, pts_warped