(续)LED局部背光算法MATLAB仿真
在上一篇博客<Local dimming algorithm in matlab>中，我们实现了对一篇论文的算法用matlab仿真。在本篇论文中，对另一篇论文进行了MATLAB仿真。
这篇论文<<A Novel Two-Dimensional Adaptive Dimming Technique of X-Y Channel Drivers for LED Backlight System in LCD TVs >>和前一篇一样，都是新的ocal dimming算法，前一篇博客的算法我们记为算法1，如下图1.2所示， 这篇论文的算法框架相对复杂很多，我们简记为算法2，如下图Fig .2所示：
当然，实现算法的前提是准确理解算法的原理，这一步是很难的，而编写代码仅仅是照图施工。一般我理解论文的算法，都是先按照上一篇博客顺序通读全文，然后重点看算法原理部分。对于算法原理部分：法一反复三遍，不光看论文的数学公式介绍，还会结合图表对照；法二对于论文的算法中的某个表达式不理解，我还会网上查阅和看对应的书籍如数字图像处理；法三若是前面两种方式还是困惑，我会看论文的参考文献和引用它的文献；通过这三种方式后，绝大多数SCI论文的算法都能准确理解。
一 论文的算法原理
论文<<A Novel Two-Dimensional Adaptive Dimming Technique of X-Y Channel Drivers for LED Backlight System in LCD TVs >>的算法流程如Fig.2所示，和图1.2的算法思路一样，即先将图像分块，通过一系列处理改变图像块的亮度，把修改亮度后的图像块组合成一幅图像输出。在硬件上实现则是把算法计算得到的亮度作为LED分块后的背光亮度。这篇论文算法的总体思路：图像像素的亮度由背光亮度(BL)和图像的透射率(T)两者相乘得到。故算法围绕这两个指标展开。
（一）背光亮度BL的计算
由上面Fig .2知，背光亮度BL由四个步骤得到：第一步求图像块的MLD；第二步求图像块的调光因子k；第三步，求图像块归一化的背光亮度；第四步求总的背光亮度。下面我们将对每一步i详细介绍：
1.图像块的MLD
(1)原理
理解这个需要费点心思，直接从论文的英语句子(maximum level data)意思是不准确的,结合下面图(b)可知，MLD表示图像块中某一行或列的最大灰度值，针对一个图像块的所有行的最大MLD和列的分别用MLDrow,m和MLD:col,m 表示，其数学表达式分别如公式(1)和(2)
根据上面对MLD的分析，MLD表示一个图像块中某行或列的最大灰度等级，故在公式(1)中：MLD:m,1为第一行的最大灰度值，MLD:m,2为则表示第m行的最大灰度值，以此类推MLD:m,N为第N行的，则公式(1)的左边MLD:row,m为图像块中所有行的最大灰度等级（其灰度等级范围为0-255）；公式（2）和公式（1）类似，表示列的最大灰度等级，不再赘述。
(2)MATLAB代码
对于算法的设计，其原理是核心。准确理解了数学表达式(1)和(2)的原理，那么MATLAB代码仅仅是工具，很简单，很难的是准确理解这个论文中的数学表达式，刚开始我仅看论文的描述错误地理解为编码值，后面结合图9(b)才理解原来是图像块的最大灰度值。对于MATLAB代码编写，数字图像就是一个二维矩阵，故某行的最大灰度值就是矩阵某行的最大值，故所有行的最大灰度值MLD:row,m为矩阵中所有行的最大值，直接命令窗help max函数，仔细查看这个max函数如何表达即可。其代码如下：
1 AA = matrix_block; 2 MLD_row_m =double( max(max(AA,[],2)));%Eq.(1) 3 MLD_col_m =double( max(max(AA,[],1)) ); %Eq.(2)
2.图像块的调节因子km,n
经过上面的第一步得到所有行或列的最大灰度值，在此基础上，通过三个数学公式(下图的456)得到调节因子，(4)和(5)为求透射率，在原论文的公式(1)中有介绍。在这里很容易犯错，即把公式(4)的y理解为gamma变换，但结合论文的Fig. 8知为欧拉常数0.5772，由于经过这个y变换后，其值小于1，而gamma变换则大于1.当然理解这个是欧拉常数，除了看原论文的介绍和Fig, 8,还需要百度搜索gamma变换的信息和看数学书籍
(2)MATLAB代码
只需知道y为欧拉常数，公式456很简单，调节因子只需一个函数几行代码：
1 function k_m_n = diming_k(matrix_block) 2 y_ol = 0.5772; %Euler's constant 3 AA = matrix_block; 4 MLD_row_m =double( max(max(AA,[],2)));%Eq.(2) 5 MLD_col_m =double( max(max(AA,[],1)) ); %Eq.(3) 6 D_row_m = (MLD_row_m/255)^y_ol;%Eq.(4) 7 D_col_m = (MLD_col_m/255)^y_ol;%Eq.(5) 8 %Eq.(6) 9 k_m_n = min(D_row_m,D_col_m); 10 end
3.求图像块归一化的背光亮度
(1)原理
理解这一步，着重点是归一化，在matlab函数中，归一化表示把某个变量的范围转换到[0 1]区间，而灰度图的范围为[0 255]，如下图公式(11)所示，f(x-xn,y-ym)表示第(m,n)块图像在亮度归一化后峰值坐标点的亮度值（在0-1区间）,结合下面图11(b),每个图像块的归一化亮度为峰值亮度，是一个常数。理解了这个f函数的含义，而调节因子km,n (已知)，则轻松得到图像块归一化的背光亮度BLm,n(x,y)
(2)MATLAB代码
在matlab函数中，有一个归一化函数normalize，其中的一个用法Nr = normalize(A,'range');就是把变量的范围转换为[0 1]区间，而这里是将图像块的亮度归一化，normalize(A,'range');输入为图像块的亮度矩阵则就得到了归一化的亮度，峰值亮度直接调用max函数得到，之后与第二步调节因子函数的输出相乘即可。
1 function peak_luminance_out = luminance_trans(Y_luminance_IN) 2 A = Y_luminance_IN; 3 A = double (A); 4 Nr = normalize(A,'range'); 5 peak_luminance_out = max(Nr,[],'all'); % the max value of matrix A 6 %end 7 end
4.总的背光亮度
(1)原理
在公式(11)中求出了一个块的归一化背光亮度，公式(12)为总的背光亮度，它等于所有块的背光亮度的和，由下面公式(12)的表达可知，总的背光亮度为M*N块的累加和
(2)MATLAB代码
由于我们已经得到了一个图像块的背光亮度，接着把每个图像块的背光亮度累加即可，在matlab中累加求和，直接用两层for循环得到：
1 mm = 10; 2 nn = 10; 3 BLmn_out = 0; 4 for i3 = 1 : mm 5 for j3 = 1 : nn 6 temp_k_mn3 = (temp_kmn1(i3,j3)); 7 temp_f_luminance = f_luminance_temp(i3,j3); 8 BLmn_out_mn = temp_k_mn3* temp_f_luminance; 9 BLmn_out = BLmn_out + BLmn_out_mn; 10 end 11 end
我们刚开始在<<一 论文算法原理>>中总结了图Fig .2的算法流程，像素的亮度= BL*T.其中BL是背光亮度，我们已经得到，接下是求T，透射率。
（二）最终透射率Tfinal的计算
1.the LC transmittance
(1)原理
透射率T和编码值(0-255)的关系符合gamma曲线，如图8所示，从这个图中，我们能避免犯错，由于T的值为[0 1]，故y为欧拉常数0.5772，他们的数学表达式如公式(1)：
(2)MTLAB代码
由于公式(1)表达式为中括号里的幂函数，而幂为0.5772，输入的编码值为变量，取值为[0 255],matlab代码如下：
1 function T = LC_trans(cv) 2 % for cv = 1:255 3 y_ol = 0.5772; %Euler's constant 4 T = (cv/255)^y_ol; 5 %end 6 end
2. The enhanced Tm,n
(1)原理
这一步的透射率与上面公式1不同，它是在第一步基础上除以调节因子km,n
(2)MATLAB代码
由于在上一步中，我们已经得到T，而且调节因子km,n在求背光亮度的第二步中做了详细介绍，故此处仅仅将两者调用即可，函数如下：
1 function TT_m_n = TT_trans(Tmn_IN,k1_m_n) 2 % for cv = 1:255 3 %y_ol = 0.5772; %Euler's constant 4 TT_m_n = (Tmn_IN/k1_m_n); 5 %end 6 end
3. The corrected code value
1.原理
在这一步中，将公式(1)的编码值进行了变形，修改为如公式(10)所示的形式：
2.MATLAB代码
直接按照公式(10)的表达式编写matlab代码即可
1 function BLD_m_n = BLD_trans(Kmn_IN) 2 % for cv = 1:255 3 y_ol = 0.5772; %Euler's constant 4 BLD_m_n = (255)*(Kmn_IN^(1/y_ol)); 5 %end 6 end
而要得到最终的透射率Tfina 将上面公式(10)的输出带入公式(1)即可。通过这一系列运算，我们得到了背光亮度和最终的透射率，之后相乘得到像素亮度为每一个图像块的亮度，直接替换原图像的亮度，再显示修改亮度后的整幅图像即可，我们将完整代码粘贴如下;
1 %matlab code 2 %%MLD:row,m 3 %matrix A express the block 4 %max(A) is the col vector of matrix A 5 close all 6 clear 7 clc; 8 %cv(m,n)begin 9 % max_gv = max(block_mn,[],'all'); 10 % min_gv = min(block,mn,[],'all'); 11 % cv_mn = max_gv - min_gv; 12 %end 13 %X = [1 2;3 4];% 14 % max(X,[],1) 15 %return the row vector of the max value of each column for X 16 % max(X,[],2) 17 %return the column vector of the max value of each row for X 18 RGB = imread('ee0_1000.jpg');%1000*1000 pixel of image%1000*1000 pixel of image 19 I =rgb2gray(RGB); 20 %I = mat2gray (I);%convert the range of intensity to [0 1] for image 21 figure 22 imshow(I); 23 % divide image into block 10*10 begin 24 temp_kmn1 = zeros(10,10); 25 t_row = 0:100:1000; % the row'coordinates of each block 26 t_row1= t_row; 27 t_row = t_row+1;%pattention-->reduced add the Extra 1 28 t_col = 0:100:1000; % the column'coordinates of each block 29 t_col1= t_col; 30 t_col = t_col+1; 31 temp1 = cell(10);% creat cell struct 32 len = 10; %the number of block in row or column 33 for i = 1 : len 34 for j = 1 : len 35 temp = I(t_row(i):t_row1(i+1), t_col(j):t_col1(j+1)); 36 temp1{i,j}=temp; 37 In = temp1{i,j}; 38 temp_kmn1(i,j) = diming_k(In); 39 40 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 41 end 42 end 43 %power reduction rate of Eq.(7) 44 power_kmn_in = temp_kmn1; 45 power_rate_out = power_rate(power_kmn_in); 46 power_rate_out 47 %block for luminance Y BEGIN 48 YCBCR = rgb2ycbcr(RGB); 49 Y = YCBCR(:,:,1); 50 tt_row = 0:100:1000; % the row'coordinates of each block 51 tt_row1= tt_row; 52 tt_row = tt_row+1;%pattention-->reduced add the Extra 1 53 tt_col = 0:100:1000; % the column'coordinates of each block 54 tt_col1= tt_col; 55 tt_col = tt_col+1; 56 Y1_temp = cell(10); %the cell struct of 10*10 57 f_luminance_temp = zeros(10);%the zero matrix of 10*10 58 len1 = 10; %the number of block in row or column 59 for i2 = 1 : len1 60 for j2 = 1 : len1 61 62 Y1_temp{i2,j2} = Y(tt_row(i2):tt_row1(i2+1), tt_col(j2):tt_col1(j2+1)); 63 f_luminance_temp(i2,j2) = luminance_trans( Y1_temp{i2,j2}); 64 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 65 end 66 end 67 %block for luminance Y END 68 69 % Eq.(12) begin 70 mm = 10; 71 nn = 10; 72 BLmn_out = 0; 73 for i3 = 1 : mm 74 for j3 = 1 : nn 75 temp_k_mn3 = (temp_kmn1(i3,j3)); 76 temp_f_luminance = f_luminance_temp(i3,j3); 77 BLmn_out_mn = temp_k_mn3* temp_f_luminance; 78 BLmn_out = BLmn_out + BLmn_out_mn; 79 end 80 end 81 %Eq.(12)end 82 83 BLmn_out 84 %Substitute formula 10 into Formula 1 begin 85 y_ol = 0.5772; %Euler's constant 86 len2 = 10; 87 CV_mn_E10p = zeros(10); 88 T_temp1 = zeros(10); 89 temp_BL_E8 = zeros(10); 90 Block_luminance_fina = zeros(10); 91 for i4 = 1 : len2 92 for j4 = 1 : len2 93 temp1_E910 = temp1{i4,j4}; 94 temp1_cv_mn = cv89_trans(temp1_E910);%cv_mn in Eq.(9) and Eq.(10) 95 temp1_Kmn_E10 = diming_k(temp1_E910); 96 Kmn_E10 = (temp1_Kmn_E10)^(1/y_ol); 97 CV_mn_E10p(i4,j4) = temp1_cv_mn/Kmn_E10; %the output of Eq.(10) 98 T_temp1(i4,j4) = LC_trans( CV_mn_E10p(i4,j4)); 99 temp_BL_E8(i4,j4) = BLD_trans(Kmn_E10); 100 T_fina = T_temp1(i4,j4); 101 BL_fina = temp_BL_E8(i4,j4); 102 Block_luminance_fina(i4,j4) = T_fina*BL_fina; 103 end 104 end 105 %the final luminance of block in-- 106 %--B. Dimming Algorithm for the Proposed LED Backlight end 107 %end 108 %luminance of the whole image begin 109 ftt_row = 0:100:1000; 110 ftt_row1= ftt_row; 111 ftt_row = ftt_row+1;%pattention-->reduced add the Extra 1 112 ftt_col = 0:100:1000; % the column'coordinates of each block 113 ftt_col1= ftt_col; 114 ftt_col = ftt_col+1; 115 len3 = 10; %the number of block in row or column 116 for i5 = 1 : len3 117 for j5 = 1 : len3 118 BLT= Block_luminance_fina(i5,j5); 119 Y(ftt_row(i5):ftt_row1(i5+1), ftt_col(j5):ftt_col1(j5+1))=BLT; 120 Y_OUT = Y; 121 %subplot(10, 10, 10*(i-1)+j); imshow(temp); 122 end 123 end 124 %end 125 YCBCR1 = rgb2ycbcr(RGB); 126 YCBCR1(:,:,1) = Y_OUT; 127 OUT = ycbcr2rgb(YCBCR1); 128 figure 129 imshow(OUT); 130 OUT1 = OUT + RGB;% comibine the original image and modified image 131 %OUT1 the output of 132 figure 133 imshow(OUT1);title('OUT1'); 134 %Eq.(13) 135 function cc_xy = ccxy_trans(cv_x_y_IN) 136 % for cv = 1:255 137 y_ol = 0.5772; %Euler's constant 138 cc_xy =cv_x_y_IN /((BLmn_out)^(1/y_ol)); 139 %end 140 end 141 % Eq.(13)end 142 %Eq.(11) f(x-xn,y-ym) 143 %normalized backlight luminance of each division block 144 function peak_luminance_out = luminance_trans(Y_luminance_IN) 145 A = Y_luminance_IN; 146 A = double (A); 147 Nr = normalize(A,'range'); 148 peak_luminance_out = max(Nr,[],'all'); % the max value of matrix A 149 %end 150 end 151 %Eq.(10) 152 function cc_v = ccv_trans(TT_m_n_IN) 153 % for cv = 1:255 154 y_ol = 0.5772; %Euler's constant 155 cc_v = 255*(TT_m_n_IN)^(1/y_ol); 156 %end 157 end 158 159 %Eq.(9) 160 function TT_m_n = TT_trans(Tmn_IN,k1_m_n) 161 % for cv = 1:255 162 %y_ol = 0.5772; %Euler's constant 163 TT_m_n = (Tmn_IN/k1_m_n); 164 %end 165 end 166 167 %cv(m,n) the gray level of block(m,n) begin 168 % in EQ.(10) and Eq.(9) 169 function cv89_m_n = cv89_trans(cv89_IN) 170 block_mn = cv89_IN; 171 max_gv = max(block_mn,[],'all'); 172 min_gv = min(block_mn,[],'all'); 173 cv89_m_n = max_gv - min_gv; 174 end 175 %cv(m,n) end 176 %Eq.(8) 177 function BLD_m_n = BLD_trans(Kmn_IN) 178 % for cv = 1:255 179 y_ol = 0.5772; %Euler's constant 180 BLD_m_n = (255)*(Kmn_IN^(1/y_ol)); 181 %end 182 end 183 %Eq.(7) 184 function power_out = power_rate(temp_kmn) 185 temp_out = 0; 186 [m,n] = size(temp_kmn); 187 for i2 = 1 : m 188 for j2 = 1 : n 189 temp_k_m_n = (temp_kmn(i2,j2)/(m*n)); 190 temp_out = temp_out+ temp_k_m_n; 191 power_out = 1 - temp_out; 192 end 193 end 194 end 195 196 197 %block end 198 199 %Eq.(1) 200 function T = LC_trans(cv) 201 % for cv = 1:255 202 y_ol = 0.5772; %Euler's constant 203 T = (cv/255)^y_ol; 204 %end 205 end 206 %the factor k_m_n from Eq.(2) to Eq.(6) 207 function k_m_n = diming_k(matrix_block) 208 y_ol = 0.5772; %Euler's constant 209 AA = matrix_block; 210 MLD_row_m =double( max(max(AA,[],2)));%Eq.(2) 211 MLD_col_m =double( max(max(AA,[],1)) ); %Eq.(3) 212 D_row_m = (MLD_row_m/255)^y_ol;%Eq.(4) 213 D_col_m = (MLD_col_m/255)^y_ol;%Eq.(5) 214 %Eq.(6) 215 k_m_n = min(D_row_m,D_col_m); 216 end
二 仿真及实验
在这次仿真实验，我们做了一组实验，并且包含算法1，如下所示：
（a1）原图 （b1）算法1的输出结果 （c1）算法2的输出结果
（d1）原图和算法1叠加后的输出 （e1）原图和算法2叠加后的输出
图1.2.1：不同算法的实验结果
在图1.2.1中，（b1）为算法1的输出结果，它作为背光亮度的图像与原图（a1）很接近，亮度略有增加，对比度有着显著提升；（c1）为算法2的输出结果，但由于亮度变化不均匀，出现了很多不同亮度级的小块，与原图相比亮度过度增强；（d1）为原图和算法1叠加后的输出，与原图（a1）相比，细节得到适度增强，比如图像中的眼睛更清晰，亮度的适度提升导致图像的对比度提高明显，视觉效果更佳；（e1）为原图和算法2叠加后的输出，与原图（a1）相比，由于亮度的过度增强，图像（e1）出现了光晕，鼻子和嘴不能呈现出来，且由于亮度变化不连续，出现了分块。
三 总结
对新的一篇论文进行了复现，并和上一篇的实验效果作为比较，本文算法的实验效果不如算法1的效果，故应用在工业上，可能优先考虑算法1。上一篇论文复现用了整整一周，这次的论文复现仅仅三天，有做过的经验和基础是很重要，万事开头难。第一天精读论文三遍，理解算法原理。第二天用matlab代码实现了一半的算法，第三天结合实现算法1的思路，直接套用在算法2上。
这次的感悟还是准确理解算法原理是最重要的，代码就是描述而已。