def test_denoise_tv_bregman_3d():
    img = checkerboard.copy()
    # add some random noise
    img += 0.5 * img.std() * np.random.rand(*img.shape)
    img = np.clip(img, 0, 1)

    out1 = restoration.denoise_tv_bregman(img, weight=10)
    out2 = restoration.denoise_tv_bregman(img, weight=5)

    # make sure noise is reduced in the checkerboard cells
    assert_(img[30:45, 5:15].std() > out1[30:45, 5:15].std())
    assert_(out1[30:45, 5:15].std() > out2[30:45, 5:15].std())

def test_denoise_tv_bregman_3d():
    img = lena
    # add some random noise
    img += 0.5 * img.std() * np.random.random(img.shape)
    img = np.clip(img, 0, 1)

    out1 = restoration.denoise_tv_bregman(img, weight=10)
    out2 = restoration.denoise_tv_bregman(img, weight=5)

    # make sure noise is reduced
    assert img.std() > out1.std()
    assert out1.std() > out2.std()

def blur_predict(model, X, type="median", filter_size=3, sigma=1.0):
  
    if type == "median":
        blured_X = np.array(list(map(lambda x: ndimage.median_filter(x, filter_size), 
                                     X)))
    elif type == "gaussian":
        blured_X = np.array(list(map(lambda x: ndimage.gaussian_filter(x, filter_size),
                                     X)))
    elif type == "f_gaussian":
        blured_X = np.array(list(map(lambda x: filters.gaussian_filter(x.reshape((28, 28)), sigma=sigma).reshape(784),
                                     X))) 
    elif type == "tv_chambolle":
        blured_X = np.array(list(map(lambda x: restoration.denoise_tv_chambolle(x.reshape((28, 28)), weight=0.2).reshape(784),
                                     X)))
    elif type == "tv_bregman":
        blured_X = np.array(list(map(lambda x: restoration.denoise_tv_bregman(x.reshape((28, 28)), weight=5.0).reshape(784),
                                     X)))
    elif type == "bilateral":
        blured_X = np.array(list(map(lambda x: restoration.denoise_bilateral(np.abs(x).reshape((28, 28))).reshape(784),
                                     X)))
    elif type == "nl_means":
        blured_X = np.array(list(map(lambda x: restoration.nl_means_denoising(x.reshape((28, 28))).reshape(784),
                                     X)))
        
    elif type == "none":
        blured_X = X 

    else:
        raise ValueError("unsupported filter type", type)

    return predict(model, blured_X)

    def __call__(self, x0, rho):
        """
        Proximal operator for the total variation denoising penalty

        Requires scikit-image be installed

        Parameters
        ----------
        x0 : array_like
            The starting or initial point used in the proximal update step

        rho : float
            Momentum parameter for the proximal step (larger value -> stays closer to x0)

        Raises
        ------
        ImportError
            If scikit-image fails to be imported

        """
        try:
            from skimage.restoration import denoise_tv_bregman
        except ImportError:
            print('Error: scikit-image not found. TVD will not work.')
            return x0

        return denoise_tv_bregman(x0, rho / self.penalty)

def tvd(x0, rho, gamma):
    """
    Proximal operator for the total variation denoising penalty

    Requires scikit-image be installed

    Parameters
    ----------
    x0 : array_like
        The starting or initial point used in the proximal update step

    rho : float
        Momentum parameter for the proximal step (larger value -> stays closer to x0)

    gamma : float
        A constant that weights how strongly to enforce the constraint

    Returns
    -------
    theta : array_like
        The parameter vector found after running the proximal update step

    Raises
    ------
    ImportError
        If scikit-image fails to be imported

    """

    return denoise_tv_bregman(x0, rho / gamma)

def make_prediction():
    import sys
    import numpy as np
    import pandas as pd

    from skimage.data import imread
    from skimage.filters import threshold_adaptive
    from skimage.restoration import denoise_tv_bregman

    from sklearn.cross_validation import train_test_split
    from sklearn.grid_search import GridSearchCV

    from sklearn.neighbors import KNeighborsClassifier
    from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier

    from model_design import model_design
    classifier = model_design()

    X, IDs = [], range(6284, 12504)
    for ID in IDs:
        original = imread('../data/testResized/' + str(ID) +'.Bmp', as_grey=True)
        denoised = denoise_tv_bregman(original, 3)
        binarilized = threshold_adaptive(denoised, block_size=13, method='gaussian')
        feature = binarilized.reshape(1,400)[0]
        X.append(feature)
    X = np.array(X)

    y = classifier.predict(X)
    result = pd.DataFrame({'Id': IDs, 'Class': y})
    result.to_csv('../result/06-09-2015_AdaBoostXTC.csv', sep=',', index=None, columns=['Id', 'Class'])

 def filter_frame(self, data):
     logging.debug("Running Denoise")
     weight = self.parameters['weight']
     max_iter = self.parameters['max_iterations']
     eps = self.parameters['error_threshold']
     isotropic = self.parameters['isotropic']
     return denoise_tv_bregman(data[0, ...], weight, max_iter=max_iter,
                               eps=eps, isotropic=isotropic)

