【16位RAW影像處理四】 HDR->LDR：Matlab中tonemapfarbman函式的解析和自我實現。

　　最近受朋友的委託，想自己實現Matlab裡的一個HDR轉LDR的函式，函式名是tonemapfarbman，乘著十一假期，稍微瀏覽下這個函式，並做了一點C++的實現和最佳化。

　　為了看到這個函式的效果，需要至少matlab R2018b及其以上的版本。

　　首先，我們下載了matlab幫助文件中提到的該演算法對應的論文：Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation.

　　通讀了下這個論文，感覺他並沒有提出一個什麼具體的演算法，而是一個大的框架，其最有用的部分為 3.1 Multi-scale edge-preserving decompositions。

　　

　　這個意思呢，就是用某一種保邊濾波器進行K+1層次的分解，但是分解是依次的，這K+1次保邊濾波器使用某一種逐漸增強的引數，使得細節資訊呢越來越少，但是邊緣資訊還是儘量儲存。這個時候，前一層的結果減去後一層結果即為細節層，最後依次的保邊濾波結果即為基礎層，此時，由基礎層再反向加上每一層的細節層即可以得到原始一模一樣的資料。

　　這個呢，和早期的16位RAW影像處理一】：基於Fast Bilateral Filtering 演算法的 High-Dynamic Range(HDR) 影像顯示技術 相比，其實最大的區別就在於用了多層保邊濾波器，而這個文章只用了一層。

　　同一些基於金字塔的的HDR相比呢，其核心區別是沒有進行任何的下采樣，而都是全圖進行的。

　　那麼要進行HDR到LDR的操作，總的來說就是一個要點，想辦法調節細節層的資訊，一個最簡單的方法就是每個細節層乘以一個係數，在文章後續還描述了一些過程，針對不同的需求有不同的實現方式，不過我們去翻看matlab針對這個函式的核心程式碼，發現m處理的是相對來說非常簡單的，我們摘取下他的核心部分程式碼：

% Compute log luminance
logLum = log(lum);

% Apply diffusion filter for multi-scale decomposition
% Compute detail layers
% Recombine the layers together while moderately boosting multiple scale detail
uPre = logLum;
comLogLum = zeros(size(HDR,1),size(HDR,2), 'like', HDR);
numItr = 5; % No. of iterations ('NumberOfIterations' of anisotropic diffusion) for the first level decomposition.
logLum(~isfinite(logLum)) = (realmax('single')/100); % replaced Inf with a high value.
rangeLum = sum([max(logLum(:)) - min(logLum(:)),eps('single')],'omitnan'); % range of the image
for scaleInd = 1:numScale
    % gradThresh is the 'GradientThreshold' parameter of the anisotropic diffusion.
    % Here, it is taken as 5% of the dynamic range of the image.
    % gradThresh is increased for each decomposition level (multiplied by scaleInd (No. of iterations))
    gradThresh = scaleInd*5*(rangeLum)/100;
    uCurr = imdiffusefilt(logLum,'NumberOfIterations',numItr*scaleInd,'GradientThreshold',gradThresh);
    detail = uPre - uCurr; % multi-scale decomposition (detail extraction)
    comLogLum = comLogLum + weight(scaleInd).*detail; % weighted summation
    uPre = uCurr;
end
comLogLum = comLogLum + uCurr;

% Convert back to RGB
% The 99.9% intensities of the compressed image are remapped to give the output image a consistent look
newI = exp(comLogLum);
sI = sort(newI(:));
mx = sI(round(length(sI) * (99.9/100)));
newI = newI/mx;

% Exposure, RangeCompression and saturation correction
expoCorrect = exposure*newI;
if size(HDR,3) == 1
    LDR = expoCorrect.^gamma;
else
    LDR = (expoCorrect .* (rgb .^ sat)).^gamma;
end
LDR = im2uint8(LDR);
end

　　第二行程式碼，將影像資料轉換到log空間，這基本上是HDR演算法第一步的標準做法。

　　從第12行到第22行是演算法的核心部分，在這個迴圈裡，使用了imdiffusefilt這個函式作為保邊濾波器，他實際上是多次各向異性濾波器的迭代版本呢，這個濾波器具有梯度閾值和迭代次數兩個引數，迴圈中，迭代次數隨著迴圈的增加線性增加，梯度閾值也在每次迭代時做相應的調整，從而得到一個逐漸模糊且保邊影像，如下圖所示：

　　　　　　　　原圖　　　　　　　　　　　　　　 GradientThreshold = 12，NumberOfIterations = 5　　　　GradientThreshold = 24，NumberOfIterations = 10　　GradientThreshold = 36，NumberOfIterations = 15

