直方图均衡化的 C++ 实现（基于 openCV）

这是数字图像处理课的大作业，完成于 2013/06/17，需要调用 openCV 库，完整源码和报告如下：

#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string>

/* 灰度级结点 */
typedef struct {
    int pixels;        // 灰度级对应像素个数
    float rate;        // 像素比例
    float accuRate;    // 累计像素比例
    int map;        // 到均衡化后的灰度级的映射
} levNode;

void histeqGray(IplImage* pGray, int levels, int argc);
IplImage* histImage(IplImage* pSrc, int histWidth, int histHeight, int nScale);

int main(int argc, char* argv[])
{
    int levels;
    std::string imgName, inTmp;
    if (argc == 3) {
        levels = atoi(argv[1]);
        imgName = argv[2];
    } 
    else if (argc == 2)
        imgName = argv[1];
    else {
        printf("usage: histeq [levels] image_name 
");
        return -1;
    }
    
    IplImage* pSrc = cvLoadImage(imgName.c_str(), CV_LOAD_IMAGE_UNCHANGED);
    int channel = pSrc->nChannels;

    IplImage* pChnl[4] = { NULL };

    for (int i = 0; i < channel; ++i)
        pChnl[i] = cvCreateImage(cvGetSize(pSrc), pSrc->depth, 1);

    cvSplit(pSrc, pChnl[0], pChnl[1], pChnl[2], pChnl[3]);
    
    for (int i = 0; i < channel; ++i)
        histeqGray(pChnl[i], levels, argc);

    IplImage* pEql = cvCreateImage(cvGetSize(pSrc), pChnl[0]->depth, pSrc->nChannels);

    cvMerge(pChnl[0], pChnl[1], pChnl[2], pChnl[3], pEql);

    inTmp = imgName + "_Eql.jpg";
    cvSaveImage(inTmp.c_str(), pEql);

    //cvNamedWindow(imgName.c_str(), CV_WINDOW_AUTOSIZE);
    cvShowImage(imgName.c_str(), pSrc);
    //cvNamedWindow(inTmp.c_str(), CV_WINDOW_AUTOSIZE);
    cvShowImage(inTmp.c_str(), pEql);

    IplImage* pSrcGray = cvCreateImage(cvGetSize(pSrc), IPL_DEPTH_8U, 1);
    if (pSrc->nChannels == 3)
        cvCvtColor(pSrc, pSrcGray, CV_BGR2GRAY);
    else
        cvCopyImage(pSrc, pSrcGray);
    IplImage* pEqlGray = cvCreateImage(cvGetSize(pEql), IPL_DEPTH_8U, 1);
    if (pSrc->nChannels == 3)
        cvCvtColor(pEql, pEqlGray, CV_BGR2GRAY);
    else
        cvCopyImage(pEql, pEqlGray);
    imgName += "_Hist.jpg";
    inTmp += "_Hist.jpg";
    int nScale = 2;
    int histWidth = /*pSrc->width * nScale*/256 * nScale;
    int histHeight = /*pSrc->height*/128;
    IplImage* pSrcGrayHist = histImage(pSrcGray, histWidth, histHeight, nScale);
    IplImage* pEqlGrayHist = histImage(pEqlGray, histWidth, histHeight, nScale);
    cvSaveImage(imgName.c_str(), pSrcGrayHist);
    cvSaveImage(inTmp.c_str(), pEqlGrayHist);
    cvShowImage(imgName.c_str(), pSrcGrayHist);
    cvShowImage(inTmp.c_str(), pEqlGrayHist);

    cvWaitKey();

    cvReleaseImage(&pEql);
    cvReleaseImage(&pEqlGray);
    for (int i = 0; i < channel; ++i)
        cvReleaseImage(&pChnl[i]);
    cvReleaseImage(&pSrc);
    cvReleaseImage(&pSrcGray);

    return 0;
}

/*
* 直方图均衡化函数 
* pGray为输入的灰度图
* levels为均衡化的灰度级
*/
void histeqGray(IplImage* pGray, int levels, int argc)
{
    int depth = pGray->depth;
    printf("%d 
", depth);
    int width = pGray->width;
    int height = pGray->height;
    int sumPixels = width * height;        // 总像素数
    printf("%d 
", sumPixels);
    int values = static_cast<int>(pow((float)2, depth)); // 根据图像深度计算像素取值范围
    if (argc == 2) levels = values;
    printf("%d 
", levels);

    int outDepth;
    /*if (levels <= 2)
        outDepth = 1;
    else*/ if (levels <= 256)
        outDepth = 8;
    else if (levels <= 65536)
        outDepth = 16;

    assert(levels <= values);
    int intervals = values / levels; // 根据像素取值范围和灰度级求每个灰度级的像素间隔
    levNode* levNodes = (levNode*)calloc(levels, sizeof(levNode)); // 生成灰度结点
    //for (int lev = 0; lev < levels; ++lev) printf("%d 
", levNodes[lev].pixels);
    //char* pValues = pGray->imageData;

    /* 统计每个灰度级的像素个数 */
    for (int y = 0; y < height; ++y)
        for (int x = 0; x < width; ++x) {
            CvScalar scal = cvGet2D(pGray, y, x);
            int val = (int)scal.val[0];
            //printf("%d 
", val);
            for (int lev = 0; lev < levels; ++lev) {
                if ( val >= intervals*lev && val < intervals*(lev+1)) {
                    ++levNodes[lev].pixels; break;
                }
            }
        }

    int sum = 0;
    for (int lev = 0; lev < levels; ++lev)
        sum += levNodes[lev].pixels;
    printf("%d 
", sum);

    /* 计算每个灰度级像素比例和累计比例 */
    levNodes[0].accuRate = levNodes[0].rate = levNodes[0].pixels / (float)sumPixels;
    levNodes[0].map = (int)(levNodes[0].accuRate * (levels - 1) + 0.5);
    printf("%d 
", levNodes[0].pixels);
    for (int lev = 1; lev < levels; ++lev) {
        levNodes[lev].rate = levNodes[lev].pixels / (float)sumPixels;
        levNodes[lev].accuRate = levNodes[lev-1].accuRate + levNodes[lev].rate;
        levNodes[lev].map = (int)(levNodes[lev].accuRate * (levels - 1) + 0.5);
    }
    printf("%f 
", levNodes[levels-1].accuRate);

    /* 生成均衡化后的图像 */
    for (int y = 0; y < height; ++y)
        for (int x = 0; x < width; ++x) {
            CvScalar scal = cvGet2D(pGray, y, x);
            int val = (int)scal.val[0];
            //printf("%d 
", val);
            for (int lev = 0; lev < levels; ++lev) {
                if (val >= intervals*lev && val < intervals*(lev+1)) {                    
                        scal.val[0] = levNodes[lev].map;
                        //printf("%f 
", scal.val[0]);
                        cvSet2D(pGray, y, x, scal);
                        break;
                }
            }
        }
    pGray->depth = outDepth;

    free(levNodes);
}

/*
* 绘制直方图函数
*/
IplImage* histImage(IplImage* pSrc, int histWidth, int histHeight, int nScale)
{
    int histSize = static_cast<int>(pow((float)2, pSrc->depth));
    CvHistogram* pHist = cvCreateHist(/*pSrc->nChannels*/1, &histSize, CV_HIST_ARRAY);
    cvCalcHist(&pSrc, pHist);

    IplImage* pHistImg = cvCreateImage(cvSize(histWidth, histHeight), IPL_DEPTH_8U, 1);
    cvRectangle(pHistImg, cvPoint(0,0), cvPoint(pHistImg->width,pHistImg->height), CV_RGB(255,255,255), CV_FILLED);

    float histMaxVal = 0;
    cvGetMinMaxHistValue(pHist, 0, &histMaxVal);

    for(int i = 0; i < histSize; i++)
    {
        float histValue= cvQueryHistValue_1D(pHist, i); // 像素为i的直方块大小
        int nRealHeight = cvRound((histValue / histMaxVal) * histHeight);  // 要绘制的高度
        cvRectangle(pHistImg,  
            cvPoint(i*nScale, histHeight - 1),  
            cvPoint((i + 1)*nScale - 1, histHeight - nRealHeight),
            cvScalar(i),   
            CV_FILLED
        );   
    }
    //cvFillConvexPoly

    cvReleaseHist(&pHist);
    return pHistImg;
}

一、直方图均衡化概述

直方图均衡化是一种图像增强方法，其基本思想是把给定图像的直方图分布改造成均匀分布的直方图，从而增加象素灰度值的动态范围，达到增强图像整体对比度的效果。由信息学的理论来解释，具有最大熵（信息量）的图像为均衡化图像。

直方图均衡化可表示为：，t为某个象素变换后的灰度级，s为该象素变换前的灰度级。

该灰度变换函数应满足如下两个条件：

1）f(s)在范围内是单值递增函数；

2）对有

条件1：保证原图各灰度级在变换后仍保持从黑到白（或从白到黑）的排列顺序；

条件2：保证变换前后灰度值动态范围的一致性。

可以证明累积分布函数（cumulative distribution function CDF）满足上述两个条件并能将s的分布转换为t的均匀分布。

事实上，s的CDF就是原始图的累积直方图，即：

其中、、k=0,1,…,L-1

根据这个公式，可以直接算出直方图均衡化后各象素的灰度值。需要取整，以满足数字图象的要求。

 

二、算法步骤

步骤	运算
1	列出原始图灰度级
2	统计原始直方图各灰度级象素数
3	计算原始直方图（像素比例）
4	计算累积直方图 （累计像素比例）
5	取整
6	确定映射对应关系（ → ）
7	计算新的直方图

三、算法测试

1、灰度图

2、彩色图

 

四、结果分析

（1）对于灰度图和彩色图，算法结果都不错，直方图显示像素分布很广、很平均。

（2）直方图均衡化的优点：自动增强整个图像的对比度。

（3）直方图均衡化的不足：具体增强效果不易控制，处理的结果总是得到全局均衡化的直方图。