def refine_worm(image, initial_area, candidate_edges):
    # find strong worm edges (roughly equivalent to the edges found by find_initial_worm,
    # which are in candidate_edges): smooth the image, do canny edge-finding, and
    # then keep only those edges near candidate_edges
    smooth_image = restoration.denoise_tv_bregman(image, 140).astype(numpy.float32)
    smoothed, gradient, sobel = canny.prepare_canny(smooth_image, 8, initial_area)
    local_maxima = canny.canny_local_maxima(gradient, sobel)
    candidate_edge_region = ndimage.binary_dilation(candidate_edges, iterations=4)
    strong_edges = local_maxima & candidate_edge_region

    # Now threshold the image to find dark blobs as our initial worm region
    # First, find areas in the initial region unlikely to be worm pixels
    mean, std = mcd.robust_mean_std(smooth_image[initial_area][::4], 0.85)
    non_worm = (smooth_image > mean - std) & initial_area
    # now fit a smoothly varying polynomial to the non-worm pixels in the initial
    # region of interest, and subtract that from the actual image to generate
    # an image with a flat illumination field
    background = polyfit.fit_polynomial(smooth_image, mask=non_worm, degree=2)
    minus_bg = smooth_image - background
    # now recalculate a threshold from the background-subtracted pixels
    mean, std = mcd.robust_mean_std(minus_bg[initial_area][::4], 0.85)
    initial_worm = (minus_bg < mean - std) & initial_area
    # Add any pixels near the strong edges to our candidate worm position
    initial_worm |= ndimage.binary_dilation(strong_edges, iterations=3)
    initial_worm = mask.fill_small_radius_holes(initial_worm, 5)

    # Now grow/shrink the initial_worm region so that as many of the strong
    # edges from the canny filter are in contact with the region edges as possible.
    ac = active_contour.EdgeClaimingAdvection(initial_worm, strong_edges,
        max_region_mask=initial_area)
    stopper = active_contour.StoppingCondition(ac, max_iterations=100)
    while stopper.should_continue():
        ac.advect(iters=1)
        ac.smooth(iters=1, depth=2)
    worm_mask = mask.fill_small_radius_holes(ac.mask, 7)

    # Now, get edges from the image at a finer scale
    smoothed, gradient, sobel = canny.prepare_canny(smooth_image, 0.3, initial_area)
    local_maxima = canny.canny_local_maxima(gradient, sobel)
    strong_sum = strong_edges.sum()
    highp = 100 * (1 - 1.5*strong_sum/local_maxima.sum())
    lowp = max(100 * (1 - 3*strong_sum/local_maxima.sum()), 0)
    low_worm, high_worm = numpy.percentile(gradient[local_maxima], [lowp, highp])
    fine_edges = canny.canny_hysteresis(local_maxima, gradient, low_worm, high_worm)

    # Expand out the identified worm area to include any of these finer edges
    closed_edges = ndimage.binary_closing(fine_edges, structure=S)
    worm = ndimage.binary_propagation(worm_mask, mask=worm_mask|closed_edges, structure=S)
    worm = ndimage.binary_closing(worm, structure=S, iterations=2)
    worm = mask.fill_small_radius_holes(worm, 5)
    worm = ndimage.binary_opening(worm)
    worm = mask.get_largest_object(worm)
    # Last, smooth the shape a bit to reduce sharp corners, but not too much to
    # sand off the tail
    ac = active_contour.CurvatureMorphology(worm, max_region_mask=initial_area)
    ac.smooth(depth=2, iters=2)
    return strong_edges, ac.mask

def denoiseTV_Bregman(imagen,isotropic):
    """
    -isotropic es el atributo para cambiar entre filtrado isotropico y anisotropico
    """
    noisy = img_as_float(imagen)

    denoise = denoise_tv_bregman(noisy, 7, 9, 0.08, isotropic)

    return denoise

def tvd(x, rho, penalty):
    """
    Total variation denoising proximal operator

    Parameters
    ----------
    penalty : float
    """

    return denoise_tv_bregman(x, rho / penalty)

