(續)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
BL(m,n)
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
BL
我們剛開始在<<一 論文演算法原理>>中總結了圖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
the transmittance of LC
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
ALL Code
二 模擬及實驗
在這次模擬實驗，我們做了一組實驗，並且包含演算法1，如下所示：
（a1）原圖 （b1）演算法1的輸出結果 （c1）演算法2的輸出結果
（d1）原圖和演算法1疊加後的輸出 （e1）原圖和演算法2疊加後的輸出
圖1.2.1：不同演算法的實驗結果
從圖1.2.1到圖1.2.3為兩種不同演算法的三組實驗，在圖1.2.1中，（b1）為演算法1的輸出結果，它作為背光亮度的影象與原圖（a1）很接近，亮度略有增加，對比度有著顯著提升；（c1）為演算法2的輸出結果，但由於亮度變化不均勻，出現了很多不同亮度級的小塊，與原圖相比亮度過度增強；（d1）為原圖和演算法1疊加後的輸出，與原圖（a1）相比，細節得到適度增強，比如影象中的眼睛更清晰，亮度的適度提升導致影象的對比度提高明顯，視覺效果更佳；（e1）為原圖和演算法2疊加後的輸出，與原圖（a1）相比，由於亮度的過度增強，影象（e1）出現了光暈，鼻子和嘴不能呈現出來，且由於亮度變化不連續，出現了分塊。
三 總結
對新的一篇論文進行了復現，並和上一篇的實驗效果作為比較，本文演算法的實驗效果不如演算法1的效果，故應用在工業上，可能優先考慮演算法1。上一篇論文復現用了整整一週，這次的論文復現僅僅三天，有做過的經驗和基礎是很重要，萬事開頭難。第一天精讀論文三遍，理解演算法原理。第二天用matlab程式碼實現了一半的演算法，第三天結合實現演算法1的思路，直接套用在演算法2上。
這次的感悟還是準確理解演算法原理是最重要的，程式碼就是描述而已。