BioMulti_L_NL_Model_Complex_Ongoing.py

#!usr/bin/env python3
# -*- coding: utf-8 -*-

'''
======================================================================
DigitialBrain Version 2.0
Copyright(c) 2020  Qiang Li
All Rights Reserved.
qiang.li@uv.es
Distributed under the (new) BSD License.
######################################################################
----------------------------------------------------------------------
Permission to use, copy, or modify this software and its documentation
for educational and research purposes only and without fee is here
granted, provided that this copyright notice and the original authors'
names appear on all copies and supporting documentation. This program
shall not be used, rewritten, or adapted as the basis of a commercial
software or hardware product without first obtaining permission of the
authors. The authors make no representations about the suitability of
this software for any purpose. It is provided "as is" without express
or implied warranty.

I would like to thank all of the open suorce contributors in the Python
open community, because open source make this script works. 

The *.ipynb version also can get when you test code.
'''

# Dependent toolbox 
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# pip3
#sudo pip3 install MotionClouds
#sudo pip3 install NeuroTools
#sudo pip3 install statsmodels==0.10.0rc2 --pre
#sudo pip3 install pyrtools
#sudo pip3 install PyWavelets
#sudo pip3 install colour-science
#sudo pip3 install colorama

# @caution: you need first run (InstallDependent.sh) for install necessary
# dependent toolbox

# Bash InstallDependent.sh

# Input the mudules and exter-packages
from __future__ import division
from __future__ import absolute_import
from __future__ import print_function
import sys
sys.path.insert(0,'/home/qiang/QiangLi/Python_Utils_Functional/Retina-LGN-V1-ongoing/PyTorchSteerablePyramid')
#@caution: if you download via pip, you need to update some funciton:
#1). hand change scipy.msic to scipy.special otherwise you can not import out factorial function.
#2). add from imageio import imread in the utils script.
from steerable.SCFpyr_NumPy import SCFpyr_NumPy
import steerable.utils as utils
from google.colab.patches import cv2_imshow
from skimage.color import rgb2lab, lab2rgb, rgb2gray, xyz2lab, rgb2xyz, rgb2yuv, yuv2rgb, gray2rgb, xyz2rgb, rgb2hsv
import glob
from imageio import imread
import math
from skimage.transform import resize, rescale
import scipy
import scipy.misc
from PIL import Image
import matplotlib.pyplot as plt
import os
from mpl_toolkits.axes_grid1 import ImageGrid
import cv2
import numpy as np
from numpy.linalg import inv
import time
from matplotlib import transforms
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib import colors
from mpl_toolkits.mplot3d import Axes3D
from skimage.util import dtype
import torch
import torch.nn.functional as F
import PIL.Image as pim
import warnings
from skimage import data, img_as_float
from skimage import exposure
from scipy import linalg
import MotionClouds as mc
import pyrtools as pt
import pywt
import colour
from pywt._doc_utils import wavedec2_keys, draw_2d_wp_basis
from mpl_toolkits import axes_grid1
from tqdm import tqdm_notebook as tqdm
from tqdm import trange
from colorama import Fore

#%tensorflow_version 1x
import tensorflow as tf

warnings.filterwarnings('ignore')
#%matplotlib inline
plt.style.use('ggplot')
plt.matplotlib.rcParams['font.size'] = 6
#%load_ext autoreload
#%autoreload 2

#Sometimes, it will erros. Why? Go to check right postion use !ls
sys.path.insert(0, '/home/qiang/QiangLi/Python_Utils_Functional/Retina-LGN-V1-ongoing/SLIP/SLIP')
from SLIP import Image

#LogGabor can be install in the first time but when you second time install it,
#It will cause error. Here for how to fix it.
#!pip3 install --upgrade pip setuptools wheel
#!sudo apt-get install libpq-dev
sys.path.insert(0, '/home/qiang/QiangLi/Python_Utils_Functional/Retina-LGN-V1-ongoing/LogGabor/LogGabor')
from LogGabor import LogGabor
parameterfile = '/home/qiang/QiangLi/Python_Utils_Functional/Retina-LGN-V1-ongoing/LogGabor/default_param.py'
lg = LogGabor(parameterfile)
print('It works on 23 April2020')

#DWT wavelet filters 
sys.path.insert(0, '/home/qiang/QiangLi/Python_Utils_Functional/Retina-LGN-V1-ongoing/')

from nt_toolbox.general import *
from nt_toolbox.signal import *
from nt_toolbox.compute_wavelet_filter import *
print('It works on 27 April2020')
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

'''
############################################################################
#                          Visual Computing
#                          Fourier Brain
############################################################################
# Human visual inspired multi-layer LNL model. In this model, the main 
# component are:
#
#    Nature Image --> VonKries Adaptation --> ATD  (Color processing phase)
#    Wavelets Transform --> Contrast sensivity function (CSF) --> Divisive
#    Normalization(DN)  --> Noise(Gaussian or Poisson)
#
# Evalute of model with TID2008 database.
#
# Redundancy redunction measure with Total Correlation(RBIG or Cortex module)
#
# This model derivated two version script： Matlab, Python. In the future, I
# want to implemented all of these code on C++ or Java. If our goal is 
# simulate of primate brain, we need to implement all everything to High 
# performance Computer(HPC) with big framework architecture(C/C++/Java).
'''

############################################################################
# Function 
############################################################################
def prepare_colorarray(arr):
    """Check the shape of the array and convert it to
    floating point representation.
    """
    arr = np.asanyarray(arr)

    if arr.ndim != 3 or arr.shape[2] != 3:
        msg = "the input array must be have a shape == (.,.,3))"
        raise ValueError(msg)

    return dtype.img_as_float(arr)
#---------------------------------------------------------------
# Matrices that define conversion between different color spaces
# I will update more convert matrix at here.
#---------------------------------------------------------------
#rgb_to_xyz = np.array([[0.412453, 0.357580, 0.180423],
#              [0.212671, 0.715160, 0.072169],
#              [0.019334, 0.119193, 0.950227]])
#xyz_to_rgb = linalg.inv(rgb_to_xyz)

xyz_to_atd = np.array([[0.297, 0.72, -0.107],
            [-0.449, 0.29, -0.077],
            [0.086, -0.59, 0.501]])
atd_to_xyz = linalg.inv(xyz_to_atd)
atd_to_xyz_updated=np.array([[0.979, 1.189, 1.232],
              [-1.535, 0.764, 1.163],
              [0.445, 0.135, 2.079]])

srgb_to_iou=np.array([[0.2814, 0.6938, 0.0638],
          [-0.0971, 0.1458, -0.0250],
          [-0.0930, -0.2529, 0.4665]])
iou_to_srgb = linalg.inv(srgb_to_iou)

rgb_to_yuv=np.array([[0.299, 0.587, 0.114],
          [-0.147, -0.288, 0.436],
          [0.615, -0.515, -0.100]])

rgb_to_yiq=np.array([[0.299, 0.587, 0.114],
          [0.596, -0.274, -0.322],
          [0.211, -0.523, 0.312]])
#---------------------------------------------------------------
# Conversion between different color spaces
#---------------------------------------------------------------
def convert(matrix, arr):
    """Do the color space conversion.
    Parameters
    ----------
    matrix : array_like
        The 3x3 matrix to use.
    arr : array_like
        The input array.
    Returns
    -------
    out : ndarray, dtype=float
        The converted array.
    """
    arr = prepare_colorarray(arr)
    arr = np.swapaxes(arr, 0, 2)
    oldshape = arr.shape
    arr = np.reshape(arr, (3, -1))
    out = np.dot(matrix, arr)
    out.shape = oldshape
    out = np.swapaxes(out, 2, 0)

    return np.ascontiguousarray(out)

def xyz2atd(xyz):
  return convert(xyz_to_atd, xyz)
def atd2xyz(atd):
  return convert(atd_to_xyz, atd)
  
def atd2xyz_updated(atd):
  return convert(atd_to_xyz_updated, atd)
 
def linsrgb_to_srgb (linsrgb):
  gamma = 1.055 * linsrgb**(1./2.4) - 0.055
  scale = linsrgb * 12.92
  return np.where (linsrgb > 0.0031308, gamma, scale)

def srgb2iou(srgb):
  return convert(srgb_to_iou, srgb)

def iou2srgb(iou):
  return convert(iou_to_srgb,iou)

def rgb2yuv_matrix(rgb):
  return convert(rgb_to_yuv,rgb)

def rgb2yiq(rgb):
  return convert(rgb_to_yiq, rgb)
#---------------------------------------------------------------
# Check max and min value in each channel
#---------------------------------------------------------------
def plot_MinMax(X, labels=["R", "G", "B"]):
  print("-------------------------------")
  for i, lab in enumerate(labels):
    mi=np.min(X[:,:,i])*255
    ma=np.max(X[:,:,i])*255
    print("{} : MIN={:8.4f}, MAX={:8.4f}".format(lab,mi,ma))
#test=np.array(images[2])    
#plot_MinMax(test, labels=["R","G","B"])
#---------------------------------------------------------------
# Visualization the ATD opponent channel with presudo color
#---------------------------------------------------------------
def visual_channel_iou(image_iou, dim):
  """ 
  Opponent channel processing theory
  PresudoColor visualization
  The input array must be ATD color space
  """
  z=np.zeros(image_iou.shape)
  if dim != 0:
    print(np.mean(image_iou[:,:,0]))
    print(np.max(image_iou[:,:,0]))
    print(np.min(image_iou[:,:,0]))
    z[:,:,0]= np.min(image_iou[:,:,0])
  z[:,:,dim] = image_iou[:,:,dim]
  z = iou2srgb(z)
  z = (255*np.clip(z,0,1)).astype('uint8') 

  return(z)
#---------------------------------------------------------------
# Implemented Weber's and Fechner's law
#---------------------------------------------------------------
def waber(lumin, lambdaValue):
  """ 
   Weber's law nspired by Physhilogy experiment.
  """
  # lambdaValue normally select 0.6
  w = lumin**lambdaValue
  #w = (255*np.clip(w,0,1)).astype('uint8') 
  return(w)
#---------------------------------------------------------------
# Implemented Entropy 
#---------------------------------------------------------------
def entropy(img):
  """ 
  Calculate the entropy of image
  """
  hist, _ = np.histogram(img)
  hist = hist[hist > 0]
  return -np.log2(hist / hist.sum()).sum()

def show_entropy(band_name, img):    
    """ 
    Plot the entropy of image
    """
    bits = entropy(img)
    per_pixel = bits / img.size
    print(f"{band_name:3s} entropy = {bits:7.2f} bits, {per_pixel:7.6f} per pixel")

def image_entropy(img):
  """
  Calculate the entropy of an image-update funciton
  """
  histogram, _ = np.histogram(img)
  histogram_length = sum(histogram)

  samples_probability = [float(h) / histogram_length for h in histogram]

  return -sum([p * math.log(p, 2) for p in samples_probability if p != 0])

def print_entropy(band_name, img):
  """ 
  plot the entropy of image with each band name  
  """
  bits = image_entropy(img)
  #per_pixel = bits / img.size
  print(f"{band_name:3s} entropy = {bits:7.2f} bits")
#---------------------------------------------------------------
# Implemented mutual_information
#---------------------------------------------------------------
def mutual_information(hgram):
     """ 
     Mutual information for joint histogram
     """
     # Convert bins counts to probability values
     pxy = hgram / float(np.sum(hgram))
     px = np.sum(pxy, axis=1) # marginal for x over y
     py = np.sum(pxy, axis=0) # marginal for y over x
     px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
     # Now we can do the calculation using the pxy, px_py 2D arrays
     nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
     return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
#---------------------------------------------------------------
# Implemented resized_img function
#---------------------------------------------------------------
def resized_img(img):
  '''
  Resized function for image size with (256, 256)
  '''
  image_resized = resize(img, (256, 256, 3), anti_aliasing=True)
  return image_resized
