def whiten(image):
    
    #return image
    
    tmp = image - image.mean()
    return tmp / numpy.std(tmp)
   
    if image.ndim > 3:
        raise TypeError('Not more than 3 dimensions supported')

    if image.ndim == 3:
        tmp = numpy.empty_like(image)
        for c in range(image.shape[2]):
            tmp[:, :, c] = whiten(image[:, :, c])

        result = numpy.zeros_like(image)
        for c1 in range(image.shape[2]):
            for c2 in range(image.shape[2]):
                if c1 == c2:
                    result[:, :, c1] += tmp[:, :, c2]
                else:
                    result[:, :, c1] -= tmp[:, :, c2]

        return result

    sigma1 = 0.2
    img1 = ndimage.gaussian_filter(image, sigma1)
    sigma2 = 5 * sigma1
    img2 = ndimage.gaussian_filter(image, sigma2)
    result = img1 - img2
    return result

def mkData():
    global data, cache, ui
    dtype = (ui.dtypeCombo.currentText(), ui.rgbCheck.isChecked())
    if dtype not in cache:
        if dtype[0] == 'uint8':
            dt = np.uint8
            loc = 128
            scale = 64
            mx = 255
        elif dtype[0] == 'uint16':
            dt = np.uint16
            loc = 4096
            scale = 1024
            mx = 2**16
        elif dtype[0] == 'float':
            dt = np.float
            loc = 1.0
            scale = 0.1
        
        if ui.rgbCheck.isChecked():
            data = np.random.normal(size=(20,512,512,3), loc=loc, scale=scale)
            data = ndi.gaussian_filter(data, (0, 6, 6, 0))
        else:
            data = np.random.normal(size=(20,512,512), loc=loc, scale=scale)
            data = ndi.gaussian_filter(data, (0, 6, 6))
        if dtype[0] != 'float':
            data = np.clip(data, 0, mx)
        data = data.astype(dt)
        cache[dtype] = data
        
    data = cache[dtype]
    updateLUT()

def test_gaussian_filter():
    # Test gaussian filter with np.float16
    # gh-8207
    data = np.array([1],dtype = np.float16)
    sigma = 1.0
    with assert_raises(RuntimeError):
        sndi.gaussian_filter(data,sigma)

def smooth_corrected(z, mask, w_smooth):
     mask1=snd.gaussian_filter(np.float32(mask), w_smooth, mode="constant", cval=0)
     ztemp=np.nan_to_num(z)
     ztemp[mask==0]=0.0
     zs=snd.gaussian_filter(ztemp, w_smooth, mode="constant", cval=0)
     zs[mask1>0]=zs[mask1>0]/mask1[mask1>0]
     return zs, mask1

    def _process(self, img):
        for c in xrange(3):
            channel = img[:, :, c]
            gaussian_filter(channel, output=channel, sigma=self._sigma)
            img[:, :, c] = channel

        return img

def density_to_sim_data(density, out=None):
    """
    Takes a 2D image input, returns a 4D SIM dataset output
    """
    if out == None:
        sim_data = np.zeros(
            (num_rotations, num_phases) + density.shape, dtype=np.float64)
    else:
        sim_data = out

    """
    Construct the illumination pattern
    """
    illumination = generate_illumination(density)
    sim_data_to_visualization(illumination, 'illumination.tif')

    """
    Simulate the imaging process: multiply by the illumination, and blur
    """
    for t in range(num_rotations):
        for p in range(num_phases):
            sim_data[t, p, :, :] = illumination[t, p, :, :] * density
            gaussian_filter(sim_data[t, p, :, :],
                            sigma=emission_sigma,
                            output=sim_data[t, p, :, :])
    return sim_data

 def loss(self, H): 
     if not hasattr(self, 'sigma'): 
         self.sigma = .05*np.array(H.shape)
     npts = nonzero(self.npoints(H))
     Hs = H/npts
     gaussian_filter(Hs, sigma=self.sigma, mode='constant', output=Hs)
     return dist2loss(Hs)

def blur_image(img):
    '''Return the blurred image that's used when sampling'''
    blur = np.zeros(list(img.shape)+[2], img.dtype)
    for z in range(img.shape[2]):
        blur[:,:,z, 0] = laplace(gaussian_filter(img[:,:,z], 3))
        blur[:,:,z, 1] = gaussian_filter(img[:,:,z], 5)
    return blur

