利用本征图像分解(Intrinsic Image Decomposition)算法，将图像分解为shading(illumination) image 和 reflectance(albedo) image，计算图像的reflectance image。
Reflectance Image 是指在变化的光照条件下能够维持不变的图像部分。
Shading Image 是反应原图像光照情况的图像部分。

本文提供的算法是比较传统的，效果虽然不SOTA，但是可以满足要求不高的朋友，可以考虑用CUDA实现加速。

//main.cpp#include "lum_retinex.h"
#include const float threshold = 0.13;int main(int argc, char* argv[])
{ cv::Mat input = cv::imread("Input.png", 1);input.convertTo(input, CV_32FC3);int w = input.cols;int h = input.rows;cv::Mat reflectance(h, w, CV_32FC3);cv::Mat shading(h, w, CV_32FC1);lum::retinex_decomp rdecomp(w, h);rdecomp.solve_rgb(threshold, (const float*)input.data, (float*)reflectance.data, (float*)shading.data);cv::imshow("input", input);cv::imshow("shading", shading);cv::imshow("reflectance", reflectance);cv::Mat out(h, w, CV_32FC3);for (int i = 0; i < h; i++) { for (int j = 0; j < w; j++) { out.at<cv::Vec3f>(i, j)[0] = reflectance.at<cv::Vec3f>(i, j)[0] * 255.0;out.at<cv::Vec3f>(i, j)[1] = reflectance.at<cv::Vec3f>(i, j)[1] * 255.0;out.at<cv::Vec3f>(i, j)[2] = reflectance.at<cv::Vec3f>(i, j)[2] * 255.0;}}cv::imwrite("img/reflectance1.png", out);cv::waitKey(0);return 0;
}

//lum_retinex.h#ifndef LUM_RETINEX_H
#define LUM_RETINEX_H#include 
#include namespace lum { // Run retinex as single function call (cannot reuse decomposition).void retinex(float threshold, const float* img, int w, int h, float* reflectance, float* shading);// First create decomposition for linear solver (slow)// Then solve for various right-hand-sides (fast)class retinex_decomp { public:retinex_decomp(int w, int h);void solve(float threshold, const float* img, float* refl, float* shading);void solve_rgb(float threshold, const float* img, float* refl, float* shading);private:int m_w;int m_h;using SpMat = Eigen::SparseMatrix < float > ;SpMat m_At;Eigen::SimplicialCholesky<SpMat> m_solver;};
}
#endif

// retinextest.cpp #include "lum_retinex.h"
#include 
#include 
#include namespace lum { // timing helperdouble get_s() { using namespace std::chrono;auto now = system_clock::now();system_clock::duration tse = now.time_since_epoch();return duration_cast<nanoseconds>(tse).count() / 1e9;}// Helper image processing routinesvoid log(const float* inimg, float* outimg, int sz) { 
#ifdef USE_MKLvsLn(sz, inimg, outimg);
#elsefor (int i = 0; i < sz; ++i) { float in = inimg[i];outimg[i] = std::logf(inimg[i] + 0.00001);}
#endif}void exp(const float* inimg, float* outimg, int sz) { 
#ifdef USE_MKLvsExp(sz, inimg, outimg);
#elsefor (int i = 0; i < sz; ++i) { outimg[i] = std::expf(inimg[i]);}
#endif}void mean(const float* inimg, float* outimg, int outsz) { for (int i = 0; i < outsz; ++i) { outimg[i] = (inimg[i * 3] + inimg[i * 3 + 1] + inimg[i * 3 + 2]) / 3;}}// helpers for image indicesclass reflshadidx { public:reflshadidx(int w, int h):m_w(w), m_h(h){ }int reflidx(int x, int y) const { return m_w * y + x;}int shadidx(int x, int y) const { return m_w * m_h + reflidx(x, y);}private:int m_w;int m_h;};class imwrap { public:imwrap(const float* img, int w, int h):m_w(w), m_h(h), m_img(img){  }float operator()(int x, int y) const { assert(x >= 0);assert(y >= 0);assert(x < m_w);assert(y < m_h);return m_img[m_w * y + x];}private:const float* m_img;int m_w;int m_h;};// forward declaration of internal functionsvoid reflect_clamp(int w, int h, float* refl_in, float* shading_in);void preprocess(int w, int h, const float* img, float* logimg);void postprocess(int w, int h, float* refl_in, float* shading_in, float* refl_out, float* shading_out);Eigen::VectorXf makeB(float threshold, const float* im, int w, int h);using Triplet = Eigen::Triplet<float>;int nconstraints(int w, int h) { return w*h + 2 * w*(h - 1) + 2 * (w - 1) * h;}int nentries(int w, int h) { return 2 * (w*h + 2 * w*(h - 1) + 2 * (w - 1) * h);}std::vector<Triplet> makeTriplets(int w, int h) { reflshadidx I(w, h);printf("Assemble matrix.
");double assemble_start = get_s();std::vector <Triplet> entries;int cit = 0;for (int y = 0; y < h; ++y) { for (int x = 0; x < w; ++x) { if (x < w - 1) { // dxR(r, c) = -lR(r, c) + lR(r, c + 1)entries.push_back(Triplet(cit, I.reflidx(x, y), -1));entries.push_back(Triplet(cit, I.reflidx(x + 1, y), +1));cit++;// dxS(r, c) = -lS(r, c) + lS(r, c + 1)entries.push_back(Triplet(cit, I.shadidx(x, y), -1));entries.push_back(Triplet(cit, I.shadidx(x + 1, y), +1));cit++;}if (y < h - 1) { entries.push_back(Triplet(cit, I.reflidx(x, y