Histogram Matching
===
[TOC]

Python version : https://stackoverflow.com/questions/32655686/histogram-matching-of-two-images-in-python-2-x

## The Concept of Histogram Matching
G1(z1), G2(z2), z1 represents Match image, G1 represents histogram eaulization of z1, z2 represents src
z3 (dst) = G1^(-1)G2(z2)

## Image read

// If don't know histogram equalization, please refer to https://codeinfo.space/imageprocessing/histogram-equalization/
Mat src = imread("lenna.jpg",0);
Mat match = imread("mountains.jpg",0);

## Get Histogram of Match Image

    for (int y = 0; y < match.rows; y++){
        for (int x = 0; x < match.cols; x++){
            uchar grey = *(match.data + y*match.cols + x);
            matchHist[grey]++;
            if(grey<=tmpMinGrey) {
                tmpMinGrey=grey;
            }
        }
    }

## Get CDF of Match Image

// Get CDF of match
    matchCDF[0] = matchHist[0];
    for (int i = 0; i <= 255; i++){
        matchCDF[i] = matchCDF[i-1] + matchHist[i];
    }

## Get Equalization ToDst table of mach image

// Get Equalization To Dst table of mach image
    tmpConst = 255.0f / (float)(match.rows*match.cols - matchHist[tmpMinGrey]);
    cout << "tmpConst = " << tmpConst << endl;
    for (int y = 0; y < match.rows; y++){
        for (int x = 0; x < match.cols; x++){
            uchar grey = *(match.data + y*match.cols + x);
            uchar equalized = (uchar)((float)(matchCDF[grey] - matchHist[tmpMinGrey])*tmpConst);
            eToDst[equalized] = grey;   // save inverse value  equalized(match)->grey(match)
        }
    }

## Get Histogram and CDF of src

// Get Histogram of src;
    for (int y = 0; y < src.rows; y++){
        for (int x = 0; x < src.cols; x++){
            uchar grey = *(src.data + y*src.cols + x);
            srcHist[grey]++;
            if (grey <= tmpMinGrey){
                tmpMinGrey = grey;
            }
        }
    }

    // Get CDF of src
    srcCDF[0] = srcHist[0];
    for (int i = 1; i <= 255; i++){
        srcCDF[i] = srcCDF[i-1] + srcHist[i];
    }

## Get Dst

   // Get dst
    tmpConst=255.0f/(float)(src.rows*src.cols-srcHist[tmpMinGrey]);
    for (int y = 0; y < dst.rows; y++){
        for (int x = 0; x < dst.cols; x++){
            uchar srcgrey = *(src.data+y*src.cols+x);
            uchar equalgrey = (uchar)((float)(srcCDF[srcgrey]-srcHist[tmpMinGrey])*tmpConst);
            *(dst.data+y*dst.cols+x) = getMatchGrey(equalgrey,eToDst);//通過查表得到結果
        }
    }

// get the matching grey 通過匹配圖逆變換的表，查到某灰度值最近接的對應值
// according to equalizetion To Dst table 
// and the input grey which is the equalization of src
uchar getMatchGrey(uchar input, uchar* table){  //input爲某灰度值，table爲逆變換的表
    uchar grey = *(table+input);
    if (grey != 0) return grey;
    // 如果沒有，就往當前值的兩邊尋找，尋找時注意灰度範圍的一頭一尾
    int right = 0, left = 0;
    uchar outNum=input;
    int i = 1, j = 1;
    while (1) {
        if (right != 0 && left != 0) break;
        if(*(table+input+i)!=0) {right=input+i;}
        if(*(table+input-j)!=0) {left=input-j;}
        i++;
        j++;
    }

    return (*(table+right) + *(table+left)) / 2;// return the colsest grey value
}




## How about RGB image?

#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;

// show histogram of RGB image with continuous curve
// 用於顯示彩色圖的直方圖
void myShowHistRGB(int psrccnt[][256],string winname,int rows=400,int cols=512)
{
    Mat hist(rows,cols,CV_8UC3,Scalar(188,188,188));

    int maxcnt=psrccnt[0][0];
    for(int k=0; k<3; k++){
        for(int i=0; i<=255; i++){
            if(psrccnt[k][i]>maxcnt) maxcnt=psrccnt[k][i];
        }
    }
    int thickness = 2;
    int lineType = 8;
    int k;
    int i;
    int start;// find where existed gery start
    // B
    k=0;
    i=0;
    start=0;
    while(psrccnt[k][i]==0)
    {
        start=i;
        i++;
    }
    for(int i=start+1; i<=255; i++) {
        int cnt=psrccnt[k][i];
        if(cnt!=0) {
            line(hist,
                Point((int)((float)(start)*(float)cols/256.0f),rows-(int)((float)psrccnt[k][start]*(float)rows/(float)maxcnt)),
                Point((int)((float)i*(float)cols/256.0f),rows-(int)((float)cnt*(float)rows/(float)maxcnt)),
                Scalar(255,0,0),thickness,lineType);
            start=i;
        }
    }
    // G
    k=1;
    i=0;
    start=0;
    while(psrccnt[k][i]==0)
    {
        start=i;
        i++;
    }
    for(int i=start+1; i<=255; i++) {
        int cnt=psrccnt[k][i];
        if(cnt!=0) {
            line(hist,
                Point((int)((float)(start)*(float)cols/256.0f),rows-(int)((float)psrccnt[k][start]*(float)rows/(float)maxcnt)),
                Point((int)((float)i*(float)cols/256.0f),rows-(int)((float)cnt*(float)rows/(float)maxcnt)),
                Scalar(0,255,0),thickness,lineType);
            start=i;
        }
    }
    // R
    k=2;
    i=0;
    start=0;
    while(psrccnt[k][i]==0)
    {
        start=i;
        i++;
    }
    for(int i=start+1; i<=255; i++) {
        int cnt=psrccnt[k][i];
        if(cnt!=0) {
            line(hist,
                Point((int)((float)(start)*(float)cols/256.0f),rows-(int)((float)psrccnt[k][start]*(float)rows/(float)maxcnt)),
                Point((int)((float)i*(float)cols/256.0f),rows-(int)((float)cnt*(float)rows/(float)maxcnt)),
                Scalar(0,0,255),thickness,lineType);
            start=i;
        }
    }
    namedWindow(winname);
    imshow(winname, hist);
    imwrite(winname+".jpg",hist);
}

