图像处理与模式识别作业一：直方图均衡与灰度拉伸

算法原理

1.直方图均衡

直方图均衡化的基本思想是把原始图像的直方图变换为均匀分布的形式。这样增加了灰度值的动态范围，从而达到增强图像整体对比度的效果。

1）计算图像f(x,y)的各灰度级中像素出现的概率p(i)。

$$p(i)=\frac{n_i}{n},i=0,1,...,L-1 (L：灰度级个数)$$

2) 计算p的累计概率函数c(i),c即为图像的累计归一化直方图

$$c(i)=\sum_{j=0}^{i}p(x_j)$$

3）将c(i)缩放至0~255范围内

$$y(i)=255*c(i)$$

2.灰度拉伸

灰度拉伸属于线性点运算的一种。灰度拉伸。也称对比度拉伸，是一种简单的线性点运算。它扩展图像的直方图，使其充满整个灰度级范围内。
设f(x,y)为输入图像，它的最小灰度级A和最大灰度级B的定义，如下：

$$A = min[f(x,y)]$$ $$B = max[f(x,y)]$$
将A和B分别线性映射到0和255，最终得到的图像g(x,y)为：
$$g(x,y)=(\frac{255}{B-A})[f(x,y)-A]$$

源码解读

1.直方图均衡

首先对全图进行遍历，统计各灰度值出现的频数

 /*扫描全图，统计各灰度级出现的概率*/
void scanimg(Mat &src){
    NUM = src.rows * src.cols;
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar *data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            counts[data[j]] ++;
        }
    }
    for (int i = 0; i < 256; i++){
        change[i] = tween(i);
    }
}

然后计算累积概率，即变换函数

/*计算累积概率*/
double tween(int k){
    int res = 0;
    for (int i = 0; i <= k; i++){
        res += counts[i];
    }
    return (double)res / NUM;
}

为了避免重复运算，我设置了一张转换表change将转换值预存

for (int i = 0; i < 256; i++){
        change[i] = tween(i);
}

然后再次遍历全图，用已经处理好的转换表对所有像素进行转换

/*灰度直方图均衡*/
void myEqualizeHist(Mat &src){
    scanimg(src);
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar * data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            data[j] = 255 * change[data[j]];
        }
    }
}

2.灰度拉伸

首先第一步直方图均衡一样，也是统计各灰度出现的频数

/*扫描全图，统计各灰度级出现的概率*/
void scanimg(Mat &src){
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar *data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            counts[data[j]] ++;
        }
    }

然后统计最低灰度值与最高灰度值，由于高光部分的灰度值并不为0但数值较小，为了让效果更明显，我认为限定了阈值为105

for (int i = 0; i < 256; i++){
        if (counts[i] > 105){
            a = i;
            break;
        }
    }
    for (int i = 0; i < 256;i++)
        if (counts[255 - i] > 105){
            b = 255 - i;
            break;
        }

由于灰度拉伸的结果是全等级，因而将原计算公式化简为255/(a-b) *(f(x,y) -a)
得到转换函数：

double  tween(int k){
    double t = 255.0 / (b - a)*(k - a);
    return t;
}

最后遍历全图对每个像素进行转换：

void myProc(Mat &src){
    scanimg(src);
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar * data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            data[j] = tween(data[j]);
        }
    }
}

3.RGB三通道拉伸

计算方法基本与灰度拉伸一样，只不过需要对三通道的像素值分别进行一次计算，统计时也需要分开单独统计，不过由于三个通道的像素值分布是不同的，因此阈值也分别根据统计情况设置为300、105和200

结果分析

1.直方图均衡

2.灰度拉伸

3.RGB拉伸

个人讨论

通过本次大作业掌握了最最基础的图像处理方法，过程中遇到了不少小问题，比如转换函数的计算要用浮点不能用整型，opencv中RGB三通道的排列是反的等等。
从结果来看，直方图均衡和灰度拉伸的效果差异非常大，灰度拉伸显得更加自然，而直方图均衡则有些失真。因此这也解释了为什么它们的应用场景也有所不同，直方图均衡多用于图像增强，显示出一些不易发现的细节，而灰度拉伸则类似于PS中的亮度提升。

完整代码

/*
灰度直方图均衡
author: Billow_Tau
编译环境：Visual Studio 2013
第三方库：opencv
*/

#include <iostream>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>

using namespace std;
using namespace cv;

int counts[256] = { 0 };
double change[256] = { 0 };
int NUM;
int a = 0, b = 128;

