IVL

IVL/Dockerfile

@@ -0,0 +1,12 @@
+# syntax=docker/dockerfile:1
+
+FROM python:3.10-slim-buster
+
+COPY requirements.txt .
+RUN pip install --no-cache -r requirements.txt
+
+RUN apt-get update
+RUN apt-get install ffmpeg libsm6 libxext6 libopenblas-dev -y
+
+COPY . /pipe24
+WORKDIR /pipe24

IVL/Grayness_Index.py

@@ -0,0 +1,118 @@
+import cv2
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+
+
+def DerivGauss(im, sigma=0.5):
+    gaussian_die_off = 0.000001
+    var = sigma ** 2
+    # compute filter width
+    width = None
+    for i in range(1, 51):
+        if np.exp(-(i ** 2) / (2 * var)) > gaussian_die_off:
+            width = i
+    if width is None:
+        width = 1
+
+    # create filter (derivative of Gaussian filter)
+    x = np.arange(-width, width + 1)
+    y = np.arange(-width, width + 1)
+    coordinates = np.meshgrid(x, y)
+    x = coordinates[0]
+    y = coordinates[1]
+    derivate_gaussian_2D = -x * \
+        np.exp(-(x * x + y * y) / (2 * var)) / (var * np.pi)
+
+    # apply filter and return magnitude
+    ax = cv2.filter2D(im, -1, derivate_gaussian_2D)
+    ay = cv2.filter2D(im, -1, np.transpose(derivate_gaussian_2D))
+    magnitude = np.sqrt((ax ** 2) + (ay ** 2))
+    return magnitude
+
+
+def GPconstancy_GI(im, gray_pixels, delta_th=10**(-4)):
+
+    # mask saturated pixels and mask very dark pixels
+    mask = np.logical_or(np.max(im, axis=2) >= 0.95,
+                         np.sum(im, axis=2) <= 0.0315)
+
+    # remove noise with mean filter
+    # mean_kernel = np.ones((7, 7), np.float32) / 7**2
+    # im = cv2.filter2D(im, -1, mean_kernel)
+
+    # decompose rgb values
+    r = im[:, :, 0]
+    g = im[:, :, 1]
+    b = im[:, :, 2]
+
+    # mask 0 elements
+    mask = np.logical_or.reduce((mask, r == 0, g == 0, b == 0))
+
+    # replace 0 values with machine epsilon
+    eps = np.finfo(np.float32).eps
+    r[r == 0] = eps
+    g[g == 0] = eps
+    b[b == 0] = eps
+    norm = r + g + b
+
+    # mask low contrast pixels
+    delta_r = DerivGauss(r)
+    delta_g = DerivGauss(g)
+    delta_b = DerivGauss(b)
+    mask = np.logical_or(mask, np.logical_and.reduce(
+        (delta_r <= delta_th, delta_g <= delta_th, delta_b <= delta_th)))
+
+    # compute colors in log domain, only red and blue
+    log_r = np.log(r) - np.log(norm)
+    log_b = np.log(b) - np.log(norm)
+
+    # mask low contrast pixels in the log domain
+    delta_log_r = DerivGauss(log_r)
+    delta_log_b = DerivGauss(log_b)
+    mask = np.logical_or.reduce(
+        (mask, delta_log_r == np.inf, delta_log_b == np.inf))
+
+    # normalize each channel in log domain
+    data = np.concatenate(
+        (np.reshape(delta_log_r, (-1, 1)), np.reshape(delta_log_b, (-1, 1))), axis=1)
+    mink_norm = 2
+    norm2_data = np.sum(data ** mink_norm, axis=1) ** (1 / mink_norm)
+    map_uniquelight = np.reshape(norm2_data, delta_log_r.shape)
+
+    # make masked pixels to max value
+    map_uniquelight[mask] = np.max(map_uniquelight)
+
+    # denoise
+    # map_uniquelight = cv2.filter2D(map_uniquelight, -1, mean_kernel)
+
+    # filter using map_uniquelight
+    gray_index_unique = map_uniquelight
+    sort_unique = np.sort(gray_index_unique.flatten())
+    gindex_unique = np.full(gray_index_unique.shape, False, dtype=bool)
+    gindex_unique[gray_index_unique <= sort_unique[gray_pixels - 1]] = True
+    choosen_pixels = im[gindex_unique]
+    mean = np.mean(choosen_pixels, axis=0)
+    result = mean / np.apply_along_axis(np.linalg.norm, 0, mean)
+    return result
+
+
+if __name__ == "__main__":
+
+    # read image and convert to 0 1
+    im_path = 'example.png'
+    im = cv2.cvtColor(cv2.imread(im_path), cv2.COLOR_BGR2RGB)
+    im = np.float32(im) / 255
+
+    tot_pixels = im.shape[0] * im.shape[1]
+    # compute number of gray pixels
+    n = 0.1  # 0.01%
+    num_gray_pixels = int(np.floor(n * tot_pixels / 100))
+
+    # compute global illuminant values
+    gt = np.array([0.6918, 0.6635, 0.2850])
+    lumTriplet = GPconstancy_GI(im, num_gray_pixels, 10**(-4))
+
+    # show results and angular error w.r.t gt
+    print(lumTriplet)
+    print(np.arccos(np.dot(lumTriplet, gt)) * 180 / np.pi)

IVL/blocks.py

