IIRLab

IIR-Lab/Dockerfile

@@ -0,0 +1,17 @@
+FROM pytorch/pytorch:2.0.1-cuda11.7-cudnn8-runtime
+
+ENV TZ=Asia
+ARG DEBIAN_FRONTEND=noninteractive
+
+RUN apt-get update && apt-get install -y \
+    libpng-dev libjpeg-dev \
+    libopencv-dev ffmpeg \
+    libgl1-mesa-glx 
+
+COPY requirements.txt .
+RUN python -m pip install --upgrade pip
+RUN pip install --no-cache -r requirements.txt
+
+COPY . /nightimage
+RUN chmod +x /nightimage/run.sh
+WORKDIR /nightimage

IIR-Lab/ISP_pipeline/.gitignore

@@ -0,0 +1,129 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/

IIR-Lab/ISP_pipeline/Dockerfile

@@ -0,0 +1,12 @@
+FROM python:3.9
+
+RUN apt-get update && apt-get install -y \
+    libpng-dev libjpeg-dev \
+    libopencv-dev ffmpeg \
+    libgl1-mesa-glx
+
+COPY requirements.txt .
+RUN python -m pip install --no-cache -r requirements.txt
+
+COPY . /nightimaging
+WORKDIR /nightimaging

IIR-Lab/ISP_pipeline/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023 Color Reproduction and Synthesis
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

IIR-Lab/ISP_pipeline/cfa_pattern_change.py