#---------------------------------------------------------------
# Implemented visualization image histogram with 3D function
#---------------------------------------------------------------
def histogram3dplot(h, e, fig=None):
  '''
  Visualization 3D histogram of image. 
  '''  
  M, N, O = h.shape
  idxR = np.arange(M)
  idxG = np.arange(N)
  idxB = np.arange(O)

  R, G, B = np.meshgrid(idxR, idxG, idxB)
  a = np.diff(e[0])[0]
  b = a/2
  R = a * R + b

  a = np.diff(e[1])[0]
  b = a/2
  G = a * G + b

  a = np.diff(e[2])[0]
  b = a/2
  B = a * B + b

  colors = np.vstack((R.flatten(), G.flatten(), B.flatten())).T/255
  h = h / np.sum(h)
  if fig is not None:
      f = plt.figure(fig)
  else:
      f = plt.gcf()
  ax = f.add_subplot(111, projection='3d')     
  mxbins = np.array([M,N,O]).max()
  ax.scatter(R.flatten(), G.flatten(), B.flatten(), s=h.flatten()*(256/mxbins)**3/2, c=colors)
  #ax.set_xlabel('Red')
  #ax.set_ylabel('Green')
  #ax.set_zlabel('Blue')
#---------------------------------------------------------------
# Calculate correlation coffecient value for image split channel
#---------------------------------------------------------------
def show_coeff(band_name, img_ch, img_chs):
  '''
  Plot the correlation coefficient between two color channel
  '''
  coeff=np.corrcoef(img_ch.ravel(), img_chs.ravel())[0, 1]
  print(f"{band_name:3s} coeff = {coeff:.2f} coeff")
#---------------------------------------------------------------
# Visualization Lab colorspace with presudo-color
#---------------------------------------------------------------
def extract_single_dim_from_LAB_convert_to_RGB(image,dim):
  '''
  Visualization presudo color in the LAB color space.
  '''
  z = np.zeros(image.shape)
  if dim != 0:
      z[:,:,0]=30 ## I need brightness to plot the image along 1st or 2nd axis
  z[:,:,dim] = image[:,:,dim]
  z = lab2rgb(z)
  return(z)
#---------------------------------------------------------------
# Calculate MSE between raw image and reference image
#---------------------------------------------------------------
def mse(reference, query):
  '''
  Calculate the mse between origin image and query image
  '''
  (ref, que) = (reference.astype('double'), query.astype('double'))
  diff = ref - que
  square = (diff ** 2)
  mean = square.mean()
  return mean
#---------------------------------------------------------------
# Calculate PSNR between raw image and reference image
#---------------------------------------------------------------
def psnr(reference, query, normal=255):
  '''
    Calculate the PSNR of original image and query image.
  '''
  normalization = float(normal)
  msev = mse(reference, query)
  if msev != 0:
      value = 10.0 * np.log10(normalization * normalization / msev)
  else:
      value = float("inf")
  return value

def show_psnr(band_name, reference, query):
  '''
  plot psnr of reference image and query image
  '''
  psnr_ = psnr(reference, query, normal=255)
  print(f"{band_name:3s} psnr_ = {psnr_:.2f} psnr_")
#---------------------------------------------------------------
# Convert RGB image to gray
#---------------------------------------------------------------
def convert_gray(im_rgb):
  '''
  Convert image gray scale, image of shape (None,None,3)
  '''
  R=im_rgb[:,:,0]
  G=im_rgb[:,:,1]
  B=im_rgb[:,:,2]
  L = R*299/1000 + G*587/1000+B*114/1000

  return(L)

def convert_to_luminance(x):
  '''
  Convert color image into luminacne
  '''
  return np.dot(x[..., :3], [0.299, 0.587, 0.144]).astype('double')
#---------------------------------------------------------------
# Implemented CSF, you can adjust fourier space size and frequency
#---------------------------------------------------------------
def make_CSF(x, nfreq):
  '''
  Contrast Sensitivity Function implemented with Delay version.
  
  The CSF measures the sensitivity of human visual system to the various frequencies of visual stimuli, 
  Here we apply an adjusted CSF model given by:

  The mathmatic equation of CSF located in my CortexComputing notebook(Random section)

  Input:
    x - Define size of domain, float
    nfreq - Fourier frequency, float

  Output:
    CSF -  Fourier Space of CSF, 2darray
  '''
  #w=0.7
  #up_bound=7.8909
  #down_bound= 0.9809
  #a=2.6
  #b=0.0192
  #c=0.114
  #e=1.1
  params=[0.7, 7.8909, 0.9809, 2.6, 0.0192, 0.114, 1.1]

  N_x=nfreq
  N_x_=np.linspace(0, N_x, x+1)-nfreq/2
  N_x_up=N_x_[:-1]

  [xplane,yplane]=np.meshgrid(N_x_up, N_x_up)
  plane=(xplane+1j*yplane)  
  radfreq=np.abs(plane)	
  
  s=(1-params[0])/2*np.cos(4*np.angle(plane))+(1+params[0])/2
  radfreq=radfreq/s
  print(radfreq.shape)
  csf = params[3]*(params[4]+params[5]*radfreq)*np.exp(-(params[5]*radfreq)**params[6])
  f = radfreq < params[1]
  csf[f] = params[2]
  return csf
