How to account for non-uniform illumination - Scenes are changing from image to image

Question

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.

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.

First series of images

Second series of images

Full analysis (histogram + image + CV + mean intensity + threshold) of the first series of the images

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