def compute_harris_response(image, eps=1e-6):
    """ compute the Harris corner detector response function
        for each pixel in the image"""

    # derivatives
    image = ndimage.gaussian_filter(image, 1)
    imx = ndimage.sobel(image, axis=0, mode='constant')
    imy = ndimage.sobel(image, axis=1, mode='constant')

    Wxx = ndimage.gaussian_filter(imx * imx, 1.5, mode='constant')
    Wxy = ndimage.gaussian_filter(imx * imy, 1.5, mode='constant')
    Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant')

    # determinant and trace
    Wdet = Wxx * Wyy - Wxy ** 2
    Wtr = Wxx + Wyy
    harris = Wdet / (Wtr + eps)

    # Non maximum filter of size 3
    harris_max = ndimage.maximum_filter(harris, 3, mode='constant')
    harris *= harris == harris_max
    # Remove the image corners
    harris[:3] = 0
    harris[-3:] = 0
    harris[:, :3] = 0
    harris[:, -3:] = 0

    return harris

def generate_surface_from_stack(stack, sd=(10, 10, 10), surface_blur_sd=5):
    """Return a 2D image encoding a height map, generated from the input stack.
    The image is generated by first blurring the stack, then taking the z index
    of the brightest point for each X, Y location. The resultant surface is
    then smoothed with a gaussian filter.

    sd: standard deviation in each direction
    surface_blur_sd: standard deviation of smoothing applied to 2D surface."""


    ydim, xdim, zdim = stack.shape
    pad_val = 5
    padding = pad_val * 2
    padded_stack = np.zeros((ydim + padding, xdim + padding, zdim + padding),
                            dtype=stack.dtype)
    start = pad_val
    end_add = pad_val
    padded_stack[start:ydim+end_add, start:xdim+end_add, start:zdim+end_add] = stack
    smoothed_stack = nd.gaussian_filter(padded_stack, sd)
    cropped_stack = smoothed_stack[start:ydim+end_add, start:xdim+end_add, start:zdim+end_add]
    raw_surface = np.argmax(cropped_stack, 2)
    raw_surface[np.logical_and(raw_surface==0, cropped_stack[:,:,0] == 0)] = zdim-1
    smoothed_surface = nd.gaussian_filter(raw_surface, surface_blur_sd)

    return smoothed_surface

def preproc_dog(bg):
    """test low level image filter
    """

    #special dependencies
    from scipy.ndimage import gaussian_filter
    from pyrankfilter import rankfilter,__version__

    bg = bg.astype(float)
    f1 = gaussian_filter(bg,sigma=14.5)
    f2 = gaussian_filter(bg,sigma=15)

    f = ((f1-f2+3)*10).astype('uint8')

    loc_max = rankfilter(f,'highest',20,infSup=[0,0])

    r = bg.copy()
    r[loc_max>0]=0

    plt.figure()
    plt.subplot(2,2,1)
    plt.imshow(bg)
    plt.subplot(2,2,2)
    plt.imshow(f)
    plt.colorbar()
    plt.subplot(2,2,3)
    plt.imshow(loc_max)
    plt.subplot(2,2,4)
    plt.imshow(r)
    plt.show()

    def tantriggs(self, x, alpha=0.1,gamma=0.2,sigma0=1,sigma1=2,tau=10.0):
        x = np.array(x, dtype=np.float32)
        x = np.power(x, gamma)
        s0 = 3*sigma0
        s1 = 3*sigma1
        if ((s0%2)==0):
            s0+=1
        if ((s1%2)==0):
            s1+=1

        x = np.asarray(
            ndimage.gaussian_filter(x, sigma0) - ndimage.gaussian_filter(x, sigma1)
            )

        x = x / np.power(
            np.mean(np.power(np.abs(x), alpha)),
            1.0 / alpha
            )
        x = x / np.power(
                np.mean(
                    np.power(
                        np.minimum(np.abs(x), tau),
                        alpha
                    )
                ),
                1.0 / alpha
            )

        x = np.tanh(x / tau) * tau
        x = cv2.normalize(x,x,-220,0,cv2.NORM_MINMAX)
        return np.array(x, np.uint8)

