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

专门针对HDR高动态范围图像设计兼容浮点型数据、保留完整动态范围避免普通滤波算法对HDR的数值截断问题。同时提供可视化对比和工业级优化建议适配影视后期、计算摄影、高动态场景去噪等需求。一、核心原理适配HDR的关键设计HDR图像通常以double/single浮点格式存储动态范围可达0~10^6甚至更高和普通LDR图像0-255 uint8的处理逻辑有本质区别禁止类型转换滤波全程保持浮点类型绝对不能转uint8否则会永久丢失动态范围。滤波核权重归一化保证滤波后整体亮度均值不变仅做局部平滑。双边滤波适配高动态值域权重需要匹配HDR的动态范围避免过强的亮暗差异被误判为边缘。滤波类型适用场景HDR下的优势/劣势高斯滤波全局模糊、预平滑、降噪速度快但会糊掉高对比度边缘如亮部文字、灯箱边缘双边滤波保边去噪、纹理平滑、HDR tone mapping前的预处理同时考虑空间位置和像素值差异完美保留高对比度边缘速度较慢二、可运行代码2.1 主脚本hdr_filtering_main.m%% HDR图像高斯/双边滤波主程序clear;clc;close all;%% 1. 读取HDR图像 % 支持格式.hdrRadiance、.exrOpenEXR需Image Processing Toolbox% 内置示例HDRMATLAB自带若无可自行替换为自己的HDR路径hdr_pathmemorial.hdr;% MATLAB自带HDR示例位于 toolbox/images/imdata/if~exist(hdr_path,file)% 若找不到内置示例提示用户指定路径[file,path]uigetfile({*.hdr;*.exr},选择HDR图像);iffile0,error(未选择HDR文件);endhdr_pathfullfile(path,file);end% 读取HDR返回double类型保留完整动态范围[~,~,ext]fileparts(hdr_path);ifstrcmpi(ext,.hdr)hdr_imghdrread(hdr_path);% Radiance .hdr格式elseifstrcmpi(ext,.exr)hdr_imgexrread(hdr_path);% OpenEXR格式需Image Processing Toolboxendfprintf(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_s2.0;% 高斯滤波空间sigma1~3值越大模糊越强params.bilateral_sigma_s2.0;% 双边滤波空间sigma1~3params.bilateral_sigma_r0.15;% 双边滤波值域sigma自适应计算见下文params.filter_window9;% 滤波窗口大小奇数最大建议11避免计算过慢params.tone_map_gamma2.2;% 色调映射gamma仅用于显示不改变原始HDR数据% 自适应计算双边滤波值域sigma避免手动调参适配不同HDR的动态范围global_stdstd(hdr_img(:));params.bilateral_sigma_r0.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_gaussiangaussian_filter_hdr(hdr_img,params.gaussian_sigma_s,params.filter_window);fprintf(高斯滤波完成用时: %.2f秒\n,toc);tic;hdr_bilateralbilateral_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_originalreinhard_tone_mapping(hdr_img,params.tone_map_gamma);tm_gaussianreinhard_tone_mapping(hdr_gaussian,params.tone_map_gamma);tm_bilateralreinhard_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_optionquestdlg(是否保存滤波后的HDR图像,保存,是,否,是);ifstrcmp(save_option,是)[save_file,save_path]uiputfile({*.hdr,Radiance HDR (*.hdr);*.exr,OpenEXR (*.exr)},保存HDR,filtered_hdr.hdr);ifsave_file~0save_path_fullfullfile(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);endend2.2 高斯滤波函数适配HDR浮点数据functionhdr_filteredgaussian_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_winfloor(window_size/2);[x,y]meshgrid(-half_win:half_win,-half_win:half_win);gaussian_kernelexp(-(x.^2y.^2)/(2*sigma_s^2));gaussian_kernelgaussian_kernel/sum(gaussian_kernel(:));% 归一化% 对每个颜色通道独立滤波HDR的RGB通道互不干扰hdr_filteredzeros(size(hdr_img),like,hdr_img);forc1:size(hdr_img,3)% imfilter默认支持double输入输出保持doublehdr_filtered(:,:,c)imfilter(hdr_img(:,:,c),gaussian_kernel,replicate);endend2.3 双边滤波函数手写无工具箱依赖functionhdr_filteredbilateral_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_winfloor(window_size/2);% 预计算空间权重所有像素共享提高效率[x,y]meshgrid(-half_win:half_win,-half_win:half_win);spatial_weightexp(-(x.^2y.^2)/(2*sigma_s^2));hdr_filteredzeros(H,W,C,like,hdr_img);% 对每个通道独立处理forc1:Cfprintf( 双边滤波处理第%d/%d通道...\n,c,C);channelhdr_img(:,:,c);% 遍历每个像素嵌套循环较慢大图像建议用GPU/积分图优化见下文扩展建议fori1:Hforj1:W% 提取当前像素的邻域i_minmax(1,i-half_win);i_maxmin(H,ihalf_win);j_minmax(1,j-half_win);j_maxmin(W,jhalf_win);% 邻域内的像素值neighbor_valschannel(i_min:i_max,j_min:j_max);% 当前像素值center_valchannel(i,j);% 截取对应的空间权重sp_winspatial_weight(...(i_min:i_max)-ihalf_win1,...(j_min:j_max)-jhalf_win1);% 计算值域权重exp(-(邻域值-中心值)²/(2*sigma_r²))range_weightexp(-(neighbor_vals-center_val).^2/(2*sigma_r^2));% 总权重 空间权重 × 值域权重归一化total_weightsp_win.*range_weight;total_weighttotal_weight/sum(total_weight(:));% 加权求和得到滤波后的值hdr_filtered(i,j,c)sum(neighbor_vals(:).*total_weight(:));endendendend2.4 色调映射函数仅用于显示functiontm_imgreinhard_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%分位数避免过曝L0.2126*hdr_img(:,:,1)0.7152*hdr_img(:,:,2)0.0722*hdr_img(:,:,3);% 亮度通道white_pointprctile(L(:),99);% Reinhard色调映射公式L_tmL.*(1L./(white_point^2))./(1L);% 亮度归一化L_tmL_tm/max(L_tm(:));% 保持颜色比例对每个通道应用相同的亮度缩放scaleL_tm./(Leps);tm_imghdr_img.*scale;% Gamma校正tm_imgtm_img.^(1/gamma);% 截断到0-1区间tm_imgmax(min(tm_img,1),0);end2.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 直接运行将所有函数保存为.m文件放在同一文件夹运行主脚本hdr_filtering_main.m程序会自动加载MATLAB自带的memorial.hdr示例也可手动选择自己的HDR文件3.2 参数调优建议参数建议范围效果gaussian_sigma_s1~3值越大全局模糊越强去噪效果越好但边缘越糊bilateral_sigma_s1~3空间邻域范围值越大平滑范围越广bilateral_sigma_r全局标准差的5%~20%值越小对像素值差异越敏感边缘保留越好但去噪效果弱值越大越接近高斯滤波filter_window5~11奇数窗口越大计算量呈平方增长建议不超过113.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图像上会非常慢以下是优化方案GPU加速将数据转为gpuArray滤波后在转回CPUhdr_img_gpugpuArray(hdr_img);% 滤波函数在GPU上运行需修改代码适配GPU数组hdr_filtered_gpubilateral_filter_hdr(hdr_img_gpu,sigma_s,sigma_r,window_size);hdr_filteredgather(hdr_filtered_gpu);积分图优化预计算积分图和平方积分图将双边滤波的时间复杂度从O(N×k²)降到O(N)适合大图像。使用内置函数若安装了Computer Vision Toolbox可直接用imbilatfilt速度比手写快10倍以上% 内置双边滤波支持HDRhdr_bilateralimbilatfilt(hdr_img,SpatialSigma,2,RangeSigma,params.bilateral_sigma_r);降采样滤波先将HDR降采样到1/2或1/4尺寸滤波后再上采样速度提升4~16倍精度损失很小。五、应用场景扩展HDR去噪双边滤波是HDR图像去噪的首选比高斯滤波保留更多细节适合夜景、高对比度场景。HDR tone mapping预处理滤波后可以减少色调映射后的伪影让亮部过渡更自然。HDR纹理平滑比如游戏HDR贴图的纹理平滑保留硬边缘的同时平滑表面纹理。高动态场景目标检测滤波后降低噪声提升检测算法对亮部和暗部目标的识别率。