Correcting a Non-Uniform Background

I created a sample image by doing the following:

%%Create a Sample Image
for i = 1:501
  x(i,:) = linspace(0,.5,501);
end

% display image
figure(1), clf; imshow(x);

Now, I want to know how I can use a polynomial fit or some other trick to correct the non-uniform background. I'm trying to model something fairly robust, so that I can use it with a real image.

Thank you!

Answer 1

[Image Analyst]

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

[Caleb]

That was an outstanding demo!

Could you have made the background image darker before overlaying the blobs back on the image without messing with the regions of interest?

[Image Analyst]

Only by masking. But then that's not really background correction by dividing by the background. That's really segmentation into foreground and background and assigning the background to something else since it wouldn't be applying the correction to the foreground that it should be.

By the way,there are two ways to do background correction, by dividing, and by subtraction. You have to understand the physics and optics of what's going on to decide which is the theoretically correct way to do it. For most situations where you have light reflecting off a scene, you'd use background division. For radiological imaging and fluorescence imaging you'd use background subtraction.

[Caleb]

That makes sense. If you're working with an image that does not easily segment the objects in the foreground from the background, then what can you do?

Take, for example, the image you worked with. What if some of those "blobs" were darker than the background in some areas due to the background not being uniform? i.e. if the background in the center of the image is actually brighter than the blobs around the edges, how can you still segment these effectively?

Also, I've been having a whale of a time with the shapes of my blobs due to some light pixels around the edges. Is there a way to smudge the region or improve the quality of the segmentation?

[Image Analyst]

You have to employ sophisticated methods. The harder it is to tell the foreground from background, the smarter the algorithm will have to be. Blogs on the edges being darker than the background in the middle can be solved in different ways depending on the size, shape, and texture of the objects. You could use a texture filter like I did in my demo, or use an edge filter followed by imfill, or you could use morphological methods like a tophat or bottomhat filter, or model your background assuming the objects' contribution won't affect the model that much and then threshold, or maybe color segmentation or something else. You could upload your image to tinypic.com if you want me to take a look at it and give opinion.

[Caleb]

I'm working with a .tif image, and I don't have the software to convert it to jpeg, so I'm not sure that I can upload it.

[Ryan]

imgur.com allows you to upload .tif images

[Jay]

@Image Analyst: What is the exact algorithm you are using in this demo. Would you be able to write down the steps of this algorithm. I just want to have the overview of the demo. thanks.

[Image Analyst]

The exact algorithm is shown. In short, it's background division.

[Parth]

@Image Analyst please comment on this. https://www.mathworks.com/matlabcentral/answers/320992-what-kind-of-bias-actually-surface-fitting-is-correcting

Answer 2

[Image Analyst]

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.

[Caleb]

You've referred me to that before, but I wasn't really sure how to use it. Do I need to modify the function in some way or does it just act on an existing image without my help?

[Image Analyst]

See my demo elsewhere on this page.

Answer 3

[Teja Muppirala]

One example that uses morphological operations is here:

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

Answer 4

[Parth]

Channel wise normalization would introduce inconsistent hue in the dataset when the lighting conditions are uncorrelated. Isn't it?

In BackgroundCorrect it is normalizing it per channel.

[Image Analyst]

If you did not have the same pattern of light in each color channel, then it's not white to begin with and there are already local variations in color. Doing normalization on each color channel independently will flat field each color channel thus making the whole image one uniform color, which might not be necessarily having red=green=blue. Since it's all one uniform color, of course the color will vary from the original non-uniform pattern of colors you started with. But at least you have a flat field. If you then want to apply color correction/matching (RGB to RGB) or color calibration (RGB to CIELAB) then you can do that as a second step.