Local dimming algorithm in matlab
LED局部背光算法的matlab仿真
最近公司接了华星光电（TCL）的一个项目LCD-BackLight-Local-Diming-Algorithm-IP ，由于没有实际的硬件，只能根据客户给的论文 算法进行调研，评估和确认。即先理解论文的算法，再用MATLAB或OpenCV仿真，再通过视觉或客观图像评价指标评估算法效果，最后通过对几种论文算法的实验仿真效果分析比较确定一种算法，作为fpga实现。最大感触是得搞清楚实际情况再行动！
一 论文算法原理
论文题目为<< Backlight Local Dimming Algorithm for High Contrast LCD-TV>>，我一般看英语论文先用excel表格记录论文编号，题目，背景，关键词，主要方法，缺点（可以创新的地方），结论等核心点，接着看文献顺序为：摘要，结论，实验，着重看论文算法（准确理解，很耗时），若很闲则会看Introduction。当然，这次对这个算法的改进是多次，每次蹦出一个想法就改善算法，主要有三点：一是多用参数化定义，提高算法移植性和通用性。二是多想办法减少重复劳动，如在一个窗口中一次显示多个参数取不同值的结果，这一点很重要，以后多用。三是，一个算法的改进周期很长，直接方法是看高手的写法，而且得多次对比和看文档书籍才有提升，而不是成谜自己的小成果。
下面直接给出算法整体框架：
从上图可以看出整个算法分为四步，LED亮度强度计算，暗区域增强，空间滤波器，时间滤波器。注意点1：算法的仿真是放在FPGA上通过对LED的亮度调节，即整个算法是对LED块的分块的亮度进行处理，由于没有LED等实际硬件，本次仿真是对图像分块后的每个图像块亮度强度进行处理。注意点2：算法的实验是进行视频仿真，而MATLAB主要是进行图像处理，故本次时间滤波器的处理，这一步包含当前帧和上一帧，通过用potplayer软件截取一帧视频图像，方法如下：
第一步是视频素材下载网站：
视频下载网站（水印）
https://www.vjshi.com/watch/6162152.html
视频素材下载（无水印）
https://mixkit.co/free-stock-video/nature/
第二步是PotPlayer的截取视频的一帧图像：在用potpalyer播放视频时，快捷键crtl+G打开连续截取画面，按照下面的12345步骤即可成功截取想要的两帧图像。
论文算法的参数如下表格所示：
上面的参数在第一步到第三步都会有所涉及，这里提前给出。
（一）LED亮度强度计算
1.算法介绍
Linit(n)为第n个LED块的初始亮度强度，第一步的输出；LMAX可允许的最大LED亮度强度水平；TL 是一个预定义参数，用来控制局部调节水平，由表1可知为0.0625；PNUM（第n个图像块的总像素数量）；
Hn(i)是第n个图像块的直方图；W(i)是预定义的权重向量，论文里为变量i的平方。
2.MATLAB代码实现
从上面公式1可知，算法整体为两项中取最小值，第二项有点繁琐，那就从第二项开始：第二项的TL，Pnum都很好计算，一个为常数0.0625，一个为图像块的总像素数，即图像的尺寸M*N，图像行列的乘积；
第二项中的图像块直方图和权重W(i)的乘积的累加和需注意，图像块的直方图怎么实现，若是调用imhist仅仅是显示直方图，而在这里明显不对，故需要查资料弄懂图像的直方图的含义，通过看<<数字图像处理>>这本书，明白图像的直方图：对应灰度级的像素数。而MATLAB怎么表示，通过在MATLAB中搜索help imhist和不断看相关函数发现了histogram函数。在MATLABhelp文档中，知道了一种用法可以表示图像某个灰度级的像素数：h = histogram(In,256);counts=h.Values ;通过这两行代码把第n个图像块的直方图搞定了，那么公式1则轻松表示出来了。故对公式1的处理关键：拆分，弄明白直方图的含义，通过经典的数字图像处理书籍和matlab官方的文档，准确理解直方图含义。由于后面的步骤要调用这个结果，故把这一步封装成函数，调用很方便。代码如下：
1 function L_init = Linit(In) 2 [m,n]=size(In); 3 num_p=m*n; 4 s=0; 5 h = histogram(In,256); 6 counts=h.Values ; 7 for i2 = 1:256 8 h = counts(i2); 9 ii= i2^2; 10 s=s+(h*ii); %求得所有像素与灰度级平方的乘积。 11 % s=s+(h*i2); 12 end 13 num_p = double (num_p); 14 s = double (s); 15 L_temp=(0.0625/num_p)*s; 16 L_init = min(63,L_temp); 17 end
LED亮度强度计算
（二）暗区域增强
通过第一步的matlab代码实现，我们可以学会MATLAB编写方法：第一对公式的每个变量准备理解，如最大值，像素总数，灰度级，权重变量。第二对公式拆分，把握关键的一点，即难一点实现的，如上面的直方图。三是查阅权威书籍<<数字图像处理>>和matlab官方文档，弄懂直方图的实现，再把算法完整编写。对于这一步的暗区域增强，实现过程相似。
1.算法介绍
Lben 是增强的LED亮度强度灰度；Lmean 是当前图像的所有图像块的 Linit的平均值；TM用来定义暗区域的阈值；TB用来控制暗区域增强水平的预定义参数；
2.MATLAB代码实现
通过上面的分析可知，Linit的使用通过调用第一步的函数输出值即可，而Lmean 的计算则需要注意，它是包含图像的所有块的Linit的均值，故需要用一个矩阵存放所有图像块的Linit的值，再用matlab的mean2函数即可，而后面的TM，TB用通过表格可知分别是5和3。整体是选择结构，用matlab的if语句实现即可。
1 function L_ben = Lben(Lini,L_mean) 2 TM = 5; 3 if(Lini<L_mean | TM<L_mean) % L_mean = 20 4 L_ben = Lini; 5 else 6 L_ben =min(63 ,Lini+3*(Lini-L_mean)); 7 end 8 end
暗区域增强Lben
（三）空间滤波器处理
1.算法介绍
LSF(n)是第n个LED块的亮度强度的空间滤波器的输出；TSF 是用来平滑滤波器的预定义参数；Φ(n)是大小为3X3的第n个块的邻域中心；这一步的关键是对3*3的领域空间的理解，以第n块为中心，共3*3总共为9个图像块的领域。这个领域用一个3*3的矩阵索引实现，即Lben(m)是领域中每一块的 Lben 值。
2.MATLAB代码实现
其实关键点就是3*3的邻域实现，我的想法是，将10*10的Lben最外围全部舍掉，即行和列范围都是2-9，然后舍掉本身即第n块的Lben，这总共8个元素，构成一个数组，取数组的最大值和第n块的Lben比较取较大值，得到空间滤波器的输出，每一块都能得到一个输出。思路的代码如下：
1 function L_SF = S_filter(LSF_ini,temp_Lben,T_SF) 2 %function L_SF = S_filter(LSF_ini,T_SF) 3 temp_Lsf = temp_Lben; 4 [row,col]=size(temp_Lben); 5 LSF_temp1=zeros(1,9); 6 for nx=2:row-1 7 for ny=2:col-1 8 % 取出对应需要进行运算区域的数据 9 LSF_temp1 = temp_Lsf(nx-1:nx+1,ny-1:ny+1); 10 % convert to 11 temp11=[LSF_temp1(1,:),LSF_temp1(2,1),LSF_temp1(2,3),LSF_temp1(3,:),]; 12 % 赋值 13 temp_max = max(temp11); 14 end 15 end 16 Len_nn = LSF_ini; 17 L_SF = max(Len_nn,temp_max-T_SF); 18 end
空间滤波器
(四) 时间滤波
1.算法介绍
LTF(n)是时间滤波器的输出；(k) 和 (k-1)指当前帧和前一帧；R (0<R ≤1)为控制IIR低通滤波器的平滑程度的参数。
其中，R由视频内容自适应决定，其表达式如下：
TTF 用来控制IIR滤波器形状的预定义参数，为0.125。Pmean整个图像的归一化平均像素值。
2.MATLAB代码实现
时间滤波由上面两个公式(4)(5)组成，包含两个部分，即第一步(5)是公式(4)中的R，故我们可以先编写matlab函数实现R，再带入公式(4)；第二步LSF为空间滤波的输出，可以封装为一个函数；在顶层调用这两部分的函数即可实现公式(4)的功能。
R的代码，第一步是计算两帧图像的归一化平均像素值，两帧图像这里为上面用potplayer软件截取的视频图像，归一化平均像素值那么需要用到归一化函数 normalize，还有平均值函数mean和最小值函数min,按照公式(4)直接敲即可：
%calculate the result of R in EQ.(4) for algorithm %the input is two frame of image %the output is R in EQ.(4) function R_S4 = cal_R(IMA,IMA1) %RGB = imread(\'../pic/EE01000_1000.jpeg\');%1000*1000 pixel of image %RGB1 = imread(\'../pic/EE4_1000.jpeg\');%1000*1000 pixel of EE4_1000.jpeg %THE RESULT OF above ima is R = 0.4626 %RGB = imread(\'ee0_1000.jpg\');%1000*1000 pixel of image I =rgb2gray(IMA); I1 =rgb2gray(IMA1); %normalize the pixel of ima I = double(I); I1 = double(I1); I_nr=normalize(I,\'range\'); I1_nr = normalize(I1,\'range\'); Pmean0 = mean(I_nr,\'all\');%the normalize mean of previous frame Pmean1 = mean(I1_nr,\'all\');%the normalize mean of current frame R_S4_temp = 0.125+Pmean0+Pmean1; R_S4 = min(1,R_S4_temp);%the parameter R of step4 for algorithm end
公式(4)还调用了第三步的公式(3)空间滤波的输出，故我们可以把第三步空间滤波的输出封装成一个函数如下所示，在时间滤波这一步直接调用即可：
1 %main function begin 2 function Lsf_out = cal_Lsf(image) 3 tic; 4 RGB = image; 5 %RGB = imread(\'ee0_1000.jpg\');%1000*1000 pixel of image 6 %parameter 7 global ima_width; 8 global ima_hight; 9 global block_width; 10 global block_hight; 11 global block_num_h; 12 global block_num_v; 13 %parameter 14 %I =rgb2gray(RGB); 15 %normalize the pixel of ima 16 %---Table 1. Thresholds and paramters---- 17 global Y_M; 18 Y_M = 230;%the most bright scene 19 global T_L; 20 % T_L = 0.025; 21 %T_L is a predefined parameter 22 %-used to control the local dimming degree 23 global T_M; 24 T_M = 18;%the most dark scene 25 global T_B; 26 T_B = 3; 27 global T_SF; 28 T_SF = 120; 29 30 %----Table 1. Thresholds and paramters--- 31 IMA_YCBCR = rgb2ycbcr(RGB); 32 Y=IMA_YCBCR(:,:,1); 33 % figure 34 % imshow(I); 35 temp_Linit = zeros(block_num_v,block_num_h); 36 % % I1 = rgb2gray(I); 37 t_row = 0:block_hight:ima_hight; % the row\'coordinates of each block 38 t_row1= t_row; 39 t_row = t_row+1;%pattention-->reduced add the Extra 1 40 t_col = 0:block_width:ima_width; % the column\'coordinates of each block 41 t_col1= t_col; 42 t_col = t_col+1; 43 % % save space to predifine variable 44 % %temp = repmat(int8(0), 100, 100);%arry cannot deliver to point 45 temp1 = cell(block_num_v,block_num_h);% creat cell struct 46 %the number of block in row or column 47 for i = 1 : block_num_v 48 for j = 1 : block_num_h 49 temp = Y(t_row(i):t_row1(i+1), t_col(j):t_col1(j+1)); 50 temp1{i,j}=temp; 51 In = temp1{i,j}; 52 temp_Linit(i,j)= Linit(In); 53 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 54 end 55 end 56 % %%cal L_ben by cycle 57 % 58 L_mean = mean2(temp_Linit); 59 temp_Lben = zeros(block_num_v,block_num_h); 60 for i1 = 1 : block_num_v 61 for j1 = 1 : block_num_h 62 Lben_In = temp_Linit(i1,j1); 63 temp_Lben(i1,j1)= Lben(Lben_In,L_mean); 64 end 65 end 66 % %%combine L_SF to arry 67 %T_SF = 20; 68 temp_Lsf = zeros(block_num_v,block_num_h); 69 for i4 = 1 : block_num_v 70 for j4 = 1 : block_num_h 71 Lsf_In = temp_Lben(i4,j4); 72 temp_Lsf(i4,j4)= S_filter(Lsf_In,temp_Lben); 73 end 74 end 75 % 76 % %YCBCR1_temp{i,j}=temp_Lsf(i,j); 77 Lsf_out = temp_Lsf; 78 end 79 %main function end 80 81 % %CONTINU TO l_init, L_ben,L_sf 82 function L_init = Linit(In) 83 global Y_M; 84 global T_L; 85 86 [m,n]=size(In); 87 num_p=m*n; 88 s=0; 89 % h = histogram(In,256); 90 %counts=h.Values ; 91 counts = imhist(In,256); 92 for i2 = 1:256 93 h = counts(i2); 94 ii= (i2-1)^2; 95 s=s+(h*ii); %求得所有像素与灰度级平方的乘积。 96 % s=s+(h*i2); 97 end 98 num_p = double (num_p); 99 s = double (s); 100 L_temp=(T_L/num_p)*s; 101 L_init = min(Y_M,L_temp); 102 end 103 % % %part 2 the enhance method 104 function L_ben = Lben(Lini,L_mean) 105 global T_M; 106 global Y_M; 107 global T_B; 108 if(Lini<L_mean | T_M<L_mean) % L_mean = 20 109 L_ben = Lini; 110 else 111 L_ben =min(Y_M ,Lini+T_B*(Lini-L_mean)); 112 end 113 end 114 % % %part 3 the Spatial filter 115 % rng(0,\'twister\');%初始化随机数生成器 116 % a = 1; 117 % b = 64; 118 % rr = (b-a).*rand(16,1) + a;%creat a rand matrix 119 % LSF=reshape(rr,4,4); 120 % %contains 16 datas the value range from 1 to 64 121 % for i = 1:4 122 % for j = 1:4 123 % LSF_i= LSF(i,j); 124 % ll = S_filter(LSF_i,i,j,20); 125 % end 126 % end 127 function L_SF = S_filter(LSF_ini,temp_Lben) 128 %function L_SF = S_filter(LSF_ini,T_SF) 129 global T_SF; 130 temp_Lsf = temp_Lben; 131 [row,col]=size(temp_Lben); 132 LSF_temp1=zeros(1,9); 133 for nx=2:row-1 134 for ny=2:col-1 135 % 取出对应需要进行运算区域的数据 136 LSF_temp1 = temp_Lsf(nx-1:nx+1,ny-1:ny+1); 137 % convert to 138 temp11=[LSF_temp1(1,:),LSF_temp1(2,1),LSF_temp1(2,3),LSF_temp1(3,:),]; 139 % 赋值 140 temp_max = max(temp11); 141 end 142 end 143 Len_nn = LSF_ini; 144 L_SF = max(Len_nn,temp_max-T_SF); 145 end 146 147 148 %设置循环，对每个像素进行八邻域中值运算 149 %设置循环条件
LSF Function
如上面代码所示，含有全局参数，这是增加代码的移植性和通用性，即简单改参数就能仿真不同输入和条件下的图像效果。
（一）时间滤波的顶层函数
由于这个LED的局部调光算法的核心思想：
根据图像的内容分块，以及对不同块使用不同的系数（图像块或LED块的亮度强度）去调节背光（达到高对比度，即亮的地方更亮，暗的地方更暗）。
上面是我对LED局部调光算法的总结。故要想实现这个算法，得搞定怎么对图像分块，用MATLAB表示第n块图像，对第n块图像实现上面的（一）（二）（三）步处理，这是难点。我刚开始直接想的是MATLAB现成的函数：B = blockproc(A,[m n],fun)函数。但这个函数的输入为原图像，输出直接和原图同大小的一幅图像，不能实现我们想要的第n块图像。我在这卡了很久，通过大量查阅官网的图像分块资料和看数字图像处理书籍，及网上查阅csdn博客和MATLAB中文论坛，终于找到办法：即数字图像就是一个多维矩阵，像素点对应矩阵元素，故图像分块就是矩阵的分块。如本文的图像分块为10*10共100块，输入图像分辨率大小为1000*1000的rgb图像，则每一图像块为100*100的矩阵。我们把这个100*100的矩阵赋值给一个10*10的元胞数组的一个元素（关键，由于矩阵不能赋值给一点），然后把元胞数组的元素赋值给临时变量temp，这个大小为100*100的矩阵temp就是第n块的图像块，再传递给上面的函数调用即可。当然，针对10*10共100块的图像块处理，会频繁用到for循环结构，幸亏上面把每一步封装成函数，我们调用即可，不然很难搞。做matlab算法和FPGA的设计一样，都是循序渐进，即模块化设计（一个功能一个函数），把每个模块仿真实现无误后，再顶层调用即可。为了提高算法的通用性，可以把图像的宽和高等参数，表示为一个参数，以后输入不同分辨率图像直接修改参数即可，本次就省略哦。LED局部调光算法的主函数如下：
1 %calculate the luminous of rgb ima by the whole algorithm 2 %display the modified luminous of rgb image 3 %input image 4 %output the modified luminous of RGB image 5 %date 08.04 6 close all 7 clear 8 clc; 9 %load(\'S1.mat\'); 10 tic; 11 % %blockproc 12 % 13 % % myfun = @(block_struct) block_struct.data; 14 % %I = imread(\'e1.jpeg\');%100*100 pixel of image 15 % RGB = imread(\'../pic/meadowP1000_1000.jpg\');%1000*1000 pixel of image previous frame 16 % RGB1 = imread(\'../pic/meadowC1000_1000.jpg\');%1000*1000 pixel 17 % RGB = imread(\'../pic/fireworksP1000_1000.jpg\');%1000*1000 pixel of image previous frame 18 % RGB1 = imread(\'../pic/fireworksC1000_1000.jpg\');%1000*1000 pixel 19 RGB = imread(\'../pic/moonP1000_1000.jpg\');%1000*1000 pixel of image previous frame 20 RGB1 = imread(\'../pic/moonC1000_1000.jpg\');%1000*1000 pixel 21 %RGB = imread(\'ee0_1000.jpg\');%1000*1000 pixel of image 22 %calculate the R of between current frame and previous frame 23 R = cal_R(RGB,RGB1); 24 %calculate the Ltf of the Nth image blcok 25 temp_Lsf0 = cal_Lsf(RGB);%Lsf of previous frame 26 temp_Lsf1 = cal_Lsf(RGB1); 27 Ltf0 = R*temp_Lsf0; 28 Ltf1 = (R*temp_Lsf1)+((1-R)*Ltf0);%Ltf of current frame 29 temp_Lsf = Ltf1; 30 % %YCBCR1_temp{i,j}=temp_Lsf(i,j); 31 tt_row = 0:100:1000; % the row\'coordinates of each block 32 tt_row1= tt_row; 33 tt_row = tt_row+1;%pattention-->reduced add the Extra 1 34 tt_col = 0:100:1000; % the column\'coordinates of each block 35 tt_col1= tt_col; 36 tt_col = tt_col+1; 37 % % save space to predifine variable 38 % %temp = repmat(int8(0), 100, 100);%arry cannot deliver to point 39 YCBCR1_out = zeros(1000); 40 YCBCR1_temp = cell(10);% creat cell struct 41 len1 = 10; %the number of block in row or column 42 for i5 = 1 : len1 43 for j5 = 1 : len1 44 YCBCR1_temp{i5,j5}=temp_Lsf(i5,j5); 45 YCBCR1_out(tt_row(i5):tt_row1(i5+1), tt_col(j5):tt_col1(j5+1))=YCBCR1_temp{i5,j5}; 46 47 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 48 end 49 end 50 % figure 51 % imshow(RGB);title(\'RGB\'); 52 YCBCR = rgb2ycbcr(RGB); 53 % figure 54 % imshow(YCBCR);title(\'YCBCR\'); 55 YCBCR(:,:,1)=YCBCR1_out; 56 OUT = ycbcr2rgb(YCBCR); 57 OUT1 = OUT + RGB;% comibine the original image and modified image 58 %OUT1 the output of 59 % figure 60 % imshow(OUT);title(\'OUT\'); 61 % figure 62 % imshow(OUT1);title(\'OUT1\'); 63 toc; 64 %离y轴距离，离x轴距离，子图宽，子图高 65 subplot(\'Position\',[0.1 0.1 0.8 0.8]);imshow(RGB);title(\'RGB\'); 66 % [left bottom width height]. 67 figure 68 subplot(1,3,1);imshow(RGB);title(\'RGB\'); 69 subplot(1,3,2); imshow(YCBCR);title(\'YCBCR\'); 70 subplot(1,3,3);imshow(OUT1);title(\'OUT1\'); 71 % % for m = 1:10 72 % % for n = 1:10 73 % % temp2=(temp1{m,n}); 74 % % subplot(10, 10, 10*(m-1)+n); imshow(temp2); 75 % % end 76 % % end 77 % %CONTINU TO l_init, L_ben,L_sf 78 79 %设置循环，对每个像素进行八邻域中值运算 80 %设置循环条件
main Code
当然，这里的代码就是公式(4)的顶层代码。
(二)高效的仿真手段
回头看看最开始的（一 论文算法原理）Table 1表格，竟然五个参数！天啊，做实验恐怖得一周，要是手动一个个改。这种体力活一直都是这样搞的，这次想偷偷懒，没必要搞重复劳动。自然想到了，要是直接一次显示多个参数取不同值的结果，多爽哈！以前一周的实验，现在几分钟万事！哈哈，这就是效率。
怎么一次显示多个参数取不同值？
思路是：直接把算法封装成一个函数，用for循环遍历不同值，结合subplot和title函数，在一个窗口显示输入不同参数的对应图像。
1.简单例子1–直接上图：一次显示一个参数取四个不同值的图像效果
上面这一个窗口，包含一个参数a取四个不同值的图像效果，这个简单的例子为图像增强的线性变换，这个线性变换的增强函数如下：
for a = 1:4
out = a*RGB+0.5;
subplot(2,2,a);imshow(out);
% title(sprintf(a,b));
%title(num2str(4*(a-1)+b));%one title display different images from 1 to 16
title(num2str(a));%Create Multiline Title
end
2 例子2为一次显示两个参数各取四个不同值的图像效果
如上图所示，一次显示参数a和b,一个为对比度增强系数，一个为亮度，两个参数在1到4之间取值的16中组合结果，代码如下：
for a = 1:4
for b = 1:4
out = a*RGB2+b;
subplot(4,4,(4*(a-1)+b));imshow(out);
% title(sprintf(a,b));
%title(num2str(4*(a-1)+b));%one title display different images from 1 to 16
title({num2str(a);num2str(b)});%Create Multiline Title
end
end
可以看出关键为subplot的坐标和title中变量的写法。直接搜索相关函数文档理解。既然可以显示两个参数取不同值，自然能够推广到多个参数。
3 输入任意参数的用法（一次显示参数取多个不同值的结果）
在上面的两个例子中，我们为了着重设计思路，将参数取为1：4这种情况，但实际你会碰到参数取值30-50或者0.01到0.05这种不是从1开始或每次增加1.既然，一个参数，两个参数都会了，那么任意参数只需要转换为上面的两个例子即可，如下图的代码所示：
上面的25张图像，为水平和垂直方向的分区数量的参数，各自取30：50到范围的仿真效果，MATLAB代码如下：
for block_num_v = 30:5:50
for block_num_h = 30:5:50
a =((block_num_v-30)/5)+1;
b =((block_num_h-30)/5)+1;
out = A1_Localdim_814(RGB,RGB1);
subplot(5, 5, 5*(a-1)+b);imshow(out);
% title(sprintf(a,b));
%title(num2str(4*(a-1)+b));%one title display different images from 1 to 16
title({num2str(block_num_v);num2str(block_num_h)});%Create Multiline Title
end
end
可以看出，对于输入多个任意参数，用到数值转换：即把输入参数减去起点30，除以公差（数列），加上1就转换为我们的线性转换例子1和2。
注意点：在调用的子函数中，不能包含有画图相关的函数，如figure, subplot, imshow等，否则就会把参数不同取值的仿真结果分别显示在多个窗不同口中，如下图所示：
那是在调用的函数GABF_DDE_top814中有下面几行代码：
figure;
subplot(1,1,1);
imshow(test_R);
相应的例子代码如下：
1 %this file is used to 2 %once display the result of Algorithm for 20 images 3 %pattention:the function can\'t have a figure 4 %--subplot(1,1,1);or imshow(test_R2); 5 clc; 6 close all; 7 clear; 8 %----a example for display 4 images---- 9 I = cell(2); 10 I{1,1}=load(\'../data/S1.mat\',\'S\'); 11 I{1,2}=load(\'../data/S2.mat\',\'S\'); 12 I{2,1}=load(\'../data/S3.mat\',\'S\'); 13 I{2,2}=load(\'../data/S4.mat\',\'S\'); 14 15 for i =1:2 16 for j = 1:2 17 temp_I = I{i,j}; 18 temp = temp_I.S; 19 out = GABF_DDE_top814(temp); 20 subplot(2,2,(2*(i-1)+j)); 21 imshow(out);title({num2str(i);num2str(j)}); 22 end 23 end 24 %----a example for display 4 images---- 25 %extend to display 20 images
GABF_DDE_top814
通过对实现一次显示参数取不同值的结果的分析，利用这个方法，我们就能极大的节省时间，偷个懒！
4 matlab的图像算法的结构
如上图所示，有三层结构，最顶层为仿真测试，其次为算法顶层，最后是功能函数。仿真测试顶层，用来做不同输入参数和不同输入图像的实验，如上面讲的参数取不同值的方法，输入不同分辨率图像，灰度图像，RGB图像等。算法顶层，即能够完整实现算法功能，核心就是能够封装成一个函数，实现算法的效果。如本次算法为LCD局部背光调节算法，就能够实现根据图像的内容，对一帧图像进行分块M*N块后，进行算法流程的四步处理后，达到预期的效果：即让亮的区域更亮，暗的区域更暗。从而提升视频图像对比度和亮度，输出高质量的视觉图像，给人更好的视觉体验。功能函数，对应于具体的算法步骤，如算法流程中第一步求初始亮度，空间滤波等，多个函数可以放在一个MATLAB文件里，也可以单独放一个文件。为了理解，冗余复杂的函数可以单独放一个文件，简单的函数多个放一起。
(一)测试顶层函数
1 %This file is used to 2 %simulate A1_Localdim_814.m or verify it 3 4 %------------1000*1000 of ima------------- 5 % RGB = imread(\'../pic/meadowP1000_1000.jpg\');%1000*1000 pixel of image previous frame 6 % RGB1 = imread(\'../pic/meadowC1000_1000.jpg\');%1000*1000 pixel 7 % RGB = imread(\'../pic/fireworksP1000_1000.jpg\');%1000*1000 pixel of image previous frame 8 % RGB1 = imread(\'../pic/fireworksC1000_1000.jpg\');%1000*1000 pixel 9 RGB = imread(\'../pic/moonP1000_1000.jpg\');%1000*1000 pixel of image previous frame 10 RGB1 = imread(\'../pic/moonC1000_1000.jpg\');%1000*1000 pixel 11 %------------1000*1000 of ima------------- 12 %------------1920*1080 of ima------------- 13 % RGB = imread(\'../pic/sunset_pre1920_1080.jpg\');%1920*1080 pixel of image previous frame 14 % RGB1 = imread(\'../pic/sunset_current1920_1080.jpg\');%1920*1080 pixel 15 % RGB = imread(\'../pic/blinking_round_lights_pre1920_1080.jpg\'); 16 % RGB1 = imread(\'../pic/blinking_round_lights_cur1920_1080.jpg\'); 17 % RGB = imread(\'../pic/volumetric_pre1920_1080.jpg\');%volumetric_pre1920_1080 18 % RGB1 = imread(\'../pic/volumetric_cur1920_1080.jpg\'); 19 RGB2 = imread(\'../pic/LYF2.jpg\'); 20 %------------1920*1080 of ima------------- 21 I = rgb2gray(RGB); 22 YCBCR = rgb2ycbcr(RGB); 23 YCBCR_Y =( YCBCR(:,:,1)); 24 % subplot(1,1,1);imshow(YCBCR_Y);title(\'YCBCR_Y\'); 25 %%test 1 26 %-------the result of different numbers of blocks for a ima------------- 27 %parameters 28 tic; 29 global block_num_h; 30 global block_num_v; 31 global T_L; 32 T_L = 0.025; 33 % %---1000*1000 of ima---- 34 % block_num_v=50; 35 % block_num_h=50; %the number of block in the horizontal direction 36 % %---1000*1000 of ima---- 37 %---1920*1080 of ima----the parameters is a multiple of 2 38 %block_num_h=64,32,16 or 8 39 %block_num_v=40,20,10 40 % block_num_v=40; 41 % block_num_h=64; %the number of block in the horizontal direction 42 %---1000*1000 of ima---- 43 %parameters 44 for block_num_v = 30:5:50 45 for block_num_h = 30:5:50 46 a =((block_num_v-30)/5)+1; 47 b =((block_num_h-30)/5)+1; 48 out = A1_Localdim_814(RGB,RGB1); 49 subplot(5, 5, 5*(a-1)+b);imshow(out); 50 % title(sprintf(a,b)); 51 %title(num2str(4*(a-1)+b));%one title display different images from 1 to 16 52 title({num2str(block_num_v);num2str(block_num_h)});%Create Multiline Title 53 end 54 end 55 %-------the result of different numbers of blocks for a ima----------- 56 %%test 2 57 %-------the result of different values of global T_L;----------- 58 % global T_L; 59 % for T_L = 0.01:0.005:0.05 60 % a =((T_L-0.01)/0.005)+1; 61 % out = A1_Localdim_814(RGB,RGB1); 62 % subplot(3,3,a);imshow(out); 63 % % title(sprintf(a,b)); 64 % %title(num2str(4*(a-1)+b));%one title display different images from 1 to 16 65 % title(num2str(T_L));%Create Multiline Title 66 % end 67 %example--once display different values of parameter at the same figure 68 % for a = 1:4 69 % for b = 1:4 70 % out = a*RGB2+b; 71 % subplot(4,4,(4*(a-1)+b));imshow(out); 72 % % title(sprintf(a,b)); 73 % %title(num2str(4*(a-1)+b));%one title display different images from 1 to 16 74 % title({num2str(a);num2str(b)});%Create Multiline Title 75 % end 76 % end 77 %example 78 toc;
test_demo814
主要包含三个测试：一是在水平和垂直方向上不同分区数量的分块时的仿真效果。二是Table 1表格中T_B取不同值的仿真效果。三是上面分析的提高算法效率的例子，一次显示多个参数取不同值的仿真效果。
(二)算法顶层
1 %calculate the luminous of rgb ima by the whole algorithm 2 %display the modified luminous of rgb image 3 %input image 4 %output the modified luminous of RGB image 5 %date 08.12 6 %main work is adjust the parameters of Table 1 7 % close all 8 % clear 9 % clc; 10 %load(\'S1.mat\'); 11 function YCBCR1_out = A1_Localdim_814(RGB,RGB1) 12 13 % %blockproc 14 % 15 % % myfun = @(block_struct) block_struct.data; 16 %------------1000*1000 of ima------------- 17 % RGB = imread(\'../pic/meadowP1000_1000.jpg\');%1000*1000 pixel of image previous frame 18 % RGB1 = imread(\'../pic/meadowC1000_1000.jpg\');%1000*1000 pixel 19 % RGB = imread(\'../pic/fireworksP1000_1000.jpg\');%1000*1000 pixel of image previous frame 20 % RGB1 = imread(\'../pic/fireworksC1000_1000.jpg\');%1000*1000 pixel 21 % RGB = imread(\'../pic/moonP1000_1000.jpg\');%1000*1000 pixel of image previous frame 22 % RGB1 = imread(\'../pic/moonC1000_1000.jpg\');%1000*1000 pixel 23 %------------1000*1000 of ima------------- 24 %------------1920*1080 of ima------------- 25 % RGB = imread(\'../pic/sunset_pre1920_1080.jpg\');%1920*1080 pixel of image previous frame 26 % RGB1 = imread(\'../pic/sunset_current1920_1080.jpg\');%1920*1080 pixel 27 % RGB = imread(\'../pic/blinking_round_lights_pre1920_1080.jpg\'); 28 % RGB1 = imread(\'../pic/blinking_round_lights_cur1920_1080.jpg\'); 29 % RGB = imread(\'../pic/volumetric_pre1920_1080.jpg\');%volumetric_pre1920_1080 30 % RGB1 = imread(\'../pic/volumetric_cur1920_1080.jpg\'); 31 %------------1920*1080 of ima------------- 32 I =rgb2gray(RGB); 33 I1 =rgb2gray(RGB1); 34 %% parameters... 35 %The horizontal direction---h 36 %The vertical direction---v 37 global ima_width; %the number of column for image 38 global ima_hight; %the number of row 39 global block_width; %the number of column for a block 40 global block_hight; %the number of row for a block 41 global block_num_h; 42 global block_num_v; 43 44 [ima_hight,ima_width]=size(I); 45 % block_num_v=40; 46 % block_num_h=64; %the number of block in the horizontal direction 47 48 block_width = floor(ima_width/block_num_h); 49 block_hight = floor(ima_hight/block_num_v); 50 %% 51 %calculate the R of between current frame and previous frame 52 global T_TF; 53 T_TF = 0.125; 54 R = cal_R(RGB,RGB1); 55 %calculate the Ltf of the Nth image blcok 56 temp_Lsf0 = cal_Lsf(RGB);%Lsf of previous frame 57 temp_Lsf1 = cal_Lsf(RGB1); 58 Ltf0 = R*temp_Lsf0; 59 Ltf1 = (R*temp_Lsf1)+((1-R)*Ltf0);%Ltf of current frame 60 temp_Lsf = Ltf1; 61 % %YCBCR1_temp{i,j}=temp_Lsf(i,j); 62 tt_row = 0:block_hight:ima_hight; % the row\'coordinates of each block 63 tt_row1= tt_row; 64 tt_row = tt_row+1;%pattention-->reduced add the Extra 1 65 tt_col = 0:block_width:ima_width; % the column\'coordinates of each block 66 tt_col1= tt_col; 67 tt_col = tt_col+1; 68 % % save space to predifine variable 69 % %temp = repmat(int8(0), 100, 100);%arry cannot deliver to point 70 YCBCR1_out = zeros(ima_hight,ima_width); 71 YCBCR1_temp = cell(block_num_v,block_num_h);% creat cell struct 72 %the number of block in row or column 73 for i5 = 1 : block_num_v 74 for j5 = 1 : block_num_h 75 YCBCR1_temp{i5,j5}=temp_Lsf(i5,j5); 76 YCBCR1_out(tt_row(i5):tt_row1(i5+1), tt_col(j5):tt_col1(j5+1))=YCBCR1_temp{i5,j5}; 77 78 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 79 end 80 end 81 % figure 82 % imshow(RGB);title(\'RGB\'); 83 84 % figure 85 % subplot(1,2,1); imshow(YCBCR1_out);title(\'YCBCR1_out\'); 86 YCBCR = rgb2ycbcr(RGB); 87 YCBCR_Y =( YCBCR(:,:,1)); 88 %------min/max luminous-------- 89 %find the max luminous of ima 90 %the min luminous 91 Y_MAX = max(YCBCR_Y,[],\'all\'); 92 Y_MIN = min(YCBCR_Y,[],\'all\'); 93 %------min/max luminous-------- 94 YCBCR1_out = uint8(YCBCR1_out); 95 % figure 96 % subplot(1,2,2); imshow(YCBCR1_out);title(\'YCBCR1_out\'); 97 % subplot(1,2,1);imshow(YCBCR_Y);title(\'YCBCR_Y\'); 98 % figure 99 % imshow(YCBCR);title(\'YCBCR\'); 100 % YCBCR(:,:,1)=YCBCR1_out; 101 102 % OUT = ycbcr2rgb(YCBCR); 103 % OUT1 = OUT + RGB;% comibine the original image and modified image 104 %OUT1 the output of 105 % figure 106 % imshow(OUT);title(\'OUT\'); 107 % figure 108 % imshow(OUT1);title(\'OUT1\'); 109 110 end 111 %离y轴距离，离x轴距离，子图宽，子图高 112 % subplot(\'Position\',[0.1 0.1 0.8 0.8]);imshow(RGB);title(\'RGB\'); 113 % [left bottom width height]. 114 % figure 115 % subplot(1,3,1);imshow(RGB);title(\'RGB\'); 116 % subplot(1,3,2); imshow(YCBCR);title(\'YCBCR\'); 117 % subplot(1,3,3);imshow(OUT1);title(\'OUT1\');
A1_Localdim_814
上面为局部背光调节算法，能够实现预期的算法效果。核心是用到全局变量，这样算法具备通用性。可以根据输入图像分辨率和分块数量，分块的大小调节对应的全局参数，而不用每次改代码，太麻烦了。
(三)功能函数
1.计算空间滤波的输出
1 %main function begin 2 function Lsf_out = cal_Lsf(image) 3 tic; 4 RGB = image; 5 %RGB = imread(\'ee0_1000.jpg\');%1000*1000 pixel of image 6 %parameter 7 global ima_width; 8 global ima_hight; 9 global block_width; 10 global block_hight; 11 global block_num_h; 12 global block_num_v; 13 %parameter 14 %I =rgb2gray(RGB); 15 %normalize the pixel of ima 16 %---Table 1. Thresholds and paramters---- 17 global Y_M; 18 Y_M = 230;%the most bright scene 19 global T_L; 20 % T_L = 0.025; 21 %T_L is a predefined parameter 22 %-used to control the local dimming degree 23 global T_M; 24 T_M = 18;%the most dark scene 25 global T_B; 26 T_B = 3; 27 global T_SF; 28 T_SF = 120; 29 30 %----Table 1. Thresholds and paramters--- 31 IMA_YCBCR = rgb2ycbcr(RGB); 32 Y=IMA_YCBCR(:,:,1); 33 % figure 34 % imshow(I); 35 temp_Linit = zeros(block_num_v,block_num_h); 36 % % I1 = rgb2gray(I); 37 t_row = 0:block_hight:ima_hight; % the row\'coordinates of each block 38 t_row1= t_row; 39 t_row = t_row+1;%pattention-->reduced add the Extra 1 40 t_col = 0:block_width:ima_width; % the column\'coordinates of each block 41 t_col1= t_col; 42 t_col = t_col+1; 43 % % save space to predifine variable 44 % %temp = repmat(int8(0), 100, 100);%arry cannot deliver to point 45 temp1 = cell(block_num_v,block_num_h);% creat cell struct 46 %the number of block in row or column 47 for i = 1 : block_num_v 48 for j = 1 : block_num_h 49 temp = Y(t_row(i):t_row1(i+1), t_col(j):t_col1(j+1)); 50 temp1{i,j}=temp; 51 In = temp1{i,j}; 52 temp_Linit(i,j)= Linit(In); 53 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 54 end 55 end 56 % %%cal L_ben by cycle 57 % 58 L_mean = mean2(temp_Linit); 59 temp_Lben = zeros(block_num_v,block_num_h); 60 for i1 = 1 : block_num_v 61 for j1 = 1 : block_num_h 62 Lben_In = temp_Linit(i1,j1); 63 temp_Lben(i1,j1)= Lben(Lben_In,L_mean); 64 end 65 end 66 % %%combine L_SF to arry 67 %T_SF = 20; 68 temp_Lsf = zeros(block_num_v,block_num_h); 69 for i4 = 1 : block_num_v 70 for j4 = 1 : block_num_h 71 Lsf_In = temp_Lben(i4,j4); 72 temp_Lsf(i4,j4)= S_filter(Lsf_In,temp_Lben); 73 end 74 end 75 % 76 % %YCBCR1_temp{i,j}=temp_Lsf(i,j); 77 Lsf_out = temp_Lsf; 78 end 79 %main function end 80 81 % %CONTINU TO l_init, L_ben,L_sf 82 function L_init = Linit(In) 83 global Y_M; 84 global T_L; 85 86 [m,n]=size(In); 87 num_p=m*n; 88 s=0; 89 % h = histogram(In,256); 90 %counts=h.Values ; 91 counts = imhist(In,256); 92 for i2 = 1:256 93 h = counts(i2); 94 ii= (i2-1)^2; 95 s=s+(h*ii); %求得所有像素与灰度级平方的乘积。 96 % s=s+(h*i2); 97 end 98 num_p = double (num_p); 99 s = double (s); 100 L_temp=(T_L/num_p)*s; 101 L_init = min(Y_M,L_temp); 102 end 103 % % %part 2 the enhance method 104 function L_ben = Lben(Lini,L_mean) 105 global T_M; 106 global Y_M; 107 global T_B; 108 if(Lini<L_mean | T_M<L_mean) % L_mean = 20 109 L_ben = Lini; 110 else 111 L_ben =min(Y_M ,Lini+T_B*(Lini-L_mean)); 112 end 113 end 114 % % %part 3 the Spatial filter 115 % rng(0,\'twister\');%初始化随机数生成器 116 % a = 1; 117 % b = 64; 118 % rr = (b-a).*rand(16,1) + a;%creat a rand matrix 119 % LSF=reshape(rr,4,4); 120 % %contains 16 datas the value range from 1 to 64 121 % for i = 1:4 122 % for j = 1:4 123 % LSF_i= LSF(i,j); 124 % ll = S_filter(LSF_i,i,j,20); 125 % end 126 % end 127 function L_SF = S_filter(LSF_ini,temp_Lben) 128 %function L_SF = S_filter(LSF_ini,T_SF) 129 global T_SF; 130 temp_Lsf = temp_Lben; 131 [row,col]=size(temp_Lben); 132 LSF_temp1=zeros(1,9); 133 for nx=2:row-1 134 for ny=2:col-1 135 % 取出对应需要进行运算区域的数据 136 LSF_temp1 = temp_Lsf(nx-1:nx+1,ny-1:ny+1); 137 % convert to 138 temp11=[LSF_temp1(1,:),LSF_temp1(2,1),LSF_temp1(2,3),LSF_temp1(3,:),]; 139 % 赋值 140 temp_max = max(temp11); 141 end 142 end 143 Len_nn = LSF_ini; 144 L_SF = max(Len_nn,temp_max-T_SF); 145 end 146 147 148 %设置循环，对每个像素进行八邻域中值运算 149 %设置循环条件
cal_Lsf
2.计算第四步时间滤波的参数R
1 %calculate the result of R in EQ.(4) for algorithm 2 %the input is two frame of image 3 %the output is R in EQ.(4) 4 function R_S4 = cal_R(IMA,IMA1) 5 %RGB = imread(\'../pic/EE01000_1000.jpeg\');%1000*1000 pixel of image 6 %RGB1 = imread(\'../pic/EE4_1000.jpeg\');%1000*1000 pixel of EE4_1000.jpeg 7 %THE RESULT OF above ima is R = 0.4626 8 %RGB = imread(\'ee0_1000.jpg\');%1000*1000 pixel of image 9 % I =rgb2gray(IMA); 10 % I1 =rgb2gray(IMA1); 11 %normalize the pixel of ima 12 global T_TF; 13 % I = double(I); 14 % I1 = double(I1); 15 IMA = double(IMA); 16 IMA1 = double(IMA1); 17 I_nr=normalize(IMA,\'range\'); 18 I1_nr = normalize(IMA1,\'range\'); 19 Pmean0 = mean(I_nr,\'all\');%the normalize mean of previous frame 20 Pmean1 = mean(I1_nr,\'all\');%the normalize mean of current frame 21 R_S4_temp = T_TF+Pmean0+Pmean1; 22 R_S4 = min(1,R_S4_temp);%the parameter R of step4 for algorithm 23 end
cal_R
(四)实验图像
上面为T_B取不同值的实验
实验现象简单分析，算法处理后，输出的图像变暗很多，即亮度改变很大，但视觉效果差，和原图差很远，即细节增强效果不好。
本文算法主要是对图像的亮度调整，故关键点2是颜色空间转换，即要实现对亮度的调整，得把RGB图像转换为YCbCr，然后Y分量即RGB图像的亮度，把算法的结果1000*1000的矩阵的点赋值给Y分量，其他两个分量不变，之后将YCbCr形式的图像转换回RGB图像显示，这时的RGB图像的亮度已经是我们算法的处理结果。关键点1是图像的分块（怎么表示第n块图像块），分块的逆问题：把分块后的图像块组合成一幅图像，通用是利用到主函数里的思路，数字图像即矩阵，即将每个图像块赋值给矩阵的某个范围的行和列即可，如100*100的图像块直接赋值给1000*1000矩阵行和列的1：100范围，把十个图像块循环赋值完，即得到一个矩阵。
第二组实验图像：
在上面的图f1中，为经过时间滤波处理后的图像，算法的输出OUT1与原图RGB相比，明显更加明亮鲜艳，整副图像的亮度提升显著，如黑色方框圈出的阳光与原图相比更加耀眼，红色圆圈标注的为暗区域，但算法处理后亮度也有提升，故算法有人工添加的痕迹，整体对比度和亮度改善明显。图g1中，为时间滤波的结果，与f1相比，几乎一致，肉眼很难区别。后面的实验与f1的分析类似，整体亮度和对比度显著增强。
实验结果讨论：
算法分析，一是从实验看分区100个与1600个，扩大16倍并无明显的视觉改善。二是有无时间滤波，对图像的影响很小。三是，用亮度分量代替RGB进行算法处理，几乎一致。
若是用到时间滤波，即 公式4，那么首先得计算两帧图像的统计量，即缓存三帧四帧，可以先RGB转换为灰度图像，若是硬件缓存2560*1600@90HZ的RGB视频，若仅仅提取亮度分量Y，一是视频传输带宽为2560*1600*90=368.64MB/S，这仅仅是最小的理论带宽，一般能够满足为千兆网，若是RGB24视频传输带宽，则为1.1GB/S。二是帧缓存，至少两帧用于计算，那么得考虑缓存三帧四帧，缓存大小=2560*1600*4=15.625兆字节。
三 论文算法复现总结
（一）清楚认识世界比你咋想重要得多（首要核心）
首先得准确理解算法，在算法的设计实现过程中犯了五个误区。
第一是如算法的处理对象，即LED块的亮度强度，得用颜色空间转换YCbCr才行，而我刚开始直接当成图像的灰度，当成和大多数图像处理算法一样的思路，惯性思维导致。
第二是，算法是对视频处理，我看到最后才发现这一点，当前帧和上一帧，真的很突然，发现自己想错了，以为搞不来视频仿真，后面才懂直接截取视频两帧图像即可。
第三是，RGB图像的R分量一开始搞错了，直接写成R = RGB(1, :, :);% R = 1*M*3. 正确写法如下：R = RGB(:, :, 1);% R = M*N.
第四是对于直方图的写法，刚开始直接调用了网上的，后面仔细分析发现统计的灰度值，和直方图代表一个灰度级的像素数差很远。
第五是，图像的像素总数搞成了图像的灰度总数，这是没有查阅权威资料和文档。这是犯了没有仔细思考，以科学事实为原则。
（二）循序渐进（拆分原则）
即设计一个算法，把他拆分成很多步骤或模块组成，例子1如在银行取钱，第一步插入银行卡，第二步输入密码，第三步输入取款金额，第三步取出银行卡离开。例子2，时钟的表针分成三部分，时针，分针，秒针。故本次算法也用拆分，由四个步骤组成，那么就一步步的实现，把第一步实现了做下一步；或者各个击破，把每一步分别实现，组合起来。正确的设计方法，是模块化设计，循序渐进，把每一个部分当成一个函数，直接在主函数调用即可，当遇到问题时，由于是一个模块，且模块间关联不大，很容易解决。很错误就是想一下子把整个功能实现，由于代码量很大，一旦出现问题，会很复杂的逻辑关系，到时候是问题解决你了，败给了时间。
（三）用多种方式去做
多查阅官方文档，多看权威书籍，多勇敢尝试不同方法，以实践为核心。
方法1：最开始想复制粘贴别人的算法，发现搜了半天，如CSDN博客，知乎上怎么查找论文的代码列出的网站，博客园，GitHub上，谷歌搜索，搜作者的主页等方式，都没找到，只能另想他法。
方法2：先问了一个同事，他的想法是可以用C++做图像处理，为我开启了一种新思路，但我一直用matlab做红外图像细节增强算法，好歹搞出一篇论文，故不想放弃MATLAB图像处理的基础。他还启发了我，对图像亮度处理，用颜色空间转换: RGB-YCbCr-RGB,让我发现了自己的坑，我设计的算法对图像灰度处理，错的离谱。
方法3：自己编代码。很耗时，同样很多问题伴随，这就是工科生，在问题中收获知识和成长！
了解世界和自己！