一、 绘制直方图

1.1 代码

#include <iostream>
#include "pch.h"
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace std;
using namespace cv;

/*  Histogram 直方图 */
void getHistogram(Mat& img) {
    // 分割通道
    vector<Mat> rgb_planes;
    split(img, rgb_planes);

    int histSize = 255;

    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 = 600;
    int hist_h = 600;
    int bin_w = cvRound((double)hist_w / histSize);

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

    // 归一化
    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(255, 0, 0), 1, 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), 1, 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(0, 0, 255), 1, 8, 0);
    }

    namedWindow("Histogram", WINDOW_AUTOSIZE);
    imshow("Histogram", histImage);
}

int main()
{
    Mat src;

    src = imread("input.jpg");

    if (!src.data)
        return -1;

    getHistogram(src);

    waitKey();
    return 0;
}

1.2 效果

二、全局直方图均衡

2.1 代码

#include <iostream>
#include "pch.h"
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/types_c.h>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace std;
using namespace cv;

/* Histogram 直方图 */
void getHistogram(Mat& img, String str) {
    Mat dst(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(dst);
    // 分割通道
    vector<Mat> rgb_planes;
    split(dst, rgb_planes);

    int histSize = 255;

    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 = 500;
    int hist_h = 500;
    int bin_w = cvRound((double)hist_w / histSize);

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

    // 归一化
    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(255, 0, 0), 1, 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), 1, 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(0, 0, 255), 1, 8, 0);
    }

    namedWindow(str, WINDOW_AUTOSIZE);
    imshow(str, histImage);
}

/* 直方图均衡 Histogram Equalization*/
Mat equalizeHist_Gray(Mat& img) {
    Mat dst(img.rows, img.cols, CV_8UC1, Scalar(0));
    img.copyTo(dst);
    if (dst.type() == CV_8UC3)
        cvtColor(dst, dst, CV_BGR2GRAY);

    // 计算灰度值
    double grayScale[256] = { 0 };
    for (int i = 0; i < dst.rows; i++)
        for (int j = 0; j < dst.cols; j++)
            grayScale[dst.at<uchar>(i, j)]++;

    // 计算灰度分布密度
    int pixels = dst.rows*dst.cols;
    double density[256] = { 0 };

    for (int i = 0; i < 256; i++)
        density[i] = grayScale[i] / (double)pixels;

    // 计算累计直方图分布
    double temp[256] = { 0 };

    temp[0] = density[0];

    for (int i = 1; i < 256; i++)
        temp[i] = temp[i - 1] + density[i];

    // 累计分布取整
    int result[256] = { 0 };

    for (int i = 0; i < 256; i++)
        result[i] = (int)(255 * temp[i] + 0.5);

    // 灰度值映射
    for (int i = 0; i < dst.rows; i++)
        for (int j = 0; j < dst.cols; j++)
            dst.at<uchar>(i, j) = result[dst.at<uchar>(i, j)];

    return dst;
}

