這是我參與8月更文挑戰的第9天，活動詳情查看： 8月更文挑戰算法

正文

原理是基於這篇論文——《Adaptive Local Tone Mapping Based on Retinex for High Dynamic Range Images》markdown

論文提出的背景

雖然能夠從不一樣曝光的照片中得到包含真實場景全動態範圍的高動態範圍 (HDR) 圖像，可是普通顯示器等低動態範圍(LDR)顯示設備沒法處理場景的全動態範圍。LDR 設備只能顯示兩個數量級的。一旦將 HDR 圖像線性映射到顯示設備，就會丟失不少信息。所以，在映射到設備以前，必須對 HDR 圖像進行壓縮。從 HDR 圖像到 LDR 顯示設備的映射技術稱爲色調映射或色調再現。app

這些色調映射技術能夠分爲兩類：第一類是全局操做符，另外一類是局部操做符。koa

全局色調映射操做符是對全部像素應用一個函數；
局部色調映射算子根據相鄰像素對每一個像素應用不一樣的函數。

什麼是中心環繞視網膜

retinex 理論最初是由 Land 定義的。說明了人類視覺系統如何提取世界上可靠的顏色信息。Jobson 等人在中心/環繞視網膜的基礎上，引入了單尺度視網膜 (SSR) 和多尺度視網膜 (MSR) 。函數

SSR 的公式oop

其中 x,y 爲圖像像素座標， $R_{i}(x,y)$ 爲視網膜輸出， $I_{i}(x,y)$ 爲圖像在第 i 個光譜波段的分佈，*爲卷積運算，F(x,y) 爲高斯環繞函數post

其中 c 是高斯環繞空間常數。K 是歸一化因子（一個小的空間常數會產生良好的動態範圍壓縮，但顏色再現卻比較糟糕。相反，一個大的常數能產生良好的顏色再現，但動態範圍壓縮卻不太好）。性能

MSR 的公式優化

其中 N 是尺度的數量, $R_{ni}(x,y)$ 表示第 n 個尺度的第 i 個組成部分， $R_{MSR_{i}}(x,y)$ 是 MSR 輸出的第 i 個光譜份量， $\omega _{n}$ 是與第 n 個尺度相關的權重。 spa

MSR 的目標是減小高對比度邊緣周圍的暈圈僞影，並與動態範圍壓縮和顏色再現保持平衡。MSR 產生了良好的動態範圍壓縮，但仍然遭受光暈僞影的影響。此外，小空間常數的 SSR 使圖像中較大的均勻區域變灰變平。

如圖所示，

左圖是 MSR 的輸出，右圖是帶小空間常量的 SSR 的輸出

做者提出的算法

爲了解決以上的缺陷，Hyunchan Ahn 等人提出了一種新的局部色調映射方法，該方法基於中心/環繞視網膜，既能保留細節，又能防止光暈產生。

在算法中，亮度值從輸入的 HDR 圖像中得到並處理。首先，應用全局色調映射做爲預處理。在此基礎上，應用了基於視網膜算法的局部色調映射。最後，對處理後的亮度值和輸入的色度值進行歸一化，獲得輸出圖像。

全局色調映射

原理不復雜，這篇博客已經詳細闡述了——blog.csdn.net/just_sort/a…，我就再也不贅述。

根據博主的樸素實現（可是博主只實現了全局自適應部分的代碼，局部自適應並無實現）

我優化的 OpenCV 代碼以下：