@@ -0,0 +1,450 @@
+import cv2
+import numpy as np
+import skimage.color as cl
+from Grayness_Index import GPconstancy_GI
+from scipy.ndimage.filters import gaussian_filter
+from skimage import exposure
+from skimage.restoration import denoise_nl_means, estimate_sigma
+
+_RGB_TO_YCBCR = np.array([[0.257, 0.504, 0.098],
+                          [-0.148, -0.291, 0.439],
+                          [0.439, -0.368, -0.071]])
+
+_YCBCR_OFF = np.array([0.063, 0.502, 0.502])
+
+
+def _mul(coeffs, image):
+
+    r = image[:, :, 0]
+    g = image[:, :, 1]
+    b = image[:, :, 2]
+
+    r0 = np.repeat(r[:, :, np.newaxis], 3, 2) * coeffs[:, 0]
+    r1 = np.repeat(g[:, :, np.newaxis], 3, 2) * coeffs[:, 1]
+    r2 = np.repeat(b[:, :, np.newaxis], 3, 2) * coeffs[:, 2]
+
+    return r0 + r1 + r2
+
+
+def rgb2ycbcr(rgb):
+    """sRGB to YCbCr conversion."""
+    clip_rgb = False
+    if clip_rgb:
+        rgb = np.clip(rgb, 0, 1)
+    return _mul(_RGB_TO_YCBCR, rgb) + _YCBCR_OFF
+
+
+def ycbcr2rgb(rgb):
+    """YCbCr to sRGB conversion."""
+    clip_rgb = False
+    rgb = _mul(np.linalg.inv(_RGB_TO_YCBCR), rgb - _YCBCR_OFF)
+    if clip_rgb:
+        rgb = np.clip(rgb, 0, 1)
+    return rgb
+
+
+def normalize(raw_image, black_level, white_level):
+    if type(black_level) is list and len(black_level) == 1:
+        black_level = float(black_level[0])
+    if type(white_level) is list and len(white_level) == 1:
+        white_level = float(white_level[0])
+    black_level_mask = black_level
+    if type(black_level) is list and len(black_level) == 4:
+        if type(black_level[0]) is Ratio:
+            black_level = ratios2floats(black_level)
+        if type(black_level[0]) is Fraction:
+            black_level = fractions2floats(black_level)
+        black_level_mask = np.zeros(raw_image.shape)
+        idx2by2 = [[0, 0], [0, 1], [1, 0], [1, 1]]
+        step2 = 2
+        for i, idx in enumerate(idx2by2):
+            black_level_mask[idx[0]::step2, idx[1]::step2] = black_level[i]
+    normalized_image = raw_image.astype(np.float32) - black_level_mask
+    # if some values were smaller than black level
+    normalized_image[normalized_image < 0] = 0
+    normalized_image = normalized_image / (white_level - black_level_mask)
+    return normalized_image
+
+
+class LCC():
+
+    def __init__(self, sigma=None):
+        super(LCC, self).__init__()
+        if sigma is None:
+            sigma = np.sqrt(512 ** 2 + 512 ** 2) * 0.01
+        self.sigma = sigma
+
+    def __call__(self, image):
+        ycbcr = cl.rgb2ycbcr(image)
+        y = (ycbcr[:, :, 0] - 16) / 219
+        cb = ycbcr[:, :, 1]
+        cr = ycbcr[:, :, 2]
+
+        blurred_y = gaussian_filter(y, sigma=self.sigma)
+        mask = 1 - blurred_y
+
+        mean_intensity = np.mean(y)
+
+        alpha_lower = np.log(mean_intensity) / np.log(0.5)
+        alpha_upper = np.log(0.5) / np.log(mean_intensity)
+
+        condition = mean_intensity < 0.5
+        alpha = np.zeros(mask.shape)
+        alpha = np.where(condition, alpha_lower, alpha_upper)
+
+        gamma = alpha ** ((0.5 - mask) / 0.5)
+
+        new_y = y ** gamma
+
+        new_y = new_y * 219 + 16
+
+        new_ycbcr = np.stack([new_y, cb, cr], 2)
+
+        im_rgb = cl.ycbcr2rgb(new_ycbcr)
+        # im_rgb = np.clip(im_rgb, 0, 1)
+
+        im_out = contrast_saturation_fix(im_rgb, image)
+
+        return im_out
+
+
+def contrast_saturation_fix(enhanced_image, input_image, mode="LCC", n_bits=8):
+
+    im_ycbcr = rgb2ycbcr(enhanced_image)
+    or_ycbcr = rgb2ycbcr(input_image)
+
+    y_new = im_ycbcr[:, :, 0];
+    cb_new = im_ycbcr[:, :, 1];
+    cr_new = im_ycbcr[:, :, 2];
+
+    y = or_ycbcr[:, :, 0];
+    cb = or_ycbcr[:, :, 1];
+    cr = or_ycbcr[:, :, 2];
+
+    # dark pixels percentage
+    mask = np.logical_and(y < (35 / 255), (((cb - 0.5) * 2 +
+                                            (cr - 0.5) * 2) / 2) < (20 / 255))
+
+    dark_pixels = mask.flatten().sum()
+
+    if dark_pixels > 0:
+
+        ipixelCount, _ = np.histogram(y.flatten(), 256, range=(0, 1))
+        cdf = np.cumsum(ipixelCount)
+        idx = np.argmin(abs(cdf - (dark_pixels * 0.3)))
+        b_input30 = idx
+
+        ipixelCount, _ = np.histogram(y_new.flatten(), 256, range=(0, 1))
+        cdf = np.cumsum(ipixelCount)
+        idx = np.argmin(abs(cdf - (dark_pixels * 0.3)))
+        b_output30 = idx
+
+        bstr = (b_output30 - b_input30)
+    else:
+
+        bstr = np.floor(np.quantile(y_new.flatten(), 0.002) * 255)
+
+    if bstr > 50:
+        bstr = 50
+
+    dark_bound = bstr / 255
+
+    bright_b = np.floor(np.quantile(y_new.flatten(), 1 - 0.002) * 255)
+
+    if (255 - bright_b) > 50:
+        bright_b = 255 - 50
+
+    bright_bound = bright_b / 255
+
+    # y_new = (y_new - dark_bound) / (bright_bound - dark_bound)
+    y_new = exposure.rescale_intensity(y_new, in_range=(
+        y_new.min(), y_new.max()), out_range=(dark_bound, bright_bound))
+    y_new = y_new.clip(0, 1)
+
+    im_ycbcr[:, :, 0] = y_new
+    im_new = ycbcr2rgb(im_ycbcr)
+
+    im_new = im_new.clip(0, 1)
+
+    # Saturation
+
+    im_tmp = input_image
+
+    r = im_tmp[:, :, 0]
+    g = im_tmp[:, :, 1]
+    b = im_tmp[:, :, 2]
+
+    r_new = 0.5 * (((y_new / (y + 1e-40)) * (r + y)) + r - y)
+    g_new = 0.5 * (((y_new / (y + 1e-40)) * (g + y)) + g - y)
+    b_new = 0.5 * (((y_new / (y + 1e-40)) * (b + y)) + b - y)
+
+    im_new[:, :, 0] = r_new
+    im_new[:, :, 1] = g_new
+    im_new[:, :, 2] = b_new
+
+    return im_new
+
+
+def gamma_correction(img, exp):
+
+    return img ** exp
+
+
+def black_stretch(img, perc=0.2):
+
+    im_hsv = cl.rgb2hsv(img.clip(0, 1))
+    v = im_hsv[:, :, 2]
+
+    dark_bound = np.quantile(v.flatten(), perc, method='closest_observation')
+
+    v_new = (v - dark_bound) / (1 - dark_bound)
+
+    im_hsv[:, :, 2] = v_new.clip(0, 1)
+
+    out = cl.hsv2rgb(im_hsv)
+
+    return out.clip(0, 1)
+
+
+def saturation_scale(img, scale=2.):
+
+    img_hsv = cl.rgb2hsv(img.clip(0, 1))
+    s = img_hsv[:, :, 1]
+    s *= scale
+    img_hsv[:, :, 1] = s
+
+    return cl.hsv2rgb(np.clip(img_hsv, 0, 1))
+
+
+def global_mean_contrast(x, beta=0.5):
+
+    x_mean = np.mean(np.mean(x, 0), 0)
+    x_mean = np.expand_dims(np.expand_dims(x_mean, 0), 0)
+    x_mean = np.repeat(np.repeat(x_mean, x.shape[1], 1), x.shape[0], 0)
+
+    # scale all channels
+    out = x_mean + beta * (x - x_mean)
+
+    return out
+
+
+def sharpening(image, sigma=2.0, scale=1):
+
+    gaussian = cv2.GaussianBlur(image, (0, 0), sigma)
+
+    unsharp_image = image + scale * (image - gaussian)
+
+    return unsharp_image.clip(0, 1)
+
+
+def illumination_parameters_estimation(current_image, illumination_estimation_option):
+    ie_method = illumination_estimation_option.lower()
+    if ie_method == "gw":
+        ie = np.mean(current_image, axis=(0, 1))
+        ie /= ie[1]
+        return ie
+    elif ie_method == "sog":
+        sog_p = 4.
+        ie = np.mean(current_image**sog_p, axis=(0, 1))**(1 / sog_p)
+        ie /= ie[1]
+        return ie
+    elif ie_method == "wp":
+        ie = np.max(current_image, axis=(0, 1))
+        ie /= ie[1]
+        return ie
+    elif ie_method == "iwp":
+        samples_count = 20
+        sample_size = 20
+        rows, cols = current_image.shape[:2]
+        data = np.reshape(current_image, (rows * cols, 3))
+        maxima = np.zeros((samples_count, 3))
+        for i in range(samples_count):
+            maxima[i, :] = np.max(data[np.random.randint(
+                low=0, high=rows * cols, size=(sample_size)), :], axis=0)
+        ie = np.mean(maxima, axis=0)
+        ie /= ie[1]
+        return ie
+    else:
+        raise ValueError(
+            'Bad illumination_estimation_option value! Use the following options: "gw", "wp", "sog", "iwp"')
+
+
+def wb(demosaic_img, as_shot_neutral):
+
+    as_shot_neutral = np.asarray(as_shot_neutral)
+    # transform vector into matrix
+    if as_shot_neutral.shape == (3,):
+        as_shot_neutral = np.diag(1. / as_shot_neutral)
+
+    assert as_shot_neutral.shape == (3, 3)
+
+    white_balanced_image = np.dot(demosaic_img, as_shot_neutral.T)
+    white_balanced_image = np.clip(white_balanced_image, 0.0, 1.0)
+
+    return white_balanced_image
+
+
+def white_balance(img, n=0.1, th=1e-4, denoise_first=False):
+
+    uint = False
+    if np.issubdtype(img.dtype, np.uint8):
+        uint = True
+
+    if uint:
+        img = img.astype(np.float32) / 255
+
+    tot_pixels = img.shape[0] * img.shape[1]
+    # compute number of gray pixels
+    num_gray_pixels = int(np.floor(n * tot_pixels / 100))
+
+    # denoise if necessary
+    if denoise_first:
+        sigma_est = 1# np.mean(estimate_sigma(img, channel_axis=-1))
+        img_ = cv2.GaussianBlur(img,(0,0),5)
+    else:
+        img_ = img
+
+    # compute global illuminant values
+    lumTriplet = GPconstancy_GI(img_, num_gray_pixels, th)
+
+    lumTriplet /= lumTriplet.max()
+    out = wb(img, lumTriplet)
+
+    if uint:
+        return (out * 255).astype(np.uint8)
+    else:
+        return out
+
+
+def scurve(img, alpha=None, lmbd=1 / 1.8, blacks=False):
+
+    x = img
+
+    if alpha is None:
+        im_hsv = cl.rgb2hsv(img.clip(0, 1))
+        v = im_hsv[:, :, 2]
+
+        alpha = np.quantile(v.flatten(), 0.02, method='closest_observation')
+
+    if not blacks:
+        out = np.where(x <= alpha,
+                       x,  # alpha - alpha * (1 - x / alpha) ** lmbd,
+                       alpha + (1 - alpha) *
+                       ((x - alpha) / (1 - alpha)).clip(min=0) ** lmbd
+                       )
+    else:
+        out = np.where(x <= alpha,
+                       alpha - alpha * (1 - x / alpha) ** lmbd,
+                       x
+                       )
+
+    # out = out.clip(0, 1)
+
+    return out
+
+
+def scurve_central(img, lmbd=1 / 1.4, blacks=False):
+
+    x = img
+
+    im_hsv = cl.rgb2hsv(img.clip(0, 1))
+    v = im_hsv[:, :, 2]
+
+    alpha1 = np.quantile(v.flatten(), 0.2, method='closest_observation')
+    alpha2 = np.quantile(v.flatten(), 0.9, method='closest_observation')
+
+    out = np.where(x <= alpha1,
+                   x,
+                   np.where(x >= alpha2,
+                            x,
+                            alpha1 + (alpha2 - alpha1) *
+                            ((x - alpha1) / (alpha2 - alpha1)).clip(min=0) ** lmbd
+                            )
+                   )
+
+    return out
+
+
+def imadjust(img, hi=0.9999, pi=0.0001):
+    '''
+    Python version of matlab imadjust
+    '''
+
+    im_hsv = cl.rgb2hsv(img.clip(0, 1))
+    v = im_hsv[:, :, 2]
+
+    hi = np.quantile(v.flatten(), hi, method='closest_observation')
+    li = np.quantile(v.flatten(), pi, method='closest_observation')
+
+    if hi < 0.7:
+        hi = np.quantile(v.flatten(), 0.995, method='closest_observation')
+
+    if hi == 1:
+        v_tmp = v.flatten()
+        v_tmp = v_tmp[v_tmp != 1]
+        hi = np.quantile(v_tmp, 0.9995, method='closest_observation')
+    if li == 0:
+        v_tmp = v.flatten()
+        v_tmp = v_tmp[v_tmp != 0]
+        li = np.quantile(v_tmp, 0.0001, method='closest_observation')
+
+    x = img
+    li = li
+    hi = hi
+
+    lo = 0
+    ho = 0.9
+    gamma = 1
+
+    out = ((x - li) / (hi - li)) ** gamma
+    out = out * (ho - lo) + lo
+
+    return out
+
+
+def denoise_raw(image, l_w=3, ch_w=20):
+
+    im_yuv = cl.rgb2yuv(image)
+
+    # Separately process luma and choma
+
+    patch_kw = dict(patch_size=5,
+                    patch_distance=6
+                    )
+    sigma_est = np.mean(estimate_sigma(im_yuv[:, :, 0]))
+
+    den_y = denoise_nl_means(im_yuv[:, :, 0], h=l_w * sigma_est, fast_mode=True,
+                             **patch_kw)
+
+    patch_kw = dict(patch_size=5,
+                    patch_distance=6,
+                    channel_axis=-1
+                    )
+    sigma_est = np.mean(estimate_sigma(im_yuv[:, :, 1:2], channel_axis=-1))
+    den_uv = denoise_nl_means(im_yuv[:, :, 1:3], h=ch_w * sigma_est, fast_mode=True,
+                              **patch_kw)
+
+    out = im_yuv
+
+    out[:, :, 0] = den_y
+    out[:, :, 1:3] = den_uv
+
+    del den_y
+    del den_uv
+
+    out = cl.yuv2rgb(out)
+
+    return out
+
+def denoise_rgb(image, l_w=3, ch_w=None):
+
+    #patch_kw = dict(patch_size=5,
+    #                patch_distance=6
+    #               )
+    patch_kw = {}
+    sigma_est = np.mean(estimate_sigma(image, channel_axis=2))
+    out = denoise_nl_means(image, h=l_w * sigma_est, fast_mode=True,
+                             **patch_kw, channel_axis=2)
+
+
+    return out

