다차원 선형 보간을위한 가중치 사전 계산

Question

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 

Answer 1

나는 (지금 현재까지) 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 

Answer 2

작업을 수행하려면 보간 된 값을 계산하는 것 외에는 전체 보간 프로세스를 수행해야합니다. 다음은 Octave c++ source에서 번역 된 솔루션입니다. 입력 형식은 v 배열을 필요로하지 않는다는 점을 제외하면 interpn 함수의 frst 시그니처와 동일합니다. 또한 X은 벡터이어야하며 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,2,[]),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 덕분에 그의 유용한 의견.