def get_gauss_points(self, ng):
     if isinstance(ng, int):
         return self.gauss_points.get(ng)
     elif isinstance(ng, list):
         if len(ng) > 3:
             raise Exception('Gauss points for 4 dimensions' +
                     'and above not supported')
         if len(ng) == 2:
             Xi1, W1 = self.get_gauss_points(ng[0])
             Xi2, W2 = self.get_gauss_points(ng[1])
             Xi1g, Xi2g = numpy.meshgrid(Xi1.flatten(), Xi2.flatten())
             Xi1 = numpy.array([Xi1g.flatten(), Xi2g.flatten()]).T
             W1g, W2g = numpy.meshgrid(W1.flatten(), W2.flatten())
             W1 = W1g.flatten() * W2g.flatten()
             return Xi1, W1
         elif len(ng) == 3:
             Xi1, W1 = self.get_gauss_points(ng[0])
             Xi2, W2 = self.get_gauss_points(ng[1])
             Xi3, W3 = self.get_gauss_points(ng[2])
             gindex = numpy.mgrid[0:ng[0], 0:ng[1], 0:ng[2]]
             gindex = numpy.array([
                 gindex[2, :, :].flatten(),
                 gindex[1, :, :].flatten(),
                 gindex[0, :, :].flatten()]).T
             Xi = numpy.array([
                 Xi1[gindex[:,0]], Xi2[gindex[:,1]], Xi3[gindex[:,2]]])[:,:,0].T
             W = numpy.array([
                 W1[gindex[:,0]], W2[gindex[:,1]], W3[gindex[:,2,]]]).T.prod(1)
             return Xi, W
         
     raise Exception('Invalid number of gauss points')
     return None, None    

def get_patches(data, var, dql, centers, patch_size=5):
    """
    Return patched data, etc, at given centers.
    """
    buff = (patch_size - 1) / 2

    # SExtractor doesnt alway put peak on brightest pixel
    # find shifts, up to 1 pixel that correct for this
    x, y = np.meshgrid(range(-2, 3),
                       range(-2, 3))
    xcs = (x[None, :, :] + centers[:,0][:, None, None]).astype(np.int)
    ycs = (y[None, :, :] + centers[:,1][:, None, None]).astype(np.int)
    
    for i in range(centers.shape[0]):
        test = data[xcs[i], ycs[i]]
        ind = (test == np.max(test))
        centers[i, 0] += x[ind]
        centers[i, 1] += y[ind]

    # Full patch indices
    x, y = np.meshgrid(range(-buff, buff + 1),
                       range(-buff, buff + 1))
    xcs = (x[None, :, :] + centers[:,0][:, None, None]).astype(np.int)
    ycs = (y[None, :, :] + centers[:,1][:, None, None]).astype(np.int)
    
    return data[xcs, ycs], var[xcs, ycs], dql[xcs, ycs], xcs, ycs, \
        centers

	def logLikelihood(self, data):
		assert self.numImages == data.numImages

		which = [np.nonzero(data.id == i)[0]\
				for i in xrange(0, self.numImages)]

		# Mean vector
		m = np.empty(data.t.size)
		for i in xrange(0, self.numImages):
			m[which[i]] = self.mag[i]

		# Covariance matrix
		[t1, t2] = np.meshgrid(data.t, data.t)
		lags = np.abs(t2 - t1)
		plt.imshow(lags)
		plt.show()
		ids = np.meshgrid(data.id, data.id)


		equal = ids[0] == ids[1]

##		try:
##			L = la.cholesky(C)
##		except:
##			return -np.inf
#		y = data.y - m
#		logDeterminant = 2.0*np.sum(np.log(np.diag(L)))
#		solution = la.cho_solve((L, True), y)
#		exponent = np.dot(y, solution)
#		logL = -0.5*data.t.size*np.log(2.0*np.pi) - 0.5*logDeterminant - 0.5*exponent
		return 0.

