Maybe wirting the class/function in c is more performant?
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
Respuestas (0)
Ver también
Categorías
Más información sobre Interpolation 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!