【C++】【Intel MKL】【OpenMP】【模糊聚类】多线程英特尔MKL实战 - 模糊C均值图像分割

最新推荐文章于 2024-06-04 16:47:17 发布
最新推荐文章于 2024-06-04 16:47:17 发布
阅读量1.2k 收藏 3
点赞数 1
分类专栏： 基本知识
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/u011861755/article/details/114025187
版权

1、英特尔MKL

面对计算量庞大的线性代数计算，英特尔为其CPU开发了MKL（Math Kernel Library）。Intel MKL能够充分利用CPU指令集和并行计算优势，较大提高矩阵计算速度。
具体参考：https://software.intel.com/content/www/us/en/develop/articles/intel-math-kernel-library-documentation.html

2、概述

本实验采用英特尔MKL与OpenMP完成模糊C均值图像分割。Intel MKL没有为矩阵逐元素乘积开发API。若要使用MKL计算矩阵逐元素乘积，步骤如下：
1、取出两个矩阵的相同一行（列），得到两个行向量A和B。
2、将A转置为列向量，得到 A T A^T AT
3、计算 A T B A^TB ATB，得到一个方阵。
4、取出方阵的对角线元素，即为行向量A和B的逐元素乘积。
然而，第四步仍然需要较多次循环，不如用循环直接完成逐元素相乘。所以本文采用OpenMP多线程循环计算矩阵逐元素乘积。

3、实验

输入图像
采用纯循环实现的FCM：
循环FCM
纯循环FCM

采用Intel MKL+OpenMP多线程
在这里插入图片描述
加速FCM
可以看到计算效率提高了5倍左右，结果没有变化，具有一定提升。

4、循环版FCM代码

#include <opencv2\opencv.hpp>
#include <stdio.h>
#include <iostream>
#include <time.h>
using namespace cv;
using namespace std;
int main()
{
	Mat input = imread("F:\\C++语言程序\\模糊聚类\\Images\\118035.jpg", 0);
	const int rows = input.rows;
	const int cols = input.cols;
	namedWindow("输入");
	imshow("输入", input);
	// 输入参数
	int m = 2, K = 3, max_iter = 100;
	double error = 0.01;
	// 为输入图像开辟内存空间，用二维数组表示
	double** input_ptr = (double**)malloc(rows * sizeof(double*));
	for (int i = 0; i < rows; i++)
	{
		input_ptr[i] = (double*)malloc(cols * sizeof(double));
	}
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			input_ptr[i][j] = (double)input.at<uchar>(i, j);
	 初始化隶属度数组U，用三维数组表示
	// 开辟U的内存
	double*** U = (double***)malloc(rows * sizeof(double**));
	for (int i = 0; i < rows; i++)
	{
		U[i] = (double**)malloc(cols * sizeof(double*));
		for (int j = 0; j < cols; j++)
			U[i][j] = (double*)malloc(K * sizeof(double));
	}
	// 随机初始化U
	srand((int)time(0));
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			for (int l = 0; l < K; l++)
				U[i][j][l] = (double)rand() / RAND_MAX;
	// 令每一个像素的U之和为1
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
		{
			double U_sum = 0;
			for (int l = 0; l < K; l++)
				U_sum += U[i][j][l];
			for (int l = 0; l < K; l++)
				U[i][j][l] /= U_sum;
		}
	// 开辟聚类中心center的内存
	double* center = (double*)malloc(K * sizeof(double));
	// 开辟目标函数J的内存
	double* J = (double*)calloc(max_iter, sizeof(double));
	// 聚类
	clock_t time_start = clock();
	int iter;
	for (iter = 0; iter < max_iter; iter++)
	{
		// 更新聚类中心
		for (int l = 0; l < K; l++)
		{
			double center_up, center_down;
			center_up = center_down = 0;
			for (int i = 0; i < rows; i++)
				for (int j = 0; j < cols; j++)
				{
					center_up += input_ptr[i][j] * pow(U[i][j][l], m);
					center_down += pow(U[i][j][l], m);
				}
			center[l] = center_up / center_down;
		}
		// 更新隶属度
		for (int i = 0; i < rows; i++)
			for (int j = 0; j < cols; j++)
			{
				double U_down = 0;
				for (int l = 0; l < K; l++)
					U_down = U_down + 1 / pow(input_ptr[i][j] - center[l], 2 / (m - 1));
				for (int l = 0; l < K; l++)
					U[i][j][l] = 1 / (pow(input_ptr[i][j] - center[l], 2 / (m - 1)) * U_down);
			}
		// 更新目标函数
		for (int i = 0; i < rows; i++)
			for (int j = 0; j < cols; j++)
				for (int l = 0; l < K; l++)
					J[iter] += pow(U[i][j][l], m) * pow(input_ptr[i][j] - center[l], 2);
		printf("第%d次迭代，J = %.6f\n", iter + 1, J[iter]);
		// 迭代条件判定
		if (abs(J[iter] - J[iter - 1]) <= error && iter >= 1)
		{ 
			cout << "目标函数已收敛\n"; 
			break; 
		}
		else if (iter == max_iter - 1 && iter >= 1)
		{ 
			cout << "目标函数未收敛，达到最大迭代次数\n"; 
			break;
		}
	}
	clock_t time_end = clock();
	printf("聚类耗时%.3f秒\n", (double)(time_end - time_start) / CLOCKS_PER_SEC);
	printf("单次迭代耗时%.3f秒\n", (double)(time_end - time_start) / CLOCKS_PER_SEC / ((double)iter + 1));
	cout << "聚类中心分别为\n";
	for (int l = 0; l < K; l++)
		cout << center[l] << "\n";
	// 计算隶属度划分系数和划分熵
	double Vpc, Vpe;
	Vpc = Vpe = 0;
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			for (int l = 0; l < K; l++)
			{
				Vpc += U[i][j][l] * U[i][j][l];
				Vpe += U[i][j][l] * log(U[i][j][l]);
			}
	Vpc = Vpc / (rows * cols) * 100;
	Vpe = -Vpe / (rows * cols) * 100;
	cout << "隶属度划分系数Vpc = " << Vpc << "%" << endl;
	cout << "隶属度划分熵Vpe = " << Vpe << "%" << endl;
	// 标签数组
	int** labels = (int**)calloc(rows, sizeof(int*));
	for (int i = 0; i < rows; i++)
		labels[i] = (int*)calloc(cols, sizeof(int));
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
		{
			double U_max = 0;
			for (int l = 0; l < K; l++)
				if (U_max <= U[i][j][l])
				{
					U_max = U[i][j][l];
					labels[i][j] = l + 1;
				}
		}
	// 标签数组着色
	Mat result = Mat(rows, cols, CV_8U);
	uchar* colors = new uchar[K];
	RNG rng(time(0));
	for (int i = 0; i < K; i++)
		*(colors + i) = (uchar)(center[i]);
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			for (int l = 1; l <= K; l++)
				if (labels[i][j] == l)
					result.at<uchar>(i, j) = *(colors + l - 1);
	imshow("聚类结果", result);
	waitKey(0);
	return 0;
}