Mat equalizeHist_RGB(Mat& img) {
    Mat temp(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(temp);
    if (temp.type() == CV_8UC1)
        cvtColor(temp, temp, CV_GRAY2RGB);

    Mat matArray[3];
    split(temp, matArray);
    for (int i = 0; i < 3; i++)
        matArray[i] = equalizeHist_Gray(matArray[i]);
    Mat dst;
    merge(matArray, 3, dst);
    return dst;
}

/* PSNR 峰值信噪比 */
double getPSNR(Mat& src, Mat& img) {
    Mat dst;
    
    // 矩阵对应元素作差绝对值
    absdiff(src, img, dst);
    dst.convertTo(dst, CV_32F);

    // 计算矩阵对应元素差的平方
    dst.mul(dst);
    
    Scalar s = sum(dst);
    double sse = s.val[0] + s.val[1] + s.val[2];

    if (sse <= 1e-10)
        return 0;
    else {
        double mse = sse / double(src.channels()*src.total());
        double psnr = 10.0*log10(255 * 255 / mse);
        return psnr;
    }
}

/* SSIM 结构相似性 */
Scalar getSSIM(const Mat& src, const Mat& img) {
    const double k1 = 0.01, k2 = 0.03, L = 255;
    double c1 = (k1*L)*(k1*L), c2 = (k2*L)*(k2*L);

    Mat I1, I2;
    src.convertTo(I1, CV_32F);
    img.convertTo(I2, CV_32F);
    Mat I1_2 = I1.mul(I1);
    Mat I2_2 = I2.mul(I2);
    Mat I1_I2 = I1.mul(I2);

    Mat u1, u2;
    GaussianBlur(I1, u1, Size(11, 11), 1.5);
    GaussianBlur(I2, u2, Size(11, 11), 1.5);

    Mat u1_2 = u1.mul(u1);
    Mat u2_2 = u2.mul(u2);
    Mat u1_u2 = u1.mul(u2);

    Mat sigma1_2, sigma2_2, sigma1_sigma2;
    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
    sigma1_2 -= u1_2;
    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
    sigma2_2 -= u2_2;
    GaussianBlur(I1_I2, sigma1_sigma2, Size(11, 11), 1.5);
    sigma1_sigma2 -= u1_u2;

    Mat t1, t2, t3;
    t1 = 2 * u1_u2 + c1;
    t2 = 2 * sigma1_sigma2 + c2;
    t3 = t1.mul(t2);

    t1 = u1_2 + u2_2 + c1;
    t2 = sigma1_2 + sigma2_2 + c2;
    t1 = t1.mul(t2);

    Mat ssim_map;
    divide(t3, t1, ssim_map);
    Scalar ssim = mean(ssim_map);
    return ssim;
}

/* 调用OpenCV直方图均衡函数 */
Mat equalizeHist_OpenCVGray(Mat& matSrc) {
    Mat matArray;
    cvtColor(matSrc, matSrc, CV_BGR2GRAY);

    Mat dst;
    equalizeHist(matSrc, dst);
    return dst;
}

Mat equalizeHist_OpenCVRGB(Mat& matSrc)
{
    Mat matArray[3];
    split(matSrc, matArray);
    // 直方图均衡化
    for (int i = 0; i < 3; i++)
        equalizeHist(matArray[i], matArray[i]);
    Mat dst;
    merge(matArray, 3, dst);
    return dst;
}

int main()
{
    Mat src, imgGray, imgRGB, imgCheck;

    src = imread("./TestImages/car.jpg");
    if (!src.data)
        return -1;

    imgGray = equalizeHist_Gray(src);
    imgRGB = equalizeHist_RGB(src);
    imgCheck = equalizeHist_OpenCVRGB(src);

    imshow("Raw", src);
    imshow("Equalized_Gray", imgGray);
    imshow("Equalized_RGB", imgRGB);

    cout << "(RAW/Equalized)PSNR:" << getPSNR(src, imgRGB) << "dB" << endl;
    cout << "(RAW/EqualizedCV)PSNR:" << getPSNR(src, imgCheck) << "dB" << endl;
    cout << "(OpenCV/Mine)PSNR:" << getPSNR(imgCheck, imgRGB) << "dB\n" << endl;

    cout << "(RAW/Equalized)SSIM:" << getSSIM(src, imgRGB) << endl;
    cout << "(RAW/EqualizedCV)SSIM:" << getSSIM(src, imgCheck) << endl;
    cout << "(OpenCV/Mine)SSIM:" << getSSIM(imgCheck, imgRGB) << endl;
    
    getHistogram(src, "RawHistogram");
    getHistogram(imgRGB, "EqualizedHistogram");

    waitKey();
    return 0;
}

1.2 效果

Sample1 暗图像均衡

Sample1 直方图对比

Sample2 亮图像均衡

Sample2 直方图对比

Sample3 低对比度均衡

Sample3 直方图对比

Sample4 高对比度均衡

Sample4 直方图对比

原图-灰度图均衡-RGB均衡

直方图对比

峰值信噪比与结构相似性

三、直方图匹配

3.1 代码

#include <iostream>
#include "pch.h"
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/types_c.h>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace std;
using namespace cv;

/* Histogram 直方图 */
void getHistogram(Mat& img, String str) {
    Mat dst(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(dst);
    // 分割通道
    vector<Mat> rgb_planes;
    split(dst, rgb_planes);

    int histSize = 255;

    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 = 500;
    int hist_h = 500;
    int bin_w = cvRound((double)hist_w / histSize);

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

    // 归一化
    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(255, 0, 0), 1, 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), 1, 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(0, 0, 255), 1, 8, 0);
    }

    namedWindow(str, WINDOW_AUTOSIZE);
    imshow(str, histImage);
}

