HDR图像高斯双边滤波MATLAB实现

专门针对HDR（高动态范围）图像设计，兼容浮点型数据、保留完整动态范围，避免普通滤波算法对HDR的数值截断问题。同时提供可视化对比和工业级优化建议，适配影视后期、计算摄影、高动态场景去噪等需求。

一、核心原理（适配HDR的关键设计）

HDR图像通常以double/single浮点格式存储，动态范围可达0~10^6甚至更高，和普通LDR图像（0-255 uint8）的处理逻辑有本质区别：

1. 禁止类型转换：滤波全程保持浮点类型，绝对不能转uint8，否则会永久丢失动态范围。
2. 滤波核权重归一化：保证滤波后整体亮度均值不变，仅做局部平滑。
3. 双边滤波适配高动态：值域权重需要匹配HDR的动态范围，避免过强的亮暗差异被误判为边缘。

二、可运行代码

2.1 主脚本hdr_filtering_main.m

%% HDR图像高斯/双边滤波主程序clear;clc;close all;%% ===== 1. 读取HDR图像 =====% 支持格式：.hdr（Radiance）、.exr（OpenEXR，需Image Processing Toolbox）% 内置示例HDR（MATLAB自带，若无可自行替换为自己的HDR路径）hdr_path='memorial.hdr';% MATLAB自带HDR示例，位于 toolbox/images/imdata/if~exist(hdr_path,'file')% 若找不到内置示例，提示用户指定路径[file,path]=uigetfile({'*.hdr;*.exr'},'选择HDR图像');iffile==0,error('未选择HDR文件');endhdr_path=fullfile(path,file);end% 读取HDR（返回double类型，保留完整动态范围）[~,~,ext]=fileparts(hdr_path);ifstrcmpi(ext,'.hdr')hdr_img=hdrread(hdr_path);% Radiance .hdr格式elseifstrcmpi(ext,'.exr')hdr_img=exrread(hdr_path);% OpenEXR格式（需Image Processing Toolbox）endfprintf('HDR图像尺寸: %d×%d×%d | 动态范围: [%.2f, %.2f]\n',...size(hdr_img,1),size(hdr_img,2),size(hdr_img,3),...min(hdr_img(:)),max(hdr_img(:)));%% ===== 2. 滤波参数配置 =====params.gaussian_sigma_s=2.0;% 高斯滤波空间sigma（1~3，值越大模糊越强）params.bilateral_sigma_s=2.0;% 双边滤波空间sigma（1~3）params.bilateral_sigma_r=0.15;% 双边滤波值域sigma（自适应计算，见下文）params.filter_window=9;% 滤波窗口大小（奇数，最大建议11，避免计算过慢）params.tone_map_gamma=2.2;% 色调映射gamma（仅用于显示，不改变原始HDR数据）% 自适应计算双边滤波值域sigma（避免手动调参，适配不同HDR的动态范围）global_std=std(hdr_img(:));params.bilateral_sigma_r=0.1*global_std;% 取全局标准差的10%作为值域sigmafprintf('滤波参数:\n 高斯sigma: %.1f | 双边空间sigma: %.1f | 双边值域sigma: %.3f\n',...params.gaussian_sigma_s,params.bilateral_sigma_s,params.bilateral_sigma_r);%% ===== 3. 执行滤波 =====tic;hdr_gaussian=gaussian_filter_hdr(hdr_img,params.gaussian_sigma_s,params.filter_window);fprintf('高斯滤波完成，用时: %.2f秒\n',toc);tic;hdr_bilateral=bilateral_filter_hdr(hdr_img,params.bilateral_sigma_s,params.bilateral_sigma_r,params.filter_window);fprintf('双边滤波完成，用时: %.2f秒\n',toc);%% ===== 4. 色调映射（仅用于可视化，不改变原始HDR数据）=====% HDR动态范围太大，无法直接显示，需要映射到0-1区间tm_original=reinhard_tone_mapping(hdr_img,params.tone_map_gamma);tm_gaussian=reinhard_tone_mapping(hdr_gaussian,params.tone_map_gamma);tm_bilateral=reinhard_tone_mapping(hdr_bilateral,params.tone_map_gamma);%% ===== 5. 结果可视化对比 =====visualize_hdr_filtering(tm_original,tm_gaussian,tm_bilateral,hdr_img,hdr_gaussian,hdr_bilateral);%% ===== 6. 保存滤波后的HDR图像 =====save_option=questdlg('是否保存滤波后的HDR图像？','保存','是','否','是');ifstrcmp(save_option,'是')[save_file,save_path]=uiputfile({'*.hdr','Radiance HDR (*.hdr)';'*.exr','OpenEXR (*.exr)'},'保存HDR','filtered_hdr.hdr');ifsave_file~=0save_path_full=fullfile(save_path,save_file);[~,~,ext]=fileparts(save_path_full);ifstrcmpi(ext,'.hdr')hdrwrite(hdr_bilateral,save_path_full);% 保存双边滤波结果（更常用）elseifstrcmpi(ext,'.exr')exrwrite(hdr_bilateral,save_path_full);endfprintf('HDR图像已保存到: %s\n',save_path_full);endend

