引言



在大部分的图像处理程序中，其中必不可少的一步就是对传入的彩图进行灰度处理，将三个通道的RGB图片转化为单通道的Gray图，而对于灰度图进行直方图统计同样是观察检测图像特征的常用方法。在OpenCV中已经有成熟的封装函数进行上述功能的实现，本文主要讲述CUDA实现，加快对大图像的处理速度。

任务要求



输入一张彩色图片，通过CUDA将其转换为灰度图，并对灰度图进行灰度直方图统计。

实现思路



关于彩色图片灰度化原理可参照此博客图像灰度化总结简单来说就是将传入的RGB图中三个通道的像素值分别乘以一定系数即可得到对应像素点灰度值。而在灰度值统计时，由于在CUDA中block内的线程同时并行，故无法采用简单的累加统计，需用到CUDA内自带的原子操作函数atomicAdd()。

实现代码

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cuda.h>
#include <device_functions.h>
#include <opencv2\opencv.hpp>
#include <iostream>
#include <math.h>
using namespace std;
using namespace cv;

//输入图像为BGR图，将其转化为gray图
__global__ void rgb2grayInCuda(uchar3 *dataIn, unsigned char *dataOut, int imgHeight, int imgWidth)
{
    int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
    int yIndex = threadIdx.y + blockIdx.y * blockDim.y;

    if (xIndex < imgWidth && yIndex < imgHeight)
    {
        uchar3 rgb = dataIn[yIndex * imgWidth + xIndex];

        dataOut[yIndex * imgWidth + xIndex] = 0.299f * rgb.x + 0.587f * rgb.y + 0.114f * rgb.z;
    }
}

//灰度直方图统计
__global__ void imHistInCuda(unsigned char *dataIn, int *hist)
{
    int threadIndex = threadIdx.x + threadIdx.y * blockDim.x;
    int blockIndex = blockIdx.x + blockIdx.y * gridDim.x;
    int index = threadIndex + blockIndex * blockDim.x * blockDim.y;

    atomicAdd(&hist[dataIn[index]], 1);

}

