def gaussian(g):
    ks = 5
    y, x = np.mgrid[-ks:ks, -ks:ks]
    h = np.exp(-(x ** 2 + y ** 2) / (2. * (ks / 2) ** 2))
    h = h / h.sum()
    g_conv = fftconvolve(g, h[:, :, np.newaxis], mode='same')
    return g_conv
              
            

              
# generate the gauss kernel
def gauss_kernel(sigma):
    ks = np.round(2*sigma)+1
    y,x = np.mgrid[-ks:ks,-ks:ks]
    h = np.exp(-(x**2+y**2)/(2.*(ks/2)**2))
    h = h/h.sum()
    return h

def fast_bilateral_gaussian_gray(g, sigma_s, sigma_r, K):
    M, N = g.shape
    gk = gauss_kernel(sigma_s)

    # discrete gray values
    U_k = np.zeros((M, N, K))

    # pre-filter images
    for k in range(0, K):
        w_k = np.exp(-(k / (K - 1) - g) ** 2 / (2 * sigma_r ** 2))
        g_k = fftconvolve(g * w_k, gk, mode='same')
        W_k = fftconvolve(w_k, gk, mode='same')
        U_k[:, :, k] = g_k / np.maximum(1e-06, W_k)

    # indices for linear interpolation
    l = np.int32(np.floor(g * (K - 1)))

    # delta l
    dl = (g * (K - 1) - l) / (K - 1)

    U = np.zeros_like(g)
    for k in range(0, K):
        u = U_k[:, :, k]
        if k < K - 1:
            up = U_k[:, :, k + 1]
        else:
            up = U_k[:, :, k]
        U[l == k] = u[l == k] + dl[l == k] * (up[l == k] - u[l == k])

    return U
              
            

              
from skimage.restoration import denoise_bilateral

def denoise_bilateral(noisy_image, sigma_s, sigma_r):
    return denoise_bilateral(g,
                             sigma_color=sigma_s,
                             sigma_spatial=sigma_r,
                             mode="symmetric",
                             channel_axis=-1)