def predict_data():
    X, IDs = [], range(6284, 12504)
    for ID in IDs:
        original = imread('../data/testResized/' + str(ID) +'.Bmp', as_grey=True)
        denoised = denoise_tv_bregman(original, 3)
        binarilized = threshold_adaptive(denoised, block_size=13, method='gaussian')
        feature = binarilized.reshape(1,400)[0]
        X.append(feature)
    X = np.array(X)
    return X

def test_denoise_tv_bregman_float_result_range():
    # lena image
    img = lena_gray
    int_lena = np.multiply(img, 255).astype(np.uint8)
    assert np.max(int_lena) > 1
    denoised_int_lena = restoration.denoise_tv_bregman(int_lena, weight=60.0)
    # test if the value range of output float data is within [0.0:1.0]
    assert denoised_int_lena.dtype == np.float
    assert np.max(denoised_int_lena) <= 1.0
    assert np.min(denoised_int_lena) >= 0.0

def test_denoise_tv_bregman_float_result_range():
    # astronaut image
    img = astro_gray.copy()
    int_astro = np.multiply(img, 255).astype(np.uint8)
    assert_(np.max(int_astro) > 1)
    denoised_int_astro = restoration.denoise_tv_bregman(int_astro, weight=60.0)
    # test if the value range of output float data is within [0.0:1.0]
    assert_(denoised_int_astro.dtype == np.float)
    assert_(np.max(denoised_int_astro) <= 1.0)
    assert_(np.min(denoised_int_astro) >= 0.0)

 def filter_frames(self, data):
     data = data[0]
     logging.debug("Running Denoise")
     weight = self.parameters['weight']
     max_iter = self.parameters['max_iterations']
     eps = self.parameters['error_threshold']
     isotropic = self.parameters['isotropic']
     data = np.nan_to_num(data[0, ...])
     result = denoise_tv_bregman(data, weight, max_iter=max_iter,
                                 eps=eps, isotropic=isotropic)
     return result

def make_step(net, xy, step_size=1.5, end='fc8', clip=True, unit=None, denoise_weight=0.1, margin=0, w=224, h=224):
    '''Basic gradient ascent step.'''

    src = net.blobs['data'] # input image is stored in Net's 'data' blob
    
    dst = net.blobs[end]
    acts = net.forward(end=end)

    if end in fc_layers:
        fc = acts[end][0]
        best_unit = fc.argmax()
        best_act = fc[best_unit]
        obj_act = fc[unit]
        # print "unit: %s [%.2f], obj: %s [%.2f]" % (best_unit, fc[best_unit], unit, obj_act)
    
    one_hot = np.zeros_like(dst.data)
    
    if end in fc_layers:
      one_hot.flat[unit] = 1.
    elif end in conv_layers:
      one_hot[:, unit, xy, xy] = 1.
    else:
      raise Exception("Invalid layer type!")

    dst.diff[:] = one_hot

    net.backward(start=end)
    g = src.diff[0]

    # Mask out gradient to limit the drawing region
    if margin != 0:
      mask = np.zeros_like(g)

      for dx in range(0 + margin, w - margin):
        for dy in range(0 + margin, h - margin):
          mask[:, dx, dy] = 1
      g *= mask

    src.data[:] += step_size/np.abs(g).mean() * g

    if clip:
        bias = net.transformer.mean['data']
        src.data[:] = np.clip(src.data, -bias, 255-bias) 

    # Run a separate TV denoising process on the resultant image
    asimg = deprocess( net, src.data[0] ).astype(np.float64)
    denoised = denoise_tv_bregman(asimg, weight=denoise_weight, max_iter=100, eps=1e-3)

    src.data[0] = preprocess( net, denoised )
    
    # reset objective for next step
    dst.diff.fill(0.)

    return best_unit, best_act, obj_act

def preprocess():
    labels = pd.read_csv('../data/trainLabels.csv', sep=',')
    X, y = [], np.array(labels.Class)

    for ID in labels.ID:
        original = imread('../data/trainResized/' + str(ID) +'.Bmp', as_grey=True)
        denoised = denoise_tv_bregman(original, 3)
        binarilized = threshold_adaptive(denoised, block_size=13, method='gaussian')
        feature = binarilized.reshape(1,400)[0]
        X.append(feature)
    X = np.array(X)
    return X, y

