算法是基于Debevec在1997年的论文Recovering High Dynamic Range Radiance Maps from Photographs
之前在博客中也给出了几个用于学习理论知识的链接CRF理论学习

下面的C++代码是基于这个修改的。包括：

1. Image Irradiance的归一化
2. 指数函数拟合ICRF
3. 输出到本地文件

C++

#include <iostream>
#include "opencv2/highgui/highgui.hpp"    
#include "opencv2/opencv.hpp"    
#include <opencv2/core/core.hpp>  
#include <algorithm>
#include <stdlib.h>
#include <numeric>
#include <vector>
#include "Eigen/Dense"
#include "Eigen/Core"
#include <unsupported/Eigen/NonLinearOptimization>
#include <unsupported/Eigen/NumericalDiff>

using namespace std;
using namespace cv;

// Generic functor

template<typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic>
struct Functor
{
	typedef _Scalar Scalar;
	enum {
		InputsAtCompileTime = NX,
		ValuesAtCompileTime = NY
	};
	typedef Eigen::Matrix<Scalar, InputsAtCompileTime, 1> InputType;
	typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, 1> ValueType;
	typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, InputsAtCompileTime> JacobianType;

	int m_inputs, m_values;

	Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
	Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}

	int inputs() const { return m_inputs; }
	int values() const { return m_values; }
};

vector<float> allIrradiance;
struct FitICRF : Functor<float>
{
	// 输出个数必须大于输入个数, 故用2不用1;
	FitICRF(void) : Functor<float>(2, 2) {}

	int operator()(const Eigen::VectorXf &x, Eigen::VectorXf &y) const
	{
		float space = 1.0 / (255 + 1);
		y(0) = 0;
		for (float crf_X = 0.0; crf_X <= 1.0; crf_X = crf_X + space) {
			float crf_Y = pow(crf_X, 1.0 / x(0));
			int xids = crf_X / space;
			y(0) += abs(crf_Y - allIrradiance[xids]);
			//cout << Icrf << endl;
			y(1) = 0;
		}

		return 0;
	}

};

float FitFinalICRF() {
	Eigen::VectorXf x(2);
	x(0) = 1.0;
	FitICRF fitcurve;
	Eigen::NumericalDiff<FitICRF> numDiff(fitcurve);
	Eigen::LevenbergMarquardt<Eigen::NumericalDiff<FitICRF>, float> lm(numDiff);
	lm.parameters.maxfev = 100;
	lm.parameters.xtol = 1.0e-8;
	Eigen::VectorXf y(2);
	fitcurve.operator()(x, y);
	int iRet = lm.minimize(x);
	float icr_gamma = x(0);

	return icr_gamma;
}

void main() {

	vector<cv::Mat> Imgs;
	for (int i = 0; i < 3; i++) {
		cv::Mat img = cv::imread(to_string(i) + ".jpg");
		Imgs.emplace_back(img);
	}

	const vector<float> times = { (float)0.001479, (float)0.008333, (float)0.066667 };

	// 计算CRF
	Mat responseDebevec;  // response if a CV_32FC3 value.
	Ptr<CalibrateDebevec> calibrateDebevec = createCalibrateDebevec();
	calibrateDebevec->process(Imgs, responseDebevec, times);

	// 采用Debevec方法计算HDR
	Mat hdr;
	Ptr<MergeDebevec> mergeDebevec = createMergeDebevec();
	mergeDebevec->process(Imgs, hdr, times, responseDebevec);
	// map to ldr for draw/showing.
	Mat ldr;
	Ptr<Tonemap> tonemap = createTonemap(2.2f);		// Debevec算法需要经过一步色调映射
	tonemap->process(hdr, ldr);

	const int width_scale = 3, height_scale = 20;
	const int width = 255 * width_scale, height = 40 * height_scale;
	Mat img = Mat::zeros(Size(width, height), CV_8UC3);
	Point b_old = Point(0, height);
	Point g_old = Point(0, height);
	Point r_old = Point(0, height);

	for (int i = 0; i < 255; ++i) {
		Vec3f v3f = responseDebevec.at<Vec3f>(i, 0);
		float ib = v3f[0];
		float ig = v3f[1];
		float ir = v3f[2];

		float gray = (ib + ig + ir) / 3.0;
		allIrradiance.emplace_back(gray);
	}
	//归一化
	for (int j = 0; j < allIrradiance.size(); j++) {
		allIrradiance[j] = (allIrradiance[j] - allIrradiance[0]) / (allIrradiance[allIrradiance.size() - 1] - allIrradiance[0]);
	}

	float gamma = FitFinalICRF();
	cout << gamma << endl;
	ofstream file("./irradiance.txt");
	for (int i = 1; i <= 1024; i++) {
		float I = pow((float)i / 1024.0, 1.0 / gamma);
		file << setiosflags(ios::fixed)<<setprecision(5)<<I<< endl;//保留5位小数
	}
}

python

import cv2
import numpy as np
import matplotlib.pyplot as plt

# load exposure images into a list
img_fn = ["data/multiexpo/0.JPG", "data/multiexpo/1.JPG", "data/multiexpo/2.JPG"]
img_list = [cv2.imread(fn) for fn in img_fn]
exposure_times = np.array([0.001479, 0.008333, 0.05], dtype=np.float32)

# Estimate camera response function (CRF)
cal_debevec = cv2.createCalibrateDebevec()
crf_debevec = cal_debevec.process(img_list, times=exposure_times)

x = np.arange(256)

plt.plot(x, crf_debevec[:, 0, 0], c='r', marker='o')
plt.plot(x, crf_debevec[:, 0, 1], c='g', marker='o')
plt.plot(x, crf_debevec[:, 0, 2], c='b', marker='o')
plt.show()