def P1(mode="a"):
    if mode=="a":
        fig = plt.figure()
        for N,p in [(8,221),(16,222),(32,223),(64,224)]:
            u = OneD(0.035,2.5,60,20,20,N,0.01,5)
            x = [i*2.5/N for i in range(N)]
            y = [j*0.01 for j in range(int(5/0.01))]
            X,Y = np.meshgrid(x,y)
            ax = fig.add_subplot(p, projection='3d')
            surf = ax.plot_surface(X,Y,u,rstride=1, cstride=1,linewidth=0, cmap=cm.coolwarm)
            ax.set_zlim(20,60.01)
            ax.set_xlabel('X Position (cm)')
            ax.set_ylabel('Time (s)')
            ax.set_zlabel('Temperature (C)')
            ax.set_title('N='+str(N))
            fig.colorbar(surf,ticks=[t for t in range(20,61,5)])
        plt.show()
    if mode == "b":
        fig = plt.figure()
        for t,p in [(0.001,221),(0.01,222),(0.1,223)]:
            u = OneD(0.035,2.5,60,20,20,32,t,5)
            x = [i*2.5/32 for i in range(32)]
            y = [j*t for j in range(int(5/t))]
            X,Y = np.meshgrid(x,y)
            ax = fig.add_subplot(p, projection='3d')
            surf = ax.plot_surface(X,Y,u,rstride=1, cstride=1,linewidth=0, cmap=cm.coolwarm)
            if t!=0.1:
                ax.set_zlim(20,60.01)
            ax.set_xlabel('X Position (cm)')
            ax.set_ylabel('Time (s)')
            ax.set_zlabel('Temperature (C)')
            ax.set_title('Timestep:'+str(t))
            if t != 0.1:
                fig.colorbar(surf,ticks=[t for t in range(20,61,5)])
        plt.show()

def makePathAndGrid(data, xsteps, ysteps, x_col='x_m', y_col='y_m'):
    x = data[x_col]
    y = data[y_col]

    xmin = min(x)
    xmax = max(x)
    ymin = min(y)
    ymax = max(y)

    # make grid
    dx = (xmax - xmin)/xsteps
    dy = (ymax - ymin)/ysteps
    XXraw, YYraw = np.meshgrid(np.arange(xmin, xmax, dx) + dx/2, np.arange(ymin, ymax, dy) + dy/2)

    x = ((x - xmin)/(xmax-xmin) * xsteps).astype(int)
    y = ((y - ymin)/(ymax-ymin) * ysteps).astype(int)

    xx = np.arange(0, xsteps)
    yy = np.arange(0, ysteps)

    M = np.zeros((len(xx),len(yy)))

    N = len(xx) * len(yy)
    idx = np.arange(N)
    np.random.shuffle(idx)
    XX, YY = np.meshgrid(xx, yy)

    XX = XX.reshape(N)
    YY = YY.reshape(N)
    XXraw = XXraw.reshape(N)
    YYraw = YYraw.reshape(N)

    return idx, np.c_[XXraw, YYraw].T, np.c_[XX, YY].T, M

