Correcting a Non-Uniform Background

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

I use John D'Errico's polyfitn: http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn to fit a 2D polynomial to some smooth but noisy background image.

One example that uses morphological operations is here:

http://www.mathworks.com/products/image/examples.html?file=/products/demos/shipping/images/ipexrice.html