from PIL import Image
import numpy as np
from scipy.signal import convolve2d
from scipy import stats

img = imread('lettera.bmp')
arr = np.asarray(img, int)
arr = arr[:0:-1,:0:-1]  # Flip it
mn = np.mean(arr)
arr[arr < mn] = -1
arr[arr >= mn] = 1
original_mean = np.mean(arr)
print original_mean

sigma = 2
noisy_arr = arr + sigma*np.random.normal(size=arr.shape)

imshow(arr)

imshow(noisy_arr)

noisy_mean = np.mean(noisy_arr)
print noisy_mean
imshow(noisy_arr > noisy_mean)

noisy_image_copy = noisy_arr.copy()
noisy_image_iterations = []

burn_in = 10
iterations = 5
for i in range(burn_in + iterations):
    for x in range(noisy_image_copy.shape[0]):
        low_range_x = max(x - 1, 0)
        high_range_x = min(x + 1, noisy_image_copy.shape[0])
        for y in range(noisy_image_copy.shape[1]):
            low_range_y = max(y - 1, 0)
            high_range_y = min(y + 1, noisy_image_copy.shape[1])
            eta = sum(noisy_image_copy[low_range_x:high_range_x+1, low_range_y:high_range_y+1]) - noisy_image_copy[x, y]
            local_evidence = log(stats.norm.pdf(noisy_arr[x, y], 1, 2) / stats.norm.pdf(noisy_arr[x, y], -1, 2))
            threshold = 2 * eta - local_evidence
            sigmoid_threshold = exp(threshold) / (exp(threshold) + 1)
            filter_matrix = np.random.random() < sigmoid_threshold
            noisy_image_copy[x, y] = 1 if filter_matrix else -1
    if i > burn_in:
        noisy_image_iterations.append(noisy_image_copy.copy())

noisy_image_copy[low_range_x:high_range_x+1, low_range_y:high_range_y+1]

imshow(noisy_image_copy)

sum(abs(arr - noisy_image_copy)) / arr.size

sum(abs(arr - np.array(noisy_image_iterations).mean(axis=0) > 0)) / float(arr.size)

# convolution image
w_conv = np.ones((3, 3)) 
w_conv[1, 1] = 0

noisy_image_copy = noisy_arr.copy()
noisy_image_iterations = []
burn_in = 10
iterations = 5
for i in range(burn_in + iterations):
    eta_arr = convolve2d(noisy_image_copy, w_conv, mode='same')
    local_evidence = log(stats.norm.pdf(noisy_arr, 1, 2) / stats.norm.pdf(noisy_arr, -1, 2))
    threshold = 2 * eta_arr - local_evidence
    sigmoid_threshold = exp(threshold) / (exp(threshold) + 1)
    filter_matrix = np.random.random(noisy_image_copy.shape) < sigmoid_threshold
    noisy_image_copy[filter_matrix] = 1
    noisy_image_copy[~filter_matrix] = -1
    if i > burn_in:
        noisy_image_iterations.append(noisy_image_copy.copy())
    

imshow(np.array(noisy_image_iterations).mean(axis=0) > 0)

imshow(noisy_image_copy)

sum(abs(arr - np.array(noisy_image_iterations).mean(axis=0) > 0)) / float(arr.size)

sum(abs(arr - noisy_image_copy)) / arr.size