How to apply element-wise operations without using for loops and without influencing the speed

Question

Suppose that I have an RGB image matrix and I want to apply some spatial filters on it. In general I want to apply element-wise operations (note that it's a college assignment and I'm not permitted to use any built-in functions available in the Image Processing toolbox). I decided to write the filters as functions and then apply bsxfun to these functions on the image.

A simple example would be this:
I want to add 50 to all gray levels of an image and then replace all gray levels with above 200 with 200. Here's my code:

a='C:\Users\sepideh\Desktop\IP_abadpour\S45C-113050518040.jpg';  
b=imread(a);  
b(:,:,1)=b(:,:,1)+50;  
b(:,:,2)=b(:,:,2)+50;  
b(:,:,3)=b(:,:,3)+50;  
c=reshape(b,[],1);  
d=bsxfun(@test,c,200);

test is a function in this form:

function Out = test(in,a)  
   if in>a
      in=200;
   end
   Out = in;
end  

This code won't work because in the second line "in > a" is a matrix having 0's and 1's (I mean all of the elements are not 1 and should not be) so the debugger won't branch into the if statement.

Could you guide me how to write this function and how to apply spatial and fourier analyses on the image, without affecting the performance and run-time speed?

Answer 1

Here's a couple of suggestions:

1. First of all, you don't need to add 50 to each layer of the RGB matrix individually. You can just do:

   b = b + 50;
   

2. Why do you reshape b before passing it to bsxfun? The size of the output of bsxfun is the same as your image's, there's really no need in reshaping here anything.

3. Regarding your test function, note what the official documentation of bsxfun states:

   A binary element-wise function of the form C = fun(A,B) accepts arrays A and B of arbitrary but equal size and returns output of the same size. Each element in the output array C is the result of an operation on the corresponding elements of A and B only. fun must also support scalar expansion, such that if A or B is a scalar, C is the result of applying the scalar to every element in the other input array.

   So bsxfun performs singleton expansion and "inflates" its two input arrays to the same size, and then applies the specified function to the inflated arrays. The element-wise function fun operates, in fact, on the arrays, not scalars. I don't see any actual gain in employing bsxfun here.

   That said, you can simplify your code as shown in Dan's suggestion, or implement it as a function:

   function out = test(in, a);
       out = in;
       out(in > a) = a;
   

   I assume that if you were using the value 210 instead of 200, you'd like to cap all gray levels with 210 as well, so you should really be using a instead of a hard-coded value 200. You could also write your function like so:

   function out = test(in, a)
       out = min(in, a);
   

   and then invoke it with:

   d = test(b, 200);
   

   instead of the more complicated d = bsxfun(@test, b, 200).

   Another alternative is to use arrayfun:

   d = arrayfun(@(x)test(x, 200), a);
   

   or

   d = arrayfun(@test, a, 200 * ones(size(a)));
   

   in which arrayfun will apply test element-wise, and the test function would need to operate only on scalars. However, arrayfun usually runs slower than loops, let alone vectorized operations.

4. For spatial analysis, check out conv2 just like Dan suggested (or implement your own 2-D convolution, for the sake of practice). For Fourier analysis, consider using the fft2 and ifft2 functions in the frequency domain.

Hope this helps!

Answer 2

So for the example you posted you can just take advantage of the fact that most operators in matlab work on matrices natively:

b=imread(a);
c = a + 50;
c(c > 200) = 200;

It's as simple as that.

For the filtering, if you are allowed, I would have a look at the conv2 function. You can do spatial filtering this way without transforming to the frequency domain (remeber, multiplication of a filter in the frequency domain is the same as convolution in the spatial domain). So for example a basic low pass filter:

lpf = ones(5)./25;
c(:,:,1) = conv2(b(:,:,1), lpf);
c(:,:,2) = conv2(b(:,:,2), lpf);
c(:,:,3) = conv2(b(:,:,3), lpf);