|
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 |
|
|
|
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_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]; |
|
|
|
|
|
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 = 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) |
|
|
|
|
|
|
|
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) |
|
|
|
|
|
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) |
|
|
|
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] |
|
|
|
num_gray_pixels = int(np.floor(n * tot_pixels / 100)) |
|
|
|
|
|
if denoise_first: |
|
sigma_est = 1 |
|
img_ = cv2.GaussianBlur(img,(0,0),5) |
|
else: |
|
img_ = img |
|
|
|
|
|
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 + (1 - alpha) * |
|
((x - alpha) / (1 - alpha)).clip(min=0) ** lmbd |
|
) |
|
else: |
|
out = np.where(x <= alpha, |
|
alpha - alpha * (1 - x / alpha) ** lmbd, |
|
x |
|
) |
|
|
|
|
|
|
|
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) |
|
|
|
|
|
|
|
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 = {} |
|
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 |
|
|