Vectorizing of interpolation code

1 visualización (últimos 30 días)
Philipp
Philipp el 14 de Jul. de 2014
Comentada: Philipp el 14 de Jul. de 2014
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
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.
Philipp
Philipp el 14 de Jul. de 2014
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.

Iniciar sesión para comentar.

Respuestas (0)

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!

Translated by