Vectorizing of interpolation code

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.

:)

[Philipp]

Maybe wirting the class/function in c is more performant?

[Oleg Komarov]

Please, provide sample inputs for X, Y and FF.

[Philipp]

ok, here we go.

 x  = 0:1:100;
 y = 0:1:300;
 [xx,yy] = ndgrid(x,y);
 F = xx/10.*exp(xx/10.^2 - yy/10.^2);

 Int = BiCubicInterpolation(x,y,F);

 xi = 1:.1:100;
 yi = xi;

 P = [xi' yi'];
 Fi = Int.Interpolate(P);

[John D'Errico]

Is there an intelligent reason why you are not simply using interp2?

[Philipp]

Yes - at least I think so. interp2 is less performant when want to interpolate only a small set of data. (no compact grid - just a list of points) So in my particular case, I do maybe 100000 calls of *.Interpolate(P), and a lot of them are requesting several points to interpolate. So overall performace with interp2 is significant lower than using my interpolation (If i made no faults). So I decided to gain speed by vectorizing this method.

[John D'Errico]

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]

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.