python类gaussian_filter()的实例源码

def get_instance(test=False):
  ge_fail, link_fail = random_fail(G_E)
  path_chg = path_changed(G_E, ge_fail, G_OBS)


  return link_fail, path_chg, ge_fail 
  # img, _x = gen_crack()
  # img = gaussian_filter(img, 1.0)
  # return img, _x

# a trace is named tuple
# (Img, S, Os) 
# where Img is the black/white image
# where S is the hidden hypothesis (i.e. label of the img)
# Os is a set of Observations which is (qry_pt, label)

def _get_smoothed_histogram(self, chain, parameter):
        data = chain.get_data(parameter)
        smooth = chain.config["smooth"]
        if chain.grid:
            bins = get_grid_bins(data)
        else:
            bins = chain.config['bins']
            bins, smooth = get_smoothed_bins(smooth, bins, data, chain.weights)

        hist, edges = np.histogram(data, bins=bins, normed=True, weights=chain.weights)
        edge_centers = 0.5 * (edges[1:] + edges[:-1])
        xs = np.linspace(edge_centers[0], edge_centers[-1], 10000)

        if smooth:
            hist = gaussian_filter(hist, smooth, mode=self.parent._gauss_mode)
        kde = chain.config["kde"]
        if kde:
            kde_xs = np.linspace(edge_centers[0], edge_centers[-1], max(200, int(bins.max())))
            ys = MegKDE(data, chain.weights, factor=kde).evaluate(kde_xs)
            area = simps(ys, x=kde_xs)
            ys = ys / area
            ys = interp1d(kde_xs, ys, kind="linear")(xs)
        else:
            ys = interp1d(edge_centers, hist, kind="linear")(xs)
        cs = ys.cumsum()
        cs /= cs.max()
        return xs, ys, cs

def draw_npcs(all_npcs, npc_locations, input_path, output_path, mark_npcs=None,
              relative=False, filter_sigma=10, sigmoid_k=8, sigmoid_c=0.1):

    map_orig = PIL.Image.open(input_path)

    # Get the overlay for NPCs we wish to plot
    if not mark_npcs:
        mark_npcs = all_npcs    
    overlay = get_overlay_for(mark_npcs, npc_locations, map_orig.size)

    # If we want a relative fraction, we also want an overlay for the entire population
    if relative:
        overlay_population = get_overlay_for(all_npcs, npc_locations, map_orig.size)
        overlay = np.nan_to_num(np.divide(overlay, gaussian_filter(overlay_population, filter_sigma)))

    # Blur and normalize it (TODO if relative, we've already blurred the overall
    # population overlay once, is that a problem?)
    overlay = gaussian_filter(overlay, filter_sigma)
    overlay /= np.max(overlay)
    overlay = sigmoid(overlay, k=sigmoid_k, c=sigmoid_c)

    # Apply a colormap and delete the alpha channel
    overlay_cm = np.delete(matplotlib.cm.Blues(overlay.T), 3, axis=2)
    overlay_im = PIL.Image.fromarray(np.uint8(overlay_cm * 255))

    map_final = PIL.Image.blend(map_orig, overlay_im, 0.7)
    map_final.save(output_path)

def blur_grid(grid):
        filtered = gaussian_filter(grid, sigma=1)
        filtered[(grid > 0.45) & (grid < 0.55)] = grid[(grid > 0.45) &
                                                       (grid < 0.55)]
        return filtered

