直方图处理

  • Post author:
  • Post category:其他


#include"stdafx.h"

Mat histogramImage(Mat srcImage); //输出直方图
int equalization(Mat srcImage); //均衡化
int gradient(Mat srcImage); //梯度模图
Mat calHist(Mat srcImage); //计算直方图
void histMap(Mat& img_src, Mat& img_obj); //单通道规定化
int guidinghua(); //三通道规定化

int main() {
	char ch;
	Mat srcImage = imread("lena.jpg", IMREAD_GRAYSCALE);
	cout << "输入对应数字执行功能:0:退出;1:输出原始图像及直方图;2:输出均衡化后图像及直方图;3:输出规定化后图像及直方图;4:输出梯度的模图像。" << endl;
	while (cin.get(ch)) {
		switch (ch) {
		case '0':
			return 0;
		case '1':
			histogramImage(srcImage);
			break;
		case '2':
			equalization(srcImage);
			break;
		case '3':
			guidinghua();
			break;
		case '4':
			gradient(srcImage);
			break;
		default:
			cout << "输入对应数字执行功能:0:退出;1:输出原始图像及直方图;2:输出均衡化后图像及直方图;3:输出规定化后图像及直方图;4:输出梯度的模图像。" << endl;
			break;
		}
	}
	
}

//画直方图
 Mat histogramImage(Mat srcImage) {
	Mat Hist(256, 1, CV_16UC1, Scalar(0));  //直方图
	int index = 0;
	int width = srcImage.cols;
	int height = srcImage.rows;
	int i, kernel;
	
	omp_set_num_threads(THREAD_NUMS);
	Mat sum(256, THREAD_NUMS, CV_16UC1, Scalar(0));     //并行直方图

	//resize(srcImage, srcImage, Size(height * 10, width * 10), 0, 0, INTER_LINEAR);
    width = srcImage.cols;
	height = srcImage.rows;

	double start = omp_get_wtime();

#pragma omp parallel private(kernel)
    {
		kernel = omp_get_thread_num();
	    for (int i = kernel; i < height; i+=THREAD_NUMS)
	    {
		    for (int j = 0; j < width; j++)
		    {
			    index = srcImage.at<uchar>(i, j);
			    sum.at<ushort>(index, kernel)++;
		    }
	    }
    }

	for (int m = 0; m < 256; m++)
	{
		int a = 0;
    #pragma omp parallel for reduction(+:a)
		for(int j=0;j<THREAD_NUMS;j++)
		   {

			a += sum.at<ushort>(m, j);
		   }
		Hist.at<ushort>(m, 0) = a;
	}

	double maxVal = 0;
	int scale = 2; //每个直方图单位宽度2
	minMaxLoc(Hist, NULL, &maxVal, NULL, NULL);   //最大元素
	int hist_height = 256; 
	Mat hist_img = Mat::zeros(hist_height, 256*scale, CV_8UC3);
	
#pragma omp parallel for
	for (int i = 0; i < 256; i++) 
	{
		float value = Hist.at<ushort>(i, 0); 
		int intensity = cvRound(value * hist_height / maxVal); //最大高度
		rectangle(hist_img, Point(i*scale, hist_height - 1), Point((i + 1)*scale-1, hist_height - intensity), Scalar(255, 255, 255));
	}

	double end = omp_get_wtime();
	cout << "绘制直方图时长为:" << end - start << "秒." << endl;
	imshow("image", srcImage);
	imshow("hist", hist_img);
	waitKey(0);
	//cout << Hist << endl;
	//cout << Hist.size << endl;

	

	return Hist;
 }

//均衡化
int equalization(Mat srcImage)
{
	int height = srcImage.rows;
	int width = srcImage.cols;
	Mat hist;  //原图直方图
	Mat hist1; //均衡化直方图

	int i, kernel;
	
	omp_set_num_threads(THREAD_NUMS);
	Mat sum(256, THREAD_NUMS, CV_16UC1, Scalar(0));

    height = srcImage.rows;
	width = srcImage.cols;
	int size = height * width;

	Mat hist_sum = Mat::zeros(256, 1, CV_64F);  //累计直方图 
	Mat hist_image = Mat::zeros(height, width, CV_64F);  //均衡化图片

	hist = calHist(srcImage);

	double start1 = omp_get_wtime();

	hist_sum.at<double>(0, 0) = hist.at<ushort>(0.0);
	//累计直方图
	for (int i = 1; i < 256; i++)
	{
		int index = 0;
		index = hist.at<ushort>(i, 0);
		hist_sum.at<double>(i, 0) = hist_sum.at<double>(i - 1, 0) + index;
	}

	hist_sum = hist_sum / size;

//均衡化
#pragma omp parallel for 
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
		{
			float value = 0.0;
			int index1 = 0;
			index1 = srcImage.at<uchar>(i, j);
			value = hist_sum.at<double>(index1, 0);
			value = round(255 * value);
			hist_image.at<double>(i,j) = value;
		}
	}
	double end1 = omp_get_wtime();
	
	cout << "均衡化过程时长为:" << end1 - start1 << "秒." << endl;

	hist_image.convertTo(hist_image, CV_8U);
	hist1 = histogramImage(hist_image);

	waitKey(0);
	return 0;

}

