def test_wavelet_denoising():
    rstate = np.random.RandomState(1234)

    # version with one odd-sized dimension
    astro_gray_odd = astro_gray[:, :-1]
    astro_odd = astro[:, :-1]

    for img, multichannel in [(astro_gray, False), (astro_gray_odd, False),
                              (astro_odd, True)]:
        sigma = 0.1
        noisy = img + sigma * rstate.randn(*(img.shape))
        noisy = np.clip(noisy, 0, 1)

        # Verify that SNR is improved when true sigma is used
        denoised = restoration.denoise_wavelet(noisy, sigma=sigma,
                                               multichannel=multichannel)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert_(psnr_denoised > psnr_noisy)

        # Verify that SNR is improved with internally estimated sigma
        denoised = restoration.denoise_wavelet(noisy,
                                               multichannel=multichannel)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert_(psnr_denoised > psnr_noisy)

        # Test changing noise_std (higher threshold, so less energy in signal)
        res1 = restoration.denoise_wavelet(noisy, sigma=2*sigma,
                                           multichannel=multichannel)
        res2 = restoration.denoise_wavelet(noisy, sigma=sigma,
                                           multichannel=multichannel)
        assert_(np.sum(res1**2) <= np.sum(res2**2))

def test_wavelet_denoising_levels():
    rstate = np.random.RandomState(1234)
    ndim = 2
    N = 256
    wavelet = 'db1'
    # Generate a very simple test image
    img = 0.2*np.ones((N, )*ndim)
    img[[slice(5, 13), ] * ndim] = 0.8

    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    denoised = restoration.denoise_wavelet(noisy, wavelet=wavelet)
    denoised_1 = restoration.denoise_wavelet(noisy, wavelet=wavelet,
                                             wavelet_levels=1)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    psnr_denoised_1 = compare_psnr(img, denoised_1)

    # multi-level case should outperform single level case
    assert_(psnr_denoised > psnr_denoised_1 > psnr_noisy)

    # invalid number of wavelet levels results in a ValueError
    max_level = pywt.dwt_max_level(np.min(img.shape),
                                   pywt.Wavelet(wavelet).dec_len)
    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
                  wavelet=wavelet, wavelet_levels=max_level+1)
    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
                  wavelet=wavelet, wavelet_levels=-1)

def test_wavelet_threshold():
    rstate = np.random.RandomState(1234)

    img = astro_gray
    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    # employ a single, user-specified threshold instead of BayesShrink sigmas
    with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
        denoised = _wavelet_threshold(noisy, wavelet='db1', method=None,
                                      threshold=sigma)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    assert_(psnr_denoised > psnr_noisy)

    # either method or threshold must be defined
    with testing.raises(ValueError):
        _wavelet_threshold(noisy, wavelet='db1', method=None, threshold=None)

    # warns if a threshold is provided in a case where it would be ignored
    with expected_warnings(["Thresholding method ",
                            PYWAVELET_ND_INDEXING_WARNING]):
        _wavelet_threshold(noisy, wavelet='db1', method='BayesShrink',
                           threshold=sigma)

def test_wavelet_denoising():
    for img, multichannel in [(astro_gray, False), (astro, True)]:
        sigma = 0.1
        noisy = img + sigma * np.random.randn(*(img.shape))
        noisy = np.clip(noisy, 0, 1)

        # Verify that SNR is improved when true sigma is used
        denoised = restoration.denoise_wavelet(noisy, sigma=sigma,
                                               multichannel=multichannel)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert psnr_denoised > psnr_noisy

        # Verify that SNR is improved with internally estimated sigma
        denoised = restoration.denoise_wavelet(noisy,
                                               multichannel=multichannel)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert psnr_denoised > psnr_noisy

        # Test changing noise_std (higher threshold, so less energy in signal)
        res1 = restoration.denoise_wavelet(noisy, sigma=2*sigma,
                                           multichannel=multichannel)
        res2 = restoration.denoise_wavelet(noisy, sigma=sigma,
                                           multichannel=multichannel)
        assert (res1.sum()**2 <= res2.sum()**2)

