我有一个沿D维度的非均匀矩形网格,一个网格上的逻辑值V矩阵,以及一个查询数据点X的矩阵 . 网格点的数量在不同维度上有所不同 .
我对同一网格G和查询X多次运行插值,但对于不同的值V.
目标是预先计算插值的索引和权重并重用它们,因为它们总是相同的 .
这是一个2维的例子,我必须在循环中每次计算索引和值,但我想在循环之前只计算一次 . 我保留了我的应用程序中的数据类型(主要是单个和逻辑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
2 回答
我找到了一种方法来做这个并在这里发布,因为(截至目前)还有两个人感兴趣 . 我只需稍微修改一下原始代码(见下文) .
为了完成任务,应该完成插值的整个过程,除了计算插值 . 这是从Octave c++ source翻译的解决方案 . 输入的格式与interpn函数的第一个签名相同,只是不需要
v
数组 . 此外X
应该是矢量,不应该是ndgrid
格式 . 输出W
(权重)和I
(位置)的大小均为(a ,b)
,其中a
是网格上某点的邻居数,b
是要插入的请求点数 .测试:
感谢@SardarUsama的测试和他的有用评论 .