Python numpy.mgrid() Examples

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 _():
        """ 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 _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 _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 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 __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 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 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 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 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 get_passive_elements(self):
        X, Y = np.mgrid[self.passive_min_x:self.passive_max_x + 1,
            self.passive_min_y:self.passive_max_y]
        pairs = np.vstack([X.ravel(), Y.ravel()]).T
        passive_to_ids = np.vectorize(lambda pair: xy_to_id(*pair,
            nelx=self.nelx - 1, nely=self.nely - 1), signature="(m)->()")
        return passive_to_ids(pairs) 

def _FSpecialGauss(size, sigma):
  """Function to mimic the 'fspecial' gaussian MATLAB function."""
  radius = size // 2
  offset = 0.0
  start, stop = -radius, radius + 1
  if size % 2 == 0:
    offset = 0.5
    stop -= 1
  x, y = np.mgrid[offset + start:stop, offset + start:stop]
  assert len(x) == size
  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))
  return g / g.sum() 

def _FSpecialGauss(size, sigma):
    """Function to mimic the 'fspecial' gaussian MATLAB function."""
    radius = size // 2
    offset = 0.0
    start, stop = -radius, radius + 1
    if size % 2 == 0:
        offset = 0.5
        stop -= 1
    x, y = np.mgrid[offset + start:stop, offset + start:stop]
    assert len(x) == size
    g = np.exp(-((x ** 2 + y ** 2) / (2.0 * sigma ** 2)))
    return g / g.sum() 

def _FSpecialGauss(size, sigma):
    """Function to mimic the 'fspecial' gaussian MATLAB function."""
    radius = size // 2
    offset = 0.0
    start, stop = -radius, radius + 1
    if size % 2 == 0:
        offset = 0.5
        stop -= 1
    x, y = np.mgrid[offset + start:stop, offset + start:stop]
    assert len(x) == size
    g = np.exp(-((x ** 2 + y ** 2) / (2.0 * sigma ** 2)))
    return g / g.sum() 

def get_dp_vertex_doubly_sparse(self, dtype=np.float64, sparseformat=lsofcsr_c, axis=0):
    """ Returns the product vertex coefficients for dominant products as a one-dimensional array of sparse matrices """
    nnz = self.get_dp_vertex_nnz()
    i1,i2,i3,data = zeros(nnz, dtype=int), zeros(nnz, dtype=int), zeros(nnz, dtype=int), zeros(nnz, dtype=dtype)
    atom2s,dpc2s,nfdp,n = self.sv.atom2s, self.dpc2s,self.dpc2s[-1],self.sv.atom2s[-1]

    inz = 0
    for s,f,sd,fd,pt,spp in zip(atom2s,atom2s[1:],dpc2s,dpc2s[1:],self.dpc2t,self.dpc2sp):
      size = self.prod_log.sp2vertex[spp].size
      lv = self.prod_log.sp2vertex[spp].reshape(size)
      dd,aa,bb = np.mgrid[sd:fd,s:f,s:f].reshape((3,size))
      i1[inz:inz+size],i2[inz:inz+size],i3[inz:inz+size],data[inz:inz+size] = dd,aa,bb,lv
      inz+=size

    for sd,fd,pt,spp in zip(dpc2s,dpc2s[1:],self.dpc2t,self.dpc2sp):
      if pt!=2: continue
      inf,(a,b),size = self.bp2info[spp],self.bp2info[spp].atoms,self.bp2info[spp].vrtx.size
      sa,fa,sb,fb = atom2s[a],atom2s[a+1],atom2s[b],atom2s[b+1]
      dd,aa,bb = np.mgrid[sd:fd,sa:fa,sb:fb].reshape((3,size))
      i1[inz:inz+size],i2[inz:inz+size],i3[inz:inz+size] = dd,aa,bb
      data[inz:inz+size] = inf.vrtx.reshape(size)
      inz+=size
      i1[inz:inz+size],i2[inz:inz+size],i3[inz:inz+size] = dd,bb,aa
      data[inz:inz+size] = inf.vrtx.reshape(size)
      inz+=size
    return sparseformat((data, (i1, i2, i3)), dtype=dtype, shape=(nfdp,n,n), axis=axis) 

def get_dp_vertex_sparse(self, dtype=np.float64, sparseformat=coo_matrix):
    """ Returns the product vertex coefficients as 3d array for dominant products, in a sparse format coo(p,ab) by default"""
    nnz = self.get_dp_vertex_nnz()                                                         #Number of non-zero elements in the dominant product vertex
    irow,icol,data = zeros(nnz, dtype=int), zeros(nnz, dtype=int), zeros(nnz, dtype=dtype) # Start to construct coo matrix

    atom2s,dpc2s,nfdp,n = self.sv.atom2s, self.dpc2s,self.dpc2s[-1],self.sv.atom2s[-1]     #self.dpc2s, self.dpc2t, self.dpc2sp: product Center -> list of the size of the basis set in this center,of center's types,of product species

    inz = 0
    for s,f,sd,fd,pt,spp in zip(atom2s,atom2s[1:],dpc2s,dpc2s[1:],self.dpc2t,self.dpc2sp):
      size = self.prod_log.sp2vertex[spp].size
      lv = self.prod_log.sp2vertex[spp].reshape(size)
      dd,aa,bb = np.mgrid[sd:fd,s:f,s:f].reshape((3,size))
      irow[inz:inz+size],icol[inz:inz+size],data[inz:inz+size] = dd, aa+bb*n, lv
      inz+=size

    for sd,fd,pt,spp in zip(dpc2s,dpc2s[1:],self.dpc2t,self.dpc2sp):
      if pt!=2: continue
      inf,(a,b),size = self.bp2info[spp],self.bp2info[spp].atoms,self.bp2info[spp].vrtx.size
      sa,fa,sb,fb = atom2s[a],atom2s[a+1],atom2s[b],atom2s[b+1]
      dd,aa,bb = np.mgrid[sd:fd,sa:fa,sb:fb].reshape((3,size))
      irow[inz:inz+size],icol[inz:inz+size] = dd,aa+bb*n
      data[inz:inz+size] = inf.vrtx.reshape(size)
      inz+=size

      irow[inz:inz+size],icol[inz:inz+size] = dd,bb+aa*n
      data[inz:inz+size] = inf.vrtx.reshape(size)
      inz+=size
      
    return sparseformat((data, (irow, icol)), dtype=dtype, shape=(nfdp,n*n)) 

def get_dp_vertex_sparse2(self, dtype=np.float64, sparseformat=coo_matrix):
    """ Returns the product vertex coefficients as 3d array for dominant products, in a sparse format coo(pa,b) by default"""
    import numpy as np
    from scipy.sparse import csr_matrix, coo_matrix
    nnz = self.get_dp_vertex_nnz()                                                           
    irow,icol,data = np.zeros(nnz, dtype=int), np.zeros(nnz, dtype=int), np.zeros(nnz, dtype=dtype)
    atom2s,dpc2s,nfdp,n = self.sv.atom2s, self.dpc2s, self.dpc2s[-1],self.sv.atom2s[-1]   

    inz = 0
    for s,f,sd,fd,pt,spp in zip(atom2s,atom2s[1:],dpc2s,dpc2s[1:],self.dpc2t,self.dpc2sp):
        size = self.prod_log.sp2vertex[spp].size
        lv = self.prod_log.sp2vertex[spp].reshape(size)
        dd,aa,bb = np.mgrid[sd:fd,s:f,s:f].reshape((3,size))
        irow[inz:inz+size],icol[inz:inz+size],data[inz:inz+size] = dd+aa*nfdp, bb, lv
        inz+=size
           
    for sd,fd,pt,spp in zip(dpc2s,dpc2s[1:],self.dpc2t, self.dpc2sp):
        if pt!=2: continue
        inf,(a,b),size = self.bp2info[spp],self.bp2info[spp].atoms,self.bp2info[spp].vrtx.size
        sa,fa,sb,fb = atom2s[a],atom2s[a+1],atom2s[b],atom2s[b+1]
        dd,aa,bb = np.mgrid[sd:fd,sa:fa,sb:fb].reshape((3,size))
        irow[inz:inz+size],icol[inz:inz+size] = dd+aa*nfdp, bb
        data[inz:inz+size] = inf.vrtx.reshape(size)
        inz+=size
     
        irow[inz:inz+size],icol[inz:inz+size] = dd+bb*nfdp, aa
        data[inz:inz+size] = inf.vrtx.reshape(size)
        inz+=size
    return coo_matrix((data, (irow, icol)), dtype=dtype, shape=(nfdp*n,n)) 

def get_gaussian_weight(self, anchor):
        """
        Args:
            anchor: coordinate of the center
        """
        ret = np.zeros(self.shape, dtype='float32')

        y, x = np.mgrid[:self.shape[0], :self.shape[1]]
        y = y.astype('float32') / ret.shape[0] - anchor[0]
        x = x.astype('float32') / ret.shape[1] - anchor[1]
        g = np.exp(-(x**2 + y ** 2) / self.sigma)
        # cv2.imshow(" ", g)
        # cv2.waitKey()
        return g 

def spatial_distance(array):
    coord = np.sqrt((np.mgrid[0:windowSize,0:windowSize]-windowSize//2)**2)
    spat_dist = np.sqrt(((0-coord[0])**2+(0-coord[1])**2))
    rel_spat_dist = spat_dist/spatImp + 1.0 # relative spatial distance
    rev_spat_dist = 1/rel_spat_dist # relative spatial distance reversed
    flat_spat_dist = np.ravel(rev_spat_dist)
    spat_dist_da = da.from_array(flat_spat_dist, chunks=flat_spat_dist.shape)
    print ("Done spatial distance!", spat_dist_da)
    
    return spat_dist_da


# Define the threshold used in the dynamic classification process 

def GridIter(self, PT = None, QM = 1.0):
        if PT is None:
            PT = self.TM
        #Get the range of values to compute grid for
        R = self.Solve3DT(np.array([0, self.w, 0, self.w]), np.array([0, 0, self.h, self.h]))
        minx, miny, _ = np.floor(R.min(axis = 0) / QM) * QM
        maxx, maxy, _ = np.ceil(R.max(axis = 0) / QM) * QM
        GP = np.zeros((int((maxx - minx) * (maxy - miny) / QM), 4))
        GP[:, 0:2] = np.mgrid[minx:maxx:QM, miny:maxy:QM].reshape(2, -1).T
        GP[:, 3] = 1
        PR = np.dot(GP, PT)
        PR = (PR[:, 0:2] / PR[:, [2]])
        ind = ((PR >= np.array([0, 0])) & (PR < np.array([self.w, self.h]))).sum(axis = 1) == 2
        return GP[ind, 0:3], PR[ind] 

def CellCorners(self):
        '''
        Gets the top left corners of the CNN prediction cells in pixels (x, y)
        '''
        return np.mgrid[self.sb[0]:(self.ss[0] + self.sb[0] + 1):self.cs[0], self.sb[1]:(self.ss[1] + self.sb[1] + 1):self.cs[1]].reshape(2, -1).T 

def live_calibrate(camera, pattern_shape, n_matches_needed):
    """ Find calibration parameters as the user moves a checkerboard in front of the camera """
    print("Looking for %s checkerboard" % (pattern_shape,))
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
    example_3d = np.zeros((pattern_shape[0] * pattern_shape[1], 3), np.float32)
    example_3d[:, :2] = np.mgrid[0 : pattern_shape[1], 0 : pattern_shape[0]].T.reshape(-1, 2)
    points_3d = []
    points_2d = []
    while len(points_3d) < n_matches_needed:
        ret, frame = camera.cap.read()
        assert ret
        gray_frame = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
        ret, corners = cv2.findCirclesGrid(
            gray_frame, pattern_shape, flags=cv2.CALIB_CB_ASYMMETRIC_GRID
        )
        cv2.imshow("camera", frame)
        if ret:
            points_3d.append(example_3d.copy())
            points_2d.append(corners)
            print("Found calibration %i of %i" % (len(points_3d), n_matches_needed))
            drawn_frame = cv2.drawChessboardCorners(frame, pattern_shape, corners, ret)
            cv2.imshow("calib", drawn_frame)
        cv2.waitKey(10)
    ret, camera_matrix, distortion_coefficients, _, _ = cv2.calibrateCamera(
        points_3d, points_2d, gray_frame.shape[::-1], None, None
    )
    assert ret
    return camera_matrix, distortion_coefficients 

def plot_mcontour(self, ndim0, ndim1, z, show_mode):
        "use mayavi.mlab to plot contour."
        if not mayavi_installed:
            self.__logger.info("Mayavi is not installed on your device.")
            return
        #do 2d interpolation
        #get slice object
        s = np.s_[0:ndim0:1, 0:ndim1:1]
        x, y = np.ogrid[s]
        mx, my = np.mgrid[s]
        #use cubic 2d interpolation
        interpfunc = interp2d(x, y, z, kind='cubic')
        newx = np.linspace(0, ndim0, 600)
        newy = np.linspace(0, ndim1, 600)
        newz = interpfunc(newx, newy)
        #mlab
        face = mlab.surf(newx, newy, newz, warp_scale=2)
        mlab.axes(xlabel='x', ylabel='y', zlabel='z')
        mlab.outline(face)
        #save or show
        if show_mode == 'show':
            mlab.show()
        elif show_mode == 'save':
            mlab.savefig('mlab_contour3d.png')
        else:
            raise ValueError('Unrecognized show mode parameter : ' +
                             show_mode)

        return 

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

def random_warp( in_image ):
    assert in_image.shape[:2] == (256,256)

    image = in_image.copy()


    scale = 5

    range_ = numpy.linspace( 128-120, 128+120, scale )
    mapx = numpy.broadcast_to( range_, (scale,scale) )
    mapy = mapx.T

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

    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[:,:,:3], interp_mapx, interp_mapy, cv2.INTER_CUBIC )

    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) )

    target_mask = target_image[:,:,3].reshape((64,64,1))
    target_image = target_image[:,:,:3]


    if len(target_image.shape)>2:
      return ( warped_image, 
               target_image, 
               target_mask )
    else:
      return ( warped_image, 
               target_image ) 

def _init_coords(self):
    idxs = np.mgrid[[slice(0, x) for x in self.canvas.image.shape]]
    sort_idx = np.argsort(self.canvas.image.flat)[::-1]
    self.coords = np.array(zip(*[idx.flat[sort_idx] for idx in idxs]))