"Recovering High Dynamic Range Radiance Maps from Photograph

# Include <opencv2/core. hpp>
# Include <opencv2/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 ++)
{
// Note that the coordinates here are width * hight.
ZMatrix. at <uchar> (I, j) = frames [j]. at <uchar> (selectedPixels [I]. y, selectedPixels [I]. x );
}
}
Imshow ("zMatrix", zMatrix );
WaitKey (1, 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 (1, 10000 );
}