def debug_Data_Augmentation(blur=False, sigma=1.0, hflip=False, vflip=False, hvsplit=False, randbright=False):
    image = cv2.imread('Dataset/young/female/180.jpg', 0)
    #image = cv2.imread('Dataset/young/female/285.jpg', 0)
    #image = cv2.resize(image, (100, 100))
    cv2.imshow('Image', image)

    # Data Augmentation:
    # Gaussian Blurred 
    if blur:
        cv2.imshow('Blur', gaussian_filter(input=image, sigma=sigma))
        #cv2.imwrite("Blur_{:1.1f}.jpg".format(sigma), 
        #           gaussian_filter(input=image, sigma=sigma))
        cv2.imwrite("../xBlur_{:1.1f}.jpg".format(sigma), 
                    gaussian_filter(input=image, sigma=sigma))
    # Flip and Rotate
    if (hflip and not vflip) or (hflip and hvsplit):
        cv2.imshow('hflip', np.fliplr(image))
        cv2.imwrite("../hflip.jpg", np.fliplr(image))
    if (vflip and not hflip) or (vflip and hvsplit):
        cv2.imshow('vflip', np.flipud(image))
        cv2.imwrite("../vflip.jpg", np.flipud(image))
    if hflip and vflip and not hvsplit:
        cv2.imshow('rot 180', np.rot90(image, k=2))
        cv2.imwrite("../rot2k.jpg", np.rot90(image, k=2))
    cv2.waitKey(0)
    cv2.destroyAllWindows()

def debug_analyse_image_texture(file, sigma=1.0):
    image = cv2.imread(file, 0)
    blur  = gaussian_filter(input=image, sigma=sigma)
    cv2.imshow('Image', image - blur)
    #analysis = ndimage.gaussian_gradient_magnitude(image, sigma=sigma)
    #cv2.imshow('Analysis', analysis * 10)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
##########################################

def distort_elastic(image, smooth=10.0, scale=100.0, seed=0):
    """
    Elastic distortion of images.

    Channel axis in RGB images will not be distorted but grayscale or
    RGB images are both valid inputs. RGB and grayscale images will be
    distorted identically for the same seed.

    Simard, et. al, "Best Practices for Convolutional Neural Networks
    applied to Visual Document Analysis",
    in Proc. of the International Conference on Document Analysis and
    Recognition, 2003.

    :param ndarayy image: Image of shape [h,w] or [h,w,c]
    :param float smooth: Smoothes the distortion.
    :param float scale: Scales the distortion.
    :param int seed: Seed for random number generator. Ensures that for the
      same seed images are distorted identically.
    :return: Distorted image with same shape as input image.
    :rtype: ndarray
    """
    # create random, smoothed displacement field
    rnd = np.random.RandomState(int(seed))
    h, w = image.shape[:2]
    dxy = rnd.rand(2, h, w, 3) * 2 - 1
    dxy = gaussian_filter(dxy, smooth, mode="constant")
    dxy = dxy / np.linalg.norm(dxy) * scale
    dxyz = dxy[0], dxy[1], np.zeros_like(dxy[0])

    # create transformation coordinates and deform image
    is_color = len(image.shape) == 3
    ranges = [np.arange(d) for d in image.shape]
    grid = np.meshgrid(*ranges, indexing='ij')
    add = lambda v, dv: v + dv if is_color else v + dv[:, :, 0]
    idx = [np.reshape(add(v, dv), (-1, 1)) for v, dv in zip(grid, dxyz)]
    distorted = map_coordinates(image, idx, order=1, mode='reflect')

    return distorted.reshape(image.shape)

def level_curves(fname, npoints, smoothing = 10, level = 0.5) :
    # Find the contour lines
    img = misc.imread(fname, flatten = True) # Grayscale
    img = img.T[:, ::-1]
    img = img / 255.
    img = gaussian_filter(img, smoothing, mode='nearest')
    lines = find_contours(img, level)

    # Compute the sampling ratio
    lengths = []
    for line in lines :
        lengths.append( arclength(line) )
    lengths = array(lengths)
    points_per_line = ceil( npoints * lengths / sum(lengths) )

    # Interpolate accordingly
    points = []
    connec = []
    index_offset = 0
    for ppl, line in zip(points_per_line, lines) :
        (p, c) = resample(line, ppl)
        points.append(p)
        connec.append(c + index_offset)
        index_offset += len(p)

    points = vstack(points)
    connec = vstack(connec)
    return Curve(points.ravel(), connec, 2) # Dimension 2 !