IVL/pipeline24.py

@@ -0,0 +1,167 @@
+###################################################
+# Night challenge 2024
+###################################################
+
+
+import argparse
+import threading
+from pathlib import Path
+import os
+import queue
+from time import time
+from tqdm import tqdm
+
+import blocks as B
+import cv2
+import numpy as np
+import skimage.color as cl
+from skimage.transform import resize
+from raw_prc_pipeline import (expected_img_ext, expected_landscape_img_height,
+                              expected_landscape_img_width)
+from raw_prc_pipeline.pipeline import RawProcessingPipelineDemo
+from tqdm import tqdm
+from utils import fraction_from_json, json_read
+
+
+def parse_args():
+    parser = argparse.ArgumentParser(
+        description='Script for processing PNG images with given metadata files.')
+    # folder params
+    parser.add_argument('-p', '--png_dir', type=Path, required=True,
+                        help='Path of the directory containing PNG images with metadata files.')
+    parser.add_argument('-o', '--out_dir', type=Path, default='./results',
+                        help='Path to the directory where processed images will be saved. Images will be saved in JPG format.')
+    # raw processing params
+    parser.add_argument('-ie', '--illumination_estimation', type=str, default='',
+                        help='Options for illumination estimation algorithms: "gw", "wp", "sog", "iwp".')
+    parser.add_argument('-tm', '--tone_mapping', type=str, default='Flash',
+                        help='Options for tone mapping algorithms: "Base", "Flash", "Storm", "Linear", "Drago", "Mantiuk", "Reinhard".')
+    # srgb processing params
+    parser.add_argument('-gc', '--gamma_correction', type=float, default=1 / 1.4,
+                        help='Global gamma correction.')
+    parser.add_argument('-dm', '--denoise_mask', type=float, default=0.6,
+                        help='Value to control denoising effect in bright regions. Should be between 0 and 1')
+    args = parser.parse_args()
+
+    if args.out_dir is None:
+        args.out_dir = args.png_dir
+
+    return args
+
+
+class PNGProcessing():
+    def __init__(self, ie_method, tone_mapping, gamma_correction, denoise_mask):
+        self.pipeline_params = {
+            'illumination_estimation': ie_method,
+            # 'tone_mapping': tone_mapping,
+            'out_landscape_width': expected_landscape_img_width,
+            'out_landscape_height': expected_landscape_img_height
+        }
+
+        self.pipeline = RawProcessingPipelineDemo(**self.pipeline_params)
+        self.gamma_correction = gamma_correction
+        self.denoise_mask = denoise_mask
+
+    def pipeline_exec(self, raw_image, metadata):
+
+        normalized_image = self.pipeline.normalize(raw_image, metadata)
+
+        demosaiced_image = self.pipeline.demosaic(normalized_image, metadata)
+        # check the original demosaicing to see if results are the same
+        
+        demosaiced_image = resize(demosaiced_image, (768, 1024), preserve_range=True, anti_aliasing=True)
+
+        wb_image = self.pipeline.white_balance(demosaiced_image, metadata)
+
+        xyz_image = self.pipeline.xyz_transform(wb_image, metadata)
+        srgb_image = self.pipeline.srgb_transform(xyz_image, metadata)
+
+        denoised_image = B.denoise_raw(
+            srgb_image, l_w=1, ch_w=7)
+            # srgb_image, l_w=4.5, ch_w=20)
+            # srgb_image, l_w=1.659923974475318, ch_w=5.459274910995606)
+
+        light_enhancer = B.LCC(2)
+        # light_enhancer = B.LCC(sigma=6.463076463115174)
+        light_image = light_enhancer(denoised_image).clip(0)
+
+        contrast_image = B.global_mean_contrast(light_image, beta=1.5).clip(0)
+        # contrast_image = B.global_mean_contrast(light_image, beta=0.8653634653721171).clip(0)
+
+        gamma_image = B.scurve(contrast_image, alpha=0, lmbd=(1 / 1.8)).clip(0)
+        # gamma_image = B.scurve(contrast_image, alpha=0.7050463096367395, lmbd=0.9740931227248038).clip(0)
+
+        black_adj_image = B.imadjust(gamma_image, 0.99).clip(0)
+        # black_adj_image = B.imadjust(gamma_image, 0.9957007298972433, 0.01697803128505186).clip(0)
+
+        im_h = cl.rgb2hsv(black_adj_image)[:, :, 2]
+        if im_h.mean() < 0.2:
+            black_adj_image = B.scurve_central(black_adj_image, lmbd=(1 / 1.8)).clip(0)
+            # black_adj_image = B.scurve_central(black_adj_image, lmbd=0.6913136563678325).clip(0)
+        elif im_h.mean() < 0.25:
+            black_adj_image = B.scurve_central(black_adj_image, lmbd=(1 / 1.4)).clip(0)
+            # black_adj_image = B.scurve_central(black_adj_image, lmbd=0.3612134419536918).clip(0)
+        elif im_h.mean() > 0.4:
+            black_adj_image = B.gamma_correction(black_adj_image, 1.6).clip(0)
+            # black_adj_image = B.gamma_correction(black_adj_image, 3.5208650132731356).clip(0)
+
+        sharp_image = B.sharpening(black_adj_image, sigma=1).clip(0)
+        # sharp_image = B.sharpening(black_adj_image, 1.5389081796026578, 0.05456721376794549).clip(0)
+
+        wb_image = B.white_balance(sharp_image, denoise_first=True).clip(0)
+        # wb_image = B.white_balance(sharp_image, 0.740831363817609, 0.004044358054560114).clip(0)
+        
+        uint8_image = self.pipeline.to_uint8(wb_image, metadata)
+        # resized_image = self.pipeline.resize(uint8_image, metadata)
+        resulted_image = self.pipeline.fix_orientation(uint8_image, metadata)
+
+
+        return resulted_image
+
+    def __call__(self, png_path: Path, out_path: Path):
+
+        # parse raw img
+        raw_image = cv2.imread(str(png_path), cv2.IMREAD_UNCHANGED)
+        # parse metadata
+        metadata = json_read(png_path.with_suffix(
+            '.json'), object_hook=fraction_from_json)
+
+        start = time()
+        output_image = self.pipeline_exec(raw_image, metadata)
+        end = time()
+        
+        # save results
+        output_image = cv2.cvtColor(output_image, cv2.COLOR_RGB2BGR)
+        cv2.imwrite(str(out_path), output_image, [
+                    cv2.IMWRITE_JPEG_QUALITY, 100])
+        return end - start
+
+
+def process(png_processor, out_dir, png_paths):
+    out_paths = [
+        out_dir / png_path.with_suffix(expected_img_ext).name for png_path in png_paths]
+    times = []
+    pbar = tqdm(total=len(png_paths), ncols=100)
+    for png_path, out_path in zip(png_paths, out_paths):
+        runtime = png_processor(png_path, out_path)
+        times.append(runtime)
+        pbar.update()
+    return times
+    
+
+def main(png_dir, out_dir, illumination_estimation, tone_mapping, gamma_correction, denoise_mask):
+    # out_dir.mkdir(exist_ok=True)
+    os.makedirs(out_dir, exist_ok=True)
+
+    png_paths = list(png_dir.glob('*.png'))
+
+    png_processor = PNGProcessing(
+        illumination_estimation, tone_mapping, gamma_correction, denoise_mask)
+
+    times = process(png_processor, out_dir, png_paths)
+    print(f'Average time: {np.mean(times)} seconds.')
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    main(**vars((args)))

