多维线性插值的预计算权重 [英] Precompute weights for multidimensional linear interpolation

问题描述

我在D维度上有一个不均匀的矩形网格，在网格上有一个逻辑值V的矩阵，以及一个查询数据点X的矩阵.网格点的数量在各个维度上都不同.

I have a non-uniform rectangular grid along D dimensions, a matrix of logical values V on the grid, and a matrix of query data points X. The number of grid points differs across dimensions.

对于相同的网格G和查询X，但对于不同的值V，我多次运行插值.

I run the interpolation multiple times for the same grid G and query X, but for different values V.

目标是为插值预先计算索引和权重并重新使用它们，因为它们始终相同.

The goal is to precompute the indexes and weights for the interpolation and to reuse them, because they are always the same.

这里是一个二维示例，其中我每次必须在循环中计算索引和值，但是我只想在循环之前计算一次.我保留了应用程序中的数据类型(主要是单个和逻辑gpuArrays).

Here is an example in 2 dimensions, in which I have to compute indexes and values every time within the loop, but I want to compute them only once before the loop. I keep the data types from my application (mostly single and logical 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

推荐答案

要完成此任务，除了计算插值，还要完成插值的整个过程.这是从 Octave c ++源代码翻译而来的解决方案.输入格式与 interpn 的第一签名相同函数，除了不需要v数组.同样，X应该是向量，并且不应该是ndgrid格式.输出W(权重)和I(位置)都具有(a ,b)的大小，其中a是网格上点的邻居数，而b是要内插的请求点数

To do the task the whole process of interpolation ,except computing the interpolated values, should be done. Here is a solution translated from the Octave c++ source. Format of the input is the same as the frst signature of the interpn function except that there is no need to the v array. Also Xs should be vectors and should not be of the ndgrid format. Both the outputs W (weights) and I (positions) have the size (a ,b) that a is the number of neighbors of a points on the grid and b is the number of requested points to be interpolated.

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的测试和他的有用评论.

Thanks to @SardarUsama for testing and his useful comments.

这篇关于多维线性插值的预计算权重的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！