/* PSNR 峰值信噪比 */
double getPSNR(Mat& src, Mat& img) {
    Mat dst;
    
    // 矩阵对应元素作差绝对值
    absdiff(src, img, dst);
    dst.convertTo(dst, CV_32F);

    // 计算矩阵对应元素差的平方
    dst.mul(dst);
    
    Scalar s = sum(dst);
    double sse = s.val[0] + s.val[1] + s.val[2];

    if (sse <= 1e-10)
        return 0;
    else {
        double mse = sse / double(src.channels()*src.total());
        double psnr = 10.0*log10(255 * 255 / mse);
        return psnr;
    }
}

/* SSIM 结构相似性 */
Scalar getSSIM(const Mat& src, const Mat& img) {
    const double k1 = 0.01, k2 = 0.03, L = 255;
    double c1 = (k1*L)*(k1*L), c2 = (k2*L)*(k2*L);

    Mat I1, I2;
    src.convertTo(I1, CV_32F);
    img.convertTo(I2, CV_32F);
    Mat I1_2 = I1.mul(I1);
    Mat I2_2 = I2.mul(I2);
    Mat I1_I2 = I1.mul(I2);

    Mat u1, u2;
    GaussianBlur(I1, u1, Size(11, 11), 1.5);
    GaussianBlur(I2, u2, Size(11, 11), 1.5);

    Mat u1_2 = u1.mul(u1);
    Mat u2_2 = u2.mul(u2);
    Mat u1_u2 = u1.mul(u2);

    Mat sigma1_2, sigma2_2, sigma1_sigma2;
    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
    sigma1_2 -= u1_2;
    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
    sigma2_2 -= u2_2;
    GaussianBlur(I1_I2, sigma1_sigma2, Size(11, 11), 1.5);
    sigma1_sigma2 -= u1_u2;

    Mat t1, t2, t3;
    t1 = 2 * u1_u2 + c1;
    t2 = 2 * sigma1_sigma2 + c2;
    t3 = t1.mul(t2);

    t1 = u1_2 + u2_2 + c1;
    t2 = sigma1_2 + sigma2_2 + c2;
    t1 = t1.mul(t2);

    Mat ssim_map;
    divide(t3, t1, ssim_map);
    Scalar ssim = mean(ssim_map);
    return ssim;
}

