colfilt() with large matrices

On my Matlab 2014a installation, colfilt() behaves rather strange with large matrices. Consider the following code:

% Define basic parameters
matSize = 160;     % must divide by 20
kernelSize = 33;

% Create checkerboard
inputMat = repmat([zeros(matSize/20), ones(matSize/20); ones(matSize/20), zeros(matSize/20)], 10, 10);

% Define "convolution" kernel which just returns the original pixel
kernel = [zeros(ceil(kernelSize^2/2 - 1), 1); 1; zeros(floor(kernelSize^2/2), 1)];

% Convolve with this kernel
outputMat = colfilt(inputMat, [kernelSize, kernelSize], 'sliding', @(x) colfilt_func(x, kernel));

% Show image
imshow(outputMat);

This code is just for demonstration purposes. It essentially creates a checkerboard-type image, defines a convolution kernel and then uses colfilt to do the convolution. The convolution kernel is defined such that it just reproduces the original image.

The code invokes a filter function which is defined as follows:

function [ result ] = colfilt_func( x, kernel )
    % Output size of x to console
    disp(['Input from colfilt: ', num2str(size(x, 1)), ' x ', num2str(size(x, 2))]);

      % Convolve
      result = sum(x .* repmat(kernel, 1, size(x, 2)));
  end

The filter function outputs some internal information to the command window so that we can monitor what's happening.

When setting the "matSize" variable to 80, we should expect 80^2 = 6400 convolutions. This is exactly what happens, hence the following output is produced:

Input from colfilt: 1089 x 6400.

By contrast, setting "matSize" to 160 should yield 16^2 = 25600 convolutions. The matrix which needs to be passed from colfilt() to the filter function is so huge that Matlab splits it into chunks in order to save memory. This is mentioned in the documentation. However, chunking is not the only thing that happens. Instead, I get this output:

Input from colfilt: 1089 x 1.

Input from colfilt: 1089 x 6400.

Input from colfilt: 1089 x 6400.

Input from colfilt: 1089 x 6400.

Input from colfilt: 1089 x 6400.

Notice the first invocation of the filter function by colfilt(). It passes a 1089x1 column vector. This vector contains the first 33x33 entries of the input matrix (the size of the filter kernel), not zero-padded. Only thereafter, the actual convolution starts and I get four invocations of the filter function with equally sized chunks. These matrices are correctly zero-padded and work perfectly fine.

My question: what is it with the first invocation? Why does this happen?

[Walter Roberson]

Sometimes, in order to pre-allocate output properly for something that returns a size that is not obvious, it is easiest to run one iteration and use the size returned as the basic template for preallocating for the computation as a whole.

Answer 1

[Blackadder]

I have now looked into the code of colfilt() and figured out the answer to my question. The crucial information is already present in the documentation. It states that colfilt() may split a convolution operation into multiple chunks to save memory. It calls the bestblk() function to decide whether and how to split the input matrix.

While this is a very reasonable thing to do, it requires an additional step during filtering. If an input matrix can be processed in a single step, the matrix is first rearranged by colfit(), then fed into the filter function and the filter function returns the complete filtered image. In this case, colfilt() never has to look at the result from the filter function but just pipe it through as its own return value.

If, however, the input matrix needs to be split into chunks, colfilt() must collect all matrices that are returned by the filter function in order to combine them into the complete filtered image later. This requires that colfilt() must know the data type of the output from the filter function, e.g., single, double, uint8, and so forth.

That's precisely what the first run with only one column is for, as suggested by Walter Roberson. It is needed for colfilt2() to determine the data type returned by the filter function. colfilt() then preallocates its return matrix for the filtered image according to the data type it just received. Only now can it start with the actual filtering.

Answer 2

[Blackadder]

Also look at my response here: https://www.mathworks.com/matlabcentral/answers/234750#answer_238217