main
原始文件为 .py 代码，本文是转换后的 Markdown 文件。
```.py
import cv2
import math
import numpy as np
from scipy import ndimage
def relinear_img(img):
result = np.copy(img)
idx = img < 0.0404482
result[idx] = img[idx] / 12.92
result[~idx] = ((img[~idx] + 0.055) / 1.055)**(2.4)
return result
def save_img(name, img):
img = np.round(img * 255).clip(0, 255).astype(np.uint8)
cv2.imwrite(name, img)
颜色的高斯核 g (color)
def compute_g(delta, sigmaColor):
# muColor = (sigmaColor)2
g = np.exp(-((delta/sigmaColor)2) / 2)
return g
# g = g / (math.sqrt(2*np.pi) * sigmaColor)
# g 加上一个很小数字，这是为了防止后续计算 g*f 时有的窗口全是 0
# return g / np.sum(g) + 1e-100
空域的高斯核 f (spatial)，通过 kernel size >= ceil(6 sigma) 计算
def compute_f(d, sigmaSpace):
if d is None:
d = 2 * math.ceil(sigmaSpace * 3)
r = d // 2
x, y = np.meshgrid(np.linspace(-r, r, d), np.linspace(-r, r, d))
f = np.exp(-(x**2 + y**2) / (sigmaSpace**2) / 2)
return f
Piecewise Bilateral Filtering
From Paper: Fast Bilateral Filtering for the Display of High-Dynamic-Range Images
def PiecewiseBilateral(img, d, sigmaColor, sigmaSpace):
# 计算最大最小值，一般都是 0 和 1
maxI, minI = img.max(), img.min()
# f = compute_f(sigmaSpace)
# # 一开始分段方式，想通过 sigmaColor 来求，最后发现直接默认 100 就行
# NB_SEGMENTS = math.ceil((maxI - minI) / sigmaColor)
# # 防止 sigma 过小，导致分段太多了
# NB_SEGMENTS = min(20, NB_SEGMENTS)
NB_SEGMENTS = 100
LEN_SEGMENTS = (maxI - minI) / NB_SEGMENTS
# 遍历的时候需要将最后一段包进去，因为遍历是 range(0, NB_SEGMENTS)，尤其是这个等于号，不考虑的话最后一段就不会遍历到
if minI + LEN_SEGMENTS * (NB_SEGMENTS - 1) <= maxI:
NB_SEGMENTS += 1
img = img.astype(float)
result = np.zeros_like(img)
for channel in range(0, 3):
I = img[..., channel]
J = np.zeros_like(I)
for j in range(0, NB_SEGMENTS):
Ij = minI + j * LEN_SEGMENTS
delta = np.abs(I - Ij)
# 对每个像素进行颜色上的处理，即 g(I-now)
Gj = compute_g(delta, sigmaColor)
# # 三种调用卷积方式，测试分别为8.3s、3.2s、0.15s，最后一种由于大核使用FFT所以很快
# # 但是最后一种有误差，在本次例子中会导致严重问题（突出点）！所以弃用
# Kj = signal.convolve2d(Gj, f, mode='same')
# Kj = signal.convolve(Gj, f, mode='same')
# Kj = ndimage.convolve(Gj, f)
# 然而事实上，直接调用 guassian 就可以了，因为高斯核是可拆分的，所以可以用两个 1D 更快
Kj = ndimage.gaussian_filter(Gj, sigmaSpace, radius=d//2)
# V1: 如果 Kj 是 0, 那其实根本就不用管这些点，这些点不可能在 Ij 左右的。
# V2: 但实际上并不如此，如果 sigmaColor 过小，对于某个像素而言，只有完全等于目标颜色才会有值
# 假设一个3*3窗口像素都是 99，假设按照 50 长度遍历，
# 当目标颜色为 50 的时候，窗口内全是 0；当目标颜色为 100 的时候，窗口内也还是 0，有问题。
# 这种时候就应该赋值为目标颜色，这种要特殊处理
zeros = Kj == 0
Kj[zeros] = 1
Hj = Gj * I
Jj = ndimage.gaussian_filter(Hj, sigmaSpace, radius=d//2)
Jj = Jj / Kj
Jj[zeros] = Ij
# 插值
mask = delta < LEN_SEGMENTS
J = J + mask * (1 - delta/LEN_SEGMENTS) * Jj
result[..., channel] = J
return result
Joint Bilateral Filtering
def JointBilateral(noflash, flash, d, sigmaColor, sigmaSpace):
maxI, minI = noflash.max(), noflash.min()
NB_SEGMENTS = 100
LEN_SEGMENTS = (maxI - minI) / NB_SEGMENTS
# 和 PiecewiseBilateral 一样，这里要检查，确保遍历的时候将最后一段包进去
if minI + LEN_SEGMENTS * (NB_SEGMENTS - 1) <= maxI:
NB_SEGMENTS += 1
result = np.zeros_like(flash)
for channel in range(0, 3):
print(channel)
F = flash[..., channel]
A = noflash[..., channel]
J = np.zeros_like(A)
for j in range(0, NB_SEGMENTS):
Ij = minI + j * LEN_SEGMENTS
# 对应 gr(Fp - Fp'), Ij 就是 Fp'
delta = np.abs(F - Ij)
Gj = compute_g(delta, sigmaColor)
Kj = ndimage.gaussian_filter(Gj, sigmaSpace, radius=d//2)
zeros = Kj == 0
Kj[zeros] = 1
# 这里是乘以 noflash
Hj = Gj * A
Jj = ndimage.gaussian_filter(Hj, sigmaSpace, radius=d//2)
Jj = Jj / Kj
Jj[zeros] = Ij
# 插值
mask = delta < LEN_SEGMENTS
J = J + mask * (1 - delta/LEN_SEGMENTS) * Jj
result[..., channel] = J
return result
def gen_mask(Alin, Flin, shad_thr=0.05):
def cv2_flood_fill(img):
contour, _ = cv2.findContours(img, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
for cnt in contour:
cv2.drawContours(img, [cnt], 0, color=255, thickness=-1)
return img
# 转成 YUV
AY = 0.299*Alin[..., 2] + 0.587*Alin[..., 1] + 0.114*Alin[..., 0]
FY = 0.299*Flin[..., 2] + 0.587*Flin[..., 1] + 0.114*Flin[..., 0]
mask_shad = np.zeros(AY.shape, dtype=np.uint8)
mask_shad[np.abs(FY-AY) < shad_thr] = 255
mask_shad = cv2.dilate(mask_shad, np.ones((3, 3)), iterations=1)
mask_shad = cv2_flood_fill(mask_shad)
mask_shad = cv2.erode(mask_shad, np.ones((3, 3)), iterations=1)
mask_shad = cv2_flood_fill(mask_shad)
mask_over = np.zeros(AY.shape, dtype=np.uint8)
mask_over[FY >= 0.95*np.max(FY)] = 255
mask_over = cv2.erode(mask_over, np.ones((3, 3)), iterations=1)
mask_over = cv2_flood_fill(mask_over)
mask_over = cv2.dilate(mask_over, np.ones((3, 3)), iterations=1)
mask_over = cv2_flood_fill(mask_over)
mask = cv2.GaussianBlur(mask_over|mask_shad, ksize=(7, 7), sigmaX=0)
return mask.astype(float) / 255
A = cv2.imread('./src/noflash.tif', -1)
F = cv2.imread('./src/flash.tif', -1)
A = np.power(A/255, 1).astype(np.float32)
F = np.power(F/255, 1).astype(np.float32)
A = A[950:1160, 1550:1550+500]
F = F[950:1160, 1550:1550+500]
delta = (F - A) / 2 + 0.5
save_img('./result/delta.png', delta)
Abase = PiecewiseBilateral(A, d=61, sigmaColor=0.08, sigmaSpace=360)
Anr = JointBilateral(A, F, d=9, sigmaColor=0.05, sigmaSpace=60)
Fbase = PiecewiseBilateral(F, d=5, sigmaColor=0.12, sigmaSpace=100)
eps = 0.02
Fdetail = (F + eps) / (Fbase + eps)
Adetail = Anr * Fdetail
Flin = relinear_img(F)
Alin = relinear_img(A)
mask = gen_mask(Alin, Flin, 0.05)
mask = np.repeat(mask[..., np.newaxis], 3, 2)
AFinal = (1 - mask) * Adetail + mask * Abase
save_img('./result/F.png', F)
save_img('./result/Abase.png', Abase)
save_img('./result/Anr.png', Anr)
save_img('./result/Fbase.png', Fbase)
save_img('./result/Fdetail.png', Fdetail)
save_img('./result/Adetail.png', Adetail)
save_img('./result/AFinal.png', AFinal)
save_img('./result/mask.png', mask)```