Vectorizing of interpolation code
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
John D'Errico
el 14 de Jul. de 2014
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.
Respuestas (0)
Ver también
Categorías
Más información sobre Interpolating Gridded Data en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!