def stmDatToArray(filename):
    with open(filename, 'r') as f:
        f.readline()
        delta = []
        for i in range(2):
            line = map(float, [entry for entry in f.readline().strip().split(' ') if entry!=''][:2])
            delta.append(line)

        shape = [entry for entry in f.readline().strip().split(' ') if len(entry)!=0]
        shape = map(int, shape[:2])

        arrays = []
        xm, ym = [delta[i][0] + np.arange(shape[i])*delta[i][1] for i in range(2)]
        xo, yo = [delta[i][0] - 0.5*delta[i][1]+ np.arange(shape[i] + 1)*delta[i][1] for i in range(2)]

        Xm, Ym = np.meshgrid(xm, ym)
        Xo, Yo = np.meshgrid(xo, yo)

        collect = []
        for line in f:
            try:
                line = [entry for entry in line.strip().split(' ') if len(entry)>0]
                line = map(float, line)
            except ValueError:
                break
            else:
                collect.extend(line)

        Z = np.reshape(collect, shape)

        return Xm, Ym, Z

 def _calc_partials(self):
     """Calculate the partial derivatives
     
     Uses sympy.utilities.lambdify to convert equation into
     a function that can be easily computed.
     
     Raises:
       MissingEquationError: if method is called before equation is
         attached.
     """
     if not self.equation:
         raise self.MissingEquationError('No equation attached')
     
     x, y = sympy.symbols('x,y')
     compute_func = sympy.utilities.lambdify((x, y), self.equation)
     X, Y = np.meshgrid(self.xrange, self.yrange)
     DX, DY = np.meshgrid(np.zeros(len(self.xrange)),
                          np.zeros(len(self.yrange)))
     
     # Iterate through grid and compute function value at each point
     # If value cannot be computed, default to 0
     # If value can be computed, scale by sqrt of the magnitude
     for i, a in enumerate(self.xrange):
         for j, b in enumerate(self.yrange):
             dx = 1
             try:
                 dy = compute_func(a, b)
                 n = sqrt(dx**2 + dy**2)
                 dy /= sqrt(n)
                 dx /= sqrt(n)
                 DX[j][i] = dx
                 DY[j][i] = dy
             except (ValueError, ZeroDivisionError):
                 pass
     return X, Y, DX, DY

    def __init__(self):
        # sift features
        self.gS = 8
        self.pS = 16
        self.nrml_thres = 1.0
        self.sigma = 0.8
        self.sift_thres = 0.2
        self.Nangles = 8
        self.Nbins = 4
        self.Nsamples = self.Nbins**2
        self.alpha = 9.0
        self.angles = np.array(range(self.Nangles))*2.0*np.pi/self.Nangles

        self.GH, self.GW = self.gen_dgauss(self.sigma)

        # compute the weight contribution map
        # weights is the contribution of each pixel to the corresponding bin center
        sample_res = self.pS / np.double(self.Nbins)
        sample_p = np.array(range(self.pS))
        sample_ph, sample_pw = np.meshgrid(sample_p,sample_p)
        sample_ph.resize(sample_ph.size)
        sample_pw.resize(sample_pw.size)
        bincenter = np.array(range(1,self.Nbins*2,2)) / 2.0 / self.Nbins * self.pS - 0.5
        bincenter_h, bincenter_w = np.meshgrid(bincenter,bincenter)
        bincenter_h.resize((bincenter_h.size,1))
        bincenter_w.resize((bincenter_w.size,1))
        dist_ph = abs(sample_ph - bincenter_h)
        dist_pw = abs(sample_pw - bincenter_w)
        weights_h = dist_ph / sample_res
        weights_w = dist_pw / sample_res
        weights_h = (1-weights_h) * (weights_h <= 1)
        weights_w = (1-weights_w) * (weights_w <= 1)
        self.weights = weights_h * weights_w