2.2 高斯滤波函数（适配HDR浮点数据）

functionhdr_filtered=gaussian_filter_hdr(hdr_img,sigma_s,window_size)% HDR图像高斯滤波（保持浮点精度，不改变动态范围）% 输入:% hdr_img: HDR图像 (H×W×3 double)% sigma_s: 空间sigma% window_size: 滤波窗口大小（奇数）% 输出:% hdr_filtered: 滤波后的HDR图像 (H×W×3 double)% 生成高斯滤波核（权重和为1，保证亮度不变）half_win=floor(window_size/2);[x,y]=meshgrid(-half_win:half_win,-half_win:half_win);gaussian_kernel=exp(-(x.^2+y.^2)/(2*sigma_s^2));gaussian_kernel=gaussian_kernel/sum(gaussian_kernel(:));% 归一化% 对每个颜色通道独立滤波（HDR的RGB通道互不干扰）hdr_filtered=zeros(size(hdr_img),'like',hdr_img);forc=1:size(hdr_img,3)% imfilter默认支持double输入，输出保持doublehdr_filtered(:,:,c)=imfilter(hdr_img(:,:,c),gaussian_kernel,'replicate');endend

2.3 双边滤波函数（手写，无工具箱依赖）

functionhdr_filtered=bilateral_filter_hdr(hdr_img,sigma_s,sigma_r,window_size)% HDR图像双边滤波（保边平滑，适配高动态范围）% 输入:% hdr_img: HDR图像 (H×W×3 double)% sigma_s: 空间sigma% sigma_r: 值域sigma（需匹配HDR动态范围）% window_size: 滤波窗口大小（奇数，建议≤11，否则极慢）% 输出:% hdr_filtered: 滤波后的HDR图像 (H×W×3 double)[H,W,C]=size(hdr_img);half_win=floor(window_size/2);% 预计算空间权重（所有像素共享，提高效率）[x,y]=meshgrid(-half_win:half_win,-half_win:half_win);spatial_weight=exp(-(x.^2+y.^2)/(2*sigma_s^2));hdr_filtered=zeros(H,W,C,'like',hdr_img);% 对每个通道独立处理forc=1:Cfprintf(' 双边滤波：处理第%d/%d通道...\n',c,C);channel=hdr_img(:,:,c);% 遍历每个像素（嵌套循环较慢，大图像建议用GPU/积分图优化，见下文扩展建议）fori=1:Hforj=1:W% 提取当前像素的邻域i_min=max(1,i-half_win);i_max=min(H,i+half_win);j_min=max(1,j-half_win);j_max=min(W,j+half_win);% 邻域内的像素值neighbor_vals=channel(i_min:i_max,j_min:j_max);% 当前像素值center_val=channel(i,j);% 截取对应的空间权重sp_win=spatial_weight(...(i_min:i_max)-i+half_win+1,...(j_min:j_max)-j+half_win+1);% 计算值域权重：exp(-(邻域值-中心值)²/(2*sigma_r²))range_weight=exp(-(neighbor_vals-center_val).^2/(2*sigma_r^2));% 总权重 = 空间权重 × 值域权重，归一化total_weight=sp_win.*range_weight;total_weight=total_weight/sum(total_weight(:));% 加权求和得到滤波后的值hdr_filtered(i,j,c)=sum(neighbor_vals(:).*total_weight(:));endendendend

2.4 色调映射函数（仅用于显示）

functiontm_img=reinhard_tone_mapping(hdr_img,gamma)% Reinhard全局色调映射（将HDR映射到0-1区间，仅用于显示）% 输入:% hdr_img: HDR图像 (H×W×3 double)% gamma: gamma校正系数（2.2为SRGB标准）% 输出:% tm_img: 色调映射后的LDR图像 (H×W×3 double，范围0-1)% 计算全局白点（取99%分位数，避免过曝）L=0.2126*hdr_img(:,:,1)+0.7152*hdr_img(:,:,2)+0.0722*hdr_img(:,:,3);% 亮度通道white_point=prctile(L(:),99);% Reinhard色调映射公式L_tm=L.*(1+L./(white_point^2))./(1+L);% 亮度归一化L_tm=L_tm/max(L_tm(:));% 保持颜色比例，对每个通道应用相同的亮度缩放scale=L_tm./(L+eps);tm_img=hdr_img.*scale;% Gamma校正tm_img=tm_img.^(1/gamma);% 截断到0-1区间tm_img=max(min(tm_img,1),0);end

