加速图像处理的神器: INTEL ISPC 编译器迁移图像旋转算法 - 从 ISPC双精度到 ISPC单精度 (四)

本文探讨了将ISPC双精度计算转换为单精度计算的方法,并对比了性能提升。通过实验发现,在不进行额外内存访问优化的情况下,单精度计算使性能提升了1.5倍。

从 ISPC双精度到 ISPC单精度

前面把原始的C代码转成了ISPC可编译的C代码,其中image_rotate_double_ispc函数里面的数据都是基于double双精度来运算的。我的电脑是支持AVX/AVX2指令集的,所以一次可以并发做4个double浮点的运算,理论上可以提升4倍的算力。

在这里插入图片描述

通过ISPC的编译,实际获得了3.74X倍的性能加速。

从上图YMM寄存器的宽度和浮点数据的宽度来看,YMM寄存器可以一次做8个单精float型数据的计算。这次就来试试把image_rotate_double_ispc函数里面的计算全部改为单精浮点运算,看看性能有多少提升。

代码改动

  1. 把所有的double替换成float
  2. 定义float常量的后缀是f, 定义double常量的后缀是d
#define M_PI_F 3.1415926535f
 
export void image_rotate_float_ispc(uniform const uint8 srcImg[], uniform uint8 dstImg[], uniform float center_x,uniform float center_y, uniform int Width, uniform int Height, uniform float RotateDegree)
{
	uniform float angle = (float)RotateDegree*M_PI_F / 180.0;
	uniform float alpha = cos(angle);
	uniform float beta = sin(angle);
	uniform float m[6];
 
	m[0] = alpha;
	m[1] = -beta;
	m[2] = (1.0 - alpha) * (float)center_x + beta * (float)center_y ;
	m[3] = beta;
	m[4] = alpha;
	m[5] = (1.0 - alpha) * (float)center_y - beta * (float)center_x;
 
	for (uniform int row = 0; row < Height; row++)
		foreach (col = 0 ... Width) {
			float x, y;
			int leftX, rightX, topY, bottomY;
			float w00, w01, w10, w11;
			float fxy;
 
			x = m[0] * (float)col + m[1] * (float)row + m[2];
			y = m[3] * (float)col + m[4] * (float)row + m[5];
 
			leftX = floor(x);
			topY = floor(y);
			rightX = leftX + 1.0;
			bottomY = topY + 1.0;
 
			w11 = abs(x - leftX)*abs(y - topY);
			w01 = abs(1.0 - (x - leftX))*abs(y - topY);
			w10 = abs(x - leftX)*abs(1 - (y - topY));
			w00 = abs(1.0 - (x - leftX))*abs(1.0 - (y - topY));
 
			if ((int)leftX >= 0 && (int)rightX < Width && (int)topY >= 0 && (int)bottomY < Height) {
				fxy = (float)srcImg[topY*Width+ leftX]*w00 +    
					  (float)srcImg[bottomY*Width+ leftX]*w01 +
					  (float)srcImg[topY*Width+ rightX]*w10 + 
					  (float)srcImg[bottomY*Width+ rightX]*w11;
 
				fxy = round(fxy);
				if (fxy < 0)
					fxy = 0;
				if (fxy > 255)
					fxy = 255;
 
				dstImg[row*Width+ col] = (uint8)(fxy);
			}
			else
				dstImg[row*Width + col] = 0;
		};
};

代码对比

image_rotate_double_ispc()和image_rotate_float_ispc()
在这里插入图片描述

运行一下单精度计算版本的代码,运行时间: 761ms

性能提升

和image_rotate_double_ispc版本的性能提升

1148ms/761ms = 1.5X

和原始C代码的性能比对

4301ms/761ms = 5.65X

5.65X虽然和理论的8X有很大的差距,但是已经很满意了,毕竟这段代码写的时候没有考虑任何内存访问的优化,现在能有这个性能提升很好很暴力!!!