void myShowHist(Mat src, MatND hist){
    int bins = 256;
    int hist_size[] = { bins };
    float range[] = { 0, 256 };
    const float* ranges[] = { range };

    int planes[] = { 0 };

    calcHist(&src, 1, planes, Mat(), hist, 1, hist_size, ranges, true, false);

    int hist_h = 256;
    double max_val;
    int scale = 2;
    minMaxLoc(hist, 0, &max_val, 0, 0);

    Mat hist_img = Mat::zeros(hist_h, bins * scale, CV_8UC3);

    for (int i = 0; i < bins; i++){
        float bin_val = hist.at<float>(i);
        int intensity = cvRound(bin_val* hist_h/ max_val);
        rectangle(hist_img, Point(i*scale, hist_h-1), Point((i+1)*scale -1, hist_h - intensity), CV_RGB(255,255,255));
    }

    imshow("GRAY", hist_img);
    //imwrite("hist_src.png",hist_img);
}

/*计算累积概率*/
double tween(int k){
    int res = 0;
    for (int i = 0; i <= k; i++){
        res += counts[i];
    }
    return (double)res / NUM;
}

/*扫描全图，统计各灰度级出现的概率*/
void scanimg(Mat &src){
    NUM = src.rows * src.cols;
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar *data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            counts[data[j]] ++;
        }
    }
    for (int i = 0; i < 256; i++){
        change[i] = tween(i);
    }
}

/*灰度直方图均衡*/
void myEqualizeHist(Mat &src){
    scanimg(src);
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar * data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            data[j] = 255 * change[data[j]];
        }
    }
}

int main(){

    Mat img;
    MatND hist2;
    img = imread("D:\\opencv\\project\\test1.jpg", 1);

    if (img.empty()){
        return -1;
    }

    cvtColor(img, img, CV_BGR2GRAY);

    //myEqualizeHist(img);

    myShowHist(img, hist2);

    //imwrite("src.png",img);

    namedWindow("MyWindow", CV_WINDOW_NORMAL);
    imshow("MyWindow", img);
    waitKey(0);

    return 0;
}

/*
灰度拉伸
author: Billow_Tau
编译环境：Visual Studio 2013
第三方库：opencv
*/

#include <iostream>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>

using namespace std;
using namespace cv;

int counts[256] = { 0 };
int a = 0, b = 255;

void myShowHist(Mat src, MatND hist){
    int bins = 256;
    int hist_size[] = { bins };
    float range[] = { 0, 256 };
    const float* ranges[] = { range };

    int planes[] = { 0 };

    calcHist(&src, 1, planes, Mat(), hist, 1, hist_size, ranges, true, false);

    int hist_h = 256;
    double max_val;
    int scale = 2;
    minMaxLoc(hist, 0, &max_val, 0, 0);

    Mat hist_img = Mat::zeros(hist_h, bins * scale, CV_8UC3);

    for (int i = 0; i < bins; i++){
        float bin_val = hist.at<float>(i);
        int intensity = cvRound(bin_val* hist_h / max_val);
        rectangle(hist_img, Point(i*scale, hist_h - 1), Point((i + 1)*scale - 1, hist_h - intensity), CV_RGB(255, 255, 255));
    }

    imshow("GRAY", hist_img);
    //imwrite("hist_dst2.png",hist_img);
}

double  tween(int k){
    double t = 255.0 / (b - a)*(k - a);
    return t;
}

/*扫描全图，统计各灰度级出现的概率*/
void scanimg(Mat &src){
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar *data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            counts[data[j]] ++;
        }
    }
    /*for (int i = 0; i<256; i++)
        cout << counts[i] << ' ';
    cout << '\n' << endl;*/
    for (int i = 0; i < 256; i++){
        if (counts[i] > 105){
            a = i;
            break;
        }
    }
    for (int i = 0; i < 256;i++)
        if (counts[255 - i] > 105){
            b = 255 - i;
            break;
        }
}

void myProc(Mat &src){
    scanimg(src);
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar * data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j++){
            data[j] = tween(data[j]);
        }
    }
}

int main(){

    Mat img;
    MatND hist2;
    img = imread("D:\\opencv\\project\\test1.jpg", 1);

    if (img.empty()){
        return -1;
    }

    cvtColor(img, img, CV_BGR2GRAY);

    myProc(img);

    myShowHist(img, hist2);

    //imwrite("dst2.png",img);

    namedWindow("MyWindow", CV_WINDOW_NORMAL);
    imshow("MyWindow", img);
    waitKey(0);

    return 0;
}

/*
RGB拉伸
author: Billow_Tau
编译环境：Visual Studio 2013
第三方库：opencv
*/

#include <iostream>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>

using namespace std;
using namespace cv;

int Rcounts[256] = { 0 };
int Gcounts[256] = { 0 };
int Bcounts[256] = { 0 };
int Ra = 0, Rb = 255;
int Ga = 0, Gb = 255;
int Ba = 0, Bb = 255;

