MATLAB: How Can I simulate non uniform background illumination

image processingImage Processing Toolbox

I would like to add non uniform background illumination to some images to test some thresholding algorithms.

Best Answer

If I understand your question correctly you would like to generate images with non-uniform background illumination in order to test thresholding algorithms.

You could try to find a test set of images which already have a non-uniformly illuminated background. For example, rice.png is an image shipped with Image Processing Toolbox with non-uniform illumination. Take a look at that example to see how you could use it:

Correcting Nonuniform Illumination

It's hard to add non-uniform background illumination to images. That would be hard even with Adobe Photoshop. It's easier to create a non-uniform background and then add the foreground. So, why not create images from the ground up using basic functions in MATLAB?

Here's a code snippet that creates a grayscale image with a non-uniform background made of a sinc wave:

    % Size of our image: N-by-N
    N = 500;
    % The xy limits of the sinc function
    % we want N points in each direction
    x = linspace(-50,50,N);
    y = linspace(-50,50,N);
    % Transform these xy coordinates into a grid pattern
    [X,Y] = meshgrid(x,y);
    % The input to the sinc function is the radius from the center
    r = sqrt(X.^2 + Y.^2);
    % Create the background as the output of sinc
    background = sin(r)./r;
    % Now we need to scale the matrix so that it is a valid image
    mmin = min(background(:));
    mmax = max(background(:));
    % Subtract the min, divide by the range, multiply by 255 -> values are in [0,255]
    background = 255*(background-mmin)/(mmax-mmin);
    % Round and cast into uint8 so that it is a valid grayscale image
    background = uint8(round(background));
    % Now display the result
    figure
    imshow(background)

    % Add objects in the foreground
    img = background;
    img(100:200,100:200) = 255;
    img(150:450,150:350) = 128;
    % Display the final image
    imshow(img)

You could also do a gradient across the image. For example:

    % Create a vector of linearly varying values in x from 100 to 200
    x = linspace(100,200,N);
    y = 1:N;
    % Transform into a grid
    [X,Y] = meshgrid(x,y);
    % Define the background
    background = uint8(round(X));
    % Display
    imshow(background)

    % Add objects in foreground
    ...

Related Solutions

MATLAB: 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

MATLAB: Thin the boundary of objects of black & white images

Use imerode.

From Matlab documentation:

IM2 = imerode(IM,SE) erodes the grayscale, binary, or packed binary image IM, returning the eroded image IM2. The argument SE is a structuring element object or array of structuring element objects returned by the strel function.

Example:

Erode a binary image with a disk structuring element.

originalBW = imread('circles.png');  
se = strel('disk',11);        
erodedBW = imerode(originalBW,se);
imshow(originalBW), figure, imshow(erodedBW)

Best Answer

Related Solutions

MATLAB: Correcting a Non-Uniform Background

MATLAB: Thin the boundary of objects of black & white images

Related Question