HDR图像生成算法详解
HDR技术概述
高动态范围（HDR）图像生成是通过合成多张不同曝光度的图像，来捕捉超出传统显示设备范围的亮度信息。以下是主要的HDR生成算法及其MATLAB实现。
HDR成像基本原理
动态范围定义
动态范围 = 最大可记录亮度 / 最小可记录亮度
- LDR图像：~100:1
- HDR图像：可达100,000:1
HDR图像生成核心算法
1. 多曝光图像合成算法
function hdr_image = generate_hdr_multi_exposure(images, exposure_times)
% 基于多曝光图像的HDR合成
% 输入:
% images - 不同曝光的图像序列 {N×1} cell array
% exposure_times - 对应的曝光时间向量
% 输出:
% hdr_image - 合成的HDR图像
num_images = length(images);
[height, width, channels] = size(images{1});
% 初始化权重矩阵和HDR图像
hdr_image = zeros(height, width, channels);
total_weight = zeros(height, width, channels);
% 相机响应函数估计（可选，如果未提供则使用默认权重）
if ~exist('camera_response', 'var')
% 使用三角权重函数
weight_function = @(z) max(0, 1 - abs(2*z - 1)); % z在[0,1]范围内
end
fprintf('开始合成HDR图像...\n');
for i = 1:num_images
fprintf('处理第%d张曝光图像...\n', i);
% 将图像转换为double类型并归一化
img = im2double(images{i});
% 计算权重
weight_map = zeros(size(img));
for c = 1:channels
weight_map(:, :, c) = weight_function(img(:, :, c));
end
% 避免过曝和欠曝区域
saturation_threshold = 0.95;
underexposure_threshold = 0.05;
% 降低过曝和欠曝区域的权重
over_exposed = img > saturation_threshold;
under_exposed = img < underexposure_threshold;
weight_map(over_exposed) = 0.01;
weight_map(under_exposed) = 0.01;
% 累加到HDR图像
hdr_image = hdr_image + weight_map .* img / exposure_times(i);
total_weight = total_weight + weight_map;
end
% 避免除零错误
total_weight(total_weight == 0) = 1e-8;
% 归一化得到最终HDR图像
hdr_image = hdr_image ./ total_weight;
fprintf('HDR图像合成完成\n');
end
2. 相机响应函数估计
function [g, lnE] = estimate_camera_response(images, exposure_times, lambda)
% 使用Debevec方法估计相机响应函数
% 输入:
% images - 不同曝光的图像序列
% exposure_times - 曝光时间
% lambda - 平滑度参数
% 输出:
% g - 估计的相机响应函数
% lnE - 场景辐照度的对数
num_images = length(images);
[height, width, channels] = size(images{1});
num_pixels = height * width;
% 采样像素点（避免过曝和欠曝）
sample_indices = sample_pixels(images, 1000);
num_samples = length(sample_indices);
% 构建线性系统
n = 256; % 量化级别 (8-bit)
A = zeros(num_samples * num_images + n + 1, num_samples + n);
b = zeros(size(A, 1), 1);
k = 1;
% 添加数据项约束
for i = 1:num_images
img = images{i};
for j = 1:num_samples
[row, col] = ind2sub([height, width], sample_indices(j));
z = round(img(row, col, 1) * (n-1)) + 1; % 使用单个通道
w_z = weighting_function(z/n);
A(k, z) = w_z;
A(k, n + j) = -w_z;
b(k) = w_z * log(exposure_times(i));
k = k + 1;
end
end
% 添加平滑度约束
for z = 2:n-1
w_z = weighting_function(z/n);
A(k, z-1) = lambda * w_z;
A(k, z) = -2 * lambda * w_z;
A(k, z+1) = lambda * w_z;
k = k + 1;
end
% 添加锚点约束（设置g(128)=0）
A(k, round(n/2)) = 1;
k = k + 1;
% 求解线性系统
x = A \ b;
g = x(1:n);
lnE = x(n+1:end);
end
function indices = sample_pixels(images, num_samples)
% 采样像素点，避免过曝和欠曝区域
[height, width, ~] = size(images{1});
% 计算像素的可靠性分数
reliability = ones(height, width);
num_images = length(images);
for i = 1:num_images
img = im2double(images{i});
% 中间曝光度最可靠
mid_tone = 0.5;
reliability = reliability .* (1 - abs(img - mid_tone));
end
% 随机采样，权重为可靠性
reliability_flat = reliability(:);
probabilities = reliability_flat / sum(reliability_flat);
indices = randsample(height * width, num_samples, true, probabilities);
end
function w = weighting_function(z)
% 权重函数：中间调高权重，过曝和欠曝低权重
z = max(0, min(1, z));
if z <= 0.5
w = z;
else
w = 1 - z;
end
end
3. 色调映射算法
function ldr_image = tone_mapping(hdr_image, method, varargin)
% 色调映射：将HDR图像映射到LDR显示范围
% 输入:
% hdr_image - HDR输入图像
% method - 色调映射方法
% 输出:
% ldr_image - 映射后的LDR图像
switch lower(method)
case 'reinhard'
ldr_image = reinhard_tone_mapping(hdr_image, varargin{:});
case 'durand'
ldr_image = durand_tone_mapping(hdr_image, varargin{:});
case 'drago'
ldr_image = drago_tone_mapping(hdr_image, varargin{:});
case 'gamma'
ldr_image = gamma_tone_mapping(hdr_image, varargin{:});
otherwise
error('不支持的色调映射方法');
end
% 确保输出在[0,1]范围内
ldr_image = max(0, min(1, ldr_image));
end
function ldr_image = reinhard_tone_mapping(hdr_image, L_white, a)
% Reinhard全局色调映射算法
% 参考: Photographic Tone Reproduction for Digital Images
if nargin < 2
L_white = max(hdr_image(:)); % 最大亮度值
end
if nargin < 3
a = 0.18; % 关键值
end
% 转换为亮度
L = 0.2126 * hdr_image(:,:,1) + 0.7152 * hdr_image(:,:,2) + 0.0722 * hdr_image(:,:,3);
% 计算场景平均亮度（对数平均）
delta = 1e-6;
L_avg = exp(mean(log(L + delta)));
% 尺度缩放
L_scaled = (a / L_avg) * L;
% Reinhard色调映射算子
L_display = L_scaled .* (1 + L_scaled / (L_white^2)) ./ (1 + L_scaled);
% 恢复颜色
ldr_image = zeros(size(hdr_image));
for c = 1:3
ldr_image(:,:,c) = (hdr_image(:,:,c) ./ L) .* L_display;
end
end
function ldr_image = durand_tone_mapping(hdr_image, base_contrast, detail_contrast)
% Durand快速双边滤波色调映射
% 分离基础层和细节层
if nargin < 2
base_contrast = 5.0;
end
if nargin < 3
detail_contrast = 1.5;
end
% 转换为亮度
L = 0.2126 * hdr_image(:,:,1) + 0.7152 * hdr_image(:,:,2) + 0.0722 * hdr_image(:,:,3);
% 对数变换
log_L = log10(L + 1e-6);
% 双边滤波分离基础层和细节层
base_layer = bilateral_filter(log_L, 0.4, 0.1);
detail_layer = log_L - base_layer;
% 压缩基础层的动态范围
log_base_compressed = (base_layer - max(base_layer(:))) * base_contrast + max(base_layer(:));
% 增强细节层
detail_enhanced = detail_layer * detail_contrast;
% 重建对数亮度
log_output = log_base_compressed + detail_enhanced;
% 指数变换回线性空间
L_output = 10.^(log_output);
% 恢复颜色
ldr_image = zeros(size(hdr_image));
for c = 1:3
ldr_image(:,:,c) = (hdr_image(:,:,c) ./ L) .* L_output;
end
end
function filtered = bilateral_filter(image, sigma_spatial, sigma_range)
% 双边滤波实现
[height, width] = size(image);
filtered = zeros(size(image));
% 空间高斯核
spatial_size = ceil(3 * sigma_spatial);
[X, Y] = meshgrid(-spatial_size:spatial_size, -spatial_size:spatial_size);
spatial_kernel = exp(-(X.^2 + Y.^2) / (2 * sigma_spatial^2));
for i = 1:height
for j = 1:width
% 计算局部窗口
i_min = max(1, i - spatial_size);
i_max = min(height, i + spatial_size);
j_min = max(1, j - spatial_size);
j_max = min(width, j + spatial_size);
% 提取局部区域
local_region = image(i_min:i_max, j_min:j_max);
% 范围核
range_kernel = exp(-(local_region - image(i,j)).^2 / (2 * sigma_range^2));
% 组合核
combined_kernel = spatial_kernel(...
(i_min:i_max)-i+spatial_size+1, ...
(j_min:j_max)-j+spatial_size+1) .* range_kernel;
% 归一化并滤波
combined_kernel = combined_kernel / sum(combined_kernel(:));
filtered(i,j) = sum(sum(local_region .* combined_kernel));
end
end
end
4. 完整的HDR生成系统
function hdr_pipeline_demo()
% 完整的HDR图像生成演示系统
fprintf('=== HDR图像生成管道演示 ===\n\n');
% 1. 生成或加载多曝光图像序列
fprintf('1. 准备多曝光图像序列...\n');
[images, exposure_times] = generate_test_exposure_sequence();
% 2. 估计相机响应函数
fprintf('2. 估计相机响应函数...\n');
[g, lnE] = estimate_camera_response(images, exposure_times, 50);
% 3. 合成HDR图像
fprintf('3. 合成HDR图像...\n');
hdr_image = generate_hdr_multi_exposure(images, exposure_times);
% 4. 色调映射
fprintf('4. 色调映射处理...\n');
ldr_results = apply_tone_mapping_methods(hdr_image);
% 5. 质量评估和结果显示
fprintf('5. 评估和显示结果...\n');
evaluate_and_display_results(images, hdr_image, ldr_results, exposure_times, g);
fprintf('HDR管道处理完成！\n');
end
function [images, exposure_times] = generate_test_exposure_sequence()
% 生成测试多曝光图像序列
% 如果没有真实数据，使用合成图像演示
fprintf('生成合成多曝光图像序列...\n');
% 创建合成HDR场景（高动态范围）
[X, Y] = meshgrid(1:512, 1:512);
center_x = 256; center_y = 256;
% 创建不同亮度的区域
hdr_scene = zeros(512, 512, 3);
% 暗区域
dark_region = ((X - 100).^2 + (Y - 100).^2) < 50^2;
hdr_scene(:,:,1) = hdr_scene(:,:,1) + 0.1 * dark_region; % 很暗的红色
% 中等亮度区域
mid_region = ((X - 256).^2 + (Y - 256).^2) < 100^2;
hdr_scene(:,:,2) = hdr_scene(:,:,2) + 1.0 * mid_region; % 中等绿色
% 高亮区域
bright_region = ((X - 400).^2 + (Y - 400).^2) < 30^2;
hdr_scene(:,:,3) = hdr_scene(:,:,3) + 10.0 * bright_region; % 很亮的蓝色
% 背景纹理
background = 0.3 + 0.1 * sin(X/20) .* cos(Y/15);
hdr_scene = hdr_scene + repmat(background, [1, 1, 3]);
% 生成不同曝光的LDR图像
exposure_times = [0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8];
num_exposures = length(exposure_times);
images = cell(num_exposures, 1);
for i = 1:num_exposures
% 模拟相机响应（简单的gamma曲线）
exposed = min(1, hdr_scene * exposure_times(i));
exposed = exposed .^ (1/2.2); % gamma校正
% 添加噪声
exposed = exposed + 0.01 * randn(size(exposed));
images{i} = max(0, min(1, exposed));
fprintf(' 曝光时间 %.4f: 动态范围 [%.3f, %.3f]\n', ...
exposure_times(i), min(exposed(:)), max(exposed(:)));
end
end
function ldr_results = apply_tone_mapping_methods(hdr_image)
% 应用不同的色调映射方法
fprintf('应用不同的色调映射算法...\n');
ldr_results = struct();
% Reinhard方法
fprintf(' - Reinhard色调映射...\n');
ldr_results.reinhard = tone_mapping(hdr_image, 'reinhard');
% Durand方法
fprintf(' - Durand色调映射...\n');
ldr_results.durand = tone_mapping(hdr_image, 'durand');
% Gamma方法（简单对比）
fprintf(' - Gamma色调映射...\n');
ldr_results.gamma = tone_mapping(hdr_image, 'gamma', 2.2, 0.5);
% 直方图均衡化方法
fprintf(' - 直方图均衡化...\n');
ldr_results.histeq = histogram_equalization_tone_mapping(hdr_image);
end
function ldr_image = gamma_tone_mapping(hdr_image, gamma_value, scale)
% 简单的gamma色调映射
if nargin < 2
gamma_value = 2.2;
end
if nargin < 3
scale = 1.0;
end
% 归一化并应用gamma校正
ldr_image = scale * hdr_image / max(hdr_image(:));
ldr_image = ldr_image .^ (1/gamma_value);
end
function ldr_image = histogram_equalization_tone_mapping(hdr_image)
% 基于直方图均衡化的色调映射
ldr_image = zeros(size(hdr_image));
for c = 1:3
channel = hdr_image(:,:,c);
% 对数变换压缩动态范围
log_channel = log10(channel + 1e-6);
log_channel = (log_channel - min(log_channel(:))) / ...
(max(log_channel(:)) - min(log_channel(:)));
% 直方图均衡化
ldr_image(:,:,c) = histeq(log_channel);
end
end
5. 质量评估和可视化
function evaluate_and_display_results(images, hdr_image, ldr_results, exposure_times, camera_response)
% 评估HDR生成质量并显示结果
fprintf('评估结果质量...\n');
figure('Position', [100, 100, 1400, 1000]);
% 1. 显示多曝光输入序列
num_inputs = min(6, length(images));
for i = 1:num_inputs
subplot(4, num_inputs, i);
imshow(images{i});
title(sprintf('曝光: %.4f', exposure_times(i)));
end
% 2. 显示相机响应函数
subplot(4, 3, 7);
if exist('camera_response', 'var') && ~isempty(camera_response)
plot(0:255, camera_response, 'b-', 'LineWidth', 2);
title('估计的相机响应函数');
xlabel('像素值');
ylabel('log曝光量');
grid on;
end
% 3. 显示HDR图像统计
subplot(4, 3, 8);
hdr_luminance = 0.2126 * hdr_image(:,:,1) + 0.7152 * hdr_image(:,:,2) + 0.0722 * hdr_image(:,:,3);
hist(log10(hdr_luminance(:) + 1e-6), 50);
title('HDR亮度分布（对数）');
xlabel('log10(亮度)');
ylabel('频率');
grid on;
% 4. 显示动态范围
subplot(4, 3, 9);
dr_info = calculate_dynamic_range(hdr_image);
bar(dr_info.range_per_channel);
title('各通道动态范围');
xlabel('颜色通道');
ylabel('动态范围 (log10)');
grid on;
% 5. 显示不同色调映射结果
methods = fieldnames(ldr_results);
for i = 1:length(methods)
subplot(4, 3, 9 + i);
imshow(ldr_results.(methods{i}));
title([methods{i} ' 色调映射']);
end
% 6. 质量评估指标
fprintf('\n=== 质量评估结果 ===\n');
evaluate_tone_mapping_quality(hdr_image, ldr_results);
end
function dr_info = calculate_dynamic_range(hdr_image)
% 计算HDR图像的动态范围
dr_info = struct();
for c = 1:3
channel = hdr_image(:,:,c);
valid_pixels = channel > 0;
if sum(valid_pixels(:)) > 0
min_val = min(channel(valid_pixels));
max_val = max(channel(valid_pixels));
dr_info.range_per_channel(c) = log10(max_val / min_val);
else
dr_info.range_per_channel(c) = 0;
end
end
dr_info.average_dynamic_range = mean(dr_info.range_per_channel);
fprintf('平均动态范围: %.2f log10单位\n', dr_info.average_dynamic_range);
end
function evaluate_tone_mapping_quality(hdr_reference, ldr_results)
% 评估色调映射结果的质量
methods = fieldnames(ldr_results);
fprintf('\n%-15s %-12s %-12s %-12s\n', ...
'方法', '对比度', '信息熵', '色彩饱和度');
fprintf('%-15s %-12s %-12s %-12s\n', ...
'------', '------', '------', '----------');
for i = 1:length(methods)
ldr_image = ldr_results.(methods{i});
% 对比度（标准差）
contrast = std(ldr_image(:));
% 信息熵
entropy_val = calculate_image_entropy(ldr_image);
% 色彩饱和度
color_saturation = calculate_color_saturation(ldr_image);
fprintf('%-15s %-12.4f %-12.4f %-12.4f\n', ...
methods{i}, contrast, entropy_val, color_saturation);
end
end
function entropy_val = calculate_image_entropy(image)
% 计算图像的信息熵
if size(image, 3) == 3
image = rgb2gray(image);
end
% 计算直方图
[counts, ~] = histcounts(image(:), 256);
probabilities = counts / sum(counts);
probabilities = probabilities(probabilities > 0);
entropy_val = -sum(probabilities .* log2(probabilities));
end
function saturation = calculate_color_saturation(rgb_image)
% 计算图像的平均色彩饱和度
hsv_image = rgb2hsv(rgb_image);
saturation = mean(hsv_image(:,:,2), 'all');
end
6. 高级HDR技术
function advanced_hdr_techniques_demo()
% 高级HDR技术演示
fprintf('=== 高级HDR技术演示 ===\n\n');
% 1. 曝光融合（不生成显式HDR）
fprintf('1. 曝光融合技术...\n');
exposure_fusion_result = exposure_fusion(images, exposure_times);
% 2. 抗鬼影处理
fprintf('2. 抗鬼影处理...\n');
ghost_free_hdr = ghost_removal_hdr(images, exposure_times);
% 3. 单图像HDR重建
fprintf('3. 单图像HDR重建...\n');
single_image_hdr = single_image_hdr_reconstruction(images{ceil(end/2)});
% 显示比较结果
figure('Position', [100, 100, 1200, 800]);
subplot(2, 2, 1);
imshow(exposure_fusion_result);
title('曝光融合结果');
subplot(2, 2, 2);
imshow(tone_mapping(ghost_free_hdr, 'reinhard'));
title('抗鬼影HDR');
subplot(2, 2, 3);
imshow(tone_mapping(single_image_hdr, 'reinhard'));
title('单图像HDR重建');
subplot(2, 2, 4);
% 显示算法比较
plot_algorithm_comparison();
end
function fused_image = exposure_fusion(images, exposure_times)
% 曝光融合：直接在LDR空间融合多曝光图像
num_images = length(images);
weight_maps = cell(num_images, 1);
% 为每张图像计算质量权重图
for i = 1:num_images
img = im2double(images{i});
% 计算对比度（使用拉普拉斯滤波）
contrast = abs(imfilter(rgb2gray(img), fspecial('laplacian')));
% 计算饱和度（颜色标准差）
saturation = std(img, [], 3);
% 计算曝光良好度（高斯函数，以0.5为中心）
well_exposed = exp(-(rgb2gray(img) - 0.5).^2 / (2 * 0.2^2));
% 组合权重图
weight_maps{i} = contrast .* saturation .* well_exposed + 1e-6;
end
% 归一化权重图
total_weight = zeros(size(weight_maps{1}));
for i = 1:num_images
total_weight = total_weight + weight_maps{i};
end
% 多分辨率融合（金字塔融合）
fused_image = pyramid_fusion(images, weight_maps);
end
function fused_image = pyramid_fusion(images, weight_maps)
% 使用拉普拉斯金字塔进行多分辨率融合
num_levels = 5;
num_images = length(images);
% 为每张图像构建拉普拉斯金字塔和权重高斯金字塔
laplacian_pyramids = cell(num_images, 1);
gaussian_pyramids = cell(num_images, 1);
for i = 1:num_images
laplacian_pyramids{i} = laplacian_pyramid(images{i}, num_levels);
gaussian_pyramids{i} = gaussian_pyramid(weight_maps{i}, num_levels);
end
% 融合金字塔
fused_pyramid = laplacian_pyramids{1};
for l = 1:num_levels
level_weights = zeros([size(gaussian_pyramids{1}{l}), num_images]);
for i = 1:num_images
level_weights(:,:,i) = gaussian_pyramids{i}{l};
end
% 归一化权重
level_weights = level_weights ./ sum(level_weights, 3);
% 加权融合
if l <= num_levels
fused_level = zeros(size(laplacian_pyramids{1}{l}));
for i = 1:num_images
weight_expanded = repmat(level_weights(:,:,i), [1, 1, size(images{i}, 3)]);
fused_level = fused_level + weight_expanded .* laplacian_pyramids{i}{l};
end
fused_pyramid{l} = fused_level;
end
end
% 重建融合图像
fused_image = reconstruct_laplacian_pyramid(fused_pyramid);
end
参考代码 HDR图像生成的算法 www.youwenfan.com/contentcnj/66127.html
使用方法
- 基本HDR生成：
% 运行完整的HDR管道
hdr_pipeline_demo();
- 自定义图像序列：
% 加载自己的多曝光图像
image_files = {'exp1.jpg', 'exp2.jpg', 'exp3.jpg'};
exposure_times = [1/30, 1/15, 1/8]; % 对应的曝光时间
images = cell(length(image_files), 1);
for i = 1:length(image_files)
images{i} = imread(image_files{i});
end
% 生成HDR
hdr_image = generate_hdr_multi_exposure(images, exposure_times);
- 高级功能：
% 运行高级HDR技术演示
advanced_hdr_techniques_demo();
浙公网安备 33010602011771号