＂Recovering High Dynamic Range Radiance Maps from Photograph

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <iostream>
#include <vector>
#define zMin 0
#define zMax 255

using namespace std;
using namespace cv;

class HDRCreator
{
private:
int lamuta;
vector<Point2i> selectedPixels;
vector<double> lnDeltaT;
float weight[256];
vector<Mat> frames;
Mat zMatrix;
Mat g,lE;
public:
HDRCreator(const vector<Point2i> & p_selectedPixels,const vector<double> & p_lnDeltaT,
const vector<Mat> & p_frames,const int & p_lamuta);
void buildPixelsMap();
void buildWeightArray();
void solveResponseFunctionInv();
void buildLnRadianceMap(Mat &lnRandianceMap);
//void displaylnRadianceMap();
//
//void buildOneExposureRadiance(const int & j,Mat & radiance);
};

#include "HDRCreator.h"

HDRCreator::HDRCreator(const vector<Point2i> & p_selectedPixels,const vector<double> & p_lnDeltaT,
const vector<Mat> & p_frames,const int & p_lamuta)
{
selectedPixels = p_selectedPixels;
lnDeltaT = p_lnDeltaT;
frames = p_frames;
lamuta = p_lamuta;
}

void HDRCreator::buildPixelsMap()
{
zMatrix.create(selectedPixels.size(),frames.size(),CV_8UC1);
int i,j;
for (i = 0;i < zMatrix.rows;i ++)
{
for (j = 0;j < zMatrix.cols;j ++)
{
//特別要注意這裡的坐標是width * hight
zMatrix.at<uchar>(i,j) = frames[j].at<uchar>(selectedPixels[i].y,selectedPixels[i].x);
}
}
imshow("zMatrix",zMatrix);
waitKey(10000);
}

void HDRCreator::buildWeightArray()
{
float middle = 1.0 / 2.0 * (zMax + zMin);
for (int z = 0;z < 256;z++)
{
if (z <= middle)
{
weight[z] = z - zMin;
}
else
{
weight[z] = zMax - z;
}
}
}

void HDRCreator::solveResponseFunctionInv()
{
int n = 256;
Mat A = Mat::zeros(zMatrix.rows * zMatrix.cols + n - 1,n + zMatrix.rows,CV_64FC1);
Mat b = Mat::zeros(A.rows,1,CV_64FC1);
Mat X = Mat::zeros(A.cols,1,CV_64FC1);

int k = 0;
for (int i = 0;i < zMatrix.rows;i ++)
{
for (int j = 0;j < zMatrix.cols;j ++)
{
uchar zij = zMatrix.at<uchar>(i,j);
float wij = weight[zij];
A.at<double>(k,zij) = wij;
A.at<double>(k,n + i) = - wij;
b.at<double>(k,1) = wij * lnDeltaT[j];
k ++;
}
}
A.at<double>(k,128) = 1;
k ++;
for (int z = 1;z < n - 1;z ++)
{
A.at<double>(k,z - 1) = lamuta * weight[z];
A.at<double>(k,z + 0) = lamuta * weight[z] * -2;
A.at<double>(k,z + 1) = lamuta * weight[z];
k ++;
}

solve(A,b,X,CV_SVD);
g.create(n,1,CV_64FC1);
lE.create(A.rows - n,1,CV_64FC1);
g = X.rowRange(0,n).clone();
//exp(g,g);
lE = X.rowRange(n,X.rows).clone();
}

void HDRCreator::buildLnRadianceMap(Mat & lnRadianceMap)
{
buildWeightArray();
buildPixelsMap();
solveResponseFunctionInv();

lnRadianceMap.create(frames[0].rows,frames[0].cols,CV_64FC1);

for (int x = 0;x < lnRadianceMap.rows;x ++)
{
for (int y = 0;y < lnRadianceMap.cols;y ++)
{
double sumDenominator = 0.0;
double sumNominator = 0.0;
for (int j = 0;j < frames.size();j ++)
{
sumDenominator += weight[frames[j].at<uchar>(x,y)];
}

for (int j = 0;j < frames.size();j ++)
{
sumNominator += weight[frames[j].at<uchar>(x,y)] * (g.at<double>(frames[j].at<uchar>(x,y),0) - lnDeltaT[j]);
}
lnRadianceMap.at<double>(x,y) = sumNominator / sumDenominator;
}
}

imshow("lnradiance",lnRadianceMap);
waitKey(10000);
}