2

在我现在正在处理的一个问题中，我计算矩阵中的一些值，x然后对于x需要的每个元素，在单调递增的向量中找到下面最近元素的索引X以及x元素与他们两边的第一个元素。（这本质上是线性插值，没有进行实际插值。）我经常这样做，所以我真的对它尽可能快地超级感兴趣。

我写了一个函数locate，我可以用一些示例数据调用它：

X = linspace(5, 300, 40)';
x = randi(310, 5, 6, 7);

[ii, weights] = locate(x, X);

我写了两个版本的locate. 第一个是为了说明，第二个是我加快计算速度的最佳尝试。对于如何进一步提高性能，您有什么建议或替代方法吗？

1. 博览会

function [ii, weights] = locate(x, X)
    % LOCATE Locate first node on grid below a given value.
    %
    %   [ii, weights] = locate(x, X) returns the first node in X that is below
    %   each element in x and the relative proximities to the two closest nodes.
    %
    %   X must be a monotonically increasing vector. x is a matrix (of any
    %   order).

    % Preallocate
    ii = ones(size(x));  % Indices of first node below (or 1 if no nodes below)
    weights = zeros([2, size(x)]);  % Relative proximity of the two closest nodes

    % Find indices and compute weights
    for ix = 1:numel(x)
        if x(ix) <= X(1)
            ii(ix) = 1;
            weights(:, ix) = [1; 0];
        elseif x(ix) >= X(end)
            ii(ix) = length(X) - 1;
            weights(:, ix) = [0; 1];
        else
            ii(ix) = find(X <= x(ix), 1, 'last');
            weights(:, ix) = ...
                [X(ii(ix) + 1) - x(ix); x(ix) - X(ii(ix))] / (X(ii(ix) + 1) - X(ii(ix)));
        end
    end
end

2. 最佳尝试

function [ii, weights] = locate(x, X)
    % LOCATE Locate first node on grid below a given value.
    %
    %   [ii, weights] = locate(x, X) returns the first node in X that is below
    %   each element in x and the relative proximities to the two closest nodes.
    %
    %   X must be a monotonically increasing vector. x is a matrix (of any
    %   order).

    % Preallocate
    ii = ones(size(x));  % Indices of first node below (or 1 if no nodes below)
    weights = zeros([2, size(x)]);  % Relative proximity of the two closest nodes

    % Find indices
    for iX = 1:length(X) - 1
        ii(X(iX) <= x) = iX;
    end

    % Find weights
    below = x <= X(1);
    weights(1, below) = 1;  % All mass on the first node
    weights(2, below) = 0;

    above = x >= X(end);
    weights(1, above) = 0;
    weights(2, above) = 1;  % All mass on the last node

    interior = ~below & ~above;
    xInterior = x(interior)';
    iiInterior = ii(interior);
    XBelow = X(iiInterior)';
    XAbove = X(iiInterior + 1)';
    weights(:, interior) = ...
        [XAbove - xInterior; xInterior - XBelow] ./ (XAbove - XBelow);
end
4

1 回答 1

1

polylineinterp在Brain2Mesh 工具箱中检查我的功能。

https://github.com/fangq/brain2mesh/blob/master/polylineinterp.m

几乎完全一样，除了polylen输入就像你的 X 的差异。

一般来说，矢量化这种操作是使用histc()，比如这条线

https://github.com/fangq/brain2mesh/blob/master/polylineinterp.m#L52

于 2020-04-04T13:44:56.203 回答