Here's my demo that uses polyfitn() to background correct an image. Save the code as BGCorrect_demo.m and run it. Requires the Image Processing Toolbox, and polyfitn().
% Demo to background correct a MATLAB demo image.
% Estimates the background from the image itself.
% Then models it to a 4th order polynomial in both directions
% using John D'Errico's File Exchange submission:
% http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn
function BGCorrect_demo()
try
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.
fontSize = 20;
% Make sure polyfitn is on my path.
addpath(genpath('C:\Program Files\MATLAB\work'));
% Read in the noisy nonuniform background image.
% Read in a standard MATLAB gray scale demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'AT3_1m4_01.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);
[rows columns numberOfColorChannels] = size(grayImage)
% Display the original background image.
subplot(2, 3, 1);
imshow(grayImage); % Display image.
title('Original Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize'));
drawnow; % Force display to update immediately.
% Let's process this to get the background, which will be the smooth areas.
filteredImage = entropyfilt(grayImage, ones(15));
subplot(2, 3, 2);
imshow(filteredImage, []); % Display image.
title('Entropy-Filtered Image', 'FontSize', fontSize);
binaryImage = filteredImage > 5;
binaryImage = imfill(binaryImage, 'holes');
subplot(2, 3, 3);
imshow(binaryImage, []); % Display image.
title('Binary Image', 'FontSize', fontSize);
% Fill in the textured areas with local background
% to get an estimate of what the background
% would be without the objects in there.
nonuniformBackgroundImage = roifill(grayImage, binaryImage);
subplot(2, 3, 4);
imshow(nonuniformBackgroundImage, []); % Display image.
title('Non-Uniform Background Image', 'FontSize', fontSize);
% Now we have an otiginal gray scale image and a background image.
% If you're able to take a "blank shot" - an image with no objects
% in there in the first place - that would be better than trying
% to estimate a background from an image with objects in it.
% Anyway, we now have our images and we're ready to begin!
% Get the modeled background image.
% Fit our actual background with a 2D polynomial.
modeledBackgroundImage = GetModelBackgroundImage(nonuniformBackgroundImage);
% Scale the model for display so that it looks more like the original.
meanOfOriginal = mean(nonuniformBackgroundImage(:))
meanOfModel = mean(modeledBackgroundImage(:))
modelImageForDisplay = modeledBackgroundImage * meanOfOriginal / meanOfModel;
meanOfDisplayModel = mean(modelImageForDisplay(:))
% Display the model background image.
subplot(2, 3, 5);
modelImageForDisplay = uint8(modelImageForDisplay);
imshow(modelImageForDisplay, []); % Display image.
title('Model Background', 'FontSize', fontSize);
% Do background correction on the raw background image.
correctedBackgroundImage = BackgroundCorrect(nonuniformBackgroundImage, modeledBackgroundImage);
subplot(2, 3, 6);
cla reset;
imshow(correctedBackgroundImage); % Display image.
caption = sprintf('Corrected background image');
title(caption, 'FontSize', fontSize);
% Bring up a new figure to show some more image.
figure;
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize'));
drawnow; % Force display to update immediately.
% Show the original images, the background image and the actual image.
subplot(2, 3, 1);
imshow(modelImageForDisplay, []); % Display image.
title('Model Background Image.', 'FontSize', fontSize);
subplot(2, 3, 2);
imshow(nonuniformBackgroundImage, []); % Display image.
title('Original Background Image.', 'FontSize', fontSize);
subplot(2, 3, 3);
imshow(grayImage, []); % Display image.
title('Original Gray Scale Image.', 'FontSize', fontSize);
% Show the corrected background image.
subplot(2, 3, 5);
imshow(correctedBackgroundImage); % Display image.
caption = sprintf('Corrected background image');
title(caption, 'FontSize', fontSize);
% Do background correction on the original gray scale image.
correctedGrayImage = BackgroundCorrect(grayImage, modeledBackgroundImage);
subplot(2, 3, 6);
imshow(correctedGrayImage, []); % Display image.
title('Corrected Gray Scale image.', 'FontSize', fontSize);
message = sprintf('Done with demo!\nBe sure to check out both figure windows.');
msgbox(message);
catch ME
errorMessage = sprintf('Error in function BGCorrect_demo().\n\nError Message:\n%s', ME.message);
fprintf(1, '%s\n', errorMessage);
uiwait(warndlg(errorMessage));
end
return; % from BGCorrect_demo
%================================================================================================================
% Divides the test image by the modeled image to get a background corrected image.
function correctedImage = BackgroundCorrect(inputImage, modeledBackgroundImage)
try
correctedImage = inputImage;
[rows columns numberOfColorChannels] = size(inputImage);
if numberOfColorChannels > 1
% It's a color image. Correct each color channel one at a time.
for colorIndex = 1 : 3
oneColorChannel = inputImage(:, :, colorIndex);
% maxNoisyValue = max(max(oneColorChannel));
% Model it to get rid of noise, scratches, smudges, etc.
noiselessBackgroundImage = modeledBackgroundImage(:, :, colorIndex);
% Divide the actual image by the modeled image.
% maxValue = max(max(noiselessBackgroundImage));
% correctedColorChannel = uint8(maxValue * single(oneColorChannel) ./ noiselessBackgroundImage);
correctedColorChannel = uint8(single(oneColorChannel) ./ noiselessBackgroundImage);
correctedImage(:, :, colorIndex) = correctedColorChannel;
end
else
% It's actually a gray scale image.
% Divide the actual image by the modeled image.
% maxValue = max(max(noiselessBackgroundImage));
correctedImage = uint8(single(inputImage) ./ modeledBackgroundImage);
end
catch ME
errorMessage = sprintf('Error in function BackgroundCorrect().\n\nError Message:\n%s', ME.message);
fprintf(1, '%s\n', errorMessage);
uiwait(warndlg(errorMessage));
end
return; % from BackgroundCorrect
%================================================================================================================
% Takes a color, noisy, non-uniform background image and fits a 2D model to it.
% Returns a noise-free, smooth color image where each color plane is independently
% normalized with a max value of 1. See lengthy comment farther down.
% modeledBackgroundImage is a double image in the range 0-1.
function modeledBackgroundImage = GetModelBackgroundImage(nonuniformBackgroundImage)
try
% Get the dimensions of the image. numberOfColorBands should be = 3.
% [rows columns numberOfColorBands] = size(nonuniformBackgroundImage);
% Preallocate output image.
modeledBackgroundImage = zeros(size(nonuniformBackgroundImage));
[rows columns numberOfColorChannels] = size(nonuniformBackgroundImage);
if numberOfColorChannels > 1
% It's a color image. Correct each color channel one at a time.
for colorIndex = 1 : 3
oneColorChannel = nonuniformBackgroundImage(:, :, colorIndex);
% Model it to get rid of noise, scratches, smudges, etc.
noiselessImage = ModelOneColorChannel(oneColorChannel);
% Divide the actual image by the modeled image.
maxValue = max(max(noiselessImage));
% Divide by the max of the modeled image so that you can get a "percentage" image
% that you can later divide by to get the color corrected image.
% Note: each color channel may have a different max value,
% but you do not want to divide by the max of all three color channel
% because if you do that you'll also be making all color channels have the same
% range, in effect "whitening" the image in addition to correcting for non-uniform
% background intensity. This is not what you want to do.
% You want to correct for background non-uniformity WITHOUT changing the color.
% That can be done as a separate step if desired.
noiselessImage = noiselessImage / maxValue; % This is a percentage image.
% Now insert that percentage image for this one color channel into the color percentage image.
modeledBackgroundImage(:, :, colorIndex) = noiselessImage;
% minValue = min(modeledBackgroundImage(:))
% maxValue = max(modeledBackgroundImage(:))
end
else
% It's a gray scale image. (Much simpler situation.)
% Model it to get rid of noise, scratches, smudges, etc.
noiselessImage = ModelOneColorChannel(nonuniformBackgroundImage);
% Divide the actual image by the modeled image.
maxValue = max(max(noiselessImage));
% Divide by the max of the modeled image so that you can get a "percentage" image
% that you can later divide by to get the color corrected image.
modeledBackgroundImage = noiselessImage / maxValue; % This is a percentage image.
% minValue = min(modeledBackgroundImage(:))
% maxValue = max(modeledBackgroundImage(:))
end
catch ME
errorMessage = sprintf('Error in function GetSmoothedBackgroundImage().\n\nError Message:\n%s', ME.message);
fprintf(1, '%s\n', errorMessage);
uiwait(warndlg(errorMessage));
end
% Free up memory. Should be automatically freed but we have had problems here.
% (Make sure you don't clear any input or output arguments or you'll get an exception.)
clear('noiselessImage', 'oneColorChannel');
return; % from GetSmoothedBackgroundImage
%================================================================================================================
% Fit a biquartic model to the image.
% You need to have a path to polyfitn(), John D'Errico's File Exchange submission.
function modeledColorChannel = ModelOneColorChannel(colorChannel)
try
modeledColorChannel = colorChannel; % Initialize.
rows = size(colorChannel, 1);
columns = size(colorChannel, 2);
% midX = columns / 2;
% midY = rows / 2;
[X,Y] = meshgrid(1:columns, 1:rows);
z = colorChannel;
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
z1d = double(reshape(z, numel(z), 1));
x = [x1d y1d];
%--------------------------------------------------------
% Get a 4th order polynomial model.
% CHANGE THE ORDER HERE IF YOU WANT TO.
%--------------------------------------------------------
polynomialOrder = 4;
p = polyfitn(x, z1d, polynomialOrder);
% Evaluate on a grid and plot:
zg = polyvaln(p, x);
modeledColorChannel = reshape(zg, [rows columns]);
catch ME
errorMessage = sprintf('Error in ModelOneColorChannel():\n%s', ME.message);
if strfind(errorMessage, 'HELP MEMORY')
memory;
errorMessage = sprintf('%s\n\nCheck console for memory information.', errorMessage);
end
uiwait(warndlg(errorMessage));
end
% Free up memory. Should be automatically freed but we have had problems here.
% (Make sure you don't clear any input or output arguments or you'll get an exception.)
clear('zg', 'x1d', 'y1d', 'z1d', 'X', 'Y', 'x', 'z', 'p');
return; % from ModelOneColorChannel