IVL/raw_prc_pipeline/__init__.py

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

IVL/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'),
+}

IVL/raw_prc_pipeline/exif_utils.py

@@ -0,0 +1,208 @@
+"""
+Manual parsing of image file directories (IFDs).
+"""
+
+
+import struct
+from fractions import Fraction
+from raw_prc_pipeline.exif_data_formats import exif_formats
+
+class Ifd:
+    def __init__(self):
+        self.offset = -1
+        self.tags = {}  # <key, tag> dict; tag number will be key.
+
+
+class Tag:
+    def __init__(self):
+        self.offset = -1
+        self.tag_num = -1
+        self.data_format = -1
+        self.num_values = -1
+        self.values = []
+
+
+def parse_exif(image_path, verbose=True):
+    """
+    Parse EXIF tags from a binary file and return IFDs.
+    Returned IFDs include EXIF SubIFDs, if any.
+    """
+
+    def print_(str_):
+        if verbose:
+            print(str_)
+
+    ifds = {}  # dict of <offset, Ifd> pairs; using offset to IFD as key.
+
+    with open(image_path, 'rb') as fid:
+        fid.seek(0)
+        b0 = fid.read(1)
+        _ = fid.read(1)
+        # byte storage direction (endian):
+        # +1: b'M' (big-endian/Motorola)
+        # -1: b'I' (little-endian/Intel)
+        endian = 1 if b0 == b'M' else -1
+        print_("Endian = {}".format(b0))
+        endian_sign = "<" if endian == -1 else ">"  # used in struct.unpack
+        print_("Endian sign = {}".format(endian_sign))
+        _ = fid.read(2)  # 0x002A
+        b4_7 = fid.read(4)  # offset to first IFD
+        offset_ = struct.unpack(endian_sign + "I", b4_7)[0]
+        i = 0
+        ifd_offsets = [offset_]
+        while len(ifd_offsets) > 0:
+            offset_ = ifd_offsets.pop(0)
+            # check if IFD at this offset was already parsed before
+            if offset_ in ifds:
+                continue
+            print_("=========== Parsing IFD # {} ===========".format(i))
+            ifd_ = parse_exif_ifd(fid, offset_, endian_sign, verbose)
+            ifds.update({ifd_.offset: ifd_})
+            print_("=========== Finished parsing IFD # {} ===========".format(i))
+            i += 1
+            # check SubIFDs; zero or more offsets at tag 0x014a
+            sub_idfs_tag_num = int('0x014a', 16)
+            if sub_idfs_tag_num in ifd_.tags:
+                ifd_offsets.extend(ifd_.tags[sub_idfs_tag_num].values)
+            # check Exif SUbIDF; usually one offset at tag 0x8769
+            exif_sub_idf_tag_num = int('0x8769', 16)
+            if exif_sub_idf_tag_num in ifd_.tags:
+                ifd_offsets.extend(ifd_.tags[exif_sub_idf_tag_num].values)
+    return ifds
+
+
+def parse_exif_ifd(binary_file, offset_, endian_sign, verbose=True):
+    """
+    Parse an EXIF IFD.
+    """
+
+    def print_(str_):
+        if verbose:
+            print(str_)
+
+    ifd = Ifd()
+    ifd.offset = offset_
+    print_("IFD offset = {}".format(ifd.offset))
+    binary_file.seek(offset_)
+    num_entries = struct.unpack(endian_sign + "H", binary_file.read(2))[0]  # format H = unsigned short
+    print_("Number of entries = {}".format(num_entries))
+    for t in range(num_entries):
+        print_("---------- Tag {} / {} ----------".format(t + 1, num_entries))
+        if t == 22:
+            ttt = 1
+        tag_ = parse_exif_tag(binary_file, endian_sign, verbose)
+        ifd.tags.update({tag_.tag_num: tag_})  # supposedly, EXIF tag numbers won't repeat in the same IFD
+    # TODO: check for subsequent IFDs by parsing the next 4 bytes immediately after the IFD
+    return ifd
+
+
+def parse_exif_tag(binary_file, endian_sign, verbose=True):
+    """
+    Parse EXIF tag from a binary file starting from the current file pointer and returns the tag values.
+    """
+
+    def print_(str_):
+        if verbose:
+            print(str_)
+
+    tag = Tag()
+
+    # tag offset
+    tag.offset = binary_file.tell()
+    print_("Tag offset = {}".format(tag.offset))
+
+    # tag number
+    bytes_ = binary_file.read(2)
+    tag.tag_num = struct.unpack(endian_sign + "H", bytes_)[0]  # H: unsigned 2-byte short
+    print_("Tag number = {} = 0x{:04x}".format(tag.tag_num, tag.tag_num))
+
+    # data format (some value between [1, 12])
+    tag.data_format = struct.unpack(endian_sign + "H", binary_file.read(2))[0]  # H: unsigned 2-byte short
+    exif_format = exif_formats[tag.data_format]
+    print_("Data format = {} = {}".format(tag.data_format, exif_format.name))
+
+    # number of components/values
+    tag.num_values = struct.unpack(endian_sign + "I", binary_file.read(4))[0]  # I: unsigned 4-byte integer
+    print_("Number of values = {}".format(tag.num_values))
+
+    # total number of data bytes
+    total_bytes = tag.num_values * exif_format.size
+    print_("Total bytes = {}".format(total_bytes))
+
+    # seek to data offset (if needed)
+    data_is_offset = False
+    current_offset = binary_file.tell()
+    if total_bytes > 4:
+        print_("Total bytes > 4; The next 4 bytes are an offset.")
+        data_is_offset = True
+        data_offset = struct.unpack(endian_sign + "I", binary_file.read(4))[0]
+        current_offset = binary_file.tell()
+        print_("Current offset = {}".format(current_offset))
+        print_("Seeking to data offset = {}".format(data_offset))
+        binary_file.seek(data_offset)
+
+    # read values
+    # TODO: need to distinguish between numeric and text values?
+    if tag.num_values == 1 and total_bytes < 4:
+        # special case: data is a single value that is less than 4 bytes inside 4 bytes, take care of endian
+        val_bytes = binary_file.read(4)
+        # if endian_sign == ">":
+        # val_bytes = val_bytes[4 - total_bytes:]
+        # else:
+        # val_bytes = val_bytes[:total_bytes][::-1]
+        val_bytes = val_bytes[:total_bytes]
+        tag.values.append(struct.unpack(endian_sign + exif_format.short_name, val_bytes)[0])
+    else:
+        # read data values one by one
+        for k in range(tag.num_values):
+            val_bytes = binary_file.read(exif_format.size)
+            if exif_format.name == 'unsigned rational':
+                tag.values.append(eight_bytes_to_fraction(val_bytes, endian_sign, signed=False))
+            elif exif_format.name == 'signed rational':
+                tag.values.append(eight_bytes_to_fraction(val_bytes, endian_sign, signed=True))
+            else:
+                tag.values.append(struct.unpack(endian_sign + exif_format.short_name, val_bytes)[0])
+        if total_bytes < 4:
+            # special case: multiple values less than 4 bytes in total, inside the 4 bytes; skip the extra bytes
+            binary_file.seek(4 - total_bytes, 1)
+
+    if verbose:
+        if len(tag.values) > 100:
+            print_("Got more than 100 values; printing first 100 only:")
+            print_("Tag values = {}".format(tag.values[:100]))
+        else:
+            print_("Tag values = {}".format(tag.values))
+    if tag.data_format == 2:
+        print_("Tag values (string) = {}".format(b''.join(tag.values).decode()))
+
+    if data_is_offset:
+        # seek back to current position to read the next tag
+        print_("Seeking back to current offset = {}".format(current_offset))
+        binary_file.seek(current_offset)
+
+    return tag
+
+
+def get_tag_values_from_ifds(tag_num, ifds):
+    """
+    Return values of a tag, if found in ifds. Return None otherwise.
+    Assuming any tag exists only once in all ifds.
+    """
+    for key, ifd in ifds.items():
+        if tag_num in ifd.tags:
+            return ifd.tags[tag_num].values
+    return None
+
+
+def eight_bytes_to_fraction(eight_bytes, endian_sign, signed):
+    """
+    Convert 8-byte array into a Fraction. Take care of endian and sign.
+    """
+    if signed:
+        num = struct.unpack(endian_sign + "l", eight_bytes[:4])[0]
+        den = struct.unpack(endian_sign + "l", eight_bytes[4:])[0]
+    else:
+        num = struct.unpack(endian_sign + "L", eight_bytes[:4])[0]
+        den = struct.unpack(endian_sign + "L", eight_bytes[4:])[0]
+    den = den if den != 0 else 1
+    return Fraction(num, den)