@@ -0,0 +1,74 @@
+import numpy as np
+import os
+import json
+import cv2
+
+def change_cfa_pattern(img):
+    raw_colors = np.asarray([3,1,2,0]).reshape((2, 2))
+    changed_raw_colors = np.asarray([0,1,2,3]).reshape((2, 2))
+    demosaiced_image = np.zeros((img.shape[0]//2, img.shape[1]//2, 4))
+    for i in range(2):
+        for j in range(2):
+            ch = raw_colors[i, j]
+            demosaiced_image[:, :, ch] = img[i::2, j::2]
+    for i in range(2):
+        for j in range(2):
+            ch1 = changed_raw_colors[i, j]
+            img[i::2, j::2] = demosaiced_image[:, :, ch1]
+
+    return img
+
+def rggb_raw(raw):
+    # pack RGGB Bayer raw to 4 channels
+    H, W = raw.shape
+    raw = raw[None, ...]
+    raw_pack = np.concatenate((raw[:, 0:H:2, 0:W:2],
+                               raw[:, 0:H:2, 1:W:2],
+                               raw[:, 1:H:2, 0:W:2],
+                               raw[:, 1:H:2, 1:W:2]), axis=0)
+    # tmp = rggb[...,0]
+    # rggb[...,0] = rggb[...,-1]
+    # rggb[...,-1] = tmp
+    return raw_pack
+
+def raw_rggb(raws):
+    # depack 4 channels raw to RGGB Bayer
+    C, H, W = raws.shape
+    output = np.zeros((H * 2, W * 2)).astype(np.uint16)
+
+    output[0:2 * H:2, 0:2 * W:2] = raws[0:1, :, :]
+    output[0:2 * H:2, 1:2 * W:2] = raws[1:2, :, :]
+    output[1:2 * H:2, 0:2 * W:2] = raws[2:3, :, :]
+    output[1:2 * H:2, 1:2 * W:2] = raws[3:4, :, :]
+
+    return output
+
+if __name__ == "__main__":
+    json_path = "/data1/02_data/Train_Data/"
+    file_name = os.listdir(json_path)
+    json_list = []
+    for file_name_all in file_name:
+        if file_name_all.endswith(".json"):
+            json_list.append(json_path+file_name_all)
+    a = []
+    for i in range(len(json_list)):
+        with open(json_list[i],'r',encoding='UTF-8') as f:
+            result = json.load(f)
+            # a,b = result["noise_profile"]
+            # black = result["white_level"]
+            cfa_pattern = result["cfa_pattern"]
+            if cfa_pattern[0] == 2:
+                a.append(json_list[i])
+    for j in range(len(a)):
+        pic_name,_ = os.path.splitext(a[j])
+        img = cv2.imread(pic_name+str(".png"), cv2.IMREAD_UNCHANGED)
+        # img1 = cv2.imread(pic_name+".png", cv2.IMREAD_UNCHANGED)
+        # test = img - img1
+        # print(test)
+        changed_img = change_cfa_pattern(img=img)
+        # cv2.imwrite(pic_name+"test1.png",changed_img)
+        np.save(pic_name+"origin.npy",img)
+        np.save(pic_name+"changed.npy",changed_img)
+    # np.save("./json_all.npy",result)
+    
+    

IIR-Lab/ISP_pipeline/debayer.py

@@ -0,0 +1,1295 @@
+import numpy as np
+import math
+import time
+import utility
+from scipy import signal
+
+# =============================================================
+# function: dbayer_mhc
+#   demosaicing using Malvar-He-Cutler algorithm
+#   http://www.ipol.im/pub/art/2011/g_mhcd/
+# =============================================================
+def debayer_mhc(raw, bayer_pattern="rggb", clip_range=[0, 65535], timeshow=False):
+
+    # convert to float32 in case it was not
+    raw = np.float32(raw)
+
+    # dimensions
+    width, height = utility.helpers(raw).get_width_height()
+
+    # number of pixels to pad
+    no_of_pixel_pad = 2
+    raw = np.pad(raw, \
+                 (no_of_pixel_pad, no_of_pixel_pad),\
+                 'reflect') # reflect would not repeat the border value
+
+    # allocate space for the R, G, B planes
+    R = np.empty( (height + no_of_pixel_pad * 2, width + no_of_pixel_pad * 2), dtype = np.float32 )
+    G = np.empty( (height + no_of_pixel_pad * 2, width + no_of_pixel_pad * 2), dtype = np.float32 )
+    B = np.empty( (height + no_of_pixel_pad * 2, width + no_of_pixel_pad * 2), dtype = np.float32 )
+
+    # create a RGB output
+    demosaic_out = np.empty( (height, width, 3), dtype = np.float32 )
+
+    # fill up the directly available values according to the Bayer pattern
+    if (bayer_pattern == "rggb"):
+
+        G[::2, 1::2]  = raw[::2, 1::2]
+        G[1::2, ::2]  = raw[1::2, ::2]
+        R[::2, ::2]   = raw[::2, ::2]
+        B[1::2, 1::2] = raw[1::2, 1::2]
+
+        # Green channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # G at Red location
+                if (((i % 2) == 0) and ((j % 2) == 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+                    2. * G[i+1, j], \
+                    -1. * R[i+2, j]])
+                # G at Blue location
+                elif (((i % 2) != 0) and ((j % 2) != 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+                    2. * G[i+1, j],\
+                    -1. * B[i+2, j]])
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Green: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+        # Red and Blue channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # Green locations in Red rows
+                if (((i % 2) == 0) and ((j % 2) != 0)):
+                    # R at Green locations in Red rows
+                    R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+                     -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+                     -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                      .5 * G[i+2, j]])
+
+                    # B at Green locations in Red rows
+                    B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                # Green locations in Blue rows
+                elif (((i % 2) != 0) and ((j % 2) == 0)):
+
+                    # R at Green locations in Blue rows
+                    R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                    # B at Green locations in Blue rows
+                    B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+                    -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+                    -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+                    -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                    .5 * G[i+2, j]])
+
+                # R at Blue locations
+                elif (((i % 2) != 0) and ((j % 2) != 0)):
+                    R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+                    2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+                    -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+                    2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+                    -1.5 * B[i+2, j]])
+
+                # B at Red locations
+                elif (((i % 2) == 0) and ((j % 2) == 0)):
+                    B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+                    2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+                    -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+                    2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+                    -1.5 * R[i+2, j]])
+
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+
+    elif (bayer_pattern == "gbrg"):
+
+        G[::2, ::2]   = raw[::2, ::2]
+        G[1::2, 1::2] = raw[1::2, 1::2]
+        R[1::2, ::2]  = raw[1::2, ::2]
+        B[::2, 1::2]  = raw[::2, 1::2]
+
+        # Green channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # G at Red location
+                if (((i % 2) != 0) and ((j % 2) == 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+                    2. * G[i+1, j], \
+                    -1. * R[i+2, j]])
+                # G at Blue location
+                elif (((i % 2) == 0) and ((j % 2) != 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+                    2. * G[i+1, j],\
+                    -1. * B[i+2, j]])
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Green: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+        # Red and Blue channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # Green locations in Red rows
+                if (((i % 2) != 0) and ((j % 2) != 0)):
+                    # R at Green locations in Red rows
+                    R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+                     -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+                     -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                      .5 * G[i+2, j]])
+
+                    # B at Green locations in Red rows
+                    B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                # Green locations in Blue rows
+                elif (((i % 2) == 0) and ((j % 2) == 0)):
+
+                    # R at Green locations in Blue rows
+                    R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                    # B at Green locations in Blue rows
+                    B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+                    -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+                    -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+                    -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                    .5 * G[i+2, j]])
+
+                # R at Blue locations
+                elif (((i % 2) == 0) and ((j % 2) != 0)):
+                    R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+                    2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+                    -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+                    2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+                    -1.5 * B[i+2, j]])
+
+                # B at Red locations
+                elif (((i % 2) != 0) and ((j % 2) == 0)):
+                    B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+                    2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+                    -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+                    2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+                    -1.5 * R[i+2, j]])
+
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+    elif (bayer_pattern == "grbg"):
+
+        G[::2, ::2]   = raw[::2, ::2]
+        G[1::2, 1::2] = raw[1::2, 1::2]
+        R[::2, 1::2]  = raw[::2, 1::2]
+        B[1::2, ::2]  = raw[1::2, ::2]
+
+        # Green channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # G at Red location
+                if (((i % 2) == 0) and ((j % 2) != 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+                    2. * G[i+1, j], \
+                    -1. * R[i+2, j]])
+                # G at Blue location
+                elif (((i % 2) != 0) and ((j % 2) == 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+                    2. * G[i+1, j],\
+                    -1. * B[i+2, j]])
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Green: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+        # Red and Blue channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # Green locations in Red rows
+                if (((i % 2) == 0) and ((j % 2) == 0)):
+                    # R at Green locations in Red rows
+                    R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+                     -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+                     -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                      .5 * G[i+2, j]])
+
+                    # B at Green locations in Red rows
+                    B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                # Green locations in Blue rows
+                elif (((i % 2) != 0) and ((j % 2) != 0)):
+
+                    # R at Green locations in Blue rows
+                    R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                    # B at Green locations in Blue rows
+                    B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+                    -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+                    -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+                    -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                    .5 * G[i+2, j]])
+
+                # R at Blue locations
+                elif (((i % 2) != 0) and ((j % 2) == 0)):
+                    R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+                    2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+                    -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+                    2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+                    -1.5 * B[i+2, j]])
+
+                # B at Red locations
+                elif (((i % 2) == 0) and ((j % 2) != 0)):
+                    B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+                    2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+                    -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+                    2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+                    -1.5 * R[i+2, j]])
+
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+    elif (bayer_pattern == "bggr"):
+
+        G[::2, 1::2]  = raw[::2, 1::2]
+        G[1::2, ::2]  = raw[1::2, ::2]
+        R[1::2, 1::2] = raw[1::2, 1::2]
+        B[::2, ::2]   = raw[::2, ::2]
+
+        # Green channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # G at Red location
+                if (((i % 2) != 0) and ((j % 2) != 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+                    2. * G[i+1, j], \
+                    -1. * R[i+2, j]])
+                # G at Blue location
+                elif (((i % 2) == 0) and ((j % 2) == 0)):
+                    G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+                    2. * G[i-1, j], \
+                    -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+                    2. * G[i+1, j],\
+                    -1. * B[i+2, j]])
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Green: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+        # Red and Blue channel
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+
+            # to display progress
+            t0 = time.process_time()
+
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # Green locations in Red rows
+                if (((i % 2) != 0) and ((j % 2) == 0)):
+                    # R at Green locations in Red rows
+                    R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+                     -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+                     -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                      .5 * G[i+2, j]])
+
+                    # B at Green locations in Red rows
+                    B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                # Green locations in Blue rows
+                elif (((i % 2) == 0) and ((j % 2) != 0)):
+
+                    # R at Green locations in Blue rows
+                    R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+                    -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+                    .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+                    -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+                    -1. * G[i+2, j]])
+
+                    # B at Green locations in Blue rows
+                    B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+                    -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+                    -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+                    -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+                    .5 * G[i+2, j]])
+
+                # R at Blue locations
+                elif (((i % 2) == 0) and ((j % 2) == 0)):
+                    R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+                    2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+                    -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+                    2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+                    -1.5 * B[i+2, j]])
+
+                # B at Red locations
+                elif (((i % 2) != 0) and ((j % 2) != 0)):
+                    B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+                    2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+                    -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+                    2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+                    -1.5 * R[i+2, j]])
+
+            if (timeshow):
+                elapsed_time = time.process_time() - t0
+                print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+                      " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+
+    else:
+        print("Invalid bayer pattern. Valid pattern can be rggb, gbrg, grbg, bggr")
+        return demosaic_out # This will be all zeros
+
+    # Fill up the RGB output with interpolated values
+    demosaic_out[0:height, 0:width, 0] = R[no_of_pixel_pad : height + no_of_pixel_pad, \
+                                           no_of_pixel_pad : width + no_of_pixel_pad]
+    demosaic_out[0:height, 0:width, 1] = G[no_of_pixel_pad : height + no_of_pixel_pad, \
+                                           no_of_pixel_pad : width + no_of_pixel_pad]
+    demosaic_out[0:height, 0:width, 2] = B[no_of_pixel_pad : height + no_of_pixel_pad, \
+                                           no_of_pixel_pad : width + no_of_pixel_pad]
+
+    demosaic_out = np.clip(demosaic_out, clip_range[0], clip_range[1])
+    return demosaic_out
+
+
+def fill_channel_directional_weight(data, bayer_pattern):
+
+    #== Calculate the directional weights (weight_N, weight_E, weight_S, weight_W.
+    # where N, E, S, W stand for north, east, south, and west.)
+    data = np.asarray(data)
+    v = np.asarray(signal.convolve2d(data, [[1],[0],[-1]], mode="same", boundary="symm"))
+    h = np.asarray(signal.convolve2d(data, [[1, 0, -1]], mode="same", boundary="symm"))
+
+    weight_N = np.zeros(np.shape(data), dtype=np.float32)
+    weight_E = np.zeros(np.shape(data), dtype=np.float32)
+    weight_S = np.zeros(np.shape(data), dtype=np.float32)
+    weight_W = np.zeros(np.shape(data), dtype=np.float32)
+
+    value_N = np.zeros(np.shape(data), dtype=np.float32)
+    value_E = np.zeros(np.shape(data), dtype=np.float32)
+    value_S = np.zeros(np.shape(data), dtype=np.float32)
+    value_W = np.zeros(np.shape(data), dtype=np.float32)
+
+    if ((bayer_pattern == "rggb") or (bayer_pattern == "bggr")):
+
+
+        # note that in the following the locations in the comments are given
+        # assuming the bayer_pattern rggb
+
+        #== CALCULATE WEIGHTS IN B LOCATIONS
+        weight_N[1::2, 1::2] = np.abs(v[1::2, 1::2]) + np.abs(v[::2, 1::2])
+
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp_h_b = np.hstack((h, np.atleast_2d(h[:, -2]).T))
+        weight_E[1::2, 1::2] = np.abs(h[1::2, 1::2]) + np.abs(temp_h_b[1::2, 2::2])
+
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp_v_b = np.vstack((v, v[-1]))
+        weight_S[1::2, 1::2] = np.abs(v[1::2, 1::2]) + np.abs(temp_v_b[2::2, 1::2])
+        weight_W[1::2, 1::2] = np.abs(h[1::2, 1::2]) + np.abs(h[1::2, ::2])
+
+        #== CALCULATE WEIGHTS IN R LOCATIONS
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp_v_r = np.delete(np.vstack((v[1], v)), -1, 0)
+        weight_N[::2, ::2] = np.abs(v[::2, ::2]) + np.abs(temp_v_r[::2, ::2])
+
+        weight_E[::2, ::2] = np.abs(h[::2, ::2]) + np.abs(h[::2, 1::2])
+
+        weight_S[::2, ::2] = np.abs(v[::2, ::2]) + np.abs(v[1::2, ::2])
+
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp_h_r = np.delete(np.hstack((np.atleast_2d(h[:, 1]).T, h)), -1, 1)
+        weight_W[::2, ::2] = np.abs(h[::2, ::2]) + np.abs(temp_h_r[::2, ::2])
+
+        weight_N = np.divide(1., 1. + weight_N)
+        weight_E = np.divide(1., 1. + weight_E)
+        weight_S = np.divide(1., 1. + weight_S)
+        weight_W = np.divide(1., 1. + weight_W)
+
+        #== CALCULATE DIRECTIONAL ESTIMATES IN B LOCATIONS
+        value_N[1::2, 1::2] = data[::2, 1::2] + v[::2, 1::2] / 2.
+
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((data, np.atleast_2d(data[:, -2]).T))
+        value_E[1::2, 1::2] = temp[1::2, 2::2] - temp_h_b[1::2, 2::2] / 2.
+
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp = np.vstack((data, data[-1]))
+        value_S[1::2, 1::2] = temp[2::2, 1::2] - temp_v_b[2::2, 1::2] / 2.
+
+        value_W[1::2, 1::2] = data[1::2, ::2] + h[1::2, ::2] / 2.
+
+        #== CALCULATE DIRECTIONAL ESTIMATES IN R LOCATIONS
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp = np.delete(np.vstack((data[1], data)), -1, 0)
+        value_N[::2, ::2] = temp[::2, ::2] + temp_v_r[::2, ::2] / 2.
+
+        value_E[::2, ::2] = data[::2, 1::2] - h[::2, 1::2] / 2.
+
+        value_S[::2, ::2] = data[1::2, ::2] - v[1::2, ::2] / 2.
+
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp = np.delete(np.hstack((np.atleast_2d(data[:, 1]).T, data)), -1, 1)
+        value_W[::2, ::2] = temp[::2, ::2] + temp_h_r[::2, ::2] / 2.
+
+        output = np.zeros(np.shape(data), dtype=np.float32)
+        output = np.divide((np.multiply(value_N, weight_N) + \
+                            np.multiply(value_E, weight_E) + \
+                            np.multiply(value_S, weight_S) + \
+                            np.multiply(value_W, weight_W)),\
+                            (weight_N + weight_E + weight_S + weight_W))
+
+        output[::2, 1::2] = data[::2, 1::2]
+        output[1::2, ::2] = data[1::2, ::2]
+
+        return output
+
+    elif ((bayer_pattern == "gbrg") or (bayer_pattern == "grbg")):
+
+        # note that in the following the locations in the comments are given
+        # assuming the bayer_pattern gbrg
+
+        #== CALCULATE WEIGHTS IN B LOCATIONS
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp_v_b = np.delete(np.vstack((v[1], v)), -1, 0)
+        weight_N[::2, 1::2] = np.abs(v[::2, 1::2]) + np.abs(temp_v_b[::2, 1::2])
+
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp_h_b = np.hstack((h, np.atleast_2d(h[:, -2]).T))
+        weight_E[::2, 1::2] = np.abs(h[::2, 1::2]) + np.abs(temp_h_b[::2, 2::2])
+
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        weight_S[::2, 1::2] = np.abs(v[::2, 1::2]) + np.abs(v[1::2, 1::2])
+        weight_W[::2, 1::2] = np.abs(h[::2, 1::2]) + np.abs(h[::2, ::2])
+
+        #== CALCULATE WEIGHTS IN R LOCATIONS
+        weight_N[1::2, ::2] = np.abs(v[1::2, ::2]) + np.abs(v[::2, ::2])
+        weight_E[1::2, ::2] = np.abs(h[1::2, ::2]) + np.abs(h[1::2, 1::2])
+
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp_v_r = np.vstack((v, v[-1]))
+        weight_S[1::2, ::2] = np.abs(v[1::2, ::2]) + np.abs(temp_v_r[2::2, ::2])
+
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp_h_r = np.delete(np.hstack((np.atleast_2d(h[:, 1]).T, h)), -1, 1)
+        weight_W[1::2, ::2] = np.abs(h[1::2, ::2]) + np.abs(temp_h_r[1::2, ::2])
+
+        weight_N = np.divide(1., 1. + weight_N)
+        weight_E = np.divide(1., 1. + weight_E)
+        weight_S = np.divide(1., 1. + weight_S)
+        weight_W = np.divide(1., 1. + weight_W)
+
+        #== CALCULATE DIRECTIONAL ESTIMATES IN B LOCATIONS
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp = np.delete(np.vstack((data[1], data)), -1, 0)
+        value_N[::2, 1::2] = temp[::2, 1::2] + temp_v_b[::2, 1::2] / 2.
+
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((data, np.atleast_2d(data[:, -2]).T))
+        value_E[::2, 1::2] = temp[::2, 2::2] - temp_h_b[::2, 2::2] / 2.
+
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        value_S[::2, 1::2] = data[1::2, 1::2] - v[1::2, 1::2] / 2.
+
+        value_W[::2, 1::2] = data[::2, ::2] + h[::2, ::2] / 2.
+
+        #== CALCULATE DIRECTIONAL ESTIMATES IN R LOCATIONS
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        value_N[1::2, ::2] = data[::2, ::2] + v[::2, ::2] / 2.
+        value_E[1::2, ::2] = data[1::2, 1::2] - h[1::2, 1::2] / 2.
+
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp = np.vstack((data, data[-1]))
+        value_S[1::2, ::2] = temp[2::2, ::2] - temp_v_r[2::2, ::2] / 2.
+
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp = np.delete(np.hstack((np.atleast_2d(data[:, 1]).T, data)), -1, 1)
+        value_W[1::2, ::2] = temp[1::2, ::2] + temp_h_r[1::2, ::2] / 2.
+
+        output = np.zeros(np.shape(data), dtype=np.float32)
+        output = np.divide((np.multiply(value_N, weight_N) + \
+                            np.multiply(value_E, weight_E) + \
+                            np.multiply(value_S, weight_S) + \
+                            np.multiply(value_W, weight_W)),\
+                            (weight_N + weight_E + weight_S + weight_W))
+
+        output[::2, ::2] = data[::2, ::2]
+        output[1::2, 1::2] = data[1::2, 1::2]
+
+        return output
+
+
+def fill_br_locations(data, G, bayer_pattern):
+
+    # Fill up the B/R values interpolated at R/B locations
+    B = np.zeros(np.shape(data), dtype=np.float32)
+    R = np.zeros(np.shape(data), dtype=np.float32)
+
+    data = np.asarray(data)
+    G = np.asarray(G)
+    d1 = np.asarray(signal.convolve2d(data, [[-1, 0, 0],[0, 0, 0], [0, 0, 1]], mode="same", boundary="symm"))
+    d2 = np.asarray(signal.convolve2d(data, [[0, 0, 1], [0, 0, 0], [-1, 0, 0]], mode="same", boundary="symm"))
+
+    df_NE = np.asarray(signal.convolve2d(G, [[0, 0, 0], [0, 1, 0], [-1, 0, 0]], mode="same", boundary="symm"))
+    df_SE = np.asarray(signal.convolve2d(G, [[-1, 0, 0], [0, 1, 0], [0, 0, 0]], mode="same", boundary="symm"))
+    df_SW = np.asarray(signal.convolve2d(G, [[0, 0, -1], [0, 1, 0], [0, 0, 0]], mode="same", boundary="symm"))
+    df_NW = np.asarray(signal.convolve2d(G, [[0, 0, 0], [0, 1, 0], [0, 0, -1]], mode="same", boundary="symm"))
+
+    weight_NE = np.zeros(np.shape(data), dtype=np.float32)
+    weight_SE = np.zeros(np.shape(data), dtype=np.float32)
+    weight_SW = np.zeros(np.shape(data), dtype=np.float32)
+    weight_NW = np.zeros(np.shape(data), dtype=np.float32)
+
+    value_NE = np.zeros(np.shape(data), dtype=np.float32)
+    value_SE = np.zeros(np.shape(data), dtype=np.float32)
+    value_SW = np.zeros(np.shape(data), dtype=np.float32)
+    value_NW = np.zeros(np.shape(data), dtype=np.float32)
+
+    if ((bayer_pattern == "rggb") or (bayer_pattern == "bggr")):
+
+        #== weights for B in R locations
+        weight_NE[::2, ::2] = np.abs(d2[::2, ::2]) + np.abs(df_NE[::2, ::2])
+        weight_SE[::2, ::2] = np.abs(d1[::2, ::2]) + np.abs(df_SE[::2, ::2])
+        weight_SW[::2, ::2] = np.abs(d2[::2, ::2]) + np.abs(df_SW[::2, ::2])
+        weight_NW[::2, ::2] = np.abs(d1[::2, ::2]) + np.abs(df_NW[::2, ::2])
+
+        #== weights for R in B locations
+        weight_NE[1::2, 1::2] = np.abs(d2[1::2, 1::2]) + np.abs(df_NE[1::2, 1::2])
+        weight_SE[1::2, 1::2] = np.abs(d1[1::2, 1::2]) + np.abs(df_SE[1::2, 1::2])
+        weight_SW[1::2, 1::2] = np.abs(d2[1::2, 1::2]) + np.abs(df_SW[1::2, 1::2])
+        weight_NW[1::2, 1::2] = np.abs(d1[1::2, 1::2]) + np.abs(df_NW[1::2, 1::2])
+
+        weight_NE = np.divide(1., 1. + weight_NE)
+        weight_SE = np.divide(1., 1. + weight_SE)
+        weight_SW = np.divide(1., 1. + weight_SW)
+        weight_NW = np.divide(1., 1. + weight_NW)
+
+        #== directional estimates of B in R locations
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp = np.delete(np.vstack((data[1], data)), -1, 0)
+        value_NE[::2, ::2] = temp[::2, 1::2] + df_NE[::2, ::2] / 2.
+        value_SE[::2, ::2] = data[1::2, 1::2] + df_SE[::2, ::2] / 2.
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp = np.delete(np.hstack((np.atleast_2d(data[:, 1]).T, data)), -1, 1)
+        value_SW[::2, ::2] = temp[1::2, ::2] + df_SW[::2, ::2] / 2.
+
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp = np.delete(np.vstack((data[1], data)), -1, 0)
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp = np.delete(np.hstack((np.atleast_2d(temp[:, 1]).T, temp)), -1, 1)
+        value_NW[::2, ::2] = temp[::2, ::2] + df_NW[::2, ::2]
+
+        #== directional estimates of R in B locations
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((data, np.atleast_2d(data[:, -2]).T))
+        value_NE[1::2, 1::2] = temp[::2, 2::2] + df_NE[1::2, 1::2] / 2.
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((data, np.atleast_2d(data[:, -2]).T))
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp = np.vstack((temp, temp[-1]))
+        value_SE[1::2, 1::2] = temp[2::2, 2::2] + df_SE[1::2, 1::2] / 2.
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp = np.vstack((data, data[-1]))
+        value_SW[1::2, 1::2] = temp[2::2, ::2] + df_SW[1::2, 1::2] / 2.
+        value_NW[1::2, 1::2] = data[::2, ::2] + df_NW[1::2, 1::2] / 2.
+
+        RB = np.divide(np.multiply(weight_NE, value_NE) + \
+                       np.multiply(weight_SE, value_SE) + \
+                       np.multiply(weight_SW, value_SW) + \
+                       np.multiply(weight_NW, value_NW),\
+                       (weight_NE + weight_SE + weight_SW + weight_NW))
+
+        if (bayer_pattern == "rggb"):
+
+            R[1::2, 1::2] = RB[1::2, 1::2]
+            R[::2, ::2] = data[::2, ::2]
+            B[::2, ::2] = RB[::2, ::2]
+            B[1::2, 1::2] = data[1::2, 1::2]
+
+        elif (bayer_pattern == "bggr"):
+            R[::2, ::2] = RB[::2, ::2]
+            R[1::2, 1::2] = data[1::2, 1::2]
+            B[1::2, 1::2] = RB[1::2, 1::2]
+            B[::2, ::2] = data[::2, ::2]
+
+
+        R[1::2, ::2] = G[1::2, ::2]
+        R[::2, 1::2] = G[::2, 1::2]
+        R = fill_channel_directional_weight(R, "gbrg")
+
+        B[1::2, ::2] = G[1::2, ::2]
+        B[::2, 1::2] = G[::2, 1::2]
+        B = fill_channel_directional_weight(B, "gbrg")
+
+
+    elif ((bayer_pattern == "grbg") or (bayer_pattern == "gbrg")):
+        #== weights for B in R locations
+        weight_NE[::2, 1::2] = np.abs(d2[::2, 1::2]) + np.abs(df_NE[::2, 1::2])
+        weight_SE[::2, 1::2] = np.abs(d1[::2, 1::2]) + np.abs(df_SE[::2, 1::2])
+        weight_SW[::2, 1::2] = np.abs(d2[::2, 1::2]) + np.abs(df_SW[::2, 1::2])
+        weight_NW[::2, 1::2] = np.abs(d1[::2, 1::2]) + np.abs(df_NW[::2, 1::2])
+
+        #== weights for R in B locations
+        weight_NE[1::2, ::2] = np.abs(d2[1::2, ::2]) + np.abs(df_NE[1::2, ::2])
+        weight_SE[1::2, ::2] = np.abs(d1[1::2, ::2]) + np.abs(df_SE[1::2, ::2])
+        weight_SW[1::2, ::2] = np.abs(d2[1::2, ::2]) + np.abs(df_SW[1::2, ::2])
+        weight_NW[1::2, ::2] = np.abs(d1[1::2, ::2]) + np.abs(df_NW[1::2, ::2])
+
+        weight_NE = np.divide(1., 1. + weight_NE)
+        weight_SE = np.divide(1., 1. + weight_SE)
+        weight_SW = np.divide(1., 1. + weight_SW)
+        weight_NW = np.divide(1., 1. + weight_NW)
+
+        #== directional estimates of B in R locations
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp = np.delete(np.vstack((data[1], data)), -1, 0)
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((temp, np.atleast_2d(temp[:, -2]).T))
+        value_NE[::2, 1::2] = temp[::2, 2::2] + df_NE[::2, 1::2] / 2.
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((data, np.atleast_2d(data[:, -2]).T))
+        value_SE[::2, 1::2] = temp[1::2, 2::2] + df_SE[::2, 1::2] / 2.
+        value_SW[::2, 1::2] = data[1::2, ::2] + df_SW[::2, 1::2] / 2.
+
+        # repeating the second row at the top of matrix so that sampling does
+        # not cause any dimension mismatch, also remove the bottom row
+        temp = np.delete(np.vstack((data[1], data)), -1, 0)
+        value_NW[::2, 1::2] = temp[::2, ::2] + df_NW[::2, 1::2]
+
+        #== directional estimates of R in B locations
+        value_NE[1::2, ::2] = data[::2, 1::2] + df_NE[1::2, ::2] / 2.
+        # repeating the column before the last to the right so that sampling
+        # does not cause any dimension mismatch
+        temp = np.hstack((data, np.atleast_2d(data[:, -2]).T))
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp = np.vstack((temp, temp[-1]))
+        value_SE[1::2, ::2] = temp[2::2, 1::2] + df_SE[1::2, ::2] / 2.
+        # repeating the row before the last row to the bottom so that sampling
+        # does not cause any dimension mismatch
+        temp = np.vstack((data, data[-1]))
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp = np.delete(np.hstack((np.atleast_2d(temp[:, 1]).T, temp)), -1, 1)
+        value_SW[1::2, ::2] = temp[2::2, ::2] + df_SW[1::2, ::2] / 2.
+        # repeating the second column at the left of matrix so that sampling
+        # does not cause any dimension mismatch, also remove the rightmost
+        # column
+        temp = np.delete(np.hstack((np.atleast_2d(data[:, 1]).T, data)), -1, 1)
+        value_NW[1::2, ::2] = temp[::2, ::2] + df_NW[1::2, ::2] / 2.
+
+        RB = np.divide(np.multiply(weight_NE, value_NE) + \
+                       np.multiply(weight_SE, value_SE) + \
+                       np.multiply(weight_SW, value_SW) + \
+                       np.multiply(weight_NW, value_NW),\
+                       (weight_NE + weight_SE + weight_SW + weight_NW))
+
+        if (bayer_pattern == "grbg"):
+
+            R[1::2, ::2] = RB[1::2, ::2]
+            R[::2, 1::2] = data[::2, 1::2]
+            B[::2, 1::2] = RB[::2, 1::2]
+            B[1::2, ::2] = data[1::2, ::2]
+
+        elif (bayer_pattern == "gbrg"):
+            R[::2, 1::2] = RB[::2, 1::2]
+            R[1::2, ::2] = data[1::2, ::2]
+            B[1::2, ::2] = RB[1::2, ::2]
+            B[::2, 1::2] = data[::2, 1::2]
+
+
+        R[::2, ::2] = G[::2, ::2]
+        R[1::2, 1::2] = G[1::2, 1::2]
+        R = fill_channel_directional_weight(R, "rggb")
+
+        B[1::2, 1::2] = G[1::2, 1::2]
+        B[::2, ::2] = G[::2, ::2]
+        B = fill_channel_directional_weight(B, "rggb")
+
+
+    return B, R
+
+# # =============================================================
+# # function: dbayer_mhc_fast
+# #   demosaicing using Malvar-He-Cutler algorithm
+# #   http://www.ipol.im/pub/art/2011/g_mhcd/
+# # =============================================================
+# def debayer_mhc_fast(raw, bayer_pattern="rggb", clip_range=[0, 65535], timeshow=False):
+#
+#     # convert to float32 in case it was not
+#     raw = np.float32(raw)
+#
+#     # dimensions
+#     width, height = utility.helpers(raw).get_width_height()
+#
+#     # allocate space for the R, G, B planes
+#     R = np.empty((height, width), dtype = np.float32)
+#     G = np.empty((height, width), dtype = np.float32)
+#     B = np.empty((height, width), dtype = np.float32)
+#
+#     # create a RGB output
+#     demosaic_out = np.empty( (height, width, 3), dtype = np.float32 )
+#
+#     # define the convolution kernels
+#     kernel_g_at_rb = [[0., 0., -1., 0., 0.],\
+#                       [0., 0., 2., 0., 0.],\
+#                       [-1., 2., 4., 2., -1.],\
+#                       [0., 0., 2., 0., 0.],\
+#                       [0., 0., -1., 0., 0.]] * .125
+#
+#     kernel_r_at_gr = [[0., 0., .5, 0., 0.],\
+#                       [0., -1., 0., -1., 0.],\
+#                       [-1., 4., 5., 4., -1.],\
+#                       [0., -1., 0., -1., 0.],\
+#                       [0., 0., .5, 0., 0.]] * .125
+#
+#     kernel_b_at_gr = [[0., 0., -1., 0., 0.],\
+#                       [0., -1., 4., -1., 0.],\
+#                       [.5., 0., 5., 0., .5],\
+#                       [0., -1., 4., -1., 0],\
+#                       [0., 0., -1., 0., 0.]] * .125
+#
+#     kernel_r_at_gb = [[0., 0., -1., 0., 0.],\
+#                       [0., -1., 4., -1., 0.],\
+#                       [.5, 0., 5., 0., .5],\
+#                       [0., -1., 4., -1., 0.],\
+#                       [0., 0., -1., 0., 0.]] * .125
+#
+#     kernel_b_at_gb = [[0., 0., .5, 0., 0.],\
+#                       [0., -1., 0., -1., 0.],\
+#                       [-1., 4., 5., 4., -1.],\
+#                       [0., -1., 0., -1., 0.],\
+#                       [0., 0., .5, 0., 0.]] * .125
+#
+#     kernel_r_at_b = [[0., 0., -1.5, 0., 0.],\
+#                      [0., 2., 0., 2., 0.],\
+#                      [-1.5, 0., 6., 0., -1.5],\
+#                      [0., 2., 0., 2., 0.],\
+#                      [0., 0., -1.5, 0., 0.]] * .125
+#
+#     kernel_b_at_r = [[0., 0., -1.5, 0., 0.],\
+#                      [0., 2., 0., 2., 0.],\
+#                      [-1.5, 0., 6., 0., -1.5],\
+#                      [0., 2., 0., 2., 0.],\
+#                      [0., 0., -1.5, 0., 0.]] * .125
+#
+#
+#
+#     # fill up the directly available values according to the Bayer pattern
+#     if (bayer_pattern == "rggb"):
+#
+#         G[::2, 1::2]  = raw[::2, 1::2]
+#         G[1::2, ::2]  = raw[1::2, ::2]
+#         R[::2, ::2]   = raw[::2, ::2]
+#         B[1::2, 1::2] = raw[1::2, 1::2]
+#
+#         # Green channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # G at Red location
+#                 if (((i % 2) == 0) and ((j % 2) == 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+#                     2. * G[i+1, j], \
+#                     -1. * R[i+2, j]])
+#                 # G at Blue location
+#                 elif (((i % 2) != 0) and ((j % 2) != 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+#                     2. * G[i+1, j],\
+#                     -1. * B[i+2, j]])
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Green: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#         # Red and Blue channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # Green locations in Red rows
+#                 if (((i % 2) == 0) and ((j % 2) != 0)):
+#                     # R at Green locations in Red rows
+#                     R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+#                      -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+#                      -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+#                      -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                       .5 * G[i+2, j]])
+#
+#                     # B at Green locations in Red rows
+#                     B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                 # Green locations in Blue rows
+#                 elif (((i % 2) != 0) and ((j % 2) == 0)):
+#
+#                     # R at Green locations in Blue rows
+#                     R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                     # B at Green locations in Blue rows
+#                     B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+#                     -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+#                     -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+#                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                     .5 * G[i+2, j]])
+#
+#                 # R at Blue locations
+#                 elif (((i % 2) != 0) and ((j % 2) != 0)):
+#                     R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+#                     2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+#                     -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+#                     2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+#                     -1.5 * B[i+2, j]])
+#
+#                 # B at Red locations
+#                 elif (((i % 2) == 0) and ((j % 2) == 0)):
+#                     B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+#                     2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+#                     -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+#                     2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+#                     -1.5 * R[i+2, j]])
+#
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#
+#     elif (bayer_pattern == "gbrg"):
+#
+#         G[::2, ::2]   = raw[::2, ::2]
+#         G[1::2, 1::2] = raw[1::2, 1::2]
+#         R[1::2, ::2]  = raw[1::2, ::2]
+#         B[::2, 1::2]  = raw[::2, 1::2]
+#
+#         # Green channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # G at Red location
+#                 if (((i % 2) != 0) and ((j % 2) == 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+#                     2. * G[i+1, j], \
+#                     -1. * R[i+2, j]])
+#                 # G at Blue location
+#                 elif (((i % 2) == 0) and ((j % 2) != 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+#                     2. * G[i+1, j],\
+#                     -1. * B[i+2, j]])
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Green: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#         # Red and Blue channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # Green locations in Red rows
+#                 if (((i % 2) != 0) and ((j % 2) != 0)):
+#                     # R at Green locations in Red rows
+#                     R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+#                      -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+#                      -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+#                      -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                       .5 * G[i+2, j]])
+#
+#                     # B at Green locations in Red rows
+#                     B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                 # Green locations in Blue rows
+#                 elif (((i % 2) == 0) and ((j % 2) == 0)):
+#
+#                     # R at Green locations in Blue rows
+#                     R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                     # B at Green locations in Blue rows
+#                     B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+#                     -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+#                     -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+#                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                     .5 * G[i+2, j]])
+#
+#                 # R at Blue locations
+#                 elif (((i % 2) == 0) and ((j % 2) != 0)):
+#                     R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+#                     2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+#                     -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+#                     2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+#                     -1.5 * B[i+2, j]])
+#
+#                 # B at Red locations
+#                 elif (((i % 2) != 0) and ((j % 2) == 0)):
+#                     B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+#                     2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+#                     -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+#                     2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+#                     -1.5 * R[i+2, j]])
+#
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#     elif (bayer_pattern == "grbg"):
+#
+#         G[::2, ::2]   = raw[::2, ::2]
+#         G[1::2, 1::2] = raw[1::2, 1::2]
+#         R[::2, 1::2]  = raw[::2, 1::2]
+#         B[1::2, ::2]  = raw[1::2, ::2]
+#
+#         # Green channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # G at Red location
+#                 if (((i % 2) == 0) and ((j % 2) != 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+#                     2. * G[i+1, j], \
+#                     -1. * R[i+2, j]])
+#                 # G at Blue location
+#                 elif (((i % 2) != 0) and ((j % 2) == 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+#                     2. * G[i+1, j],\
+#                     -1. * B[i+2, j]])
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Green: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#         # Red and Blue channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # Green locations in Red rows
+#                 if (((i % 2) == 0) and ((j % 2) == 0)):
+#                     # R at Green locations in Red rows
+#                     R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+#                      -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+#                      -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+#                      -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                       .5 * G[i+2, j]])
+#
+#                     # B at Green locations in Red rows
+#                     B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                 # Green locations in Blue rows
+#                 elif (((i % 2) != 0) and ((j % 2) != 0)):
+#
+#                     # R at Green locations in Blue rows
+#                     R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                     # B at Green locations in Blue rows
+#                     B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+#                     -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+#                     -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+#                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                     .5 * G[i+2, j]])
+#
+#                 # R at Blue locations
+#                 elif (((i % 2) != 0) and ((j % 2) == 0)):
+#                     R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+#                     2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+#                     -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+#                     2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+#                     -1.5 * B[i+2, j]])
+#
+#                 # B at Red locations
+#                 elif (((i % 2) == 0) and ((j % 2) != 0)):
+#                     B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+#                     2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+#                     -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+#                     2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+#                     -1.5 * R[i+2, j]])
+#
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#     elif (bayer_pattern == "bggr"):
+#
+#         G[::2, 1::2]  = raw[::2, 1::2]
+#         G[1::2, ::2]  = raw[1::2, ::2]
+#         R[1::2, 1::2] = raw[1::2, 1::2]
+#         B[::2, ::2]   = raw[::2, ::2]
+#
+#         # Green channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # G at Red location
+#                 if (((i % 2) != 0) and ((j % 2) != 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * R[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * R[i, j-2], 2. * G[i, j-1], 4. * R[i,j], 2. * G[i, j+1], -1. * R[i, j+2],\
+#                     2. * G[i+1, j], \
+#                     -1. * R[i+2, j]])
+#                 # G at Blue location
+#                 elif (((i % 2) == 0) and ((j % 2) == 0)):
+#                     G[i, j] = 0.125 * np.sum([-1. * B[i-2, j], \
+#                     2. * G[i-1, j], \
+#                     -1. * B[i, j-2], 2. * G[i, j-1], 4. * B[i,j], 2. * G[i, j+1], -1. * B[i, j+2], \
+#                     2. * G[i+1, j],\
+#                     -1. * B[i+2, j]])
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Green: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#         # Red and Blue channel
+#         for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+#
+#             # to display progress
+#             t0 = time.process_time()
+#
+#             for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+#
+#                 # Green locations in Red rows
+#                 if (((i % 2) != 0) and ((j % 2) == 0)):
+#                     # R at Green locations in Red rows
+#                     R[i, j] = 0.125 * np.sum([.5 * G[i-2, j],\
+#                      -1. * G[i-1, j-1], -1. * G[i-1, j+1], \
+#                      -1. * G[i, j-2], 4. * R[i, j-1], 5. * G[i,j], 4. * R[i, j+1], -1. * G[i, j+2], \
+#                      -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                       .5 * G[i+2, j]])
+#
+#                     # B at Green locations in Red rows
+#                     B[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * B[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * B[i+1,j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                 # Green locations in Blue rows
+#                 elif (((i % 2) == 0) and ((j % 2) != 0)):
+#
+#                     # R at Green locations in Blue rows
+#                     R[i, j] = 0.125 * np.sum([-1. * G[i-2, j], \
+#                     -1. * G[i-1, j-1], 4. * R[i-1, j], -1. * G[i-1, j+1], \
+#                     .5 * G[i, j-2], 5. * G[i,j], .5 * G[i, j+2], \
+#                     -1. * G[i+1, j-1], 4. * R[i+1, j],  -1. * G[i+1, j+1], \
+#                     -1. * G[i+2, j]])
+#
+#                     # B at Green locations in Blue rows
+#                     B[i, j] = 0.125 * np.sum([.5 * G[i-2, j], \
+#                     -1. * G [i-1, j-1], -1. * G[i-1, j+1], \
+#                     -1. * G[i, j-2], 4. * B[i, j-1], 5. * G[i,j], 4. * B[i, j+1], -1. * G[i, j+2], \
+#                     -1. * G[i+1, j-1], -1. * G[i+1, j+1], \
+#                     .5 * G[i+2, j]])
+#
+#                 # R at Blue locations
+#                 elif (((i % 2) == 0) and ((j % 2) == 0)):
+#                     R[i, j] = 0.125 * np.sum([-1.5 * B[i-2, j], \
+#                     2. * R[i-1, j-1], 2. * R[i-1, j+1], \
+#                     -1.5 * B[i, j-2], 6. * B[i,j], -1.5 * B[i, j+2], \
+#                     2. * R[i+1, j-1], 2. * R[i+1, j+1], \
+#                     -1.5 * B[i+2, j]])
+#
+#                 # B at Red locations
+#                 elif (((i % 2) != 0) and ((j % 2) != 0)):
+#                     B[i, j] = 0.125 * np.sum([-1.5 * R[i-2, j], \
+#                     2. * B[i-1, j-1], 2. * B[i-1, j+1], \
+#                     -1.5 * R[i, j-2], 6. * R[i,j], -1.5 * R[i, j+2], \
+#                     2. * B[i+1, j-1], 2. * B[i+1, j+1], \
+#                     -1.5 * R[i+2, j]])
+#
+#             if (timeshow):
+#                 elapsed_time = time.process_time() - t0
+#                 print("Red/Blue: row index: " + str(i-1) + " of " + str(height) + \
+#                       " | elapsed time: " + "{:.3f}".format(elapsed_time) + " seconds")
+#
+#     else:
+#         print("Invalid bayer pattern. Valid pattern can be rggb, gbrg, grbg, bggr")
+#         return demosaic_out # This will be all zeros
+#
+#     # Fill up the RGB output with interpolated values
+#     demosaic_out[0:height, 0:width, 0] = R[no_of_pixel_pad : height + no_of_pixel_pad, \
+#                                            no_of_pixel_pad : width + no_of_pixel_pad]
+#     demosaic_out[0:height, 0:width, 1] = G[no_of_pixel_pad : height + no_of_pixel_pad, \
+#                                            no_of_pixel_pad : width + no_of_pixel_pad]
+#     demosaic_out[0:height, 0:width, 2] = B[no_of_pixel_pad : height + no_of_pixel_pad, \
+#                                            no_of_pixel_pad : width + no_of_pixel_pad]
+#
+#     demosaic_out = np.clip(demosaic_out, clip_range[0], clip_range[1])
+#     return demosaic_out

