基于C++的本征图像分解(Intrinsic Image Decomposition)

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

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

//main.cpp#include "lum_retinex.h"
#include <opencv2/opencv.hpp>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 <Eigen/SparseCore>
#include <Eigen/SparseCholesky>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 <Eigen/SparseCore>
#include <Eigen/SparseCholesky>
#include <chrono>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.\n");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), -1));entries.push_back(Triplet(cit, I.reflidx(x, y + 1), +1));cit++;entries.push_back(Triplet(cit, I.shadidx(x, y), -1));entries.push_back(Triplet(cit, I.shadidx(x, y + 1), +1));cit++;// dyR(r, c) = -lR(r, c) + lR(r + 1)// dyS(r, c) = -lS(r, c) + lS(r + 1)}// reflectance plus shading (log space) == final imageentries.push_back(Triplet(cit, I.reflidx(x, y), 1));entries.push_back(Triplet(cit, I.shadidx(x, y), 1));cit++;}}assert(entries.size() == nentries(w, h));double assemble_end = get_s();printf("Makemtx took %.1fms\n", (assemble_end- assemble_start) * 1000);return entries;}// to log domainvoid preprocess(int w, int h, const float* img, float* logimg) {printf("Start Preprocessing.\n");double preprocess_start = get_s();log(img, logimg, w*h);double preprocess_b_end = get_s();printf("Preprocess took %.1fms\n", (preprocess_b_end - preprocess_start) * 1000);}// back to linearvoid postprocess(int w, int h, float* refl_in, float* shading_in, float* refl_out, float* shading_out) {printf("Post process start.\n");double postprocess_start = get_s();reflect_clamp(w, h, refl_in, shading_in);exp(refl_in, refl_out, w*h);exp(shading_in, shading_out, w*h);double postprocess_end = get_s();printf("Postprocess took %.1fms\n", (postprocess_end - postprocess_start)*1000 );}Eigen::VectorXf makeB(float threshold, const float* im, int w, int h) {Eigen::VectorXf b(nconstraints(w, h));double assemble_b_start = get_s();imwrap I(im, w, h);int cit = 0;for (int y = 0; y < h; ++y) {for (int x = 0; x < w; ++x) {if (x < w - 1) {float dx = -I(x, y) + I(x + 1, y);float dxR;float dxS;if (std::abs(dx) > threshold) {dxR = dx;dxS = 0;} else {dxR = 0;dxS = dx;}// dxR(r, c) = -lR(r, c) + lR(r, c + 1)b(cit++) = dxR;// dxS(r, c) = -lS(r, c) + lS(r, c + 1)b(cit++) = dxS;}if (y < h - 1) {float dy = -I(x, y) + I(x, y + 1);float dyR;float dyS;if (std::abs(dy) > threshold) {dyR = dy;dyS = 0;} else {dyR = 0;dyS = dy;}b(cit++) = dyR;b(cit++) = dyS;// dyR(r, c) = -lR(r, c) + lR(r + 1)// dyS(r, c) = -lS(r, c) + lS(r + 1)}// reflectance plus shading (log space) == final imageb(cit++) = I(x, y);}}return b;}// operates on log-reflectance// makes sure that log-reflectane is less than 0void reflect_clamp(int w, int h, float* refl_in, float* shading_in) {float max_reflectance = -FLT_MIN;float min_reflectance = FLT_MAX;int nancount = 0;int infcount = 0;for (int i = 0; i < w * h; ++i) {if (refl_in[i] > max_reflectance) {max_reflectance = refl_in[i];}if (refl_in[i] < min_reflectance) {min_reflectance = refl_in[i];}}if (max_reflectance > 0) {for (int i = 0; i < w * h; ++i) {refl_in[i] -= max_reflectance;}for (int i = 0; i < w * h; ++i) {shading_in[i] += max_reflectance;}}}void retinex(float threshold, const float* im, int w, int h, float* reflectance, float* shading) {assert(reflectance);assert(shading);assert(im);retinex_decomp rdec(w, h);rdec.solve(threshold, im, reflectance, shading);}/* I would prefer solving direction Ax = b.However, this doesn't work with Eigen's QR decomp (why?).I solve A'Ax = A'b, using Cholesky, instead.*/retinex_decomp::retinex_decomp(int w, int h): m_w(w), m_h(h) {std::vector <Triplet> entries = makeTriplets(w, h);SpMat A(nconstraints(w, h), w * h * 2);A.setFromTriplets(entries.begin(), entries.end());m_At = A.transpose();{printf("factorize %d-by-%d matrix\n", A.cols(), A.cols());double decompose_start = get_s();m_solver.compute(m_At * A);double decompose_end = get_s();printf("factorize took %.1fms\n", (decompose_end - decompose_start) * 1000);}}void retinex_decomp::solve(float threshold, const float* im, float* reflectance, float* shading) {assert(reflectance);assert(shading);assert(im);std::vector<float> logimg(m_w*m_h);preprocess(m_w, m_h, im, logimg.data());Eigen::VectorXf b = makeB(threshold, logimg.data(), m_w, m_h);double solve_start = get_s();Eigen::VectorXf x = m_solver.solve(m_At * b);double solve_end = get_s();printf("solve took %.1fms\n", (solve_end - solve_start) * 1000);postprocess(m_w, m_h, x.data(), x.data() + m_w * m_h, reflectance, shading);}// does greyscale retinex with post processing.void retinex_decomp::solve_rgb(float threshold, const float* im, float* reflectance, float* shading) {assert(reflectance);assert(shading);assert(im);const int sz = m_w * m_h;std::vector<float> grayimg(sz);mean(im, grayimg.data(), grayimg.size());std::vector<float> graylog(sz);double before = get_s();log(grayimg.data(), graylog.data(), sz);double after = get_s();Eigen::VectorXf b = makeB(threshold, graylog.data(), m_w, m_h);double solve_start = get_s();Eigen::VectorXf x = m_solver.solve(m_At * b);double solve_end = get_s();float* log_shading = x.data() + sz;reflect_clamp(m_w, m_h, x.data(), log_shading);std::vector<float> rgb_logimg(3*sz);log(im, rgb_logimg.data(), rgb_logimg.size());// do gray to rgb conversion// log R = log I - log Sfor (int i = 0; i < sz; ++i) {rgb_logimg[i * 3 + 0] = rgb_logimg[i * 3 + 0] - log_shading[i];rgb_logimg[i * 3 + 1] = rgb_logimg[i * 3 + 1] - log_shading[i];rgb_logimg[i * 3 + 2] = rgb_logimg[i * 3 + 2] - log_shading[i];}exp(rgb_logimg.data(), reflectance, 3 * sz);exp(log_shading, shading, sz);}
}

