Mat GrayHistMatching(Mat I, Mat R)
{
/* Histogram Matching of a gray image with a reference*/
// accept two images I (input image) and R (reference image)
Mat Result; // The Result image
int L = 256; // Establish the number of bins
if(I.channels()!=1)
{
cout<<"Please use Gray image"<<endl;
return Mat::zeros(I.size(),CV_8UC1);
}
Mat G,S,F; //G is the reference CDF, S the CDF of the equlized given image, F is the map from S->G
if(R.cols>1)
{
if(R.channels()!=1)
{
cout<<"Please use Gray reference"<<endl;
return Mat::zeros(I.size(),CV_8UC1);
}
Mat R_hist, Rp_hist; //R_hist the counts of pixels for each level, Rp_hist is the PDF of each level
/// Set the ranges ( for B,G,R) )
float range[] = { 0, 256 } ;
const float* histRange = { range };
bool uniform = true; bool accumulate = false;
//calcHist(const Mat* images, int nimages, const int* channels, InputArray mask, OutputArray hist, int dims, const int* histSize, const float** ranges, bool uniform=true, bool accumulate=false )
calcHist( &R, 1, 0, Mat(), R_hist, 1, &L, &histRange, uniform, accumulate );
//Calc PDF of the image
Rp_hist=R_hist/(I.rows*I.cols);
//calc G=(L-1)*CDF(p)
Mat CDF=cumsum(Rp_hist);
G=(L-1)*CDF;
for(int i=0;i<G.rows;i++)
G.at<Point2f>(i,0).x=(float)cvRound(G.at<Point2f>(i,0).x);//round G
}
else
{
//else, the given R is the reference PDF
Mat CDF=cumsum(R);
G=(L-1)*CDF;
for(int i=0;i<G.rows;i++)
G.at<Point2f>(i,0).x=(float)cvRound(G.at<Point2f>(i,0).x);//round G
}
/// Establish the number of bins
Mat S_hist, Sp_hist; //S_hist the counts of pixels for each level, Sp_hist is the PDF of each level
/// Set the ranges ( for B,G,R) )
float range[] = { 0, 256 } ;
const float* histRange = { range };
bool uniform = true; bool accumulate = false;
//calcHist(const Mat* images, int nimages, const int* channels, InputArray mask, OutputArray hist, int dims, const int* histSize, const float** ranges, bool uniform=true, bool accumulate=false )
calcHist( &I, 1, 0, Mat(), S_hist, 1, &L, &histRange, uniform, accumulate );
//Calc PDF of the image
Sp_hist=S_hist/(I.rows*I.cols);
//calc s=(L-1)*CDF(p)
Mat CDF=cumsum(Sp_hist);
S=(L-1)*CDF;
for(int i=0;i<S.rows;i++)
S.at<Point2f>(i,0).x=(float)cvRound(S.at<Point2f>(i,0).x);//round S
F=Mat::zeros(S.size(),CV_32F);
int minIndex=-1;
double T,min=100000;
for(int i=0;i<S.rows;i++)
{
for(int j=0;j<G.rows;j++)
{
T=abs(double(S.at<Point2f>(i,0).x)-double(G.at<Point2f>(j,0).x));
if (T==0)
{
minIndex=j;
break;
}
else
if(T<min)
{
minIndex=j;
min=T;
}
}
F.at<Point2f>(i,0).x=(float)minIndex;
minIndex=-1;
min=1000000;
}
uchar table[256];
for(int i=0;i<256;i++)
{
table[i]=(int)F.at<Point2f>(i,0).x;
}
Result= ScanImageAndReduceC(I, table);
return Result;
}