int HDR(const cv::Mat &input_img, cv::Mat &out_img)
{
    int rows = input_img.rows;
    int cols = input_img.cols;

    cv::Mat src;
    input_img.convertTo(src, CV_32FC3);

    // -----------------------------------------------------------------------
    //Timer t;

    //Lw_max
    float Lw_max = 0.0;

    //Lw
    cv::Mat Lw = cv::Mat(rows, cols, CV_32FC1, cv::Scalar(0));

    for (int i = 0; i < rows; i++)
    {
        cv::Vec3f *ptr_vec3 = src.ptr<cv::Vec3f>(i);
        float *ptr = Lw.ptr<float>(i);

        for (int j = 0; j < cols; j++)
        {
            float tmp_0 = 0.299f * ptr_vec3[j][2];
            float tmp_1 = 0.587f * ptr_vec3[j][1];
            float tmp_2 = 0.114f * ptr_vec3[j][0];

            float val = tmp_0 + tmp_1 + tmp_2;

#if 1
            if (fabs(val) < 1e-5)
            {
                ptr[j] = 1;
            }
            else
            {
                ptr[j] = val;
            }
#else
            Lw.ptr<float>(i)[j] = val == 0 ? 1 : val;
#endif

            Lw_max = max(val, Lw_max); // 公式 4 中的 Lw_max
        }
    }
    //std::cout << "Lw_max; " << t.elapsed() << std::endl;

    // -----------------------------------------------------------------------

    //t.restart();
    // Lw_sum
    float Lw_sum = 0;
    for (int i = 0; i < rows; i++)
    {
        float *ptr = Lw.ptr<float>(i);
        for (int j = 0; j < cols; j++)
        {
            // .001f 其做用主要是爲了不對純黑色像素進行 log 計算時數值溢出
            float val = log(0.001f + ptr[j]);
            Lw_sum += val; // 公式 5 的求和部分
        }
    }

    //std::cout << "sum; " << t.elapsed() << std::endl;

    // -----------------------------------------------------------------------

    //t.restart();
    // Lwaver 公式 5 的倒數
    float inv_Lwaver = 1.f  / exp(Lw_sum / (rows * cols));

    for (int i = 0; i < rows; i++)
    {
        float *ptr = Lw.ptr<float>(i);
        cv::Vec3f *ptr_vec3 = src.ptr<cv::Vec3f>(i);

        //#pragma omp parallel for schedule (dynamic)
        for (int j = 0; j < cols; j++)
        {
            float val = ptr[j];
            float top = log(val * inv_Lwaver + 1);
            float bottom = log(Lw_max * inv_Lwaver + 1);
            float Lg = top / bottom;

            // 低照度的亮度部分比高照度的部分要能獲得更大程度的提高，因此對於低照度圖，該公式能起到很好的加強做用
            // 式中使用了全局的對數平均值，這就有了必定的自適應性。
            float gain = Lg / val;

            ptr_vec3[j][0] *= gain;
            ptr_vec3[j][1] *= gain;
            ptr_vec3[j][2] *= gain;
        }

    }
    //std::cout << "公式 5; " << t.elapsed() << std::endl;

    src.convertTo(out_img, CV_8UC3, 255.f);

    return 0;
}
複製代碼

執行的效果爲（圖片大小爲 894*1080，屏幕大小有限，截圖只截取了一部分）：

OpenCV 的版本執行（100次）時間平均爲 48.14ms 左右

通過 CUDA 優化以後的效果基本一致

可是執行時間大大減小，減小爲僅需 8.2 ms（執行 100 次，取平均，包括 host 和 device 之間的拷貝時間）

其餘圖片（也放縮成高度爲 1080，寬度保持縱橫比）效果：

7.55 ms

10.80 ms

基於視網膜算法的局部色調映射（個人補充）

基於視網膜理論的局部適應是在全局適應過程以後進行的。做者發現位於亮得多的像素值附近的輸入 $L_{l}(x,y)$ 的輸出 $R_{i}(x,y)$ 會變得很是暗，這將致使光暈僞影，使結果看起來不天然。經過引入邊緣保持濾波器，能夠減小僞影。所以做者採用導向濾波器代替高斯濾波器。導向濾波器是一種與雙邊濾波器相似的保邊濾波器，其權值不只取決於歐幾里德距離，還取決於亮度差。這些濾波器的性能類似，但導向濾波器在邊緣附近有更好的性能。此外，它的計算複雜度是線性時間的，沒有近似，與核大小無關。

局部適應方程可表示爲

式中， $L_{l}(x,y)$ 爲局部自適應的輸出， $H_{g}(x,y)$ 爲導向濾波器對 $L_{g}(x,y)$ 的輸出，

$\xi _{x}, \xi _{y}$ 鄰域像素座標, $\omega (x, y)$ 是一個（以像素爲 (x, y) 中心，以 r 爲半徑的）局部方形窗口，|ω| 是在 $\omega (x, y)$ 內像素的數量 $a(\xi _{x}, \xi _{y})$ ) 和 $b(\xi _{x}, \xi _{y})$ 是一些線性係數。

$\mu (\xi _{x}, \xi _{y})$ 和 $\sigma ^{2}(\xi_{x}, \xi _{y})$ 分別是 $L_{g}$ 在 $\omega (\xi _{x}, \xi _{y})$ 窗口下的的均值和方差， $\mu_{2}(\xi x, \xi y)$ 是 $L^{2}_{g}$ 在 $\omega (\xi _{x}, \xi _{y})$ 窗口下的均值，ε 是正則化參數。【由於在做者的算法中傳進導向濾波器的導向圖像與輸入圖像是一致的，因此公式看起來和何凱明的原公式有點不太同樣，實際上是一致的】

應用該濾波器後，光暈僞影明顯減小，但因爲其全局對比度較低，輸出的總體外觀仍然使人不夠滿意。

爲了防止濾波產生的平面現象，提升咱們的方法的性能，做者引入了兩個重要因素。

①、對比度加強因子 α

②、根據場景內容變化的自適應非線性偏移量 β

其中 η 表示對比控制參數， $L_{gmax}$ 是全局自適應輸出的最大亮度值。

其中 λ 是非線性控制參數, $\overline{L}_{g}$ 是全局自適應輸出的亮度對數均值。

最後獲得局部自適應的公式爲：

最終的效果

先通過全局自適應，再通過局部自適應後的結果以下所示

一種基於視神經網絡的高動態範圍（HDR）圖像自適應局部色調映射的實現【OpenCV】【CUDA】