MATLAB: Correcting a Non-Uniform Background

image analysisimage processingImage Processing Toolboxnon-uniform illuminationpolyfit

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!

Best Answer

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

Related Solutions

MATLAB: Crop dark area and display the face only

Convert to grayscale and then threshold and find the bounding box

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. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% Read in the color image.
folder = pwd;
baseFileName = 'm2.bmp';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % 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.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
rgbImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows, columns, numberOfColorBands] = size(rgbImage);
% Display the original color image.
subplot(2, 3, 1);
imshow(rgbImage, []);
axis on;
title('Original Color Image', 'FontSize', fontSize);
drawnow;
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Histogram the red channel, but there is a huge spike at 255 so let's not count those.
subplot(2, 3, 2);
grayImage = rgb2gray(rgbImage);
imshow(grayImage, []);
title('Gray Scale version', 'FontSize', fontSize);
subplot(2, 3, 3);
histogram(grayImage(grayImage<255))
grid on;
title('Histogram of Gray Scale Image', 'FontSize', fontSize);
someThresholdValue = 150;
binaryImage = grayImage > someThresholdValue;
% Take largest blob
binaryImage = bwareafilt(binaryImage, 1);
subplot(2, 3, 4);
imshow(binaryImage);
grid on;
title('Binary Image', 'FontSize', fontSize);
% Label the binary image.
labeledImage = bwlabel(binaryImage);
props = regionprops(labeledImage, 'BoundingBox');
bbox = props.BoundingBox
croppedImageRGB = imcrop(rgbImage, bbox);
subplot(2, 3, 5);
imshow(croppedImageRGB);
grid on;
title('Cropped RGB Image', 'FontSize', fontSize);
% Gray scale version of cropped image
croppedImageGray = imcrop(grayImage, bbox);
subplot(2, 3, 6);
imshow(croppedImageGray);
grid on;
title('Cropped Gray Scale Image', 'FontSize', fontSize);

MATLAB: Normalize image orientation and scale

To apply PCA to an image, see this link: http://www.mathworks.com/matlabcentral/answers/167938#answer_163231 I don't think it will work, though you can try. Make sure you pass it the location of the black pixels, not the image gray levels.

I'd probably use a radon transform. See my attached demo where I use the radon transform to figure out the angle and rotate the image.

% Demo to use the radon transform to determine the angle to rotate an image to straighten it.
% Initialization / clean-up code.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% Read in the color demo image.
[rgbImage, colorMap] = imread('football.jpg');
subplot(2, 3, 1);
imshow(rgbImage, colorMap);
axis on;
title('Original Color Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Extract the red channel and display it.
grayImage = rgbImage(:, :, 1);
subplot(2, 3, 2);
imshow(grayImage, colorMap);
axis on;
title('Red Channel Image', 'FontSize', fontSize);
% Do the Radon transform.
theta = 0:180;
[R,xp] = radon(grayImage,theta);
% Find the location of the peak of the radon transform image.
maxR = max(R(:));
[rowOfMax, columnOfMax] = find(R == maxR)
% Display the Radon Transform image.
h3 = subplot(2, 3, [3,6]);
imshow(R,[],'Xdata',theta,'Ydata',xp,...
            'InitialMagnification','fit')
axis on;
% Plot a blue circle over the max.
hold on;
plot(h3, columnOfMax, xp(rowOfMax), 'bo', 'MarkerSize', 30, 'LineWidth', 3);
line([columnOfMax, columnOfMax], [xp(end), xp(rowOfMax)+15], 'Color', 'b', 'LineWidth', 3);
caption = sprintf('Radon Transform.  Max at angle %.1f', columnOfMax);
title(caption, 'FontSize', fontSize);
xlabel('\theta (degrees)', 'FontSize', fontSize)
ylabel('x''', 'FontSize', fontSize)
colormap(h3, hot(256));
colorbar;
% The column of the max is the angle of the football --
% the angle that the projected sum (profile) will have the highest sum.
% Rotate it by minus that angle to straighten it.
rotatedImage = imrotate(rgbImage, -columnOfMax);
% Display the rotated image.
subplot(2, 3, 4);
imshow(rotatedImage);
axis on;
title('Rotated Color Image', 'FontSize', fontSize);
% Rotate perpendicular to that angle and display that rotation.
rotatedImage = imrotate(rgbImage, -columnOfMax+90);
subplot(2, 3, 5);
imshow(rotatedImage);
axis on;
title('Rotated Color Image', 'FontSize', fontSize);

Best Answer

Related Solutions

MATLAB: Crop dark area and display the face only

MATLAB: Normalize image orientation and scale

Related Question