【OpenCV】【openMP】使用openMP对sobel边缘检测算子加速研究

Sobel 算子常用于图像的边缘检测，计算公式如下

OpenMP加速方法计算尺寸为2304X2304，8位灰度掌纹图像的梯度图（或自行选自其他图像），计算采用OpenMP带来的加速比。
图像信息：

原图像：

步骤：
（1）读取图像，转化为Mat矩阵，src为原图像

原图像：

（2）对原图像进行横向运算，找出纵向边缘

Gx图像：

（3）对原图像进行纵向运算，找出横向边缘

Gy图像：

（4）求G，判决门限选择50

G图像：

Threshold=70 Threshold=60

（5）比较运行时间
由于求Gx和Gy时，运算均与中心点（i，j）周围的8个点有关，不能private（tepe），又因为不是累加，不能用reduction，只是在最外围循环使用#pragma ompparallel for，比较时间后，时间并没有大幅度减少，可见并行运算对卷积形式的运算没有多少帮助。
在求G时，由于num变量只和i，j有关，所以可以使用private（num），运算时间相对于前两者运算时间减少的更多，但是不明显。
未并行运算时间：

优化后：

附程序：

#include <opencv2/core/core.hpp>    
#include <opencv2/highgui/highgui.hpp>   
#include <opencv2\opencv.hpp>  
#include <iostream>  
#include <string>  
#include <sstream>  
#include <windows.h>
//#include <atlimage.h> 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using namespace cv;
using namespace std;

#define width 2304
#define height  2304
//时间节点
clock_t t1, t2,t3,t4,t5,t6,tx,ty,txy;
int main()
{
	t1 = clock();//时间节点
	Mat src_Mat,Gx,Gy,Gxy;
	Mat src,dst,temp_G,temp_Gy,temp_Gxy;
	size_t image_size = width * height;
	Size re_size(512,512);
	unsigned char *Data = new unsigned char[image_size];
	FILE *file;
	fopen_s(&file,"palm.raw", "rb+");
	fread(Data, sizeof(unsigned char), image_size, file);
	fclose(file);
	cv::Mat temp(height, width, CV_8UC1, Data);    //单通道的Mat raw数据
	src = temp.clone();
	//初始化矩阵
	Gx = Mat::zeros(width, height, src.type());
	Gy = Mat::zeros(width, height, src.type());
	Gxy = Mat::zeros(width, height, src.type());
	resize(temp, temp, re_size, INTER_LINEAR);
	imshow("src_image", temp);

	int tepe=0;
	t2 = clock();//时间节点
#pragma ompparallel for //并行计算
	//Sobel横向边缘检测
	for (int i = 0; i < src.rows; i++) {
		for (int j = 0; j < src.cols; j++) {
			//原图像与横向矩阵相乘相加
			if ((i - 1 >= 0) && (j - 1 >= 0) && (i + 1) < src.rows && (j + 1) < src.cols) {
				tepe = src.at<uchar>(i - 1, j - 1)*(-1) + src.at<uchar>(i - 1, j)*(-2) +
					src.at<uchar>(i - 1, j + 1)*(-1) + src.at<uchar>(i + 1, j - 1) * 1 
					+ src.at<uchar>(i + 1, j + 1) * 1 +
					src.at<uchar>(i + 1, j) * 2;
				if (tepe > 255)
					tepe = 255;
				if (tepe < 0)
					tepe = 0;
				Gx.at<uchar>(i, j) = (uchar)tepe;
			}
			else {//边缘赋值
				Gx.at<uchar>(i, j) = src.at<uchar>(i, j);
			}
			tepe = 0;
		}
	}
	tx = clock();//时间节点
	
	temp_G = Gx.clone();
	//裁剪图像为512*512，方便显示
	resize(temp_G, temp_G, re_size, INTER_LINEAR);
	imshow("Gx", temp_G);//显示图像

	t3 = clock();//时间节点
#pragma ompparallel for //并行计算
	//Sobel纵向边缘检测
	for (int i = 0; i < src.rows; i++) {
		for (int j = 0; j < src.cols; j++) {
			//原图像与纵向矩阵相乘相加
			if ((i - 1) >= 0 && (j - 1) >= 0 && (i + 1) < src.rows && (j + 1) < src.cols) {
				tepe = src.at<uchar>(i - 1, j - 1)*(-1) + src.at<uchar>(i, j - 1)*(-2) +
					src.at<uchar>(i - 1, j + 1) * 1 + src.at<uchar>(i + 1, j - 1)*(-1) +
					src.at<uchar>(i + 1, j + 1) * 1 +
					src.at<uchar>(i, j + 1) * 2;
				if (tepe > 255)
					tepe = 255;
				if (tepe < 0)
					tepe = 0;
				Gy.at<uchar>(i, j) = (uchar)tepe;
			}
			else {//边缘赋值
				Gy.at<uchar>(i, j) = src.at<uchar>(i, j);
			}
			tepe = 0;
		}
	}
	ty = clock();//时间节点
	
	temp_G = Gy.clone();
	//裁剪图像为512*512，方便显示
	resize(temp_G, temp_G, re_size, INTER_LINEAR);
	imshow("Gy", temp_G);//显示图像

	t4 = clock();//时间节点
	double num;
	//计算G
	for (int i = 0; i < src.rows; i++){
#pragma ompparallel for private(num)//并行计算
		for (int j = 0; j < src.cols; j++) {
			//计算G
			num = sqrt(Gx.at<uchar>(i, j)*Gx.at<uchar>(i, j) + 
				Gy.at<uchar>(i, j)*Gy.at<uchar>(i, j));
			//阈值为50，二值化
			if (num < 50) {
				Gxy.at<uchar>(i, j) = 0;
			}
			else {
				Gxy.at<uchar>(i, j) = 255;
			}
			num = 0.0;			
		}
	}
	txy = clock();	//时间节点
	temp_G = Gxy.clone();
	//裁剪图像为512*512，方便显示
	resize(temp_G, temp_G, re_size, INTER_LINEAR);
	imshow("SobelDetect", temp_G);//显示图像
	t5 = clock();//时间节点
	//计算时间
	printf("t1 =%d\n", t1);
	printf("t2 =%d\n", t2);
	printf("t3 =%d\n", t3);
	printf("t4 =%d\n", t4);
	printf("t5 =%d\n", t5);
	//Gx运算所需时间
	printf("SobelDetect_rows time =%d\n", tx - t2);
	//Gy运算所需时间
	printf("SobelDetect_cols time =%d\n", ty - t3);
	//Gxy运算所需时间
	printf("SobelDetect_sqrt time =%d\n", txy - t4);
	//程序运行时间
	printf("Run time =%d\n", t5 - t1);
	waitKey();
	return 0;
}

大家如果有什么疑问，欢迎回复提问。