def nan_gaussian_filter(X, sigma, keep_nans=False, gf_kwargs=dict()):
    """Equivalent to scipy.ndimage.gaussian_filter, but allows NaNs

    For inputs see original function

    if keep_nans is True, then output has NaNs everywhere input had nans

    http://stackoverflow.com/questions/18697532/
    gaussian-filtering-a-image-with-nan-in-python"""
    nanX = np.isnan(X)

    X1 = X.copy()
    X1[nanX] = 0
    X1_filtered = gaussian_filter(X1, sigma, **gf_kwargs)

    X2 = np.ones_like(X)
    X2[nanX] = 0
    X2_filtered = gaussian_filter(X2, sigma, **gf_kwargs)

    out = X1_filtered/X2_filtered

    if keep_nans:
        out[nanX] = np.nan

    return out

def sharpen(Xs,ys,sigma1=3,sigma2=1,alpha=30,ratio=0.6,t=None):
    print "sharpening images..."

    if ratio > 1.0:
        print "Do you really want a ratio of %f for sharpen?" % ratio
        print "Every images produced will always be similar"

    rand_list = randomindex(Xs.shape[0]*ratio,Xs.shape[0])
    Xs = Xs[rand_list]
    ys = ys[rand_list]

    tx = []
    ty = []

    for X, y in zip(Xs,ys):
        blurred_l = ndimage.gaussian_filter(X, sigma1)

        filter_blurred_l = ndimage.gaussian_filter(blurred_l, sigma2)
        sharpened = blurred_l + alpha * (blurred_l - filter_blurred_l)

        tx.append(sharpened)
        ty.append(y)

        if t: t.print_update(1)

    return np.array(tx), np.array(ty)

def PlotProfile(label1, label2, data1, data2, filename, smooth, binsz):
    plt.clf()
    plt.figure(figsize=(5,4))
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif')
    
    counts1 = np.append(data1['FLUX'].data/(EXPOSURE*PIXEL_SA), 0)
    glats1 = np.append(data1['GLAT_MIN'][0], data1['GLAT_MAX'].data)

    profile1 = np.histogram(glats1, bins=np.sort(glats1), weights=counts1)
    xstep=binsz
    y_smooth_1 = gaussian_filter(profile1[0], smooth / xstep)
    print y_smooth_1.mean()
    
    counts2 = np.append(data2['FLUX'].data/(EXPOSURE*PIXEL_SA*(data2['FLUX'].data.size())), 0)
    glats2 = np.append(data2['GLAT_MIN'][0], data2['GLAT_MAX'].data)

    profile2 = np.histogram(glats2, bins=np.sort(glats2), weights=counts2)
    xstep=binsz
    y_smooth_2 = gaussian_filter(profile2[0], smooth / xstep)
    print y_smooth_2.mean()
    
    x1 = 0.5 * (glats1[1:] + glats1[:-1])
    plt.plot(x1, y_smooth_1, label='{0}'.format(label1))
    plt.plot(x1, y_smooth_2, label='{0}'.format(label2))
    plt.hlines(0, LatLow+binsz, LatHigh-binsz)
    plt.xlabel(r'Galactic Latitude/$deg$', fontsize=10)
    plt.ylabel(r'Surface Brightness/ph cm$^{-2}$ s$^{-1} sr^{-1}$', fontsize=10)
    plt.xlim([LatLow+binsz, LatHigh-binsz])
    plt.tick_params(axis='x', labelsize=10)
    plt.grid(b=True, which='major', color='0.75', linewidth=0.5)
    plt.legend(prop={'size':8})
    plt.savefig(filename)

    def _plot_path_energy(self, experiment_results):
        Tmax = 6000
        for repetition in experiment_results:
            nodes = experiment_results.nodes_have_metric("routeEnergy")
            for node in nodes:
                data = experiment_results.get_tuple_metric_per_node("routeEnergy", node, repetition)

                T = [float(pair[0]) for pair in data]
                R = [float(pair[1]) for pair in data]
                bw = 0.5

                trange = [0, Tmax]
                bins = 5000

                dx = (trange[1] - trange[0]) / bins

                # compute sum_i K(x - x_i) y_i
                hist_R, edges = np.histogram(T, range=trange, bins=bins, weights=R)
                kde_R = gaussian_filter(hist_R, bw / dx)

                # compute sum_i K(x - x_i)
                hist_T, edges = np.histogram(T, range=trange, bins=bins)
                kde_T = gaussian_filter(hist_T, bw / dx)

                # compute the Nadaraya-Watson estimate
                interpolated_R = kde_R / kde_T

                # computer x-axis
                domain = (edges[1:] + edges[:-1]) / 2.0

                if self.draw:
                    plt.plot(domain, interpolated_R)
                    file_name = "alternative_path-energy_node-" + str(node)
                    plt.savefig(os.path.join(self.location, self.scenario + "_" + file_name + ".png"))
                    plt.close()

 def __init__(self, location, fig, ax):
     self.location = location
     self.fig = fig
     self.ax = ax
     # get alignment data for given location #
     try:
         self.afm_data_exists = 1
         self.afm_data = afm_alignment_data(location)
     except IOError:
         self.afm_data_exists = 0
     if self.afm_data_exists:
         # smooth data #
         self.afm_data.amplitude = ndimage.gaussian_filter(self.afm_data.amplitude, 0.7)
         self.afm_data.phase = ndimage.gaussian_filter(self.afm_data.phase, 0.7)
         
         # make afm data subplots #
         gs = gridspec.GridSpecFromSubplotSpec(1, 2, self.ax.get_subplotspec(), wspace=0.3)
         r_ax = plt.subplot(gs[0])
         theta_ax = plt.subplot(gs[1])
         # plot on subplot axes #
         self.add_alignment_subplot(self.afm_data.amplitude, r_ax, x=None, y=None, name='$r$')
         self.add_alignment_subplot(self.afm_data.phase, theta_ax, x=None, y=None, name=r'$\theta$')
         for tl in theta_ax.get_yticklabels():
             tl.set_visible(False)
         r_ax.set_ylabel('y position (nm)')
         r_ax.set_xlabel('x position (nm)')

    def get_M(self, X, mode='reflect', cval=0.):
        '''Get the Harris-Laplace scale-adapted second moment matrix'''
        # Compute the gaussian smoothed image
        # N.B.: using ndimage for speed; could replace with gaussian filter
        # class above
        G = ndimage.gaussian_filter(X, self.sigma_d, mode=mode, cval=cval)

        # Compute derivatives of gaussian-smoothed image in x and y directions
        Lx, Ly = self.get_grads_xy(G, mode=mode, cval=cval)

        # Compute second derivatives from first derivatives in x and y
        # directions
        Lxx, Lxy = self.get_grads_xy(Lx, mode=mode, cval=cval)
        _, Lyy = self.get_grads_xy(Ly, mode=mode, cval=cval)

        # Convolve each second derivative matrix with the integration gaussian
        # N.B.: using ndimage for speed; could replace with gaussian filter
        # class above
        Lxx = ndimage.gaussian_filter(Lxx, self.sigma_i, mode=mode, cval=cval)
        Lxy = ndimage.gaussian_filter(Lxy, self.sigma_i, mode=mode, cval=cval)
        Lyy = ndimage.gaussian_filter(Lyy, self.sigma_i, mode=mode, cval=cval)

        # Get an empty matrix and store the second derivatives
        M = np.empty((*X.shape, 2, 2))
        M[:, :, 0, 0] = Lxx
        M[:, :, 1, 0] = Lxy
        M[:, :, 0, 1] = Lxy
        M[:, :, 1, 1] = Lyy

        # Apply scale correction
        M = (self.sigma_d ** 2.) * M
        return M