5、Intel MKL+OpenMP版FCM

#include <stdio.h>
#include <iostream>
#include <mkl.h>
#include <opencv2/opencv.hpp>
#include <omp.h>
using namespace std;
using namespace cv;
void power(double** A, double** result, int& rows, int& cols, double& p);
void power(double** A, int& rows, int& cols, double& p);
void equal(double* A, double* result, int& N);
void equal(double& number, double* result, int& N);
void print(double** A, int& rows, int& cols);
void print(double* A, int& N);
void reshape(short* A, short** result, int& N, int& row_num, int& col_num);
void vectorProduct(double* A, double* B, double* result, int& N)
{
#pragma omp parallel for
	for (int i = 0; i < N; i++)
		result[i] = A[i] * B[i];
}
void power(double** A, double** result, int& rows, int& cols, double& p)
{
#pragma omp parallel for collapse(2)
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			result[i][j] = pow(A[i][j], p);
}
void power(double** A, int& rows, int& cols, double& p)
{
#pragma omp parallel for collapse(2)
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			A[i][j] = pow(A[i][j], p);
}
void equal(double* A, double* result, int& N)
{
#pragma omp parallel for
	for (int i = 0; i < N; i++)
		result[i] = A[i];
}
void equal(double& number, double* result, int& N)
{
#pragma omp parallel for
	for (int i = 0; i < N; i++)
		result[i] = number;
}
void print(double** A, int& rows, int& cols)
{
	int rows_dhow = rows;
	int cols_dhow = cols;
	if (rows > 10)
		rows_dhow = 10;
	if (cols > 10)
		cols_dhow = 10;
	cout << '\t';
	for (int j = 0; j < cols_dhow; j++)
		cout << j + 1 << '\t';
	cout << endl;
	for (int i = 0; i < rows_dhow; i++)
	{
		cout << i + 1 << '\t';
		for (int j = 0; j < cols_dhow; j++)
			printf("%.2f\t", A[i][j]);
		cout << endl;
	}
}
void print(double* A, int& N)
{
	int cols_dhow = N;
	if (N > 10)
		cols_dhow = 10;
	for (int j = 0; j < cols_dhow; j++)
		cout << j + 1 << '\t';
	cout << endl;
	for (int j = 0; j < cols_dhow; j++)
		printf("%.2f\t", A[j]);
	cout << endl;
}
void reshape(short* A, short** result, int& N, int& row_num, int& col_num)
{
	if (row_num * col_num != N)
	{
		cout << "元素数量发生改变，函数退出\n";
		return;
	}
	for (int i = 0; i < row_num; i++)
		for (int j = 0; j < col_num; j++)
			result[i][j] = A[i * col_num + j];
}
int main()
{
	int K = 3, max_iter = 100;
	double m = 2, error = 0.01;
	Mat f = imread("F:\\C++语言程序\\模糊聚类\\Images\\118035.jpg", 0);
	imshow("原始图像", f);
	int rows = f.rows;
	int cols = f.cols;
	int N = rows * cols;
	// 将图像一维化
	double* input = new double[N];
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			input[i * cols + j] = (double)f.at<uchar>(i, j);
	f.release();
	double* input_temp = new double[N];
	// 初始化聚类隶属度U
	double** U = new double*[K];
	double** U_m = new double*[K];
	double* U_down = new double[N];
	for (int i = 0; i < K; i++)
	{
		U[i] = new double[N];
		U_m[i] = new double[N];
	}
	srand((int)time(0));
	for (int i = 0; i < K; i++)
		for (int j = 0; j < N; j++)
			U[i][j] = (double)rand() / RAND_MAX;
	for (int i = 0; i < N; i++)
	{
		double U_dum = 0;
		for (int j = 0; j < K; j++)
			U_dum += *(*(U + j) + i);
		for (int j = 0; j < K; j++)
			*(*(U + j) + i) /= U_dum;
	}
	power(U, U_m, K, N, m);
	// 初始化聚类中心center
	double* center = new double[K];
	// 初始化距离
	double** d = new double* [K];
	double** d_rec = new double* [K];
	for (int i = 0; i < K; i++)
	{
		d[i] = new double[N];
		d_rec[i] = new double[N];
	}
	// 初始化目标函数值
	double J = 0, J_old = 0;
	double time_start = dsecnd();
	// 聚类
	int iter;
	for (iter = 0; iter < max_iter; iter++)
	{
		for (int i = 0; i < K; i++)
		{
			// 更新聚类中心
			*(center + i) = cblas_ddot(N, U_m[i], 1, input, 1) / cblas_dasum(N, U_m[i], 1);
			equal(*(center + i), *(d + i), N); // d[:][i] = center[i]
			equal(input, input_temp, N); // input_temp = input
			// 更新距离
			cblas_drot(N, input_temp, 1, *(d + i), 1, -1, -1); // d[i][:] = input_temp - d[i][:]
		}
		equal(input, input_temp, N);
		// 距离平方，得到像素与聚类中心的欧氏距离
		power(d, K, N, m);
		// 计算隶属度
		double p = -1;
		power(d, d_rec, K, N, p);
		double init = 0;
		equal(init, U_down, N);
		for (int i = 0; i < K; i++)
			cblas_drot(N, U_down, 1, d_rec[i], 1, 1, 1); // U_down = U_down + d_rec[i]
		for (int i = 0; i < K; i++)
			vectorProduct(d[i], U_down, U[i], N); // U[i][:] = d[i][:] .* U_down[:]
		power(U, K, N, p);
		power(U, U_m, K, N, m);
		// 计算目标函数
		double J_temp;
		J = 0;
		for (int i = 0; i < K; i++)
		{
			vectorProduct(U_m[i], d[i], d[i], N);
			J_temp = cblas_dasum(N, d[i], 1);
			J += J_temp;
		}
		printf("第%d次聚类, J = %.5f\n", iter + 1, J);
		// 迭代条件判定
		if (iter >= 1)
		{
			if (abs(J - J_old) <= error)
			{
				cout << "目标函数已收敛\n";
				break;
			}
			else if (iter == max_iter - 1)
			{
				cout << "目标函数未收敛，达到最大迭代次数\n";
				break;
			}
		}
		J_old = J;
	}
	double time_end = dsecnd();
	printf("聚类耗时%.3f秒\n", (time_end - time_start));
	printf("单次迭代耗时%.3f秒\n", (time_end - time_start) / ((double)iter + 1));
	// 清理聚类过程中的二维临时数组
	for (int i = 0; i < K; i++)
	{
		delete[]U_m[i];
		U_m[i] = NULL;
		delete[]d[i];
		d[i] = NULL;
		delete[]d_rec[i];
		d_rec[i] = NULL;
	}
	delete[]U_m;
	U_m = NULL;
	delete[]d;
	d = NULL;
	delete[]d_rec;
	d_rec = NULL;
	// 清理聚类过程中的一维临时数组
	delete[]input_temp;
	input_temp = NULL;
	delete[]input;
	input = NULL;
	delete[]U_down;
	U_down = NULL;
	cout << "聚类中心分别为\n";
	print(center, K);
	// 计算隶属度划分系数和划分熵
	double Vpc, Vpe;
	Vpc = Vpe = 0;
	for (int i = 0; i < K; i++)
		for (int j = 0; j < N; j++)
			{
				Vpc += U[i][j] * U[i][j];
				Vpe += U[i][j] * log(U[i][j]);
			}
	Vpc = Vpc / N * 100;
	Vpe = -Vpe / N * 100;
	cout << "隶属度划分系数Vpc = " << Vpc << "%" << endl;
	cout << "隶属度划分熵Vpe = " << Vpe << "%" << endl;
	// 标签数组
	short* labels = new short[N];
	for (int i = 0; i < N; i++)
	{
		double U_max = 0;
		for (int j = 0; j < K; j++)
		{
			if (U_max <= U[j][i])
			{
				U_max = U[j][i];
				labels[i] = (short)j;
			}
		}
	}
	for (int i = 0; i < K; i++)
	{
		delete[]U[i];
		U[i] = NULL;
	}
	delete[]U;
	U = NULL;
	short** labels_reshape = new short* [rows];
	for (int i = 0; i < rows; i++)
		labels_reshape[i] = new short[cols];
	reshape(labels, labels_reshape, N, rows, cols);
	delete[]labels;
	labels = NULL;
	// 标签数组着色
	Mat result = Mat(rows, cols, CV_8U);
	uchar* colors = new uchar[K];
	RNG rng(time(0));
	for (int i = 0; i < K; i++)
		*(colors + i) = (uchar)(center[i]);
	delete[]center;
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < cols; j++)
			for (int l = 0; l < K; l++)
				if (labels_reshape[i][j] == (short)l)
					result.at<uchar>(i, j) = *(colors + l);
	imshow("聚类结果", result);
	waitKey(0);
	return 0;
}
寒霜雨刃
关注
  • 1
    点赞
  • 踩
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
在CLion中使用MKL(C++)
weixin_36665857的博客
10-13 1454
MKL 安装
关于MKL库的多线程设置
wh15542311366的博客
08-09 3872
对于多核程序，多线程对于程序的性能至关重要。 下面，我们将对Intel MKL 有关多线程方面的设置做一些介绍： 我们提到MKL 支持多线程，它包括的两个概念： 1>MKL 是线程安全的： MKL在设计时，就保证它是一个线程安全的库函数。 也就是说，无论是在单个线程中调用MKL函数，还是在多个线程中同时使用Intel MKL 函数，都能够确保函数有正确的计算结果。 2>MKL函数内部实现了多线程优化。许多MKL的函数，已经包括内部多线程的实现。用户调这些函数时，只需设置多线程的数目，MKL
linux多进程编程计算圆周率,多线程并行程序性能分析方法综述之MKL 数学核心函数库...
weixin_39717318的博客
04-30 497
2.3 MKL数学核心函数库Intel® 性能库提供了高度优化的函数，这些函数可充分利用Intel® 多核处理器，从而能够最大限度地获得应用程序的性能并减少开发时间。Intel® 数学核心函数库(MKL)是一套经过高度优化的、线程安全的数学例程、函数，并在此基础上加入了在Linux集群环境下的ScaLAPACK(Scalable LAPACK – http://www.netlib.org/sca...
mkl和openmp同时使用时需要注意的坑
最新发布
myth52066的博客
06-04 140
在单独的线程中,需要调用mkl库的mkl_set_num_threads_local函数设置并行的线程数量:int policy;while (1)i < 10;++i)
MKL -C++基本操作
04-18
MKL基本操作的C++实现
MKL学习——基本操作C++实现
热门推荐
风翼冰舟的博客
04-17 1万+
前言前面介绍了各种向量-向量，矩阵-向量，矩阵-矩阵的函数简介。根据自身目前状况，主要使用实数域的操作，也就是说关注单精度float类型的s和双精度double类型的d。还有就是用的基本都是全矩阵，没有经过压缩，也不是对称、三角、带状的某一种情况。所以主要还是总结一般的乘法、加法操作。【注】代码都以单精度float的情况书写，主要流程要记住，使用mkl_malloc申请内存，使用mkl_free释放
【C++】【图像分割】【模糊聚类】FLICM算法C++实现，OpenCV负责图像读取与显示
【原创】【论文复现】【Matlab】【图像处理】基于三维直方图降维重建与高斯分布拟合的Otsu全局阈值分割算法的优化
02-19 1008
S. Krinidis, V. Chatzis. A Robust Fuzzy Local Information C-Means Clustering Algorithm[J]. IEEE Transactions on Image Processing, 19(5), 2010: 1328-1337. 该算法推导有误，使目标函数无法最小化。 正确的推导请参考： T. Celik, H. K. Lee. Comments on “A Robust Fuzzy Local Information C-Mea
Openmp和MKL的一点兼容性问题：VS2010+MKL10.0+OpenMP2.0
战斗蜗牛的专栏
11-18 5604
发现这个问题网上比较少碰到，只有一个帖子ti'j
intel-mkl-src:重新分配英特尔MKL板条箱
05-15
英特尔-mkl-src 箱crate.io 描述英特尔-mkl-src 英特尔MKL的源板条箱英特尔-mkl-sys 用于Intel-MKL 和FFI 英特尔-mkl工具 用于重新分发Intel-MKL的CLI实用程序将英特尔MKL重新包装为板条箱。 在Linux，macOS和...
MKL-数据并行C++开发人员参考
09-13
《Intel oneAPI Math Kernel Library - 数据并行C++开发人员参考》 Intel oneAPI Math Kernel Library（oneMKL）是Intel公司提供的一款高效能数学核心函数库，它旨在提升大型计算问题的软件应用性能。oneMKL包含了...
numpy-1.19.1+mkl-cp39-cp39-win_amd64.whl
08-11
numpy-1.19.1+mkl-cp39-cp39-win_amd64.whl python第三方库 64位操作系统
MKL入门教程
04-12
MKL初级tutorial
MKL官方文档
04-12
MKL2017官方开发文档，可以用于当做工具书
VC++用MKL进行最优化计算的示例源程序
12-08
一个用Visual C++ 2012 , intel MKL 进行无约束最优化计算的示例源程序。
mkl-2018.1-1.4.1-API文档-中英对照版.zip
06-18
赠送jar包：mkl-2018.1-1.4.1.jar； 赠送原API文档：mkl-2018.1-1.4.1-javadoc.jar； 赠送源代码：mkl-2018.1-1.4.1-sources.jar； 赠送Maven依赖信息文件：mkl-2018.1-1.4.1.pom； 包含翻译后的API文档：mkl-2018.1...
Tensorflow2.6.0-MKL for C++
wd1603926823的专栏
11-05 1690
在intel官网看到Optimize AI Applications with Intel® oneAPI Deep Neural Network... 查了一下，好像即使只要有intel，以后CPU预测速度也可以提高。所以决定试试，去年我编译的libtensorflow.so好像没有加入这一块！ 试一下，看是否比去年的快。 ...
C++ 科学计算矩阵数学库: Intel MKL+ Dlib +Armadillo 编程环境配置安装以及使用matlab将M文件编译成链接库
wongpp的博客
11-13 2986
dlib +matlab+ QT +armadillo, C++
C++ & Intel MKL 混合编程
Santorinisu Blogs
05-10 7394
1. 软件版本信息Windows 10Visual Studio 2015 ProfessionalIntel MKL2. 软件来源链接Intel MKL下载链接：方式 1：直接进入Intel Software官网，进行注册后下载Intel Math Kernel Library(MKL) Package.下载链接：https://software.intel.com/en-us/performa...
关于C/C++以及Fortran使用MKL函数库编译环境的配置方案
crysky7ye的博客
01-13 2567
关于C/C++以及Fortran使用MKL函数库编译环境的配置方案 安装vs2019方法自行百度。 安装Intel Parallel Studio XE 2019方法自行百度，强调一下该软件自带MKL函数库。 C/C++方案配置： 点击工程属性：Test222属性（P） C/C++调用MKL函数库的配置方案如下： 此处设置是选用对应的MKL配置模式，由程序自动设置包含的头文件...
Intel MKL官方编程手册
08-05
11.MKL 的发展历史：Intel MKL 自 2003 年发布以来，已经历经了多个版本的更新和改进，提供了更多的功能和工具。 12.MKL 的支持资源：Intel MKL 提供了多种支持资源，包括文档、示例代码、社区支持等，以帮助开发者...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值