% 主函数:处理GAF、MTF以及融合 function main() % 设置输入和输出目录 inputDir = '.'; % 当前目录作为输入目录 gafOutputDir = 'GAF_Images'; mtfOutputDir = 'MTF_Images'; fusionOutputDir = 'Fused_Images'; % 创建输出目录 if ~exist(gafOutputDir, 'dir') mkdir(gafOutputDir); end if ~exist(mtfOutputDir, 'dir') mkdir(mtfOutputDir); end if ~exist(fusionOutputDir, 'dir') mkdir(fusionOutputDir); end % 线路名称列表 line_names = {'L1', 'L2', 'L3', 'L4', 'L5'}; % 设置参数 num_bins = 8; % MTF分箱数量 max_points = 1000; % 最大点数(降采样) % 启动并行池(加速处理) if isempty(gcp('nocreate')) parpool(); end fprintf('已启动并行池,使用 %d 个工作进程\n', gcp().NumWorkers); % 使用并行循环处理线路 parfor i = 1:length(line_names) line_name = line_names{i}; filename = [line_name '.mat']; % 检查文件是否存在 if ~exist(filename, 'file') warning('文件不存在: %s', filename); continue; end % 加载数据 data_struct = load(filename); % 获取数据(假设变量名与文件名相同) if ~isfield(data_struct, line_name) warning('文件 %s 中未找到变量 %s', filename, line_name); continue; end ts_obj = data_struct.(line_name); % 处理timeseries对象 if isa(ts_obj, 'timeseries') ts_data = ts_obj.Data; time_vector = ts_obj.Time; else ts_data = ts_obj; time_vector = []; end % 确保数据为一维 ts_data = ts_data(:)'; % 转换为行向量 % 降采样(减少计算量) n = length(ts_data); if n > max_points factor = ceil(n / max_points); ts_data = downsample(ts_data, factor); if ~isempty(time_vector) time_vector = downsample(time_vector, factor); end fprintf('%s: 降采样 %d -> %d 点\n', line_name, n, length(ts_data)); end % 数据预处理:归一化到[-1,1]区间 normalized_data = normalizeData(ts_data); % 计算GAF矩阵 gasf = computeGASF_optimized(normalized_data); gadf = computeGADF_optimized(normalized_data); % 计算MTF矩阵 mtf_matrix = computeMTF(normalized_data, num_bins); % 保存GAF、MTF图像 saveGAFImage(line_name, gasf, gafOutputDir); saveMTFImage(line_name, mtf_matrix, time_vector, mtfOutputDir); % 融合GAF和MTF fused_image = fuseGAFandMTF(gasf, mtf_matrix, 'weighted'); saveFusedImage(line_name, fused_image, fusionOutputDir); fprintf('已完成: %s\n', line_name); end fprintf('\n所有线路处理完成!\n'); fprintf('GAF图像保存在: %s\n', fullfile(pwd, gafOutputDir)); fprintf('MTF图像保存在: %s\n', fullfile(pwd, mtfOutputDir)); fprintf('融合图像保存在: %s\n', fullfile(pwd, fusionOutputDir)); % 尝试打开输出目录 try if ispc winopen(gafOutputDir); winopen(mtfOutputDir); winopen(fusionOutputDir); elseif ismac system(['open "' gafOutputDir '"']); system(['open "' mtfOutputDir '"']); system(['open "' fusionOutputDir '"']); else system(['xdg-open "' gafOutputDir '"']); system(['xdg-open "' mtfOutputDir '"']); system(['xdg-open "' fusionOutputDir '"']); end catch fprintf('请手动访问目录:\n'); fprintf('GAF图像目录: %s\n', gafOutputDir); fprintf('MTF图像目录: %s\n', mtfOutputDir); fprintf('融合图像目录: %s\n', fusionOutputDir); end end % 数据归一化函数 function normalized = normalizeData(data) data = double(data(:)); % 转换为双精度列向量 min_val = min(data); max_val = max(data); if min_val == max_val normalized = zeros(size(data)); else normalized = (data - min_val) / (max_val - min_val); end end % 计算GASF的优化函数 function gasf = computeGASF_optimized(data) data = double(data(:)); phi = acos(data); gasf = cos(bsxfun(@plus, phi, phi')); end % 计算GADF的优化函数 function gadf = computeGADF_optimized(data) data = double(data(:)); phi = acos(data); gadf = sin(bsxfun(@minus, phi, phi')); end % 计算MTF的函数 function mtf = computeMTF(data, num_bins) % 数据归一化 data = normalizeData(data); % 计算分位数边界 quantiles = quantile(data, num_bins-1); bin_edges = [-inf, quantiles, inf]; % 离散化数据 [~, ~, bin_indices] = histcounts(data, bin_edges); % 计算马尔可夫转移矩阵 n = length(bin_indices); transition_matrix = zeros(num_bins, num_bins); for k = 1:(n-1) i = bin_indices(k); j = bin_indices(k+1); transition_matrix(i, j) = transition_matrix(i, j) + 1; end % 归一化转移矩阵 row_sums = sum(transition_matrix, 2); row_sums(row_sums == 0) = 1; % 避免除以零 transition_matrix = transition_matrix ./ row_sums; % 构建MTF矩阵 mtf = zeros(n); for i = 1:n for j = 1:n state_i = bin_indices(i); state_j = bin_indices(j); mtf(i, j) = transition_matrix(state_i, state_j); end end end % 保存GAF图像的函数(灰度图像) function saveGAFImage(line_name, gasf, outputDir) % 数据归一化到[0,1]范围 gasf = normalizeImage(gasf); % 创建图像 fig = figure('Visible', 'off'); imagesc(gasf); axis square; colormap(gray); % 使用灰度色图 colorbar; title(sprintf('%s - GAF', line_name), 'FontSize', 10); % 保存图像 outputPath = fullfile(outputDir, sprintf('%s_GAF.png', line_name)); saveas(fig, outputPath); close(fig); fprintf('已保存: %s\n', outputPath); end % 保存MTF图像的函数(灰度图像) function saveMTFImage(line_name, mtf_matrix, time_vector, outputDir) % 数据归一化到[0,1]范围 mtf_matrix = normalizeImage(mtf_matrix); % 创建图像 fig = figure('Visible', 'off'); if ~isempty(time_vector) imagesc(time_vector, time_vector, mtf_matrix); xlabel('Time (s)', 'FontSize', 10); ylabel('Time (s)', 'FontSize', 10); else imagesc(mtf_matrix); xlabel('Time Index', 'FontSize', 10); ylabel('Time Index', 'FontSize', 10); end axis square; colormap(gray); % 使用灰度色图 colorbar; title(sprintf('%s - MTF', line_name), 'FontSize', 10); % 保存图像 outputPath = fullfile(outputDir, sprintf('%s_MTF.png', line_name)); saveas(fig, outputPath); close(fig); fprintf('已保存: %s\n', outputPath); end % 融合GAF和MTF的函数 function fused = fuseGAFandMTF(gasf, mtf, method) % 确保数据在相同范围内 gasf = normalizeImage(gasf); mtf = normalizeImage(mtf); % 应用选择的融合方法 switch lower(method) case 'average' fused = (gasf + mtf) / 2; case 'product' fused = gasf .* mtf; case 'max' fused = max(gasf, mtf); case 'min' fused = min(gasf, mtf); case 'absdiff' fused = abs(gasf - mtf); case 'weighted' % 加权融合 - 给GAF更高权重 fused = 0.7 * gasf + 0.3 * mtf; otherwise error('未知融合方法: %s', method); end % 归一化到[0,1]范围 fused = normalizeImage(fused); end % 图像归一化函数 function normalized = normalizeImage(img) min_val = min(img(:)); max_val = max(img(:)); if max_val == min_val normalized = zeros(size(img)); else normalized = (img - min_val) / (max_val - min_val); end end % 保存融合图像的函数(灰度图像) function saveFusedImage(line_name, fused_image, outputDir) % 创建灰度图像(单通道) gray_image = mat2gray(fused_image); % 确保在[0,1]范围 % 将图像数据转换为uint8类型 gray_image = im2uint8(gray_image); % 保存图像 outputPath = fullfile(outputDir, sprintf('%s_Fused_GAF_MTF.png', line_name)); % 确保文件路径有效 if ~ischar(outputPath) && ~isstring(outputPath) error('文件路径必须为字符向量或字符串标量。'); end imwrite(gray_image, outputPath, 'Quality', 100); fprintf('已保存: %s\n', outputPath); end % 运行主函数 main();请分析这段代码
最新发布
08-08
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值