Vectorizing of interpolation code

1 次查看（过去 30 天）
Philipp
Philipp 2014-7-14
评论： Philipp 2014-7-14
Hello Matlab Community,
I have written the following interpolation class that performs a bicubic interpolation. I would very appreciate if somebody could give me some concrete hints how to avoid the particular 'for loop' in the 'interpolate method'
% classdef BiCubicInterpolation
%UNTITLED2 Summary of this class goes here
% Detailed explanation goes here
properties (Constant)
%1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
PolynomCoefficient = ...
[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 %1
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 %2
-3 3 0 0 -2 -1 0 0 0 0 0 0 0 0 0 0 %3
2 -2 0 0 1 1 0 0 0 0 0 0 0 0 0 0 %4
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 %5
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 %6
0 0 0 0 0 0 0 0 -3 3 0 0 -2 -1 0 0 %7
0 0 0 0 0 0 0 0 2 -2 0 0 1 1 0 0 %8
-3 0 3 0 0 0 0 0 -2 0 -1 0 0 0 0 0 %9
0 0 0 0 -3 0 3 0 0 0 0 0 -2 0 -1 0 %10
9 -9 -9 9 6 3 -6 -3 6 -6 3 -3 4 2 2 1 %11
-6 6 6 -6 -3 -3 3 3 -4 4 -2 2 -2 -2 -1 -1 %12
2 0 -2 0 0 0 0 0 1 0 1 0 0 0 0 0
0 0 0 0 2 0 -2 0 0 0 0 0 1 0 1 0
-6 6 6 -6 -4 -2 4 2 -3 3 -3 3 -2 -1 -2 -1
4 -4 -4 4 2 2 -2 -2 2 -2 2 -2 1 1 1 1];
PolynomExponent = ...
[0 1 2 3
0 1 2 3
0 1 2 3
0 1 2 3];
end
properties (SetAccess = private)
VecX
VecY
nx
ny
dx
dy
F
Fdx
Fdy
Fdxdy
end
methods (Access = private)
function Filter = CentralDiffOperatorKernel(h)
Filter = 1/2/h*[0 0 0
-1 0 1
0 0 0];
end
end
methods (Access = public)
function obj = BiCubicInterpolation(X,Y,FF)
% Constructor
if(nargin == 3)
tic
obj.VecX = X;
obj.VecY = Y;
obj.nx = numel(obj.VecX);
obj.ny = numel(obj.VecY);
obj.F = [FF zeros(obj.nx,1) ; zeros(1,obj.ny+1)];
obj.dx = obj.VecX(2) - obj.VecX(1);
obj.dy = obj.VecY(2) - obj.VecY(1);
% Calculation of derivatives and crossderivatives
obj.Fdx = conv2(obj.F, CentralDiffOperatorKernel(obj.dx), 'same');
obj.Fdy = conv2(obj.F, CentralDiffOperatorKernel(obj.dy)', 'same');
obj.Fdxdy = conv2(obj.Fdx, CentralDiffOperatorKernel(obj.dy)', 'same');
toc
disp('Constructor')
end
end
function Int = Interpolate(obj, P)
% Evaluate index
xi = round((P(:,1) - mod(P(:,1), obj.dx)).*10)./10+1;
yi = round((P(:,2) - mod(P(:,2), obj.dy)).*10)./10+1;
% Normalize x[0 1] y[0 1]
x = mod(P(:,1), obj.dx)./obj.dx;
y = mod(P(:,2), obj.dy)./obj.dy;
ll = numel(P(:,1));
Int = zeros(ll,1);
xi(xi < 1) = 1;
yi(yi < 1) = 1;
xi(xi > obj.nx) = obj.nx;
yi(yi > obj.ny) = obj.ny;
for kk = 1:1:ll
if(isnan(xi(kk))||isnan(yi(kk)) )
yy = nan(16,1);
else
yy = [obj.F(xi(kk) ,yi(kk)) obj.F(xi(kk) + 1 ,yi(kk)) obj.F(xi(kk),yi(kk) + 1) obj.F(xi(kk) + 1,yi(kk) + 1)...
obj.Fdx(xi(kk) ,yi(kk)) obj.Fdx(xi(kk) + 1 ,yi(kk)) obj.Fdx(xi(kk),yi(kk) + 1) obj.Fdx(xi(kk) + 1,yi(kk) + 1)...
obj.Fdy(xi(kk) ,yi(kk)) obj.Fdy(xi(kk) + 1 ,yi(kk)) obj.Fdy(xi(kk),yi(kk) + 1) obj.Fdy(xi(kk) + 1,yi(kk) + 1)...
obj.Fdxdy(xi(kk) ,yi(kk)) obj.Fdxdy(xi(kk) + 1 ,yi(kk)) obj.Fdxdy(xi(kk),yi(kk) + 1) obj.Fdxdy(xi(kk) + 1,yi(kk) + 1) ]';
end
ap = obj.PolynomCoefficient *yy;
ap = (reshape(ap,4,4));
Int(kk) = sum(sum(ap.*x(kk).^(obj.PolynomExponent)'.*y(kk).^(obj.PolynomExponent)));
end
Int((P(:,1)>(obj.VecX(end)))|(P(:,1)<(obj.VecX(1)))|(P(:,2)>(obj.VecY(end)))|(P(:,2)<(obj.VecY(1)))) = nan;
end
end
end
Any help would be great.
:)
  7 个评论
显示 5更早的评论
John D'Errico
John D'Errico 2014-7-14
Are you doing 100000 calls with interp2? It is already vectorized, so ONE call would suffice.
There will be some loss of speed with interp2, since it checks for errors, and it allows a general non-uniform grid spacing.
Philipp
Philipp 2014-7-14
Missunderstanding. Evey call of the 100000 calls is requesting an array of points. The lentgth of this array goes from 10 to 1000. The spacing of the points is non uniformly.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Interpolating Gridded Data 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by