(續)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上。
這次的感悟還是准確理解算法原理是最重要的，代碼就是描述而已。