IVL/raw_prc_pipeline/fs.py

@@ -0,0 +1,43 @@
+import cv2
+import numpy as np
+
+
+def perform_flash(source, a=5, target=-1, perform_gamma_correction=True):
+    rows, cols, _ = source.shape
+
+    v = np.max(source, axis=2)
+    vd = np.copy(v)
+    vd[vd == 0] = 1e-9
+    result = source / (a * np.exp(np.mean(np.log(vd))) + np.tile(np.expand_dims(vd, axis=2), (1, 1, 3)))
+
+    if perform_gamma_correction:
+        result **= 1.0 / 2.2
+
+    if target >= 0:
+        result *= target / np.mean((0.299 * result[:, :, 2] + 0.587 * result[:, :, 1] + 0.114 * result[:, :, 0]))
+    else:
+        result *= 255.0 / np.max(result)
+
+    return result
+
+
+def perform_storm(source, a=5, target=-1, kernels=(1, 4, 16, 64, 256), perform_gamma_correction=True):
+    rows, cols, _ = source.shape
+
+    v = np.max(source, axis=2)
+    vd = np.copy(v)
+    vd[vd == 0] = 1e-9
+    lv = np.log(vd)
+    result = sum([source / np.tile(
+        np.expand_dims(a * np.exp(cv2.boxFilter(lv, -1, (int(min(rows // kernel, cols // kernel)),) * 2)) + vd, axis=2),
+        (1, 1, 3)) for kernel in kernels])
+
+    if perform_gamma_correction:
+        result **= 1.0 / 2.2
+
+    if target >= 0:
+        result *= target / np.mean((0.299 * result[:, :, 2] + 0.587 * result[:, :, 1] + 0.114 * result[:, :, 0]))
+    else:
+        result *= 255.0 / np.max(result)
+
+    return result

IVL/raw_prc_pipeline/pipeline.py

@@ -0,0 +1,211 @@
+"""
+Demo raw processing pipeline and pipeline executor.
+"""
+
+import numpy as np
+from raw_prc_pipeline.pipeline_utils import *
+
+
+class RawProcessingPipelineDemo:
+    """
+    Demonstration pipeline of raw image processing.
+
+    This pipeline is a baseline pipeline to process raw image.
+    The public methods of this class are successive steps of raw image processing pipeline.
+    The declaration order of the public methods must correspond to the order in which these methods (steps) are supposed to be called when processing raw image.
+
+    It is assumed that each public method has 2 parameters:
+    raw_img : ndarray
+        Array with images data.
+    img_meta : Dict
+        Some metadata of image.
+    
+    Also each such public method must return an image (ndarray) as the result of processing.
+    """
+    def __init__(self, illumination_estimation='', denoise_flg=True, tone_mapping='Flash', out_landscape_width=None, out_landscape_height=None):
+        """
+        RawProcessingPipelineDemo __init__ method.
+
+        Parameters
+        ----------
+        illumination_estimation : str, optional
+            Options for illumination estimation algorithms: '', 'gw', 'wp', 'sog', 'iwp', by default ''.
+        denoise_flg : bool, optional
+            Denoising flag, by default True.
+            If True, resulted images will be denoised with some predefined parameters.
+        tone_mapping : str, optional
+            Options for tone mapping methods, defined in function `apply_tone_map` from `pipeline_utils` module.
+            By default 'Flash'.
+        out_landscape_width : int, optional
+            The width of output image (when orientation is landscape). If None, the image resize will not be performed.
+            By default None.
+        out_landscape_height : int, optional
+            The height of output image (when orientation is landscape). If None, the image resize will not be performed.
+            By default None.
+        """
+        self.params = locals()
+        del self.params['self']
+
+    # Linearization not handled.
+    def linearize_raw(self, raw_img, img_meta):
+        return raw_img
+
+    def normalize(self, linearized_raw, img_meta):
+        return normalize(linearized_raw, img_meta['black_level'], img_meta['white_level'])
+
+    def demosaic(self, normalized, img_meta):
+        return simple_demosaic(normalized, img_meta['cfa_pattern'])
+
+    def denoise(self, demosaic, img_meta):
+        if not self.params['denoise_flg']:
+            return demosaic
+        return denoise_image(demosaic)
+
+    def white_balance(self, demosaic, img_meta):
+        if self.params['illumination_estimation'] == '':
+            wb_params = img_meta['as_shot_neutral']
+        else:
+            wb_params = illumination_parameters_estimation(
+                demosaic, self.params['illumination_estimation'])
+
+        white_balanced = white_balance(demosaic, wb_params)
+        return white_balanced
+
+    def xyz_transform(self, white_balanced, img_meta):
+        # in case of absence of color matrix we use mean color matrix
+        if "color_matrix_1" not in img_meta.keys():
+            ccm_default = [1.06835938, -0.29882812, -0.14257812, 
+                           -0.43164062,  1.35546875,  0.05078125, 
+                           -0.1015625, 0.24414062, 0.5859375]
+            img_meta["color_matrix_1"] = ccm_default
+            img_meta["color_matrix_2"] = ccm_default
+        return apply_color_space_transform(white_balanced, img_meta['color_matrix_1'], img_meta['color_matrix_2'])
+
+    def srgb_transform(self, xyz, img_meta):
+        return transform_xyz_to_srgb(xyz)
+
+    def tone_mapping(self, srgb, img_meta):
+        if self.params['tone_mapping'] is None:
+            return apply_tone_map(srgb, 'Base')
+        return apply_tone_map(srgb, self.params['tone_mapping'])
+
+    def gamma_correct(self, srgb, img_meta):
+        return apply_gamma(srgb)
+
+    def autocontrast(self, srgb, img_meta):
+        # return autocontrast(srgb)
+        return autocontrast_using_pil(srgb)
+
+    def to_uint8(self, srgb, img_meta):
+        return (srgb*255).astype(np.uint8)
+    
+    def resize(self, img, img_meta):
+        if self.params['out_landscape_width'] is None or self.params['out_landscape_height'] is None:
+            return img
+        return resize_using_pil(img, self.params['out_landscape_width'], self.params['out_landscape_height'])
+    
+    def fix_orientation(self, img, img_meta):
+        return fix_orientation(img, img_meta['orientation'])
+
+
+class PipelineExecutor:
+    """
+    Pipeline executor class.
+
+    This class can be used to successively execute the steps of some image pipeline class (for example `RawProcessingPipelineDemo`).
+    The declaration order of the public methods of pipeline class must correspond to the order in which these methods (steps) are supposed to be called when processing image.
+
+    It is assumed that each public method of the pipeline class has 2 parameters:
+    raw_img : ndarray
+        Array with images data.
+    img_meta : Dict
+        Some meta data of image.
+    
+    Also each such public method must return an image (ndarray) as the result of processing.
+    """
+    def __init__(self, img, img_meta, pipeline_obj, first_stage=None, last_stage=None):
+        """
+        PipelineExecutor __init__ method.
+
+        Parameters
+        ----------
+        img : ndarray
+            Image that should be processed by pipeline.
+        img_meta : Dict
+            Some image metadata.
+        pipeline_obj : pipeline object
+            Some pipeline object such as RawProcessingPipelineDemo.
+        first_stage : str, optional
+            The name of first public method of pipeline object that should be called by PipelineExecutor.
+            If None, the first public method from defined in pipeline object will be considered as `first_stage` method.
+            By default None.
+        last_stage : str, optional
+            The name of last public method of pipeline object that should be called by PipelineExecutor.
+            If None, the last public method from defined in pipeline object will be considered as `last_stage` method.
+            By default None.
+        """
+        self.pipeline_obj = pipeline_obj
+        self.stages_dict = self._init_stages()
+        self.stages_names, self.stages = list(
+            self.stages_dict.keys()), list(self.stages_dict.values())
+
+        if first_stage is None:
+            self.next_stage_indx = 0
+        else:
+            assert first_stage in self.stages_names, f"Invalid first_stage={first_stage}. Try use the following stages: {self.stages_names}"
+            self.next_stage_indx = self.stages_names.index(first_stage)
+
+        if last_stage is None:
+            self.last_stage_indx = len(self.stages_names) - 1
+        else:
+            assert last_stage in self.stages_names, f"Invalid last_stage={last_stage}. Try use the following stages: {self.stages_names}"
+            self.last_stage_indx = self.stages_names.index(last_stage)
+            if self.next_stage_indx > self.last_stage_indx:
+                print(f'Warning: the specified first_stage={first_stage} follows the specified last_stage={last_stage}, so using __call__ no image processing will be done.')
+
+        self.current_image = img
+        self.img_meta = img_meta
+
+    def _init_stages(self):
+        stages = {func: getattr(self.pipeline_obj, func) for func in self.pipeline_obj.__class__.__dict__ if callable(
+            getattr(self.pipeline_obj, func)) and not func.startswith("_")}
+        return stages
+
+    @property
+    def next_stage(self):
+        if self.next_stage_indx < len(self.stages):
+            return self.stages_names[self.next_stage_indx]
+        else:
+            return None
+
+    @property
+    def last_stage(self):
+        return self.stages_names[self.last_stage_indx]
+
+    def __iter__(self):
+        return self
+
+    def __next__(self):
+        if self.next_stage_indx < len(self.stages):
+            stage_func = self.stages[self.next_stage_indx]
+            self.current_image = stage_func(self.current_image, self.img_meta)
+            self.next_stage_indx += 1
+            return self.current_image
+        else:
+            raise StopIteration
+
+    def __call__(self):
+        """
+        PipelineExecutor __call__ method.
+
+        This method will sequentially execute the methods defined in the pipeline object from the `first_stage` to the `last_stage` inclusive.
+
+        Returns
+        -------
+        ndarray
+            Resulted processed raw image.
+        """
+        for current_image in self:
+            if self.next_stage_indx > self.last_stage_indx:
+                return current_image
+        return self.current_image

IVL/raw_prc_pipeline/pipeline_utils.py

@@ -0,0 +1,493 @@
+"""
+Camera pipeline utilities.
+"""
+
+import os
+from fractions import Fraction
+
+import cv2
+import numpy as np
+import exifread
+# from exifread import Ratio
+from exifread.utils import Ratio
+import rawpy
+from scipy.io import loadmat
+from raw_prc_pipeline.exif_utils import parse_exif, get_tag_values_from_ifds
+from raw_prc_pipeline.fs import perform_storm, perform_flash
+from PIL import Image, ImageOps
+from skimage.restoration import denoise_bilateral
+from skimage.transform import resize as skimage_resize
+
+
+def get_visible_raw_image(image_path):
+    raw_image = rawpy.imread(image_path).raw_image_visible.copy()
+    # raw_image = rawpy.imread(image_path).raw_image.copy()
+    return raw_image
+
+
+def get_image_tags(image_path):
+    with open(image_path, 'rb') as f:
+        tags = exifread.process_file(f)
+    return tags
+
+
+def get_image_ifds(image_path):
+    ifds = parse_exif(image_path, verbose=False)
+    return ifds
+
+
+def get_metadata(image_path):
+    metadata = {}
+    tags = get_image_tags(image_path)
+    ifds = get_image_ifds(image_path)
+    metadata['linearization_table'] = get_linearization_table(tags, ifds)
+    metadata['black_level'] = get_black_level(tags, ifds)
+    metadata['white_level'] = get_white_level(tags, ifds)
+    metadata['cfa_pattern'] = get_cfa_pattern(tags, ifds)
+    metadata['as_shot_neutral'] = get_as_shot_neutral(tags, ifds)
+    color_matrix_1, color_matrix_2 = get_color_matrices(tags, ifds)
+    metadata['color_matrix_1'] = color_matrix_1
+    metadata['color_matrix_2'] = color_matrix_2
+    metadata['orientation'] = get_orientation(tags, ifds)
+    # isn't used
+    metadata['noise_profile'] = get_noise_profile(tags, ifds)
+    # ...
+    # fall back to default values, if necessary
+    if metadata['black_level'] is None:
+        metadata['black_level'] = 0
+        print("Black level is None; using 0.")
+    if metadata['white_level'] is None:
+        metadata['white_level'] = 2 ** 16
+        print("White level is None; using 2 ** 16.")
+    if metadata['cfa_pattern'] is None:
+        metadata['cfa_pattern'] = [0, 1, 1, 2]
+        print("CFAPattern is None; using [0, 1, 1, 2] (RGGB)")
+    if metadata['as_shot_neutral'] is None:
+        metadata['as_shot_neutral'] = [1, 1, 1]
+        print("AsShotNeutral is None; using [1, 1, 1]")
+    if metadata['color_matrix_1'] is None:
+        metadata['color_matrix_1'] = [1] * 9
+        print("ColorMatrix1 is None; using [1, 1, 1, 1, 1, 1, 1, 1, 1]")
+    if metadata['color_matrix_2'] is None:
+        metadata['color_matrix_2'] = [1] * 9
+        print("ColorMatrix2 is None; using [1, 1, 1, 1, 1, 1, 1, 1, 1]")
+    if metadata['orientation'] is None:
+        metadata['orientation'] = 0
+        print("Orientation is None; using 0.")
+    # ...
+    return metadata
+
+
+def get_linearization_table(tags, ifds):
+    possible_keys = ['Image Tag 0xC618', 'Image Tag 50712',
+                     'LinearizationTable', 'Image LinearizationTable']
+    return get_values(tags, possible_keys)
+
+
+def get_black_level(tags, ifds):
+    possible_keys = ['Image Tag 0xC61A', 'Image Tag 50714',
+                     'BlackLevel', 'Image BlackLevel']
+    vals = get_values(tags, possible_keys)
+    if vals is None:
+        # print("Black level not found in exifread tags. Searching IFDs.")
+        vals = get_tag_values_from_ifds(50714, ifds)
+    return vals
+
+
+def get_white_level(tags, ifds):
+    possible_keys = ['Image Tag 0xC61D', 'Image Tag 50717',
+                     'WhiteLevel', 'Image WhiteLevel']
+    vals = get_values(tags, possible_keys)
+    if vals is None:
+        # print("White level not found in exifread tags. Searching IFDs.")
+        vals = get_tag_values_from_ifds(50717, ifds)
+    return vals
+
+
+def get_cfa_pattern(tags, ifds):
+    possible_keys = ['CFAPattern', 'Image CFAPattern']
+    vals = get_values(tags, possible_keys)
+    if vals is None:
+        # print("CFAPattern not found in exifread tags. Searching IFDs.")
+        vals = get_tag_values_from_ifds(33422, ifds)
+    return vals
+
+
+def get_as_shot_neutral(tags, ifds):
+    possible_keys = ['Image Tag 0xC628', 'Image Tag 50728',
+                     'AsShotNeutral', 'Image AsShotNeutral']
+    return get_values(tags, possible_keys)
+
+
+def get_color_matrices(tags, ifds):
+    possible_keys_1 = ['Image Tag 0xC621', 'Image Tag 50721',
+                       'ColorMatrix1', 'Image ColorMatrix1']
+    color_matrix_1 = get_values(tags, possible_keys_1)
+    possible_keys_2 = ['Image Tag 0xC622', 'Image Tag 50722',
+                       'ColorMatrix2', 'Image ColorMatrix2']
+    color_matrix_2 = get_values(tags, possible_keys_2)
+    #print(f'Color matrix 1:{color_matrix_1}')
+    #print(f'Color matrix 2:{color_matrix_2}')
+    #print(np.sum(np.abs(np.array(color_matrix_1) - np.array(color_matrix_2))))
+    return color_matrix_1, color_matrix_2
+
+
+def get_orientation(tags, ifds):
+    possible_tags = ['Orientation', 'Image Orientation']
+    return get_values(tags, possible_tags)
+
+
+def get_noise_profile(tags, ifds):
+    possible_keys = ['Image Tag 0xC761', 'Image Tag 51041',
+                     'NoiseProfile', 'Image NoiseProfile']
+    vals = get_values(tags, possible_keys)
+    if vals is None:
+        # print("Noise profile not found in exifread tags. Searching IFDs.")
+        vals = get_tag_values_from_ifds(51041, ifds)
+    return vals
+
+
+def get_values(tags, possible_keys):
+    values = None
+    for key in possible_keys:
+        if key in tags.keys():
+            values = tags[key].values
+    return values
+
+
+def normalize(raw_image, black_level, white_level):
+    if type(black_level) is list and len(black_level) == 1:
+        black_level = float(black_level[0])
+    if type(white_level) is list and len(white_level) == 1:
+        white_level = float(white_level[0])
+    black_level_mask = black_level
+    if type(black_level) is list and len(black_level) == 4:
+        if type(black_level[0]) is Ratio:
+            black_level = ratios2floats(black_level)
+        if type(black_level[0]) is Fraction:
+            black_level = fractions2floats(black_level)
+        black_level_mask = np.zeros(raw_image.shape)
+        idx2by2 = [[0, 0], [0, 1], [1, 0], [1, 1]]
+        step2 = 2
+        for i, idx in enumerate(idx2by2):
+            black_level_mask[idx[0]::step2, idx[1]::step2] = black_level[i]
+    normalized_image = raw_image.astype(np.float32) - black_level_mask
+    # if some values were smaller than black level
+    normalized_image[normalized_image < 0] = 0
+    normalized_image = normalized_image / (white_level - black_level_mask)
+    return normalized_image
+
+
+def ratios2floats(ratios):
+    floats = []
+    for ratio in ratios:
+        floats.append(float(ratio.num) / ratio.den)
+    return floats
+
+
+def fractions2floats(fractions):
+    floats = []
+    for fraction in fractions:
+        floats.append(float(fraction.numerator) / fraction.denominator)
+    return floats
+
+
+def illumination_parameters_estimation(current_image, illumination_estimation_option):
+    ie_method = illumination_estimation_option.lower()
+    if ie_method == "gw":
+        ie = np.mean(current_image, axis=(0, 1))
+        ie /= ie[1]
+        return ie
+    elif ie_method == "sog":
+        sog_p = 4.
+        ie = np.mean(current_image**sog_p, axis=(0, 1))**(1 / sog_p)
+        ie /= ie[1]
+        return ie
+    elif ie_method == "wp":
+        ie = np.max(current_image, axis=(0, 1))
+        ie /= ie[1]
+        return ie
+    elif ie_method == "iwp":
+        samples_count = 20
+        sample_size = 20
+        rows, cols = current_image.shape[:2]
+        data = np.reshape(current_image, (rows * cols, 3))
+        maxima = np.zeros((samples_count, 3))
+        for i in range(samples_count):
+            maxima[i, :] = np.max(data[np.random.randint(
+                low=0, high=rows * cols, size=(sample_size)), :], axis=0)
+        ie = np.mean(maxima, axis=0)
+        ie /= ie[1]
+        return ie
+    else:
+        raise ValueError(
+            'Bad illumination_estimation_option value! Use the following options: "gw", "wp", "sog", "iwp"')
+
+
+def white_balance(demosaic_img, as_shot_neutral):
+    if type(as_shot_neutral[0]) is Ratio:
+        as_shot_neutral = ratios2floats(as_shot_neutral)
+
+    as_shot_neutral = np.asarray(as_shot_neutral)
+    # transform vector into matrix
+    if as_shot_neutral.shape == (3,):
+        as_shot_neutral = np.diag(1. / as_shot_neutral)
+
+    assert as_shot_neutral.shape == (3, 3)
+
+    white_balanced_image = np.dot(demosaic_img, as_shot_neutral.T)
+    white_balanced_image = np.clip(white_balanced_image, 0.0, 1.0)
+
+    return white_balanced_image
+
+
+def simple_demosaic(img, cfa_pattern):
+    raw_colors = np.asarray(cfa_pattern).reshape((2, 2))
+    demosaiced_image = np.zeros((img.shape[0] // 2, img.shape[1] // 2, 3))
+    for i in range(2):
+        for j in range(2):
+            ch = raw_colors[i, j]
+            if ch == 1:
+                demosaiced_image[:, :, ch] += img[i::2, j::2] / 2
+            else:
+                demosaiced_image[:, :, ch] = img[i::2, j::2]
+    return demosaiced_image
+
+
+def denoise_image(demosaiced_image):
+    current_image = denoise_bilateral(
+        demosaiced_image, sigma_color=None, sigma_spatial=0.01, channel_axis=2, mode='reflect')
+    return current_image
+
+
+def apply_color_space_transform(demosaiced_image, color_matrix_1, color_matrix_2):
+    if isinstance(color_matrix_1[0], Fraction):
+        color_matrix_1 = fractions2floats(color_matrix_1)
+    if isinstance(color_matrix_2[0], Fraction):
+        color_matrix_2 = fractions2floats(color_matrix_2)
+
+    xyz2cam1 = np.reshape(np.asarray(color_matrix_1), (3, 3))
+    xyz2cam2 = np.reshape(np.asarray(color_matrix_2), (3, 3))
+
+    # normalize rows (needed?)
+
+    xyz2cam1 = xyz2cam1 / np.sum(xyz2cam1, axis=1, keepdims=True)
+    xyz2cam2 = xyz2cam2 / np.sum(xyz2cam1, axis=1, keepdims=True)
+
+    # inverse
+    cam2xyz1 = np.linalg.inv(xyz2cam1)
+    cam2xyz2 = np.linalg.inv(xyz2cam2)
+
+    # cam2xyz1 = cam2xyz1 * 0.9 + cam2xyz2 * 0.1
+    # for now, use one matrix  # TODO: interpolate btween both
+
+    # simplified matrix multiplication
+    # xyz_image = cam2xyz1[np.newaxis, np.newaxis, :, :] * \
+    #     demosaiced_image[:, :, np.newaxis, :]
+    # xyz_image = np.sum(xyz_image, axis=-1)
+
+    xyz_image = np.einsum('kc,ijc', cam2xyz1, demosaiced_image)
+
+    xyz_image = np.clip(xyz_image, 0.0, 1.0)
+
+    return xyz_image
+
+
+def srgb2xyz(xyz_image):
+    srgb2xyz = np.array([[0.4124564, 0.3575761, 0.1804375],
+                         [0.2126729, 0.7151522, 0.0721750],
+                         [0.0193339, 0.1191920, 0.9503041]])
+
+    out = srgb2xyz[np.newaxis, np.newaxis,
+                   :, :] * xyz_image[:, :, np.newaxis, :]
+    out = np.sum(out, axis=-1)
+    out = np.clip(out, 0.0, 1.0)
+    return out
+
+
+def transform_xyz_to_srgb(xyz_image):
+    # srgb2xyz = np.array([[0.4124564, 0.3575761, 0.1804375],
+    #                      [0.2126729, 0.7151522, 0.0721750],
+    #                      [0.0193339, 0.1191920, 0.9503041]])
+
+    # xyz2srgb = np.linalg.inv(srgb2xyz)
+
+    xyz2srgb = np.array([[2.0413690, -0.5649464, -0.3446944],
+                         [-0.9692660, 1.8760108, 0.0415560],
+                         [0.0134474, -0.1183897, 1.0154096]])
+    # xyz2srgb = np.array([[3.2404542, -1.5371385, -0.4985314],
+    #                      [-0.9692660, 1.8760108, 0.0415560],
+    #                      [0.0556434, -0.2040259, 1.0572252]])
+
+    # normalize rows (needed?)
+    xyz2srgb = xyz2srgb / np.sum(xyz2srgb, axis=-1, keepdims=True)
+
+    srgb_image = xyz2srgb[np.newaxis, np.newaxis,
+                          :, :] * xyz_image[:, :, np.newaxis, :]
+    srgb_image = np.sum(srgb_image, axis=-1)
+    srgb_image = np.clip(srgb_image, 0.0, 1.0)
+    return srgb_image
+
+
+def reverse_orientation(image, orientation):
+    # 1 = Horizontal(normal)
+    # 2 = Mirror horizontal
+    # 3 = Rotate 180
+    # 4 = Mirror vertical
+    # 5 = Mirror horizontal and rotate 270 CW
+    # 6 = Rotate 90 CW
+    # 7 = Mirror horizontal and rotate 90 CW
+    # 8 = Rotate 270 CW
+    rev_orientations = np.array([1, 2, 3, 4, 5, 8, 7, 6])
+    return fix_orientation(image, rev_orientations[orientation - 1])
+
+
+def apply_gamma(x):
+    # return x ** (1.0 / 2.2)
+    x = x.copy()
+    idx = x <= 0.0031308
+    x[idx] *= 12.92
+    x[idx == False] = (x[idx == False] ** (1.0 / 2.4)) * 1.055 - 0.055
+    return x
+
+
+def apply_tone_map(x, tone_mapping='Base'):
+    if tone_mapping == 'Flash':
+        return perform_flash(x, perform_gamma_correction=0) / 255.
+    elif tone_mapping == 'Storm':
+        return perform_storm(x, perform_gamma_correction=0) / 255.
+    elif tone_mapping == 'Drago':
+        tonemap = cv2.createTonemapDrago()
+        return tonemap.process(x.astype(np.float32))
+    elif tone_mapping == 'Mantiuk':
+        tonemap = cv2.createTonemapMantiuk()
+        return tonemap.process(x.astype(np.float32))
+    elif tone_mapping == 'Reinhard':
+        tonemap = cv2.createTonemapReinhard()
+        return tonemap.process(x.astype(np.float32))
+    elif tone_mapping == 'Linear':
+        return np.clip(x / np.sort(x.flatten())[-50000], 0, 1)
+    elif tone_mapping == 'Base':
+        # return 3 * x ** 2 - 2 * x ** 3
+        # tone_curve = loadmat('tone_curve.mat')
+        tone_curve = loadmat(os.path.join(os.path.dirname(
+            os.path.realpath(__file__)), 'tone_curve.mat'))
+        tone_curve = tone_curve['tc']
+        x = np.round(x * (len(tone_curve) - 1)).astype(int)
+        tone_mapped_image = np.squeeze(tone_curve[x])
+        return tone_mapped_image
+    else:
+        raise ValueError(
+            'Bad tone_mapping option value! Use the following options: "Base", "Flash", "Storm", "Linear", "Drago", "Mantiuk", "Reinhard"')
+
+
+def autocontrast(output_image, cutoff_prcnt=2, preserve_tone=False):
+    if preserve_tone:
+        min_val, max_val = np.percentile(
+            output_image, [cutoff_prcnt, 100 - cutoff_prcnt])
+        output_image = (output_image - min_val) / (max_val - min_val)
+    else:
+        channels = [None] * 3
+        for ch in range(3):
+            min_val, max_val = np.percentile(
+                output_image[..., ch], [cutoff_prcnt, 100 - cutoff_prcnt])
+            channels[ch] = (output_image[..., ch] - min_val) / \
+                (max_val - min_val)
+        output_image = np.dstack(channels)
+    output_image = np.clip(output_image, 0, 1)
+    return output_image
+
+
+def autocontrast_using_pil(img, cutoff=2):
+    img_uint8 = np.clip(255 * img, 0, 255).astype(np.uint8)
+    img_pil = Image.fromarray(img_uint8)
+    img_pil = ImageOps.autocontrast(img_pil, cutoff=cutoff)
+    output_image = np.array(img_pil).astype(np.float32) / 255
+    return output_image
+
+
+def raw_rgb_to_cct(rawRgb, xyz2cam1, xyz2cam2):
+    """Convert raw-RGB triplet to corresponding correlated color temperature (CCT)"""
+    pass
+    # pxyz = [.5, 1, .5]
+    # loss = 1e10
+    # k = 1
+    # while loss > 1e-4:
+    #     cct = XyzToCct(pxyz)
+    #     xyz = RawRgbToXyz(rawRgb, cct, xyz2cam1, xyz2cam2)
+    #     loss = norm(xyz - pxyz)
+    #     pxyz = xyz
+    #     fprintf('k = %d, loss = %f\n', [k, loss])
+    #     k = k + 1
+    # end
+    # temp = cct
+
+
+def resize_using_skimage(img, width=1296, height=864):
+    out_shape = (height, width) + img.shape[2:]
+    if img.shape == out_shape:
+        return img
+    out_img = skimage_resize(
+        img, out_shape, preserve_range=True, anti_aliasing=True)
+    out_img = out_img.astype(np.uint8)
+    return out_img
+
+
+def resize_using_pil(img, width=1296, height=864):
+    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.ANTIALIAS)
+    out_img = np.array(out_img)
+    return out_img
+
+
+def fix_orientation(image, orientation):
+    # 1 = Horizontal(normal)
+    # 2 = Mirror horizontal
+    # 3 = Rotate 180
+    # 4 = Mirror vertical
+    # 5 = Mirror horizontal and rotate 270 CW
+    # 6 = Rotate 90 CW
+    # 7 = Mirror horizontal and rotate 90 CW
+    # 8 = Rotate 270 CW
+
+    map_ = {'Horizontal (normal)': 1,
+            'Mirror horizontal': 2,
+            'Rotate 180': 3,
+            'Mirror vertical': 4,
+            'Mirror horizontal and rotate 270 CW': 5,
+            'Rotate 90 CW': 6,
+            'Mirror horizontal and rotate 90 CW': 7,
+            'Rotate 270 CW': 8
+
+            }
+
+    if type(orientation) is list:
+        orientation = orientation[0]
+
+    orientation = map_[orientation]
+
+    if orientation == 1:
+        pass
+    elif orientation == 2:
+        image = cv2.flip(image, 0)
+    elif orientation == 3:
+        image = cv2.rotate(image, cv2.ROTATE_180)
+    elif orientation == 4:
+        image = cv2.flip(image, 1)
+    elif orientation == 5:
+        image = cv2.flip(image, 0)
+        image = cv2.rotate(image, cv2.ROTATE_90_COUNTERCLOCKWISE)
+    elif orientation == 6:
+        image = cv2.rotate(image, cv2.ROTATE_90_CLOCKWISE)
+    elif orientation == 7:
+        image = cv2.flip(image, 0)
+        image = cv2.rotate(image, cv2.ROTATE_90_CLOCKWISE)
+    elif orientation == 8:
+        image = cv2.rotate(image, cv2.ROTATE_90_COUNTERCLOCKWISE)
+    else:
+        raise NotImplementedError('Orientation not defined')
+
+    return image

IVL/requirements.txt

@@ -0,0 +1,9 @@
+ExifRead==2.3.2
+matplotlib==3.5.1
+numpy==1.24.2
+opencv_python==4.5.5.62
+Pillow==10.2.0
+rawpy==0.17.0
+scipy==1.9.1
+scikit-image==0.20.0
+tqdm==4.62.3

IVL/run.sh

@@ -0,0 +1,3 @@
+#!/usr/bin/env bash
+
+python pipeline24.py -p /data/ -o /data/

IVL/utils/__init__.py

@@ -0,0 +1,36 @@
+from fractions import Fraction
+from pathlib import Path
+from json import JSONEncoder
+from .utils import *
+
+
+def rmtree(path: Path):
+    if path.is_file():
+        path.unlink()
+    else:
+        for ch in path.iterdir():
+            rmtree(ch)
+        path.rmdir()
+
+
+def safe_save(fpath, data, save_fun, rewrite=False, error_msg='File {fpath} exists! To rewite it use `--rewrite` flag', **kwargs):
+    if not fpath.is_file() or rewrite:
+        save_fun(str(fpath), data, **kwargs)
+    else:
+        raise FileExistsError(error_msg.format(fpath=fpath))
+
+
+class FractionJSONEncoder(JSONEncoder):
+    def default(self, o):
+        if isinstance(o, Fraction):
+            return {'Fraction': [o.numerator, o.denominator]}
+        else:
+            return o.__dict__
+
+
+def fraction_from_json(json_object):
+    if 'Fraction' in json_object:
+        return Fraction(*json_object['Fraction'])
+    return json_object
+
+

IVL/utils/utils.py

@@ -0,0 +1,56 @@
+from PIL import Image
+import json
+import os
+
+def json_read(fname, **kwargs):
+    with open(fname) as j:
+        data = json.load(j, **kwargs)
+    return data
+
+
+def json_save(fname, data, indent_len=4, **kwargs):
+    with open(fname, "w") as f:
+        s = json.dumps(data, sort_keys=True, ensure_ascii=False,
+                       indent=" " * indent_len, **kwargs)
+        f.write(s)
+
+
+def process_wb_from_txt(txt_path):
+    with open(txt_path, 'r') as fh:
+        txt = [line.rstrip().split() for line in fh]
+
+    txt = [[float(k) for k in row] for row in txt]
+
+    assert len(txt) in [1, 3]
+
+    if len(txt) == 1:
+        # wb vector
+        txt = txt[0]
+
+    return txt
+
+
+def process_ids_from_txt(txt_path):
+    with open(txt_path, 'r') as fh:
+        temp = fh.read().splitlines()
+    return temp
+
+
+def save_txt(p, s):
+    with open(p, 'w') as text_file:
+        text_file.write(s)
+
+
+def downscale_jpg(img_path, new_shape, quality_perc=100):
+    img = Image.open(img_path)
+    if (img.size[0], img.size[1]) != new_shape:
+        new_img = img.resize(new_shape, Image.ANTIALIAS)
+        new_img.save(img_path[:-len('.jpg')] + '.jpg',
+                     'JPEG', quality=quality_perc)
+
+
+def rename_img(img_path):
+    if img_path.lower().endswith('jpeg'):
+        os.rename(img_path, img_path[:-len('jpeg')] + 'jpg')
+    else:
+        os.rename(img_path, img_path[:-len('JPG')] + 'jpg')