def plot_beam(dirname, input_beam, Rho, Phi, ref_output_fluence):
    filename = lambda name: os.path.join(dirname, name)
    
    if len(Rho) > 1:
        vfluence = np.vectorize(input_beam.fluence)
        ref_input_fluence = vfluence(*np.meshgrid(Rho, Phi)).T
        norm_input_fluence = ref_input_fluence / input_beam.ref_fluence
        norm_output_fluence = ref_output_fluence / input_beam.ref_fluence
        max_output_fluence = np.amax(norm_output_fluence)
        
        n_ref = -1
        for n, phi in enumerate(Phi):
            if n_ref < 0 or abs(phi - input_beam.phi_ref) < abs(Phi[n_ref] - input_beam.phi_ref):
                n_ref = n
        rholim = (Rho[0], Rho[-1])
        plot.plot_data(filename("fluences"), "Input and Output Fluence", ((Rho,)*2, None, rholim, rho_label), ((norm_input_fluence[:, n_ref], norm_output_fluence[:, n_ref]), None, None, fluence_rel_label), ("input beam", "output beam"))
        plot.plot_data(filename("fluences_norm"), "Normalized Input and Output Fluence", ((Rho,)*2, None, rholim, rho_label), ((norm_input_fluence[:, n_ref], norm_output_fluence[:, n_ref] / max_output_fluence), None, None, fluence_norm_rel_label), ("input beam", "output beam"))
        
        if len(Phi) > 1:
            FR, RF = np.meshgrid(Phi, Rho)
            XY, YX = RF * np.cos(FR), RF * np.sin(FR)
            stride_rho = max(len(Rho) // params.out_count_rho, 1)
            stride_phi = max(len(Phi) // params.out_count_phi, 1)
            plot.plot_projection(filename("fluence_in"), "Input Fluence", (XY, None, x_label), (YX, None, y_label), (norm_input_fluence, None, fluence_rel_label), (30, -60), (stride_rho, stride_phi))
            plot.plot_projection(filename("fluence_out"), "Output Fluence", (XY, None, x_label), (YX, None, y_label), (norm_output_fluence, None, fluence_rel_label), (30, -60), (stride_rho, stride_phi))

def test_label_propagation_closed_form():
    n_classes = 2
    X, y = make_classification(n_classes=n_classes, n_samples=200,
                               random_state=0)
    y[::3] = -1
    Y = np.zeros((len(y), n_classes + 1))
    Y[np.arange(len(y)), y] = 1
    unlabelled_idx = Y[:, (-1,)].nonzero()[0]
    labelled_idx = (Y[:, (-1,)] == 0).nonzero()[0]

    clf = label_propagation.LabelPropagation(max_iter=10000,
                                             gamma=0.1)
    clf.fit(X, y)
    # adopting notation from Zhu et al 2002
    T_bar = clf._build_graph()
    Tuu = T_bar[tuple(np.meshgrid(unlabelled_idx, unlabelled_idx,
                      indexing='ij'))]
    Tul = T_bar[tuple(np.meshgrid(unlabelled_idx, labelled_idx,
                                  indexing='ij'))]
    Y = Y[:, :-1]
    Y_l = Y[labelled_idx, :]
    Y_u = np.dot(np.dot(np.linalg.inv(np.eye(Tuu.shape[0]) - Tuu), Tul), Y_l)

    expected = Y.copy()
    expected[unlabelled_idx, :] = Y_u
    expected /= expected.sum(axis=1)[:, np.newaxis]

    assert_array_almost_equal(expected, clf.label_distributions_, 4)

    def graph(self, x_axis, y_axis, image, params, result = None):
        #plot the sample data
        from matplotlib import pyplot
        pyplot.figure()
        pyplot.subplot(2,1,1)
        pyplot.contourf(x_axis, y_axis, image, alpha = 0.5)
        #plot the fit
        #sample the fit with more precision

        x_axis_fit = np.linspace(x_axis.min(), x_axis.max(), x_axis.size * 5)
        y_axis_fit = np.linspace(y_axis.min(), y_axis.max(), y_axis.size * 5)
        xx, yy = np.meshgrid(x_axis_fit, y_axis_fit)
        params['background_level'].value = 0                                                                                
        fit_1d = self.ion_model(params, xx, yy, use_1d=True)
        xx, yy = np.meshgrid(x_axis, y_axis)
        pyplot.subplot(2,1,2)
        ximage = image.sum(axis=0)
        ximage = ximage - ximage[0]
        a1d = np.ones_like(self.allions_1d)
        #low_values_indices = self.allions_1d - bg < (0.5)*params['amplitude'].value  # Where values are low
        low_values_indices = self.allions_1d < self.model_threshold  # Where values are low
        a1d[low_values_indices] = 0  # All low values set to 0
        pyplot.plot(x_axis, a1d * ximage.max())
        #pyplot.plot(x_axis, ximage * a1d, 'o', markersize=10)
        pyplot.plot(x_axis, ximage)
        pyplot.plot(x_axis_fit, fit_1d)
        pyplot.tight_layout()
        pyplot.show()

    def __init__(self, gridSpacing, patchSize, nrml_thres=1.0, sigma_edge=0.8, sift_thres=0.2):
        """
gridSpacing: the spacing for sampling dense descriptors
patchSize: the size for each sift patch
nrml_thres: low contrast normalization threshold
sigma_edge: the standard deviation for the gaussian smoothing
before computing the gradient
sift_thres: sift thresholding (0.2 works well based on
Lowe's SIFT paper)
"""
        self.gS = gridSpacing
        self.pS = patchSize
        self.nrml_thres = nrml_thres
        self.sigma = sigma_edge
        self.sift_thres = sift_thres
        # compute the weight contribution map
        sample_res = self.pS / np.double(Nbins)
        sample_p = np.array(range(self.pS))
        sample_ph, sample_pw = np.meshgrid(sample_p, sample_p)
        sample_ph.resize(sample_ph.size)
        sample_pw.resize(sample_pw.size)
        bincenter = np.array(range(1, Nbins * 2, 2)) / 2.0 / Nbins * self.pS - 0.5
        bincenter_h, bincenter_w = np.meshgrid(bincenter, bincenter)
        bincenter_h.resize((bincenter_h.size, 1))
        bincenter_w.resize((bincenter_w.size, 1))
        dist_ph = abs(sample_ph - bincenter_h)
        dist_pw = abs(sample_pw - bincenter_w)
        weights_h = dist_ph / sample_res
        weights_w = dist_pw / sample_res
        weights_h = (1 - weights_h) * (weights_h <= 1)
        weights_w = (1 - weights_w) * (weights_w <= 1)
        # weights is the contribution of each pixel to the corresponding bin center
        self.weights = weights_h * weights_w

    def _make_carts_dict(self):
        """ Return a carts dictionary, distances in meters. """
        az_deg, range_km, el_deg = self._calc_geometry()
        # simple calculation involving 4/3 earth radius
        nsweeps = self.master_header['max_nz']
        nrays = self.master_header['max_ny']
        ngates = self.master_header['max_nx']
        xx = np.empty([nsweeps, nrays, ngates], dtype=np.float32)
        yy = np.empty([nsweeps, nrays, ngates], dtype=np.float32)
        zz = np.empty([nsweeps, nrays, ngates], dtype=np.float32)

        if self.projection == 'rhi':
            rg, ele = np.meshgrid(range_km, el_deg)
            rg = np.array(rg, dtype=np.float64)
            ele = np.array(ele, dtype=np.float64)
            for aznum in range(nsweeps):
                azg = np.ones(rg.shape, dtype=np.float64) * az_deg[aznum]
                x, y, z = antenna_to_cartesian(rg, azg, ele)
                zz[aznum, :, :] = z
                xx[aznum, :, :] = x
                yy[aznum, :, :] = y

        elif self.projection == 'ppi':
            rg, azg = np.meshgrid(range_km, az_deg)
            rg = np.array(rg, dtype=np.float64)
            azg = np.array(azg, dtype=np.float64)
            for elnum in range(nsweeps):
                ele = np.ones(rg.shape, dtype=np.float64) * el_deg[elnum]
                x, y, z = antenna_to_cartesian(rg, azg, ele)
                zz[elnum, :, :] = z
                xx[elnum, :, :] = x
                yy[elnum, :, :] = y

        return {'x': xx, 'y': yy, 'z': zz}

    def find_intersections(self, A, B):
        ''' (matrix, matrix -> bool'''
        amin = lambda x1, x2: np.where(x1<x2, x1, x2)
        amax = lambda x1, x2: np.where(x1>x2, x1, x2)
        aall = lambda abools: np.dstack(abools).all(axis=2)
        slope = lambda line: (lambda d: d[:,1]/d[:,0])(np.diff(line, axis=0))

        x11, x21 = np.meshgrid(A[:-1, 0], B[:-1, 0])
        x12, x22 = np.meshgrid(A[1:, 0], B[1:, 0])
        y11, y21 = np.meshgrid(A[:-1, 1], B[:-1, 1])
        y12, y22 = np.meshgrid(A[1:, 1], B[1:, 1])

        m1, m2 = np.meshgrid(slope(A), slope(B))
        m1inv, m2inv = 1/m1, 1/m2

        yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
        xi = (yi - y21)*m2inv + x21

        xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12), 
                  amin(x21, x22) < xi, xi <= amax(x21, x22) )
        yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
                  amin(y21, y22) < yi, yi <= amax(y21, y22) )
        if xi[aall(xconds)] and yi[aall(yconds)]:
            return  xi[aall(xconds)], yi[aall(yconds)], True # intersection, don't go!
        return  xi[aall(xconds)], yi[aall(yconds)], False

