Vectorizing of interpolation code

Question

0 个投票

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更早的评论隐藏 5更早的评论

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 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.

请先登录，再进行评论。

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

Follow Question