%% Step 2: Image subtraction 

% Author: Eleanor Fadely, University of California Davis
% Last updated: 01/27/2026

% Key steps:

% (1) Load cropped/registered color brightfield images from mineralization
% phase. 
% (2) Make a binary mask of each color brightfield image.
% (3) Mask reference image and color brightfield images to isolate
% biofilm-biomineral region of interest at each timepoint.
% (4) Subtract all masked color brightfield images from masked reference 
% image to capture decrease in pixel intensity due to Mn oxide
% precipitation.
% (5) Sum each subtracted matrix to measure the total change in intensity
% at each timepoint.
% (6) Save variables for use in subsequent scripts.

% --------------- Set up file paths --------------------------------------

clear all; clc; close all

% For PC
addpath('F:\');

% -------------- Inputs (need to modify) ---------------------------------

% Specific experiment (CHANGE THIS)
ExperimentName = '240923_Constit_mCh';
FolderPath = 'Oxidation';

% --------- Set up directories and image retrieval from folder -----------

% Root directory is external hard drive, will never change
RD = 'F:\';

% Assign pathways for variable storage and retrieval
VDG = [RD ExperimentName '/Growth/Protocol_Variables'];
VDO = [RD ExperimentName '/Oxidation/Protocol_Variables'];

PIO = [VDO '/Processed_Images'];
GM = [VDG '/Grain_Mask'];
IS = [VDO '/Image_Subtraction'];

%% Load brightfield images from mineralization phase and grain mask

cd(PIO);
load('CroppedOxidationPhase.mat');
  
cd(IS);
load('ReferenceImageNumber.mat');
load('EndingImageNumber.mat');

cd(GM);
load('GrainMask.mat');
cd(VDO);

'Done loading...'

%% Make mask of each brightfield image to capture mineral accumulation

% Loop through brightfield images and binarize
for ii = 1:EndNum
    GrayIm = rgb2gray(ProcOxidationBF{1,ii});
    TR = adaptthresh(GrayIm, 0.55, 'ForegroundPolarity', 'dark');
    IM = imbinarize(GrayIm, TR);
    Cleaned = imcomplement(GM_Final) .* imcomplement(IM);
    ImageMask{1,ii} = logical(Cleaned);

    'Image masked...'

end

%% Mask all mineralization phase images and subtract subsequent timepoints from reference

% Identify reference image
RefImage = ProcOxidationBF{1,RefNum};

% Set up arrays to populate
dim = EndNum - RefNum;
dROIInt = zeros(3,dim);
dGrayInt = zeros(1,dim);

for m = RefNum:EndNum

    % Assign mask and image to loop through
    TimepointMask16 = uint16(ImageMask{1,m});
    Image = ProcOxidationBF{1,m};

    % Apply timpeoint mask to image
    RedROI = TimepointMask16 .* Image(:,:,1);
    GreenROI = TimepointMask16 .* Image(:,:,2);
    BlueROI = TimepointMask16 .* Image(:,:,3);
 
    % Multiply reference image by timepoint mask
    RefROI = TimepointMask16 .* RefImage;

    % Split masked reference image into 3 color channels
    RefR = RefROI(:,:,1);
    RefG = RefROI(:,:,2);
    RefB = RefROI(:,:,3);

    % Subtract masked timepoint image from masked reference image
    SubR = imsubtract(RefR, RedROI);
    SubG = imsubtract(RefG, GreenROI);
    SubB = imsubtract(RefB, BlueROI);

    % Convert separate color channels back to 3D image --> rgbImage = cat(3, R, G, B)
    subtractedImage = cat(3, SubR, SubG, SubB); 

    % Convert color image to grayscale (1D)
    subtractedGray{m} = rgb2gray(subtractedImage);

    % Sum change in pixel intensity at each timepoint 
    dROIInt(1,m-(RefNum-1)) = sum(SubR, 'all'); % red channel
    dROIInt(2,m-(RefNum-1)) = sum(SubG, 'all'); % green channel
    dROIInt(3,m-(RefNum-1)) = sum(SubB, 'all'); % blue channel
    dGrayInt(1,m-(RefNum-1)) = sum(subtractedGray{m}, 'all'); % grayscale composite

    'Image masked and subtracted...'
   
end

%% Plot accumulation of Mn oxide product over time 

gc = [0.5020 0.5020 0.5020]; % color for grayscale composite data plotting
TP = 0:0.5:((EndNum/2)-0.5); % time vector for mineralization phase

f1 = figure(1);
set(gcf,'Position',[100 100 800 700]);
hold on;
plot(TP(1,RefNum+1:end), dROIInt(1,2:end), '-o', 'Color', 'r', 'MarkerSize', 8, 'LineWidth', 1.5, 'DisplayName', 'Red channel');
plot(TP(1,RefNum+1:end), dROIInt(2,2:end), '-o', 'Color', 'g', 'MarkerSize', 8, 'LineWidth', 1.5, 'DisplayName', 'Green channel');
plot(TP(1,RefNum+1:end), dROIInt(3,2:end), '-o', 'Color', 'b', 'MarkerSize', 8, 'LineWidth', 1.5, 'DisplayName', 'Blue channel');
plot(TP(1,RefNum+1:end), dGrayInt(1,2:end), '-o', 'Color', gc, 'MarkerSize', 8, 'LineWidth', 1.5, 'DisplayName', 'Grayscale composite');
hold off;
xlabel('Time (h)', 'FontSize', 20);
ylabel('Relative pixel intensity', 'FontSize', 20);
title('Mn oxides produced over time');
legend SHOW;
legend('Location', 'NorthWest');
box on;
set(gca, 'FontSize', 16, 'LineWidth', 1.25);
xlim([0 30]);

%% Save variables 

cd(IS);
save('Brightfield_Image_Masks.mat', 'ImageMask', '-v7.3');      % masks of all brightfield image timepoints
save('Subtracted_Images_Ref.mat', 'subtractedGray', '-v7.3');   % matrix of subtracted images 
save('RGB_Total_Amount.mat', 'dROIInt');                        % sum of relative intensity in R, G, B channels
save('Grayscale_Total_Amount.mat', 'dGrayInt');                 % sum of relative intensity in grayscale composite