#---------------------------------------------------------------
# Implemented Gabor function with control contrast, luminacne
# velocity, orienttion and sigma of Gaussian function. 
#---------------------------------------------------------------
def Gabor(size, L0, c, f, theta, sigma, center=None):
  """ Make a Gabor functional with control parameters.
  :Input:
    -size: array, float
    -L0: luminance, float
    -c: contrast, float
    -theta: orientation, float
    -sigma: fwhm of gaussian, float
    -center: center coordinate of spatial map, float
  :output:
    -Gabor
  """
  x=np.arange(0, size, 1, float)
  y = x[:,np.newaxis]

  if center is None:
    x=y=int(size//2)
  else:
    x=center[0]
    y=center[1]
  [x,y]=np.meshgrid(range(-x,x),range(-y,y))  
  L=L0*(1.0+c*np.sin(2*np.pi*f*(x*np.cos(theta)+y*np.sin(theta))) * np.exp(-(x**2+y**2)/2*sigma**2))
  return L
#---------------------------------------------------------------
# Global contrast normalization
#---------------------------------------------------------------
def global_contrast_normalization(X, s, lmda, epsilon):  
  '''
    Global contrast normalization
  '''
  for i in range(len(X)):
    X_average = np.mean(X[i])
    print('Mean: ', X_average)
    x = X[i] - X_average
    # `su` is here the mean, instead of the sum
    contrast = np.sqrt(lmda + np.mean(x**2))
    x = s * x / max(contrast, epsilon)
    #global_contrast_normalization(images, 1, 10, 0.000000001)
#---------------------------------------------------------------
# Implemented gaussian filter
#---------------------------------------------------------------
def gaussian_filter(kernel_shape):
  '''
    Gaussian filter 
  '''
  x = np.zeros(kernel_shape, dtype='float32')

  def gauss(x, y, sigma=2.0):
      Z = 2 * np.pi * sigma ** 2
      return  1. / Z * np.exp(-(x ** 2 + y ** 2) / (2. * sigma ** 2))
  mid = np.floor(kernel_shape[-1] / 2.)
  for kernel_idx in range(0, kernel_shape[1]):
      for i in range(0, kernel_shape[2]):
          for j in range(0, kernel_shape[3]):
              x[0, kernel_idx, i, j] = gauss(i - mid, j - mid)
  return x / np.sum(x)
#---------------------------------------------------------------
# Local mean and local divisive
# local contrast normalization for increase feature
# which inspired by Divisive Normalization
#---------------------------------------------------------------
def DivisiveNormalization(image,radius=9):
  """
  image: torch.Tensor , .shape => (1,channels,height,width) 

  radius: Gaussian filter size (int), odd
  """
  if radius%2 == 0:
      radius += 1
  def get_gaussian_filter(kernel_shape):
      x = np.zeros(kernel_shape, dtype='float64')

      def gauss(x, y, sigma=2.0):
          Z = 2 * np.pi * sigma ** 2
          return  1. / Z * np.exp(-(x ** 2 + y ** 2) / (2. * sigma ** 2))

      mid = np.floor(kernel_shape[-1] / 2.)
      for kernel_idx in range(0, kernel_shape[1]):
          for i in range(0, kernel_shape[2]):
              for j in range(0, kernel_shape[3]):
                  x[0, kernel_idx, i, j] = gauss(i - mid, j - mid)

      return x / np.sum(x)

  n,c,h,w = image.shape[0],image.shape[1],image.shape[2],image.shape[3]

  gaussian_filter = torch.Tensor(get_gaussian_filter((1,c,radius,radius)))
  filtered_out = F.conv2d(image,gaussian_filter,padding=radius-1)
  mid = int(np.floor(gaussian_filter.shape[2] / 2.))
  ### Subtractive Normalization
  centered_image = image - filtered_out[:,:,mid:-mid,mid:-mid]

  ## Variance Calc
  sum_sqr_image = F.conv2d(centered_image.pow(2),gaussian_filter,padding=radius-1)
  s_deviation = sum_sqr_image[:,:,mid:-mid,mid:-mid].sqrt()
  per_img_mean = s_deviation.mean()

  ## Divisive Normalization
  divisor = np.maximum(per_img_mean.numpy(),s_deviation.numpy())
  divisor = np.maximum(divisor, 1e-4)
  new_image = centered_image / torch.Tensor(divisor)
  return new_image

#---------------------------------------------------------------
# Plot image and hist
#---------------------------------------------------------------
def plot_img_and_hist(image, axes, bins=256):
  """
  The script modifies from skimage. 
  Plot an image along with its histogram and cumulative histogram.
  """
  image = img_as_float(image)
  ax_img, ax_hist = axes
  ax_cdf = ax_hist.twinx()

  # Display image
  ax_img.imshow(image, cmap=plt.cm.gray)
  ax_img.set_axis_off()

  # Display histogram
  ax_hist.hist(image.ravel(), bins=bins, histtype='step', color='black')
  ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
  ax_hist.set_xlabel('Pixel intensity')
  ax_hist.set_xlim(0, 1)
  ax_hist.set_yticks([])

  # Display cumulative distribution
  img_cdf, bins = exposure.cumulative_distribution(image, bins)
  ax_cdf.plot(bins, img_cdf, 'r')
  ax_cdf.set_yticks([])

  return ax_img, ax_hist, ax_cdf
#---------------------------------------------------------------
# Adjust Gamma
#---------------------------------------------------------------
def adjust_gamma(image, gamma=0.8):  

  # The nonlinearlity convert between human visual and display screen 

  invGamma = 1.0 / gamma
  table = np.array([((i / 255.0) ** invGamma) * 255 for i in np.arange(0, 256)]).astype("uint8")
  return cv2.LUT(image, table)
#---------------------------------------------------------------
# Visualization wavelet transform coefficient with verious scale
#---------------------------------------------------------------
def plot_wavelet_ql(fW, Jmin=0):
  """
      plot_wavelet - plot wavelets coefficients.

      U = plot_wavelet(fW, Jmin):
  """
  def rescaleWav(A):
      v = abs(A).max()
      B = A.copy()
      if v > 0:
          B = .5 + .5 * A / v
      return B
  
  n = fW[0,0].shape[1]
  Jmax = int(np.log2(n)) - 1
  U = fW[0,0].copy()
  for j in np.arange(Jmax, Jmin - 1, -1):
      U[:2 ** j:,    2 ** j:2 **
          (j + 1):] = rescaleWav(U[:2 ** j:, 2 ** j:2 ** (j + 1):])
      U[2 ** j:2 ** (j + 1):, :2 **
        j:] = rescaleWav(U[2 ** j:2 ** (j + 1):, :2 ** j:])
      U[2 ** j:2 ** (j + 1):, 2 ** j:2 ** (j + 1):] = (
          rescaleWav(U[2 ** j:2 ** (j + 1):, 2 ** j:2 ** (j + 1):]))
  U[:2 ** j:, :2 ** j:] = nt.rescale(U[:2 ** j:, :2 ** j:])
  imageplot(U)
  for j in np.arange(Jmax, Jmin - 1, -1):
      plt.plot([0, 2 ** (j + 1)], [2 ** j, 2 ** j], 'r')
      plt.plot([2 ** j, 2 ** j], [0, 2 ** (j + 1)], 'r')
  plt.plot([0, n], [0, 0], 'r')
  plt.plot([0, n], [n, n], 'r')
  plt.plot([0, 0], [0, n], 'r')
  plt.plot([n, n], [0, n], 'r')
  return U
#---------------------------------------------------------------
# Wavelet transform coefficient with Daubechies filter
# Visualization coefficients
#---------------------------------------------------------------
def wavelet_transform(I, Jmin=1, h = compute_wavelet_filter("Daubechies",6)) :
    """
    2D-Wavelet decomposition, using Mallat's algorithm.
    By default, the convolution filters are those proposed by
    Ingrid Daubechies in her landmark 1988-1992 papers.
    """
    wI = perform_wavortho_transf(I, Jmin, + 1, h)
    return wI

def iwavelet_transform(wI, Jmin=1, h = compute_wavelet_filter("Daubechies",6)) :
    """
    Invert the Wavelet decomposition by rolling up the operations above.
    """
    I = perform_wavortho_transf(wI, Jmin, - 1, h)
    return I

def display(im):  
    """
    Displays an image using the methods of the 'matplotlib' library.
    """
    plt.figure(figsize=(8,8))                   
    plt.imshow( im, cmap="gray", vmin=0, vmax=1)
    plt.axis("off") 

#---------------------------------------------------------------
# Wavelet threshold coefficient with Daubechies filter
# Visualization coefficients
#---------------------------------------------------------------
def Wavelet_threshold(wI, threshold) :   
  """
  Re-implement a thresholding routine and create a copy of the Wavelet transform
  Remove all the small coefficients then Invert the new transform
  """
  wI_thresh = wI.copy()                  
  wI_thresh[ abs(wI) < threshold ] = 0   
  I_thresh = iwavelet_transform(wI_thresh)  
  return I_thresh

def thresh_hard(u,t):
  'Hard threshold for Ortho-wavelets'
  return u*(np.abs(u)>t)
  
def thresh_soft(u,t):
  'Soft threshold for Ortho-wavelets'
  return np.maximum(1-t/np.abs(u), 0)*u

#---------------------------------------------------------------
# Visualization sp3_filter coefficients
#--------------------------------------------------------------
def make_grid_coeff_ql(coeff, normalize=True):
  '''
  Visualization function for building a large image that contains the
  low-pass, high-pass and all intermediate levels in the steerable pyramid. 
  For the complex intermediate bands, the real part is visualized.

  Args:
      coeff (list): complex pyramid stored as list containing all levels
      normalize (bool, optional): Defaults to True. Whether to normalize each band

  Returns:
      np.ndarray: large image that contains grid of all bands and orientations
  '''

  M, N = coeff[0,1].shape
  Norients = pyr.num_scales
  #out = np.zeros((M * 2 - coeff[3, 3].shape[0], Norients * N))
  out=np.zeros((258,514))
  currentx, currenty = 0, 0

  for i in range(1, pyr.num_scales):
    for j in range(4):
        tmp = coeff[i,j].real
        m, n = tmp.shape
        if normalize:
            tmp = 255 * tmp/tmp.max()
        tmp[m-1,:] = 255
        tmp[:,n-1] = 255
        out[currentx:currentx+m,currenty:currenty+n] = tmp
        currenty += n
    currentx += coeff[i, 0].shape[0]
    currenty = 0

  a, b = coeff[3,3].shape
  out[currentx: currentx+a, currenty: currenty+b] = 255 * coeff[3,3]/coeff[3,3].max()
  out[0,:] = 255
  out[:,0] = 255
  return out.astype(np.uint8)

#---------------------------------------------------------------
# Visualization entropy 2D map
#--------------------------------------------------------------
def entropy_2d(signal):
  '''
  function returns entropy of a signal
  signal must be a 1-D numpy array
  '''
  lensig=signal.size
  symset=list(set(signal))
  numsym=len(symset)
  propab=[np.size(signal[signal==i])/(1.0*lensig) for i in symset]
  ent=np.sum([p*np.log2(1.0/p) for p in propab])
  return ent
#---------------------------------------------------------------
# Add colorbar with right size
#--------------------------------------------------------------
def add_colorbar(im, aspect=20, pad_fraction=0.5, **kwargs):
  """Add a vertical color bar to an image plot."""
  divider = axes_grid1.make_axes_locatable(im.axes)
  width = axes_grid1.axes_size.AxesY(im.axes, aspect=1./aspect)
  pad = axes_grid1.axes_size.Fraction(pad_fraction, width)
  current_ax = plt.gca()
  cax = divider.append_axes("right", size=width, pad=pad)
  plt.sca(current_ax)
  return im.axes.figure.colorbar(im, cax=cax, **kwargs)
#---------------------------------------------------------------
# View image and Fourier space in the same time and
# Plot the centerl response 
#--------------------------------------------------------------
def viewimage(img):
  '''
  In the current figure window, show the image and its
  central rows and columns (round((end+1)/2)), as well 
  as its amplitude spectrum and its central rows and columns.

  Input: 2d array (img)

  Output: 2d array (img and correspond its Fourier space)

  '''
  img2 = 20*np.log10(np.fft.fftshift(np.abs(np.fft.fft2(img))))
  plt.figure(figsize=(10, 10))
  plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
  plt.subplot(2,3,1)
  plt.imshow(img, cmap='gray'), plt.axis('off'), plt.title('image')
  plt.subplot(2,3,2)
  plt.plot(img[np.int64((img.shape[0]+1)/2) ,:], 'r.-' ) 
  plt.axis('off')
  plt.title('central row of image')
  plt.subplot(2,3,3)
  plt.plot(img[:, np.int64((img.shape[0]+1)/2)], 'r.-')
  plt.axis('off')
  plt.title('central column of image')
  plt.subplot(2,3,4)
  plt.imshow(img2, cmap='gray')
  plt.axis('off')
  plt.title('magnitude spectrum')
  plt.subplot(2,3,5)
  plt.plot(img2[:, np.int64((img2.shape[0]+1)/2)], 'b.-')
  plt.axis('off')
  plt.title('central row of magnitude spectrum')
  plt.subplot(2,3,6)
  plt.plot(img2[np.int64((img2.shape[0]+1)/2), :], 'b.-')
  plt.axis('off')
  plt.title('central column of magnitude spectrum')
  plt.tight_layout()
  plt.show()
#---------------------------------------------------------------
# Saturation is an element-wise exponential function.
#--------------------------------------------------------------
def saturation_f(x,g,xm,epsilon,sizeT):
  '''
  SATURATION_F is an element-wise exponential function (saturation). 
  It is good to (1) model the saturation in Wilson-Cowan recurrent networks, 
  and (2) as a crude (fixed) approximation to the luminance-brightness transform. 

  This saturation is normalized and modified to have these two good properties:
  (a) Some specific input, xm (e.g. the median, the average) maps into itself: xm = f(xm).

      f(x) = sign(x)*K*|x|^g  , where the constant K=xm^(1-g)

  (b) Plain exponential is modified close to the origin to avoid the singularity of
  the derivative of saturating exponentials at zero.
  This problem is solved by imposing a parabolic behavior below a certain
  threshold (epsilon) and imposing the continuity of the parabola and its
  first derivative at epsilon.

      f(x) = sign(x)*K*|x|^g             for |x| > epsilon
            sign(x)*(a*|x|+b*|x|^2)     for |x| <= epsilon

  with:
            a = (2-g)*K*epsilon^(g-1)
            b = (g-1)*K*epsilon^(g-2)

  The derivative is (of course) signal dependent:

      df/dx = g*K*|x|^(g-1)   for |x| > epsilon   [bigger with xm and decreases with signal]
              a + 2*b*|x|     for |x| <= epsilon  [bigger with xm and decreases with signal (note that b<0)]

  In the end, the slope at the origin depends on the constant xm (bigger for bigger xm). 

  The program gives the function and the derivative. For the inverse see INV_SATURATION_F.M

  For vector/matrix inputs x, the vector/matrix with anchor points, xm, has to be the same size as x.

  USE:    [f,dfdx] = saturation_f(x,gamma,xm,epsilon);

  x     = n*m matrix with the values 
  gamma = exponent (scalar)
  xm    = n*m matrix with the anchor values (in wavelet representations typically anchors will be different for different subbands)
  epsilon = threshold (scalar). It can also be a matrix the same size as x (again different epsilons per subband, e.g. epsilon = 1e-3*x_average)

  '''
  K = tf.pow(xm, tf.scalar_mul(1 - g, tf.ones([sizeT,1])))
  K = tf.where(tf.math.is_nan(K), tf.zeros_like(K), K)

  a = (2 - g) * K*(epsilon**(g - 1))
  a = tf.where(tf.math.is_nan(a), tf.zeros_like(a), a)
  b = (g-1) * K*(epsilon**(g - 2))
  b = tf.where(tf.math.is_nan(b), tf.zeros_like(b), b)

  pG = tf.math.greater(x, tf.ones([sizeT,1]) * epsilon)
  pG_zeros = tf.count_nonzero(pG)

  pp1 = tf.math.less_equal(x, tf.ones([sizeT,1]) * epsilon)
  pp2 = tf.math.greater_equal(x, tf.zeros([sizeT,1]))
  pp1_zeros = tf.count_nonzero(pp1)
  pp2_zeros = tf.count_nonzero(pp2)

  nG = tf.math.less(x, -tf.ones([sizeT,1]) * epsilon)
  np1 = tf.math.greater(x, -tf.ones([sizeT,1]) * epsilon)
  np2 = tf.math.less_equal(x, tf.zeros([sizeT,1]))
  nG_zeros = tf.count_nonzero(nG)
  np1_zeros = tf.count_nonzero(np1)
  np2_zeros = tf.count_nonzero(np2)

  f = x

  def f1(): return tf.where(pG, K*tf.pow(x, (g * tf.ones([sizeT,1]))), f)
  def f2(): return f
  f1 = tf.cond(tf.math.greater(pG_zeros, 0), f1, f2)

  def f3(): return tf.where(nG, -K*tf.pow(tf.abs(x), (g * tf.ones([sizeT,1]))), f1)
  def f4(): return f1
  f2 = tf.cond(tf.math.greater(nG_zeros, 0), f3, f4)

  def f5(): return tf.where(tf.math.logical_and(pp1,pp2), a * tf.abs(x) + b *tf.pow(x, 2 * tf.ones([sizeT,1])), f2)
  def f6(): return f2
  f3 = tf.cond(tf.math.greater(pp1_zeros + pp2_zeros, 1), f5, f6)

  def f7(): return tf.where(tf.math.logical_and(np1,np2), -(a * tf.abs(x) + b * tf.pow(x, 2 * tf.ones([sizeT,1]))), f3)
  def f8(): return f3
  f4 = tf.cond(tf.math.greater(np1_zeros + np2_zeros, 1), f7, f8)

  return f4

############################################################################
# Main section 
############################################################################
if __name__ is '__main__':

  ############################################################################
  #Setting parameters
  ############################################################################  
  wavelets_types=['SCFpyr_NumPy_SteerbalePyramid', 'LogGabor', 'sp3_filters', 'Daubechies', 'DWT']
  wavelets_type=wavelets_types[3]
  noise_types=['Gaussian', 'Possion']
  noise_type=noise_types[1]
  threshold=False
  adaptation_gain_control=False
  statis_wavlets=True
  saturation =False
  dim = (256, 256)
  g_sa = 0.5
  epsilon_sa = 0.1
  k_sa = 1
  deltat_sa = 1e-5
  g_sa = np.array(g_sa).astype('float32')
  epsilon_sa = np.array(epsilon_sa).astype('float32')
  ############################################################################
  #Load image and preprocess
  ############################################################################  
  image=cv2.imread('drive/My Drive/PyTorchSteerablePyramid/assets/lena.jpg', cv2.IMREAD_UNCHANGED)
  print('Original Dimensions : ',image.shape)
  im= cv2.resize(image, dim, interpolation = cv2.INTER_AREA)
  print('Resized Dimensions : ',im.shape)
  im=im/255
  b,g,r = cv2.split(im)      
  rgb_img = cv2.merge([r,g,b])     
  
  plt.figure(figsize=(4,5))
  plt.imshow(rgb_img)
  plt.axis('off')
  plt.title("Orignal")
  plt.tight_layout()
  plt.show()

  print('########################')
  print('Step 1')
  print('Load img test passed!!!')
  print('########################')
  #---------------------------------------------------------------
  # Gamma correction with monit and eye nonlinearlity
  #---------------------------------------------------------------
  im_gamma=exposure.adjust_gamma(im, gamma=0.6)
  b,g,r = cv2.split(im_gamma)      
  img_gamma_rgb = cv2.merge([r,g,b])     
  
  plt.figure(figsize=(4,5))
  plt.imshow(img_gamma_rgb)
  plt.axis("off")
  plt.tight_layout()
  plt.title("GammaCorrection")
  plt.show()

  print('########################')
  print('Step 2')
  print('Gammma correction test passed!!!')
  print('########################')
  #---------------------------------------------------------------
  # VonKries adaptation
  #---------------------------------------------------------------
  if adaptation_gain_control==True:
    M = np.asarray([[.40024,0.7076,-.08081], [-.2263,1.16532,0.0457],[0,0,.91822]])
    int_M = inv(M)
    XYZ2BGR = np.asarray([[3.240479,-1.53715,-0.498535], [-0.969256,1.875991,0.041556],[0.055648,-0.204043,1.057311]])
    BGR2XYZ = inv(XYZ2BGR)

    L2 = 0.6
    M2 = 0.6
    S2 = 0.6

    width = im_gamma.shape[0]
    height = im_gamma.shape[1]
    img_LMS = np.zeros(im_gamma.shape)
    img_XYZ = np.zeros(im_gamma.shape)
    img_XYZ_corrected = np.zeros(im_gamma.shape)
    img_corrected = np.zeros(im_gamma.shape)

    for x in trange(width, desc='Running'):
        for y in trange(height, desc='Running'):
            img_XYZ[x,y,:] = np.matmul(BGR2XYZ,im_gamma[x,y,:]/255)
            
    for x in trange(width, desc='Running'):
        for y in trange(height, desc='Running'):
            img_LMS[x,y,:] = np.matmul(M,img_XYZ[x,y,:])

    L_max = img_LMS[:,:,0].max()
    M_max = img_LMS[:,:,1].max()
    S_max = img_LMS[:,:,2].max()

    img_LMS[:,:,0] = img_LMS[:,:,0]/L_max
    img_LMS[:,:,1] = img_LMS[:,:,1]/M_max
    img_LMS[:,:,2] = img_LMS[:,:,2]/S_max

    img_LMS[:,:,0] = img_LMS[:,:,0]*L2
    img_LMS[:,:,1] = img_LMS[:,:,1]*M2
    img_LMS[:,:,2] = img_LMS[:,:,2]*S2

    for x in trange(width, desc='Running'):
        for y in trange(height, desc='Running'):
            img_XYZ_corrected[x,y,:] = np.matmul(int_M,img_LMS[x,y,:])
      
    for x in trange(width, desc='Running'):
        for y in trange(height, desc='Running'):
            img_corrected[x,y,:] = np.matmul(XYZ2BGR,img_XYZ_corrected[x,y,:])

    results_RGB = np.zeros(img_XYZ.shape)
    results_RGB[:,:,0] = img_corrected[:,:,2]
    results_RGB[:,:,1] = img_corrected[:,:,1]
    results_RGB[:,:,2] = img_corrected[:,:,0]
    VonKries=(255*np.clip(results_RGB,0,1)).astype('uint8') 

    plt.figure(figsize=(4,5))
    plt.imshow((VonKries),vmin=results_RGB.min(),vmax=results_RGB.max())
    plt.axis("off")
    plt.title("VonKries")
    plt.tight_layout()
    plt.show()

    print('########################')
    print('Step 3')
    print('VonKries test passed!!!')
    print('########################')
  
  #----------------------------------------------------------------
  # Chromatic adaptation with different methods
  #----------------------------------------------------------------
  else:
    print(sorted(colour.CHROMATIC_ADAPTATION_TRANSFORMS))
    # Test view condition
    XYZ_w = np.array([234, 240, 220])
    # Reference view condition
    XYZ_wr = np.array([224, 187, 70]) 
    method = 'Von Kries'
    VonKries=colour.adaptation.chromatic_adaptation_VonKries(rgb2xyz(im_gamma), XYZ_w, XYZ_wr, method)
    try:
      print(VonKries.shape)
    except ValueError:
      print('Dimentional error!') 

    fig=plt.figure(figsize=(4,5))
    fig.add_subplot(111)
    plt.imshow(VonKries)
    plt.title('VonKriesAdaptation')
    plt.tight_layout()
    plt.axis('off')
    plt.show()
    #---------------------------------------------------------------
  # Colorsapce conversation- destation ATD
  # Visualization channel with fake color
  #---------------------------------------------------------------
  plt.figure(figsize=(10,8))
  plt.subplots_adjust(wspace=0, hspace=0)

  xyz_=rgb2xyz(VonKries)
  IOU=xyz2atd(xyz_)

  plt.subplot(1, 3, 1)
  plt.imshow(IOU[:,:,0],cmap='gray')
  plt.axis('off')
  plt.title("A")

  plt.subplot(1, 3, 2)
  plt.imshow(IOU[:,:,1],cmap='gray')
  plt.axis('off')
  plt.title("T")

  plt.subplot(1, 3, 3)
  plt.imshow(IOU[:,:,2],cmap='gray')
  plt.axis('off')
  plt.title("D")
  plt.tight_layout()
  plt.show()
 
  print('########################')
  print('Step 4')
  print('ATD test passed!!!')
  print('########################')
  #---------------------------------------------------------------
  # Weber law application
  #---------------------------------------------------------------
  waber_img=waber(IOU[:,:,0], lambdaValue=0.6)
  plt.figure(figsize=(4,5))
  plt.imshow(waber_img, cmap='gray')
  plt.axis('off')
  plt.title('brightness')
  plt.tight_layout()
  plt.show()
 
  print('########################')
  print('Step 5')
  print('Weber law test passed!!!')
  print('########################')
  
  # 1) -- Complex steerable pyramid
  ############################################################################
  # Complex steerable pyramid
  ############################################################################
  #---------------------------------------------------------------
  # Implemented Steerable Wavelets with Pytorch
  # I will extend more options at here in the future.
  # 1. logGabor(including Wavelets)
  # 2. pyrtools, it ready have but I need to upgrade
  #    the visualization function. I will do it later.
  #---------------------------------------------------------------
  if wavelets_type== 'SCFpyr_NumPy_SteerbalePyramid':
    # Number of pyramid levels
    pyr_height = 5
    # Number of orientation bands
    pyr_nbands = 4
    # Tolerance for error checking
    tolerance = 1e-3

    pyr_numpy = SCFpyr_NumPy(pyr_height, pyr_nbands, scale_factor=1.7)
    coeff_numpy = pyr_numpy.build(waber_img)  
    v = reverse( np.sort(abs(coeff_numpy[1][0]).ravel()) )
    plt.figure(figsize=(4,6))
    plt.plot(np.log10(v))
    plt.axis('tight')
    plt.xlabel('length')
    plt.ylabel('coeffs')
    plt.show()
    reconstruction_numpy = pyr_numpy.reconstruct(coeff_numpy)
    reconstruction_numpy = (np.clip(reconstruction_numpy,0,1)*255).astype(np.uint8)
    
    grid = utils.make_grid_coeff(coeff_numpy, normalize=True)
    
    for level, _ in tqdm(enumerate(coeff_numpy), desc='Running'):

      print('Pyramid Level {level}'.format(level=level))
      coeff_level = coeff_numpy[level]
      
      if isinstance(coeff_level, np.ndarray):
          print('  NumPy.   min = {min:.3f}, max = {max:.3f},'
                ' mean = {mean:.3f}, std = {std:.3f}'.format(
                    min=np.min(coeff_level), max=np.max(coeff_level), 
                    mean=np.mean(coeff_level), std=np.std(coeff_level)))
      
      elif isinstance(coeff_level, list):

          for band, _ in tqdm(enumerate(coeff_level), desc='Running'):
              band_level = coeff_level[band]
              print('  Orientation Band {}'.format(band))
              print('    NumPy.   min = {min:.3f}, max = {max:.3f},'
                    ' mean = {mean:.3f}, std = {std:.3f}'.format(
                        min=np.min(band_level), max=np.max(band_level), 
                        mean=np.mean(band_level), std=np.std(band_level)))
              
    #cv2.imwrite('assets/coeff.png', grid)
    #cv2_imshow(np.rot90(grid.T))
    cv2_imshow(np.ascontiguousarray(grid))
    cv2_imshow(reconstruction_numpy.astype(np.uint8))
    print('########################################################')
    print('Step 6')
    print('SteerablePyramid build and reconstructure test passed!!!')
    print('########################################################')
    
    ############################################################################
    # Flatten all decomposition frequency sub-band 
    ############################################################################
    imgList=[]
    
    band_=coeff_numpy[0]
    imgList.append(band_)
    
    for s in trange(3, desc='Running'):
      for r in trange(4, desc='Running'):
        band_s = coeff_numpy[s+1][r]
        imgList.append(band_s)

    band_po=coeff_numpy[4]
    imgList.append(band_po)
    
    plt.figure(figsize=(8,6))
    plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
    plt.margins(0,0)
    for i in trange(len(imgList[1:13]), desc='Running'):
      plt.subplot(3,4,i+1)
      plt.imshow(imgList[i].real, cmap='gray')
      plt.axis('off')
    plt.tight_layout()
    plt.show()

    print('############################')
    print('Step 6')
    print('decomposition test passed!!!')
    print('############################')
    ############################################################################
    # CSF -  operate in the Fourier space then convert back to spatial domain
    #     -  all wavelets domain 
    ############################################################################
    CSF_wavelets_vis=[]
    CSF_wavelets=[]
    CSF_wavelet_l=[]
    CSF_wavelet_m=[]
    CSF_wavelet_h=[]

    res_LL=make_CSF(x=256, nfreq=32)
    Img_CSF_LL = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgList[0].real))* res_LL))))
    CSF_wavelets_vis.append(Img_CSF_LL)
    
    for r in trange(len(imgList[1:5]), desc='Running'):
      res_L=make_CSF(x=256, nfreq=32)
      Img_CSF_L = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgList[r].real))* res_L))))
      CSF_wavelet_l.append(Img_CSF_L)
      CSF_wavelets_vis.append(Img_CSF_L)
    for t in tqdm(range(len(imgList))[5:9], desc='Running'):  
      res_M=make_CSF(x=128, nfreq=32)
      Img_CSF_M = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgList[t].real))* res_M))))
      CSF_wavelet_m.append(Img_CSF_M)
      CSF_wavelets_vis.append(Img_CSF_M)
    for y in tqdm(range(len(imgList))[9:13], desc='Running'):  
      res_H=make_CSF(x=64, nfreq=32)
      Img_CSF_H = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgList[y].real))* res_H))))
      CSF_wavelet_h.append(Img_CSF_H)
      CSF_wavelets_vis.append(Img_CSF_H)
    
    res_PP=make_CSF(x=32, nfreq=32)
    Img_CSF_PP = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgList[13].real))* res_PP))))
    CSF_wavelets_vis.append(Img_CSF_PP)
    
    CSF_wavelets=[Img_CSF_LL, CSF_wavelet_l, CSF_wavelet_m, CSF_wavelet_h, Img_CSF_PP]

    if len(CSF_wavelets)==5:
      print('CSF filter was right!!')
    else:
      print('CSF filter test wrong, go back check!!')
    
    plt.figure(figsize=(8,6))
    plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
    plt.margins(0,0)
    for j in trange(len(CSF_wavelets_vis[1:13]), desc='Running'):
      plt.subplot(3,4,j+1)
      plt.imshow(CSF_wavelets_vis[j], cmap='gray')
      plt.axis('off')
    plt.tight_layout()
    plt.show()
    print('##################################')
    print('Step 7')
    print('CSF filter wavelets test passed!!!')
    print('##################################')

    ############################################################################
    # Upsampling 
    ############################################################################
    reconstruction_num = pyr_numpy.reconstruct(CSF_wavelets)
    reconstruction_num = (np.clip(reconstruction_num, 0,1)*255).astype(np.uint8)
    cv2_imshow(reconstruction_num)
    print('#######################################')
    print('Step 7')
    print('after CSF Reconstructure test passed!!!')
    print('#######################################')
    
    ############################################################################
    # Saturation of each sub-band of steerbale wavelets after CSF filter
    ############################################################################
    if saturation == True:  
      sizeT_subL=CSF_wavelets_vis[0].shape[0]
      sizeT_subM=CSF_wavelets_vis[5].shape[0]
      sizeT_subS=CSF_wavelets_vis[9].shape[0]
      sizeT_subR=CSF_wavelets_vis[13].shape[0]
      xm_l = np.array(np.ones([256,256])).astype('float32') 
      xm_m = np.array(np.ones([128,128])).astype('float32') 
      xm_s = np.array(np.ones([64,64])).astype('float32') 
      xm_r = np.array(np.ones([32,32])).astype('float32') 

      CsfWaveletSa=[]
      CSF_wavelets_sa=[]
      CSF_wavelet_l_sa=[]
      CSF_wavelet_m_sa=[]
      CSF_wavelet_s_sa=[]

      with tf.Session() as sess:
        for sub_l in tqdm(range(len(CSF_wavelets_vis))[:5], desc='Running'):
          f_l = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis[sub_l]).astype('float32'), g_sa, xm_l, epsilon_sa, sizeT_subL)), dtype=tf.float32))
          CsfWaveletSa.append(f_l)
          CSF_wavelet_l_sa.append(f_l)
      
      with tf.Session() as sess:
        for sub_m in tqdm(range(len(CSF_wavelets_vis))[5:9], desc='Running'):
          f_m = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis[sub_m]).astype('float32'), g_sa, xm_m, epsilon_sa, sizeT_subM)), dtype=tf.float32))
          CsfWaveletSa.append(f_m)
          CSF_wavelet_m_sa.append(f_m)
      
      with tf.Session() as sess:
        for sub_s in tqdm(range(len(CSF_wavelets_vis))[9:13], desc='Running'):
          f_s = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis[sub_s]).astype('float32'), g_sa, xm_s, epsilon_sa, sizeT_subS)), dtype=tf.float32))
          CsfWaveletSa.append(f_s)
          CSF_wavelet_s_sa.append(f_s)
      
      with tf.Session() as sess:
        f_r = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis[13]).astype('float32'), g_sa, xm_r, epsilon_sa, sizeT_subR)), dtype=tf.float32))
        CsfWaveletSa.append(f_r)
      
      CSF_wavelets_sa=[CSF_wavelet_l_sa[0], CSF_wavelet_l_sa[1:], CSF_wavelet_m_sa, CSF_wavelet_s_sa, f_r]

      plt.figure(figsize=(8,6))
      plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
      plt.margins(0,0)
      for sat in trange(len(CsfWaveletSa[1:13]), desc='Running'):
        plt.subplot(3,4,sat+1)
        plt.imshow(CsfWaveletSa[sat], cmap='gray')
        plt.axis('off')
      plt.tight_layout()
      plt.show()
      print('##################################')
      print('Step 8')
      print('Saturation test passed!!!')
      print('##################################')
      ############################################################################
      # Upsampling for saturation 
      ############################################################################
      reconstruction_num_sa = pyr_numpy.reconstruct(CSF_wavelets_sa)
      reconstruction_num_sa = (np.clip(reconstruction_num_sa, 0,1)*255).astype(np.uint8)
      cv2_imshow(reconstruction_num_sa)
      print('#######################################')
      print('Step 8')
      print('Saturation Reconstructure test passed!!!')
      print('#######################################')
      ############################################################################
      # Neural noise - 1) Gaussian  2) Possion 
      ############################################################################
      #Gaussian
      if noise_type=='Gaussian':
        sigma=.5
        f_gaussian = reconstruction_num_sa + sigma*np.random.standard_normal(reconstruction_num_sa.shape)
        cv2_imshow((np.clip(f_gaussian, 0,1)*255).astype('uint8'))

      #Possion
      if noise_type=='Possion':
        f_possion = np.random.poisson(reconstruction_num_sa, reconstruction_num_sa.shape)
        cv2_imshow((np.clip(f_possion, 0,1)*255).astype('uint8'))
      
      print('####################################')
      print('Step 9')
      print('Simulate neural noise test passed!!!')
      print('####################################')

      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('Retina-LGN-V1 Visual Computing Finished')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
    ############################################################################
    # Nonlinear Divisive Normalizaton model(DN)
    # DN model inspired by Human Visual Cortex
    ############################################################################
    else:
      DN_band=[]
      reconstruction_num_level_LL = np.reshape(CSF_wavelets_vis[0], (256, 256, 1))
      DN_LL = DivisiveNormalization(torch.Tensor([reconstruction_num_level_LL.transpose((2,0,1))]),radius=9) 
      ret_LL = DN_LL[0].numpy().transpose((1,2,0)) 
      scaled_ret_LL = (ret_LL - ret_LL.min())/(ret_LL.max() - ret_LL.min()) 
      scaled_ret_LL = np.reshape(scaled_ret_LL, (256, 256))  
      DN_band.append(scaled_ret_LL)
      for j in trange(1, 5, desc='Running'):
        reconstruction_num_level_l = np.reshape(CSF_wavelets_vis[j], (256, 256, 1))
        DN_l = DivisiveNormalization(torch.Tensor([reconstruction_num_level_l.transpose((2,0,1))]),radius=9) 
        ret_l = DN_l[0].numpy().transpose((1,2,0)) 
        scaled_ret_l = (ret_l - ret_l.min())/(ret_l.max() - ret_l.min()) 
        scaled_ret_l = np.reshape(scaled_ret_l, (256, 256))
        DN_band.append(scaled_ret_l)
      for k in trange(5, 9, desc='Running'):
        reconstruction_num_level_m = np.reshape(CSF_wavelets_vis[k], (128, 128, 1))
        DN_m = DivisiveNormalization(torch.Tensor([reconstruction_num_level_m.transpose((2,0,1))]),radius=9) 
        ret_m = DN_m[0].numpy().transpose((1,2,0)) 
        scaled_ret_m = (ret_m - ret_m.min())/(ret_m.max() - ret_m.min()) 
        scaled_ret_m = np.reshape(scaled_ret_m, (128, 128))
        DN_band.append(scaled_ret_m)
      for p in trange(9, 13, desc='Running'):
        reconstruction_num_level_h = np.reshape(CSF_wavelets_vis[p], (64, 64, 1))
        DN_h = DivisiveNormalization(torch.Tensor([reconstruction_num_level_h.transpose((2,0,1))]),radius=9) 
        ret_h = DN_h[0].numpy().transpose((1,2,0)) 
        scaled_ret_h = (ret_h - ret_h.min())/(ret_h.max() - ret_h.min()) 
        scaled_ret_h = np.reshape(scaled_ret_h, (64, 64))
        DN_band.append(scaled_ret_h)
      reconstruction_num_level_HH = np.reshape(CSF_wavelets_vis[13], (32, 32, 1))
      DN_HH = DivisiveNormalization(torch.Tensor([reconstruction_num_level_HH.transpose((2,0,1))]),radius=9) 
      ret_HH = DN_HH[0].numpy().transpose((1,2,0)) 
      scaled_ret_HH = (ret_HH - ret_HH.min())/(ret_HH.max() - ret_HH.min()) 
      scaled_ret_HH = np.reshape(scaled_ret_HH, (32, 32))  
      DN_band.append(scaled_ret_HH)
      
      try:
        print(len(DN_band))
        if len(DN_band)==14:
          print('DN model works!')

          plt.figure(figsize=(8,6))
          plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
          plt.margins(0,0)
          for j in trange(len(DN_band[1:13]), desc='Running'):
            plt.subplot(3,4,j+1)
            plt.imshow(DN_band[j], cmap='gray')
            plt.axis('off')
          plt.tight_layout()
          plt.show()
          print('#############################')
          print('Step 8')
          print('DN model test passed!!!')
          print('#############################')

      except ValueError:
        print('DN model errors')
      ############################################################################
      # Upsampling 
      ############################################################################
      reconstruction_num = np.reshape(reconstruction_num, (256, 256, 1))
      DN = DivisiveNormalization(torch.Tensor([reconstruction_num.transpose((2,0,1))]),radius=9) 
      ret = DN[0].numpy().transpose((1,2,0)) 
      scaled_ret = (ret - ret.min())/(ret.max() - ret.min()) 
      scaled_ret = np.reshape(scaled_ret, (256, 256))
      #plt.imsave('{}_lcn.png'.format(i),scaled_ret)
      plt.figure(figsize=(4,6))
      plt.imshow(scaled_ret, cmap='gray')
      plt.axis('off')
      plt.show()
      print('################################')
      print('Step 8')
      print('DN Reconstructure test passed!!!')
      print('################################')
      ############################################################################
      # Neural noise - 1) Gaussian  2) Possion 
      ############################################################################
      #Gaussian
      if noise_type=='Gaussian':
        sigma=.5
        f_gaussian = scaled_ret + sigma*np.random.standard_normal(scaled_ret.shape)
        cv2_imshow((np.clip(f_gaussian, 0,1)*255).astype('uint8'))

      #Possion
      if noise_type=='Possion':
        f_possion = np.random.poisson(scaled_ret, scaled_ret.shape)
        cv2_imshow((np.clip(f_possion, 0,1)*255).astype('uint8'))
      
      print('####################################')
      print('Step 9')
      print('Simulate neural noise test passed!!!')
      print('####################################')

      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('Retina-LGN-V1 Visual Computing Finished')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')

  # 2) -- LogGabor 
  ############################################################################
  # Build Pyramid with LogGabor
  # Reference: Self-Invertible 2D LogGabor Wavelets, 2007  
  ############################################################################
  if wavelets_type== 'LogGabor':
    phi = (np.sqrt(5)+1) / 2
    fig = plt.figure(figsize=(8, 6))
    border = 0.2
    a = fig.add_axes((border, border*.5, .9-border, 1.-border*.5), facecolor='w')
    a.axis(c='b', lw=0)
    N_X_gabor = 128
    lg.set_size((N_X_gabor, N_X_gabor))
    samples = 8
    image_gabor = np.zeros((N_X_gabor*samples, N_X_gabor*samples))
    for i in trange(samples, desc='Running'):
        for j in trange(samples, desc='Running'):
            params = {'sf_0':1./16, 'B_sf':lg.pe.B_sf, 'theta':i*np.pi/samples, 'B_theta':lg.pe.B_theta}
            FT_lg = lg.loggabor(N_X_gabor/2, N_X_gabor/2, **params)
            image_temp = lg.invert(FT_lg * np.exp(-1j*2*np.pi*j/samples))
            image_gabor[(i*N_X_gabor):(i*N_X_gabor+N_X_gabor), (j*N_X_gabor):(j*N_X_gabor+N_X_gabor)] = image_temp
    image_gabor /= np.abs(image_gabor).max()
    lg.set_size(image_gabor)
    fig, a = lg.imshow(image_gabor, fig=fig, ax=a)
    a.grid(False)
    plt.show()

    fig, ax = plt.subplots(figsize=(6, 6))
    env_total = np.zeros_like(image_gabor)
    for sf_0_ in lg.sf_0:
        for theta_ in lg.theta:
            for sym in [0, np.pi]:
                FT_lg = lg.loggabor(0, 0, sf_0=sf_0_, theta=theta_+sym, B_sf=lg.pe.B_sf,  B_theta=lg.pe.B_theta)
                env = np.absolute(FT_lg)**2
                env_total += env
                ax.contour(env, levels=[env.max()/2], lw=1)
    fig.suptitle('Tiling of Fourier space using the pyramid')
    ax.set_xlabel(r'$f_Y$')
    ax.set_ylabel(r'$f_X$')
    plt.tight_layout()
    plt.axis('off')
    plt.show()

    lg.pe.base_levels = phi
    lg.set_size(waber_img)
    lg.init()
    C = lg.linear_pyramid(waber_img)
    fig, axs = lg.golden_pyramid(C, mask=True, spiral=True, fig_width=13)
    plt.show()
    print('##########################################')
    print('Step 6')
    print('Build pyramid with LogGabor test passed!!!')
    print('##########################################')

  # 3) -- sp3——filters 
  ############################################################################
  # Build Pyramid with sp3_filters from pyrtools
  ############################################################################
  if wavelets_type== 'sp3_filters':
    filt = 'sp3_filters' 
    pyr = pt.pyramids.SteerablePyramidSpace(waber_img, height=4, order=3)
    imgCoeffs = []
    for j in trange(pyr.num_scales, desc='Running'):
        for s in trange(pyr.num_scales, desc='Running'):
            band = pyr.pyr_coeffs[j,s]
            imgCoeffs.append(band)

    plt.figure(figsize=(8,6))
    for m in trange(len(imgCoeffs), desc='Running'):
      plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
      plt.subplot(len(imgCoeffs)/4, len(imgCoeffs)/4, m+1 )
      plt.imshow(imgCoeffs[m], cmap='gray')
      plt.axis('off')
      plt.tight_layout()

    coeff_=pyr.pyr_coeffs
    plt.figure(figsize = (6,6))
    plot_wavelet_ql(coeff_, 1)
    plt.show()
    print('#############################################')
    print('Step 6')
    print('Build pyramid with sp3_filters test passed!!!')
    print('#############################################')
    ############################################################################
    # CSF -  operate in the Fourier space then convert back to spatial domain
    #     -  all wavelets domain 
    ############################################################################
    CSF_wavelets_vis_sp=[]
    CSF_wavelets_sp=[]
    CSF_wavelet_l_sp=[]
    CSF_wavelet_m_sp=[]
    CSF_wavelet_h_sp=[]
    CSF_wavelet_s_sp=[]
    for r_sp in trange(len(imgCoeffs[0:4]), desc='Running'):
      res_L_sp=make_CSF(x=256, nfreq=32)
      Img_CSF_L_sp = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgCoeffs[r_sp].real))* res_L_sp))))
      CSF_wavelet_l_sp.append(Img_CSF_L_sp)
      CSF_wavelets_vis_sp.append(Img_CSF_L_sp)
    for t_sp in tqdm(range(len(imgCoeffs))[4:8], desc='Running'):  
      res_M_sp=make_CSF(x=128, nfreq=32)
      Img_CSF_M_sp = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgCoeffs[t_sp].real))* res_M_sp))))
      CSF_wavelet_m_sp.append(Img_CSF_M_sp)
      CSF_wavelets_vis_sp.append(Img_CSF_M_sp)
    for y_sp in tqdm(range(len(imgCoeffs))[8:12], desc='Running'):  
      res_H_sp=make_CSF(x=64, nfreq=32)
      Img_CSF_H_sp = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgCoeffs[y_sp].real))* res_H_sp))))
      CSF_wavelet_h_sp.append(Img_CSF_H_sp)
      CSF_wavelets_vis_sp.append(Img_CSF_H_sp)
    for e_sp in tqdm(range(len(imgCoeffs))[12:16], desc='Running'):  
      res_S_sp=make_CSF(x=32, nfreq=32)
      Img_CSF_S_sp = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(imgCoeffs[e_sp].real))* res_S_sp))))
      CSF_wavelet_s_sp.append(Img_CSF_S_sp)
      CSF_wavelets_vis_sp.append(Img_CSF_S_sp)
    
    CSF_wavelets_sp=[CSF_wavelet_l_sp[0], CSF_wavelet_l_sp, CSF_wavelet_m_sp, CSF_wavelet_h_sp, CSF_wavelet_s_sp[0]]

    if len(CSF_wavelets_sp)==5:
      print('CSF filter was right!!')
    else:
      print('CSF filter test wrong, go back check!!')
    
    plt.figure(figsize=(8,6))
    plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
    plt.margins(0,0)
    for j in trange(len(CSF_wavelets_vis_sp), desc='Running'):
      plt.subplot(4,4,j+1)
      plt.imshow(CSF_wavelets_vis_sp[j], cmap='gray')
      plt.axis('off')
    plt.tight_layout()
    plt.show()
    print('##################################')
    print('Step 7')
    print('CSF filter wavelets test passed!!!')
    print('##################################')
    ############################################################################
    # Upsampling 
    ############################################################################
    pyr_height=5
    pyr_nbands=4
    pyr_numpy = SCFpyr_NumPy(pyr_height, pyr_nbands, scale_factor=1.7)
    reconstruction_num_sp = pyr_numpy.reconstruct(CSF_wavelets_sp)
    reconstruction_num_sp=np.reshape(reconstruction_num_sp, (256, 256))
    plt.figure(figsize=(4,5))
    plt.imshow(reconstruction_num_sp, cmap='gray')
    plt.axis('off')
    plt.tight_layout()
    plt.show()
    print('#######################################')
    print('Step 7')
    print('after CSF Reconstructure test passed!!!')
    print('#######################################')
    
    ############################################################################
    # Saturation of each sub-band of steerbale wavelets after CSF filter
    ############################################################################
    if saturation == True:  
      sizeT_subL_sp=CSF_wavelets_vis_sp[0].shape[0]
      sizeT_subM_sp=CSF_wavelets_vis_sp[4].shape[0]
      sizeT_subS_sp=CSF_wavelets_vis_sp[8].shape[0]
      sizeT_subR_sp=CSF_wavelets_vis_sp[12].shape[0]
      
      xm_l_sp = np.array(np.ones([256,256])).astype('float32') 
      xm_m_sp = np.array(np.ones([128,128])).astype('float32') 
      xm_s_sp = np.array(np.ones([64,64])).astype('float32') 
      xm_r_sp = np.array(np.ones([32,32])).astype('float32') 

      CsfWaveletSa_sp=[]
      CSF_wavelet_l_sp=[]
      CSF_wavelet_m_sp=[]
      CSF_wavelet_s_sp=[]
      CSF_wavelet_r_sp=[]
      
      with tf.Session() as sess:
        for sub_l in tqdm(range(len(CSF_wavelets_vis_sp))[0:4], desc='Running'):
          f_l = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis_sp[sub_l]).astype('float32'), g_sa, xm_l_sp, epsilon_sa, sizeT_subL_sp)), dtype=tf.float32))
          CsfWaveletSa_sp.append(f_l)
          CSF_wavelet_l_sp.append(f_l)
      
      with tf.Session() as sess:
        for sub_m in tqdm(range(len(CSF_wavelets_vis_sp))[4:8], desc='Running'):
          f_m = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis_sp[sub_m]).astype('float32'), g_sa, xm_m_sp, epsilon_sa, sizeT_subM_sp)), dtype=tf.float32))
          CsfWaveletSa_sp.append(f_m)
          CSF_wavelet_m_sp.append(f_m)
      
      with tf.Session() as sess:
        for sub_s in tqdm(range(len(CSF_wavelets_vis_sp))[8:12], desc='Running'):
          f_s = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis_sp[sub_s]).astype('float32'), g_sa, xm_s_sp, epsilon_sa, sizeT_subS_sp)), dtype=tf.float32))
          CsfWaveletSa_sp.append(f_s)
          CSF_wavelet_s_sp.append(f_s)
      
      with tf.Session() as sess:
        for sub_r in tqdm(range(len(CSF_wavelets_vis_sp))[12:16], desc='Running'):
          f_r = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(CSF_wavelets_vis_sp[sub_r]).astype('float32'), g_sa, xm_r_sp, epsilon_sa, sizeT_subR_sp)), dtype=tf.float32))
          CsfWaveletSa_sp.append(f_r)
          CSF_wavelet_r_sp.append(f_r)
        
      CSF_wavelets_sp=[CSF_wavelet_l_sp[0], CSF_wavelet_l_sp, CSF_wavelet_m_sp, CSF_wavelet_s_sp, CSF_wavelet_r_sp[0]]

      plt.figure(figsize=(8,6))
      plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
      plt.margins(0,0)
      for sat in trange(len(CsfWaveletSa_sp), desc='Running'):
        plt.subplot(4,4,sat+1)
        plt.imshow(CsfWaveletSa_sp[sat], cmap='gray')
        plt.axis('off')
      plt.tight_layout()
      plt.show()
      print('##################################')
      print('Step 8')
      print('Saturation test passed!!!')
      print('##################################')
      ############################################################################
      # Upsampling for saturation 
      ############################################################################
      reconstruction_num_sp = pyr_numpy.reconstruct(CSF_wavelets_sp)
      reconstruction_num_sp = (np.clip(reconstruction_num_sp, 0,1)*255).astype(np.uint8)
      cv2_imshow(reconstruction_num_sp)
      print('#######################################')
      print('Step 8')
      print('Saturation Reconstructure test passed!!!')
      print('#######################################')  
      ############################################################################
      # Neural noise - 1) Gaussian  2) Possion 
      ############################################################################
      #Gaussian
      if noise_type=='Gaussian':
        sigma=.1
        f_gaussian_sp = reconstruction_num_sp + sigma*np.random.standard_normal(reconstruction_num_sp.shape)
        cv2_imshow((np.clip(f_gaussian_sp, 0,1)*255).astype('uint8'))

      #Possion
      if noise_type=='Possion':
        f_possion_sp = np.random.poisson(reconstruction_num_sp, scaled_ret_sp.shape)
        cv2_imshow((np.clip(f_possion_sp, 0,1)*255).astype('uint8'))
      
      print('####################################')
      print('Step 9')
      print('Simulate neural noise test passed!!!')
      print('####################################')

      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('Retina-LGN-V1 Visual Computing Finished')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
        
    ############################################################################
    # Nonlinear Divisive Normalizaton model(DN)
    # DN model inspired by Human Visual Cortex
    ############################################################################
    else:
      DN_band_sp=[]
      for j in tqdm(range(0, 4), desc='Running'):
        reconstruction_num_level_l_sp = np.reshape(CSF_wavelets_vis_sp[j], (256, 256, 1))
        DN_l_sp = DivisiveNormalization(torch.Tensor([reconstruction_num_level_l_sp.transpose((2,0,1))]),radius=9) 
        ret_l_sp = DN_l_sp[0].numpy().transpose((1,2,0)) 
        scaled_ret_l_sp = (ret_l_sp - ret_l_sp.min())/(ret_l_sp.max() - ret_l_sp.min()) 
        scaled_ret_l_sp = np.reshape(scaled_ret_l_sp, (256, 256))
        DN_band_sp.append(scaled_ret_l_sp)
      for k in tqdm(range(4,8), desc='Running'):
        reconstruction_num_level_m_sp = np.reshape(CSF_wavelets_vis_sp[k], (128, 128, 1))
        DN_m_sp = DivisiveNormalization(torch.Tensor([reconstruction_num_level_m_sp.transpose((2,0,1))]),radius=9) 
        ret_m_sp = DN_m_sp[0].numpy().transpose((1,2,0)) 
        scaled_ret_m_sp = (ret_m_sp - ret_m_sp.min())/(ret_m_sp.max() - ret_m_sp.min()) 
        scaled_ret_m_sp = np.reshape(scaled_ret_m_sp, (128, 128))
        DN_band_sp.append(scaled_ret_m_sp)
      for p in tqdm(range(8,12), desc='Running'):
        reconstruction_num_level_h_sp = np.reshape(CSF_wavelets_vis_sp[p], (64, 64, 1))
        DN_h_sp = DivisiveNormalization(torch.Tensor([reconstruction_num_level_h_sp.transpose((2,0,1))]),radius=9) 
        ret_h_sp = DN_h_sp[0].numpy().transpose((1,2,0)) 
        scaled_ret_h_sp = (ret_h_sp - ret_h_sp.min())/(ret_h_sp.max() - ret_h_sp.min()) 
        scaled_ret_h_sp = np.reshape(scaled_ret_h_sp, (64, 64))
        DN_band_sp.append(scaled_ret_h_sp)
      for q in tqdm(range(12,16), desc='Running'):
        reconstruction_num_level_s_sp = np.reshape(CSF_wavelets_vis_sp[q], (32, 32, 1))
        DN_s_sp = DivisiveNormalization(torch.Tensor([reconstruction_num_level_s_sp.transpose((2,0,1))]),radius=9) 
        ret_s_sp = DN_s_sp[0].numpy().transpose((1,2,0)) 
        scaled_ret_s_sp = (ret_s_sp - ret_s_sp.min())/(ret_s_sp.max() - ret_s_sp.min()) 
        scaled_ret_s_sp = np.reshape(scaled_ret_s_sp, (32, 32))
        DN_band_sp.append(scaled_ret_s_sp)

      try:
        print(len(DN_band_sp))
        if len(DN_band_sp)==16:
          print('DN model works!')

          plt.figure(figsize=(8,6))
          plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
          plt.margins(0,0)
          for j in trange(len(DN_band_sp), desc='Running'):
            plt.subplot(4,4,j+1)
            plt.imshow(DN_band_sp[j], cmap='gray')
            plt.axis('off')
          plt.tight_layout()
          plt.show()
          print('##################################')
          print('Step 8')
          print('DN model test passed!!!')
          print('##################################')

      except ValueError:
        print('DN model errors')
      ############################################################################
      # Upsampling 
      ############################################################################
      reconstruction_num_sp = np.reshape(reconstruction_num_sp, (256, 256, 1))
      DN_sp = DivisiveNormalization(torch.Tensor([reconstruction_num_sp.transpose((2,0,1))]),radius=9) 
      ret_sp = DN_sp[0].numpy().transpose((1,2,0)) 
      scaled_ret_sp = (ret_sp - ret_sp.min())/(ret_sp.max() - ret_sp.min()) 
      scaled_ret_sp = np.reshape(scaled_ret_sp, (256, 256))
      plt.figure(figsize=(4,6))
      plt.imshow(scaled_ret_sp, cmap='gray')
      plt.axis('off')
      plt.show()
      print('################################')
      print('Step 8')
      print('DN Reconstructure test passed!!!')
      print('################################')
      ############################################################################
      # Neural noise - 1) Gaussian  2) Possion 
      ############################################################################
      #Gaussian
      if noise_type=='Gaussian':
        sigma=.1
        f_gaussian_sp = scaled_ret_sp + sigma*np.random.standard_normal(scaled_ret_sp.shape)
        cv2_imshow((np.clip(f_gaussian_sp, 0,1)*255).astype('uint8'))

      #Possion
      if noise_type=='Possion':
        f_possion_sp = np.random.poisson(scaled_ret_sp, scaled_ret_sp.shape)
        cv2_imshow((np.clip(f_possion_sp, 0,1)*255).astype('uint8'))
      
      print('####################################')
      print('Step 9')
      print('Simulate neural noise test passed!!!')
      print('####################################')

      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('Retina-LGN-V1 Visual Computing Finished')
      print('+++++++++++++++++++++++++++++++++++++++')
      print('+++++++++++++++++++++++++++++++++++++++')
    
  # 4) -- Daubechies convolution filter
  ############################################################################
  # Build Pyramid with Daubechies convolution filter
  ############################################################################
  if wavelets_type== 'Daubechies':
    if threshold==True:
      wI = wavelet_transform(waber_img)     
      thesh_=Wavelet_threshold(wI, .005)  
      plt.figure(figsize = (8,8))        
      plot_wavelet(wI, 1)
      plt.axis('off')              
      plt.title('threshold Wavelet coefficients')
      plt.show()
    else:
      wI = wavelet_transform(waber_img)     
      plt.figure(figsize = (8,8))        
      plot_wavelet(wI, 1)              
      plt.show()
      ############################################################################
      # Staistics
      ############################################################################
      V_wI=wI[1:128, 128+1:256]
      viewimage(V_wI)
      plt.figure(figsize = (6,4))
      plt.subplots_adjust(wspace=0.5)        
      plt.subplot(121)
      plt.imshow(V_wI, cmap='gray')
      plt.axis('off')              
      plt.subplot(122)
      plt.hist(V_wI.ravel(), 128, histtype='step', log=True)
      plt.axis('tight')
      plt.show()
      #Energy-to-Shannon entropy ratio
      S=0
      Etot=0
      for d in wI:
          E=d**2
          P=E/np.sum(E)
          S+=-np.sum(P*np.log(P))
          Etot+=np.sum(E)
      ratio=Etot/S
      print('Energy-Shannon Ratio:{:.3f}'.format(ratio))
      ############################################################################
      #  Hard/Soft thresholding
      ############################################################################
      alpha = np.linspace(-6,6,1000)
      plt.figure(figsize=(6,4))
      plt.subplot(1,2,1)
      plt.plot(alpha, thresh_hard(alpha,1),'b-')   
      plt.axis('equal')
      plt.title('Hard threshold')
      plt.subplot(1,2,2)
      plt.plot(alpha, thresh_soft(alpha,1),'r-')
      plt.axis('equal')
      plt.title('Soft threshold')

      #reconstructure
      I_thresh = iwavelet_transform(wI)
      cv2_imshow((np.clip(I_thresh,0,1)*255).astype('uint8'))
      print('############################################')
      print('Step 6')
      print('Build pyramid with Daubechies test passed!!!')
      print('############################################')
      ############################################################################
      # CSF 
      ############################################################################
      csf_ortho=make_CSF(x=256, nfreq=32)
      csf_ortho_ = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(wI))*csf_ortho))))
      plt.figure(figsize = (8,8))        
      plot_wavelet(csf_ortho_, 1)
      plt.axis('off')
      plt.show()
      print('#########################')
      print('Step 7')
      print('CSF filter test passed!!!')
      print('#########################') 
      I_thresha = iwavelet_transform(csf_ortho_)
      cv2_imshow((np.clip(I_thresha,0,1)*255).astype('uint8'))
      e=mse(I_thresh, I_thresha)
      print('MSE:{:.3f}'.format(e))
      print('############################')
      print('Recon CSF_filter test passed')
      print('############################')
      ############################################################################
      # Saturation of each sub-band of steerbale wavelets after CSF filter
      ############################################################################
      if saturation == True:
        sizeT_subL_db=wI.shape[0]
        xm_l_db = np.array(np.ones([256,256])).astype('float32') 
      
        with tf.Session() as sess:
          f_l_sn = sess.run(tf.cast(tf.squeeze(saturation_f(np.array(csf_ortho_).astype('float32'), g_sa, xm_l_db, epsilon_sa, sizeT_subL_db)), dtype=tf.float32))
        plt.figure(figsize = (8,8))        
        plot_wavelet(f_l_sn, 1)
        plt.axis('off')
        plt.show()
        print('#########################')
        print('Step 8')
        print('Saturation test passed!!!')
        print('#########################') 
        ############################################################################
        # Neural noise - 1) Gaussian  2) Possion 
        ############################################################################
        #Gaussian
        if noise_type=='Gaussian':
          sigma=.1
          f_gaussian_ortho = np.abs(f_l_sn) + sigma*np.random.standard_normal(f_l_sn.shape)
          #cv2_imshow((np.clip(f_gaussian_ortho, 0,1)*255).astype('uint8'))
          plt.figure(figsize = (8,8))        
          plot_wavelet(f_gaussian_ortho, 1)
          plt.axis('off')
          plt.show()
        
        #Possion
        if noise_type=='Possion':
          f_possion_ortho = np.random.poisson(np.abs(f_l_sn), f_l_sn.shape)
          #cv2_imshow((np.clip(f_possion_ortho, 0,1)*255).astype('uint8'))
          plt.figure(figsize = (8,8))        
          plot_wavelet(f_possion_ortho, 1)
          plt.axis('off')
          plt.show()
        
        print('####################################')
        print('Step 9')
        print('Simulate neural noise test passed!!!')
        print('####################################')

        print('+++++++++++++++++++++++++++++++++++++++')
        print('+++++++++++++++++++++++++++++++++++++++')
        print('Retina-LGN-V1 Visual Computing Finished')
        print('+++++++++++++++++++++++++++++++++++++++')
        print('+++++++++++++++++++++++++++++++++++++++')
      ############################################################################
      # DN 
      ############################################################################
      else:
        csf_ortho_dn=np.reshape(I_thresha, (256, 256, 1))
        DN_ortho = DivisiveNormalization(torch.Tensor([csf_ortho_dn.transpose((2,0,1))]),radius=9) 
        ret_ortho = DN_ortho[0].numpy().transpose((1,2,0)) 
        scaled_ret_ortho = (ret_ortho - ret_ortho.min())/(ret_ortho.max() - ret_ortho.min()) 
        scaled_ret_ortho = np.reshape(scaled_ret_ortho, (256, 256))
        plt.figure(figsize=(8,8))
        plt.imshow(scaled_ret_ortho, cmap='gray')
        plt.axis('off')
        plt.show()
        print('##########################')
        print('Step 9')
        print('DN model test passed!!!')
        DN_ortho_recon = iwavelet_transform(scaled_ret_ortho)
        cv2_imshow((np.clip(DN_ortho_recon,0,1)*255).astype('uint8'))
        print('##########################')
        print('Recon DN test passed')
        print('##########################')
        ############################################################################
        # Neural noise - 1) Gaussian  2) Possion 
        ############################################################################
        #Gaussian
        if noise_type=='Gaussian':
          sigma=.1
          f_gaussian_ortho = np.abs(DN_ortho_recon) + sigma*np.random.standard_normal(DN_ortho_recon.shape)
          cv2_imshow((np.clip(f_gaussian_ortho, 0,1)*255).astype('uint8'))

        #Possion
        if noise_type=='Possion':
          f_possion_ortho = np.random.poisson(np.abs(DN_ortho_recon), DN_ortho_recon.shape)
          cv2_imshow((np.clip(f_possion_ortho, 0,1)*255).astype('uint8'))
        
        print('####################################')
        print('Step 9')
        print('Simulate neural noise test passed!!!')
        print('####################################')

        print('+++++++++++++++++++++++++++++++++++++++')
        print('+++++++++++++++++++++++++++++++++++++++')
        print('Retina-LGN-V1 Visual Computing Finished')
        print('+++++++++++++++++++++++++++++++++++++++')
        print('+++++++++++++++++++++++++++++++++++++++')
      
  # 5) -- DWT
  ############################################################################
  # Build Pyramid with DWT
  ############################################################################
  if wavelets_type== 'DWT':
    max_lev = 3
    level=max_lev   
    shape=waber_img.shape
    c = pywt.wavedec2(waber_img, 'db2', mode='periodization', level=level)
    c[0] /= np.abs(c[0]).max()
    for detail_level in trange(level, desc='Running'):
        c[detail_level + 1] = [d/np.abs(d).max() for d in c[detail_level + 1]]

    arr, slices = pywt.coeffs_to_array(c)
    plt.figure(figsize=(8,6))
    plt.imshow(arr, cmap='gray')
    plt.axis('off')
    plt.tight_layout()
    plt.show()
    
    print('######################################')
    print('Step 6')
    print('Build pyramid with DWT test passed!!!')
    print('#####################################')
    ##################################################
    #Raw image entropy
    ##################################################
    if statis_wavlets==True:
      start=time.time()
      N=5
      colorIm=rgb_img
      greyIm=convert_to_luminance(colorIm)
      colorIm=np.array(colorIm)
      greyIm=np.array(greyIm)

      Sql=greyIm.shape
      Eql=np.array(greyIm)
      for row in trange(Sql[0], desc='Running'):
        for col in trange(Sql[1], desc='Running'):
          Lx=np.max([0,col-N])
          Ux=np.min([Sql[1],col+N])
          Ly=np.max([0,row-N])
          Uy=np.min([Sql[0],row+N])
          region=greyIm[Ly:Uy,Lx:Ux].flatten()
          Eql[row,col]=entropy(region)

      Img_init=[colorIm, greyIm, Eql]
      plt.figure(figsize=(12,10))
      plt.subplots_adjust(wspace=0.1)
      for i in tqdm(range(len(Img_init)), desc='Running'):
        plt.subplot(1,3,i+1)
        plt.imshow(Img_init[i], cmap=plt.cm.jet)
        plt.axis('off')
      plt.show()
      ##################################################
      # Wavelet energy map
      ##################################################
      Energ_levels=[]
      for i in tqdm(range(1, max_lev+1), desc='Running'):
        for j in tqdm(range(max_lev), desc='Running'):
          Energ_level=np.abs(c[i][j])**2
          Energ_levels.append(Energ_level)
      #---------------------------------------------------
      #Energy map with grid show
      #---------------------------------------------------
      Init_en=np.abs(c[0])**2
      Eng_grid=[Init_en, Energ_levels[:3], Energ_levels[3:6], Energ_levels[6:9]]

      Eng_grid[0] /= np.abs(Eng_grid[0]).max()
      for detail_level in trange(level, desc='Running'):
        Eng_grid[detail_level + 1] = [d/np.abs(d).max() for d in Eng_grid[detail_level + 1]]
      arr_, slices_ = pywt.coeffs_to_array(Eng_grid)
      plt.figure(figsize=(8,6))
      plt.imshow(arr_, cmap='gray')
      plt.axis('off')
      plt.tight_layout()
      plt.show()
      #---------------------------------------------------
      #Energy map with flatten map
      #---------------------------------------------------
      fin=Energ_levels[0]+Energ_levels[1]+Energ_levels[2]
      fina=Energ_levels[3]+Energ_levels[4]+Energ_levels[5]
      finb=Energ_levels[6]+Energ_levels[7]+Energ_levels[8]

      plt.figure(figsize=(12,10))
      plt.subplots_adjust(wspace=0.1)
      ener_q=[fin, fina, finb]
      for i in tqdm(range(len(ener_q)), desc='Running'):
        plt.subplot(1,3,i+1)
        plt.imshow(ener_q[i], cmap='gray')
        plt.axis('off')
      plt.show()
      #------------------------------------------------------------
      #Entropy map(sum of all different level) with 10*10 neigh 
      # IT EXIST BUG HERE, BE ATTENATION
      #------------------------------------------------------------
      Entropy_maps=[]
      for i in tqdm(range(len(Energ_levels)), desc='Running'):
        Sty=Energ_levels[i].shape
        Ety=np.array(Energ_levels[i])
        for row in tqdm(range(Sty[0]), desc='Running'):
          for col in tqdm(range(Sty[1]), desc='Running'):
            Lx=np.max([0,col-N])
            Ux=np.min([Sty[1],col+N])
            Ly=np.max([0,row-N])
            Uy=np.min([Sty[0],row+N])
            region=(Energ_levels[i])[Ly:Uy,Lx:Ux].flatten()
            Ety[row,col]=entropy(region)
            Entropy_maps.append(Ety)

      e_fin=Entropy_maps[0]+Entropy_maps[1]+Entropy_maps[2]
      e_fina=Entropy_maps[3]+Entropy_maps[4]+Entropy_maps[5]
      e_finb=Entropy_maps[6]+Entropy_maps[7]+Entropy_maps[8]

      plt.figure(figsize=(12,10))
      plt.subplots_adjust(wspace=0.1)
      e_ener_we=[e_fin, e_fina, e_finb]
      for i in tqdm(range(len(e_ener_we)), desc='Running'):
        plt.subplot(1,3,i+1)
        plt.imshow(e_ener_we[i], cmap=plt.cm.jet)
        plt.axis('off')
      plt.show()
      #---------------------------------------------------------
      # Entropy map with grid show
      #---------------------------------------------------------
      Entropy_levelss=[]
      for i in tqdm(range(1, max_lev+1), desc='Running'):
        for j in tqdm(range(max_lev), desc='Running'):
          eS=(c[i][j]).shape
          eE=np.array(c[i][j])
          for row in tqdm(range(eS[0]), desc='Running'):
            for col in trange(eS[1], desc='Running'):
              Lx=np.max([0,col-N])
              Ux=np.min([eS[1],col+N])
              Ly=np.max([0,row-N])
              Uy=np.min([eS[0],row+N])
              region=(c[i][j])[Ly:Uy,Lx:Ux].flatten()
              eE[row,col]=entropy(region)
              Entropy_levelss.append(eE)

      Init_entr=c[0]
      eSql=Init_entr.shape
      eEql=np.array(Init_entr)
      for row in trange(eSql[0], desc='Running'):
        for col in trange(eSql[1], desc='Running'):
          Lx=np.max([0,col-N])
          Ux=np.min([eSql[1],col+N])
          Ly=np.max([0,row-N])
          Uy=np.min([eSql[0],row+N])
          region=Init_entr[Ly:Uy,Lx:Ux].flatten()
          eEql[row,col]=entropy(region)

      eEng_grid=[eEql, Entropy_levelss[:3], Entropy_levelss[3:6], Entropy_levelss[6:9]]

      eEng_grid[0] /= np.abs(eEng_grid[0]).max()
      for detail_level in trange(level, desc='Running'):
        eEng_grid[detail_level + 1] = [d/np.abs(d).max() for d in eEng_grid[detail_level + 1]]
      arr_e, slices_e = pywt.coeffs_to_array(eEng_grid)
      plt.figure(figsize=(8,6))
      plt.imshow(arr_e, cmap=plt.cm.jet)
      plt.axis('off')
      plt.tight_layout()
      plt.show()

      elapsed=time.time()-start
      print('Total need Time:{:.3f}'.format(elapsed))
    ############################################################################
    # CSF  
    ############################################################################
    CSF_wavelets_vis_dwt=[]
    CSF_wavelet_l_dwt=[]
    CSF_wavelet_m_dwt=[]
    CSF_wavelet_h_dwt=[]
    CSF_wavelet_s_dwt=[]
    
    res_L_dwt=make_CSF(x=32, nfreq=32)
    Img_CSF_L_dwt = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(c[0].real))* res_L_dwt))))
    CSF_wavelet_l_dwt.append(Img_CSF_L_dwt)
    CSF_wavelets_vis_dwt.append(Img_CSF_L_dwt)
    for t_dwt in tqdm(range(0,3), desc='Running'):
      res_M_dwt=make_CSF(x=32, nfreq=32)
      Img_CSF_M_dwt = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(c[1][t_dwt].real))* res_M_dwt))))
      CSF_wavelet_m_dwt.append(Img_CSF_M_dwt)
      CSF_wavelets_vis_dwt.append(Img_CSF_M_dwt)
    for y_dwt in tqdm(range(0,3), desc='Running'):  
      res_H_dwt=make_CSF(x=64, nfreq=32)
      Img_CSF_H_dwt = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(c[2][y_dwt].real))* res_H_dwt))))
      CSF_wavelet_h_dwt.append(Img_CSF_H_dwt)
      CSF_wavelets_vis_dwt.append(Img_CSF_H_dwt)
    for e_dwt in tqdm(range(0,3), desc='Running'):  
      res_S_dwt=make_CSF(x=128, nfreq=32)
      Img_CSF_S_dwt = np.double (np.real(np.fft.ifft2(np.fft.ifftshift(np.fft.fftshift(np.fft.fft2(c[3][e_dwt].real))* res_S_dwt))))
      CSF_wavelet_s_dwt.append(Img_CSF_S_dwt)
      CSF_wavelets_vis_dwt.append(Img_CSF_S_dwt)

    CSF_wavelets_dwt=[CSF_wavelet_l_dwt[0], CSF_wavelet_m_dwt, CSF_wavelet_h_dwt, CSF_wavelet_s_dwt]

    if len(CSF_wavelets_dwt)==4:
      print('CSF filter was right!!')
    else:
      print('CSF filter test wrong, go back check!!')
    
    plt.figure(figsize=(10,10))
    plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
    plt.margins(0,0)
    for j in tqdm(range(len(CSF_wavelets_vis_dwt[1:])), desc='Running'):
      plt.subplot(3,3,j+1)
      plt.imshow(CSF_wavelets_vis_dwt[j], cmap='gray')
      plt.axis('off')
    plt.tight_layout()
    plt.show()
    #-----------------------------------------------------------------
    # CSF grid show
    #------------------------------------------------------------------
    CSF_wavelets_dwt[0] /= np.abs(CSF_wavelets_dwt[0]).max()
    for detail_level in trange(level, desc='Running'):
      CSF_wavelets_dwt[detail_level + 1] = [d/np.abs(d).max() for d in CSF_wavelets_dwt[detail_level + 1]]
    arr_dwt, slices_dwt = pywt.coeffs_to_array(CSF_wavelets_dwt)
    plt.figure(figsize=(8,6))
    plt.imshow(arr_dwt, cmap='gray')
    plt.axis('off')
    plt.tight_layout()
    plt.show()
    print('##############################################')
    print('Step 7')
    print('CSF filter wavelets filter test passed!!!')
    print('#############################################')
    
    #-----------------------------------------------------------------
    # CSF wavelets reconstructure
    #------------------------------------------------------------------
    CSF_idwt=pywt.waverec2(CSF_wavelets_dwt, 'db2', 'periodization')
    plt.figure(figsize=(4,6))
    plt.imshow(CSF_idwt, cmap='gray')
    plt.axis('off')
    plt.tight_layout()
    plt.show()
    print('##############################################')
    print('Step 8')
    print('CSF filter wavelets upsampling test passed!!!')
    print('#############################################')
    ############################################################################
    # DN 
    ############################################################################
    csf_dwt_dn=np.reshape(CSF_idwt,(256, 256, 1))
    DN_dwt = DivisiveNormalization(torch.Tensor([csf_dwt_dn.transpose((2,0,1))]),radius=9) 
    ret_dwt = DN_dwt[0].numpy().transpose((1,2,0)) 
    scaled_ret_dwt = (ret_dwt - ret_dwt.min())/(ret_dwt.max() - ret_dwt.min()) 
    scaled_ret_dwt = np.reshape(scaled_ret_dwt, (256, 256))
    plt.figure(figsize=(6,6))
    plt.imshow(scaled_ret_dwt, cmap='gray')
    plt.axis('off')
    plt.show()
    print('########################')
    print('Step 9')
    print('DN model test passed!!!')
    ############################################################################
    # Neural noise - 1) Gaussian  2) Possion 
    ############################################################################
    #Gaussian
    if noise_type=='Gaussian':
      sigma=.1
      f_gaussian_dwt = scaled_ret_dwt + sigma*np.random.standard_normal(scaled_ret_dwt.shape)
      viewimage((np.clip(f_gaussian_dwt, 0,1)*255).astype('uint8'))

    #Possion
    if noise_type=='Possion':
      f_possion_dwt = np.random.poisson(scaled_ret_dwt, scaled_ret_dwt.shape)
      viewimage((np.clip(f_possion_dwt, 0,1)*255).astype('uint8'))
    
    print('########################')
    print('Step 10')
    print('Simulate neural noise test passed!!!')
    print('########################')

    print('+++++++++++++++++++++++++++++++++++++++')
    print('+++++++++++++++++++++++++++++++++++++++')
    print('Retina-LGN-V1 Visual Computing Finished')
    print('+++++++++++++++++++++++++++++++++++++++')
    print('+++++++++++++++++++++++++++++++++++++++')