【HDR】曝光融合（Exposure Fusion）
文章目录
-
0 前言
-
1 算法细节
-
-
1.1 Naive
-
- 1.1.1 主要思想
- 1.1.2 权重计算
- 1.1.3 融合
-
1.2 Multi-resolution
-
-
2 实验
-
3 参考
0 前言
当时提出的一种技术，在未采用现代先进的计算能力时显得较为繁琐。生成高动态范围图像的过程通常分为两个主要步骤：首先将多幅不同曝光度的低动态范围图像融合生成高动态范围图像（例如Debevec提出的加权融合方法），其中通常采用12位或16位深度；其次通过对高动态范围图像进行tone mapping压缩处理来适应低dynamic range显示设备的需求（例如Durand提出的基于双边滤波的tone mapping算法），结果通常压缩到8位深度。
该算法的优势体现在无需对相机响应曲线进行标定点校，并省去了tone mapping步骤，在生成用于显示用途下的HDR图像方面表现出色。
1 算法细节
1.1 Naive
1.1.1 主要思想
在多曝光图像序列中，在每一幅图像中都应选取其像素重要性程度较高的区域进行合成处理；例如，在曝光时间较长的图像中往往能够获得细节较为丰富且信噪比较高的区域，则这些区域即为重要的部分；显然地，则需建立一种评估机制来量化每一幅图像内各像素的重要性程度；随后通过对每一幅图片内各像素的重要性指标赋予权重值并结合加权融合方法最终完成高动态范围（HDR）的重建过程
1.1.2 权重计算
从对比度、饱和度和亮度三个维度对像素的价值进行评估：
-
对比度指标
这里的对比度指标指的是图像的梯度。具体而言，在处理边缘、纹理等重要信息时赋予了较大的权重。这一过程通过将图像转换为灰度图并进行拉普斯滤波处理后得到的结果取绝对值得到的结果作为对比度指标 C(I_k)。 -
饱和度
若RGB三个通道间的差异较大，则该区域可被视作饱和度较高的区域；反地，则过曝或欠曝区域的RGB三个通道值趋于一致，其饱和度较低。因此可以采用RGB三个通道的标准差作为衡量该区域饱和度的指标。
S_{k,i,j}=\sqrt{\frac{(I_{k,i,j}^{B}-M_{k,i,j})^2+(I_{k,i,j}^{G}-M_{k,i,j})^2+(I_{k,i,j}^{R}-M_{k,i,j})^2}{3}}
M_{k,i,j}=\frac{I_{k,i,j}^{B}+I_{k,i,j}^{G}+I_{k,i,j}^{R}}{3}
对于归一化后取值范围限定在0到1之间的图像，在其像素中具有明暗程度处于中等水平（约在0.5附近）的部分应赋予较高的权重；而亮度接近极端值（接近于0或1）的部分则应赋予较低的权重。具体而言，在红、绿、蓝三个通道上每个像素与其对应的权重之间均遵循均值为0.5的标准差为σ的高斯分布关系：
E_{k,i,j}=e^{-\frac{(I_{k,i,j}^{B}-0.5)^2}{2\sigma^2}} \cdot e^{-\frac{(I_{k,i,j}^{G}-0.5)^2}{2\sigma^2}} \cdot e^{-\frac{(I_{k,i,j}^{R}-0.5)^2}{2\sigma^2}}
基于以上三个指标的获取后（即通过这三个指标的数据），可以计算出每个像素所对应的权重值）。其中，默认情况下w_c=w_s=w_e=1；此外，在图像数量维度上对权重进行归一化处理以确保同一位置的权重和为1）。数学表达式如下所示：\hat{W}_{{ij,k}}=\frac{{W}_{{ij,k}}}{\sum_{k^{\prime}=1}^{N}{W}_{{ij,k^{\prime}}}}
1.1.3 融合
通过预设的权重系数对原始图像执行加权叠加运算，从而获得融合后的图像：
H_{ij} = \sum_{k=1}^{N} \hat{W}_{{ij,k}} \cdot I_{{ij,k}}
在进行这种粗糙融合时会出现的问题是在权重突变区域。因为各图像的曝光时间不一致,其绝对强度水平也存在差异,这会导致融合后的灰度变化过于剧烈,并且会使得图像呈现大量黑色和白色斑块的状态,这一现象与噪声分布非常相似。
在权重突变的区域出现处理困难的情况下, 可以适当减小过渡程度, 这自然联想到使用高斯滤波器进行模糊处理. 作者将权重图经过高斯滤波后再次融合合成, 尽管消除了斑点噪声到一定程度, 但合成图像在边界区域可能出现模糊现象.
光晕现象源于边缘区域权重的平滑处理。因此，在图像处理中，建议采用保边滤波来取代高斯滤波。
1.2 Multi-resolution
由于简单版本的融合方式不仅未能消除黑白斑点问题，并且还带来了光晕等新问题。因此，在此背景下，作者建议采用了拉普拉斯金字塔融合的方法来进行图像合并处理
具体而言，在不同曝光条件下的原始图像中提取出拉普拉斯金字塔的同时，在对应的权重图中也提取出了高斯金字塔。随后，在各个尺度上进行融合处理，并将融合后的结果作为新的拉普拉斯细节图输出。最后从这个细节图的顶层开始进行反向映射操作，在每个层次上叠加相应的细节信息，并逐步递归至最高分辨率层以上的位置处完成最终结果的合成过程。（需要注意的是，在拉普拉斯金字塔结构中其顶层与原始图像的高斯金字塔顶层是相同的）
在处理图像特征时，我们采取了区分处理策略。具体而言，在这一阶段的处理中，我们采用了高斯滤波和平移下采样的方法来生成权重图。其变化主要集中在较大的纹理区域。基于此，在构建拉普拉斯金字塔时，高频细节被优先考虑。这种分离策略使得在处理尖锐过渡区域时不会干扰到整体的平滑区域特性。
2 实验
import os
import sys
import glob
import numpy as np
import cv2
import argparse
def show_image(message, src):
cv2.namedWindow(message, 0)
cv2.imshow(message, src)
cv2.waitKey(0)
cv2.destroyAllWindows()
def gauss_curve(src, mean, sigma):
dst = np.exp(-(src - mean)**2 / (2 * sigma**2))
return dst
class ExposureFusion(object):
def __init__(self, sequence, best_exposedness=0.5, sigma=0.2, eps=1e-12, exponents=(1.0, 1.0, 1.0), layers=7):
self.sequence = sequence # [N, H, W, 3], (0..1), float32
self.img_num = sequence.shape[0]
self.best_exposedness = best_exposedness
self.sigma = sigma
self.eps = eps
self.exponents = exponents
self.layers = layers
@staticmethod
def cal_contrast(src):
gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)
laplace_kernel = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], dtype=np.float32)
contrast = cv2.filter2D(gray, -1, laplace_kernel, borderType=cv2.BORDER_REPLICATE)
return np.abs(contrast)
@staticmethod
def cal_saturation(src):
mean = np.mean(src, axis=-1)
channels = [(src[:, :, c] - mean)**2 for c in range(3)]
saturation = np.sqrt(np.mean(channels, axis=0))
return saturation
@staticmethod
def cal_exposedness(src, best_exposedness, sigma):
exposedness = [gauss_curve(src[:, :, c], best_exposedness, sigma) for c in range(3)]
exposedness = np.prod(exposedness, axis=0)
return exposedness
def cal_weight_map(self):
weights = []
for idx in range(self.sequence.shape[0]):
contrast = self.cal_contrast(self.sequence[idx])
saturation = self.cal_saturation(self.sequence[idx])
exposedness = self.cal_exposedness(self.sequence[idx], self.best_exposedness, self.sigma)
weight = np.power(contrast, self.exponents[0]) * np.power(saturation, self.exponents[1]) * np.power(exposedness, self.exponents[2])
# Gauss Blur
# weight = cv2.GaussianBlur(weight, (21, 21), 2.1)
weights.append(weight)
weights = np.stack(weights, 0) + self.eps
# normalize
weights = weights / np.expand_dims(np.sum(weights, axis=0), axis=0)
return weights
def naive_fusion(self):
weights = self.cal_weight_map() # [N, H, W]
weights = np.stack([weights, weights, weights], axis=-1) # [N, H, W, 3]
naive_fusion = np.sum(weights * self.sequence * 255, axis=0)
naive_fusion = np.clip(naive_fusion, 0, 255).astype(np.uint8)
return naive_fusion
def build_gaussian_pyramid(self, high_res):
gaussian_pyramid = [high_res]
for idx in range(1, self.layers):
gaussian_pyramid.append(cv2.GaussianBlur(gaussian_pyramid[-1], (5, 5), 0.83)[::2, ::2])
return gaussian_pyramid
def build_laplace_pyramid(self, gaussian_pyramid):
laplace_pyramid = [gaussian_pyramid[-1]]
for idx in range(1, self.layers):
size = (gaussian_pyramid[self.layers - idx - 1].shape[1], gaussian_pyramid[self.layers - idx - 1].shape[0])
upsampled = cv2.resize(gaussian_pyramid[self.layers - idx], size, interpolation=cv2.INTER_LINEAR)
laplace_pyramid.append(gaussian_pyramid[self.layers - idx - 1] - upsampled)
laplace_pyramid.reverse()
return laplace_pyramid
def multi_resolution_fusion(self):
weights = self.cal_weight_map() # [N, H, W]
weights = np.stack([weights, weights, weights], axis=-1) # [N, H, W, 3]
image_gaussian_pyramid = [self.build_gaussian_pyramid(self.sequence[i] * 255) for i in range(self.img_num)]
image_laplace_pyramid = [self.build_laplace_pyramid(image_gaussian_pyramid[i]) for i in range(self.img_num)]
weights_gaussian_pyramid = [self.build_gaussian_pyramid(weights[i]) for i in range(self.img_num)]
fused_laplace_pyramid = [np.sum([image_laplace_pyramid[n][l] *
weights_gaussian_pyramid[n][l] for n in range(self.img_num)], axis=0) for l in range(self.layers)]
result = fused_laplace_pyramid[-1]
for k in range(1, self.layers):
size = (fused_laplace_pyramid[self.layers - k - 1].shape[1], fused_laplace_pyramid[self.layers - k - 1].shape[0])
upsampled = cv2.resize(result, size, interpolation=cv2.INTER_LINEAR)
result = upsampled + fused_laplace_pyramid[self.layers - k - 1]
result = np.clip(result, 0, 255).astype(np.uint8)
return result
if __name__ == '__main__':
root_path = sys.argv[1]
sequence_path = [os.path.join(root_path, fname) for fname in os.listdir(root_path)]
sequence = np.stack([cv2.imread(path) for path in sequence_path], axis=0)
mef = ExposureFusion(sequence.astype(np.float32) / 255.0)
naive_fusion_result = mef.naive_fusion()
multi_res_fusion = mef.multi_resolution_fusion()
show_image('muti-resolution', multi_res_fusion)
代码解读3 参考
Mertens T等. 暖光融合技术[C]//第十五届亚洲计算机图形及其应用会议（PG'07）. 北京: IEEE计算机协会出版社, 2007: 382-390.