void myShowHist(Mat src){
    /// 分割成3个单通道图像 ( R, G 和 B )
    vector<Mat> rgb_planes;
    split(src, rgb_planes);

    /// 设定bin数目
    int histSize = 255;

    /// 设定取值范围 ( R,G,B) )
    float range[] = { 0, 255 };
    const float* histRange = { range };

    bool uniform = true; bool accumulate = false;

    Mat r_hist, g_hist, b_hist;

    /// 计算直方图:
    calcHist(&rgb_planes[0], 1, 0, Mat(), r_hist, 1, &histSize, &histRange, uniform, accumulate);
    calcHist(&rgb_planes[1], 1, 0, Mat(), g_hist, 1, &histSize, &histRange, uniform, accumulate);
    calcHist(&rgb_planes[2], 1, 0, Mat(), b_hist, 1, &histSize, &histRange, uniform, accumulate);

    // 创建直方图画布
    int hist_w = 400; int hist_h = 400;
    int bin_w = cvRound((double)hist_w / histSize);

    Mat histImage(hist_w, hist_h, CV_8UC3, Scalar(0, 0, 0));

    /// 将直方图归一化到范围 [ 0, histImage.rows ]
    normalize(r_hist, r_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
    normalize(g_hist, g_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
    normalize(b_hist, b_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());

    /// 在直方图画布上画出直方图
    for (int i = 1; i < histSize; i++)
    {
        line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(r_hist.at<float>(i - 1))),
            Point(bin_w*(i), hist_h - cvRound(r_hist.at<float>(i))),
            Scalar(0, 0, 255), 2, 8, 0);
        line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(g_hist.at<float>(i - 1))),
            Point(bin_w*(i), hist_h - cvRound(g_hist.at<float>(i))),
            Scalar(0, 255, 0), 2, 8, 0);
        line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(b_hist.at<float>(i - 1))),
            Point(bin_w*(i), hist_h - cvRound(b_hist.at<float>(i))),
            Scalar(255, 0, 0), 2, 8, 0);
    }

    /// 显示直方图
    namedWindow("calcHist Demo", CV_WINDOW_AUTOSIZE);
    imshow("calcHist Demo", histImage);
    //imwrite("hist_rgb_dst.png", histImage);
}

double  tween(int k, int select){
    double t;
    switch (select){
    case 2:
        t = 255.0 / (Rb - Ra)*(k - Ra);
        break;
    case  1:
        t = 255.0 / (Gb - Ga)*(k - Ga);
        break;
    case 0:
        t = 255.0 / (Bb - Ba)*(k - Ba);
        break;
    }
    return t;
}

/*扫描全图，统计各通道出现的概率*/
void scanimg(Mat &src){
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar *data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j+=3){
            Rcounts[data[j + 2]] ++;
            Gcounts[data[j + 1]]++;
            Bcounts[data[j]]++;
        }
    }
    for (int i = 0; i<256; i++)
        cout << Rcounts[i] << ' ';
    cout << '\n' << endl;
    for (int i = 0; i<256; i++)
        cout << Gcounts[i] << ' ';
    cout << '\n' << endl;
    for (int i = 0; i<256; i++)
    cout << Bcounts[i] << ' ';
    cout << '\n' << endl;
    for (int i = 0; i < 256; i++)
        if (Rcounts[i] > 300){Ra = i;break;}
    for (int i = 0; i < 256; i++)
        if (Gcounts[i] > 105){ Ga = i; break; }
    for (int i = 0; i < 256; i++)
        if (Bcounts[i] > 200){ Ba = i; break; }
    for (int i = 0; i < 256; i++)
        if (Rcounts[255 - i] > 300){Rb = 255 - i;break;}
    for (int i = 0; i < 256; i++)
        if (Gcounts[255 - i] > 105){ Gb = 255 - i; break; }
    for (int i = 0; i < 256; i++)
        if (Bcounts[255 - i] > 200){ Bb = 255 - i; break; }
}

void myProc(Mat &src){
    scanimg(src);
    int nr = src.rows;
    int nc = src.cols * src.channels();
    for (int i = 0; i < nr; i++){
        uchar * data = src.ptr<uchar>(i);
        for (int j = 0; j < nc; j+=3){
            data[j+2] = tween(data[j+2],2);
            data[j + 1] = tween(data[j + 1], 1);
            data[j] = tween(data[j], 0);
        }
    }
}

int main(){

    Mat img;
    img = imread("D:\\opencv\\project\\test1.jpg", 1);

    if (img.empty()){
        return -1;
    }

    //cvtColor(img, img, CV_BGR2GRAY);

    myProc(img);

    myShowHist(img);

    //imwrite("dst_rgb.png",img);

    namedWindow("MyWindow", CV_WINDOW_NORMAL);
    imshow("MyWindow", img);
    waitKey(0);

    return 0;
}