def test_wavelet_denoising_nd():
    rstate = np.random.RandomState(1234)
    for method in ['VisuShrink', 'BayesShrink']:
        for ndim in range(1, 5):
            # Generate a very simple test image
            if ndim < 3:
                img = 0.2*np.ones((128, )*ndim)
            else:
                img = 0.2*np.ones((16, )*ndim)
            img[(slice(5, 13), ) * ndim] = 0.8

            sigma = 0.1
            noisy = img + sigma * rstate.randn(*(img.shape))
            noisy = np.clip(noisy, 0, 1)

            # Mark H. 2018.08:
            #   The issue arises because when ndim in [1, 2]
            #   ``waverecn`` calls ``_match_coeff_dims``
            #   Which includes a numpy 1.15 deprecation.
            #   for larger number of dimensions _match_coeff_dims isn't called
            #   for some reason.
            anticipated_warnings = (PYWAVELET_ND_INDEXING_WARNING
                                    if ndim < 3 else None)
            with expected_warnings([anticipated_warnings]):
                # Verify that SNR is improved with internally estimated sigma
                denoised = restoration.denoise_wavelet(noisy, method=method)
            psnr_noisy = compare_psnr(img, noisy)
            psnr_denoised = compare_psnr(img, denoised)
            assert_(psnr_denoised > psnr_noisy)

def test_PSNR_dynamic_range_and_data_range():
    # Tests deprecation of "dynamic_range" in favor of "data_range"
    out1 = compare_psnr(cam/255., cam_noisy/255., data_range=1)
    with expected_warnings(
            '`dynamic_range` has been deprecated in favor of '
            '`data_range`. The `dynamic_range` keyword argument '
            'will be removed in v0.14'):
        out2 = compare_psnr(cam/255., cam_noisy/255., dynamic_range=1)
    assert_equal(out1, out2)

def test_denoise_nl_means_multichannel():
    # for true 3D data, 3D denoising is better than denoising as 2D+channels
    img = np.zeros((13, 10, 8))
    img[6, 4:6, 2:-2] = 1.
    sigma = 0.3
    imgn = img + sigma * np.random.randn(*img.shape)
    denoised_wrong_multichannel = restoration.denoise_nl_means(
        imgn, 3, 4, 0.6 * sigma, fast_mode=True, multichannel=True)
    denoised_ok_multichannel = restoration.denoise_nl_means(
        imgn, 3, 4, 0.6 * sigma, fast_mode=True, multichannel=False)
    psnr_wrong = compare_psnr(img, denoised_wrong_multichannel)
    psnr_ok = compare_psnr(img, denoised_ok_multichannel)
    assert_(psnr_ok > psnr_wrong)

def test_wavelet_threshold():
    rstate = np.random.RandomState(1234)

    img = astro_gray
    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    # employ a single, uniform threshold instead of BayesShrink sigmas
    denoised = _wavelet_threshold(noisy, wavelet='db1', threshold=sigma)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    assert_(psnr_denoised > psnr_noisy)

def test_wavelet_denoising():
    rstate = np.random.RandomState(1234)

    # version with one odd-sized dimension
    astro_gray_odd = astro_gray[:, :-1]
    astro_odd = astro[:, :-1]

    for img, multichannel, convert2ycbcr in [(astro_gray, False, False),
                                             (astro_gray_odd, False, False),
                                             (astro_odd, True, False),
                                             (astro_odd, True, True)]:
        sigma = 0.1
        noisy = img + sigma * rstate.randn(*(img.shape))
        noisy = np.clip(noisy, 0, 1)

        # Verify that SNR is improved when true sigma is used
        with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
            denoised = restoration.denoise_wavelet(noisy, sigma=sigma,
                                                   multichannel=multichannel,
                                                   convert2ycbcr=convert2ycbcr)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert_(psnr_denoised > psnr_noisy)

        # Verify that SNR is improved with internally estimated sigma
        with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
            denoised = restoration.denoise_wavelet(noisy,
                                                   multichannel=multichannel,
                                                   convert2ycbcr=convert2ycbcr)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert_(psnr_denoised > psnr_noisy)

        # SNR is improved less with 1 wavelet level than with the default.
        denoised_1 = restoration.denoise_wavelet(noisy,
                                                 multichannel=multichannel,
                                                 wavelet_levels=1,
                                                 convert2ycbcr=convert2ycbcr)
        psnr_denoised_1 = compare_psnr(img, denoised_1)
        assert_(psnr_denoised > psnr_denoised_1)
        assert_(psnr_denoised_1 > psnr_noisy)

        # Test changing noise_std (higher threshold, so less energy in signal)
        with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
            res1 = restoration.denoise_wavelet(noisy, sigma=2 * sigma,
                                               multichannel=multichannel)
        with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
            res2 = restoration.denoise_wavelet(noisy, sigma=sigma,
                                               multichannel=multichannel)
        assert_(np.sum(res1**2) <= np.sum(res2**2))

