2

This is a bit more of a general question, but, no matter how many times I read the description of MATLAB's im2col function, I cannot fully understand it. I need it for the computational efficiency because MATLAB is awful with nested for loops. Here's what I'm attempting to do, but using nested for loops:

 [TRIMMED]=TM_FILTER(IMAGE, FILTER_SIZE, PERCENT)
    Takes a 2-D array and returns the array, filtered with a
    square trimed mean filter with length/width equal to FILTER_SIZE and percent equal to PERCENT.

    %}
    function [trimmed]=tm_filter(image, filter_size, percent)
    if rem(filter_size, 2)==0                            %make sure filter has a center pixel
        error('filter size must be odd numbered');       %error and return if number is odd
        return 
    end
    if percent > 100 || percent < 0
        error('Percentage must be ? [0, 100]');
        return
    end

    [rows, columns]=size(image);              %figure out pixels needed
    n=(filter_size-1)/2;                      %n is pixel distance from center pixel to boundaries  
    padded=(padarray(image, [n,n],128));      %padding on boundaries so center pixel always has neighborhood

for i=1+n:rows                            %rows from first non-padded entry to last nonpadded entry
    for j=1+n:columns                     %colums from first non-padded entry to last nonpadded entry
    subimage=padded(i-n:i+n,j-n:j+n);     %neighborhood same size as filter
    average=trimmean(trimmean(subimage, percent), percent);         %computes trimmed mean of neighborhood as trimmed mean of vector of trimmed means
    trimmed(i-n, j-n)=average;         %stores averaged pixel in new array
    end
end
trimmed=uint8(trimmed);             %converts image to gray levels from 0-255
Shai
  • 111,146
  • 38
  • 238
  • 371
Jordon Birk
  • 480
  • 1
  • 9
  • 28
  • So what is the question/problem here? – Schorsch Jun 20 '13 at 18:17
  • That I need to do the same procedure, but using im2col (or some other form of block processing in order to change the structure to not involve nested for loops), and I have no idea how. – Jordon Birk Jun 20 '13 at 18:25

2 Answers2

3

Here is the code you want: note the entire nested loop was replaced with a single statement.

 [TRIMMED]=TM_FILTER(IMAGE, FILTER_SIZE, PERCENT)
    Takes a 2-D array and returns the array, filtered with a
    square trimed mean filter with length/width equal to FILTER_SIZE and percent equal to PERCENT.

    %}
    function [trimmed]=tm_filter(image, filter_size, percent)
    if rem(filter_size, 2)==0                            %make sure filter has a center pixel
        error('filter size must be odd numbered');       %error and return if number is odd
        return 
    end
    if percent > 100 || percent < 0
        error('Percentage must be ? [0, 100]');
        return
    end

    trimmed = (uint8)trimmean(im2col(image, filter_size), percent);

Explanation:

the im2col function turns each region of filter_size into a column. Your trimmean function can then operate on each of the regions (columns) in a single operation - much more efficient than extracting each shape in turn. Also note this requires only a single application of trimmean - in your original you first do it on the columns, then again on the rows, which will actually cause a more severe trim than I think you intended (exclude 50% first time, then 50% again - feels like excluding 75%. Not exactly true but you get my point). Also you will find that changing the order of operations (row, then column vs column, then row) would change the result because the filter is nonlinear.

For example

im = reshape(1:9, [3 3]);
disp(im2col(im,[2 2])

results in

1  2  4  5
2  3  5  6
4  5  7  8
5  6  8  9

since you took each of the 4 possible blocks of 2x2 from this matrix:

1  4  7
2  5  8
3  6  9

and turned them into columns

Note - with this technique (applied to the unpadded image) you do lose some pixels on the edge; your method added some padding so that every pixel (even ones on the edge) has a complete neighborhood, and as such the filter returns an image that is the same size as the original (but it's not clear what the effect of padding/filtering will be near the margin, and especially the corner: you have almost 75% percent of pixels fixed at 128 and that is likely to dominate the behavior in the corner).

Floris
  • 45,857
  • 6
  • 70
  • 122
2
  1. why im2col? why not nlfilter?

    >> trimmed = nlfilter( image, [filter_size filter_size],...
                           @(x) treimmean( trimmean(x, percent), percent ) );
    
  2. Are you sure you process the entire image?
    i and j only goes up to rows and columns respectively. However, when you update trimmed you access i-n and j-n. What about the last n rows and columns?

  3. Why do you apply trimmean twice for each block? Isn't it more appropriate to process the block at once, as in trimmean( x(:), percent)?
    I believe the results of trimmean( trimmean(x, percent), percent) will be different than those of trimmean( x(:), percent). Have you give it a thought?

  4. A small remark, it is best not to use i and j as variable names in matlab.

Community
  • 1
  • 1
Shai
  • 111,146
  • 38
  • 238
  • 371
  • 1
    "A small remark, it is best not to use i and j as variable names in matlab" because they have the built-in value `sqrt(-1)` until they are used otherwise. If you have any code that relies on their original meaning, it will break. Not a good habit. – Floris Jun 20 '13 at 18:29
  • 1. I'm not sure why. The documentation says nlfilter can be slow and that colfilt can be faster, but mot of the images I'm working with aren't exceptionally large. I do have to have these images processed in close to real time, eventually, so it might be worth considering...but, but then, I'll probably need to redo everything in C to get the speed. 2. Yes. It runs a trimmed mean filter, which is a sort of generalization of a mean and median filter that varies based on parameters. It's main goal is to weed out speckle noise in an image while preserving borders. 3. Oh yeah. – Jordon Birk Jun 20 '13 at 18:34
  • In actuality, I just based this off of a moving average filter that worked. It probably does not work for trimmed mean filters, especially since trimmed mean is a nonlinear operator. – Jordon Birk Jun 20 '13 at 18:42
  • The edits and comments that appeared while I was writing my answer suggests that we are mostly in complete agreement... – Floris Jun 20 '13 at 18:45
  • Two comments: first, cool that you "asked the question about i and j"; second, I think (after re-reading) that after the padding operation the region filtered is in fact correct (although the padding adds a bias towards 128 near the edge and especially the corner). – Floris Jun 20 '13 at 21:53