def cont(imagen, depth=2**16, gaussian=3, screenpercent=0.7,t=0):

    imagen = gaussian_filter(imagen, gaussian)

    if t==0:
        otsu = threshold_otsu(imagen, depth)
    elif t==1:
        otsu = filters.threshold_isodata(imagen, depth)
    else:
        otsu = filters.threshold_li(imagen)
    imagen = binarizar(imagen, otsu)
    imagen = gaussian_filter(imagen, gaussian)

    contours = measure.find_contours(imagen, 1)
    centro = np.asanyarray([1280*0.5, 960*0.5])
    while len(contours) > 1:
        if sum(np.abs(centro - contours[1].mean(axis=0)) < [1280*screenpercent*0.5, 960*screenpercent*0.5]) != 2:
            del contours[1]
        elif sum(np.abs(centro - contours[0].mean(axis=0)) < [1280*screenpercent*0.5, 960*screenpercent*0.5]) != 2:
            del contours[0]
        else:
            if contours[1].size < contours[0].size:
                del contours[1]
            else:
                del contours[0]
    return imagen, contours[0]

    def __sharpner(self, image_data):
        blurred_f = ndimage.gaussian_filter(image_data, 3)
        filter_blurred_f = ndimage.gaussian_filter(blurred_f, 1)
        alpha = 30
        image_sharpened = blurred_f + alpha * (blurred_f - filter_blurred_f)

        return image_sharpened