filter that uses elements from two arrays at the same time

Question

Let's assume we have two arrays of the same size - A and B.

Now, we need a filter that, for a given mask size, selects elements from A, but removes the central element of the mask, and inserts there corresponding element from B.

So the 3x3 "pseudo mask" will look similar to this:

A A A
A B A
A A A

Doing something like this for averaging filter is quite simple. We can compute the mean value for elements from A without the central element, and then combine it with a proper proportion with elements from B:

h = ones(3,3);
h(2,2) =0; 
h = h/sum(h(:));
A_ave = filter2(h, A);
C = (8/9) * A_ave + (1/9) * B;

But how to do something similar for median filter (medfilt2 or even better for ordfilt2)

Answer 1

The way to solve this is to find a way to combine the information from A and B so that the filtering itself becomes easy.

The first thing I thought of was to catenate A and B along the third dimension and to pass with a filter mask that would take 8 elements from the 'A-slice' and the center element from the 'B-slice'. This is, unfortunately, not supported by Matlab.

While nlfilter only works on 2D images, it does allow you to specify any function for filtering. Thus, you could create a function that somehow is able to look up the right values of A and B. Thus I came to my first solution.

You create a new array, C, that contains the element index at each element, i.e. the first element is 1, the second element is 2, etc. Then, you run nlfilter, which takes a 3x3 sliding window and passes the values of C inside the window to the filtering function, ffn. ffn is an anonymous function, that calls crazyFilter, and that has been initialized so that A and B get passed at each call. CrazyFunction takes the values from the sliding window of C, which are nothing but indices into A and B, and collects the values from A and B from them.

The second solution is exactly the same, except that instead of moving a sliding window, you create a new array that, in every column, has the contents of the sliding window at every possible location. With an overlapping window, the column array gets larger than the original array. Again, you then just need to use the values of the column array, C, which are indices into A and B, to look up the values of A and B at the relevant locations.

EDIT If you have enough memory, im2col and col2im can speed up the process a lot

%# define A,B
A = randn(100);
B = rand(100);

%# pad A, B - you may want to think about how you want to pad
Ap = padarray(A,[1,1]);
Bp = padarray(B,[1,1]);

#% EITHER -- the more more flexible way
%# create a pseudo image that has indices instead of values
C = zeros(size(Ap));
C(:) = 1:numel(Ap);
%# convert to 'column image', where each column represents a block
C = im2col(C,[3,3]);
%# read values from A
data = Ap(C);
%# replace centers with values from B
data(5,:) = Bp(C(5,:));

%# OR -- the more efficient way
%# reshape A directly into windows and fill in B
data = im2col(Ap,[3,3]);
data(5,:) = B(:);

% median and reshape
out = reshape(median(data,1),size(A));

Old version (uses less memory, may need padding)

%# define A,B
A = randn(100);
B = rand(100);

%# define the filter function
ffun = @(x)crazyFilter(x,A,B);

%# create a pseudo image that has indices instead of values
C = zeros(size(A));
C(:) = 1:numel(A);

%# filter
filteredImage = nlfilter(C,[3,3],ffun);




%# filter function
function out = crazyFilter(input,A,B)
%#CRAZYFILTER takes the median of a 3x3 mask defined by input, taking 8 elements from A and 1 from B

%# read data from A
data = A(input(:));
%# replace center element with value from B
data(5) = B(input(5));

%# return the median
out = median(data);

Answer 2

Here's a solution that will work if your data is an unsigned integer type (like a typical grayscale image of type uint8). You can combine your two matrices A and B into a single matrix of a larger integer type, with the data from one matrix stored in the lower bits and the data from the other matrix stored in the higher bits. You can then use NLFILTER to apply a filtering function that extracts the appropriate bits of data in order to collect the necessary matrix values.

The following example applies a median filter of the form you describe above (a 3-by-3 array of elements from A with the center element from B) to two unsigned 8-bit matrices of random values:

%# Initialize some variables:

A = randi([0 255],[3 3],'uint8');  %# One random matrix of uint8 values
B = randi([0 255],[3 3],'uint8');  %# Another random matrix of uint8 values

C = uint16(A)+bitshift(uint16(B),8);  %# Convert to uint16 and place the values
                                      %#   of A in the lowest 8 bits and the
                                      %#   values of B in the highest 8 bits

C = padarray(C,[1 1],'symmetric');  %# Pad the array edges

%# Make the median filtering function for each 3-by-3 block:

medFcn = @(x) median([bitand(x(1:4),255) ...  %# Get the first four A values
                      bitshift(x(5),-8) ...   %# Get the fifth B value
                      bitand(x(6:9),255)]);   %# Get the last four A values

%# Perform the filtering:

D = nlfilter(C,[3 3],medFcn);
D = uint8(D(2:end-1,2:end-1));  %# Remove the padding and convert to uint8

Here are additional links for some of the key functions used above: PADARRAY, BITAND, BITSHIFT.