Matlab: Faster finding of 1D linear interpolation nodes and weights for each element in ND matrix

Question

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.

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);

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. Exposition

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. Best attempt

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

Answer 1

checkout my polylineinterp function in the Brain2Mesh toolbox.

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

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

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

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