def find_intersections(contour1, contour2):
    """
    Vectorized code to find intersections between contours. All
    successive duplicate vertices along the input contours must be
    removed to help avoid division-by-zero errors.

    There are cases were no intersection exists (eg. parallel lines)
    where division by zero and invalid value exceptions occur. These
    exceptions should be caught as warnings: these edge cases are
    unavoidable with this algorithm and do not indicate that the
    output is erroneous.
    """
    amin = lambda x1, x2: np.where(x1<x2, x1, x2)       # Elementwise min selection
    amax = lambda x1, x2: np.where(x1>x2, x1, x2)       # Elementwise max selection
    aall = lambda abools: np.dstack(abools).all(axis=2) # dstacks, checks True depthwise
    # Uses delta (using np.diff) to find successive slopes along path
    slope = lambda line: (lambda d: d[:,1]/d[:,0])(np.diff(line, axis=0))
    # Meshgrids between both paths (x and y). One element sliced off end/beginning
    x11, x21 = np.meshgrid(contour1[:-1, 0], contour2[:-1, 0])
    x12, x22 = np.meshgrid(contour1[1:, 0], contour2[1:, 0])
    y11, y21 = np.meshgrid(contour1[:-1, 1], contour2[:-1, 1])
    y12, y22 = np.meshgrid(contour1[1:, 1], contour2[1:, 1])
    # Meshgrid of all slopes for both paths
    m1, m2 = np.meshgrid(slope(contour1), slope(contour2))
    m2inv = 1/m2 # m1inv was not used.
    yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
    xi = (yi - y21)*m2inv + x21 # (xi, yi) is intersection candidate
    # Bounding box type conditions for intersection candidates
    xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12),
              amin(x21, x22) < xi, xi <= amax(x21, x22) )
    yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
              amin(y21, y22) < yi, yi <= amax(y21, y22) )
    return xi[aall(xconds)], yi[aall(yconds)]

