diff --git a/data/.gitignore b/data/.gitignore new file mode 100644 index 0000000..e7f6a90 --- /dev/null +++ b/data/.gitignore @@ -0,0 +1 @@ +fivek/ diff --git a/data/.gitkeep b/data/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/eiuie/base_model.py b/eiuie/base_model.py new file mode 100644 index 0000000..faa2ee8 --- /dev/null +++ b/eiuie/base_model.py @@ -0,0 +1,113 @@ +from abc import ABC, abstractmethod + +import numpy as np + + +def BGR2HSI(image: np.ndarray) -> np.ndarray: + """ + Convert image from BGR to HSI. + + Parameters + ---------- + image : np.ndarray + Image to be converted, as BGR. + + Returns + ------- + np.ndarray + Converted image, as HSI. + """ + + # Normalize BGR values to [0,1] + bgr = image.astype(np.float32) / 255.0 + blue = bgr[:, :, 0] + green = bgr[:, :, 1] + red = bgr[:, :, 2] + + # Compute intensity + I = (blue + green + red) / 3.0 + + # Compute saturation + min_val = np.minimum(np.minimum(blue, green), red) + S = 1 - 3.0 / (blue + green + red + 1e-6) * min_val + + # Compute hue + num = 0.5 * ((red - green) + (red - blue)) + den = np.sqrt((red - green) ** 2 + (red - blue) * (green - blue)) + theta = np.arccos(num / (den + 1e-6)) + H = theta + H[blue > green] = 2 * np.pi - H[blue > green] + H /= 2 * np.pi + + hsi = np.stack([H, S, I], axis=2) + return hsi + + +def HSI2BGR(image: np.ndarray) -> np.ndarray: + """ + Convert image from HSI to BGR. + + Parameters + ---------- + image : np.ndarray + Image to be converted, as HSI. + + Returns + ------- + np.ndarray + Converted image, as BGR. + """ + + H = image[:, :, 0] * 2 * np.pi + S = image[:, :, 1] + I = image[:, :, 2] + + R = np.zeros_like(H) + G = np.zeros_like(H) + B = np.zeros_like(H) + + # RG sector + cond = np.logical_and(0 <= H, H < 2 * np.pi / 3) + B[cond] = I[cond] * (1 - S[cond]) + R[cond] = I[cond] * (1 + S[cond] * np.cos(H[cond]) / np.cos(np.pi / 3 - H[cond])) + G[cond] = 3 * I[cond] - (R[cond] + B[cond]) + + # GB sector + cond = np.logical_and(2 * np.pi / 3 <= H, H < 4 * np.pi / 3) + H[cond] = H[cond] - 2 * np.pi / 3 + R[cond] = I[cond] * (1 - S[cond]) + G[cond] = I[cond] * (1 + S[cond] * np.cos(H[cond]) / np.cos(np.pi / 3 - H[cond])) + B[cond] = 3 * I[cond] - (R[cond] + G[cond]) + + # BR sector + cond = np.logical_and(4 * np.pi / 3 <= H, H < 2 * np.pi) + H[cond] = H[cond] - 4 * np.pi / 3 + G[cond] = I[cond] * (1 - S[cond]) + B[cond] = I[cond] * (1 + S[cond] * np.cos(H[cond]) / np.cos(np.pi / 3 - H[cond])) + R[cond] = 3 * I[cond] - (G[cond] + B[cond]) + + bgr = np.stack([B, G, R], axis=2) + return (bgr * 255).astype(np.uint8) + + +class BaseModel(ABC): + """ + Abstract class for all models. + """ + + @abstractmethod + def process_image(self, image: np.ndarray) -> np.ndarray: + """ + Process image using the model. + + Parameters + ---------- + image : np.ndarray + Image to be processed, as BGR. + + Returns + ------- + np.ndarray + Processed image, as BGR. + """ + ... diff --git a/eiuie/download.py b/eiuie/download.py new file mode 100644 index 0000000..8c3d993 --- /dev/null +++ b/eiuie/download.py @@ -0,0 +1,91 @@ +import os +import requests +from pathlib import Path +import rawpy +from PIL import Image +from tqdm import tqdm +from joblib import Parallel + + +class ProgressParallel(Parallel): + def __init__(self, use_tqdm=True, total=None, *args, **kwargs): + self._use_tqdm = use_tqdm + self._total = total + super().__init__(*args, **kwargs) + + def __call__(self, *args, **kwargs): + with tqdm(disable=not self._use_tqdm, total=self._total) as self._pbar: + return Parallel.__call__(self, *args, **kwargs) + + def print_progress(self): + if self._total is None: + self._pbar.total = self.n_dispatched_tasks + self._pbar.n = self.n_completed_tasks + self._pbar.refresh() + + +def dng_to_jpg(img_path: str): + raw = rawpy.imread(img_path) + rgb = raw.postprocess(use_camera_wb=True) + new_name = img_path[:-3] + "jpg" + Image.fromarray(rgb).save(new_name, optimize=True) + os.remove(img_path) + + +def tif_to_jpg(img_path: str): + img = Image.open(img_path) + img.save(img_path[:-3] + "jpg") + os.remove(img_path) + + +def download_img(img_name: str, data_dir: Path): + try: + dng_url = f"https://data.csail.mit.edu/graphics/fivek/img/dng/{img_name}.dng" + dng_file = requests.get(dng_url, allow_redirects=True) + dng_path = data_dir / "raw" / (img_name + ".dng") + open(dng_path, "wb").write(dng_file.content) + dng_to_jpg(str(dng_path)) + + tif_base = "https://data.csail.mit.edu/graphics/fivek/img/tiff16" + for expert in ["a", "b", "c", "d", "e"]: + url = f"{tif_base}_{expert}/{img_name}.tif" + path = data_dir / expert / (img_name + ".tif") + tif_file = requests.get(url, allow_redirects=True) + open(path, "wb").write(tif_file.content) + tif_to_jpg(str(path)) + + except Exception as e: + print(img_name, e) + + +def download_dataset(store_dir: Path, n_jobs: int = 8): + # * Create folders + dng_dir = store_dir / "raw" + tif_dirs = [store_dir / s for s in ["a", "b", "c", "d", "e"]] + + dng_dir.mkdir(parents=True, exist_ok=True) + for path in tif_dirs: + path.mkdir(parents=True, exist_ok=True) + + # * Get image info + f1 = requests.get( + "https://data.csail.mit.edu/graphics/fivek/legal/filesAdobe.txt" + ).text.split("\n") + f2 = requests.get( + "https://data.csail.mit.edu/graphics/fivek/legal/filesAdobeMIT.txt" + ).text.split("\n") + names = [x for x in set(f1 + f2) if x != ""] + + # * Download imgs + ProgressParallel(n_jobs=n_jobs, total=len(names))( + download_img(name, store_dir) for name in names + ) + + +def main(): + store_dir = Path("data/fivek") + download_dataset(store_dir, os.cpu_count() or 1) + + +if __name__ == "__main__": + main() diff --git a/eiuie/homomorphic_filtering.py b/eiuie/homomorphic_filtering.py new file mode 100644 index 0000000..aaf1ee7 --- /dev/null +++ b/eiuie/homomorphic_filtering.py @@ -0,0 +1,99 @@ +import numpy as np +import cv2 + +import base_model as bm + + +class HomomorphicFiltering(bm.BaseModel): + """Homomorphic Filtering""" + + ksize: int + sigma: float + gamma_1: float + gamma_2: float + rho: float + + def __init__(self, ksize: int = 3, sigma: float = 1.0): + """ + Parameters + ---------- + ksize : int, optional + Kernel size, by default 3 + sigma : float, optional + Gaussian kernel standard deviation, by default 1.0 + """ + self.ksize = ksize + self.sigma = sigma + self.gamma_1 = 0.8 + self.gamma_2 = 1.8 + self.rho = 100.0 + + def filter(self, value): + return self.gamma_1 - self.gamma_2 * ( + 1 / (1 + 2.415 * np.power(value / self.rho, 4)) + ) + + def _process_image(self, hsi: np.ndarray) -> np.ndarray: + """ + Process image using the model. + + Parameters + ---------- + image : np.ndarray + Image to be processed, as HSI. + + Returns + ------- + np.ndarray + Processed image, as HSI. + """ + i = hsi[:, :, 2] + i_log = np.log2(i + 1.0) + i_log_fft_shifted = np.fft.fftshift(np.fft.fft2(i_log)) + i_log_fft_shifted_filtered = np.zeros_like(i_log_fft_shifted) + for i in range(i_log_fft_shifted.shape[0]): + for j in range(i_log_fft_shifted.shape[1]): + i_log_fft_shifted_filtered[i, j] = i_log_fft_shifted[ + i, j + ] * self.filter(np.sqrt(i**2 + j**2)) + i_log_filtered = np.real( + np.fft.ifft2(np.fft.ifftshift(i_log_fft_shifted_filtered)) + ) + i_filtered = np.exp2(i_log_filtered) - 1.0 + hsi_filtered = hsi.copy() + hsi_filtered[:, :, 2] = i_filtered + return hsi_filtered + + def process_image(self, image: np.ndarray) -> np.ndarray: + """ + Process image using the model. + + Parameters + ---------- + image : np.ndarray + Image to be processed, as BGR. + + Returns + ------- + np.ndarray + Processed image, as BGR. + """ + og = image.copy() + image = image.astype(np.float32) + hsi = bm.BGR2HSI(image) + hsi = self._process_image(hsi) + image = bm.HSI2BGR(hsi) + image = np.clip(image, 0, 255) + image = image.astype(np.uint8) + image = cv2.cvtColor(image, cv2.COLOR_BGR2HSV) + image[:, :, 2] = cv2.equalizeHist(image[:, :, 2]) + image = cv2.cvtColor(image, cv2.COLOR_HSV2BGR) + # show difference image + diff = np.abs(og.astype(np.float32) - image.astype(np.float32)) + diff = diff.astype(np.uint8) + cv2.imshow("diff", diff) + cv2.waitKey(0) + # show og + cv2.imshow("og", og) + cv2.waitKey(0) + return image diff --git a/eiuie/main.py b/eiuie/main.py new file mode 100644 index 0000000..f52bc82 --- /dev/null +++ b/eiuie/main.py @@ -0,0 +1,65 @@ +import argparse + +import cv2 + +import download as download_module +import base_model as bm +import unsharp_masking +import retinex +import homomorphic_filtering + + +def main(): + parser = argparse.ArgumentParser(description="Entry point for the eiuie CLI.") + + parser.add_argument( + "command", + type=str, + choices=["download", "single"], + help="Command to run", + ) + + # --method=xyz + parser.add_argument( + "--method", + type=str, + default="unsharp_masking", + choices=["unsharp_masking", "retinex", "homomorphic_filtering"], + help="Filter method to use", + ) + + # --file=xyz + parser.add_argument( + "--file", + type=str, + help="Path to image file to process", + ) + + args = parser.parse_args() + + method: bm.BaseModel + match args.method: + case "unsharp_masking": + method = unsharp_masking.UnsharpMasking() + case "retinex": + method = retinex.Retinex() + case "homomorphic_filtering": + method = homomorphic_filtering.HomomorphicFiltering() + case _: + raise ValueError(f"Unknown method: {args.method}") + + match args.command: + case "download": + download_module.main() + case "single": + image = cv2.imread(args.file) + processed_image = method.process_image(image) + # show image + cv2.imshow("image", processed_image) + cv2.waitKey() + case _: + raise ValueError(f"Unknown command: {args.command}") + + +if __name__ == "__main__": + main() diff --git a/eiuie/retinex.py b/eiuie/retinex.py new file mode 100644 index 0000000..fa9ff70 --- /dev/null +++ b/eiuie/retinex.py @@ -0,0 +1,180 @@ +import numpy as np +import cv2 + +import base_model as bm + + +def get_ksize(sigma): + # opencv calculates ksize from sigma as + # sigma = 0.3*((ksize-1)*0.5 - 1) + 0.8 + # then ksize from sigma is + # ksize = ((sigma - 0.8)/0.15) + 2.0 + + return int(((sigma - 0.8) / 0.15) + 2.0) + + +def get_gaussian_blur(img, ksize=0, sigma=5): + # if ksize == 0, then compute ksize from sigma + if ksize == 0: + ksize = get_ksize(sigma) + + # Gaussian 2D-kernel can be seperable into 2-orthogonal vectors + # then compute full kernel by taking outer product or simply mul(V, V.T) + sep_k = cv2.getGaussianKernel(ksize, sigma) + + # if ksize >= 11, then convolution is computed by applying fourier transform + return cv2.filter2D(img, -1, np.outer(sep_k, sep_k)) + + +def ssr(img, sigma): + # Single-scale retinex of an image + # SSR(x, y) = log(I(x, y)) - log(I(x, y)*F(x, y)) + # F = surrounding function, here Gaussian + + return np.log10(img) - np.log10(get_gaussian_blur(img, ksize=0, sigma=sigma) + 1.0) + + +def msr(img, sigma_scales=[15, 80, 250]): + # Multi-scale retinex of an image + # MSR(x,y) = sum(weight[i]*SSR(x,y, scale[i])), i = {1..n} scales + + msr = np.zeros(img.shape) + # for each sigma scale compute SSR + for sigma in sigma_scales: + msr += ssr(img, sigma) + + # divide MSR by weights of each scale + # here we use equal weights + msr = msr / len(sigma_scales) + + # computed MSR could be in range [-k, +l], k and l could be any real value + # so normalize the MSR image values in range [0, 255] + msr = cv2.normalize(msr, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8UC3) + + return msr + + +def color_balance(img, low_per, high_per): + """Contrast stretch img by histogram equilization with black and white cap""" + + tot_pix = img.shape[1] * img.shape[0] + # no.of pixels to black-out and white-out + low_count = tot_pix * low_per / 100 + high_count = tot_pix * (100 - high_per) / 100 + + # channels of image + ch_list = [] + if len(img.shape) == 2: + ch_list = [img] + else: + ch_list = cv2.split(img) + + cs_img = [] + # for each channel, apply contrast-stretch + for i in range(len(ch_list)): + ch = ch_list[i] + # cummulative histogram sum of channel + cum_hist_sum = np.cumsum(cv2.calcHist([ch], [0], None, [256], (0, 256))) + + # find indices for blacking and whiting out pixels + li, hi = np.searchsorted(cum_hist_sum, (low_count, high_count)) + if li == hi: + cs_img.append(ch) + continue + # lut with min-max normalization for [0-255] bins + lut = np.array( + [ + 0 if i < li else (255 if i > hi else round((i - li) / (hi - li) * 255)) + for i in np.arange(0, 256) + ], + dtype="uint8", + ) + # constrast-stretch channel + cs_ch = cv2.LUT(ch, lut) + cs_img.append(cs_ch) + + if len(cs_img) == 1: + return np.squeeze(cs_img) + elif len(cs_img) > 1: + return cv2.merge(cs_img) + return None + + +def msrcr( + img, + sigma_scales=[15, 80, 250], + alpha=125, + beta=46, + G=192, + b=-30, + low_per=1, + high_per=1, +): + # Multi-scale retinex with Color Restoration + # MSRCR(x,y) = G * [MSR(x,y)*CRF(x,y) - b], G=gain and b=offset + # CRF(x,y) = beta*[log(alpha*I(x,y) - log(I'(x,y))] + # I'(x,y) = sum(Ic(x,y)), c={0...k-1}, k=no.of channels + + img = img.astype(np.float64) + 1.0 + # Multi-scale retinex and don't normalize the output + msr_img = msr(img, sigma_scales, apply_normalization=False) + # Color-restoration function + crf = beta * (np.log10(alpha * img) - np.log10(np.sum(img, axis=2, keepdims=True))) + # MSRCR + msrcr = G * (msr_img * crf - b) + # normalize MSRCR + msrcr = cv2.normalize(msrcr, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8UC3) + # color balance the final MSRCR to flat the histogram distribution with tails on both sides + msrcr = color_balance(msrcr, low_per, high_per) + + return msrcr + + +def msrcp(img, sigma_scales=[15, 80, 250], low_per=1, high_per=1): + # Multi-scale retinex with Color Preservation + # Int(x,y) = sum(Ic(x,y))/3, c={0...k-1}, k=no.of channels + # MSR_Int(x,y) = MSR(Int(x,y)), and apply color balance + # B(x,y) = MAX_VALUE/max(Ic(x,y)) + # A(x,y) = max(B(x,y), MSR_Int(x,y)/Int(x,y)) + # MSRCP = A*I + + # Intensity image (Int) + int_img = (np.sum(img, axis=2) / img.shape[2]) + 1.0 + # Multi-scale retinex of intensity image (MSR) + msr_int = msr(int_img, sigma_scales) + # color balance of MSR + msr_cb = color_balance(msr_int, low_per, high_per) + + # B = MAX/max(Ic) + B = 256.0 / (np.max(img, axis=2) + 1.0) + # BB = stack(B, MSR/Int) + BB = np.array([B, msr_cb / int_img]) + # A = min(BB) + A = np.min(BB, axis=0) + # MSRCP = A*I + msrcp = np.clip(np.expand_dims(A, 2) * img, 0.0, 255.0) + + return msrcp.astype(np.uint8) + + +class Retinex(bm.BaseModel): + def __init__( + self, + sigma_scales=[15, 80, 250], + alpha=125, + beta=46, + G=192, + b=-30, + low_per=1, + high_per=1, + ): + self.sigma_scales = sigma_scales + self.alpha = alpha + self.beta = beta + self.G = G + self.b = b + self.low_per = low_per + self.high_per = high_per + + def process_image(self, img): + return msrcp(img) diff --git a/eiuie/unsharp_masking.py b/eiuie/unsharp_masking.py new file mode 100644 index 0000000..55c1612 --- /dev/null +++ b/eiuie/unsharp_masking.py @@ -0,0 +1,41 @@ +import numpy as np +import cv2 + +import base_model as bm + + +class UnsharpMasking(bm.BaseModel): + """Unsharp Masking""" + + ksize: int + sigma: float + + def __init__(self, ksize: int = 33, sigma: float = 20.0): + """ + Parameters + ---------- + ksize : int, optional + Kernel size, by default 33 + sigma : float, optional + Gaussian kernel standard deviation, by default 20.0 + """ + self.ksize = ksize + self.sigma = sigma + + def process_image(self, image: np.ndarray) -> np.ndarray: + """ + Process image using the model. + + Parameters + ---------- + image : np.ndarray + Image to be processed, as BGR. + + Returns + ------- + np.ndarray + Processed image, as BGR. + """ + blurred = cv2.GaussianBlur(image, (self.ksize, self.ksize), self.sigma) + sharpened = cv2.addWeighted(image, 2, blurred, -1, 0) + return sharpened diff --git a/pyproject.toml b/pyproject.toml index fe9811b..5d50bc7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,14 +10,23 @@ classifiers = [ "Operating System :: OS Independent", ] requires = [ + "torch", + "tdqm", + "mpire", + "numpy", + "opencv-python", ] [build-system] requires = ["hatchling"] build-backend = "hatchling.build" +[project.scripts] +eiui = "eiuie.main:main" + [tool.hatch.build.targets.wheel] -only-include = ["eiuie"] +only-include = ["eiuie", "data"] [tool.hatch.build.targets.wheel.sources] "eiuie" = "eiuie" +"data" = "eiuie/data"