MATLB 优化多维矩阵数据的二维区间划分

问题描述 投票:0回答:1

有一个四维矩阵数据,以及每个元素对应的X坐标data_X和Y坐标data_Y。我想根据edge_X和edge_Y进行二维区间分割,或者说网格划分,得到一个新的数组data_class。当前的循环执行速度非常慢。如何提高计算速度?

clc;clear
rng("default")
data = randi([0, 100], 40,30,20, 10);
data_X =  randi([0, 100], 40,30,20, 10);
data_Y =  randi([0, 100], 40,30,20, 10);

% Define the edges for histograms in X and Y directions
edge_X = [0:10:100];
edge_Y = [0:20:100];

sz = size(data);
for ix = 1:sz(1)
for iy = 1:sz(2)
    for ii = 1:sz(3)
        for ij = 1:sz(4)
           % Extract the current data, X, and Y values
            g_data = data(ix, iy, ii, ij);
            g_X = data_X(ix, iy, ii, ij);
            g_Y = data_Y(ix, iy, ii, ij);

            % Compute the 2D histogram count  
            g_id = histcounts2(g_X, g_Y, edge_X, edge_Y);
            g_new = g_id * g_data;
            data_class  (ix,iy,ii,ij,:,:) = g_new; 
        end
    end
end
end

由于data的每个元素都是单独处理的,所以g_data在一个bin中只会被记为1,然后

g_new = g_id * g_data; 

g_new是重分类后g_data的位置

当前循环执行速度非常慢。如何提高计算速度?

matlab histogram
1个回答
0
投票

您应该做的第一件事是查看编辑器中关于

data_class
随着每次循环迭代而增长的橙色下划线警告。在所有循环之前初始化它可以立即节省大部分运行时间

在循环之前预分配

data_class

szout = [sz, [numel(edge_X),numel(edge_Y)]-1];
data_class = NaN(szout);

但主要瓶颈是对每个值单独调用

histcounts2
,您可以通过将其置于
ii/ij
循环之外来节省大量时间,而不是使用第一个输入作为其所在容器的
1
计数在中,使用第四个和第五个输出来存储每个值应进入的容器。在您的情况下,这是行索引和列索引。

以下是一些时间安排:

  • 基线代码:20.53秒
  • 预分配
    data_class
    7.75秒(-62%总运行时间)
  • 解除嵌套
    histcounts2
    0.53秒(-97%总运行时间)
  • 仅分配非零值:0.10秒(-99.5%总运行时间)

解巢

histcounts2

sz = size(data);
szout = [sz, [numel(edge_X),numel(edge_Y)]-1];
data_class = NaN(szout);
g_new = zeros(szout(end-1:end));
for ix = 1:sz(1)
    for iy = 1:sz(2)

        g_X = data_X(ix, iy, :, :);
        g_Y = data_Y(ix, iy, :, :);
        [~,~,~,bx,by] = histcounts2(g_X, g_Y, edge_X, edge_Y);
        bx = squeeze(bx);
        by = squeeze(by);
        g_data = squeeze( data(ix, iy, :, :) );
            
        for ii = 1:sz(3)
            for ij = 1:sz(4)
                g_new(:) = 0;        
                g_new(bx(ii,ij),by(ii,ij)) = g_data(ii,ij);
                data_class(ix,iy,ii,ij,:,:) = g_new; 
            end
        end
    end
end

此外,您不需要仅为每个切片初始化零数组

g_new
来设置一个值。相反,将
data_class
初始化为零矩阵,然后覆盖您感兴趣的特定值:

仅分配非零值

sz = size(data);
szout = [sz, [numel(edge_X),numel(edge_Y)]-1];
data_class = zeros(szout);
for ix = 1:sz(1)
    for iy = 1:sz(2)

        g_X = data_X(ix, iy, :, :);
        g_Y = data_Y(ix, iy, :, :);
        [~,~,~,bx,by] = histcounts2(g_X, g_Y, edge_X, edge_Y);
        bx = squeeze(bx);
        by = squeeze(by);
        g_data = squeeze( data(ix, iy, :, :) );
            
        for ii = 1:sz(3)
            for ij = 1:sz(4)
                data_class(ix,iy,ii,ij,bx(ii,ij),by(ii,ij)) = g_data(ii,ij); 
            end
        end
    end
end

您可以使用

sub2ind
做一些事情来完全摆脱
ii/ij
循环,但在快速测试中,与仅使用循环相比,没有显示出任何性能。

© www.soinside.com 2019 - 2024. All rights reserved.