def test_wavelet_denoising_nd():
    for ndim in range(1, 5):
        # Generate a very simple test image
        img = 0.2*np.ones((16, )*ndim)
        img[[slice(5, 13), ] * ndim] = 0.8

        sigma = 0.1
        noisy = img + sigma * np.random.randn(*(img.shape))
        noisy = np.clip(noisy, 0, 1)

        # Verify that SNR is improved with internally estimated sigma
        denoised = restoration.denoise_wavelet(noisy)
        psnr_noisy = compare_psnr(img, noisy)
        psnr_denoised = compare_psnr(img, denoised)
        assert psnr_denoised > psnr_noisy

def test_PSNR_float():
    p_uint8 = compare_psnr(cam, cam_noisy)
    p_float64 = compare_psnr(cam / 255., cam_noisy / 255.,
                             data_range=1)
    assert_almost_equal(p_uint8, p_float64, decimal=5)

    # mixed precision inputs
    p_mixed = compare_psnr(cam / 255., np.float32(cam_noisy / 255.),
                           data_range=1)
    assert_almost_equal(p_mixed, p_float64, decimal=5)

    # mismatched dtype results in a warning if data_range is unspecified
    with expected_warnings(['Inputs have mismatched dtype']):
        p_mixed = compare_psnr(cam / 255., np.float32(cam_noisy / 255.))
    assert_almost_equal(p_mixed, p_float64, decimal=5)

def test_denoise_nl_means_3d():
    img = np.zeros((12, 12, 8))
    img[5:-5, 5:-5, 2:-2] = 1.
    sigma = 0.3
    imgn = img + sigma * np.random.randn(*img.shape)
    psnr_noisy = compare_psnr(img, imgn)
    for s in [sigma, 0]:
        denoised = restoration.denoise_nl_means(imgn, 3, 4, h=0.75 * sigma,
                                                fast_mode=True,
                                                multichannel=False, sigma=s)
        # make sure noise is reduced
        assert_(compare_psnr(img, denoised) > psnr_noisy)
        denoised = restoration.denoise_nl_means(imgn, 3, 4, h=0.75 * sigma,
                                                fast_mode=False,
                                                multichannel=False, sigma=s)
        # make sure noise is reduced
        assert_(compare_psnr(img, denoised) > psnr_noisy)

def test_wavelet_denoising_levels():
    rstate = np.random.RandomState(1234)
    ndim = 2
    N = 256
    wavelet = 'db1'
    # Generate a very simple test image
    img = 0.2*np.ones((N, )*ndim)
    img[(slice(5, 13), ) * ndim] = 0.8

    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
        denoised = restoration.denoise_wavelet(noisy, wavelet=wavelet)
    denoised_1 = restoration.denoise_wavelet(noisy, wavelet=wavelet,
                                             wavelet_levels=1)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    psnr_denoised_1 = compare_psnr(img, denoised_1)

    # multi-level case should outperform single level case
    assert_(psnr_denoised > psnr_denoised_1 > psnr_noisy)

    # invalid number of wavelet levels results in a ValueError or UserWarning
    max_level = pywt.dwt_max_level(np.min(img.shape),
                                   pywt.Wavelet(wavelet).dec_len)
    if Version(pywt.__version__) < '1.0.0':
        # exceeding max_level raises a ValueError in PyWavelets 0.4-0.5.2
        with testing.raises(ValueError):
            with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
                restoration.denoise_wavelet(
                    noisy, wavelet=wavelet, wavelet_levels=max_level + 1)
    else:
        # exceeding max_level raises a UserWarning in PyWavelets >= 1.0.0
        with expected_warnings([
                'all coefficients will experience boundary effects']):
            restoration.denoise_wavelet(
                noisy, wavelet=wavelet, wavelet_levels=max_level + 1)

    with testing.raises(ValueError):
        with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
            restoration.denoise_wavelet(
                noisy,
                wavelet=wavelet, wavelet_levels=-1)

