Local dimming algorithm in matlab

• 2021 年 7 月 3 日
• 筆記
• 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主要是进行图像处理，故本次省略了最后一步时间滤波器的处理，这一步包含当前帧和上一帧。

论文算法的参数如下表格所示：

 上面的参数在第一步到第三步都会有所涉及，这里提前给出。

（一）LED亮度强度计算

1.算法介绍

 Linit(n)为第n个LED块的初始亮度强度，第一步的输出；LMAX可允许的最大LED亮度强度水平；TL 是一个预定义参数，用来控制局部调节水平，由表1可知为0.0625；P_NUM（第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亮度强度灰度；L_mean 是当前图像的所有图像块的 L_init的平均值；T_M用来定义暗区域的阈值；T_B用来控制暗区域增强水平的预定义参数；

2.MATLAB代码实现

     通过上面的分析可知，L_init的使用通过调用第一步的函数输出值即可，而L_mean 的计算则需要注意，它是包含图像的所有块的L_init的均值，故需要用一个矩阵存放所有图像块的L_init的值，再用matlab的mean2函数即可，而后面的T_M，T_B用通过表格可知分别是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

空间滤波器

二 主函数的设计

   由于这个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 close all
 2 clear
 3 clc;
 4 %load('S1.mat');
 5 
 6 % %blockproc  
 7 % 
 8 % % myfun = @(block_struct) block_struct.data;
 9 % %I = imread('e1.jpeg');%100*100 pixel of image
10 RGB = imread('ee2_1000.jpg');%1000*1000 pixel of image
11 I =rgb2gray(RGB);
12 % figure 
13 % imshow(I);
14  temp_Linit = zeros(10,10);
15 % % I1 = rgb2gray(I);
16 t_row = 0:100:1000;   % the row'coordinates of each block
17 t_row1= t_row;
18 t_row = t_row+1;%pattention-->reduced add the Extra 1
19 t_col = 0:100:1000;   % the column'coordinates of each block
20 t_col1= t_col;
21 t_col = t_col+1;
22 % % save space to predifine variable
23 % %temp = repmat(int8(0), 100, 100);%arry cannot deliver to point
24 temp1 = cell(10);% creat cell struct
25 len = 10;          %the number of block in row or column
26 for i = 1 : len
27     for j = 1 : len
28         temp = I(t_row(i):t_row1(i+1), t_col(j):t_col1(j+1));
29         temp1{i,j}=temp;
30         In =  temp1{i,j};
31         temp_Linit(i,j)= Linit(In);
32  %subplot(10, 10, 10*(i-1)+j); imshow(temp);
33     end
34 end
35 % %%cal L_ben by cycle
36 % 
37 L_mean = mean2(temp_Linit);
38  temp_Lben = zeros(10,10);
39 for i1 = 1 : 10
40     for j1 = 1 : 10
41         Lben_In =   temp_Linit(i1,j1);
42         temp_Lben(i1,j1)= Lben(Lben_In,L_mean);
43     end
44 end
45 % %%combine L_SF to arry
46 %T_SF = 20;
47  temp_Lsf = zeros(10,10);
48 for i4 = 1 : 10
49     for j4 = 1 : 10
50         Lsf_In =   temp_Lben(i4,j4);
51         temp_Lsf(i4,j4)= S_filter(Lsf_In,temp_Lben,20);
52     end
53 end
54 % 
55 % %YCBCR1_temp{i,j}=temp_Lsf(i,j);
56 tt_row = 0:100:1000;   % the row'coordinates of each block
57 tt_row1= tt_row;
58 tt_row = tt_row+1;%pattention-->reduced add the Extra 1
59 tt_col = 0:100:1000;   % the column'coordinates of each block
60 tt_col1= tt_col;
61 tt_col = tt_col+1;
62 % % save space to predifine variable
63 % %temp = repmat(int8(0), 100, 100);%arry cannot deliver to point
64 YCBCR1_out = zeros(1000);
65 YCBCR1_temp = cell(10);% creat cell struct
66 len1 = 10;          %the number of block in row or column
67 for i5 = 1 : len1
68     for j5 = 1 : len1
69         YCBCR1_temp{i5,j5}=temp_Lsf(i5,j5);
70    YCBCR1_out(tt_row(i5):tt_row1(i5+1), tt_col(j5):tt_col1(j5+1))=YCBCR1_temp{i5,j5};
71 
72  %subplot(10, 10, 10*(i-1)+j); imshow(temp);
73     end
74 end
75 
76  YCBCR = rgb2ycbcr(RGB);
77   figure
78  imshow(YCBCR);
79 YCBCR(:,:,1)=YCBCR1_out;
80 RGB1 = ycbcr2rgb(YCBCR);
81 figure
82 imshow(RGB);
83 figure
84 imshow(RGB1);

main code

实验效果如下：

实验现象简单分析，算法处理后，输出的图像变暗很多，即亮度改变很大，但视觉效果差，和原图差很远，即细节增强效果不好。

本文算法主要是对图像的亮度调整，故关键点2是颜色空间转换，即要实现对亮度的调整，得把RGB图像转换为YCbCr，然后Y分量即RGB图像的亮度，把算法的结果1000*1000的矩阵的点赋值给Y分量，其他两个分量不变，之后将YCbCr形式的图像转换回RGB图像显示，这时的RGB图像的亮度已经是我们算法的处理结果。关键点1是图像的分块（怎么表示第n块图像块），分块的逆问题：把分块后的图像块组合成一幅图像，通用是利用到主函数里的思路，数字图像即矩阵，即将每个图像块赋值给矩阵的某个范围的行和列即可，如100*100的图像块直接赋值给1000*1000矩阵行和列的1：100范围，把十个图像块循环赋值完，即得到一个矩阵。

三 论文算法复现总结

（一）清楚认识世界比你咋想重要得多（首要核心）

       首先得准确理解算法，在算法的设计实现过程中犯了五个误区。

       第一是如算法的处理对象，即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：自己编代码。很耗时，同样很多问题伴随，这就是工科生，在问题中收获知识和成长！

了解世界和自己！