C / C++ HDR rgbe 讀檔與色調映射(toneMapping) 實作範例代碼

以前做的作業整理起來順便發表上來，使用的方法是最簡單的能看能主的方法，並沒有任何優化。

參考論文：

細節與對比強化之高動態範圍影像顯示方法
High Dynamic Range Image Display with Detail and Contrast Enhancement

本篇不是實現上述論文，實現的是他提到的先前的方法
基於這個公式做色調映射：

色調映射有分全域與局域(本文實作的是全域)，全域就是同一個公式對整張圖的每個點做，比較簡單迅速，缺點是可能會有某些地方表現不好，通常是最亮或是最暗的地方。需要比較好看的話要做局域色調映射，針對特別亮或暗的地方做不同處理。

使用函式庫：

rgbe
OpenBMP

流程

讀取 rgbe 轉換為 rgb
轉換彩色模型至 Yxy
對亮度通道做色調映射(tone mapping)
色彩模型轉回 rgb
進行 gama 校正 (微軟 Win10 設置為 2.2)
完成 global tone mapping

讀取 rgbe 轉換為 rgb

使用函式庫 rgbe 讀取，如果要重刻過程其實還蠻麻煩的~

使用的資料結構

struct basic_rgbeData {
    int width;
    int height;
    rgbe_header_info info;
    vector<float> img;
};

讀取的程式

void rgbeData_read(basic_rgbeData& hdr, string name) {
    FILE* hdrFile;
    fopen_s(&hdrFile, name.c_str(), "rb");

    // 讀檔頭
    RGBE_ReadHeader(hdrFile, &hdr.width, &hdr.height, &hdr.info);
    // 要求空間
    hdr.img.resize(hdr.width*hdr.height*3);
    // 讀rgbe
    RGBE_ReadPixels_RLE(hdrFile, hdr.img.data(), hdr.width, hdr.height);

    fclose(hdrFile);
}

讀進來之後就可以獲得圖像的長寬與 rgb 值了，這裡的rgb值會是小數點介於 0~1 之間，可以直接全部*255，然後把圖片顯示出來就可以看到黑黑的圖像了。

大於255的就變成255，小於0的就變成0

轉換彩色模型至 Yxy

把色彩模型從 rgb 轉換至 Yxy，實際上轉了兩次 rgbe -> yxz -> Yxy，因為公式可以合併並不衝突，我把它寫在一起了。(Yxy只是將ZXY轉換到坐標系上)
這個做是為了分離出光源與彩度，需要做處理的是光源，彩度則維持原本的樣子。

XYZ公式可以參考維基百科：https://www.wikiwand.com/zh-tw/CIE1931%E8%89%B2%E5%BD%A9%E7%A9%BA%E9%97%B4

void rgb2Yxy(const float* src, float* dst, int size) {
//#pragma omp parallel for
    for(int i = 0; i < size; ++i) {
        float a, b, c;
        a = (0.412453) * src[i*3 + 0]+
            (0.357580) * src[i*3 + 1]+
            (0.180423) * src[i*3 + 2];

        b = (0.212671) * src[i*3 + 0]+
            (0.715160) * src[i*3 + 1]+
            (0.072169) * src[i*3 + 2];

        c = (0.019334) * src[i*3 + 0]+
            (0.119193) * src[i*3 + 1]+
            (0.950227) * src[i*3 + 2];
        dst[i*3 + 0] = b;
        dst[i*3 + 1] = a / (a+b+c);
        dst[i*3 + 2] = b / (a+b+c);
    }
}

#pragma omp parallel for 是平行算的標記，下面接著的for迴圈會被展開來平行運算。

對亮度通道做色調映射(tone mapping)

這裡就是全文的重點了，來自於上面介紹的公式。

其中的 dmax 與 b 來源是依據論文所提的建議設置
size是圖像的指標dst的總長度

#define HDR_dmax 100.0
#define HDR_b 0.85