def test(model):
    
    print('Start to test on {}'.format(args.test_dir))
    out_dir = save_dir + args.test_dir.split('/')[-1] + '/'
    if not os.path.exists(out_dir):
            os.mkdir(out_dir)
            
    name = []
    psnr = []
    ssim = []
    file_list = glob.glob('{}/*.png'.format(args.test_dir))
    for file in file_list:
        # read image
        img_clean = np.array(Image.open(file), dtype='float32') / 255.0
        img_test = img_clean + np.random.normal(0, args.sigma/255.0, img_clean.shape)
        img_test = img_test.astype('float32')
        # predict
        x_test = img_test.reshape(1, img_test.shape[0], img_test.shape[1], 1) 
        y_predict = model.predict(x_test)
        # calculate numeric metrics
        img_out = y_predict.reshape(img_clean.shape)
        img_out = np.clip(img_out, 0, 1)
        psnr_noise, psnr_denoised = compare_psnr(img_clean, img_test), compare_psnr(img_clean, img_out)
        ssim_noise, ssim_denoised = compare_ssim(img_clean, img_test), compare_ssim(img_clean, img_out)
        psnr.append(psnr_denoised)
        ssim.append(ssim_denoised)
        # save images
        filename = file.split('/')[-1].split('.')[0]    # get the name of image file
        name.append(filename)
        img_test = Image.fromarray((img_test*255).astype('uint8'))
        img_test.save(out_dir+filename+'_sigma'+'{}_psnr{:.2f}.png'.format(args.sigma, psnr_noise))
        img_out = Image.fromarray((img_out*255).astype('uint8')) 
        img_out.save(out_dir+filename+'_psnr{:.2f}.png'.format(psnr_denoised))
    
    psnr_avg = sum(psnr)/len(psnr)
    ssim_avg = sum(ssim)/len(ssim)
    name.append('Average')
    psnr.append(psnr_avg)
    ssim.append(ssim_avg)
    print('Average PSNR = {0:.2f}, SSIM = {1:.2f}'.format(psnr_avg, ssim_avg))
    
    pd.DataFrame({'name':np.array(name), 'psnr':np.array(psnr), 'ssim':np.array(ssim)}).to_csv(out_dir+'/metrics.csv', index=True)

def test_wavelet_denoising_nd():
    rstate = np.random.RandomState(1234)
    for method in ['VisuShrink', 'BayesShrink']:
        for ndim in range(1, 5):
            # Generate a very simple test image
            if ndim < 3:
                img = 0.2*np.ones((128, )*ndim)
            else:
                img = 0.2*np.ones((16, )*ndim)
            img[[slice(5, 13), ] * ndim] = 0.8

            sigma = 0.1
            noisy = img + sigma * rstate.randn(*(img.shape))
            noisy = np.clip(noisy, 0, 1)

            # Verify that SNR is improved with internally estimated sigma
            denoised = restoration.denoise_wavelet(noisy, method=method)
            psnr_noisy = compare_psnr(img, noisy)
            psnr_denoised = compare_psnr(img, denoised)
            assert_(psnr_denoised > psnr_noisy)

def test_denoise_nl_means_2d_multichannel():
    # reduce image size because nl means is slow
    img = np.copy(astro[:50, :50])
    img = np.concatenate((img, ) * 2, )  # 6 channels

    # add some random noise
    sigma = 0.1
    imgn = img + sigma * np.random.standard_normal(img.shape)
    imgn = np.clip(imgn, 0, 1)
    for fast_mode in [True, False]:
        for s in [sigma, 0]:
            for n_channels in [2, 3, 6]:
                psnr_noisy = compare_psnr(img[..., :n_channels],
                                          imgn[..., :n_channels])
                denoised = restoration.denoise_nl_means(imgn[..., :n_channels],
                                                        3, 5, h=0.75 * sigma,
                                                        fast_mode=fast_mode,
                                                        multichannel=True,
                                                        sigma=s)
                psnr_denoised = compare_psnr(denoised[..., :n_channels],
                                             img[..., :n_channels])
                # make sure noise is reduced
                assert_(psnr_denoised > psnr_noisy)