/* Histogram matching 直方图匹配 */
Mat matchHistogram_Gray(Mat& img1, Mat& img2) {
    Mat src, dst;
    img1.copyTo(src);
    img2.copyTo(dst);

    double grayScaleSrc[256] = { 0 };
    double grayScaleDst[256] = { 0 };

    // 统计源图像灰度值
    for (int i = 0; i < src.rows; i++) 
        for (int j = 0; j < src.cols; j++) 
            grayScaleSrc[src.at<uchar>(i, j)]++;
        
    // 统计目标图像灰度值
    for (int i = 0; i < dst.rows; i++) 
        for (int j = 0; j < dst.cols; j++) 
            grayScaleDst[dst.at<uchar>(i, j)]++;
    
    // 计算灰度分布密度
    double densitySrc[256] = { 0 }, densityDst[256] = { 0 };
    for (int i = 0; i < 255; i++) {
        densitySrc[i] = grayScaleSrc[i] / (double)(src.rows*src.cols);
        densityDst[i] = grayScaleDst[i] / (double)(dst.rows*dst.cols);
    }
    
    // 计算累计直方图分布
    double tempSrc[256] = { 0 }, tempDst[256] = { 0 };
    tempSrc[0] = densitySrc[0];
    tempDst[0] = densityDst[0];
    for (int i = 1; i < 256; i++) {
        tempSrc[i] = tempSrc[i - 1] + densitySrc[i];
        tempDst[i] = tempDst[i - 1] + densityDst[i];
    }

    int srcMap[256] = { 0 }, dstMap[256] = { 0 };
    for (int i = 0; i < 256; i++) {
        srcMap[i] = (int)(255 * tempSrc[i] + 0.5);
        dstMap[i] = (int)(255 * tempDst[i] + 0.5);
    }

    // 构建直方图匹配灰度映射
    int grayScaleMap[256] = { 0 };
    for (int i = 0; i < 256; i++) {
        // value记录映射后的灰度值
        int value = 0;
        // value_1记录没有对应灰度值时的最接近灰度值
        int value_1 = 0;
        int sum = 0;
        int temp = srcMap[i];
        for (int j = 0; j < 256; j++) {
            if (dstMap[j] == temp) {
                value += j;
                sum++;
            }
            if (dstMap[j] > temp)
            {
                value_1 = j;
                break;
            }
        }
        if (sum == 0) {
            value = value_1;
            sum = 1;
        }
        grayScaleMap[i] = value / sum;
    }

    for (int i = 0; i < src.rows; i++)
        for (int j = 0; j <src.cols; j++) 
            src.at<uchar>(i, j) = grayScaleMap[src.at<uchar>(i, j)];
    return src;
}

Mat matchHistogram_RGB(Mat& img1, Mat& img2) {
    Mat tempSrc(img1.rows, img1.cols, CV_8UC3, Scalar(0));
    Mat tempDst(img2.rows, img2.cols, CV_8UC3, Scalar(0));
    img1.copyTo(tempSrc);
    img2.copyTo(tempDst);
    if (tempSrc.type() == CV_8UC1)
        cvtColor(tempSrc, tempSrc, CV_GRAY2RGB);
    if (tempDst.type() == CV_8UC1)
        cvtColor(tempDst, tempDst, CV_GRAY2RGB);

    Mat matArraySrc[3], matArrayDst[3], result[3];
    split(tempSrc, matArraySrc);
    split(tempDst, matArrayDst);
    
    for (int i = 0; i < 3; i++)
        result[i] = matchHistogram_Gray(matArraySrc[i], matArrayDst[i]);
    Mat dst;
    merge(result, 3, dst);
    return dst;
}

int main()
{
    Mat src, imgGray, imgRGB, imgCheck, dst;

    src = imread("./TestImages/dump.png");
    dst = imread("./TestImages/template_parking.jpg");
    if (!src.data && !dst.data)
        return -1;

    imgRGB = matchHistogram_RGB(src, dst);
    imshow("Raw", src);
    imshow("result", imgRGB);

    getHistogram(src, "RawHistogram");
    getHistogram(dst, "TemplateHistogram");
    getHistogram(imgRGB, "MatchedHistogram");

    cout << "PSNR:" << getPSNR(src, imgRGB) << "dB" << endl;
    cout << "SSIM:" << getSSIM(src, imgRGB) << endl;
    waitKey();
    return 0;
}

3.2 效果

原图-模板图-输出

直方图对比

峰值信噪比与结构相似性

四、局部直方图均衡

4.1 代码

#include <iostream>
#include "pch.h"
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/types_c.h>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace std;
using namespace cv;

