Предварительные вычисления весов для многомерной линейной интерполяции

Question

У меня есть неоднородная прямоугольная сетка вдоль измерений D, матрица логических значений V на сетке и матрица точек данных запроса X. Количество точек сетки различается в зависимости от размеров.

Я запускаю интерполяцию несколько раз для одной и той же сетки G и запроса X, но для разных значений V.

Цель состоит в том, чтобы предварительно вычислить индексы и веса для интерполяции и использовать их повторно, потому что они всегда одинаковы.

Вот пример в двух измерениях, в котором мне приходится вычислять индексы и значения каждый раз в цикле, но я хочу вычислить их только один раз перед циклом. Я сохраняю типы данных из своего приложения (в основном одиночные и логические gpuArrays).

% Define grid
G{1} = single([0; 1; 3; 5; 10]);
G{2} = single([15; 17; 18; 20]);

% Steps and edges are reduntant but help make interpolation a bit faster
S{1} = G{1}(2:end)-G{1}(1:end-1);
S{2} = G{2}(2:end)-G{2}(1:end-1);

gpuInf = 1e10;
% It's my workaround for a bug in GPU version of discretize in Matlab R2017a.
% It throws an error if edges contain Inf, realmin, or realmax. Seems fixed in R2017b prerelease.
E{1} = [-gpuInf; G{1}(2:end-1); gpuInf];
E{2} = [-gpuInf; G{2}(2:end-1); gpuInf];

% Generate query points
n = 50; X = gpuArray(single([rand(n,1)*14-2, 14+rand(n,1)*7]));

[G1, G2] = ndgrid(G{1},G{2});

for i = 1 : 4
    % Generate values on grid
    foo = @(x1,x2) (sin(x1+rand) + cos(x2*rand))>0;
    V = gpuArray(foo(G1,G2));

    % Interpolate
    V_interp = interpV(X, V, G, E, S);

    % Plot results
    subplot(2,2,i);
    contourf(G1, G2, V); hold on;
    scatter(X(:,1), X(:,2),50,[ones(n,1), 1-V_interp, 1-V_interp],'filled', 'MarkerEdgeColor','black'); hold off;
end

function y = interpV(X, V, G, E, S)
y = min(1, max(0, interpV_helper(X, 1, 1, 0, [], V, G, E, S) ));
end

function y = interpV_helper(X, dim, weight, curr_y, index, V, G, E, S)
if dim == ndims(V)+1
    M = [1,cumprod(size(V),2)];
    idx = 1 + (index-1)*M(1:end-1)';
    y = curr_y + weight .* single(V(idx));
else
    x = X(:,dim); grid = G{dim}; edges = E{dim}; steps = S{dim};
    iL = single(discretize(x, edges));
    weightL = weight .* (grid(iL+1) - x) ./ steps(iL);
    weightH = weight .* (x - grid(iL)) ./ steps(iL);
    y = interpV_helper(X, dim+1, weightL, curr_y, [index, iL  ], V, G, E, S) +...
        interpV_helper(X, dim+1, weightH, curr_y, [index, iL+1], V, G, E, S);
end
end

Answer 1

Для выполнения задачи должен быть выполнен весь процесс интерполяции, кроме вычисления интерполированных значений. Вот решение, переведенное из источника Octave C++. Формат ввода такой же, как и у первой сигнатуры функции interpn, за исключением того, что нет необходимости v массив. Также Xs должны быть векторами и не должны быть ndgrid формат. Оба выхода W (веса) и I (позиции) имеют размер (a ,b) тот a это число соседей точек на сетке и b количество запрошенных точек для интерполяции.

