Sliding window minimum and maximum filter

I'm trying to apply a sliding window minimum and maximum filter to an image of a certain window size. Actually, I'm trying to find the optimum window size for it. But I really haven't gotten the hang of it. I presume that I should be using blockproc to implement the sliding window, but not really sure how to find the maximum and minimum filter. As to the implementation itself, should I use loops to slide the window across the entire area of the image ?

[Image Analyst]

blockproc() moves in "jumps" of the window size. Imerode and imdilate slide over by one pixel.

[Matt J]

In fact, blockproc can be made to move in smaller jumps than the window size by using the BorderSize parameter

[doudou]

@Matt J plz i have a similaire problem, how can i do to make blockproc() to move in smaller jumps (one pixel) and this is a link that describe what i have to do: http://www.fp.ucalgary.ca/mhallbey/texture_calculations.htm

[Image Analyst]

From the description "A two-element vector, [V H], specifying the amount of border pixels to add to each block. The function adds V rows above and below each block and H columns left and right of each block." it looks like you'd have a windowSize of 1 so that it moves over in jumps of 1 pixel each time, but the 'BorderSize' would be half your window width so that it includes more than just the one pixel it's centered on. So if you wanted a 5 by 5 window, you'd set a windowSize of 1 and a BorderSize of 2. I think - I haven't tried it.

[doudou]

@image Analyst thanks for your answer

Answer 1

[Image Analyst]

If you have the Image Processing Toolbox, you're in luck!. The sliding max filter is called imdilate() and the sliding min filter is imerode(). These are called "morphological operations." No loops needed:

localMinImage = imerode(grayImage, true(3));
localMaxImage = imdilate(grayImage, true(3));

[Lab Rat]

In fact, I do have the Image Processing Toolbox. How do I make sure that I use the minimum value of the pixels in the sliding window/image, for the pixels close to the image boundaries ?

[Image Analyst]

I did, above. At least I showed how to get them. How you USE them afterwards is up to you. If you need a more complete demo, see below. In that, I just "use" them by simply displaying them in the image.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;

% Change the current folder to the folder of this m-file.
if(~isdeployed)
  cd(fileparts(which(mfilename)));
end
  
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

% Read in a standard MATLAB gray scale demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'cameraman.tif';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Give a name to the title bar.
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off') 

