%% Step 1: Registration of time-lapse microscopy images

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

% Key steps:

% (1) Load first timelapse brightfield image, crop, and
% straighten/register relative to grain mask from CAD file.
% (2) Apply registration transform to all brightfield images from
% growth and mineralization phases.
% (3) Identify reference image (first timepoint where biofilms are
% fixed in place and not growing).
% (4) Save variables for use in subsequent scripts.

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

clear all; clc; close all
addpath('F:\');

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

% Specific experiment (CHANGE THIS)
ExperimentName = '240923_Constit_mCh';
FolderPathGR = 'Growth/Export_Color';
FolderPathOX = 'Oxidation/Export_Color';

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

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

% Assign and enter working directory
GrowthD = [RD ExperimentName '/' FolderPathGR];
OxidationD = [RD ExperimentName '/' FolderPathOX];

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

%% Read in CAD file of grain geometry 

cd(RD);
CI = imread('PoreSpace_Mask.tif');
cd(GrowthD);

% Resize CAD file to reactor dimensions (1 pixel = 1.62 um)
GrainMask = imresize(CI, [(10000/1.62), (15000/1.62)]); 
figure; imshow(GrainMask);

%% Read in first timelapse image (empty pore space) to assign registration

cd(GrowthD);

% Retrieve image names
ims = dir('*.tif');
ims = ims(~startsWith({ims.name},'._'));
ImageName = {ims.name};

RotationImage = imread(ImageName{1,1});
figure; imshow(RotationImage);

%% Crop rotation image around pore space to remove inlet/outlet/background

% INTERACTIVE STEP: Opens Crop Image GUI, need to export coordinates when done
imcrop(RotationImage);

% Update this by pasting cropping coordinates from GUI
CropCoord1 = [444.5 220.5 10770 7618];

%% Process cropped rotation image to make a rough grain mask 

% Crop image
RI = imcrop(RotationImage, CropCoord1);

% Sharpen image
RISharp = imsharpen(RI, Radius=2, Amount=2); 
figure; imshow(RISharp);

% Binarize sharpened image using an adaptive threshold that prioritizes dark features (grain boundaries)
T = adaptthresh(rgb2gray(RISharp), 0.58, 'ForegroundPolarity', 'dark');
BW = imbinarize(rgb2gray(RISharp), T);
BWC = imcomplement(BW);

% Fill grains roughly (not all grains will be filled)
Filled = imfill(BWC, 'holes');
figure; imshow(Filled);

% Erode image to remove background noise, then dilate to replace missing pixels
SE = strel('square', 2);
FilledE = imerode(Filled, SE); 
FilledERD = imdilate(FilledE, SE);
figure; imshow(FilledERD);

%% Use SIFT algorithm to identify matching features for image registration
    % See Matlab documentation for more details
    % https://www.mathworks.com/help/vision/ug/local-feature-detection-and-extraction.html 

close all;

% Detect features in CAD file (ground truth) and roughly binarized microfluidic image (approximation)
ptsOriginal = detectSIFTFeatures(GrainMask);
ptsDistorted = detectSIFTFeatures(FilledERD);

% Extract features
[featuresOriginal, validPtsOriginal] = extractFeatures(GrainMask, ptsOriginal);
[featuresDistorted, validPtsDistorted] = extractFeatures(FilledERD, ptsDistorted);

% Compare features to evaluate potential match
indexPairs = matchFeatures(featuresOriginal, featuresDistorted);
matchedOriginal = validPtsOriginal(indexPairs(:,1));
matchedDistorted = validPtsDistorted(indexPairs(:,2));

figure;  
showMatchedFeatures(GrainMask, FilledERD, matchedOriginal, matchedDistorted);
title('Candidate matched points (including outliers)');

% Remove false match points 
[tform,inlierIdx] = estgeotform2d(matchedDistorted, matchedOriginal,'affine');
inlierDistorted = matchedDistorted(inlierIdx,:);
inlierOriginal = matchedOriginal(inlierIdx,:);

figure;
showMatchedFeatures(GrainMask, FilledERD, inlierOriginal, inlierDistorted);
title('Matching points (inliers only)');
legend('ptsOriginal','ptsDistorted');

% Transform roughly binarized microfluidic image
outputView = imref2d(size(GrainMask));
recovered = imwarp(FilledERD, tform, OutputView=outputView);

% View overlay of roughly binarized microfluidic image/CAD file
figure; imshowpair(GrainMask, recovered);

% Test registration on initial brightfield image
recovered2 = imwarp(RI, tform, OutputView=outputView);
figure; imshowpair(GrainMask, recovered2);

%% Manually select matching points for better local fit 
    % ~ 40 points, focusing on areas where grain mask does not fit well

% INTERACTIVE STEP: Opens Control Point Selection GUI
    % Lock ratio, zoom to at least 200%
    % Export movingPoints and fixedPoints to workspace from GUI
cpselect(recovered2, GrainMask);  

%% Apply local registration function (uses linear weighted mean)

tform2 = fitgeotform2d(movingPoints, fixedPoints, "lwm", 40);
recovered3 = imwarp(recovered2, tform2, OutputView=outputView);
figure; imshowpair(GrainMask, recovered3);

%% Crop image a second time to remove registration artifacts

% INTERACTIVE STEP: Opens Crop Image GUI
imcrop(recovered3);

% Update this by pasting cropping coordinates from GUI
CropCoord2 = [10.5 9.5 9236 6158];

%% Crop and resize final registered image and grain mask

% Crop both
RecFin = imcrop(recovered3, CropCoord2);
GMFin = imcrop(GrainMask, CropCoord2);

% Resize both (to known dimensions of 1 x 1.5 cm)
RecRes = imresize(RecFin, [(10000/1.62), (15000/1.62)]);
GMRes = imresize(GMFin, [(10000/1.62), (15000/1.62)]);
figure; imshowpair(GMRes, RecRes);

% Dilate grain mask to better remove grain boundaries in subsequent processing
SE_GM = strel('square', 3);
GM_Final = imdilate(GMRes, SE_GM);
figure; imshowpair(GM_Final, RecRes);

%% Register all brightfield images from growth phase of experiment

close all;

tic

% Load all images, register, crop, store into structure
for i = 1:length(ImageName)

    Image = imread(ImageName{1,i});                                     % read in image
    ImageInitCrop = imcrop(Image, CropCoord1);                          % rougly crop image
    ImRegister1 = imwarp(ImageInitCrop, tform, OutputView=outputView);  % transform image with affine
    ImRegister2 = imwarp(ImRegister1, tform2, OutputView=outputView);   % transform image with lwm
    ImStraightCrop = imcrop(ImRegister2, CropCoord2);                   % crop straightened image to pore space
    ResizedBF = imresize(ImStraightCrop, [(10000/1.62), (15000/1.62)]);   % resize image to correct dimensions
    ProcGrowthBF{1,i} = ResizedBF;

    'Image read, straightened, and cropped...' 
 
end

toc

%% Save growth phase image structure for future use

cd(VDG);
mkdir('Processed_Images');
PIG = [VDG '/Processed_Images'];
cd(PIG);
save('CroppedGrowthPhase.mat', 'ProcGrowthBF', '-v7.3');

'Processed image matrix saved...'

cd(OxidationD);

%% Register all brightfield images from mineralization phase of experiment

clearvars ims ImageName Image ImageInitCrop ImRegister1 ImRegister2 ImStraightCrop ResizedBF 
ims = dir('*.tif');
ims = ims(~startsWith({ims.name},'._'));
ImageName = {ims.name};
EndNum = length(ImageName);

tic

% Load all images, straighten, crop, store into structure
for i = 1:EndNum

    Image = imread(ImageName{1,i});                                         % read in image
    ImageInitCrop = imcrop(Image, CropCoord1);                              % rougly crop image
    ImRegister1 = imwarp(ImageInitCrop, tform, OutputView=outputView);      % transform image with affine
    ImRegister2 = imwarp(ImRegister1, tform2, OutputView=outputView);       % transform image with lwm
    ImStraightCrop = imcrop(ImRegister2, CropCoord2);                       % crop straightened image to pore space
    ResizedBF = imresize(ImStraightCrop, [(10000/1.62), (15000/1.62)]);     % resize image to correct dimensions
    ProcOxidationBF{i} = ResizedBF;

    'Image read, straightened, and cropped...'   
 
end

toc

%% Save structure of cropped images for future reference

cd(VDO);
mkdir('Processed_Images');
PIO = [VDO '/Processed_Images'];
cd(PIO);
save('CroppedOxidationPhase.mat', 'ProcOxidationBF', '-v7.3');

'Processed image matrix saved...'

cd(VDG);

% Save image registration variables
mkdir('Image_Registration');
IR = [VDG '/Image_Registration'];
cd(IR);
save('CroppingCoordinates_Step1.mat', 'CropCoord1');
save('CroppingCoordinates_Step2.mat', 'CropCoord2');
save('outputView.mat', 'outputView');
save('TForm.mat', 'tform');
save('TForm2.mat', 'tform2');
save('Moving_Points.mat', 'movingPoints');
save('Fixed_Points.mat', 'fixedPoints');
cd(VDG);

% Save grain mask
mkdir('Grain_Mask');
GM = [VDG '/Grain_Mask'];
cd(GM);
save('GrainMask.mat', 'GM_Final');
cd(VDG);

'Variables saved...'

%% Produce video from timelapse images to determine onset of stationary phase

close all;

% Link to directory for saving video
cd(VDO);
mkdir('Stationary_Identification');
SID = [VDO '/Stationary_Identification'];
cd(SID);

% Create labels for timepoints
TP_Growth = 0:1:15;   % growth phase
TP_Ox = 15:0.5:43.5;  % mineralization phase

Title = strcat('Stationary_Check.mp4');
video = VideoWriter(Title, 'MPEG-4');
video.FrameRate = 1;
open(video);

% Add timepoints to images and view
tic

% Growth phase
for ii = 1:length(TP_Growth)
    ImLab = ProcGrowthBF{ii};
    TimeLab = strcat('Time (h):', {' '}, num2str(TP_Growth(ii)));
    
    I = ImLab;

    I = imshow(ImLab);
    hold on;
    text(8000, 300, TimeLab{1},'FontSize', 14, 'FontWeight', 'bold', 'Color', 'k');
    hold off;
    
    F = getframe(gcf);
    writeVideo(video,F); 
    
    'Image labeled and frame acquired...'
    
end

% Mineralization phase
for ii = 1:length(TP_Ox)
    ImLab = ProcOxidationBF{ii};
    TimeLab = strcat('Time (h):', {' '}, num2str(TP_Ox(ii)));
    
    I = ImLab;

    I = imshow(ImLab);
    hold on;
    text(8000, 300, TimeLab{1},'FontSize', 14, 'FontWeight', 'bold', 'Color', 'k');
    hold off;
    
    F = getframe(gcf);
    writeVideo(video,F); 
    
    'Image labeled and frame acquired...'
    
end

toc

%% Identify and assign reference image manually (from video)

close all;

% Manually select reference image for each experiment based on motion
% observed during video
    % First timepoint at which biofilms are fixed in place and not changing
    % their footprint

% Update this manually! (Based on video)
RefNum = 8;

% Confirm that reference is a good fit for all subsequent images
% Overlay initial mask and final mask
    % Signal that only appears in initial mask is red
    % Signal that only appears in final mask is green
    % Signal that overlaps in both masks is yellow

BWInit = imcomplement(imbinarize(rgb2gray(ProcOxidationBF{1,1})));
BWRef = imcomplement(imbinarize(rgb2gray(ProcOxidationBF{1,RefNum})));
BWFin = imcomplement(imbinarize(rgb2gray(ProcOxidationBF{1,EndNum})));

IFOverlay = imfuse(BWInit, BWFin, 'falsecolor', 'ColorChannels', [1 2 0]);
drawnow;
f1 = figure(1); 
imshow(IFOverlay);
title('Overlay of initial mask and final mask');

% Overlay reference mask and final mask
    % Signal that only appears in reference mask is red
    % Signal that only appears in final mask is green
    % Signal that overlaps in both masks is yellow

IROverlay = imfuse(BWRef, BWFin, 'falsecolor', 'ColorChannels', [1 2 0]);
drawnow;
f2 = figure(2); 
imshow(IROverlay);
title('Overlay of reference mask and final mask');

%% Save remaining variables

cd(VDO);
mkdir('Image_Subtraction');
IS = [VDO '/Image_Subtraction'];
cd(IS);
save('ReferenceImageNumber.mat', 'RefNum');
save('EndingImageNumber.mat', 'EndNum');
cd(OxidationD);