def model_design(run_as_main=False):

    from skimage.data import imread
    from skimage.filters import threshold_adaptive
    from skimage.restoration import denoise_tv_bregman

    from sklearn.cross_validation import train_test_split, StratifiedKFold

    from sklearn.svm import SVC
    from sklearn.linear_model import LogisticRegression
    from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
    from sklearn.ensemble import AdaBoostClassifier, BaggingClassifier

    labels = pd.read_csv('../data/trainLabels.csv', sep=',')
    X, y = [], np.array(labels.Class)


    for ID in labels.ID:
        original = imread('../data/trainResized/' + str(ID) +'.Bmp', as_grey=True)
        denoised = denoise_tv_bregman(original, 3)
        binarilized = threshold_adaptive(denoised, block_size=13, method='gaussian')
        feature = binarilized.reshape(1,400)[0]
        X.append(feature)
    X = np.array(X)

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

    clf = AdaBoostClassifier(base_estimator=
                             ExtraTreesClassifier(
                                 n_estimators=500,
                                 criterion='entropy',
                                 class_weight='auto',
                                 n_jobs=-1
                             ), n_estimators=50)

    # clf = AdaBoostClassifier(base_estimator=
    #                          RandomForestClassifier(
    #                              n_estimators=500,
    #                              criterion='entropy',
    #                              class_weight='auto',
    #                              n_jobs=-1
    #                          ), n_estimators=20)

    clf.fit(X_train, y_train)
    print clf.score(X_test, y_test)

cameraman_uni = sk.img_as_float(io.imread('../../cameraman_uni.bmp'))
cameraman_nor = sk.img_as_float(io.imread('../../cameraman_nor.bmp'))
cameraman_ric = sk.img_as_float(io.imread('../../cameraman_ric.bmp'))
cameraman_sp = sk.img_as_float(io.imread('../../cameraman_sp.bmp'))

baboon = sk.img_as_float(io.imread('../../baboon.bmp'))
baboon_uni = sk.img_as_float(io.imread('../../baboon_uni.bmp'))
baboon_nor = sk.img_as_float(io.imread('../../baboon_nor.bmp'))
baboon_ric = sk.img_as_float(io.imread('../../baboon_ric.bmp'))
baboon_sp = sk.img_as_float(io.imread('../../baboon_sp.bmp'))

# Define weight
w = 12

# Apply filters
lena_uni_f = sk.img_as_float(res.denoise_tv_bregman(lena_uni,12))
lena_nor_f = sk.img_as_float(res.denoise_tv_bregman(lena_nor,12))
lena_ric_f = sk.img_as_float(res.denoise_tv_bregman(lena_ric,12))
lena_sp_f = sk.img_as_float(res.denoise_tv_bregman(lena_sp,12))

cameraman_uni_f = sk.img_as_float(res.denoise_tv_bregman(cameraman_uni,12))
cameraman_nor_f = sk.img_as_float(res.denoise_tv_bregman(cameraman_nor,12))
cameraman_ric_f = sk.img_as_float(res.denoise_tv_bregman(cameraman_ric,12))
cameraman_sp_f = sk.img_as_float(res.denoise_tv_bregman(cameraman_sp,12))

baboon_uni_f = sk.img_as_float(res.denoise_tv_bregman(baboon_uni,12))
baboon_nor_f = sk.img_as_float(res.denoise_tv_bregman(baboon_nor,12))
baboon_ric_f = sk.img_as_float(res.denoise_tv_bregman(baboon_ric,12))
baboon_sp_f = sk.img_as_float(res.denoise_tv_bregman(baboon_sp,12))

# Calculate MSE

img = color.rgb2gray(data.astronaut())
img += 0.5 * img.std() * np.random.randn(*img.shape)
img.clip(0., 1., img)
# %%
tic = time()
denoise_wie = wiener(img)
print('Wiener {:.1f} sec'.format(time()-tic))
tic = time()
denoise_med = medfilt2d(img)
print('Median {:.1f} sec'.format(time()-tic))
tic = time()
denoise_tvc = skres.denoise_tv_chambolle(img)
print('TV Chambolle {:.1f} sec'.format(time()-tic))
tic = time()
denoise_tvb = skres.denoise_tv_bregman(img, weight=10.)
print('TV Bregman {:.1f} sec'.format(time()-tic))
tic = time()
denoise_bil = skres.denoise_bilateral(img, multichannel=False)
print('Bilateral {:.1f} sec'.format(time()-tic))
# %%
fg, axs = subplots(2, 3, sharex=True, sharey=True, figsize=(12, 10))
axs = axs.ravel()
axs[0].imshow(img, cmap='gray')
axs[1].imshow(denoise_wie, cmap='gray')
axs[2].imshow(denoise_med, cmap='gray')
axs[3].imshow(denoise_tvc, cmap='gray')
axs[4].imshow(denoise_tvb, cmap='gray')
axs[5].imshow(denoise_bil, cmap='gray')

axs[0].set_title('Original noisy')