% Let's compute and display the histogram.
[pixelCount grayLevels] = imhist(grayImage);
subplot(2, 2, 2); 
bar(pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Compute the local max and min images:
localMinImage = imerode(grayImage, true(3));
localMaxImage = imdilate(grayImage, true(3));
% Display them.
subplot(2, 2, 3);
imshow(localMinImage, []);
title('Local Min (eroded) Image', 'FontSize', fontSize);
subplot(2, 2, 4);
imshow(localMaxImage, []);
title('Local Max (dilated) Image', 'FontSize', fontSize);

[Lab Rat]

Basically, I am attempting to correct an image with uneven illumination by applying a sliding window minimum filter on the image and then finding the illumination or shading function, by applying sliding window maximum filter on on the output of the minimum filter. Later, I'm trying to find the 'perfect image' by dividing the original image with the illumination function.

[Image Analyst]

Doing a max filter on the output of the min filter is called a morphological opening and is done by imopen(). It is not a way to get the illumination pattern that I've ever heard of. I do background correction all the time and this would not be a good way of getting the background unless your background was uniform and had small light-colored dust on it. But even then, there are flaws with that method such that you will increase the noise and have quantization/posterization errors. Do you have a uniform background that you can image? Or all you have is the image with your scene and illumination combined? If, unfortunately, you have the latter case, you can either try homomorphic filtering (basically assuming the light pattern is a very low pass filtered version of your scene) or try to find "holes" in your scene that is pure background and then try to fit a nice smooth model, like a 2D polynomial, to them. Finally, division is the right way to do it in most cases. For certain other cases (radiology and fluorescence) you should do subtraction. Do the math and you'll see why. Upload an image if you want further advice.

[Lab Rat]

Its an assignment that I've been asked to do and that's the method that has been specified. It states that an image g(x,y) can be modeled as the product of a "perfect image" denoted by f(x,y) and an illumination function l(x,y). i.e g(x,y) = f(x,y) * l(x,y) So the assignment is to find l(x,y) using the method that I've stated in my earlier response. To implement a minimum sliding window filter on g(x,y) and then implementing a sliding window maximum filter on the output of the minimum filter. This way, I'll find l(x,y) and then find f(x,y) by dividing g(x,y) by l(x,y). I've attached the image which is basically a modified rice.tif image. https://dl.dropbox.com/u/793037/lab2.tif

[Image Analyst]

If the rice image is like the MATLAB rice image with a gradient background, then doing an opening would give you just the dark gradient background because the erosion eats away the rice and then the closing gets the dark things back to their original size. You can do that, but that is exactly the situation I was saying would give you amplified image noise and posterization artifacts. It might be okay in certain situations, but not those where you need high accuracy. You can do this assignment in two lines. One call to imopen(), and another to divide the images. Then of course some more stuff to fancy it up like display the image(s), give them titles, etc.

[Lab Rat]

I tried it out, but I'm getting a white image once I'm done dividing the result with the original image.

I=imread('lab2.tif');
imshow(I);
im=im2double(I);
se=strel('square',5);
lxy=imopen(I,se);
fl=im2double(lxy);
f=imdivide(im,fl);
imshow(f);

[Image Analyst]

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 13;

% Read in a standard MATLAB gray scale demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'rice.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
% Display the original gray scale image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Do an opening to get the background.
% This is not a good method for this image.
openedImage = imopen(grayImage, true(13));
subplot(2, 2, 2);
imshow(openedImage, []);
title('Background Image', 'FontSize', fontSize);

% Divide them
backgroundCorrected = grayImage ./ openedImage;
subplot(2, 2, 3);
imshow(backgroundCorrected, []);
title('Background-corrected Image', 'FontSize', fontSize);

[Lab Rat]

Thanks for that. It works exactly like I want it to. But suppose, if both my original image and illumination function are double precision values and if I'd want my output image also to be a double precision value, wouldn't my resultant image have minimum and maximum values out of range rather than between [0 1] ?

[Lab Rat]

Nvm, I think I figured it out. Thanks a lot for your help. I really appreciate it.

[Royi Avital]

What would be the fastest way to do so without the Image Processing Toolbox?

[Image Analyst]

Did you see Matt J's answer?

[Royi Avital]

Why would it be faster than mine (See down)?

[Image Analyst]

Your method using colfilt() requires the Image Processing Toolbox and above you asked specifically about one that does not use that toolbox. That's why I said Matt's answer. Disregarding that, I have not timed the various algorithms so I don't know how they compare.

[Royi Avital]

Yes, I know. But Matt's method creates a "Shift" in the image which requires padding which is again part of the IP toolbox.

I wonder if there's a code which can be faster than imdilate and require no toolbox.

[Matt J]

Padding doesn't require the IP toolbox, e.g.,

 p=2; %padding
 Apad=zeros(size(A)+2*p,class(A));
 Apad(p+1:end-p,p+1:end-p)=A;

[Image Analyst]

What do you mean by shift? As long as the window size is odd, e.g. 3 or 5, there will be no shift. If the window size is even, e.g. 2 or 4, then there will be a half pixel shift.

Answer 2

[Matt J]

I would just use a for-loop to do 2 separable passes. BLOCKPROC can't take advantage of the separable nature of the max/min filter:

    A=rand(100);
    window=3;

    [m,n]=size(A);
    B=A;

    for ii=1:m+1-window
       B(ii,:)=max(A(ii:ii+window-1,:),[],1); 
    end

    for ii=1:n+1-window
       B(:,ii)=max(B(:,ii:ii+window-1),[],2); 
    end

[Royi Avital]

This is nice and fast. Yet it shift the matrix to the right.

How would replicate the results of `imerode` or `imdilate` with Image Processing Toolbox (no `padaaray`) most efficiently?

Thaks.

Answer 3

[Royi Avital]

I would use:

localMaxImage = colfilt(inputImage, [winLength winLength], 'sliding', @max);
localMinImage = colfilt(inputImage, [winLength winLength], 'sliding', @min);

Though it is still requires patience. I wonder what would be the fastest way to do so without Image Processing Toolbox.

Answer 4

[tilak tenneti]

i want to apply a sliding window minimum filter on input image I and obtain Imin and also apply sliding window maximum filter on Imin to obtain Imax : the following is code

    N=1;
    Imin=ordfilt2(I, 1, true(N));

    N=10;
    Imax = ordfilt2(Imin, N*N, true(N));

here i assume N as window size... but i am confused as how should i take N for minimum filter and again N for maximum filter?

Answer 5

[Dan]

function [minVals,maxVals] = minmaxfilt1(vector,nhoodSz)
vector = vector(:);
if nhoodSz < 3 || ~floor(mod(nhoodSz,2)) 
    error('nhoodSz must be odd scalar');
end

minVals = min(conv2(vector,eye(nhoodSz)),[],2);
maxVals = max(conv2(vector,eye(nhoodSz)),[],2);

minVals = minVals(((nhoodSz-1)/2)+1: end- ((nhoodSz-1)/2));
maxVals = maxVals(((nhoodSz-1)/2)+1: end - ((nhoodSz-1)/2));
end