IIR-Lab/ISP_pipeline/docker_guidelines.md

@@ -0,0 +1,29 @@
+# Final submission docker guidelines
+
+All final submissions should be submitted as a docker image. The docker image should be built from the dockerfile in the root of the repository. The docker image should be built with the following command:
+
+```bash
+    docker build -t <your_team_name> .
+```
+
+The docker image should be run with the following command:
+
+```bash
+    docker run -it --rm -v $(pwd)/data:/data <your_team_name> ./run.sh
+```
+
+As output, the docker image should produce images in `JPEG` format in the `/data` directory. All produced files should be named as the input files, but with the `.jpg` extension. The filenames should be the same as the RAW input filenames in `/data`. Make sure that your code does not create any other folders in the `/data` directory. Docker should contain all the necessary dependencies to run the code. It also should include the `run.sh` script as the entrypoint. Take into account that inside the docker image, the `/data` directory will be mounted to the `$(pwd)/data` directory of the host machine. This means that the docker image should be able to read the input files from the `/data` directory and write the output files to the `/data` directory.
+
+## Example
+
+We providing an example of a docker image that can be used as a reference. It can be found in our [github repository](https://github.com/createcolor/nightimaging23)
+
+Your dockerfile may look like this:
+
+```dockerfile
+FROM tensorflow/tensorflow:2.3.0
+WORKDIR /opt/app
+COPY . .
+RUN pip install -r /app/requirements.txt
+CMD ["./run.sh"]
+```

IIR-Lab/ISP_pipeline/imaging.py

@@ -0,0 +1,1293 @@
+# Note:
+#   The functions try to operate in float32 data precision
+
+# =============================================================
+# Import the libraries
+# =============================================================
+import numpy as np  # array operations
+import math         # basing math operations
+from matplotlib import pylab as plt
+import time         # measure runtime
+import utility
+import debayer
+import sys          # float precision
+from scipy import signal        # convolutions
+from scipy import interpolate   # for interpolation
+
+
+
+# =============================================================
+# class: ImageInfo
+#   Helps set up necessary information/metadata of the image
+# =============================================================
+class ImageInfo:
+    def __init__(self, name = "unknown", data = -1, is_show = False):
+        self.name   = name
+        self.data   = data
+        self.size   = np.shape(self.data)
+        self.is_show = is_show
+        self.color_space = "unknown"
+        self.bayer_pattern = "unknown"
+        self.channel_gain = (1.0, 1.0, 1.0, 1.0)
+        self.bit_depth = 0
+        self.black_level = (0, 0, 0, 0)
+        self.white_level = (1, 1, 1, 1)
+        self.color_matrix = [[1., .0, .0],\
+                             [.0, 1., .0],\
+                             [.0, .0, 1.]] # xyz2cam
+        self.min_value = np.min(self.data)
+        self.max_value = np.max(self.data)
+        self.data_type = self.data.dtype
+
+        # Display image only isShow = True
+        if (self.is_show):
+            plt.imshow(self.data)
+            plt.show()
+
+    def set_data(self, data):
+        # This function updates data and corresponding fields
+        self.data = data
+        self.size = np.shape(self.data)
+        self.data_type = self.data.dtype
+        self.min_value = np.min(self.data)
+        self.max_value = np.max(self.data)
+
+    def get_size(self):
+        return self.size
+
+    def get_width(self):
+        return self.size[1]
+
+    def get_height(self):
+        return self.size[0]
+
+    def get_depth(self):
+        if np.ndim(self.data) > 2:
+            return self.size[2]
+        else:
+            return 0
+
+    def set_color_space(self, color_space):
+        self.color_space = color_space
+
+    def get_color_space(self):
+        return self.color_space
+
+    def set_channel_gain(self, channel_gain):
+        self.channel_gain = channel_gain
+
+    def get_channel_gain(self):
+        return self.channel_gain
+
+    def set_color_matrix(self, color_matrix):
+        self.color_matrix = color_matrix
+
+    def get_color_matrix(self):
+        return self.color_matrix
+
+    def set_bayer_pattern(self, bayer_pattern):
+        self.bayer_pattern = bayer_pattern
+
+    def get_bayer_pattern(self):
+        return self.bayer_pattern
+
+    def set_bit_depth(self, bit_depth):
+        self.bit_depth = bit_depth
+
+    def get_bit_depth(self):
+        return self.bit_depth
+
+    def set_black_level(self, black_level):
+        self.black_level = black_level
+
+    def get_black_level(self):
+        return self.black_level
+
+    def set_white_level(self, white_level):
+        self.white_level = white_level
+
+    def get_white_level(self):
+        return self.white_level
+
+    def get_min_value(self):
+        return self.min_value
+
+    def get_max_value(self):
+        return self.max_value
+
+    def get_data_type(self):
+        return self.data_type
+
+    def __str__(self):
+        return "Image " + self.name + " info:" + \
+                          "\n\tname:\t" + self.name + \
+                          "\n\tsize:\t" + str(self.size) + \
+                          "\n\tcolor space:\t" + self.color_space + \
+                          "\n\tbayer pattern:\t" + self.bayer_pattern + \
+                          "\n\tchannel gains:\t" + str(self.channel_gain) + \
+                          "\n\tbit depth:\t" + str(self.bit_depth) + \
+                          "\n\tdata type:\t" + str(self.data_type) + \
+                          "\n\tblack level:\t" + str(self.black_level) + \
+                          "\n\tminimum value:\t" + str(self.min_value) + \
+                          "\n\tmaximum value:\t" + str(self.max_value)
+
+
+# =============================================================
+# function: black_level_correction
+#   subtracts the black level channel wise
+# =============================================================
+def black_level_correction(raw, black_level, white_level, clip_range):
+
+    print("----------------------------------------------------")
+    print("Running black level correction...")
+
+    # make float32 in case if it was not
+    black_level = np.float32(black_level)
+    white_level = np.float32(white_level)
+    raw = np.float32(raw)
+
+    # create new data so that original raw data do not change
+    data = np.zeros(raw.shape)
+
+    # bring data in range 0 to 1
+    data[::2, ::2]   = (raw[::2, ::2] - black_level[0]) / (white_level[0] - black_level[0])
+    data[::2, 1::2]  = (raw[::2, 1::2] - black_level[1]) / (white_level[1] - black_level[1])
+    data[1::2, ::2]  = (raw[1::2, ::2] - black_level[2]) / (white_level[2] - black_level[2])
+    data[1::2, 1::2] = (raw[1::2, 1::2]- black_level[3]) / (white_level[3] - black_level[3])
+
+    # bring within the bit depth range
+    data = data * clip_range[1]
+
+    # clip within the range
+    data = np.clip(data, clip_range[0], clip_range[1]) # upper level not necessary
+    data = np.float32(data)
+
+    return data
+
+
+# =============================================================
+# function: channel_gain_white_balance
+#   multiply with the white balance channel gains
+# =============================================================
+def channel_gain_white_balance(data, channel_gain):
+
+    print("----------------------------------------------------")
+    print("Running channel gain white balance...")
+
+    # convert into float32 in case they were not
+    data = np.float32(data)
+    channel_gain = np.float32(channel_gain)
+
+    # multiply with the channel gains
+    data[::2, ::2]   = data[::2, ::2] * channel_gain[0]
+    data[::2, 1::2]  = data[::2, 1::2] * channel_gain[1]
+    data[1::2, ::2]  = data[1::2, ::2] * channel_gain[2]
+    data[1::2, 1::2] = data[1::2, 1::2] * channel_gain[3]
+
+    # clipping within range
+    data = np.clip(data, 0., None) # upper level not necessary
+
+    return data
+
+
+# =============================================================
+# function: bad_pixel_correction
+#   correct for the bad (dead, stuck, or hot) pixels
+# =============================================================
+def bad_pixel_correction(data, neighborhood_size):
+
+    print("----------------------------------------------------")
+    print("Running bad pixel correction...")
+
+    if ((neighborhood_size % 2) == 0):
+        print("neighborhood_size shoud be odd number, recommended value 3")
+        return data
+
+    # convert to float32 in case they were not
+    # Being consistent in data format to be float32
+    data = np.float32(data)
+
+    # Separate out the quarter resolution images
+    D = {} # Empty dictionary
+    D[0] = data[::2, ::2]
+    D[1] = data[::2, 1::2]
+    D[2] = data[1::2, ::2]
+    D[3] = data[1::2, 1::2]
+
+    # number of pixels to be padded at the borders
+    no_of_pixel_pad = math.floor(neighborhood_size / 2.)
+
+    for idx in range(0, len(D)): # perform same operation for each quarter
+
+        # display progress
+        print("bad pixel correction: Quarter " + str(idx+1) + " of 4")
+
+        img = D[idx]
+        width, height = utility.helpers(img).get_width_height()
+
+        # pad pixels at the borders
+        img = np.pad(img, \
+                     (no_of_pixel_pad, no_of_pixel_pad),\
+                     'reflect') # reflect would not repeat the border value
+
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # save the middle pixel value
+                mid_pixel_val = img[i, j]
+
+                # extract the neighborhood
+                neighborhood = img[i - no_of_pixel_pad : i + no_of_pixel_pad+1,\
+                                   j - no_of_pixel_pad : j + no_of_pixel_pad+1]
+
+                # set the center pixels value same as the left pixel
+                # Does not matter replace with right or left pixel
+                # is used to replace the center pixels value
+                neighborhood[no_of_pixel_pad, no_of_pixel_pad] = neighborhood[no_of_pixel_pad, no_of_pixel_pad-1]
+
+                min_neighborhood = np.min(neighborhood)
+                max_neighborhood = np.max(neighborhood)
+
+                if (mid_pixel_val < min_neighborhood):
+                    img[i,j] = min_neighborhood
+                elif (mid_pixel_val > max_neighborhood):
+                    img[i,j] = max_neighborhood
+                else:
+                    img[i,j] = mid_pixel_val
+
+        # Put the corrected image to the dictionary
+        D[idx] = img[no_of_pixel_pad : height + no_of_pixel_pad,\
+                     no_of_pixel_pad : width + no_of_pixel_pad]
+
+    # Regrouping the data
+    data[::2, ::2]   = D[0]
+    data[::2, 1::2]  = D[1]
+    data[1::2, ::2]  = D[2]
+    data[1::2, 1::2] = D[3]
+
+    return data
+
+
+# =============================================================
+# class: demosaic
+# =============================================================
+class demosaic:
+    def __init__(self, data, bayer_pattern="rggb", clip_range=[0, 65535], name="demosaic"):
+        self.data = np.float32(data)
+        self.bayer_pattern = bayer_pattern
+        self.clip_range = clip_range
+        self.name = name
+
+    def mhc(self, timeshow=False):
+
+        print("----------------------------------------------------")
+        print("Running demosaicing using Malvar-He-Cutler algorithm...")
+
+        return debayer.debayer_mhc(self.data, self.bayer_pattern, self.clip_range, timeshow)
+
+    def post_process_local_color_ratio(self, beta):
+        # Objective is to reduce high chroma jump
+        # Beta is controlling parameter, higher gives more effect,
+        # however, too high does not make any more change
+
+        print("----------------------------------------------------")
+        print("Demosaicing post process using local color ratio...")
+
+        data = self.data
+
+        # add beta with the data to prevent divide by zero
+        data_beta = self.data + beta
+
+        # convolution kernels
+        # zeta1 averages the up, down, left, and right four values of a 3x3 window
+        zeta1 = np.multiply([[0., 1., 0.], [1., 0., 1.], [0., 1., 0.]], .25)
+        # zeta2 averages the four corner values of a 3x3 window
+        zeta2 = np.multiply([[1., 0., 1.], [0., 0., 0.], [1., 0., 1.]], .25)
+
+        # average of color ratio
+        g_over_b = signal.convolve2d(np.divide(data_beta[:, :, 1], data_beta[:, :, 2]), zeta1, mode="same", boundary="symm")
+        g_over_r = signal.convolve2d(np.divide(data_beta[:, :, 1], data_beta[:, :, 0]), zeta1, mode="same", boundary="symm")
+        b_over_g_zeta2 = signal.convolve2d(np.divide(data_beta[:, :, 2], data_beta[:, :, 1]), zeta2, mode="same", boundary="symm")
+        r_over_g_zeta2 = signal.convolve2d(np.divide(data_beta[:, :, 0], data_beta[:, :, 1]), zeta2, mode="same", boundary="symm")
+        b_over_g_zeta1 = signal.convolve2d(np.divide(data_beta[:, :, 2], data_beta[:, :, 1]), zeta1, mode="same", boundary="symm")
+        r_over_g_zeta1 = signal.convolve2d(np.divide(data_beta[:, :, 0], data_beta[:, :, 1]), zeta1, mode="same", boundary="symm")
+
+        # G at B locations and G at R locations
+        if self.bayer_pattern == "rggb":
+            # G at B locations
+            data[1::2, 1::2, 1] = -beta + np.multiply(data_beta[1::2, 1::2, 2], g_over_b[1::2, 1::2])
+            # G at R locations
+            data[::2, ::2, 1] = -beta + np.multiply(data_beta[::2, ::2, 0], g_over_r[::2, ::2])
+            # B at R locations
+            data[::2, ::2, 2] = -beta + np.multiply(data_beta[::2, ::2, 1], b_over_g_zeta2[::2, ::2])
+            # R at B locations
+            data[1::2, 1::2, 0] = -beta + np.multiply(data_beta[1::2, 1::2, 1], r_over_g_zeta2[1::2, 1::2])
+            # B at G locations
+            data[::2, 1::2, 2] = -beta + np.multiply(data_beta[::2, 1::2, 1], b_over_g_zeta1[::2, 1::2])
+            data[1::2, ::2, 2] = -beta + np.multiply(data_beta[1::2, ::2, 1], b_over_g_zeta1[1::2, ::2])
+            # R at G locations
+            data[::2, 1::2, 0] = -beta + np.multiply(data_beta[::2, 1::2, 1], r_over_g_zeta1[::2, 1::2])
+            data[1::2, ::2, 0] = -beta + np.multiply(data_beta[1::2, ::2, 1], r_over_g_zeta1[1::2, ::2])
+
+        elif self.bayer_pattern == "grbg":
+            # G at B locations
+            data[1::2, ::2, 1] = -beta + np.multiply(data_beta[1::2, ::2, 2], g_over_b[1::2, 1::2])
+            # G at R locations
+            data[::2, 1::2, 1] = -beta + np.multiply(data_beta[::2, 1::2, 0], g_over_r[::2, 1::2])
+            # B at R locations
+            data[::2, 1::2, 2] = -beta + np.multiply(data_beta[::2, 1::2, 1], b_over_g_zeta2[::2, 1::2])
+            # R at B locations
+            data[1::2, ::2, 0] = -beta + np.multiply(data_beta[1::2, ::2, 1], r_over_g_zeta2[1::2, ::2])
+            # B at G locations
+            data[::2, ::2, 2] = -beta + np.multiply(data_beta[::2, ::2, 1], b_over_g_zeta1[::2, ::2])
+            data[1::2, 1::2, 2] = -beta + np.multiply(data_beta[1::2, 1::2, 1], b_over_g_zeta1[1::2, 1::2])
+            # R at G locations
+            data[::2, ::2, 0] = -beta + np.multiply(data_beta[::2, ::2, 1], r_over_g_zeta1[::2, ::2])
+            data[1::2, 1::2, 0] = -beta + np.multiply(data_beta[1::2, 1::2, 1], r_over_g_zeta1[1::2, 1::2])
+
+        elif self.bayer_pattern == "gbrg":
+            # G at B locations
+            data[::2, 1::2, 1] = -beta + np.multiply(data_beta[::2, 1::2, 2], g_over_b[::2, 1::2])
+            # G at R locations
+            data[1::2, ::2, 1] = -beta + np.multiply(data_beta[1::2, ::2, 0], g_over_r[1::2, ::2])
+            # B at R locations
+            data[1::2, ::2, 2] = -beta + np.multiply(data_beta[1::2, ::2, 1], b_over_g_zeta2[1::2, ::2])
+            # R at B locations
+            data[::2, 1::2, 0] = -beta + np.multiply(data_beta[::2, 1::2, 1], r_over_g_zeta2[::2, 1::2])
+            # B at G locations
+            data[::2, ::2, 2] = -beta + np.multiply(data_beta[::2, ::2, 1], b_over_g_zeta1[::2, ::2])
+            data[1::2, 1::2, 2] = -beta + np.multiply(data_beta[1::2, 1::2, 1], b_over_g_zeta1[1::2, 1::2])
+            # R at G locations
+            data[::2, ::2, 0] = -beta + np.multiply(data_beta[::2, ::2, 1], r_over_g_zeta1[::2, ::2])
+            data[1::2, 1::2, 0] = -beta + np.multiply(data_beta[1::2, 1::2, 1], r_over_g_zeta1[1::2, 1::2])
+
+        elif self.bayer_pattern == "bggr":
+            # G at B locations
+            data[::2, ::2, 1] = -beta + np.multiply(data_beta[::2, ::2, 2], g_over_b[::2, ::2])
+            # G at R locations
+            data[1::2, 1::2, 1] = -beta + np.multiply(data_beta[1::2, 1::2, 0], g_over_r[1::2, 1::2])
+            # B at R locations
+            data[1::2, 1::2, 2] = -beta + np.multiply(data_beta[1::2, 1::2, 1], b_over_g_zeta2[1::2, 1::2])
+            # R at B locations
+            data[::2, ::2, 0] = -beta + np.multiply(data_beta[::2, ::2, 1], r_over_g_zeta2[::2, ::2])
+            # B at G locations
+            data[::2, 1::2, 2] = -beta + np.multiply(data_beta[::2, 1::2, 1], b_over_g_zeta1[::2, 1::2])
+            data[1::2, ::2, 2] = -beta + np.multiply(data_beta[1::2, ::2, 1], b_over_g_zeta1[1::2, ::2])
+            # R at G locations
+            data[::2, 1::2, 0] = -beta + np.multiply(data_beta[::2, 1::2, 1], r_over_g_zeta1[::2, 1::2])
+            data[1::2, ::2, 0] = -beta + np.multiply(data_beta[1::2, ::2, 1], r_over_g_zeta1[1::2, ::2])
+
+
+        return np.clip(data, self.clip_range[0], self.clip_range[1])
+
+
+    def directionally_weighted_gradient_based_interpolation(self):
+        # Reference:
+        # http://www.arl.army.mil/arlreports/2010/ARL-TR-5061.pdf
+
+        print("----------------------------------------------------")
+        print("Running demosaicing using directionally weighted gradient based interpolation...")
+
+        # Fill up the green channel
+        G = debayer.fill_channel_directional_weight(self.data, self.bayer_pattern)
+
+        B, R = debayer.fill_br_locations(self.data, G, self.bayer_pattern)
+
+        width, height = utility.helpers(self.data).get_width_height()
+        output = np.empty((height, width, 3), dtype=np.float32)
+        output[:, :, 0] = R
+        output[:, :, 1] = G
+        output[:, :, 2] = B
+
+        return np.clip(output, self.clip_range[0], self.clip_range[1])
+
+
+    def post_process_median_filter(self, edge_detect_kernel_size=3, edge_threshold=0, median_filter_kernel_size=3, clip_range=[0, 65535]):
+        # Objective is to reduce the zipper effect around the edges
+        # Inputs:
+        #   edge_detect_kernel_size: the neighborhood size used to detect edges
+        #   edge_threshold: the threshold value above which (compared against)
+        #                   the gradient_magnitude to declare if it is an edge
+        #   median_filter_kernel_size: the neighborhood size used to perform
+        #                               median filter operation
+        #   clip_range: used for scaling in edge_detection
+        #
+        # Output:
+        #   output: median filtered output around the edges
+        #   edge_location: a debug image to see where the edges were detected
+        #                   based on the threshold
+
+
+        # detect edge locations
+        edge_location = utility.edge_detection(self.data).sobel(edge_detect_kernel_size, "is_edge", edge_threshold, clip_range)
+
+        # allocate space for output
+        output = np.empty(np.shape(self.data), dtype=np.float32)
+
+        if (np.ndim(self.data) > 2):
+
+            for i in range(0, np.shape(self.data)[2]):
+                output[:, :, i] = utility.helpers(self.data[:, :, i]).edge_wise_median(median_filter_kernel_size, edge_location[:, :, i])
+
+        elif (np.ndim(self.data) == 2):
+            output = utility.helpers(self.data).edge_wise_median(median_filter_kernel_size, edge_location)
+
+        return output, edge_location
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: lens_shading_correction
+#   Correct the lens shading / vignetting
+# =============================================================
+class lens_shading_correction:
+    def __init__(self, data, name="lens_shading_correction"):
+        # convert to float32 in case it was not
+        self.data = np.float32(data)
+        self.name = name
+
+    def flat_field_compensation(self, dark_current_image, flat_field_image):
+        # dark_current_image:
+        #       is captured from the camera with cap on
+        #       and fully dark condition, several images captured and
+        #       temporally averaged
+        # flat_field_image:
+        #       is found by capturing an image of a flat field test chart
+        #       with certain lighting condition
+        # Note: flat_field_compensation is memory intensive procedure because
+        #       both the dark_current_image and flat_field_image need to be
+        #       saved in memory beforehand
+        print("----------------------------------------------------")
+        print("Running lens shading correction with flat field compensation...")
+
+        # convert to float32 in case it was not
+        dark_current_image = np.float32(dark_current_image)
+        flat_field_image = np.float32(flat_field_image)
+        temp = flat_field_image - dark_current_image
+        return np.average(temp) * np.divide((self.data - dark_current_image), temp)
+
+    def approximate_mathematical_compensation(self, params, clip_min=0, clip_max=65535):
+        # parms:
+        #       parameters of a parabolic model y = a*(x-b)^2 + c
+        #       For example, params = [0.01759, -28.37, -13.36]
+        # Note: approximate_mathematical_compensation require less memory
+        print("----------------------------------------------------")
+        print("Running lens shading correction with approximate mathematical compensation...")
+        width, height = utility.helpers(self.data).get_width_height()
+
+        center_pixel_pos = [height/2, width/2]
+        max_distance = utility.distance_euclid(center_pixel_pos, [height, width])
+
+        # allocate memory for output
+        temp = np.empty((height, width), dtype=np.float32)
+
+        for i in range(0, height):
+            for j in range(0, width):
+                distance = utility.distance_euclid(center_pixel_pos, [i, j]) / max_distance
+                # parabolic model
+                gain = params[0] * (distance - params[1])**2 + params[2]
+                temp[i, j] = self.data[i, j] * gain
+
+        temp = np.clip(temp, clip_min, clip_max)
+        return temp
+
+    def __str__(self):
+        return "lens shading correction. There are two methods: " + \
+                "\n (1) flat_field_compensation: requires dark_current_image and flat_field_image" + \
+                "\n (2) approximate_mathematical_compensation:"
+
+
+# =============================================================
+# class: lens_shading_correction
+#   Correct the lens shading / vignetting
+# =============================================================
+class bayer_denoising:
+    def __init__(self, data, name="bayer_denoising"):
+        # convert to float32 in case it was not
+        self.data = np.float32(data)
+        self.name = name
+
+    def utilize_hvs_behavior(self, bayer_pattern, initial_noise_level, hvs_min, hvs_max, threshold_red_blue, clip_range):
+        # Objective: bayer denoising
+        # Inputs:
+        #   bayer_pattern:  rggb, gbrg, grbg, bggr
+        #   initial_noise_level:
+        # Output:
+        #   denoised bayer raw output
+        # Source: Based on paper titled "Noise Reduction for CFA Image Sensors
+        #   Exploiting HVS Behaviour," by Angelo Bosco, Sebastiano Battiato,
+        #   Arcangelo Bruna and Rosetta Rizzo
+        #   Sensors 2009, 9, 1692-1713; doi:10.3390/s90301692
+
+        print("----------------------------------------------------")
+        print("Running bayer denoising utilizing hvs behavior...")
+
+        # copy the self.data to raw and we will only work on raw
+        # to make sure no change happen to self.data
+        raw = self.data
+        raw = np.clip(raw, clip_range[0], clip_range[1])
+        width, height = utility.helpers(raw).get_width_height()
+
+        # First make the bayer_pattern rggb
+        # The algorithm is written only for rggb pattern, thus convert all other
+        # pattern to rggb. Furthermore, this shuffling does not affect the
+        # algorithm output
+        if (bayer_pattern != "rggb"):
+            raw = utility.helpers(self.data).shuffle_bayer_pattern(bayer_pattern, "rggb")
+
+        # fixed neighborhood_size
+        neighborhood_size = 5 # we are keeping this fixed
+                              # bigger size such as 9 can be declared
+                              # however, the code need to be changed then
+
+        # pad two pixels at the border
+        no_of_pixel_pad = math.floor(neighborhood_size / 2)   # number of pixels to pad
+
+        raw = np.pad(raw, \
+                     (no_of_pixel_pad, no_of_pixel_pad),\
+                     'reflect') # reflect would not repeat the border value
+
+        # allocating space for denoised output
+        denoised_out = np.empty((height, width), dtype=np.float32)
+
+        texture_degree_debug = np.empty((height, width), dtype=np.float32)
+        for i in range(no_of_pixel_pad, height + no_of_pixel_pad):
+            for j in range(no_of_pixel_pad, width + no_of_pixel_pad):
+
+                # center pixel
+                center_pixel = raw[i, j]
+
+                # signal analyzer block
+                half_max = clip_range[1] / 2
+                if (center_pixel <= half_max):
+                    hvs_weight = -(((hvs_max - hvs_min) * center_pixel) / half_max) + hvs_max
+                else:
+                    hvs_weight = (((center_pixel - clip_range[1]) * (hvs_max - hvs_min))/(clip_range[1] - half_max)) + hvs_max
+
+                # noise level estimator previous value
+                if (j < no_of_pixel_pad+2):
+                    noise_level_previous_red   = initial_noise_level
+                    noise_level_previous_blue  = initial_noise_level
+                    noise_level_previous_green = initial_noise_level
+                else:
+                    noise_level_previous_green = noise_level_current_green
+                    if ((i % 2) == 0): # red
+                        noise_level_previous_red = noise_level_current_red
+                    elif ((i % 2) != 0): # blue
+                        noise_level_previous_blue = noise_level_current_blue
+
+                # Processings depending on Green or Red/Blue
+                # Red
+                if (((i % 2) == 0) and ((j % 2) == 0)):
+                    # get neighborhood
+                    neighborhood = [raw[i-2, j-2], raw[i-2, j], raw[i-2, j+2],\
+                                    raw[i, j-2], raw[i, j+2],\
+                                    raw[i+2, j-2], raw[i+2, j], raw[i+2, j+2]]
+
+                    # absolute difference from the center pixel
+                    d =  np.abs(neighborhood - center_pixel)
+
+                    # maximum and minimum difference
+                    d_max = np.max(d)
+                    d_min = np.min(d)
+
+                    # calculate texture_threshold
+                    texture_threshold = hvs_weight + noise_level_previous_red
+
+                    # texture degree analyzer
+                    if (d_max <= threshold_red_blue):
+                        texture_degree = 1.
+                    elif ((d_max > threshold_red_blue) and (d_max <= texture_threshold)):
+                        texture_degree = -((d_max - threshold_red_blue) / (texture_threshold - threshold_red_blue)) + 1.
+                    elif (d_max > texture_threshold):
+                        texture_degree = 0.
+
+                    # noise level estimator update
+                    noise_level_current_red = texture_degree * d_max + (1 - texture_degree) * noise_level_previous_red
+
+                # Blue
+                elif (((i % 2) != 0) and ((j % 2) != 0)):
+
+                    # get neighborhood
+                    neighborhood = [raw[i-2, j-2], raw[i-2, j], raw[i-2, j+2],\
+                                    raw[i, j-2], raw[i, j+2],\
+                                    raw[i+2, j-2], raw[i+2, j], raw[i+2, j+2]]
+
+                    # absolute difference from the center pixel
+                    d =  np.abs(neighborhood - center_pixel)
+
+                    # maximum and minimum difference
+                    d_max = np.max(d)
+                    d_min = np.min(d)
+
+                    # calculate texture_threshold
+                    texture_threshold = hvs_weight + noise_level_previous_blue
+
+                    # texture degree analyzer
+                    if (d_max <= threshold_red_blue):
+                        texture_degree = 1.
+                    elif ((d_max > threshold_red_blue) and (d_max <= texture_threshold)):
+                        texture_degree = -((d_max - threshold_red_blue) / (texture_threshold - threshold_red_blue)) + 1.
+                    elif (d_max > texture_threshold):
+                        texture_degree = 0.
+
+                    # noise level estimator update
+                    noise_level_current_blue = texture_degree * d_max + (1 - texture_degree) * noise_level_previous_blue
+
+                # Green
+                elif ((((i % 2) == 0) and ((j % 2) != 0)) or (((i % 2) != 0) and ((j % 2) == 0))):
+
+                    neighborhood = [raw[i-2, j-2], raw[i-2, j], raw[i-2, j+2],\
+                                    raw[i-1, j-1], raw[i-1, j+1],\
+                                    raw[i, j-2], raw[i, j+2],\
+                                    raw[i+1, j-1], raw[i+1, j+1],\
+                                    raw[i+2, j-2], raw[i+2, j], raw[i+2, j+2]]
+
+                    # difference from the center pixel
+                    d = np.abs(neighborhood - center_pixel)
+
+                    # maximum and minimum difference
+                    d_max = np.max(d)
+                    d_min = np.min(d)
+
+                    # calculate texture_threshold
+                    texture_threshold = hvs_weight + noise_level_previous_green
+
+                    # texture degree analyzer
+                    if (d_max == 0):
+                        texture_degree = 1
+                    elif ((d_max > 0) and (d_max <= texture_threshold)):
+                        texture_degree = -(d_max / texture_threshold) + 1.
+                    elif (d_max > texture_threshold):
+                        texture_degree = 0
+
+                    # noise level estimator update
+                    noise_level_current_green = texture_degree * d_max + (1 - texture_degree) * noise_level_previous_green
+
+                # similarity threshold calculation
+                if (texture_degree == 1):
+                    threshold_low = threshold_high = d_max
+                elif (texture_degree == 0):
+                    threshold_low = d_min
+                    threshold_high = (d_max + d_min) / 2
+                elif ((texture_degree > 0) and (texture_degree < 1)):
+                    threshold_high = (d_max + ((d_max + d_min) / 2)) / 2
+                    threshold_low = (d_min + threshold_high) / 2
+
+                # weight computation
+                weight = np.empty(np.size(d), dtype=np.float32)
+                pf = 0.
+                for w_i in range(0, np.size(d)):
+                    if (d[w_i] <= threshold_low):
+                        weight[w_i] = 1.
+                    elif (d[w_i] > threshold_high):
+                        weight[w_i] = 0.
+                    elif ((d[w_i] > threshold_low) and (d[w_i] < threshold_high)):
+                        weight[w_i] = 1. + ((d[w_i] - threshold_low) / (threshold_low - threshold_high))
+
+                    pf += weight[w_i] * neighborhood[w_i] + (1. - weight[w_i]) * center_pixel
+
+                denoised_out[i - no_of_pixel_pad, j-no_of_pixel_pad] = pf / np.size(d)
+                # texture_degree_debug is a debug output
+                texture_degree_debug[i - no_of_pixel_pad, j-no_of_pixel_pad] = texture_degree
+
+        if (bayer_pattern != "rggb"):
+            denoised_out = utility.shuffle_bayer_pattern(denoised_out, "rggb", bayer_pattern)
+
+        return np.clip(denoised_out, clip_range[0], clip_range[1]), texture_degree_debug
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: color_correction
+#   Correct the color in linaer domain
+# =============================================================
+class color_correction:
+    def __init__(self, data, color_matrix, color_space="srgb", illuminant="d65", name="color correction", clip_range=[0, 65535]):
+        # Inputs:
+        #   data:   linear rgb image before nonlinearity/gamma
+        #   xyz2cam: 3x3 matrix found from the camera metedata, specifically
+        #            color matrix 2 from the metadata
+        #   color_space: output color space
+        #   illuminance: the illuminant of the lighting condition
+        #   name: name of the class
+        self.data = np.float32(data)
+        self.xyz2cam = np.float32(color_matrix)
+        self.color_space = color_space
+        self.illuminant = illuminant
+        self.name = name
+        self.clip_range = clip_range
+
+    def get_rgb2xyz(self):
+        # Objective: get the rgb2xyz matrix dependin on the output color space
+        #            and the illuminant
+        # Source: http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
+        if (self.color_space == "srgb"):
+            if (self.illuminant == "d65"):
+                return [[.4124564,  .3575761,  .1804375],\
+                        [.2126729,  .7151522,  .0721750],\
+                        [.0193339,  .1191920,  .9503041]]
+            elif (self.illuminant == "d50"):
+                return [[.4360747,  .3850649,  .1430804],\
+                        [.2225045,  .7168786,  .0606169],\
+                        [.0139322,  .0971045,  .7141733]]
+            else:
+                print("for now, color_space must be d65 or d50")
+                return
+
+        elif (self.color_space == "adobe-rgb-1998"):
+            if (self.illuminant == "d65"):
+                return [[.5767309,  .1855540,  .1881852],\
+                        [.2973769,  .6273491,  .0752741],\
+                        [.0270343,  .0706872,  .9911085]]
+            elif (self.illuminant == "d50"):
+                return [[.6097559,  .2052401,  .1492240],\
+                        [.3111242,  .6256560,  .0632197],\
+                        [.0194811,  .0608902,  .7448387]]
+            else:
+                print("for now, illuminant must be d65 or d50")
+                return
+        else:
+            print("for now, color_space must be srgb or adobe-rgb-1998")
+            return
+
+    def calculate_cam2rgb(self):
+        # Objective: Calculates the color correction matrix
+
+        # matric multiplication
+        rgb2cam = np.dot(self.xyz2cam, self.get_rgb2xyz())
+
+        # make sum of each row to be 1.0, necessary to preserve white balance
+        # basically divice each value by its row wise sum
+        rgb2cam = np.divide(rgb2cam, np.reshape(np.sum(rgb2cam, 1), [3, 1]))
+
+        # - inverse the matrix to get cam2rgb.
+        # - cam2rgb should also have the characteristic that sum of each row
+        # equal to 1.0 to preserve white balance
+        # - check if rgb2cam is invertible by checking the condition of
+        # rgb2cam. If rgb2cam is singular it will give a warning and
+        # return an identiry matrix
+        if (np.linalg.cond(rgb2cam) < (1 / sys.float_info.epsilon)):
+            return np.linalg.inv(rgb2cam) # this is cam2rgb / color correction matrix
+        else:
+            print("Warning! matrix not invertible.")
+            return np.identity(3, dtype=np.float32)
+
+    def apply_cmatrix(self):
+        # Objective: Apply the color correction matrix (cam2rgb)
+
+        print("----------------------------------------------------")
+        print("running color correction...")
+
+        # check if data is 3 dimensional
+        if (np.ndim(self.data) != 3):
+            print("data need to be three dimensional")
+            return
+
+        # get the color correction matrix
+        cam2rgb = self.calculate_cam2rgb()
+
+        # get width and height
+        width, height = utility.helpers(self.data).get_width_height()
+
+        # apply the matrix
+        R = self.data[:, :, 0]
+        G = self.data[:, :, 1]
+        B = self.data[:, :, 2]
+
+        color_corrected = np.empty((height, width, 3), dtype=np.float32)
+        color_corrected[:, :, 0] = R * cam2rgb[0, 0] + G * cam2rgb[0, 1] + B * cam2rgb[0, 2]
+        color_corrected[:, :, 1] = R * cam2rgb[1, 0] + G * cam2rgb[1, 1] + B * cam2rgb[1, 2]
+        color_corrected[:, :, 2] = R * cam2rgb[2, 0] + G * cam2rgb[2, 1] + B * cam2rgb[2, 2]
+
+        return np.clip(color_corrected, self.clip_range[0], self.clip_range[1])
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: nonlinearity
+#   apply gamma or degamma
+# =============================================================
+class nonlinearity:
+    def __init__(self, data, name="nonlinearity"):
+        self.data = np.float32(data)
+        self.name = name
+
+    def luma_adjustment(self, multiplier, clip_range=[0, 65535]):
+        # The multiplier is applied only on luma channel
+        # by a multipler in log10 scale:
+        #   multipler 10 means multiplied by 1.
+        #   multipler 100 means multiplied by 2. as such
+
+        print("----------------------------------------------------")
+        print("Running brightening...")
+
+        return np.clip(np.log10(multiplier) * self.data, clip_range[0], clip_range[1])
+
+    def by_value(self, value, clip_range):
+
+        print("----------------------------------------------------")
+        print("Running nonlinearity by value...")
+
+        # clip within the range
+        data = np.clip(self.data, clip_range[0], clip_range[1])
+        # make 0 to 1
+        data = data / clip_range[1]
+        # apply nonlinearity
+        return np.clip(clip_range[1] * (data**value), clip_range[0], clip_range[1])
+
+    def by_table(self, table, nonlinearity_type="gamma", clip_range=[0, 65535]):
+
+        print("----------------------------------------------------")
+        print("Running nonlinearity by table...")
+
+        gamma_table = np.loadtxt(table)
+        gamma_table = clip_range[1] * gamma_table / np.max(gamma_table)
+        linear_table = np.linspace(clip_range[0], clip_range[1], np.size(gamma_table))
+
+        # linear interpolation, query is the self.data
+        if (nonlinearity_type == "gamma"):
+            # mapping is from linear_table to gamma_table
+            return np.clip(np.interp(self.data, linear_table, gamma_table), clip_range[0], clip_range[1])
+        elif (nonlinearity_type == "degamma"):
+            # mapping is from gamma_table to linear_table
+            return np.clip(np.interp(self.data, gamma_table, linear_table), clip_range[0], clip_range[1])
+
+    def by_equation(self, a, b, clip_range):
+
+        print("----------------------------------------------------")
+        print("Running nonlinearity by equation...")
+
+        # clip within the range
+        data = np.clip(self.data, clip_range[0], clip_range[1])
+        # make 0 to 1
+        data = data / clip_range[1]
+
+        # apply nonlinearity
+        return np.clip(clip_range[1] * (a * np.exp(b * data) + data + a * data - a * np.exp(b) * data - a), clip_range[0], clip_range[1])
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: tone_mapping
+#   improve the overall tone of the image
+# =============================================================
+class tone_mapping:
+    def __init__(self, data, name="tone mapping"):
+        self.data = np.float32(data)
+        self.name = name
+
+    def nonlinear_masking(self, strength_multiplier=1.0, gaussian_kernel_size=[5, 5], gaussian_sigma=1.0, clip_range=[0, 65535]):
+        # Objective: improves the overall tone of the image
+        # Inputs:
+        #   strength_multiplier: >0. The higher the more aggressing tone mapping
+        #   gaussian_kernel_size: kernel size for calculating the mask image
+        #   gaussian_sigma: spread of the gaussian kernel for calculating the
+        #                   mask image
+        #
+        # Source:
+        # N. Moroney, “Local color correction using non-linear masking”,
+        # Proc. IS&T/SID 8th Color Imaging Conference, pp. 108-111, (2000)
+        #
+        # Note, Slight changes is carried by mushfiqul alam, specifically
+        # introducing the strength_multiplier
+
+        print("----------------------------------------------------")
+        print("Running tone mapping by non linear masking...")
+
+        # convert to gray image
+        if (np.ndim(self.data) == 3):
+            gray_image = utility.color_conversion(self.data).rgb2gray()
+        else:
+            gray_image = self.data
+
+        # gaussian blur the gray image
+        gaussian_kernel = utility.create_filter().gaussian(gaussian_kernel_size, gaussian_sigma)
+
+        # the mask image:   (1) blur
+        #                   (2) bring within range 0 to 1
+        #                   (3) multiply with strength_multiplier
+        mask = signal.convolve2d(gray_image, gaussian_kernel, mode="same", boundary="symm")
+        mask = strength_multiplier * mask / clip_range[1]
+
+        # calculate the alpha image
+        temp = np.power(0.5, mask)
+        if (np.ndim(self.data) == 3):
+            width, height = utility.helpers(self.data).get_width_height()
+            alpha = np.empty((height, width, 3), dtype=np.float32)
+            alpha[:, :, 0] = temp
+            alpha[:, :, 1] = temp
+            alpha[:, :, 2] = temp
+        else:
+            alpha = temp
+
+        # output
+        return np.clip(clip_range[1] * np.power(self.data/clip_range[1], alpha), clip_range[0], clip_range[1])
+
+    def dynamic_range_compression(self, drc_type="normal", drc_bound=[-40., 260.], clip_range=[0, 65535]):
+
+        ycc = utility.color_conversion(self.data).rgb2ycc("bt601")
+        y = ycc[:, :, 0]
+        cb = ycc[:, :, 1]
+        cr = ycc[:, :, 2]
+
+        if (drc_type == "normal"):
+            edge = y
+        elif (drc_type == "joint"):
+            edge = utility.edge_detection(y).sobel(3, "gradient_magnitude")
+
+        y_bilateral_filtered = utility.special_function(y).bilateral_filter(edge)
+        detail = np.divide(ycc[:, :, 0], y_bilateral_filtered)
+
+        C = drc_bound[0] * clip_range[1] / 255.
+        temp = drc_bound[1] * clip_range[1] / 255.
+        F = (temp * (C + clip_range[1])) / (clip_range[1] * (temp - C))
+        y_bilateral_filtered_contrast_reduced = F * (y_bilateral_filtered - (clip_range[1] / 2.)) + (clip_range[1] / 2.)
+
+        y_out = np.multiply(y_bilateral_filtered_contrast_reduced, detail)
+
+        ycc_out = ycc
+        ycc_out[:, :, 0] = y_out
+        rgb_out = utility.color_conversion(ycc_out).ycc2rgb("bt601")
+
+        return np.clip(rgb_out, clip_range[0], clip_range[1])
+
+
+# =============================================================
+# class: sharpening
+#   sharpens the image
+# =============================================================
+class sharpening:
+    def __init__(self, data, name="sharpening"):
+        self.data = np.float32(data)
+        self.name = name
+
+    def unsharp_masking(self, gaussian_kernel_size=[5, 5], gaussian_sigma=2.0,\
+                        slope=1.5, tau_threshold=0.05, gamma_speed=4., clip_range=[0, 65535]):
+        # Objective: sharpen image
+        # Input:
+        #   gaussian_kernel_size:   dimension of the gaussian blur filter kernel
+        #
+        #   gaussian_sigma:         spread of the gaussian blur filter kernel
+        #                           bigger sigma more sharpening
+        #
+        #   slope:                  controls the boost.
+        #                           the amount of sharpening, higher slope
+        #                           means more aggresssive sharpening
+        #
+        #   tau_threshold:          controls the amount of coring.
+        #                           threshold value till which the image is
+        #                           not sharpened. The lower the value of
+        #                           tau_threshold the more frequencies
+        #                           goes through the sharpening process
+        #
+        #   gamma_speed:            controls the speed of convergence to the slope
+        #                           smaller value gives a little bit more
+        #                           sharpened image, this may be a fine tuner
+
+        print("----------------------------------------------------")
+        print("Running sharpening by unsharp masking...")
+
+        # create gaussian kernel
+        gaussian_kernel = utility.create_filter().gaussian(gaussian_kernel_size, gaussian_sigma)
+
+        # convolove the image with the gaussian kernel
+        # first input is the image
+        # second input is the kernel
+        # output shape will be the same as the first input
+        # boundary will be padded by using symmetrical method while convolving
+        if np.ndim(self.data > 2):
+            image_blur = np.empty(np.shape(self.data), dtype=np.float32)
+            for i in range(0, np.shape(self.data)[2]):
+                image_blur[:, :, i] = signal.convolve2d(self.data[:, :, i], gaussian_kernel, mode="same", boundary="symm")
+        else:
+            image_blur = signal.convolove2d(self.data, gaussian_kernel, mode="same", boundary="symm")
+
+        # the high frequency component image
+        image_high_pass = self.data - image_blur
+
+        # soft coring (see in utility)
+        # basically pass the high pass image via a slightly nonlinear function
+        tau_threshold = tau_threshold * clip_range[1]
+
+        # add the soft cored high pass image to the original and clip
+        # within range and return
+        return np.clip(self.data + utility.special_function(\
+                   image_high_pass).soft_coring(\
+                   slope, tau_threshold, gamma_speed), clip_range[0], clip_range[1])
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: noise_reduction
+#   reduce noise of the nonlinear image (after gamma)
+# =============================================================
+class noise_reduction:
+    def __init__(self, data, clip_range=[0, 65535], name="noise reduction"):
+        self.data = np.float32(data)
+        self.clip_range = clip_range
+        self.name = name
+
+    def sigma_filter(self, neighborhood_size=7, sigma=[6, 6, 6]):
+
+        print("----------------------------------------------------")
+        print("Running noise reduction by sigma filter...")
+
+        if np.ndim(self.data > 2): # if rgb image
+            output = np.empty(np.shape(self.data), dtype=np.float32)
+            for i in range(0, np.shape(self.data)[2]):
+                output[:, :, i] = utility.helpers(self.data[:, :, i]).sigma_filter_helper(neighborhood_size, sigma[i])
+            return np.clip(output, self.clip_range[0], self.clip_range[1])
+        else: # gray image
+            return np.clip(utility.helpers(self.data).sigma_filter_helper(neighborhood_size, sigma), self.clip_range[0], self.clip_range[1])
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: distortion_correction
+#   correct the distortion
+# =============================================================
+class distortion_correction:
+    def __init__(self, data, name="distortion correction"):
+        self.data = np.float32(data)
+        self.name = name
+
+
+    def empirical_correction(self, correction_type="pincushion-1", strength=0.1, zoom_type="crop", clip_range=[0, 65535]):
+        #------------------------------------------------------
+        # Objective:
+        #   correct geometric distortion with the assumption that the distortion
+        #   is symmetric and the center is at the center of of the image
+        # Input:
+        #   correction_type:    which type of correction needed to be carried
+        #                       out, choose one the four:
+        #                       pincushion-1, pincushion-2, barrel-1, barrel-2
+        #                       1 and 2 are difference between the power
+        #                       over the radius
+        #
+        #   strength:           should be equal or greater than 0.
+        #                       0 means no correction will be done.
+        #                       if negative value were applied correction_type
+        #                       will be reversed. Thus,>=0 value expected.
+        #
+        #   zoom_type:          either "fit" or "crop"
+        #                       fit will return image with full content
+        #                       in the whole area
+        #                       crop will return image will 0 values outsise
+        #                       the border
+        #
+        #   clip_range:         to clip the final image within the range
+        #------------------------------------------------------
+
+        if (strength < 0):
+            print("Warning! strength should be equal of greater than 0.")
+            return self.data
+
+        print("----------------------------------------------------")
+        print("Running distortion correction by empirical method...")
+
+        # get half_width and half_height, assume this is the center
+        width, height = utility.helpers(self.data).get_width_height()
+        half_width = width / 2
+        half_height = height / 2
+
+        # create a meshgrid of points
+        xi, yi = np.meshgrid(np.linspace(-half_width, half_width, width),\
+                             np.linspace(-half_height, half_height, height))
+
+        # cartesian to polar coordinate
+        r = np.sqrt(xi**2 + yi**2)
+        theta = np.arctan2(yi, xi)
+
+        # maximum radius
+        R = math.sqrt(width**2 + height**2)
+
+        # make r within range 0~1
+        r = r / R
+
+        # apply the radius to the desired transformation
+        s = utility.special_function(r).distortion_function(correction_type, strength)
+
+        # select a scaling_parameter based on zoon_type and k value
+        if ((correction_type=="barrel-1") or (correction_type=="barrel-2")):
+            if (zoom_type == "fit"):
+                scaling_parameter = r[0, 0] / s[0, 0]
+            elif (zoom_type == "crop"):
+                scaling_parameter = 1. / (1. + strength * (np.min([half_width, half_height])/R)**2)
+        elif ((correction_type=="pincushion-1") or (correction_type=="pincushion-2")):
+            if (zoom_type == "fit"):
+                scaling_parameter = 1. / (1. + strength * (np.min([half_width, half_height])/R)**2)
+            elif (zoom_type == "crop"):
+                scaling_parameter = r[0, 0] / s[0, 0]
+
+        # multiply by scaling_parameter and un-normalize
+        s = s * scaling_parameter * R
+
+        # convert back to cartesian coordinate and add back the center coordinate
+        xt = np.multiply(s, np.cos(theta))
+        yt = np.multiply(s, np.sin(theta))
+
+        # interpolation
+        if np.ndim(self.data == 3):
+
+            output = np.empty(np.shape(self.data), dtype=np.float32)
+
+            output[:, :, 0] = utility.helpers(self.data[:, :, 0]).bilinear_interpolation(xt + half_width, yt + half_height)
+            output[:, :, 1] = utility.helpers(self.data[:, :, 1]).bilinear_interpolation(xt + half_width, yt + half_height)
+            output[:, :, 2] = utility.helpers(self.data[:, :, 2]).bilinear_interpolation(xt + half_width, yt + half_height)
+
+        elif np.ndim(self.data == 2):
+
+            output = utility.helpers(self.data).bilinear_interpolation(xt + half_width, yt + half_height)
+
+        return np.clip(output, clip_range[0], clip_range[1])
+
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: memory_color_enhancement
+#   enhance memory colors such as sky, grass, skin color
+# =============================================================
+class memory_color_enhancement:
+    def __init__(self, data, name="memory color enhancement"):
+        self.data = np.float32(data)
+        self.name = name
+
+    def by_hue_squeeze(self, target_hue, hue_preference, hue_sigma, is_both_side, multiplier, chroma_preference, chroma_sigma, color_space="srgb", illuminant="d65", clip_range=[0, 65535], cie_version="1931"):
+
+        # RGB to xyz
+        data = utility.color_conversion(self.data).rgb2xyz(color_space, clip_range)
+        # xyz to lab
+        data = utility.color_conversion(data).xyz2lab(cie_version, illuminant)
+        # lab to lch
+        data = utility.color_conversion(data).lab2lch()
+
+        # hue squeezing
+        # we are traversing through different color preferences
+        width, height = utility.helpers(self.data).get_width_height()
+        hue_correction = np.zeros((height, width), dtype=np.float32)
+        for i in range(0, np.size(target_hue)):
+
+            delta_hue = data[:, :, 2] - hue_preference[i]
+
+            if is_both_side[i]:
+                weight_temp = np.exp( -np.power(data[:, :, 2] - target_hue[i], 2) / (2 * hue_sigma[i]**2)) + \
+                              np.exp( -np.power(data[:, :, 2] + target_hue[i], 2) / (2 * hue_sigma[i]**2))
+            else:
+                weight_temp = np.exp( -np.power(data[:, :, 2] - target_hue[i], 2) / (2 * hue_sigma[i]**2))
+
+            weight_hue = multiplier[i] * weight_temp / np.max(weight_temp)
+
+            weight_chroma = np.exp( -np.power(data[:, :, 1] - chroma_preference[i], 2) / (2 * chroma_sigma[i]**2))
+
+            hue_correction = hue_correction + np.multiply(np.multiply(delta_hue, weight_hue), weight_chroma)
+
+        # correct the hue
+        data[:, :, 2] = data[:, :, 2] - hue_correction
+
+        # lch to lab
+        data = utility.color_conversion(data).lch2lab()
+        # lab to xyz
+        data = utility.color_conversion(data).lab2xyz(cie_version, illuminant)
+        # xyz to rgb
+        data = utility.color_conversion(data).xyz2rgb(color_space, clip_range)
+
+        return data
+
+
+    def __str__(self):
+        return self.name
+
+
+# =============================================================
+# class: chromatic_aberration_correction
+#   removes artifacts similar to result from chromatic
+#   aberration
+# =============================================================
+class chromatic_aberration_correction:
+    def __init__(self, data, name="chromatic aberration correction"):
+        self.data = np.float32(data)
+        self.name = name
+
+    def purple_fringe_removal(self, nsr_threshold, cr_threshold, clip_range=[0, 65535]):
+        # --------------------------------------------------------------
+        # nsr_threshold: near saturated region threshold (in percentage)
+        # cr_threshold: candidate region threshold
+        # --------------------------------------------------------------
+
+        width, height = utility.helpers(self.data).get_width_height()
+
+        r = self.data[:, :, 0]
+        g = self.data[:, :, 1]
+        b = self.data[:, :, 2]
+
+        ## Detection of purple fringe
+        # near saturated region detection
+        nsr_threshold = clip_range[1] * nsr_threshold / 100
+        temp = (r + g + b) / 3
+        temp = np.asarray(temp)
+        mask = temp > nsr_threshold
+        nsr = np.zeros((height, width)).astype(int)
+        nsr[mask] = 1
+
+        # candidate region detection
+        temp = r - b
+        temp1 = b - g
+        temp = np.asarray(temp)
+        temp1 = np.asarray(temp1)
+        mask = (temp < cr_threshold) & (temp1 > cr_threshold)
+        cr = np.zeros((height, width)).astype(int)
+        cr[mask] = 1
+
+        # quantization
+        qr = utility.helpers(r).nonuniform_quantization()
+        qg = utility.helpers(g).nonuniform_quantization()
+        qb = utility.helpers(b).nonuniform_quantization()
+
+        g_qr = utility.edge_detection(qr).sobel(5, "gradient_magnitude")
+        g_qg = utility.edge_detection(qg).sobel(5, "gradient_magnitude")
+        g_qb = utility.edge_detection(qb).sobel(5, "gradient_magnitude")
+
+        g_qr = np.asarray(g_qr)
+        g_qg = np.asarray(g_qg)
+        g_qb = np.asarray(g_qb)
+
+        # bgm: binary gradient magnitude
+        bgm = np.zeros((height, width), dtype=np.float32)
+        mask = (g_qr != 0) | (g_qg != 0) | (g_qb != 0)
+        bgm[mask] = 1
+
+        fringe_map = np.multiply(np.multiply(nsr, cr), bgm)
+        fring_map = np.asarray(fringe_map)
+        mask = (fringe_map == 1)
+
+        r1 = r
+        g1 = g
+        b1 = b
+        r1[mask] = g1[mask] = b1[mask] = (r[mask] + g[mask] + b[mask]) / 3.
+
+        output = np.empty(np.shape(self.data), dtype=np.float32)
+        output[:, :, 0] = r1
+        output[:, :, 1] = g1
+        output[:, :, 2] = b1
+
+        return np.float32(output)
+
+
+    def __str__(self):
+        return self.name

IIR-Lab/ISP_pipeline/lsc_table_r_gr_gb_b_2.npy

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0932e4b3ed5feb111880988eadae26a3cc77a5bd448150087b0983d8a62bb2b7
+size 402653312

IIR-Lab/ISP_pipeline/process_pngs_isp.py

@@ -0,0 +1,276 @@
+import sys
+sys.path.append('ISP_pipeline')
+from raw_prc_pipeline.pipeline import RawProcessingPipelineDemo
+import cv2
+import numpy as np
+import json
+import PIL.Image as Image
+import os,sys
+from raw_prc_pipeline import io
+from copy import deepcopy
+import torch
+
+def resize_using_pil(img, width=1024, height=768):
+    img_pil = Image.fromarray(img)
+    out_size = (width, height)
+    if img_pil.size == out_size:
+        return img
+    out_img = img_pil.resize(out_size, Image.LANCZOS)
+    # out_img = img_pil
+    out_img = np.array(out_img)
+    return out_img
+
+def fix_orientation(image, orientation):
+
+    if type(orientation) is list:
+        orientation = orientation[0]
+
+    if orientation == 'Horizontal(normal)':
+        pass
+    elif orientation == "Mirror horizontal":
+        image = cv2.flip(image, 0)
+    elif orientation == "Rotate 180":
+        image = cv2.rotate(image, cv2.ROTATE_180)
+    elif orientation == "Mirror vertical":
+        image = cv2.flip(image, 1)
+    elif orientation == "Mirror horizontal and rotate 270 CW":
+        image = cv2.flip(image, 0)
+        image = cv2.rotate(image, cv2.ROTATE_90_COUNTERCLOCKWISE)
+    elif orientation == "Rotate 90 CW":
+        image = cv2.rotate(image, cv2.ROTATE_90_CLOCKWISE)
+    elif orientation == "Mirror horizontal and rotate 90 CW":
+        image = cv2.flip(image, 0)
+        image = cv2.rotate(image, cv2.ROTATE_90_CLOCKWISE)
+    elif orientation == "Rotate 270 CW":
+        image = cv2.rotate(image, cv2.ROTATE_90_COUNTERCLOCKWISE)
+
+    return image
+
+def isp_night_imaging(data, meta_data, iso,
+                        do_demosaic = True,  # H/2 W/2
+
+                        do_channel_gain_white_balance = True,
+                        do_xyz_transform = True,
+                        do_srgb_transform = True,
+
+                        do_gamma_correct = True,  # con
+
+                        do_refinement = True,   # 32 bit
+                        do_to_uint8 = True,
+
+                        do_resize_using_pil = True,   # H/8, W/8
+                        do_fix_orientation = True
+                    ):
+    
+    pipeline_params = {
+        'tone_mapping': 'Flash', # options: Flash, Storm, Base, Linear, Drago, Mantiuk, Reinhard
+        'illumination_estimation': 'gw', # ie algorithm, options: "gw", "wp", "sog", "iwp"
+        'denoise_flg': True,
+        'out_landscape_width': 1024,
+        'out_landscape_height': 768,
+        "color_matrix": [ 1.06835938, -0.29882812, -0.14257812, 
+                        -0.43164062,  1.35546875,  0.05078125, 
+                        -0.1015625,   0.24414062,  0.5859375]
+    }
+
+    pipeline_demo = RawProcessingPipelineDemo(**pipeline_params)
+    
+    # ===================================
+    # Demosacing
+    # ===================================
+    if do_demosaic:
+        data = torch.stack((data[0,:,:], (data[1,:,:]+data[2,:,:])/2, data[3,:,:]), dim=0)
+        data = data.permute(1, 2, 0).contiguous()
+        # torch.cuda.empty_cache()
+    else:
+        pass
+    
+    # ===================================
+    # Channel gain for white balance
+    # ===================================
+    if do_channel_gain_white_balance:
+        data = pipeline_demo.white_balance(data, img_meta=meta_data)
+
+    else:
+        pass
+
+    # ===================================
+    # xyz_transform
+    # ===================================
+    if do_xyz_transform:
+        data = pipeline_demo.xyz_transform(data,img_meta=meta_data)   # CCM
+    else:
+        pass
+
+    # ===================================
+    # srgb_transform
+    # ===================================
+    if do_srgb_transform:
+        data = pipeline_demo.srgb_transform(data, img_meta=meta_data)   # fix ccm
+    else:
+        pass
+    
+    # ===================================
+    # gamma_correct
+    # ===================================
+    if do_gamma_correct:
+        data = pipeline_demo.gamma_correct(data, img_meta=meta_data)
+    else:
+        pass
+    
+    # ===================================
+    # refinement
+    # ===================================
+    if do_refinement:
+        if iso < 1000:
+            pth1 = "Rendering_models/low_iso.pth"
+            data = pipeline_demo.do_refinement(data, "csrnet", pth1)
+        else:
+            pth1 = "Rendering_models/high_iso.pth"
+            data = pipeline_demo.do_refinement(data, "csrnet", pth1)
+        torch.cuda.empty_cache()
+        
+    else:
+        pass
+
+    # ===================================
+    # to_uint8
+    # ===================================
+    if do_to_uint8:
+        data = pipeline_demo.to_uint8(data, img_meta=meta_data)
+        torch.cuda.empty_cache()
+    else:
+        pass
+
+    # ===================================
+    # resize_using_pil
+    # ===================================
+    if do_resize_using_pil:
+        data = resize_using_pil(data, pipeline_demo.params["out_landscape_width"], pipeline_demo.params["out_landscape_height"])
+
+    else:
+        pass
+
+    # ===================================
+    # fix_orientation
+    # ===================================
+    if do_fix_orientation:
+        data = fix_orientation(data, meta_data["orientation"])
+    else:
+        pass
+
+    return data
+
+def readjson(json_path,):
+    with open(json_path,'r',encoding='UTF-8') as f:
+        result = json.load(f)
+        
+    return result
+
+def get_smooth_kernel_size(factor):
+        if factor == 1:
+            return (5, 5)
+        elif factor == 0.5:
+            return (3, 3)
+        elif factor == 0.375:
+            return (3, 3)
+        elif factor in [0.2, 0.25]:
+            return (5, 5)
+        elif factor == 0.125:
+            return (7, 7)
+        else:
+            raise Exception('Unknown factor')
+
+def read_rawpng(path, metadata):
+    
+    raw = cv2.imread(str(path), cv2.IMREAD_UNCHANGED)
+
+    if raw.shape[0] == 4:
+        return raw * 959
+    raw = (raw.astype(np.float32) - 256.) / (4095.- 256.)
+    
+    raw = bayer2raw(raw, metadata)
+    raw = np.clip(raw, 0., 1.)
+    return raw
+
+def bayer2raw(raw, metadata):
+    # pack RGGB Bayer raw to 4 channels
+    H, W = raw.shape
+    raw = raw[None, ...]
+    if metadata['cfa_pattern'][0] == 0:
+        # RGGB
+        raw_pack = np.concatenate((raw[:, 0:H:2, 0:W:2],
+                                raw[:, 0:H:2, 1:W:2],
+                                raw[:, 1:H:2, 0:W:2],
+                                raw[:, 1:H:2, 1:W:2]), axis=0)
+    else :
+        # BGGR
+        raw_pack = np.concatenate((raw[:, 1:H:2, 1:W:2],
+                                raw[:, 0:H:2, 1:W:2],
+                                raw[:, 1:H:2, 0:W:2],
+                                raw[:, 0:H:2, 0:W:2]), axis=0)
+    return raw_pack
+
+def raw_rggb_float32(raws):
+    # depack 4 channels raw to RGGB Bayer
+    C, H, W = raws.shape
+    output = np.zeros((H * 2, W * 2)).astype(np.float32)
+
+    output[0:2 * H:2, 0:2 * W:2] = raws[0:1, :, :]
+    output[0:2 * H:2, 1:2 * W:2] = raws[1:2, :, :]
+    output[1:2 * H:2, 0:2 * W:2] = raws[2:3, :, :]
+    output[1:2 * H:2, 1:2 * W:2] = raws[3:4, :, :]
+
+    return output
+
+def json_read(pth):
+    with open(pth) as j:
+        data = json.load(j)
+    return data
+
+def linear_insert_1color(img_dt, resize, fx=128, fy=128):
+    pos_0_0, pos_0_1, pos_1_1, pos_1_0, m, n = insert_linear_pos(img_dt=img_dt, resize=resize, x_scale=fx, y_scale=fy)
+    a = (pos_1_0 - pos_0_0)
+    b = (pos_0_1 - pos_0_0)
+    c = pos_1_1 + pos_0_0 - pos_1_0 - pos_0_1
+    return np.round(a * n + b * m + c * n * m + pos_0_0).astype(int)
+
+def insert_linear_pos(img_dt, resize, x_scale=128, y_scale=128):
+    m_, n_ = img_dt.shape
+    # 获取新的图像的大小
+    if resize is None:
+        n_new, m_new  =  np.round(x_scale * n_).astype(int), np.round(y_scale * m_).astype(int)
+    else:
+        n_new, m_new  = resize
+
+    n_scale, m_scale = n_ / n_new, m_ / m_new # src_with/dst_with, Src_height/dst_heaight
+    # 一、获取位置对应的四个点
+    # 1-1- 初始化位置
+    m_indxs = np.repeat(np.arange(m_new), n_new).reshape(m_new, n_new)
+    n_indxs = np.array(list(range(n_new))*m_new).reshape(m_new, n_new)
+    # 1-2- 初始化位置
+    m_indxs_c = (m_indxs + 0.5 ) * m_scale - 0.5
+    n_indxs_c = (n_indxs + 0.5 ) * n_scale - 0.5
+    ### 将小于零的数处理成0 
+    m_indxs_c[np.where(m_indxs_c < 0)] = 0.0
+    n_indxs_c[np.where(n_indxs_c < 0)] = 0.0
+
+    # 1-3 获取正方形顶点坐标
+    m_indxs_c_down = m_indxs_c.astype(int)
+    n_indxs_c_down = n_indxs_c.astype(int)
+    m_indxs_c_up = m_indxs_c_down + 1
+    n_indxs_c_up = n_indxs_c_down + 1
+    ### 溢出部分修正
+    m_max = m_ - 1
+    n_max = n_ - 1
+    m_indxs_c_up[np.where(m_indxs_c_up > m_max)] = m_max
+    n_indxs_c_up[np.where(n_indxs_c_up > n_max)] = n_max
+
+    # 1-4 获取正方形四个顶点的位置
+    pos_0_0 = img_dt[m_indxs_c_down, n_indxs_c_down].astype(int)
+    pos_0_1 = img_dt[m_indxs_c_up, n_indxs_c_down].astype(int)
+    pos_1_1 = img_dt[m_indxs_c_up, n_indxs_c_up].astype(int)
+    pos_1_0 = img_dt[m_indxs_c_down, n_indxs_c_up].astype(int)
+    # 1-5 获取浮点位置
+    m, n = np.modf(m_indxs_c)[0], np.modf(n_indxs_c)[0]
+    return pos_0_0, pos_0_1, pos_1_1, pos_1_0, m, n

IIR-Lab/ISP_pipeline/raw_prc_pipeline/__init__.py

@@ -0,0 +1,3 @@
+expected_img_ext = '.jpg'
+expected_landscape_img_height = 866
+expected_landscape_img_width = 1300

IIR-Lab/ISP_pipeline/raw_prc_pipeline/arch_util.py

@@ -0,0 +1,626 @@
+import torch
+import torch.nn as nn
+import torch.nn.init as init
+import torch.nn.functional as F
+
+
+def initialize_weights(net_l, scale=1):
+    if not isinstance(net_l, list):
+        net_l = [net_l]
+    for net in net_l:
+        for m in net.modules():
+            if isinstance(m, nn.Conv2d):
+                init.kaiming_normal_(m.weight, a=0, mode='fan_in')
+                m.weight.data *= scale  # for residual block
+                if m.bias is not None:
+                    m.bias.data.zero_()
+            elif isinstance(m, nn.Linear):
+                init.kaiming_normal_(m.weight, a=0, mode='fan_in')
+                m.weight.data *= scale
+                if m.bias is not None:
+                    m.bias.data.zero_()
+            elif isinstance(m, nn.BatchNorm2d):
+                init.constant_(m.weight, 1)
+                init.constant_(m.bias.data, 0.0)
+
+
+def make_layer(block, n_layers):
+    layers = []
+    for _ in range(n_layers):
+        layers.append(block())
+    return nn.Sequential(*layers)
+
+
+class ResidualBlock_noBN(nn.Module):
+    '''Residual block w/o BN
+    ---Conv-ReLU-Conv-+-
+     |________________|
+    '''
+
+    def __init__(self, nf=64):
+        super(ResidualBlock_noBN, self).__init__()
+        self.conv1 = nn.Conv2d(nf, nf, 3, 1, 1, bias=True)
+        self.conv2 = nn.Conv2d(nf, nf, 3, 1, 1, bias=True)
+
+        # initialization
+        initialize_weights([self.conv1, self.conv2], 0.1)
+
+    def forward(self, x):
+        identity = x
+        out = F.relu(self.conv1(x), inplace=True)
+        out = self.conv2(out)
+        return identity + out
+
+
+def flow_warp(x, flow, interp_mode='bilinear', padding_mode='zeros'):
+    """Warp an image or feature map with optical flow
+    Args:
+        x (Tensor): size (N, C, H, W)
+        flow (Tensor): size (N, H, W, 2), normal value
+        interp_mode (str): 'nearest' or 'bilinear'
+        padding_mode (str): 'zeros' or 'border' or 'reflection'
+
+    Returns:
+        Tensor: warped image or feature map
+    """
+    assert x.size()[-2:] == flow.size()[1:3]
+    B, C, H, W = x.size()
+    # mesh grid
+    grid_y, grid_x = torch.meshgrid(torch.arange(0, H), torch.arange(0, W))
+    grid = torch.stack((grid_x, grid_y), 2).float()  # W(x), H(y), 2
+    grid.requires_grad = False
+    grid = grid.type_as(x)
+    vgrid = grid + flow
+    # scale grid to [-1,1]
+    vgrid_x = 2.0 * vgrid[:, :, :, 0] / max(W - 1, 1) - 1.0
+    vgrid_y = 2.0 * vgrid[:, :, :, 1] / max(H - 1, 1) - 1.0
+    vgrid_scaled = torch.stack((vgrid_x, vgrid_y), dim=3)
+    output = F.grid_sample(x, vgrid_scaled, mode=interp_mode, padding_mode=padding_mode)
+    return output
+
+"""
+Copyright (c) 2022 Samsung Electronics Co., Ltd.
+
+Author(s):
+Luxi Zhao ([email protected]; [email protected])
+Abdelrahman Abdelhamed ([email protected])
+
+Licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) License, (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at https://creativecommons.org/licenses/by-nc-sa/4.0
+Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an
+"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and limitations under the License.
+For conditions of distribution and use, see the accompanying LICENSE.md file.
+
+"""
+
+import numpy as np
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+def utils_get_image_stats(image_shape, grid_size):
+    """
+    Information about the cropped image.
+    :return: grid size, tile size, sizes of the 4 margins, meshgrids.
+    """
+
+    grid_rows = grid_size[0]
+    grid_cols = grid_size[1]
+
+    residual_height = image_shape[0] % grid_rows
+    residual_width = image_shape[1] % grid_cols
+
+    tile_height = image_shape[0] // grid_rows
+    tile_width = image_shape[1] // grid_cols
+
+    margin_top = tile_height // 2
+    margin_left = tile_width // 2
+
+    margin_bot = tile_height + residual_height - margin_top
+    margin_right = tile_width + residual_width - margin_left
+
+    return tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right
+
+def apply_ltm_lut(imgs, luts):
+
+        imgs = (imgs - .5) * 2.
+
+        grids = imgs.unsqueeze(0).unsqueeze(0)
+        luts = luts.unsqueeze(0)
+
+        outs = F.grid_sample(luts, grids,
+        mode='bilinear', padding_mode='border', align_corners=True)
+
+        return outs.squeeze(0).squeeze(1).permute(1,2,0)
+
+
+def apply_ltm(image, tone_curve, num_curves):
+    """
+    Apply tone curve to an image (patch).
+    :param image: (h, w, 3) if num_curves == 3, else (h, w)
+    :param tone_curve: (num_curves, control_points)
+    :param num_curves: 3 for 1 curve per channel, 1 for 1 curve for all channels.
+    :return: tone-mapped image.
+    """
+    
+    if image.shape[-1] == 3:
+        if type(image) == np.ndarray:
+            r = tone_curve[0][image[..., 0]]
+            g = tone_curve[1][image[..., 1]]
+            b = tone_curve[2][image[..., 2]]
+            new_image = np.stack((r, g, b), axis=-1)
+        else:
+            r = tone_curve[0][image[..., 0].reshape(-1).long()].reshape(image[..., 0].shape)
+            g = tone_curve[1][image[..., 1].reshape(-1).long()].reshape(image[..., 1].shape)
+            b = tone_curve[2][image[..., 2].reshape(-1).long()].reshape(image[..., 2].shape)
+            new_image = torch.stack((r, g, b), axis=-1)
+            # new_image = np.stack((r, g, b), axis=-1)
+    else:
+        new_image = tone_curve[0][image[..., 0].reshape(-1).long()].reshape(image[..., 0].shape).unsqueeze(dim=2)
+        #tone_curve[0][image[..., 0].reshape(-1).long()].reshape(image[..., 0].shape)
+    
+    return new_image
+
+
+def apply_gtm(image, tone_curve, num_curves):
+    """
+    Apply a single tone curve to an image.
+    :param image: (h, w, 3) if num_curves == 3, else (h, w)
+    :param tone_curve: (1, num_curves, control_points)
+    :param num_curves: 3 for 1 curve per channel, 1 for 1 curve for all channels.
+    :return: tone-mapped image.
+    """
+    tone_curve = tone_curve[0]
+    out = apply_ltm(image, tone_curve, num_curves)
+    return out
+
+
+def apply_ltm_center(image, tone_curves, stats, num_curves):
+    """
+    Apply tone curves to the center region of an image.
+    :param image: the original image.
+    :param tone_curves: a list of all tone curves in row scan order.
+    :return: interpolated center region of an image.
+    """
+    grid_rows, grid_cols, tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right, meshgrids = stats
+    xs_tl, ys_tl, xs_br, ys_br = meshgrids['center']
+
+    if torch.cuda.is_available():
+        device = torch.device("cuda")
+    else:
+        device = torch.device("cpu")
+    xs_tl = xs_tl.to(device)
+    ys_tl = ys_tl.to(device)
+    xs_br = xs_br.to(device)
+    ys_br = ys_br.to(device)
+
+
+    # Get neighbourhoods
+    neighbourhoods = []
+    for y in range(margin_top, image.shape[0]-margin_bot, tile_height):
+        for x in range(margin_left, image.shape[1]-margin_right, tile_width):
+            neighbourhoods.append(image[y:y + tile_height, x:x + tile_width, :])
+
+    assert len(neighbourhoods) == (grid_rows-1) * (grid_cols-1)
+
+    # Get indices for all 4-tile neighbourhoods
+    tile_ids = []
+    for i in range(grid_rows - 1):
+        for j in range(grid_cols - 1):
+            start = i * grid_cols + j
+            tile_ids.append([start, start + 1, start + grid_cols, start + grid_cols + 1])
+
+    # Apply LTM and interpolate
+    new_ns = []
+    for i, n in enumerate(neighbourhoods):
+        n_tile_ids = tile_ids[i]  # ids of the 4 tone curves (tiles) of the neighbourhood
+        # n_4versions = [apply_ltm(n, tone_curves[j], num_curves) for j in n_tile_ids]  # tl, tr, bl, br
+        n_4versions = [apply_ltm_lut(n, tone_curves[j])for j in n_tile_ids]
+        out = ys_br * xs_br * n_4versions[0] + ys_br * xs_tl * n_4versions[1] + ys_tl * xs_br * n_4versions[2] + ys_tl * xs_tl * n_4versions[3]
+        out /= (tile_height-1) * (tile_width-1)
+
+        new_ns.append(out)
+
+    # Stack the interpolated neighbourhoods together
+    rows = []
+    for i in range(grid_rows - 1):
+        cols = [new_ns[i * (grid_cols - 1) + j] for j in range(grid_cols - 1)]
+        row = torch.cat(cols, dim=1)
+        rows.append(row)
+    out = torch.cat(rows, dim=0)
+    return out
+
+
+def apply_ltm_border(image, tone_curves, stats, num_curves=3):
+    """
+    Apply tone curves to the border, not including corner areas.
+    :param image: the original image.
+    :param tone_curves: a list of all tone curves in row scan order.
+    :return: interpolated border regions of the image. In order of top, bottom, left, right.
+    """
+    grid_rows, grid_cols, tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right, meshgrids = stats
+    (top_xs_l, top_xs_r), (bot_xs_l, bot_xs_r), (left_ys_t, left_ys_b), (right_ys_t, right_ys_b) = meshgrids['border']
+
+    if torch.cuda.is_available():
+        device = torch.device("cuda")
+    else:
+        device = torch.device("cpu")
+
+    top_xs_l = top_xs_l.to(device) 
+    top_xs_r = top_xs_r.to(device)
+    bot_xs_l = bot_xs_l.to(device)
+    bot_xs_r = bot_xs_r.to(device) 
+
+    left_ys_t = left_ys_t.to(device)
+    left_ys_b = left_ys_b.to(device)  
+    right_ys_t = right_ys_t.to(device)
+    right_ys_b = right_ys_b.to(device)
+    # top, bottom, left, right neighbourhoods to be interpolated
+    ntop = []
+    nbot = []
+    nleft = []
+    nright = []
+
+    for x in range(margin_left, image.shape[1] - margin_right, tile_width):
+        ntop.append(image[:margin_top, x:x + tile_width, :])
+        nbot.append(image[-margin_bot:, x:x + tile_width, :])
+
+    for y in range(margin_top, image.shape[0] - margin_bot, tile_height):
+        nleft.append(image[y:y + tile_height, :margin_left, :])
+        nright.append(image[y:y + tile_height, -margin_right:, :])
+
+    def apply_ltm_two_tiles(tc1, tc2, meshgrid1, meshgrid2, nbhd, interp_length, num_curves):
+        """
+        Apply tone curve to, and interpolate a two-tile neighbourhood, either horizontal or vertical
+        :param tc1: left / top tone curves
+        :param tc2: right / bottom tone curves
+        :param meshgrid1: left / top meshgrids (leftmost / topmost positions are 0)
+        :param meshgrid2: right / bottom meshgrids (rightmost / bottommost positions are 0)
+        :param nbhd: neighbourhood to interpolate
+        :param interp_length: normalizing factor of the meshgrid.
+               Example: if xs = np.meshgrid(np.arange(10)), then interp_length = 9
+        :return: interpolated neighbourhood
+        """
+
+        # new_nbhd1 = apply_ltm(nbhd, tc1, num_curves)
+        # new_nbhd2 = apply_ltm(nbhd, tc2, num_curves)
+
+        new_nbhd1 = apply_ltm_lut(nbhd, tc1)
+        new_nbhd2 = apply_ltm_lut(nbhd, tc2)
+
+        out = meshgrid1 * new_nbhd2 + meshgrid2 * new_nbhd1
+        out /= interp_length
+        return out
+
+    new_ntop = [apply_ltm_two_tiles(tone_curves[i],  # left tone curve
+                                    tone_curves[i + 1],  # right tone curve
+                                    top_xs_l, top_xs_r,
+                                    n, tile_width - 1, num_curves) for i, n in enumerate(ntop)]
+
+    new_nbot = [apply_ltm_two_tiles(tone_curves[(grid_rows - 1) * grid_cols + i],  # left tone curve
+                                    tone_curves[(grid_rows - 1) * grid_cols + i + 1],  # right tone curve
+                                    bot_xs_l, bot_xs_r,
+                                    n, tile_width - 1, num_curves) for i, n in enumerate(nbot)]
+
+    new_nleft = [apply_ltm_two_tiles(tone_curves[i * grid_cols],  # top tone curve
+                                     tone_curves[(i + 1) * grid_cols],  # bottom tone curve
+                                     left_ys_t, left_ys_b,
+                                     n, tile_height - 1, num_curves) for i, n in enumerate(nleft)]
+
+    new_nright = [apply_ltm_two_tiles(tone_curves[(i + 1) * grid_cols - 1],  # top tone curve
+                                      tone_curves[(i + 2) * grid_cols - 1],  # bottom tone curve
+                                      right_ys_t, right_ys_b,
+                                      n, tile_height - 1, num_curves) for i, n in enumerate(nright)]
+
+    new_ntop = torch.cat(new_ntop, dim=1)
+    new_nbot = torch.cat(new_nbot, dim=1)
+    new_nleft = torch.cat(new_nleft, dim=0)
+    new_nright = torch.cat(new_nright, dim=0)
+    return new_ntop, new_nbot, new_nleft, new_nright
+
+
+def apply_ltm_corner(image, tone_curves, stats, num_curves=3):
+    """
+    tone_curves: a list of all tone curves in row scan order.
+    return: interpolated corner tiles in the order of top left, top right, bot left, bot right
+    """
+    grid_rows, grid_cols, tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right, _ = stats
+
+    corner_ids = [0, grid_cols - 1, -grid_rows, -1]
+    tl_tile = image[:margin_top, :margin_left]
+    tr_tile = image[:margin_top, -margin_right:]
+    bl_tile = image[-margin_bot:, :margin_left]
+    br_tile = image[-margin_bot:, -margin_right:]
+
+    corner_tiles = [tl_tile, tr_tile, bl_tile, br_tile]
+    corner_tcs = [tone_curves[i] for i in corner_ids]  # tcs: (grid_size, num_curves, control_points)
+    #new_tiles = [apply_ltm(corner_tiles[i], corner_tcs[i], num_curves) for i in range(len(corner_tcs))]
+    new_tiles = [apply_ltm_lut(corner_tiles[i],corner_tcs[i]) for i in range(len(corner_tcs))]
+
+    return new_tiles[0], new_tiles[1], new_tiles[2], new_tiles[3]
+
+
+# def get_meshgrids(height, width):
+#     """
+#     Get two meshgrids of size (height, width). One with top left corner being (0, 0),
+#     the other with bottom right corner being (0, 0).
+#     :return: top left xs, ys, bottom right xs, ys
+#     """
+#     xs, ys = np.meshgrid(np.arange(width), np.arange(height))
+#     newys, newxs = torch.meshgrid(torch.arange(height, dtype=torch.int32), torch.arange(width, dtype=torch.int32))
+#     # mesh grid for top left corner
+#     xs_tl = np.tile(np.abs(xs)[..., np.newaxis], 3)  # [0, 1, 2, ..., tile_width-1]
+#     ys_tl = np.tile(np.abs(ys)[..., np.newaxis], 3)
+#     new_xs_tl = newxs[..., None].abs().repeat(1, 1, 3)
+#     new_ys_tl = newys[..., None].abs().repeat(1, 1, 3)
+#     # mesh grid for bottom right corner
+#     xs_br = np.tile(np.abs(xs - width + 1)[..., np.newaxis], 3)  # [-(tile_width-1), ..., -2, -1, 0]
+#     ys_br = np.tile(np.abs(ys - height + 1)[..., np.newaxis], 3)
+
+#     new_xs_br = (newxs - width + 1).abs()[..., None].repeat(1, 1, 3)
+#     new_ys_br = (newys - width + 1).abs()[..., None].repeat(1, 1, 3)
+#     # return xs_tl, ys_tl, xs_br, ys_br
+#     return new_xs_tl, new_ys_tl, new_xs_br, new_ys_br
+def get_meshgrids(height, width):
+    """
+    Get two meshgrids of size (height, width). One with top left corner being (0, 0),
+    the other with bottom right corner being (0, 0).
+    :return: top left xs, ys, bottom right xs, ys
+    """
+    xs, ys = np.meshgrid(np.arange(width), np.arange(height))
+    # mesh grid for top left corner
+    xs_tl = np.tile(np.abs(xs)[..., np.newaxis], 3)  # [0, 1, 2, ..., tile_width-1]
+    ys_tl = np.tile(np.abs(ys)[..., np.newaxis], 3)
+    # mesh grid for bottom right corner
+    xs_br = np.tile(np.abs(xs - width + 1)[..., np.newaxis], 3)  # [-(tile_width-1), ..., -2, -1, 0]
+    ys_br = np.tile(np.abs(ys - height + 1)[..., np.newaxis], 3)
+    
+    return torch.tensor(xs_tl), torch.tensor(ys_tl), torch.tensor(xs_br), torch.tensor(ys_br)
+
+
+
+def get_meshgrid_center(tile_height, tile_width):
+    return get_meshgrids(tile_height, tile_width)
+
+
+def get_meshgrid_border(tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right):
+    """
+    :return: meshgrids for the 4 border regions, in the order of top, bottom, left, right
+    """
+    # top
+    top_xs_l, _, top_xs_r, _ = get_meshgrids(margin_top, tile_width)
+
+    # bottom
+    bot_xs_l, _, bot_xs_r, _ = get_meshgrids(margin_bot, tile_width)
+
+    # left
+    _, left_ys_t, _, left_ys_b = get_meshgrids(tile_height, margin_left)
+
+    # right
+    _, right_ys_t, _, right_ys_b = get_meshgrids(tile_height, margin_right)
+
+    return (top_xs_l, top_xs_r), (bot_xs_l, bot_xs_r), (left_ys_t, left_ys_b), (right_ys_t, right_ys_b)
+
+
+def get_image_stats(image, grid_size):
+    """
+    Information about the cropped image.
+    :param image: the original image
+    :return: grid size, tile size, sizes of the 4 margins, meshgrids.
+    """
+
+    grid_rows = grid_size[0]
+    grid_cols = grid_size[1]
+
+    tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right = utils_get_image_stats(image.shape,
+                                                                                                       grid_size)
+
+    meshgrid_center = get_meshgrid_center(tile_height, tile_width)
+    meshgrid_border = get_meshgrid_border(tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right)
+
+    meshgrids = {
+        'center': meshgrid_center,
+        'border': meshgrid_border
+    }
+
+    return grid_rows, grid_cols, tile_height, tile_width, margin_top, margin_left, margin_bot, margin_right, meshgrids
+
+
+#自己构造image 1 * 512 * 512 *3,tone_curve 1 * 64 * 3 * 256 维度得tf tensor 然后只在这里debug  
+def do_interpolation_lut(image, tone_curves, grid_size, num_curves=3):
+    """
+    Perform tone mapping and interpolation on an image.
+    Center region: bilinear interpolation.
+    Border region: linear interpolation.
+    Corner region: no interpolation.
+    :param num_curves: 3 -> 1 curve for each R,G,B channel, 1 -> 1 curve for all channels
+    :param image: input int8
+    :param tone_curves: (grid_size, num_curves, control_points)
+    :param grid_size: (ncols, nrows)
+    :return: image: float32, between [0-1]
+    """
+    if grid_size[0] == 1 and grid_size[1] == 1:
+        return apply_gtm(image, tone_curves, num_curves).astype(np.float64)
+
+    # get image statistics
+    stats = get_image_stats(image, grid_size)
+
+
+
+    # Center area:
+    center = apply_ltm_center(image, tone_curves, stats, num_curves)
+
+    # Border area:
+    b_top, b_bot, b_left, b_right = apply_ltm_border(image, tone_curves, stats, num_curves)
+
+    # Corner area:
+    tlc, trc, blc, brc = apply_ltm_corner(image, tone_curves, stats, num_curves)
+
+    # stack the corners, borders, and center together
+    row_t = torch.cat([tlc, b_top, trc], dim=1)
+    row_c = torch.cat([b_left, center, b_right], dim=1)
+    row_b = torch.cat([blc, b_bot, brc], dim=1)
+    out = torch.cat([row_t, row_c, row_b], dim=0)
+
+    assert out.shape == image.shape
+
+    return out
+
+
+
+class Model(nn.Module):
+    def __init__(self):
+        super(Model, self).__init__()
+        # self.conv1 = nn.Conv2d(3, 4, kernel_size=3, padding=1)  
+        # self.pool1 = nn.MaxPool2d(2)
+        
+        # self.conv2 = nn.Conv2d(4, 8, kernel_size=3, padding=1)
+        # self.pool2 = nn.MaxPool2d(2) 
+        
+        # self.conv3 = nn.Conv2d(8, 16, kernel_size=3, padding=1)
+        # self.pool3 = nn.MaxPool2d(2)
+
+        # self.conv4 = nn.Conv2d(16, 32, kernel_size=3, padding=1) 
+        # self.pool4 = nn.MaxPool2d(2)
+
+        # self.conv5 = nn.Conv2d(32, 64, kernel_size=3, padding=1)
+        # self.pool5 = nn.MaxPool2d(2)
+        
+        # self.conv6 = nn.Conv2d(64, 768, kernel_size=3, padding=1)
+        # self.pool6 = nn.MaxPool2d(2)
+
+        self.layer_1 = nn.Sequential(
+            nn.Conv2d(3, 4, kernel_size=3, padding=1),
+            nn.BatchNorm2d(4),  
+            nn.ReLU(),
+            nn.MaxPool2d(2)
+            )
+        self.layer_2 = nn.Sequential(
+            nn.Conv2d(4, 8, kernel_size=3, padding=1),
+            nn.BatchNorm2d(8),  
+            nn.ReLU(),
+            nn.MaxPool2d(2)
+        )
+        self.layer_3 = nn.Sequential(
+            nn.Conv2d(8, 16, kernel_size=3, padding=1),
+            nn.BatchNorm2d(16),  
+            nn.ReLU(),
+            nn.MaxPool2d(2)
+        )
+        self.layer_4 = nn.Sequential(
+            nn.Conv2d(16, 32, kernel_size=3, padding=1),
+            nn.BatchNorm2d(32),  
+            nn.ReLU(),
+            nn.MaxPool2d(2)
+        )
+        self.layer_5 = nn.Sequential(
+            nn.Conv2d(32, 64, kernel_size=3, padding=1),
+            nn.BatchNorm2d(64),  
+            nn.ReLU(),
+            nn.MaxPool2d(2)
+        )
+        self.layer_6 = nn.Sequential(
+            nn.Conv2d(64, 768, kernel_size=3, padding=1),
+            nn.BatchNorm2d(768),  
+            nn.Sigmoid(),
+            nn.MaxPool2d(2)
+        )
+
+        
+    def forward(self, x):
+        
+        '''
+
+        original = x
+        x = self.conv1(x)
+        x = self.pool1(x)
+        
+        x = self.conv2(x) 
+        x = self.pool2(x)
+
+        x = self.conv3(x)
+        x = self.pool3(x)
+        
+        x = self.conv4(x)
+        x = self.pool4(x)
+
+        x = self.conv5(x) 
+        x = self.pool5(x)
+
+        x = self.conv6(x)
+        x = self.pool6(x)
+        oldres = x
+        x = original
+        '''
+
+        x = self.layer_1(x)
+
+        x = self.layer_2(x)
+
+        x = self.layer_3(x)
+
+        x = self.layer_4(x)
+
+        x = self.layer_5(x)
+
+        x = self.layer_6(x)
+
+        x = x.reshape(x.shape[0], x.shape[2] * x.shape[3], 3, int(x.shape[1] / 3))
+        return x
+
+
+def _lut_transform(imgs, luts):
+    # img (b, 3, h, w), lut (b, c, m, m, m)
+    if imgs.shape[1]==1:
+
+        #for gray image pro-processs
+        luts = luts.expand(1,1,64,64,64)
+        # normalize pixel values
+        imgs = (imgs - .5) * 2.
+        grids = (imgs.unsqueeze(4)).repeat(1,1,1,1,3)
+    else:
+        # normalize pixel values
+        imgs = (imgs - .5) * 2.
+        # reshape img to grid of shape (b, 1, h, w, 3)
+        # imgs = imgs.permute(2,0,1).unsqueeze(dim=0)
+        # grids = imgs.permute(0, 2, 3, 1).unsqueeze(1)
+        grids = imgs.unsqueeze(0).unsqueeze(0)
+        luts = luts.unsqueeze(0)
+        # after gridsampling, output is of shape (b, c, 1, h, w)
+    outs = F.grid_sample(luts, grids,
+        mode='bilinear', padding_mode='border', align_corners=True)
+    return outs.squeeze(2)
+
+
+if __name__ == '__main__':
+
+    import torch
+    import cv2
+
+    grid_size = [8,8]
+    
+    np.random.seed(42)
+    rand_img  = np.random.random((512, 512, 3))
+    luts_np = np.random.random((64, 3, 9))
+
+    img_torch = torch.tensor(rand_img, dtype=torch.float32).cuda()
+    luts_torch = torch.tensor(luts_np, dtype=torch.float32).cuda()
+
+
+    iluts = []
+    for i in range(luts_torch.shape[0]):
+        iluts.append(torch.stack(
+            torch.meshgrid(*(luts_torch[i].unbind(0)[::-1])),
+            dim=0).flip(0))
+    iluts = torch.stack(iluts, dim=0)
+ 
+
+    result = do_interpolation_lut(img_torch, iluts, grid_size)
+    print(result)
+
+
+
+
+
+

IIR-Lab/ISP_pipeline/raw_prc_pipeline/color.py

@@ -0,0 +1,306 @@
+import numpy as np
+
+
+def rgb2gray(data):
+    return 0.299 * data[:, :, 0] + \
+           0.587 * data[:, :, 1] + \
+           0.114 * data[:, :, 2]
+
+
+def rgb2ycc(data, rule="bt601"):
+    # map to select kr and kb
+    kr_kb_dict = {"bt601": [0.299, 0.114],
+                  "bt709": [0.2126, 0.0722],
+                  "bt2020": [0.2627, 0.0593]}
+
+    kr = kr_kb_dict[rule][0]
+    kb = kr_kb_dict[rule][1]
+    kg = 1 - (kr + kb)
+
+    output = np.empty(np.shape(data), dtype=np.float32)
+    output[:, :, 0] = kr * data[:, :, 0] + \
+                      kg * data[:, :, 1] + \
+                      kb * data[:, :, 2]
+    output[:, :, 1] = 0.5 * ((data[:, :, 2] - output[:, :, 0]) / (1 - kb))
+    output[:, :, 2] = 0.5 * ((data[:, :, 0] - output[:, :, 0]) / (1 - kr))
+
+    return output
+
+
+def ycc2rgb(data, rule="bt601"):
+    # map to select kr and kb
+    kr_kb_dict = {"bt601": [0.299, 0.114],
+                  "bt709": [0.2126, 0.0722],
+                  "bt2020": [0.2627, 0.0593]}
+
+    kr = kr_kb_dict[rule][0]
+    kb = kr_kb_dict[rule][1]
+    kg = 1 - (kr + kb)
+
+    output = np.empty(np.shape(data), dtype=np.float32)
+    output[:, :, 0] = 2. * data[:, :, 2] * (1 - kr) + data[:, :, 0]
+    output[:, :, 2] = 2. * data[:, :, 1] * (1 - kb) + data[:, :, 0]
+    output[:, :, 1] = (data[:, :, 0] - kr * output[:, :, 0] - kb * output[:, :, 2]) / kg
+
+    return output
+
+
+def degamma_srgb(data, clip_range=[0, 65535]):
+    # bring data in range 0 to 1
+    data = np.clip(data, clip_range[0], clip_range[1])
+    data = np.divide(data, clip_range[1])
+
+    data = np.asarray(data)
+    mask = data > 0.04045
+
+    # basically, if data[x, y, c] > 0.04045, data[x, y, c] = ( (data[x, y, c] + 0.055) / 1.055 ) ^ 2.4
+    #            else, data[x, y, c] = data[x, y, c] / 12.92
+    data[mask] += 0.055
+    data[mask] /= 1.055
+    data[mask] **= 2.4
+
+    data[np.invert(mask)] /= 12.92
+
+    # rescale
+    return np.clip(data * clip_range[1], clip_range[0], clip_range[1])
+
+
+def degamma_adobe_rgb_1998(data, clip_range=[0, 65535]):
+    # bring data in range 0 to 1
+    data = np.clip(data, clip_range[0], clip_range[1])
+    data = np.divide(data, clip_range[1])
+
+    data = np.power(data, 2.2)  # originally raised to 2.19921875
+
+    # rescale
+    return np.clip(data * clip_range[1], clip_range[0], clip_range[1])
+
+
+def rgb2xyz(data, color_space="srgb", clip_range=[0, 255]):
+    # input rgb in range clip_range
+    # output xyz is in range 0 to 1
+    if color_space == "srgb":
+        # degamma / linearization
+        data = degamma_srgb(data, clip_range)
+        data = np.float32(data)
+        data = np.divide(data, clip_range[1])
+
+        # matrix multiplication`
+        output = np.empty(np.shape(data), dtype=np.float32)
+        output[:, :, 0] = data[:, :, 0] * 0.4124 + data[:, :, 1] * 0.3576 + data[:, :, 2] * 0.1805
+        output[:, :, 1] = data[:, :, 0] * 0.2126 + data[:, :, 1] * 0.7152 + data[:, :, 2] * 0.0722
+        output[:, :, 2] = data[:, :, 0] * 0.0193 + data[:, :, 1] * 0.1192 + data[:, :, 2] * 0.9505
+    elif color_space == "adobe-rgb-1998":
+        # degamma / linearization
+        data = degamma_adobe_rgb_1998(data, clip_range)
+        data = np.float32(data)
+        data = np.divide(data, clip_range[1])
+
+        # matrix multiplication
+        output = np.empty(np.shape(data), dtype=np.float32)
+        output[:, :, 0] = data[:, :, 0] * 0.5767309 + data[:, :, 1] * 0.1855540 + data[:, :, 2] * 0.1881852
+        output[:, :, 1] = data[:, :, 0] * 0.2973769 + data[:, :, 1] * 0.6273491 + data[:, :, 2] * 0.0752741
+        output[:, :, 2] = data[:, :, 0] * 0.0270343 + data[:, :, 1] * 0.0706872 + data[:, :, 2] * 0.9911085
+    elif color_space == "linear":
+        # matrix multiplication`
+        output = np.empty(np.shape(data), dtype=np.float32)
+        data = np.float32(data)
+        data = np.divide(data, clip_range[1])
+        output[:, :, 0] = data[:, :, 0] * 0.4124 + data[:, :, 1] * 0.3576 + data[:, :, 2] * 0.1805
+        output[:, :, 1] = data[:, :, 0] * 0.2126 + data[:, :, 1] * 0.7152 + data[:, :, 2] * 0.0722
+        output[:, :, 2] = data[:, :, 0] * 0.0193 + data[:, :, 1] * 0.1192 + data[:, :, 2] * 0.9505
+    else:
+        print("Warning! color_space must be srgb or adobe-rgb-1998.")
+        return
+
+    return output
+
+
+def gamma_srgb(data, clip_range=[0, 65535]):
+    # bring data in range 0 to 1
+    data = np.clip(data, clip_range[0], clip_range[1])
+    data = np.divide(data, clip_range[1])
+
+    data = np.asarray(data)
+    mask = data > 0.0031308
+
+    # basically, if data[x, y, c] > 0.0031308, data[x, y, c] = 1.055 * ( var_R(i, j) ^ ( 1 / 2.4 ) ) - 0.055
+    #            else, data[x, y, c] = data[x, y, c] * 12.92
+    data[mask] **= 0.4167
+    data[mask] *= 1.055
+    data[mask] -= 0.055
+
+    data[np.invert(mask)] *= 12.92
+
+    # rescale
+    return np.clip(data * clip_range[1], clip_range[0], clip_range[1])
+
+
+def gamma_adobe_rgb_1998(data, clip_range=[0, 65535]):
+    # bring data in range 0 to 1
+    data = np.clip(data, clip_range[0], clip_range[1])
+    data = np.divide(data, clip_range[1])
+
+    data = np.power(data, 0.4545)
+
+    # rescale
+    return np.clip(data * clip_range[1], clip_range[0], clip_range[1])
+
+
+def xyz2rgb(data, color_space="srgb", clip_range=[0, 255]):
+    # input xyz is in range 0 to 1
+    # output rgb in clip_range
+
+    # allocate space for output
+    output = np.empty(np.shape(data), dtype=np.float32)
+
+    if color_space == "srgb":
+        # matrix multiplication
+        output[:, :, 0] = data[:, :, 0] * 3.2406 + data[:, :, 1] * -1.5372 + data[:, :, 2] * -0.4986
+        output[:, :, 1] = data[:, :, 0] * -0.9689 + data[:, :, 1] * 1.8758 + data[:, :, 2] * 0.0415
+        output[:, :, 2] = data[:, :, 0] * 0.0557 + data[:, :, 1] * -0.2040 + data[:, :, 2] * 1.0570
+
+        # gamma to retain nonlinearity
+        output = gamma_srgb(output * clip_range[1], clip_range)
+    elif color_space == "adobe-rgb-1998":
+        # matrix multiplication
+        output[:, :, 0] = data[:, :, 0] * 2.0413690 + data[:, :, 1] * -0.5649464 + data[:, :, 2] * -0.3446944
+        output[:, :, 1] = data[:, :, 0] * -0.9692660 + data[:, :, 1] * 1.8760108 + data[:, :, 2] * 0.0415560
+        output[:, :, 2] = data[:, :, 0] * 0.0134474 + data[:, :, 1] * -0.1183897 + data[:, :, 2] * 1.0154096
+
+        # gamma to retain nonlinearity
+        output = gamma_adobe_rgb_1998(output * clip_range[1], clip_range)
+    elif color_space == "linear":
+
+        # matrix multiplication
+        output[:, :, 0] = data[:, :, 0] * 3.2406 + data[:, :, 1] * -1.5372 + data[:, :, 2] * -0.4986
+        output[:, :, 1] = data[:, :, 0] * -0.9689 + data[:, :, 1] * 1.8758 + data[:, :, 2] * 0.0415
+        output[:, :, 2] = data[:, :, 0] * 0.0557 + data[:, :, 1] * -0.2040 + data[:, :, 2] * 1.0570
+
+        # gamma to retain nonlinearity
+        output = output * clip_range[1]
+    else:
+        print("Warning! color_space must be srgb or adobe-rgb-1998.")
+        return
+
+    return output
+
+
+def get_xyz_reference(cie_version="1931", illuminant="d65"):
+    if cie_version == "1931":
+        xyz_reference_dictionary = {"A": [109.850, 100.0, 35.585],
+                                    "B": [99.0927, 100.0, 85.313],
+                                    "C": [98.074, 100.0, 118.232],
+                                    "d50": [96.422, 100.0, 82.521],
+                                    "d55": [95.682, 100.0, 92.149],
+                                    "d65": [95.047, 100.0, 108.883],
+                                    "d75": [94.972, 100.0, 122.638],
+                                    "E": [100.0, 100.0, 100.0],
+                                    "F1": [92.834, 100.0, 103.665],
+                                    "F2": [99.187, 100.0, 67.395],
+                                    "F3": [103.754, 100.0, 49.861],
+                                    "F4": [109.147, 100.0, 38.813],
+                                    "F5": [90.872, 100.0, 98.723],
+                                    "F6": [97.309, 100.0, 60.191],
+                                    "F7": [95.044, 100.0, 108.755],
+                                    "F8": [96.413, 100.0, 82.333],
+                                    "F9": [100.365, 100.0, 67.868],
+                                    "F10": [96.174, 100.0, 81.712],
+                                    "F11": [100.966, 100.0, 64.370],
+                                    "F12": [108.046, 100.0, 39.228]}
+    elif cie_version == "1964":
+        xyz_reference_dictionary = {"A": [111.144, 100.0, 35.200],
+                                    "B": [99.178, 100.0, 84.3493],
+                                    "C": [97.285, 100.0, 116.145],
+                                    "D50": [96.720, 100.0, 81.427],
+                                    "D55": [95.799, 100.0, 90.926],
+                                    "D65": [94.811, 100.0, 107.304],
+                                    "D75": [94.416, 100.0, 120.641],
+                                    "E": [100.0, 100.0, 100.0],
+                                    "F1": [94.791, 100.0, 103.191],
+                                    "F2": [103.280, 100.0, 69.026],
+                                    "F3": [108.968, 100.0, 51.965],
+                                    "F4": [114.961, 100.0, 40.963],
+                                    "F5": [93.369, 100.0, 98.636],
+                                    "F6": [102.148, 100.0, 62.074],
+                                    "F7": [95.792, 100.0, 107.687],
+                                    "F8": [97.115, 100.0, 81.135],
+                                    "F9": [102.116, 100.0, 67.826],
+                                    "F10": [99.001, 100.0, 83.134],
+                                    "F11": [103.866, 100.0, 65.627],
+                                    "F12": [111.428, 100.0, 40.353]}
+    else:
+        print("Warning! cie_version must be 1931 or 1964.")
+        return
+    return np.divide(xyz_reference_dictionary[illuminant], 100.0)
+
+
+def xyz2lab(data, cie_version="1931", illuminant="d65"):
+    xyz_reference = get_xyz_reference(cie_version, illuminant)
+
+    data = data
+    data[:, :, 0] = data[:, :, 0] / xyz_reference[0]
+    data[:, :, 1] = data[:, :, 1] / xyz_reference[1]
+    data[:, :, 2] = data[:, :, 2] / xyz_reference[2]
+
+    data = np.asarray(data)
+
+    # if data[x, y, c] > 0.008856, data[x, y, c] = data[x, y, c] ^ (1/3)
+    # else, data[x, y, c] = 7.787 * data[x, y, c] + 16/116
+    mask = data > 0.008856
+    data[mask] **= 1. / 3.
+    data[np.invert(mask)] *= 7.787
+    data[np.invert(mask)] += 16. / 116.
+
+    data = np.float32(data)
+    output = np.empty(np.shape(data), dtype=np.float32)
+    output[:, :, 0] = 116. * data[:, :, 1] - 16.
+    output[:, :, 1] = 500. * (data[:, :, 0] - data[:, :, 1])
+    output[:, :, 2] = 200. * (data[:, :, 1] - data[:, :, 2])
+
+    return output
+
+
+def lab2xyz(data, cie_version="1931", illuminant="d65"):
+    output = np.empty(np.shape(data), dtype=np.float32)
+
+    output[:, :, 1] = (data[:, :, 0] + 16.) / 116.
+    output[:, :, 0] = (data[:, :, 1] / 500.) + output[:, :, 1]
+    output[:, :, 2] = output[:, :, 1] - (data[:, :, 2] / 200.)
+
+    # if output[x, y, c] > 0.008856, output[x, y, c] ^ 3
+    # else, output[x, y, c] = ( output[x, y, c] - 16/116 ) / 7.787
+    output = np.asarray(output)
+    mask = output > 0.008856
+    output[mask] **= 3.
+    output[np.invert(mask)] -= 16 / 116
+    output[np.invert(mask)] /= 7.787
+
+    xyz_reference = get_xyz_reference(cie_version, illuminant)
+
+    output = np.float32(output)
+    output[:, :, 0] = output[:, :, 0] * xyz_reference[0]
+    output[:, :, 1] = output[:, :, 1] * xyz_reference[1]
+    output[:, :, 2] = output[:, :, 2] * xyz_reference[2]
+
+    return output
+
+
+def lab2lch(data):
+    output = np.empty(np.shape(data), dtype=np.float32)
+
+    output[:, :, 0] = data[:, :, 0]  # L transfers directly
+    output[:, :, 1] = np.power(np.power(data[:, :, 1], 2) + np.power(data[:, :, 2], 2), 0.5)
+    output[:, :, 2] = np.arctan2(data[:, :, 2], data[:, :, 1]) * 180 / np.pi
+
+    return output
+
+
+def lch2lab(data):
+    output = np.empty(np.shape(data), dtype=np.float32)
+
+    output[:, :, 0] = data[:, :, 0]  # L transfers directly
+    output[:, :, 1] = np.multiply(np.cos(data[:, :, 2] * np.pi / 180), data[:, :, 1])
+    output[:, :, 2] = np.multiply(np.sin(data[:, :, 2] * np.pi / 180), data[:, :, 1])
+
+    return output

IIR-Lab/ISP_pipeline/raw_prc_pipeline/csrnet_network.py

@@ -0,0 +1,76 @@
+import functools
+import math
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+
+
+class Condition(nn.Module):
+    def __init__(self, in_nc=3, nf=32):
+        super(Condition, self).__init__()
+        stride = 2
+        pad = 0
+        self.pad = nn.ZeroPad2d(1)
+        self.conv1 = nn.Conv2d(in_nc, nf, 7, stride, pad, bias=True)
+        self.conv2 = nn.Conv2d(nf, nf, 3, stride, pad, bias=True)
+        self.conv3 = nn.Conv2d(nf, nf, 3, stride, pad, bias=True)
+        self.act = nn.ReLU(inplace=True)
+
+    def forward(self, x):
+        conv1_out = self.act(self.conv1(self.pad(x)))
+        conv2_out = self.act(self.conv2(self.pad(conv1_out)))
+        conv3_out = self.act(self.conv3(self.pad(conv2_out)))
+        out = torch.mean(conv3_out, dim=[2, 3], keepdim=False)
+
+        return out
+
+
+# 3layers with control
+class CSRNet(nn.Module):
+    def __init__(self, in_nc=3, out_nc=3, base_nf=48, cond_nf=24):
+        super(CSRNet, self).__init__()
+
+        self.base_nf = base_nf
+        self.out_nc = out_nc
+
+        self.cond_net = Condition(in_nc=in_nc, nf=cond_nf)
+      
+        self.cond_scale1 = nn.Linear(cond_nf, base_nf, bias=True)
+        self.cond_scale2 = nn.Linear(cond_nf, base_nf,  bias=True)
+        self.cond_scale3 = nn.Linear(cond_nf, 3, bias=True)
+
+        self.cond_shift1 = nn.Linear(cond_nf, base_nf, bias=True)
+        self.cond_shift2 = nn.Linear(cond_nf, base_nf, bias=True)
+        self.cond_shift3 = nn.Linear(cond_nf, 3, bias=True)
+
+        self.conv1 = nn.Conv2d(in_nc, base_nf, 1, 1, bias=True) 
+        self.conv2 = nn.Conv2d(base_nf, base_nf, 1, 1, bias=True)
+        self.conv3 = nn.Conv2d(base_nf, out_nc, 1, 1, bias=True)
+
+        self.act = nn.ReLU(inplace=True)
+
+
+    def forward(self, x):
+        cond = self.cond_net(x)
+
+        scale1 = self.cond_scale1(cond)
+        shift1 = self.cond_shift1(cond)
+
+        scale2 = self.cond_scale2(cond)
+        shift2 = self.cond_shift2(cond)
+
+        scale3 = self.cond_scale3(cond)
+        shift3 = self.cond_shift3(cond)
+
+        out = self.conv1(x)
+        out = out * scale1.view(-1, self.base_nf, 1, 1) + shift1.view(-1, self.base_nf, 1, 1) + out
+        out = self.act(out)
+        
+
+        out = self.conv2(out)
+        out = out * scale2.view(-1, self.base_nf, 1, 1) + shift2.view(-1, self.base_nf, 1, 1) + out
+        out = self.act(out)
+
+        out = self.conv3(out)
+        out = out * scale3.view(-1, self.out_nc, 1, 1) + shift3.view(-1, self.out_nc, 1, 1) + out
+        return out

IIR-Lab/ISP_pipeline/raw_prc_pipeline/exif_data_formats.py

@@ -0,0 +1,22 @@
+class ExifFormat:
+    def __init__(self, id, name, size, short_name):
+        self.id = id
+        self.name = name
+        self.size = size
+        self.short_name = short_name  # used with struct.unpack()
+
+
+exif_formats = {
+    1: ExifFormat(1, 'unsigned byte', 1, 'B'),
+    2: ExifFormat(2, 'ascii string', 1, 's'),
+    3: ExifFormat(3, 'unsigned short', 2, 'H'),
+    4: ExifFormat(4, 'unsigned long', 4, 'L'),
+    5: ExifFormat(5, 'unsigned rational', 8, ''),
+    6: ExifFormat(6, 'signed byte', 1, 'b'),
+    7: ExifFormat(7, 'undefined', 1, 'B'),  # consider `undefined` as `unsigned byte`
+    8: ExifFormat(8, 'signed short', 2, 'h'),
+    9: ExifFormat(9, 'signed long', 4, 'l'),
+    10: ExifFormat(10, 'signed rational', 8, ''),
+    11: ExifFormat(11, 'single float', 4, 'f'),
+    12: ExifFormat(12, 'double float', 8, 'd'),
+}