/* Local histogram equalization 局部直方图均衡*/
void getCDF(double* s, const double *const c) {
    for (int i = 1; i < 256; i++) {
        s[i] = s[i - 1] + c[i];
    }
}

inline void refresh(double* grayScale, int dim, int value, bool flag) {
    if (flag)
        grayScale[value] += 1.0 / (dim*dim);
    else
        grayScale[value] -= 1.0 / (dim*dim);
}

Mat equalizeHistLocal_Gray(Mat& img) {
    Mat dst(img.rows, img.cols, CV_8UC1, Scalar(0));
    int dim = 5;
    img.copyTo(dst);
    if (dst.type() == CV_8UC3)
        cvtColor(dst, dst, CV_BGR2GRAY);
        
    double grayScale[256] = { 0 };

    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {
            grayScale[dst.at<uchar>(i, j)] += 1.0 / (dim*dim);
        }
    }

    double s[256] = { 0 };
    s[0] = grayScale[0];
    getCDF(s, grayScale);
    int index = dst.at<uchar>(dim / 2, dim / 2);
    dst.at<uchar>(dim / 2, dim / 2) = (int)(s[index] * 255 + 0.5);
    refresh(grayScale, dim, index, false);
    refresh(grayScale, dim, dst.at<uchar>(dim / 2, dim / 2), true);

    for (int i = dim / 2; i < dst.rows - dim / 2; i++) {
        for (int j = dim / 2 + 1; j < dst.cols - dim / 2; j++) {
            bool dirty = true;

            for (int k = 0; k < dim; k++) {
                int old = dst.at<uchar>(i - dim / 2 + k, j - dim / 2 - 1);
                int add = dst.at<uchar>(i - dim / 2 + k, j + dim / 2);
                if (old != add) {
                    refresh(grayScale, dim, old, false);
                    refresh(grayScale, dim, add, true);
                    dirty = false;
                }
            }
            if (!dirty) {
                getCDF(s, grayScale);
                index = dst.at<uchar>(i, j);
                refresh(grayScale, dim, index, false);
                dst.at<uchar>(i, j) = (int)(grayScale[index] * 255 + 0.5);
                refresh(grayScale, dim, dst.at<uchar>(i, j), true);
            }
        }
    }
    return dst;
}

Mat equalizeHistLocal_RGB(Mat& img) {
    Mat temp(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(temp);
    if (temp.type() == CV_8UC1)
        cvtColor(temp, temp, CV_GRAY2RGB);

    Mat matArray[3];
    split(temp, matArray);
    for (int i = 0; i < 3; i++)
        matArray[i] = equalizeHistLocal_Gray(matArray[i]);
    Mat dst;
    merge(matArray, 3, dst);
    return dst;
}

int main()
{
    Mat src, imgGray, imgRGB, imgCheck, dst, t, a;
    src = imread("./TestImages/square.tif");
    dst = equalizeHistLocal_RGB(src);
    imshow("Raw", src);
    imshow("Equalized", dst);
    waitKey();
    return 0;
}

4.2 效果

五、基于直方图统计的局部图像增强

5.1 代码

#include <iostream>
#include "pch.h"
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/types_c.h>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace std;
using namespace cv;

/* Histogram 直方图 */
void getHistogram(Mat& img, String str) {
    Mat dst(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(dst);
    // 分割通道
    vector<Mat> rgb_planes;
    split(dst, rgb_planes);

    int histSize = 255;

    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 = 500;
    int hist_h = 500;
    int bin_w = cvRound((double)hist_w / histSize);

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

    // 归一化
    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(255, 0, 0), 1, 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), 1, 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(0, 0, 255), 1, 8, 0);
    }

    namedWindow(str, WINDOW_AUTOSIZE);
    imshow(str, histImage);
}

