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；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
（二）暗區域增強
通過第一步的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
（三）空間濾波器處理
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);
實驗效果如下：
實驗現象簡單分析，演算法處理後，輸出的影像變暗很多，即亮度改變很大，但視覺效果差，和原圖差很遠，即細節增強效果不好。
本文演算法主要是對影像的亮度調整，故關鍵點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：自己編程式碼。很耗時，同樣很多問題伴隨，這就是工科生，在問題中收穫知識和成長！
瞭解世界和自己！