def generate_ind_pairs(img_num_per_cam_person, person_inds, cam_inds):
    positive_pair_a = numpy.asarray([], dtype='int32')
    positive_pair_b = numpy.asarray([], dtype='int32')
    negative_pair_a = numpy.asarray([], dtype='int32')
    negative_pair_b = numpy.asarray([], dtype='int32')

    past_person_inds = numpy.asarray([], dtype='int32')
    for i in person_inds:
        past_person_inds = numpy.append(past_person_inds, i)
        past_cam_inds = numpy.asarray([], dtype='int32')
        for j in cam_inds:
            past_cam_inds = numpy.append(past_cam_inds, j)
            target_abs_inds = get_abs_ind_one(img_num_per_cam_person, i, j)
            pos_cand_abs_inds = get_abs_ind_per_cam_person(img_num_per_cam_person,
                                                           numpy.asarray([i], dtype='int32'),
                                                           numpy.setdiff1d(cam_inds, past_cam_inds))
            neg_cand_abs_inds = get_abs_ind_per_cam_person(img_num_per_cam_person,
                                                       numpy.setdiff1d(person_inds, past_person_inds),
                                                       numpy.setdiff1d(cam_inds, j))

            [tmp_a, tmp_b] = numpy.meshgrid(pos_cand_abs_inds, target_abs_inds)
            positive_pair_a = numpy.append(positive_pair_a, tmp_b.flatten())
            positive_pair_b = numpy.append(positive_pair_b, tmp_a.flatten())
            [tmp_a, tmp_b] = numpy.meshgrid(neg_cand_abs_inds, target_abs_inds)
            negative_pair_a = numpy.append(negative_pair_a, tmp_b.flatten())
            negative_pair_b = numpy.append(negative_pair_b, tmp_a.flatten())

    return [positive_pair_a, positive_pair_b, negative_pair_a, negative_pair_b]

def mask_polar_to_cart(mask, center, min_radius, max_radius, output_shape, zoom_factor=1):
    '''Converts a polar binary mask to Cartesian and places in an image of zeros'''

    # Account for upsampling
    if zoom_factor != 1:
        center = (center[0]*zoom_factor + zoom_factor/2, center[1]*zoom_factor + zoom_factor/2)
        min_radius = min_radius * zoom_factor
        max_radius = max_radius * zoom_factor
        output_shape = map(lambda a: a * zoom_factor, output_shape)

    # new image
    image = np.zeros(output_shape)

    # coordinate conversion
    theta, r = np.meshgrid(np.linspace(0, 2*np.pi, mask.shape[1]),
                           np.arange(0, max_radius))
    x, y = coord_polar_to_cart(r, theta, center)
    x, y = np.round(x), np.round(y)
    x, y = x.astype(int), y.astype(int)

    x = np.clip(x, 0, image.shape[0]-1)
    y = np.clip(y, 0, image.shape[1]-1)
    ix,iy = np.meshgrid(np.arange(0,mask.shape[1]), np.arange(0,mask.shape[0]))
    image[x,y] = mask

    # downsample image
    if zoom_factor != 1:
        zf = 1/float(zoom_factor)
        image = zoom(image, (zf, zf), order=4)

    # ensure image remains a filled binary mask
    image = (image > 0.5).astype(int)
    image = binary_fill_holes(image)
    return image

    def _init_grid(self):
        """ Initialize (x,y,z) coordinate grid + (re)draw plot."""
        phi_UC = np.linspace(0, 2*pi, 400, endpoint=True) # angles for unit circle
        self.xy_UC = np.exp(1j * phi_UC) # x,y coordinates of unity circle

        steps = 100              # number of steps for x, y, r, phi
        #
        self.xmin = -1.5; self.xmax = 1.5  # cartesian range limits
        self.ymin = -1.5; self.ymax = 1.5

        rmin = 0;    rmax = self.xmin  # polar range limits

        # Calculate grids for 3D-Plots
        dr = rmax / steps * 2 # grid size for polar range
        dx = (self.xmax - self.xmin) / steps
        dy = (self.ymax - self.ymin) / steps # grid size cartesian range

        if self.chkPolar.isChecked(): # # Plot circular range in 3D-Plot
            [r, phi] = np.meshgrid(np.arange(rmin, rmax, dr),
                            np.linspace(0, 2 * pi, steps, endpoint=True))
            self.x = r * cos(phi)
            self.y = r * sin(phi)
        else: # cartesian grid
            [self.x, self.y] = np.meshgrid(np.arange(self.xmin, self.xmax, dx),
                                            np.arange(self.ymin, self.ymax, dy))

        self.z = self.x + 1j*self.y # create coordinate grid for complex plane

        self.draw() # initial plot