/* 直方图均衡 Histogram Equalization*/
Mat equalizeHist_Gray(Mat& img) {
    Mat dst(img.rows, img.cols, CV_8UC1, Scalar(0));
    img.copyTo(dst);
    if (dst.type() == CV_8UC3)
        cvtColor(dst, dst, CV_BGR2GRAY);

    // 计算灰度值
    double grayScale[256] = { 0 };
    for (int i = 0; i < dst.rows; i++)
        for (int j = 0; j < dst.cols; j++)
            grayScale[dst.at<uchar>(i, j)]++;

    // 计算灰度分布密度
    int pixels = dst.rows*dst.cols;
    double density[256] = { 0 };

    for (int i = 0; i < 256; i++)
        density[i] = grayScale[i] / (double)pixels;

    // 计算累计直方图分布
    double temp[256] = { 0 };

    temp[0] = density[0];

    for (int i = 1; i < 256; i++)
        temp[i] = temp[i - 1] + density[i];

    // 累计分布取整
    int result[256] = { 0 };

    for (int i = 0; i < 256; i++)
        result[i] = (int)(255 * temp[i] + 0.5);

    // 灰度值映射
    for (int i = 0; i < dst.rows; i++)
        for (int j = 0; j < dst.cols; j++)
            dst.at<uchar>(i, j) = result[dst.at<uchar>(i, j)];

    return dst;
}

