1

这里的这个功能在我的运行中消耗了很多时间。但是看到的是大部分时间都在内置函数polyarea中。可以将此代码矢量化以提高性能吗？

探查器报告 -

  time   calls
                  1 function [S S_area] = Polygons_intersection_Compute_area(S)
                  2 % Guillaume JACQUENOT
                  3 % guillaume at jacquenot at gmail dot com
                  4 % 2007_10_08
                  5 % 2009_06_16
                  6 % Compute area of each polygon of in S.
                  7 % Results are stored as a field in S
                  8 
  0.50   51945    9 S_area = struct('A', {}); 
  0.20   51945   10 for i=1:numel(S) 
  0.28  103890   11     S(i).area = 0; 
  1.34  103890   12     S_area(i).A = zeros(1,numel(S(i).P)); 
  0.69  103890   13     for j=1:numel(S(i).P) 
  9.24  103890   14         S_area(i).A(j) = polyarea(S(i).P(j).x,S(i).P(j).y); 
  0.28  103890   15         S(i).area      = S(i).area + (1-2*S(i).P(j).hole) * S_area(i).A(j);         
  0.01  103890   16     end 
  0.08  103890   17 end
4

1 回答 1

5

我看到4个问题。我将按潜在性能增益的递增顺序讨论它们。

第一：你使用iandj作为循环变量名。这些也是 Matlab 中虚数单位的名称，这意味着 Matlab 将不得不花一些时间来查找您所指的虚数单位。问题是，如果循环不是 JIT的，它必须在每次迭代中都这样做（你的不是，我会讲到的）。

第二：索引多维结构需要的时间比你想象的要多。多维结构在这方面有点臭名昭著，最好避免对它们进行过多的索引操作。通常制作一个元素的简单副本，对该副本执行所有操作，然后将副本写回结构可以大大提高性能。

第三：您没有S_area以最有效的方式预先分配。您甚至不预先分配结构，而是在您预先分配时在第一个循环中增长它S_area(i).A。这一切都可以改进（见下文）。

第四：polyarea不是内置函数，所以这个双循环不会被 JIT 处理。如果您在您或 Mathworks 用 M 语言（而不是 C）编写的循环内调用任何函数，JIT 编译器将无法编译您的循环。这是迄今为止JIT 框架中最烦人的（也是可改进的）限制，而 JIT 循环的运行速度通常比非 JIT 循环快 100 倍或更多。

唯一的解决方案通常是在循环体中“内联”一个非内置函数。在 Matlab 中，这意味着：将函数体的全部内容复制粘贴到循环中，并对在该函数体中调用的所有非内置函数递归地执行此操作。

以上所有导致此版本的代码：

% pre-allocate S_area
S_area(numel(S)).A = [];
As = cellfun(@(x) zeros(numel(x),1), {S.P}, 'UniformOutput', false);
[S_area.A] = deal(As{:});

% number of polygons for all S(ii)
numPolys = cellfun(@numel, {S.P});

% enter loop
for ii = 1:numel(S)
    % extract S(ii) only once
    Sii = S(ii);

    Sii.area = 0;
    Aii = S_area(ii).A;        
    for jj = 1:numPolys(ii)

        p = Sii.P(jj);  % extract polygon only once
        x = p.x; % and its x and y components
        y = p.y;            
        sz = size(p);

        % NOTE: core of polyarea. Note that all checks and flexibility, and 
        % therefore user-friendliness, is GONE. Very little has to go wrong 
        % here before a hard-to-understand error is issued. 
        Area = reshape(abs(sum( (x([2:sz(1) 1],:) - x(:,:)).* ...
            (y([2:sz(1) 1],:) + y(:,:)))/2),[1 sz(2:end)]);

        Aii(jj) = Area;
        Sii.area = Sii.area + Area*(1-2*p.hole);
    end

    % place copies back into the strucure
    S_area(ii).A = Aii;
    S(ii).area = Sii.area;

end

我无法尽可能正确地对此进行测试，因此如果您发现一些错误，请告诉我，我会尽力纠正它们。

于 2012-10-30T10:58:26.450 回答