def limitedfluxdensity(simulation, instrumentname, wavelengths, cmask, filterobject, fwhm, fluxlimit):

    # get the path for the data cube corresponding to this instrument
    fitspaths = filter(lambda fn: ("_"+instrumentname+"_") in fn, simulation.totalfitspaths())
    if len(fitspaths) != 1: return None

    # get the data cube and convert it to per-wavelength units
    cube = pyfits.getdata(arch.openbinary(fitspaths[0])).T
    cube = simulation.convert(cube, to_unit='W/m3/sr', quantity='surfacebrightness', wavelength=wavelengths)

    # convolve the data cube to a single frame, and convert back to per-frequency units
    frame = filterobject.convolve(wavelengths[cmask], cube[:,:,cmask])
    frame = simulation.convert(frame, from_unit='W/m3/sr', to_unit='MJy/sr', wavelength=filterobject.pivotwavelength())

    # get information on the simulated pixels (assume all frames are identical and square)
    sim_pixels = simulation.instrumentshape()[0]
    sim_pixelarea = simulation.angularpixelarea()   # in sr
    sim_pixelwidth = np.sqrt(sim_pixelarea) * 648000 / np.pi   # in arcsec

    # convolve the frame with a Gaussian of the appropriate size
    frame = gaussian_filter(frame, sigma=fwhm/sim_pixelwidth/2.35482, mode='constant')

    # get information on the observational instrument's pixels (assume pixel width ~ fwhm/3)
    obs_pixelwidth = fwhm/3   # in arcsec

    # calculate the bin size to obtain simulated pixels similar to observed pixels
    bin_pixels = find_nearest_divisor(sim_pixels, obs_pixelwidth/sim_pixelwidth)
    bin_pixelarea = sim_pixelarea * bin_pixels**2

    # rebin the frame
    frame = frame.reshape((sim_pixels//bin_pixels,bin_pixels,sim_pixels//bin_pixels,bin_pixels)).mean(axis=3).mean(axis=1)

    # integrate over the frame to obtain the total flux density and convert from MJy to Jy
    fluxdensity = frame[frame>fluxlimit].sum() * bin_pixelarea * 1e6
    return fluxdensity

# This function returns the integer divisor of the first (integer) argument that is nearest to the second (float) argument

def limitedfluxdensity(simulation, instrumentname, wavelengths, cmask, filterobject, fwhm, fluxlimit):

    # get the path for the data cube corresponding to this instrument
    fitspaths = filter(lambda fn: ("_"+instrumentname+"_") in fn, simulation.totalfitspaths())
    if len(fitspaths) != 1: return None

    # get the data cube and convert it to per-wavelength units
    cube = pyfits.getdata(arch.openbinary(fitspaths[0])).T
    cube = simulation.convert(cube, to_unit='W/m3/sr', quantity='surfacebrightness', wavelength=wavelengths)

    # convolve the data cube to a single frame, and convert back to per-frequency units
    frame = filterobject.convolve(wavelengths[cmask], cube[:,:,cmask])
    frame = simulation.convert(frame, from_unit='W/m3/sr', to_unit='MJy/sr', wavelength=filterobject.pivotwavelength())

    # get information on the simulated pixels (assume all frames are identical and square)
    sim_pixels = simulation.instrumentshape()[0]
    sim_pixelarea = simulation.angularpixelarea()   # in sr
    sim_pixelwidth = np.sqrt(sim_pixelarea) * 648000 / np.pi   # in arcsec

    # convolve the frame with a Gaussian of the appropriate size
    frame = gaussian_filter(frame, sigma=fwhm/sim_pixelwidth/2.35482, mode='constant')

    # get information on the observational instrument's pixels (assume pixel width ~ fwhm/3)
    obs_pixelwidth = fwhm/3   # in arcsec

    # calculate the bin size to obtain simulated pixels similar to observed pixels
    bin_pixels = find_nearest_divisor(sim_pixels, obs_pixelwidth/sim_pixelwidth)
    bin_pixelarea = sim_pixelarea * bin_pixels**2

    # rebin the frame
    frame = frame.reshape((sim_pixels//bin_pixels,bin_pixels,sim_pixels//bin_pixels,bin_pixels)).mean(axis=3).mean(axis=1)

    # integrate over the frame to obtain the total flux density and convert from MJy to Jy
    fluxdensity = frame[frame>fluxlimit].sum() * bin_pixelarea * 1e6
    return fluxdensity

# This function returns the integer divisor of the first (integer) argument that is nearest to the second (float) argument

def fastFilter(arr, ksize=30, every=None, resize=True, fn='median',
                  interpolation=cv2.INTER_LANCZOS4,
                  smoothksize=0,
                  borderMode=cv2.BORDER_REFLECT):
    '''
    fn['nanmean', 'mean', 'nanmedian', 'median']

    a fast 2d filter for large kernel sizes that also 
    works with nans
    the computation speed is increased because only 'every'nsth position 
    within the median kernel is evaluated
    '''
    if every is None:
        every = max(ksize//3, 1)
    else:
        assert ksize >= 3*every
    s0,s1 = arr.shape[:2]

    ss0 = s0//every
    every = s0//ss0
    ss1 = s1//every

    out = np.full((ss0+1,ss1+1), np.nan)

    c = {'median':_calcMedian,
         'nanmedian':_calcNanMedian,
         'nanmean':_calcNanMean,
         'mean':_calcMean,
         }[fn]
    ss0,ss1 = c(arr, out, ksize, every)
    out = out[:ss0,:ss1]

    if smoothksize:
        out = gaussian_filter(out, smoothksize)


    if not resize:
        return out
    return cv2.resize(out, arr.shape[:2][::-1],
               interpolation=interpolation)

def flatField(closeDist_img=None, inPlane_img=None,
              closeDist_bg=None, inPlane_bg=None,
              vignetting_model='different_objects',
              interpolation_method='kangWeiss',
              inPlane_scale_factor=None):

    # 1. Pixel sensitivity:
    if closeDist_img is not None:
        # TODO: find better name
        ff1 = flatFieldFromCalibration(closeDist_img, closeDist_bg)
    else:
        ff1 = 0

    # 2. Vignetting from in-plane measurements:
    if inPlane_img is not None:
        bg = gaussian_filter(median_filter(ff1, 3), 9)
        ff1 -= bg

        ff2, mask = VIGNETTING_MODELS[vignetting_model](inPlane_img,
                                                        inPlane_bg, inPlane_scale_factor)
#         import pylab as plt
#         plt.imshow(mask)
#         plt.show()
        ff2smooth = INTERPOLATION_METHODS[interpolation_method](ff2, mask)
        if isinstance(ff1, np.ndarray) and ff1.shape != ff2smooth.shape:
            ff2smooth = resize(ff2smooth, ff1.shape, mode='reflect')
    else:
        ff2 = 0
        ff2smooth = 0

    return ff1 + ff2smooth, ff2

def gauss_fltr(dem, sigma=1):
    print("Applying gaussian smoothing filter with sigma %s" % sigma)
    #Note, ndimage doesn't properly handle ma - convert to nan
    from scipy.ndimage.filters import gaussian_filter
    dem_filt_gauss = gaussian_filter(dem.filled(np.nan), sigma)
    #Now mask all nans
    #dem = np.ma.array(dem_filt_gauss, mask=dem.mask)
    out = np.ma.fix_invalid(dem_filt_gauss, copy=False, fill_value=dem.fill_value)
    out.set_fill_value(dem.fill_value)
    return out

def _gkern2(kernlen=21, nsig=3):
    """Returns a 2D Gaussian kernel array."""

    # create nxn zeros
    inp = np.zeros((kernlen, kernlen))
    # set element at the middle to one, a dirac delta
    inp[kernlen // 2, kernlen // 2] = 1
    # gaussian-smooth the dirac, resulting in a gaussian filter mask
    return fi.gaussian_filter(inp, nsig)

def elastic_transform(image, alpha, sigma, alpha_affine, random_state=None):
    """Elastic deformation of images as described in [Simard2003]_ (with modifications).
    .. [Simard2003] Simard, Steinkraus and Platt, "Best Practices for
         Convolutional Neural Networks applied to Visual Document Analysis", in
         Proc. of the International Conference on Document Analysis and
         Recognition, 2003.

     Based on https://gist.github.com/erniejunior/601cdf56d2b424757de5
    """
    if random_state is None:
        random_state = np.random.RandomState(None)

    shape = image.shape
    shape_size = shape[:2]

    # Random affine
    center_square = np.float32(shape_size) // 2
    square_size = min(shape_size) // 3
    pts1 = np.float32([center_square + square_size, [center_square[0]+square_size, center_square[1]-square_size], center_square - square_size])
    pts2 = pts1 + random_state.uniform(-alpha_affine, alpha_affine, size=pts1.shape).astype(np.float32)
    M = cv2.getAffineTransform(pts1, pts2)
    image = cv2.warpAffine(image, M, shape_size[::-1], borderMode=cv2.BORDER_REFLECT_101)

    dx = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma) * alpha
    dy = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma) * alpha
    dz = np.zeros_like(dx)

    x, y, z = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]), np.arange(shape[2]))
    indices = np.reshape(y+dy, (-1, 1)), np.reshape(x+dx, (-1, 1)), np.reshape(z, (-1, 1))

    return map_coordinates(image, indices, order=1, mode='reflect').reshape(shape)

def cylinder(self,deltaT=300.0,radius=200.0,start_depth=100.0):
      '''
      creates a cylindrical temperature anomaly with a gaussian blur

      usage: First read a temperature snapshot.

      Params:
      deltaT : cylinder excess temperature
      radius : cylinder radius (im km)
      start_depth: starting depth of the cylinder (default = 100.0)
      '''
      T_ref = self.T_adiabat[::-1]
      T_here = np.zeros((self.npts_rad,self.npts_theta))

      for i in range(0,self.npts_rad):
         print T_ref[i]
         km_per_degree = self.rad_km[i]*2*np.pi/360.0 

         for j in range(0,self.npts_theta):

            cyl_th_here  = radius / km_per_degree
            th_here      = self.theta[j]
            depth_here   = 6371.0 - self.rad_km[i]
            if th_here <= cyl_th_here and depth_here > 100.0:
               T_here[(self.npts_rad-1)-i,j] = T_ref[i] + deltaT
            else:
               T_here[(self.npts_rad-1)-i,j] = T_ref[i]

      filtered = gaussian_filter(T_here,sigma=[0,3])
      self.T = filtered
      #self.T = T_here

def blur(img, sigma=1):
    """Blur an image with a Gaussian kernel"""
    return filters.gaussian_filter(img, sigma)

def smooth_image(img, sigma = 1.0):
  if sigma is not None:
    return filters.gaussian_filter(np.asarray(img, float), sigma);
  else:
    return img;

def apply_laplacian_filter(array, alpha=30):
    """
    Laplacian is approximated with difference of Gaussian
    :param array:
    :param alpha:
    :return:
    """
    blurred_f = ndimage.gaussian_filter(array, 3)
    filter_blurred_f = ndimage.gaussian_filter(blurred_f, 1)
    return blurred_f + alpha * (blurred_f - filter_blurred_f)

def blur_image(self):
        H_blur = gaussian_filter(self.H,sigma=self.blur)
        H_blur = H_blur/np.max(H_blur)
        self.H_blur = H_blur