　　第18行使用detail = uPre - uCurr，即前一次的保邊結果-本次的保邊濾波結果得到這一層的細節資訊，然後第19行

　　　　　　comLogLum = comLogLum + weight(scaleInd).*detail;

　　把detail資訊乘以使用者指定的增強係數(大於1，增加本層的細節，小於1，減少本層的細節），在累加到之前的細節中去。

　　第22行 comLogLum = comLogLum + uCurr; 中，此時的uCurr中儲存了最後一次保邊濾波器的結算結果，所以把他加入到前面的細節資訊中接得到我們處理後的結果。

　　第26行把對數空間的資料透過exp指令再次恢復到正常的空間，其實此時配合上im2uint8就應該能得到最後的LDR影像了，但是實際上這個時候影像的細節資訊基本已經得到了增強，但是整體的可視度或者視覺效果是很一般的，所以後面再透過曝光度和Gamma校正兩個引數適當調整輸出的效果。

　　第27到29行主要是取恢復後的資料前99.9%的值作為閾值，把那些過渡曝光的點消除掉，後續程式碼再進行exposure和gamma調整。

　　整個流程也沒有什麼特別複雜的地方。

　　翻譯這個函式成C++並不是一件很複雜的事情，主要是imdiffusefilt的翻譯，對於二維的資料，這個函式呼叫了anisotropicDiffusion2D函式，具體檢視這個函式呢，他有4領域和8領域的方式，以領域為例，其核心程式碼如下：

 % DiffusionRate is fixed to 1/4 because we considered nearest neighbour
        % differences in 4 directions(East,West,North,South)
        diffusionRate = 1/4;
        diffImgNorth = paddedImg(1:end-1,2:end-1) - paddedImg(2:end,2:end-1);
        diffImgEast = paddedImg(2:end-1,2:end) - paddedImg(2:end-1,1:end-1);
        switch conductionMethod
            % Conduction coefficients
            case 'exponential'
                conductCoeffNorth = exp(-(abs(diffImgNorth)/gradientThreshold).^2);
                conductCoeffEast = exp(-(abs(diffImgEast)/gradientThreshold).^2);
            case 'quadratic'
                conductCoeffNorth = 1./(1+(abs(diffImgNorth)/gradientThreshold).^2);
                conductCoeffEast = 1./(1+(abs(diffImgEast)/gradientThreshold).^2);
        end
        fluxNorth = conductCoeffNorth .* diffImgNorth;
        fluxEast =  conductCoeffEast .* diffImgEast;
        
        % Discrete PDE solution
        I = I + diffusionRate * (fluxNorth(1:end-1,:) - fluxNorth(2:end,:) + ...
            fluxEast(:,2:end) - fluxEast(:,1:end-1));

　　裡面的大部分計算還是exp函式，然後就是涉及到了3*3的臨域，這個建議不要直接翻譯，因為matlab程式碼的向量化很厲害，要先理解他的意思，然後再重新寫。

　　我載入一副1700*3700左右的單通道16點陣圖像，在matlab中測試，使用預設引數（3層），處理的時間大概需要0.6s，個人認為這個速度相對來說是非常快的，因為這個演算法內部涉及到了太多浮點計算，特別是exp函式，我初步的C++版本的速度要比matlab的慢很多，後面經過SSE指令最佳化後，也需要1100ms，後續測試發現matlab裡的程式碼使用了多執行緒，而我這個是單執行緒的版本，如果統計用多執行緒，我這個可以做到300ms。

　　進一步的最佳化手段有，修改exp的實現，用近似的版本，比如使用快速exp演算法 這裡的迭代版本，可以由如下的程式碼實現：

inline __m128 _mm_myexp_ps(__m128 x)
{
    __m128 T = _mm_add_ps(_mm_set1_ps(1.0f), _mm_mul_ps(x, _mm_set1_ps(1.0f / 256)));
    __m128 TT = _mm_mul_ps(T, T);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    return TT;
}

　　這個要比標準的exp快很多，而精度基本差不多，但是單執行緒情況速度基本就可以做到250ms了。

　　我把這個演算法也整合到我的DEMO中了，引數介面如下所示：

　　

　　個人感覺這個函式也不是特別通用，對部分圖還是要仔細調整引數才能得到較為合理的結果。

　　關於16位RAW影像，本人開發了一個簡易的增強和處理程式，可在 https://files.cnblogs.com/files/Imageshop/Optimization_Demo_16.rar下載測試。

　　如果想時刻關注本人的最新文章，也可關注公眾號：

翻譯

搜尋

複製