def polyval2d_general(c, x, y, function="polynomial", minx=None, maxx=None, miny=None, maxy=None):
    if function == "polynomial":
        xx, yy = np.meshgrid(x, y)
        return np.polynomial.polynomial.polyval2d(xx, yy, c)
    elif function in ["legendre", "chebyshev"]:
        # Scale x-direction
        if minx is None or maxx is None:
            if np.size(x) == 1:
                xmin, xmax = -1.0, 1.0
            else:
                xmin, xmax = np.min(x), np.max(x)
        else:
            xmin, xmax = minx, maxx
        xv = 2.0 * (x-xmin)/(xmax-xmin) - 1.0
        # Scale y-direction
        if miny is None or maxy is None:
            if np.size(y) == 1:
                ymin, ymax = -1.0, 1.0
            else:
                ymin, ymax = np.min(y), np.max(y)
        else:
            ymin, ymax = miny, maxy
        yv = 2.0 * (y-ymin)/(ymax-ymin) - 1.0
        xx, yy = np.meshgrid(xv, yv)
        if function == "legendre":
            return np.polynomial.legendre.legval2d(xx, yy, c)
        elif function == "chebyshev":
            return np.polynomial.chebyshev.chebval2d(xx, yy, c)
    else:
        msgs.error("Function {0:s} has not yet been implemented".format(function))
    return None

    def pflux(self, p, x):
        Egam = x * 1e-9  # in TeV

        if self.EG is None:
            self.EG, self.EP = np.meshgrid(Egam, self.Ep)
            self.Fgam = self.F_gamma(self.EG, self.EP)
        elif np.array_equal(Egam, self.EG[0, :]) is False:
            print(Egam.shape, self.EG.shape)
            print("Warning different internal vectors. Recomputing.")
            self.EG, self.EP = np.meshgrid(Egam, self.Ep)
            self.Fgam = self.F_gamma(self.EG, self.EP)

        # proton_spec = p[3]*np.exp(-self.EP*p[1]*1e-3)*(self.EP/p[0])**(-p[2])
        proton_spec = np.exp(-self.EP * p[1] * 1e-3) * (self.EP / p[0]) ** (-p[2])

        if self.Fgam is None:
            self.Fgam = self.F_gamma(self.EG, self.EP)

        res = self.Fgam * proton_spec
        integral = 3e10 * res.sum(0) * self.lbsize / np.log(10.)

        # Normalize with proton energy content at 1 kpc for n_H = 1 cm^-3
        pe_spec = self.Ep ** 2 * np.exp(-self.Ep * p[1] * 1e-3) * (self.Ep / p[0]) ** (-p[2])
        pe_spec = pe_spec[(self.Ep >= p[4]) * (self.Ep <= p[5])]
        norm = 4 * np.pi * (3.1e21) ** 2 * 1.602e9 * pe_spec.sum() * self.lbsize / np.log(10.) / 1e50

        #        return 3e10*res.sum(0)*self.lbsize/np.log(10.)*1e-9   # per TeV not keV
        return integral / norm * p[3]