Reduce loops to Vectors in Viola Jones Algorithm

For one of my methods for the Viola Jones algorithm I used 4 loops. I noticed that this really hinders performance, so I wanted to attempt to vectorize it. Does anybody know how I would approach this? I'm sorry for the long piece of code.

function featureResponse = evaluate(intImg, feature, stepsize, scale)
%EVALUATE Evaluates the haar-like feature along some step size against the
%integral image
%   intImg: The integral image
%   feature: The HaarFeature object
%   stepsize: length 2 array with y and x step sizes (in that order) along which
%       to position the feature
%   scale: The additional scale (y, y) to apply to HaarFeature.WINDOW_SIZE 
%       and the contained feature rectangles
%   featureResponse: Returns the 2d feature response with one pixel for
%       each placement of the haarFeature along the grid defined by
%       stepsize
 
    import haar.*
    windowSize = HaarFeature.WINDOW_SIZE .* scale;
    [imgSizeY, imgSizeX] = size(intImg);
    dim = num2cell(fix((size(intImg) - windowSize)./stepsize));
    featureResponse = zeros(dim{:});
    %----------------------------------
    % actual implementation here
    %----------------------------------
    
    % PLACEHOLDER:
 
    p = (feature.pos - 1) .* scale + 1; % position of the rectangle within the window
    numSubRects = size(feature.type);
    s = feature.scale .* scale; % size of subrectangle in the scaled window
        
    yIdx = 1;
    y = p(1); % top pivot location of the feature
    
    
    
    while yIdx < dim{1}
        xIdx = 1;
        x = p(2); % left pivot location of the feature
        while xIdx < dim{2}
            result = 0;
            for fy = 1:numSubRects(1)
                for fx = 1:numSubRects(2)
                    result = result + feature.type(fy, fx) *...
                    (...
                            intImg(y  + fy   *s(1), x  + fx   *s(2)) ... % lower right
                        -   intImg(y  +(fy-1)*s(1), x  + fx   *s(2)) ... % upper right
                        -   intImg(y  + fy   *s(1), x  +(fx-1)*s(2)) ... % lower left
                        +  (intImg(y  +(fy-1)*s(1), x  +(fx-1)*s(2)))... % upper left
                    );
                end
            end
             
            featureResponse(yIdx, xIdx) = result;
            x = x + stepsize(2);
            xIdx = xIdx + 1;
        end
 
        y = y + stepsize(1);
        yIdx = yIdx + 1;
    end
end

Answer 1

[dpb]

I don't have the time right now, but the two inner loops should be replaceable by filter2 or conv2. NB: filter2 has an internal rot90 in comparison to conv2 as to array orientation. Confusing, yes.

Also look at blockproc if have Image Processing TB (which I don't)...