Vectorization three for loops

Question

I'm new to Matlab, and I need help speeding up some of the code. I am writing a Matlab application that convolves three-dimensional matrices, but unlike standard convolution, the kernel is not constant, it must be calculated for each pixel in the image.

So far I have had working code, but incredibly slow:

function result = calculateFilteredImages(images, T)

% images - matrix [480,360,10] of 10 grayscale images of height=480 and width=360
% reprezented as a value in a range [0..1] 
% i.e. images(10,20,5) = 0.1231;

% T - some matrix [480,360,10, 3,3] of double values, calculated earlier 

    kerN = 5;               %kernel size
    mid=floor(kerN/2);      %half the kernel size
    offset=mid+1;           %kernel offset
    [h,w,n] = size(images);
    %add padding so as not to get IndexOutOfBoundsEx during summation: 
    %[i.e. changes [1 2 3...10] to [0 0 1 2 ... 10 0 0]]
    images = padarray(images,[mid, mid, mid]);

    result(h,w,n)=0;           %preallocate, faster than zeros(h,w,n)
    kernel(kerN,kerN,kerN)=0;  %preallocate

    % the three parameters below are not important in this problem 
    % (are used to calculate sigma in x,y,z direction inside the loop) 
    sigMin=0.5;
    sigMax=3;
    d = 3;

    for a=1:n;
        tic;
        for b=1:w;
            for c=1:h;
                M(:,:)=T(c,b,a,:,:); % M is now a 3x3 matrix
                [R D] = eig(M); %get eigenvectors and eigenvalues - R and D are now 3x3 matrices     

                % eigenvalues
                l1 = D(1,1);
                l2 = D(2,2);
                l3 = D(3,3);

                sig1=sig( l1 , sigMin, sigMax, d);
                sig2=sig( l2 , sigMin, sigMax, d);
                sig3=sig( l3 , sigMin, sigMax, d);

                % calculate kernel
                for i=-mid:mid
                    for j=-mid:mid
                        for k=-mid:mid
                            x_new = [i,j,k] * R; %calculate new [i,j,k]
                            kernel(offset+i, offset+j, offset+k) = exp(- (((x_new(1))^2 )/(sig1^2) + ((x_new(2))^2)/(sig2^2) + ((x_new(3))^2)/(sig3^2)) /2);
                        end
                    end
                end
                % normalize
                kernel=kernel/sum(kernel(:));

                %perform summation
                xm_sum=0;
                for i=-mid:mid
                    for j=-mid:mid
                        for k=-mid:mid
                            xm_sum = xm_sum + kernel(offset+i, offset+j, offset+k) * images(c+mid+i, b+mid+j, a+mid+k);
                        end
                    end
                end
                result(c,b,a)=xm_sum;

            end
        end
        toc;
    end
end

I tried to replace part of the "computational core" with

sigma=[sig1 sig2 sig3]
[x,y,z] = ndgrid(-mid:mid,-mid:mid,-mid:mid);
k2 = arrayfun(@(x, y, z) exp(-(norm([x,y,z]*R./sigma)^2)/2), x,y,z);

but it turned out to be even slower than the cycle. I looked through several articles and vectorization guides, but I am completely stuck with this. Can it be vectorized or somehow accelerated using something else? I'm new to Matlab, maybe there are some built-in functions that could help in this case?

Update Profiling Result:

, :
T.mat
grayImages.mat

Answer 1

, , , . , , ? "" Matlab . * . ^, , , . http://www.mathworks.com/help/matlab/ref/power.html

:

sigma=[sig1 sig2 sig3]
[x,y,z] = ndgrid(-mid:mid,-mid:mid,-mid:mid);
k2 = arrayfun(@(x, y, z) exp(-(norm([x,y,z]*R./sigma)^2)/2), x,y,z);

. , k2.

EDIT: matrix_to_norm x(:) . . ( )

:

% R & mid my test variables
R = [1 2 3; 4 5 6; 7 8 9];
mid = 5;
[x,y,z] = ndgrid(-mid:mid,-mid:mid,-mid:mid);
% meshgrid is also a possibility, check that you are getting the order you want
% Going to break the equation apart for now for clarity
% Matrix operation, should already be fast.
matrix_to_norm = [x(:) y(:) z(:)]*R/sig1
% Ditto
matrix_normed = norm(matrix_to_norm)
% Note the .^ - I believe you want element-by-element exponentiation, this will 
% vectorize it.
k2 = exp(-0.5*(matrix_normed.^2))