int main()
{
    //传入图片
    Mat srcImg = imread("out1.bmp");

    int imgHeight = srcImg.rows;
    int imgWidth = srcImg.cols;
    Mat grayImg(imgHeight, imgWidth, CV_8UC1, Scalar(0));//输出灰度图
    int hist[256];//灰度直方图统计数组
    memset(hist, 0, 256 * sizeof(int));

    //在GPU中开辟输入输出空间
    uchar3 *d_in;
    unsigned char *d_out;
    int *d_hist;

    cudaMalloc((void**)&d_in, imgHeight * imgWidth * sizeof(uchar3));
    cudaMalloc((void**)&d_out, imgHeight * imgWidth * sizeof(unsigned char));
    cudaMalloc((void**)&d_hist, 256 * sizeof(int));

    //将图像数据传入GPU中
    cudaMemcpy(d_in, srcImg.data, imgHeight * imgWidth * sizeof(uchar3), cudaMemcpyHostToDevice);
    cudaMemcpy(d_hist, hist, 256 * sizeof(int), cudaMemcpyHostToDevice);

    dim3 threadsPerBlock(32, 32);
    dim3 blocksPerGrid((imgWidth + threadsPerBlock.x - 1) / threadsPerBlock.x, (imgHeight + threadsPerBlock.y - 1) / threadsPerBlock.y);

    //灰度化
    rgb2grayInCuda << <blocksPerGrid, threadsPerBlock >> >(d_in, d_out, imgHeight, imgWidth);

    //灰度直方图统计
    imHistInCuda << <blocksPerGrid, threadsPerBlock >> >(d_out, d_hist);

    //将数据从GPU传回CPU
    cudaMemcpy(hist, d_hist, 256 * sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(grayImg.data, d_out, imgHeight * imgWidth * sizeof(unsigned char), cudaMemcpyDeviceToHost);


    cudaFree(d_in);
    cudaFree(d_out);
    cudaFree(d_hist);

    return 0;
}



实现结果



输入图片 



输出图片 



通过比对发现CUDA输出结果与OpenCV输出结果一致~

 

目录
1.Objectives:
2.Experiment Content：
3.Experiment Principle:
4.Experiment Steps Result and Conlusion：
5.Appendix（programs and images）：

1.Objectives:
1、 熟悉灰度直方图的概念及计算方法；

2、 熟悉直方图均衡化的计算过程；

3、 计算并绘制图像直方图，实现直方图均衡化。
2.Experiment Content：
学习使用函数 imhist, histeq, bar, stem, plot, imadjust，及title, axis, set等描述图像工具。
3.Experiment Principle:
1.直方图：直方图定义就是一个离散函数 h（rk）=nkrk是区间内的第k级灰度，nk为图像中出现这种灰度级的像素数，可以用bar stem plot三种方式显示出来‚
2.直方图均衡：其实均衡就是一个累计分布函数（CDF）



其中Sk输出图像的灰度值 rk是输入的灰度值，nj是在某个灰度值的像素点的个数，n是整个像素点的个数。
下面我用一个例子说明如何进行直方图的均衡的
假设有9个像素点的一幅图



我们从上面的表格中可以发现 原来的20 30 40 的灰度值相隔很近，不容易用肉眼辨别 均衡后 变成了85 170 255 拉大了灰度值的差别 从而增加了对比度
4.Experiment Steps Result and Conlusion：
（1）直方图显示



图1 直方图表示

Con：从图中看出，整体的图片相对比较暗，因此灰度主要集中在100值的一片
（2）分别用bar 和stem 函数显示直方图
用plot 函数显示直方图 由 h=imhist 获得直方图；分别用bar 和stem 显示直方图h，并通过参数调整，改变直方图的显示方式。用axis 设置轴的最大、最小值（例如：axis([0 255 0 15000]);），用set 设置显示坐标的间隔（例如： set(gca,’xtick’, 0:50:255）)



图2 用bar，stem，plot显示直方图

Con：这里的直方图统计分成了25块，因此不像上次图那么密，可以明显看到100值处的像素点特别多，因为主要都集中在那里，并且有把集中在那里的大概10个灰度值的像素点都加在了一起，因此数据量特别大。坐标的设置都是用了一些matlab中的函数工具。
（3）用imadjust 函数调整图像对比度，并用imhist 查看调整前后直方图的变化。



图3用imadjust 函数调整图像对比度
两张图都是对图片的灰度进行了处理 设置的变化范围都是0-1

对于第一张图 gamma<1，图像变亮 后面的直方图也汇聚在高亮度值150附近

对于第一张图 gamma>1，图像变暗 后面的直方图也汇聚在高亮度值45附近

（4）用histeq 进行直方图均衡化，并用imhist 查看均衡化前后直方图的变化。



图4 用histeq 进行直方图均衡化
Con：均衡的目的是提高的对比度，亮的更亮，暗的更暗，从结果来看，爱恩斯坦背后的书更清晰了，说明对比度的提高能帮我们“看清”东西。但对人的处理还不如之前，说明还是要根据实际情况来进行相应的处理的。另外看直方图这边，原来低亮度 高亮度的地方都没有像素点 均衡化之后，几乎布满了整个灰度区域，说明对原来图像的亮度范围进行了拉伸。像素点的亮度集中在高亮度区域。
（5）用模糊算法进行亮度均衡，利用相应函数获得以下结果



图5 用模糊算法进行亮度均衡
Con：使用addpath添加相应的.p文件，便于之后使用fuzzysysfcn函数，明显可以看出模糊处理后的图像效果很好，更加清晰真实。
5.Appendix（programs and images）：
代码：
clear all;clc;
f=imread('E:\数字图像处理\程序与图像\图像库\EINSTEIN.tif');
figure(1);
subplot(1,2,1);
imshow(f);
title('EINSTEIN');
subplot(1,2,2);
imhist(f);%画出直方图
axis([0 255 0 40000]);%为直方图设置坐标轴  

figure(2);
subplot(2,2,1);
h=imhist(f,25);%这个25参量的意思是共有25个容器 就是把255个灰度值分成25个来进行统计
horz=linspace(0,255,25);%生成一个范围是0到255的向量，向量中有25个点，分别均匀排列
bar(horz,h);%y用bar方式显示直方图，h是直方图的信息 
axis([0 255 0 120000]);%设置xy轴标记的刻度 四个参数分别是x轴的起始刻度，终止刻度 轴的起始刻度，终止刻度
set(gca,'xtick',0:50:255) %gca是获得当前轴，然后在x轴上 每隔50个值标记刻度
set(gca,'ytick',0:20000:120000);%同上，不过是在y轴
title('bar 默认宽度');
subplot(2,2,2);
bar(horz,h,0.5);%y用bar方式显示直方图，h是直方图的信息 0.5这个参数是可以调整bar的宽度的信息 
axis([0 255 0 120000]);
set(gca,'xtick',0:50:255);
set(gca,'ytick',0:20000:120000);%同上
title('bar 0.5的宽度');
subplot(2,2,3);
stem(horz,h);%用杆状图绘制直方图
axis([0 255 0 120000]);
set(gca,'xtick',0:50:255)
set(gca,'ytick',0:20000:120000);%同上
title('stem');
subplot(2,2,4);plot(horz,h);%用plot绘制曲线图
axis([0 255 0 120000]);
set(gca,'xtick',0:50:255) 
set(gca,'ytick',0:20000:120000);%同上
title('plot'); 

figure(3);
subplot(2,2,1);
f1=imadjust(f,[0 1],[0 1],0.5);%f是输入的图像 第二、三个矩阵分别是变化前 变化后的范围 最后的值gamma是决定变化的函数是上凸还是下凹%大于1下凹，图像变暗，大于1上凸，图像变亮
imshow(f1);
title('gamma>1，图像变暗');%映射图像下凹，相应灰度级降低
subplot(2,2,2);
imhist(f1);
subplot(2,2,3);
f2=imadjust(f,[0 1],[0 1],2);
imshow(f2);
title('gamma<1,图像变亮');
subplot(2,2,4);
imhist(f2);

figure(4);
subplot(2,2,1);
imshow(f);
subplot(2,2,2);
imhist(f);
title('原图直方图');
subplot(2,2,3);
g=histeq(f,256);%求出原图直方图均衡后的图像中每个像素点的灰度数据
imshow(g);%显示原图直方图均衡后的图像
subplot(2,2,4);
imhist(g);%绘制出原图直方图均衡后的直方图
title('直方图均衡');  

addpath('E:\数字图像处理\程序与图像\dipum_toolbox_2.0.2(only run)');%添加相应的.p文件，便于之后使用fuzzysysfcn函数 

udark=@(z)1-sigmamf(z, 0.35, 0.5);%进行模糊处理的隶属度函数
ugray=@(z)triangmf(z, 0.35, 0.5, 0.65);
ubright=@(z)sigmamf(z, 0.5, 0.65);
udarker=@(z)bellmf(z, 0.0, 0.1);
umidgray=@(z)bellmf(z, 0.4, 0.5);
ubrighter=@(z)bellmf(z, 0.8, 0.9); 
outmf={udark; ugray; ubright};%表示输出模糊函数
rules={udarker; umidgray; ubrighter}; %表示输入
F=fuzzysysfcn(outmf, rules, [0,1]); %使用
fuzzysysfcn函数及逆行模糊处理
z=linspace(0, 1, 256);%产生0-1范围内的指定数量点数，等分为256份，并返回一个行向量。
T=F(z);
I=f;
g=intrans(I, 'specified', T);
figure(5), imshow(g);%模糊处理的结果
figure(6), imshow(I);%原图
figure(7),plot(T);%灰度分布
结果：

图像灰度直方图
灰度直方图是图像灰度级的函数，用来描述每个灰度级在图像矩阵中的像素个数或者占有率。直方图显示图像数据时会以左暗右亮的分布曲线形式呈现出来。横坐标是灰度级，纵坐标是该灰度级出现的频率。图像的对比度是通过灰度级范围来度量的，而灰度级范围可通过观察灰度直方图得到，灰度级范围越大代表对比度越高；反之对比度越低，低对比度的图像在视觉上给人的感觉是看起来不够清晰，所以通过算法调整图像的灰度值，从而调整图像的对比度是有必要的。
获得图像的灰度直方图def calGrayHist(img):
    h,w = img.shape[:2]
gray_hist = np.zeros([256], np.uint64)
for i in range(h):
for j in range(w):
            gray_hist[img[i][j]] += 1
return gray_hist
不过通常会将纵坐标归一化到[0,1][0,1]区间内，也就是将灰度级出现的频率（像素个数）除以图像中像素的总数，得到归一化灰度直方图，计算公式如下：
其中，
是像素的灰度级，
是具有灰度
的像素的个数，MN是图像中总的像素个数。
直方图均衡化 Histogram Equalization
假如图像的灰度分布不均匀，其灰度分布集中在较窄的范围内，使图像的细节不够清晰，对比度较低。通常采用直方图均衡化及直方图规定化两种变换，使图像的灰度范围拉开或使灰度均匀分布，从而增大反差，使图像细节清晰，以达到增强的目的。
直方图均衡化，对图像进行非线性拉伸，重新分配图像的灰度值，使一定范围内图像的灰度值大致相等。这样，原来直方图中间的峰值部分对比度得到增强，而两侧的谷底部分对比度降低，输出图像的直方图是一个较为平坦的直方图。
直方图的均衡化实际也是一种灰度的变换过程，通常均衡化选择的变换函数是灰度的累积概率
全局直方图均衡化def equalizeHist_map(img_gray):#输入必须是灰度图像
    h, w = img_gray.shape
#计算灰度直方图
    gray_hist = calGrayHist(img_gray)
#计算累加灰度直方图
cumulative_hist = np.zeros([256], np.uint64)
for p in range(256):
if p == 0:
cumulative_hist[p] = gray_hist[p]
else:
cumulative_hist[p] = cumulative_hist[p-1] + gray_hist[p]
#根据累加灰度直方图得到输入灰度级和输出灰度级之间的映射关系
relative_map = np.zeros([256], np.uint8)
coff = 256.0/(h*w)
for p in range(256):
q = coff * float(cumulative_hist[p]) -1
if q < 0:
relative_map[p] = 0
else:
relative_map[p] = math.floor(q)
return relative_map
img = cv2.imread('test.png')
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
relative_map = equalizeHist_map(img_gray)
img_result = np.zeros(img.shape, np.uint8)
h,w = img.shape[:2]
for i in range(h):
for j in range(w):
        img_result[i][j] = relative_map[img[i][j]]
虽然全局直方图均衡化方法对提高对比度很有效，但是均衡化处理以后暗区域的噪声可能会被放大，变得清晰可见，而亮区域可能会损失信息。为了解决该问题，提出了自适应直方图均衡化(Aptive Histogram Equalization)方法，自适应直方图均衡化首先将图像划分为不重叠的区域块，然后对每一个块分别进行直方图均衡化。
线性变换图像增强的缺点：需要根据不同的应用及图像自身的信息进行合理的选择，可能需要进行多次测试，所以选择合适的参数是相当麻烦的。有时候，需要图像具有某一特定的直方图形状（也就是灰度分布），而不是均匀分布的直方图。直方图正规化也叫做直方图匹配，用于将图像变换为某一特定的灰度分布，也就是其目的的灰度直方图是已知的。这其实和均衡化很类似，均衡化后的灰度直方图也是已知的，是一个均匀分布的直方图；而规定化后的直方图可以随意的指定，直方图正规化就是基于当前图像情况自动选取a和b的值的方法。
————————————————
部分参考自：https://blog.csdn.net/victory_tc/article/details/97361911和https://www.cnblogs.com/wangguchangqing/p/7098213.html

1灰度直方图
       1.1 概念

灰度直方图是关于灰度级分布的函数，是对图像中灰度级分布的统计。
灰度直方图是将数字图像中的所有像素，按照灰度值的大小，统计其出现的频率。
灰度直方图是灰度级的函数，它表示图像中具有某种灰度级的像素的个数，反映了图像中某种灰度出现的频率。
  
如果将图像总像素亮度（灰度级别）看成是一个随机变量，则其分布情况就反映了图像的统计特性，这可用probability density function （PDF）来刻画和描述，表现为灰度直方图。



       1.2 opencv表示
       opencv提供了calcHist函数来计算图像直方图。

         其中C++的函数原型如下：
void calcHist(const Mat* arrays, int narrays, const int* channels, InputArray mask, OutputArray hist, int dims, const int* histSize, const float** ranges, bool
 uniform=true, bool accumulate=false )；
         参数解释：
        arrays:输入的图像的指针，可以是多幅图像，所有的图像必须有同样的深度（CV_8U or CV_32F）。同时一副图像可以有多个channes。
        narrays:输入的图像的个数。
        channels:用来计算直方图的channes的数组。比如输入是2副图像，第一副图像有0，1，2共三个channel，第二幅图像只有0一个channel，那么输入就一共有4个channes，如果int channels[3] = {3, 2, 0}，那么就表示是使用第二副图像的第一个通道和第一副图像的第2和第0个通道来计算直方图。
        mask:掩码。如果mask不为空，那么它必须是一个8位（CV_8U）的数组，并且它的大小的和arrays[i]的大小相同，值为1的点将用来计算直方图。
        hist:计算出来的直方图
        dims:计算出来的直方图的维数。
        histSize:在每一维上直方图的个数。简单把直方图看作一个一个的竖条的话，就是每一维上竖条的个数。
        ranges:用来进行统计的范围。比如 float rang1[] = {0, 20}; float rang2[] = {30, 40};  const float *rangs[] = {rang1, rang2};那么就是对0，20和30，40范围的值进行统计。
       uniform:每一个竖条的宽度是否相等。
       accumulate:  是否累加。如果为true，在下次计算的时候不会首先清空hist。




2 直方图均衡化
      1.1 概念


          如果一副图像的像素占有很多的灰度级而且分布均匀，那么这样的图像往往有高对比度和多变的灰度色调。
         直方图均衡化就是一种能仅靠输入图像直方图信息自动达到这种效果的变换函数。
         它的基本思想是对图像中像素个数多的灰度级进行展宽，而对图像中像素个数少的灰度进行压缩，从而扩展像原取值的动态范围，提高了对比度和灰度色调的变化，使图像更加清晰。
        直方图均衡化处理的“中心思想”是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。直方图均衡化就是对图像进行非线性拉伸，重新分配图像像素值，使一定灰度范围内的像素数量大致相同。直方图均衡化就是把给定图像的直方图分布改变成“均匀”分布直方图分布。
        直方图均衡化的基本思想是把原始图的直方图变换为均匀分布的形式，这样就增加了像素灰度值的动态范围从而可达到增强图像整体对比度的效果。
    2.2 opencv表示

               cv::equalizeHist(img_gray,equalize_Hist)


3.例程
int img_Hist(Mat& image)
{
    if(!image.data)
    {
        cout << "fail to load image" << endl;
        return 0;
    }
    Mat img_gray;
    //GRAY
    if(image.channels()==3)
    {
    	cvtColor(image, img_gray, CV_BGR2GRAY);
    }
    else
    {
    	image.copyTo(img_gray);
    }
    cv::imwrite("img_gray.jpg",img_gray);

    MatND hist;       // 在cv中用CvHistogram *hist = cvCreateHist
    int dims = 1;
    float hranges[] = {0, 255};
    const float *ranges[] = {hranges};   // 这里需要为const类型
    int size = 256;
    int channels = 0;
    // 计算图像的直方图
    calcHist(&img_gray, 1, &channels, Mat(), hist, dims, &size, ranges);    // cv 中是cvCalcHist
    int scale = 1;
    Mat imageShow(size * scale, size, CV_8U, Scalar(0));
    // 获取最大值和最小值
    double minVal = 0;
    double maxVal = 0;
    minMaxLoc(hist,&minVal, &maxVal, 0, 0);  //  cv中用的是cvGetMinMaxHistValue
    //显示直方图的图像
    int hpt = saturate_cast<int>(0.9 * size);

    for(int i = 0; i < 256; i++)
    {
        float value = hist.at<float>(i);           //   注意hist中是float类型    cv中用cvQueryHistValue_1D
        int realValue = saturate_cast<int>(value * hpt/maxVal);
        rectangle(imageShow,Point(i*scale, size - 1), Point((i+1)*scale - 1, size - realValue), Scalar(255));
    }
    namedWindow("Hist");
    imshow("Hist", imageShow);
    cv::imwrite("hist.jpg",imageShow);
    Mat equalize_Hist;
    cv::equalizeHist(img_gray,equalize_Hist);

    namedWindow("equalize_Hist");
    imshow("equalize_Hist", equalize_Hist);
    cv::imwrite("equalize_Hist.jpg",equalize_Hist);
    // 计算图像的直方图
       calcHist(&equalize_Hist, 1, &channels, Mat(), hist, dims, &size, ranges);    // cv 中是cvCalcHist
       Mat imageShow_equal(size * scale, size, CV_8U, Scalar(0));
       // 获取最大值和最小值
       minMaxLoc(hist,&minVal, &maxVal, 0, 0);  //  cv中用的是cvGetMinMaxHistValue
       //显示直方图的图像
       hpt = saturate_cast<int>(0.9 * size);
       for(int i = 0; i < 256; i++)
       {
           float value = hist.at<float>(i);           //   注意hist中是float类型    cv中用cvQueryHistValue_1D
           int realValue = saturate_cast<int>(value * hpt/maxVal);
           rectangle(imageShow_equal,Point(i*scale, size - 1), Point((i+1)*scale - 1, size - realValue), Scalar(255));
       }

       namedWindow("Hist_equalize");
       imshow("Hist_equalize", imageShow_equal);
       cv::imwrite("Hist_equalize.jpg",imageShow_equal);
    waitKey(0);
    return 0;
}
int main (int args, char** argv)
{
	Mat image = imread("/home/odroid/TEST/111.jpg", 1);    // 这里也可以是BGR 但是想想提取轮廓 效果是一样的
 	imshow("original", image);
 	img_Hist(image);
	waitKey();
    return 0;
}





             
        



                            原始图                                                                                                                       灰度图


                                                                         
  



灰度直方图                                                                                                                             均衡化后的灰度直方图



                                                                    








                        直方图均衡化后的灰度图