non-uniform illumination - Scenes are changing!

Hi

I am trying to calculate/analyze coefficient of variance (CV) for multiple series of images and compare them at the end. Please see the code below.

The foreground (various shades of green that is fluorescent fluid) is separated from the background (dark blue) via Otsu. All the pixels that fall below the threshold are assigned with 0 intensity. The intensity of those that are higher than the threshold, is preserved. mean intensity is then calculated, around which the standard deviation (STD) is calculated. The STD is then divided by the mean intensity of the image. This process repeats for each image in a series (15 images). The area coverage and pattern of the images are different from one another. I have attached a full series analysis where CV, histogram, mean intensity, threshold value, and the image is presented. I have also attached a few images from another series so that the difference in illumination becomes clearer.

However, the illumination is non-uniform ( in both vertical and horizontal directions) within each series and from image to image. This causes irregularities in CV measurement. I have looked at "imhistmatch", white balance correction, gray level correction, homomorphic filters, and median filters. Honestly, I cannot think of any other method. None of them works because the scenes are changing from image to image. I truly appreciate any feedback/guidance on the matter.

P.s: No background-only or foreground-only image was taken. All images have both foreground and background.

fileList = dir('*.jpg');
fileList = {fileList.name};
thr = [];
im = imread(fileList{1});
img = im(:,:,2);
refSize = size(img);
CV = [];
num_of_soil_pixels = [];
num_of_foam_pixels = [];
mean_collect = [];
thr_collect = [];
sum_foam_intensities = [];
std_only =[];

%%%%% Defining circular areas over which calcs occur %%%%%%%%%%%%%
diam = mean(refSize);
%%%%% Iteration process to compute and store CV values for each image %%%%%
for i = 1:length(fileList)
    im = imread(fileList{i});
    fileList{i};
    % Resize images with respect to the first image
    if size(im(:,:,2))~= refSize
            im_resized = imresize(im, refSize);
    else
        im_resized = im;
    end
    % Convert to grayscale
    img = rgb2gray(im_resized);
    % Threshold the image using Otsu
    thr = graythresh(img);
    segmented_img = imbinarize(img, thr);
    out_img = double(img) .* segmented_img;
    % Normalize the image via dividing by the range
    out_img = 255*(out_img - min(min(out_img)))/(max(max(out_img))-min(min(out_img)));
    % Create a circular mask to obtain the region of interest (circle)
    mask = zeros(refSize);
    [xgrid, ygrid] = meshgrid(1:refSize(2), 1:refSize(1));
    mask = ((xgrid-(refSize(1)/2)).^2 + (ygrid-(refSize(1)/2)).^2) <= (diam/2.2).^2;
    
    % compare thresholded image with the RGB version - validate how well
    % Otsu detects the green regions (foreground)
    %fh2 = figure(23);
    %fh2.WindowState = 'maximized';
    %subplot(5,3,i)
    %imshowpair(im_resized,uint8(out_img),'montage')

    % mask the grayscale image
    values = out_img(mask);
    % Calculate the mean
    mymean = mean(values);
    
    % store various variables for plotting purposes
    mean_collect = [mean_collect;mymean];
    thr_collect = [thr_collect;thr*255];
    CV = [stdev; std(values)/mymean];
    std_only = [std_only;std(values)];
    num_of_soil_pixels = [num_of_soil_pixels;sum(values==0)];
    num_of_foam_pixels = [num_of_foam_pixels;nnz(values)];
    sum_foam_intensities =[sum_foam_intensities;sum(values(values~=0))];
end 

figure (1)
plot(CV, 'ok','MarkerfaceColor','g')
ylim([0,4])
xlim([0,15])
grid on

Answer 1

[Image Analyst]

@thelabmaster, to get the mean and standard deviation only in the masked region, try this:

% Demo by Image Analyst, November 24, 2021.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clearvars;
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
fprintf('Beginning to run %s.m ...\n', mfilename);
%-----------------------------------------------------------------------------------------------------------------------------------
% Read in image.
folder = [];
baseFileName = '03.jpg';
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~isfile(fullFileName)
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, '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
    fullFileName = fullFileNameOnSearchPath;