def get_denoise_metrics(input, output, report):
    image_file_name1 = input
    image_file_name2 = output

    image_name1 = io.imread(image_file_name1)
    image_name2 = io.imread(image_file_name2)

    # estimate the standard deiviation of the images

    std_1 = numpy.std(numpy.std(numpy.array(image_name1)))
    std_2 = numpy.std(numpy.std(numpy.array(image_name2)))

    print("std is %2.10f" % std_1)

    # print ("Standard deviation of the images are"%(std_1,std_2))

    # estimate the peak signal to noise ratio (PSNR) between the image

    peak_signal_to_noise_ratio = measure.compare_psnr(image_name1, image_name2)

    print("Peak signal to noise ratio is %s" % peak_signal_to_noise_ratio)

    # estimate the mean square error between the images

    mse = measure.compare_mse(image_name1, image_name2)

    print("Mean square error between the images is %s" % mse)

    # estimate the normalised root mean square error between the images

    rmse = measure.compare_nrmse(image_name1, image_name2)

    print("Normalised root mean square error between the images is %s" % rmse)

    resp = open(report, 'w')
    resp.write("std1 is %2.10f \n" % std_1)
    resp.write("std2 is %2.10f \n" % std_2)
    resp.write(
        "Peak signal to noise ratio is %s \n" %
        peak_signal_to_noise_ratio)
    resp.write("Mean square error between the images is %s \n" % mse)
    resp.write(
        "Normalised root mean squre error between the images is %s \n" %
        rmse)
    resp.close()

def run_metrics(image_file_name1,image_file_name2 ):
    image_name1 = io.imread(image_file_name1)
    image_name2 = io.imread(image_file_name2)
    peak_signal_to_noise_ratio = measure.compare_psnr (image_name1,image_name2)
    print ("PSNR Peak signal to noise ratio is %s"%peak_signal_to_noise_ratio)
    mse = measure.compare_mse(image_name1,image_name2)
    print  ("MSE Mean square error between the images is %s"%mse)
    rmse = measure.compare_nrmse(image_name1,image_name2)
    print  ("RMSE Normalised root mean square error between the images is %s"%rmse)
    ssim = measure.compare_ssim(image_name1,image_name2, multichannel=True)
    print ("SSIM Structural Similarity Index is %s"%ssim)
    #[M3,M4] = minkowski_distance(image_name1,image_name2)
    #print ("Minkowski distance is %s %s"%(M3,M4))
    #AD = average_difference(image_name1,image_name2)
    #print ("AD Average difference is %s"%AD)
    #SC = structural_content(image_name1,image_name2)
    #print ("SC Structural Content is %s"%SC)
    #NK = normalised_cross_correlation(image_name1,image_name2)
    #print ("NK normalised cross correlation is %s"%NK)
    #MD = maximum_difference(image_name1,image_name2)
    #print ("Maximum difference is %s"%MD)
    return {'peaktonoise':peak_signal_to_noise_ratio ,'mse': mse, 'rmse': rmse, 'ssim':ssim,'score':peak_signal_to_noise_ratio}

def test_PSNR_vs_IPOL():
    # Tests vs. imdiff result from the following IPOL article and code:
    # http://www.ipol.im/pub/art/2011/g_lmii/
    p_IPOL = 22.4497
    p = compare_psnr(cam, cam_noisy)
    assert_almost_equal(p, p_IPOL, decimal=4)

print image_name1.shape
print image_name2.shape

#estimate the standard deiviation of the images

std_1 = numpy.std (numpy.std (numpy.array(image_name1)))
std_2 = numpy.std (numpy.std (numpy.array(image_name2)))

print ("std is %2.10f"%std_1)

#print ("Standard deviation of the images are"%(std_1,std_2))

#estimate the peak signal to noise ratio (PSNR) between the image

peak_signal_to_noise_ratio = measure.compare_psnr (image_name1,image_name2)

print ("Peak signal to noise ratio is %s"%peak_signal_to_noise_ratio)

# estimate the mean square error between the images

mse = measure.compare_mse(image_name1,image_name2)

print  ("Mean square error between the images is %s"%mse)

# estimate the normalised root mean square error between the images

rmse = measure.compare_nrmse(image_name1,image_name2)
ssim = measure.compare_ssim(image_name1,image_name2)

print  ("Normalised root mean squre error between the images is %s"%rmse)