//图像梯度、模、模的总和
int gradient(Mat srcImage)
{

	omp_set_num_threads(8);  //设置运行核数
	int height = srcImage.cols;
	int width = srcImage.rows;
	srcImage.convertTo(srcImage, CV_16SC1, 1, 0);

	height = srcImage.cols;
	width = srcImage.rows;

	Mat x_gradient(height, width, CV_16SC1);
	Mat y_gradient(height, width, CV_16SC1);
	Mat abs_gradient(height, width, CV_32FC1);

	srcImage.row(height - 1).copyTo(x_gradient.row(height - 1));
	srcImage.col(width - 1).copyTo(y_gradient.col(width - 1));


	double start1 = omp_get_wtime();

#pragma omp parallel for 
	for (int i = 0; i < height - 1; i++)
	{
		for (int j = 0; j < width; j++)
		{
			x_gradient.at<short>(i, j) = srcImage.at<short>(i + 1, j) - srcImage.at<short>(i, j);
		}
	}


#pragma omp parallel for 
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width - 1; j++)
		{
			y_gradient.at<short>(i, j) = srcImage.at<short>(i, j + 1) - srcImage.at<short>(i, j);
		}
	}

#pragma omp parallel for 
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
			abs_gradient.at<float>(i, j) = sqrt(pow(x_gradient.at<short>(i, j), 2) + pow(y_gradient.at<short>(i, j), 2));
	}


	float sum = 0.0;

#pragma omp parallel for  reduction(+:sum)
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
			sum += abs_gradient.at<float>(i, j);
	}

	double end1 = omp_get_wtime();



	//int n = omp_get_num_procs();
	//cout << n << endl;

	x_gradient.convertTo(x_gradient, CV_8U);
	y_gradient.convertTo(y_gradient, CV_8U);
	abs_gradient.convertTo(abs_gradient, CV_8U);

	cout << "总时间为  " << end1 - start1 << endl;
	cout << "梯度模总和为:" << sum << endl;

	imshow("y_gradient", y_gradient);
	imshow("x_gradient", x_gradient);
	imshow("abs_gradient", abs_gradient);
	waitKey(0);
	//cout<<height<<width;
	
	return 0;
}

//对单通道进行规定化
void histMap(Mat& img_src, Mat& img_obj)
{
	Mat gray_hist; //直方图计算的结果
	int i, j; //循环变量
	double gray_temp = 0; //中间结果,用于计算累计直方图
	double totalPixel; //像素总数

	//计算原图像直方图,并归一化到(0, 1)
	gray_hist = calHist(img_src);
	//cout << typeid(gray_hist).name() << endl;
	//cout << gray_hist.size << endl;
	totalPixel = img_src.rows * img_src.cols;
	double srcHist[256];
	for (i = 0; i < 256; i++)
	{
		srcHist[i] = gray_hist.at<ushort>(i) / totalPixel;
	}

	//计算原图像直方图的累计概率 0 ~ 1
	double srcCumHist[256];
	for (i = 0; i < 256; i++)
	{
		gray_temp = 0;
		for (j = 0; j <= i; j++)
		{
			gray_temp += srcHist[j];
		}
		srcCumHist[i] = gray_temp;
	}

	//计算目标图像直方图
	gray_hist = calHist(img_obj);
	totalPixel = img_obj.rows * img_obj.cols;
	double objHist[256];
	for (i = 0; i < 256; i++)
	{
		objHist[i] = gray_hist.at<ushort>(i) / totalPixel;
	}

	//计算目标图像直方图的累计概率 0 ~ 1
	double objCumHist[256];
	for (i = 0; i < 256; i++)
	{
		gray_temp = 0;
		for (j = 0; j <= i; j++)
		{
			gray_temp += objHist[j];
		}
		objCumHist[i] = gray_temp;
	}


	Mat diffCumHist(256, 256, CV_32FC1); //两幅图累计概率的差
	for (int i = 0; i < 256; i++)
	{
		for (int j = 0; j < 256; j++)
			diffCumHist.at<float>(i, j) = fabs(srcCumHist[i] - objCumHist[j]);
	}

	Mat map(256, 1, CV_8U);    //映射关系

#pragma omp parallel for
	for (int i = 0; i < 256; i++)
	{
		float min = diffCumHist.at<float>(i, 0); //每行最小值
		int index = 0; //最小值的索引
		for (int j = 0; j < 256; j++)
		{
			if (min > diffCumHist.at<float>(i, j))
			{
				min = diffCumHist.at<float>(i, j);
				index = j;
			}
		}
		map.at<uchar>(i, 0) = index;
	}

	int height = img_src.cols;
	int width = img_src.rows;

//#pragma omp parallel for
	for (i = 0; i < height; i++)
	{
		for (j = 0; j < width; j++)
		{
			int index = img_src.at<uchar>(i, j);
			img_src.at<uchar>(i, j) = map.at<uchar>(index, 0);
		}

	}
}