end
rgbImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Display the image.
subplot(2, 2, 1);
imshow(rgbImage, []);
axis('on', 'image');
caption = sprintf('Original Image : "%s"', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
hFig1 = gcf;
hFig1.Units = 'Normalized';
hFig1.WindowState = 'maximized';
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
hFig1.Name = 'Demo by Image Analyst';
%--------------------------------------------------------------------------------------------------------
% Segment (mask) the image.
[mask, maskedRGBImage] = createMask(rgbImage);
% Display the mask image.
subplot(2, 2, 2);
imshow(mask, []);
axis('on', 'image');
caption = sprintf('Mask Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Display the masked image.
subplot(2, 2, 3);
imshow(maskedRGBImage);
axis('on', 'image');
caption = sprintf('Masked Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Compute intensity image
grayImage = rgb2gray(maskedRGBImage);
% Compute mean in the mask region.
meanGL = mean(grayImage(mask))
% Compute standard deviation in the mask region.
stdGL = std(double(grayImage(mask)))


%====================================================================================================
function [BW,maskedRGBImage] = createMask(RGB)
%createMask  Threshold RGB image using auto-generated code from colorThresholder app.
%  [BW,MASKEDRGBIMAGE] = createMask(RGB) thresholds image RGB using
%  auto-generated code from the colorThresholder app. The colorspace and
%  range for each channel of the colorspace were set within the app. The
%  segmentation mask is returned in BW, and a composite of the mask and
%  original RGB images is returned in maskedRGBImage.
% Auto-generated by colorThresholder app on 24-Nov-2021
%------------------------------------------------------
% Convert RGB image to chosen color space
I = rgb2hsv(RGB);
% Define thresholds for channel 1 based on histogram settings
channel1Min = 0.247;
channel1Max = 0.565;
% Define thresholds for channel 2 based on histogram settings
channel2Min = 0.272;
channel2Max = 1.000;
% Define thresholds for channel 3 based on histogram settings
channel3Min = 0.236;
channel3Max = 1.000;
% Create mask based on chosen histogram thresholds
sliderBW = (I(:,:,1) >= channel1Min ) & (I(:,:,1) <= channel1Max) & ...
    (I(:,:,2) >= channel2Min ) & (I(:,:,2) <= channel2Max) & ...
    (I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max);
BW = sliderBW;
% Initialize output masked image based on input image.
maskedRGBImage = RGB;
% Set background pixels where BW is false to zero.
maskedRGBImage(repmat(~BW,[1 1 3])) = 0;
end

You'll see in the command window:

meanGL =

141.995641629468

stdGL =

42.7072939153801

[thelabmaster]

Thank you, @Image Analyst. I appreciate your help on this. While I am using your script on the images, I wanted to ask if you have suggestions on how to correct illumination irregularities that pervail among images.

[Image Analyst]

If your UV illumination is not uniform then you'll have to get a map of the illumination and divide out the intensity. For example, if only 70% as many UV photos hit some pixel in your image, it will look darker than it should obviously. So how to boost the signal there as if it had the full 100% of UV light available? Well, divide by 0.7. So if the intensity at a place with 70% as much light had a gray level of 140, then if it had the full lighting at 100% it would be 140/0.7 = 200.

So to get a map of the illumination you need to put down some kind of uniform fluorescent standard covering your entire field of view. These exist, just Google for them. Then take a photo of it. Then find the max intensity and divide that image by the max intensity. This will turn the image into a percentage (0 to 100%) that each pixel in the image has of the full light intensity. Then when you take a test shot with your sample in there, divide that image by the percentage image.

[thelabmaster]

Thank you for explaining, @Image Analyst. Seems crucial to correct images with respect to a reference image; however, no uniformly-distributed UV image was taken as a reference (and no more experiments are expected to run). Therefore, the max intensity cannot be established against which other images are to be compared (if I understand correctly).

[Image Analyst]

What makes you think the illumination is not uniform then? Maybe the non-uniform signal is truly due to the amount of fluorescent dye in the sample and not the lighting. Or maybe it's due to lens shading. But I'm sure you're correcting for lens shading, right?

[thelabmaster]

Because there is an illumination gradient that is typically max at the top and decreases towards the bottom of the image. To a lesser degree, the same statement holds true from left to right. Also the UV light posts were placed manually each time so that introduces some error, e.g., their locations might slightly be different.

The intensity of fluorescent dye certainly changes from image to image. I normalize the gamma and that is about it. Did I understand you correctly about correcting for lens shading?

[Image Analyst]

You can correct for lens shading just by taking a picture of a white light image that's totally blank white and uniform and illuminated not by UV but just with broad source illumination, like a lihgting dome, or maybe even your overhead lab lighting. Then determine the percentage of light at each point by dividing the blank image by the max brightness of the image. Then divide any test images by this percentage image.

So why can't you buy a fluorescent standard and put it in the field of view and take a photo of it to determine the uniformity of illumination? Maybe just try a really white sheet of photocopy paper that has lots of optical brightener in it.

[thelabmaster]

Understood. However, the experiments are done a while ago and no more photos can be taken at this point.

Answer 2

[yanqi liu]

sir，may be do some image registration to register images

[thelabmaster]

Thank you, Yanqi. Researched it but not sure how registration may help overcome bi-directional irregularities of illumination. Can you please elaborate?

[yanqi liu]

yes，sir，may be use sift or surf points to make image match