Improve vectorization without using for loops

Hello,

I'm going crazy over a problem that I have, and it's just about improving performance in a part of a code I have here: the goal is to calculate the velocity vector for each pixel in two consecutive images.

The basic code I wrote is this one: (hope is understandable, if not ask me)

    % im1, im2 - two consecutive images with the same size
    % N - pixels neighbourhood (2*N+1) x (2*N+1)
    function [ L ] = my_ctOpticFlux( im1, im2, N )
    [height, width] = size(im1);

    % Spatiotemporal gradients of the images
    % Ix, Iy, It have the same size as im1, im2
    [Ix, Iy, It] = my_gradients_xyt(im1, im2);

    % List L with
    % [x y vx vy] in each line
    % (x,y) pixel location; (vx,vy) velocity components in (x,y)
    L = zeros((width-2*N)*(height-2*N), 4);

    n = 1;
    % For each image pixel, (Px,Py), ... (no border pixels)
    for Px = 1+N:height-N
        for Py = 1+N:width-N
  % determine the neighbourhood of the current pixel in (x,y)
  % and take the correspondent values of the gradients
            xlims = Px-N:Px+N; ylims = Py-N:Py+N;
            ix = Ix(xlims, ylims);
            iy = Iy(xlims, ylims);
            it = It(xlims, ylims);

            % ... obtain the matrixs A,b for that pixel neighbourhood;
            A = [ix(:), iy(:)];
            b = -it(:);

            % ... compute velocities vector v for that neighbourhood center;
            v = (pinv(A)*b)';

            % Save results on L
            L(n,:) = [Px Py v(1) v(2)];
            n = n + 1;
        end
    end

A little modification:

    % Calculate the limits of computable image on x and y
    regionlims = [(1+N) (height-N) (1+N) (width-N)];

    [x, y] = meshgrid(regionlims(1):regionlims(2),regionlims(3):regionlims(4));
    x = x(:); y = y(:);

Any help on this? I tried to use meshgrids, but no luck.

EDIT1: Some typo in the code.

Answer 1

[Jan]

Indexing is faster, when you insert the generator of the vector:

            ix = Ix(Px-N:Px+N, Py-N:Py+N, t);
            iy = Iy(Px-N:Px+N, Py-N:Py+N, t);
            it = It(Px-N:Px+N, Py-N:Py+N, t);

See https://www.mathworks.com/matlabcentral/answers/35676-why-not-use-square-brackets#answer_252096 .

Now search the bottleneck of the code using the profiler. I guess it is pinv, but it might be my_gradients_xyt also. In the first case, I do not see a chance to accelerate this beside parfor, in the 2nd case posting the code of my_gradients_xyt would help.

[Joel Möllering]

Even doing the same thing three times? Instead of doing one time and use the result 3x later?