//三通道规定化
int guidinghua()
{
	Mat img_src = imread("lena.jpg");
	Mat img_obj = imread("flower.jpeg");
	Mat imgOutput; //规定化后的输出图像

	//分割原图像通道
	vector<Mat> src_channels;
	Mat src_blue, src_green, src_red;
	split(img_src, src_channels);
	src_blue = src_channels.at(0);
	src_green = src_channels.at(1);
	src_red = src_channels.at(2);
	//分割目标图像通道
	vector<Mat> obj_channels;
	Mat obj_blue, obj_green, obj_red;
	split(img_obj, obj_channels);
	obj_blue = obj_channels.at(0);
	obj_green = obj_channels.at(1);
	obj_red = obj_channels.at(2);

	//分别对BGR通道做直方图规定化
	histMap(src_blue, obj_blue);
	histMap(src_green, obj_green);
	histMap(src_red, obj_red);
	//合并通道,输出结果
	merge(src_channels, imgOutput);

	//显示图像
	imshow("img_src", img_src);
	imshow("img_obj", img_obj);
	imshow("imgOutput", imgOutput);
	waitKey(0);

	return 0;
}

//计算直方图
Mat calHist(Mat srcImage) {
	Mat Hist(256, 1, CV_16UC1, Scalar(0));  //直方图
	int index = 0;
	int width = srcImage.cols;
	int height = srcImage.rows;
	int i, kernel;

	omp_set_num_threads(THREAD_NUMS);
	Mat sum(256, THREAD_NUMS, CV_16UC1, Scalar(0));     //并行直方图

	//resize(srcImage, srcImage, Size(height * 10, width * 10), 0, 0, INTER_LINEAR);
	width = srcImage.cols;
	height = srcImage.rows;

	double start = omp_get_wtime();

#pragma omp parallel private(kernel)
	{
		kernel = omp_get_thread_num();
		for (int i = kernel; i < height; i += THREAD_NUMS)
		{
			for (int j = 0; j < width; j++)
			{
				index = srcImage.at<uchar>(i, j);
				sum.at<ushort>(index, kernel)++;
			}
		}
	}

	for (int m = 0; m < 256; m++)
	{
		int a = 0;
#pragma omp parallel for reduction(+:a)
		for (int j = 0; j < THREAD_NUMS; j++)
		{

			a += sum.at<ushort>(m, j);
		}
		Hist.at<ushort>(m, 0) = a;
	}

	double maxVal = 0;
	int scale = 2; //每个直方图单位宽度2
	minMaxLoc(Hist, NULL, &maxVal, NULL, NULL);   //最大元素
	int hist_height = 256;
	Mat hist_img = Mat::zeros(hist_height, 256 * scale, CV_8UC3);

#pragma omp parallel for
	for (int i = 0; i < 256; i++)
	{
		float value = Hist.at<ushort>(i, 0);
		int intensity = cvRound(value * hist_height / maxVal); //最大高度
		rectangle(hist_img, Point(i * scale, hist_height - 1), Point((i + 1) * scale - 1, hist_height - intensity), Scalar(255, 255, 255));
	}

	double end = omp_get_wtime();
	//cout << Hist << endl;
	//cout << Hist.size << endl;
	cout << "绘制图像直方图时长为:" << end - start << "秒." << endl;


	return Hist;
}
#include<iostream>
#include<omp.h>
#include<opencv2/opencv.hpp>
#include<math.h>
#include<vector>

using namespace std;
using namespace cv;

#define THREAD_NUMS 8




版权声明:本文为qq_40922398原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。