void globalToneMapping(float* dst, int size, float dmax, float b)
{
    constexpr int dim = 3; // 幾個通道
    constexpr int rgb = 0; // 選擇哪個通道

    float maxLum = dst[rgb];
    for(unsigned i = 1; i < size; ++i) {
        if(dst[i*dim+rgb] > maxLum)
            maxLum = dst[i*dim+rgb];
    }

    float logSum = 0.0;
    for(int i = 0; i < size; ++i) 
        logSum += log(dst[i*dim+rgb]);

    const float logAvgLum = logSum/size;
    const float avgLum = exp(logAvgLum);
    const float maxLumW = (maxLum / avgLum);
    const float coeff = (dmax*float(0.01)) / log10(maxLumW+1.0);

#pragma omp parallel for
    for(int i = 0; i < size; ++i) {
        auto& p = dst[i*dim+rgb];
        p /= avgLum;
        p = log(p+1.0) / log(2.0 + pow((p/maxLumW),(log(b)/log(0.5)))*8.0);
        p *= coeff;
    }
}

色彩模型轉回 rgb

跟剛剛的公式一樣只是反轉回來RGB模型

void Yxy2rgb(const float* src, float* dst, int size){
//#pragma omp parallel for
    for(int i = 0; i < size; ++i) {
        float a, b, c, newW;
        newW = src[i*3 + 0] / src[i*3 + 2];
        a = src[i*3 + 0];
        b = newW * src[i*3 + 1];
        c = newW - b - a;

        dst[i*3 + 0]  =  float(3.240479)*b;
        dst[i*3 + 0] += -float(1.537150)*a;
        dst[i*3 + 0] += -float(0.498535)*c;

        dst[i*3 + 1]  = -float(0.969256)*b;
        dst[i*3 + 1] +=  float(1.875992)*a;
        dst[i*3 + 1] +=  float(0.041556)*c;

        dst[i*3 + 2]  =  float(0.055648)*b;
        dst[i*3 + 2] += -float(0.204043)*a;
        dst[i*3 + 2] +=  float(1.057311)*c;
    }
}

轉回來之後就已經映射完畢可以輸出來看看，與一開始相比過亮與過暗的地方已經修正完畢可以看到細節，但是整體偏暗許多需要再進行下一步修正。

進行 gama 校正

最後要將亮度調整到作業系統適合的參數，圖像一般都在 1.0 ，而Windwos系統使用的是 2.2 要將它修正到正確的參數即可。

gama參數這個資訊可以從結構裡面我們一直沒用到參數 rgbe_header_info info; 裡獲得。

#define HDR_gama 2.2

void gama_fix(float* dst, int size, float gam) {
    const float fgamma = (0.45/gam)*2.0;
    float slope = 4.5;
    float start = 0.018;
    // 判定係數
    if(gam >= float(2.1)) {
        start /= ((gam - 2) * float(7.5));
        slope *= ((gam - 2) * float(7.5));
    } else if(gam <= float(1.9)) {
        start *= ((2 - gam) * float(7.5));
        slope /= ((2 - gam) * float(7.5));
    }
    // 校正像素
#pragma omp parallel for
    for (int i = 0; i < size*3; i++) {
        if(dst[i] <= start) {
            dst[i] = dst[i]*slope;
        } else {
            dst[i] = float(1.099)*pow(dst[i], fgamma) - float(0.099);
        }
    }
}

完成 global tone mapping

再來結構的裡的圖像就可以使用摟，把它印出來。
轉出來的數值還是 0~1 要把他們 *255

然後大於255的就變成255，小於0的就變成0

hdr.bmp("resultIMG/HDR_IMG.bmp");

專案原始碼

未封裝測試代碼可以參考 OpenHRD.cpp 中的 testMapping() 如下，比較容易讀懂在寫什麼。

rgbeData_writeBMP 是將 float* 輸出成 *.bmp 圖檔到硬碟上使用函式庫 OpenBMP

void testMapping(string name) {
    // read file
    basic_rgbeData hdr;
    rgbeData_read(hdr, name);
    rgbeData_info(hdr);
    //rgbeData_writeBMP(hdr, "resultIMG\HDR_non.bmp");

    // Mpping
    int imgSize = rgbeData_size(hdr);
    vector<float> Yxy(imgSize*3);

    rgb2Yxy(hdr.img.data(), Yxy.data(), imgSize);
    globalToneMapping(Yxy.data(), imgSize);
    Yxy2rgb(Yxy.data(), hdr.img.data(), imgSize);
    //rgbeData_writeBMP(hdr, "resultIMG\HDR_mapping.bmp");

    cout << hdr.info.gamma << endl;
    gama_fix(hdr.img.data(), imgSize, 2.2);
    rgbeData_writeBMP(hdr, "resultIMG/HDR_IMG.bmp");
}

這裡是完整的專案已經封裝成完整的物件：OpenHDR for VS2017

CHG

2018年6月6日星期三