Mat equalizeHist_RGB(Mat& img) {
    Mat temp(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(temp);
    if (temp.type() == CV_8UC1)
        cvtColor(temp, temp, CV_GRAY2RGB);

    Mat matArray[3];
    split(temp, matArray);
    for (int i = 0; i < 3; i++)
        matArray[i] = equalizeHist_Gray(matArray[i]);
    Mat dst;
    merge(matArray, 3, dst);
    return dst;
}

/* PSNR 峰值信噪比 */
double getPSNR(Mat& src, Mat& img) {
    Mat dst;

    // 矩阵对应元素作差绝对值
    absdiff(src, img, dst);
    dst.convertTo(dst, CV_32F);

    // 计算矩阵对应元素差的平方
    dst.mul(dst);

    Scalar s = sum(dst);
    double sse = s.val[0] + s.val[1] + s.val[2];

    if (sse <= 1e-10)
        return 0;
    else {
        double mse = sse / double(src.channels()*src.total());
        double psnr = 10.0*log10(255 * 255 / mse);
        return psnr;
    }
}

/* SSIM 结构相似性 */
Scalar getSSIM(const Mat& src, const Mat& img) {
    const double k1 = 0.01, k2 = 0.03, L = 255;
    double c1 = (k1*L)*(k1*L), c2 = (k2*L)*(k2*L);

    Mat I1, I2;
    src.convertTo(I1, CV_32F);
    img.convertTo(I2, CV_32F);
    Mat I1_2 = I1.mul(I1);
    Mat I2_2 = I2.mul(I2);
    Mat I1_I2 = I1.mul(I2);

    Mat u1, u2;
    GaussianBlur(I1, u1, Size(11, 11), 1.5);
    GaussianBlur(I2, u2, Size(11, 11), 1.5);

    Mat u1_2 = u1.mul(u1);
    Mat u2_2 = u2.mul(u2);
    Mat u1_u2 = u1.mul(u2);

    Mat sigma1_2, sigma2_2, sigma1_sigma2;
    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
    sigma1_2 -= u1_2;
    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
    sigma2_2 -= u2_2;
    GaussianBlur(I1_I2, sigma1_sigma2, Size(11, 11), 1.5);
    sigma1_sigma2 -= u1_u2;

    Mat t1, t2, t3;
    t1 = 2 * u1_u2 + c1;
    t2 = 2 * sigma1_sigma2 + c2;
    t3 = t1.mul(t2);

    t1 = u1_2 + u2_2 + c1;
    t2 = sigma1_2 + sigma2_2 + c2;
    t1 = t1.mul(t2);

    Mat ssim_map;
    divide(t3, t1, ssim_map);
    Scalar ssim = mean(ssim_map);
    return ssim;
}

/* 基于直方图统计的局部增强 */
// 计算概率均值
double getMean(double* pdf) {
    double mean = 0;
    for (int i = 0; i < 256; i++)
        mean += pdf[i] * i;
    return mean;
}

// 计算方差
double getVariance(double* pdf, double mean) {
    double n = 0;
    for (int i = 0; i < 256; i++)
        if (pdf[i] != 0)
            n += (pow((i - mean), 2)*pdf[i]);
    return n;
}

void init(double* pdf) {
    for (int i = 0; i < 256; i++) 
        pdf[i] = 0;
}

Mat statisticsHistogram_Gray(Mat& img, double E, double k0, double k1, double k2) {
    Mat dst(img.rows, img.cols, CV_8UC1, Scalar(0));
    img.copyTo(dst);
    if (dst.type() == CV_8UC3)
        cvtColor(dst, dst, CV_BGR2GRAY);
    int dim = 3;

    double grayScale[256] = { 0 };
    for (int i = 0; i < img.rows; i++)
        for (int j = 0; j < img.cols; j++)
            grayScale[img.at<uchar>(i, j)] += 1.0 / (img.rows*img.cols);

    double mean = getMean(grayScale);

    double SigmaG = sqrt(getVariance(grayScale, mean));

    double Mxy = 0, Sigmaxy = 0;
    for (int i = dim / 2; i < img.rows - dim / 2; i++) {
        for (int j = dim / 2; j < img.cols - dim / 2; j++) {
            init(grayScale);
            for (int m = i - dim / 2; m < i + dim / 2 + 1; m++) {
                for (int n = j - dim / 2; n < j + dim / 2 + 1; n++) {
                    grayScale[img.at<uchar>(m, n)] += 1.0 / (dim*dim);
                }
            }
            
            Mxy = getMean(grayScale);
            Sigmaxy = sqrt(getVariance(grayScale, Mxy));

            if (Mxy <= k0 * mean && k1*SigmaG <= Sigmaxy && Sigmaxy <= k2 * SigmaG) 
                dst.at<uchar>(i, j) = img.at<uchar>(i, j)*E;
        }
    }
    return dst;
}

Mat statisticsHistogram_RGB(Mat& img, double E, double k0, double k1, double k2) {
    Mat temp(img.rows, img.cols, CV_8UC3, Scalar(0));
    img.copyTo(temp);
    if (temp.type() == CV_8UC1)
        cvtColor(temp, temp, CV_GRAY2RGB);

    Mat matArray[3];
    split(temp, matArray);
    for (int i = 0; i < 3; i++)
        matArray[i] = statisticsHistogram_Gray(matArray[i], E, k0, k1, k2);
    Mat dst;
    merge(matArray, 3, dst);
    return dst;
}

int main()
{
    Mat src, imgGray, imgRGB, imgCheck, dst, t, a;
    src = imread("./TestImages/Fig0327(a)(tungsten_original).tif");
    t = equalizeHist_RGB(src);
    dst = statisticsHistogram_RGB(src, 20, 0.4, 0.02, 0.4);
    imshow("Raw", src);
    imshow("Enhanced", dst);
    imshow("Equalized", t);
    cout << "(RAW/Equalized)PSNR:" << getPSNR(src, t) << "dB" << endl;
    cout << "(RAW/Enhanced)PSNR:" << getPSNR(src, dst) << "dB" << endl;
    cout << "(RAW/Equalized)SSIM:" << getSSIM(src, t) << endl;
    cout << "(RAW/Enhanced)SSIM:" << getSSIM(src, dst) << endl;

    waitKey();
    return 0;
}

5.2 效果

原图-全局直方图均衡-基于直方图统计的局部增强

峰值信噪比与结构相似性

基于直方图统计的局部增强

参考链接：
http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/imgproc/histograms/histogram_calculation/histogram_calculation.html
http://www.cnblogs.com/brucemu/archive/2013/10/17/3374558.html
https://blog.csdn.net/u013263891/article/details/82987417