// get the matching grey 
// according to equalizetion To Dst table 
// and the input grey which is the equalization of src
uchar getMatchGrey(uchar input, uchar *table) {
    uchar grey=table[input];
    if(grey!=0) return grey;// if there is one corresponding value ,return it.
    uchar outNum=input;
    for (int i=1;i<=255;i++) {
        if((input+i)>=255) {outNum=255; break;}
        if((input-i)<=0) {outNum=0; break;}
        if(table[input+i]!=0) {outNum=input+i;break;}
        if(table[input-i]!=0) {outNum=input-i;break;}
    }
    if(outNum==255) return 255;
    if(outNum==0) return 0;
    return table[outNum];// return the colsest grey value

}

int main( void)
{
    Mat src = imread("lenna.jpg",1);
    Mat match = imread("fruits.jpg",1);
    if(!src.data||!match.data) return -1;

    Mat dst = Mat::zeros(src.size(),src.type());

    uchar eToDst[3][256]={0};//Equalized to dst

    uchar tmpMinGrey[3]={255,255,255};
    float tmpConst[3]={0};


    /* match */
    int matchHist[3][256]={0};
    int matchCDF[3][256]={0};
    for(int n=0;n<3;n++){
    //get match histogram
        for(int y=0; y < match.rows; y++) {
            for(int x=0; x < match.cols; x++) {
                uchar grey=*(match.data+y*match.cols*3+x*3+n);
                (matchHist[n][grey])++;
                if(grey<=tmpMinGrey[n]) {
                    tmpMinGrey[n]=grey;
                }
            }
        }

    // get match CDF
        matchCDF[n][0]=matchHist[n][0];
        for(int i=1; i<=255; i++){
            matchCDF[n][i]=matchCDF[n][i-1]+matchHist[n][i];
        }

    // get equalizetion To Dst table of match 
        tmpConst[n]=255.0f/(float)(match.rows*match.cols-matchHist[n][tmpMinGrey[n]]);
        for(int y=0; y < match.rows; y++) {
            for(int x=0; x < match.cols; x++) {
                uchar grey=*(match.data+y*match.cols*3+x*3+n);
                uchar equalized=(uchar)((float)(matchCDF[n][grey]-matchHist[n][tmpMinGrey[n]])*tmpConst[n]);
                eToDst[n][equalized]=grey;
            }
        }
    }
    /* src */
    int srcHist[3][256]={0};
    int dstHist[3][256]={0};
    int srcCDF[3][256]={0};
    tmpMinGrey[0]=255;
    tmpMinGrey[1]=255;
    tmpMinGrey[2]=255;

    for(int n=0; n<3; n++){
        // get src Histogram
        for(int y=0; y < src.rows; y++) {
            for(int x=0; x < src.cols; x++) {
                uchar grey=*(src.data+y*src.cols*3+x*3+n);
                (srcHist[n][grey])++;
                if(grey<=tmpMinGrey[n]) {
                    tmpMinGrey[n]=grey;
                }
            }
        }
        // get src CDF
        srcCDF[n][0]=srcHist[n][0];
        for(int i=1; i<=255; i++){
            srcCDF[n][i]=srcCDF[n][i-1]+srcHist[n][i];
        }
        // get dst
        tmpConst[n]=255.0f/(float)(src.rows*src.cols-srcHist[n][tmpMinGrey[n]]);
        for(int y=0; y < dst.rows; y++) {
            for(int x=0; x < dst.cols; x++) {
                uchar srcgrey=*(src.data+y*src.cols*3+x*3+n);
                uchar equalgrey= (uchar)((float)(srcCDF[n][srcgrey]-srcHist[n][tmpMinGrey[n]])*tmpConst[n]);
                *(dst.data+y*dst.cols*3+x*3+n) = getMatchGrey(equalgrey,eToDst[n]);
            }
        }

        // get dst Histogram
        for(int y=0; y < dst.rows; y++) {
            for(int x=0; x < dst.cols; x++) {
                uchar grey=*(dst.data+y*dst.cols*3+x*3+n);
                (dstHist[n][grey])++;
            }
        }
    }
    // Show src
    namedWindow("Src");
    namedWindow("Match");
    namedWindow("Dst");

    imshow("Src", src);
    imshow("Match",match);
    imshow("Dst", dst);

    myShowHistRGB(srcHist,"SrcHist");
    myShowHistRGB(matchHist,"MatchHist");
    myShowHistRGB(dstHist,"DstHist");

    //save 
    imwrite("dst.jpg",dst);
    // Wait until user press some key
    waitKey();
    return 0;
}