2.5 可视化对比函数

functionvisualize_hdr_filtering(tm_original,tm_gaussian,tm_bilateral,hdr_original,hdr_gaussian,hdr_bilateral)figure('Color','w','Position',[1001001600900]);% 1. 整体效果对比subplot(2,3,1);imshow(tm_original);title('原图（色调映射后）');axis image off;subplot(2,3,2);imshow(tm_gaussian);title('高斯滤波后');axis image off;subplot(2,3,3);imshow(tm_bilateral);title('双边滤波后');axis image off;% 2. 局部放大对比（突出边缘保留效果）% 选择高对比度区域（如亮部边缘）zoom_rect=[400,200,150,150];% [x, y, width, height]，可根据图像调整subplot(2,3,4);imshow(imcrop(tm_original,zoom_rect));title('原图局部放大');axis image off;subplot(2,3,5);imshow(imcrop(tm_gaussian,zoom_rect));title('高斯滤波：边缘模糊');axis image off;subplot(2,3,6);imshow(imcrop(tm_bilateral,zoom_rect));title('双边滤波：边缘清晰');axis image off;% 3. 动态范围对比（亮部像素值分布）figure('Color','w','Position',[1001001200400]);subplot(1,3,1);histogram(hdr_original(:),100,'Normalization','pdf');title('原图像素值分布');xlabel('亮度');ylabel('概率密度');grid on;subplot(1,3,2);histogram(hdr_gaussian(:),100,'Normalization','pdf');title('高斯滤波后分布');xlabel('亮度');grid on;subplot(1,3,3);histogram(hdr_bilateral(:),100,'Normalization','pdf');title('双边滤波后分布');xlabel('亮度');grid on;sgtitle('HDR滤波前后动态范围对比（高斯滤波不改变分布形态，仅平滑局部）','FontSize',14,'FontWeight','bold');end

三、运行说明

3.1 直接运行

1. 将所有函数保存为.m文件，放在同一文件夹
2. 运行主脚本hdr_filtering_main.m
3. 程序会自动加载MATLAB自带的memorial.hdr示例，也可手动选择自己的HDR文件

3.2 参数调优建议

3.3 预期效果

HDR图像尺寸: 768×512×3 | 动态范围: [0.00, 1023.98] 滤波参数: 高斯sigma: 2.0 | 双边空间sigma: 2.0 | 双边值域sigma: 12.345 高斯滤波完成，用时: 0.12秒 双边滤波：处理第1/3通道... 双边滤波：处理第2/3通道... 双边滤波：处理第3/3通道... 双边滤波完成，用时: 45.67秒

• 高斯滤波：亮部边缘（如窗户边框）明显模糊
• 双边滤波：亮部边缘清晰，暗部噪声被有效平滑

参考代码 对输入的HDR图像执行高斯滤波或双边滤波www.youwenfan.com/contentcsw/82271.html

四、工业级优化建议

双边滤波的嵌套循环在4K HDR图像上会非常慢，以下是优化方案：

1. GPU加速：将数据转为gpuArray，滤波后在转回CPU：

   hdr_img_gpu=gpuArray(hdr_img);% 滤波函数在GPU上运行（需修改代码适配GPU数组）hdr_filtered_gpu=bilateral_filter_hdr(hdr_img_gpu,sigma_s,sigma_r,window_size);hdr_filtered=gather(hdr_filtered_gpu);

2. 积分图优化：预计算积分图和平方积分图，将双边滤波的时间复杂度从O(N×k²)降到O(N)，适合大图像。

3. 使用内置函数：若安装了Computer Vision Toolbox，可直接用imbilatfilt，速度比手写快10倍以上：

   % 内置双边滤波（支持HDR）hdr_bilateral=imbilatfilt(hdr_img,'SpatialSigma',2,'RangeSigma',params.bilateral_sigma_r);

4. 降采样滤波：先将HDR降采样到1/2或1/4尺寸，滤波后再上采样，速度提升4~16倍，精度损失很小。

五、应用场景扩展

1. HDR去噪：双边滤波是HDR图像去噪的首选，比高斯滤波保留更多细节，适合夜景、高对比度场景。
2. HDR tone mapping预处理：滤波后可以减少色调映射后的伪影，让亮部过渡更自然。
3. HDR纹理平滑：比如游戏HDR贴图的纹理平滑，保留硬边缘的同时平滑表面纹理。
4. 高动态场景目标检测：滤波后降低噪声，提升检测算法对亮部和暗部目标的识别率。