Matlab:更快找到ND矩阵中每个元素的一维线性插值节点和权重 [英] Matlab: Faster finding of 1D linear interpolation nodes and weights for each element in ND matrix

问题描述

在我正在处理的问题中，我在矩阵x中计算一些值，然后对于x中的每个元素，需要在单调递增的向量以及x元素到其任一侧的第一个元素的相对距离. (这实际上是线性插值，没有进行实际的插值.)我一直在进行此操作，因此我对使它尽可能快地表现出了极大的兴趣.

In a problem I'm working on now, I compute some values in a matrix x and I then for each element in x need to find the index of the closest element below in a monotonically increasing vector X as well as the relative proximity of the x elements to the first elements on their either side. (This is essentially linear interpolation without doing the actual interpolation.) I'm doing this maaaany times so I really super extra interested in it being as fast as possible.

我已经编写了一个函数locate，可以使用一些示例数据进行调用:

I have written a function locate that I can call with some example data:

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

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

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

I have written two versions of locate. The first is for exposition and the second is my best attempt at speeding up the computations. Do you have any suggestions or alternative approaches for how I could accelerate performance further?

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

推荐答案

在 Brain2Mesh工具箱.

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

几乎可以做到这一点，除了polylen输入就像X的差异一样.

does almost exactly that, except the polylen input is like the diff of your X.

通常，矢量化这种操作的方法是使用histc()，就像这行

in general, vectorizing this kind of operation is to use histc(), like this line

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

这篇关于Matlab:更快找到ND矩阵中每个元素的一维线性插值节点和权重的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！