基于C++的本征图像分解(Intrinsic Image Decomposition)相关推荐

【论文笔记】—本征图像分解—Unsupervised—USI^3^D—2020-CVPR
[论文介绍] 提出了第一个基于物理的单图像无监督学习用于本征图像分解网络USI3D(Unsupervised Single Image Decomposition) 本征图像,是指将一幅图像分解成两个 ...
本征图像分解：Retinex理论【转载】
本征图像分解:Retinex理论鲜橙关注 0.1742019.07.29 00:42:51字数 1,504阅读 2,531 我们通过眼睛观察到或者相机拍摄到的物体颜色主要由两方面因素决定,第一是物体 ...
本征时间尺度分解（Intrinsic Time-Scale Decomposition）
1.在诊断领域应用. 本征时间尺度分解能够很好地处理振动信号,分解出振动信号最具代表性的特征信号,在如滚动轴承故障诊断上,采取本征时间尺度分解作为信号处理手段来进行故障诊断能够大大提高诊断性能, 2. ...
深度学习图像融合_基于深度学习的图像超分辨率最新进展与趋势【附PDF】
因PDF资源在微信公众号关注公众号:人工智能前沿讲习回复"超分辨"获取文章PDF 1.主题简介图像超分辨率是计算机视觉和图像处理领域一个非常重要的研究问题,在医疗图像分析.生物特 ...
基于深度学习的图像超分辨率方法总结
基于深度学习的SR方法懒得总结,就从一篇综述中选取了一部分基于深度学习的图像超分辨率方法. 原文:基于深度学习的图像超分辨率复原研究进展作者:孙旭李晓光李嘉锋卓力北京工业大学信号与信息处理 ...
基于强化学习的图像配准 - Image Registration: Reinforcement Learning Approaches
配准定义给定参考图像 I_f 和浮动图像 I_m ,所谓的配准就是寻找一个图像变换T,将浮动图像I_m变换到和 I_f 相同的坐标空间下,使得两个图像中对应的点处于同一坐标下,从而达到信息聚合的目的 ...
基于小波的图像边缘检测,小波变换边缘检测原理
1.什么是"小波神经网络"?能干什么用呀小波神经网络(Wavelet Neural Network, WNN)是在小波分析研究获得突破的基础上提出的一种人工神经网络.它是基于小波 ...
基于深度学习的图像语义分析及其应用
本文转自"火光摇曳"博客:语义分析的一些方法(三),主要论述了基于深度学习方法的图像语义分析,包括图片分类.图片搜索.图片标注(image2text.image2sentence ...
一种基于卷积神经网络的图像去雾研究-含matlab代码
目录一.绪论二.去雾卷积网络 2.1 特征提取 2.2 多尺度映射 2.3 局部均值 2.4 非线性回归三.实验与分析四.Matlab代码获取一.绪论雾是一种常见的大气现象,空气中悬浮的水 ...

基于C++的本征图像分解(Intrinsic Image Decomposition)

基于C++的本征图像分解(Intrinsic Image Decomposition)相关推荐

最新文章

热门文章