# The noise in this piece of code is taken from the accepted answer in https://stackoverflow.com/questions/22937589/how-to-add-noise-gaussian-salt-and-pepper-etc-to-image-in-python-with-opencv
from skimage.measure import compare_ssim
import numpy as np
import cv2
import math
from PIL import Image
def distance(a,b):
    return math.sqrt((a[1]-b[1])**2 + (a[0]-b[0])**2)
def ti_radgrad(d):
    RG = np.empty((d,d,3), np.uint8)
    mp = d//2 #midpoint / maximum path
    for y in range(d):
        for x in range(d):
            di = abs(distance((y,x),(mp,mp))/mp)
            dis = di**8
            g = int(max(0,min(255,255*dis)))
            RG[y,x,:] = g
    return RG
def noisy(noise_typ,image):
   if noise_typ == "gauss":
      row,col,ch= image.shape
      mean = 0
      var = 0.1
      sigma = var**0.5
      gauss = np.random.normal(mean,sigma,(row,col,ch))
      gauss = gauss.reshape(row,col,ch)
      noisy = image + gauss
      return noisy
   elif noise_typ == "s&p":
      row,col,ch = image.shape
      s_vs_p = 0.5
      amount = 0.004
      out = np.copy(image)
      # Salt mode
      num_salt = np.ceil(amount * image.size * s_vs_p)
      coords = [np.random.randint(0, i - 1, int(num_salt))
              for i in image.shape]
      out[coords] = 1

      # Pepper mode
      num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
      coords = [np.random.randint(0, i - 1, int(num_pepper))
              for i in image.shape]
      out[coords] = 0
      return out
   elif noise_typ == "poisson":
      vals = len(np.unique(image))
      vals = 2 ** np.ceil(np.log2(vals))
      noisy = np.random.poisson(image * vals) / float(vals)
      return noisy
   elif noise_typ =="speckle":
      row,col,ch = image.shape
      gauss = np.random.randn(row,col,ch)
      gauss = gauss.reshape(row,col,ch)
      noisy = image + image * gauss
      return noisy
def mark(im,c):
    im[20:40,20:40,c] = 255
    return im

image = ti_radgrad(256) # 256x256 RGB (numpy array uint8)
image = noisy("speckle", image)
image = cv2.blur(image, (5,5))
Image.fromarray(mark(image,0).astype(np.uint8)).show("original") # red square in top-left
noise = noisy("speckle", image.copy())
n_obs = 5
observations = np.zeros((5,256,256,3),dtype=np.uint8) # error: OpenCV(3.4.2) /tmp/build/80754af9/opencv-suite_1535558553474/work/modules/core/src/arithm.cpp:656: error: (-209:Sizes of input arguments do not match) The operation is neither 'array op array' (where arrays have the same size and the same number of channels), nor 'array op scalar', nor 'scalar op array' in function 'arithm_op'
for n in range(n_obs):
    ni = cv2.addWeighted(image,0.1+n*0.1,noised,0.9-n*0.1,0).astype(np.uint8)
    print(ni.dtype, ni.shape, type(ni))
    observations[n,:,:,:] = ni
#observations = np.dstack(observations)
print(observations.dtype, observations.shape, type(observations))
Image.fromarray(mark(observations[3,:,:,:],1)).show("noise") # green square
result = observations[4,:,:,:].copy().astype(np.uint8)
result_diff = result.copy()
print(result.dtype, result.shape, type(result))
#result = [noised.copy()]*n_obs # TypeError: result is not a numpy array, neither a scalar
for chan in range(3):
    cv2.denoise_TVL1(observations[2,:,:,chan], result[:,:,chan], 0.1, niters=200)
#Image.fromarray(mark(result,2).astype(np.uint8)).show("denoise") # cyan square
Image.fromarray(result.astype(np.uint8)).show("denoise")
diff = compare_ssim(result, result_diff, win_size=9, multichannel=True)
print(diff)

======== OUTPUT ============
uint8 (256, 256, 3) <class 'numpy.ndarray'>
uint8 (256, 256, 3) <class 'numpy.ndarray'>
uint8 (256, 256, 3) <class 'numpy.ndarray'>
uint8 (256, 256, 3) <class 'numpy.ndarray'>
uint8 (256, 256, 3) <class 'numpy.ndarray'>
uint8 (5, 256, 256, 3) <class 'numpy.ndarray'>
uint8 (256, 256, 3) <class 'numpy.ndarray'>
1.0