function [W , I] = lininterpnw(varargin)
% [W I] = lininterpnw(X1,X2,...,Xn,Xq1,Xq2,...,Xqn)
    n     = numel(varargin)/2;
    x     = varargin(1:n);
    y     = varargin(n+1:end);
    sz    = cellfun(@numel,x);
    scale = [1 cumprod(sz(1:end-1))];
    Ni    = numel(y{1});
    index = zeros(n,Ni);
    x_before = zeros(n,Ni);
    x_after = zeros(n,Ni);
    for ii = 1:n
        jj = interp1(x{ii},1:sz(ii),y{ii},'previous');
        index(ii,:) = jj-1;
        x_before(ii,:) = x{ii}(jj);
        x_after(ii,:) = x{ii}(jj+1);
    end
    coef(2:2:2*n,1:Ni) = (vertcat(y{:}) - x_before) ./ (x_after - x_before);
    coef(1:2:end,:)    = 1 - coef(2:2:2*n,:);
    bit = permute(dec2bin(0:2^n-1)=='1', [2,3,1]);
    %I = reshape(1+scale*bsxfun(@plus,index,bit), Ni, []).';  %Octave
    I = reshape(1+sum(bsxfun(@times,scale(:),bsxfun(@plus,index,bit))), Ni, []).';
    W = squeeze(prod(reshape(coef(bsxfun(@plus,(1:2:2*n).',bit),:).',Ni,n,[]),2)).';
end

Тестирование:

x={[1 3 8 9],[2 12 13 17 25]};
v = rand(4,5);
y={[1.5 1.6 1.3 3.5,8.1,8.3],[8.4,13.5,14.4,23,23.9,24.2]};

[W I]=lininterpnw(x{:},y{:});

sum(W.*v(I))
interpn(x{:},v,y{:})

Спасибо @SardarUsama за тестирование и его полезные комментарии.

Answer 2

Я нашел способ сделать это и опубликовать его здесь, потому что (на данный момент) еще два человека заинтересованы. Требуется лишь небольшое изменение моего исходного кода (см. Ниже).

% Define grid
G{1} = single([0; 1; 3; 5; 10]);
G{2} = single([15; 17; 18; 20]);

% Steps and edges are reduntant but help make interpolation a bit faster
S{1} = G{1}(2:end)-G{1}(1:end-1);
S{2} = G{2}(2:end)-G{2}(1:end-1);

gpuInf = 1e10;
% It's my workaround for a bug in GPU version of discretize in Matlab R2017a.
% It throws an error if edges contain Inf, realmin, or realmax. Seems fixed in R2017b prerelease.
E{1} = [-gpuInf; G{1}(2:end-1); gpuInf];
E{2} = [-gpuInf; G{2}(2:end-1); gpuInf];

% Generate query points
n = 50; X = gpuArray(single([rand(n,1)*14-2, 14+rand(n,1)*7]));

[G1, G2] = ndgrid(G{1},G{2});

[W, I] = interpIW(X, G, E, S); % Precompute weights W and indexes I

for i = 1 : 4
    % Generate values on grid
    foo = @(x1,x2) (sin(x1+rand) + cos(x2*rand))>0;
    V = gpuArray(foo(G1,G2));

    % Interpolate
    V_interp = sum(W .* single(V(I)), 2);

    % Plot results
    subplot(2,2,i);
    contourf(G1, G2, V); hold on;
    scatter(X(:,1), X(:,2), 50,[ones(n,1), 1-V_interp, 1-V_interp],'filled', 'MarkerEdgeColor','black'); hold off;
end

function [W, I] = interpIW(X, G, E, S)
global Weights Indexes
Weights=[]; Indexes=[];
interpIW_helper(X, 1, 1, [], G, E, S, []);
W = Weights; I = Indexes;
end

function [] = interpIW_helper(X, dim, weight, index, G, E, S, sizeV)
global Weights Indexes
if dim == size(X,2)+1
    M = [1,cumprod(sizeV,2)];
    Weights = [Weights, weight];
    Indexes = [Indexes, 1 + (index-1)*M(1:end-1)'];
else
    x = X(:,dim); grid = G{dim}; edges = E{dim}; steps = S{dim};
    iL = single(discretize(x, edges));
    weightL = weight .* (grid(iL+1) - x) ./ steps(iL);
    weightH = weight .* (x - grid(iL)) ./ steps(iL);
    interpIW_helper(X, dim+1, weightL, [index, iL  ], G, E, S, [sizeV, size(grid,1)]);
    interpIW_helper(X, dim+1, weightH, [index, iL+1], G, E, S, [sizeV, size(grid,1)]);
end
end