%% trackLoftedObjects Description of Code (DoC)
% Created by Charles T Pett
% Date created: 2022-08-02

% *** Information within a section in this format is essential to the *** %
% *** proper operation of that section.                               *** %

% This code reads in TIFF images from the clump lofting experiment,
% processes them, and analyzes the trajectory of clumps.

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read in Images, Select a Clump of Interest (COI), & Crop for Analysis
% *** Manual user input required.                                     *** %
function out = trackLoftedObjects(imDir)
% Look through image sequence, find clump lofting events, select a single
% COI to analyze, obtain suitable dimensions for the crop tool such that
% the COI remains in the image over the lifetime of its trajectory (but try
% to cut out as much as possible), extract start and end images for a given
% COI trajectory

% 'Current folder' needs to contain the folder we're grabbing the images
% from (i.e. we shouldn't be inside the image folder)
imFolderInfoMicro = dir(imDir);

% Reverse the direction of the axis such that the origin of the images is
% in the bottom left corner instead of the top left corner
% set(gca,'YDir','normal')

% Get rid of first few rows of imFolerInfo that have no image data
count = 1;
while imFolderInfoMicro(count).bytes == 0
    count = count + 1;
end

% Rewrite nonzero image data to new imFolerInfo
imFolderInfoMicro = imFolderInfoMicro(count:end,:);

% Count number of images in '0.2-0.3 mm silicon microspheres' directory
nMicro = numel(imFolderInfoMicro);

answer1 = 2;
while answer1 == 2

    % Saturate raw image to see the information and select start frame for
    % clump of interest
    prompt1 = 'Which frame in the sequence do you want to begin with? ';
    imStart = input(prompt1);
    rawStart = imadjust(flipud(imread(imFolderInfoMicro(imStart).name)));
    figure; subplot(1,2,1); imshow(rawStart); title('Start image')
    clc

    % Saturate raw image to see the information and select end frame for
    % clump of interest
    prompt2 = 'Which frame in the sequence do you want to end with? ';
    imEnd = input(prompt2);
    rawEnd = imadjust(flipud(imread(imFolderInfoMicro(imEnd).name)));
    subplot(1,2,2); imshow(rawEnd); title('End image')
    clc

    % Crop adjusted image in order to obtain the field of view that includes
    % the lifetime of the clump (i.e. its entire trajectory)
    % Coordinates are 4-element vector: [xmin ymin width_along_x height_along_y]
    prompt3 = 'What are the coordinates of the rectangular region to crop? ';
    cropDim = input(prompt3);
    figure; subplot(1,2,1); imshow(imcrop(rawStart,cropDim)); title('Start image')
    subplot(1,2,2); imshow(imcrop(rawEnd,cropDim)); title('End image')
    clc

    prompt4 = ['Are you happy with your crop?' ...
        ' [1=yes, 2=no] '];
    answer1 = input(prompt4);
    clc

    close all

end

% Output crop stats
cropStats = [nMicro,imStart,imEnd,cropDim];

% Move the code along when ready. Make sure you press enter while in the
% Command Window
input('Press ''Enter'' to continue...','s');
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SELECT THRESHOLD REGIME, CLEAR BORDERS, BINARIZE, INITIALIZE COI CENTROID
% *** Manual user input required.                                     *** %
% *** This section and those that follow are specific to a single COI *** %

close all

% Stop frame in for loop for single clump tracking
imStop = nMicro - imEnd - 1;

answer3 = 2;
while answer3 == 2

    close all

    prompt5 = ['Do you want to continue with the default ''imadjust'' or would ' newline 'you like to edit the threshold settings? [1=default, 2=edit settings] '];
    answer2 = input(prompt5);
    clc;

    if answer2 == 1
        % Saturate raw image to see the information. Toggle/uncomment between the
        % default and advanced 'imadjust' function options
        raw = imadjust(flipud(imread(imFolderInfoMicro(imStart).name)));
        imadjustStats = {'default'};
    else
        % 'Imadjust' thresholding settings if default does not suffice
        prompt6 = 'Low_in value (between 0 and 1)? '; low_in = input(prompt6);
        prompt7 = 'High_in value (between 0 and 1)? '; high_in = input(prompt7);
        prompt8 = 'Low_out value (between 0 and 1)? '; low_out = input(prompt8);
        prompt9 = 'High_out value (between 0 and 1)? '; high_out = input(prompt9);
        prompt10 = 'Gamma value (>= 0)? '; gamma = input(prompt10);

        raw = imadjust(flipud(imread(imFolderInfoMicro(imStart).name)),...
            [low_in high_in],[low_out high_out],gamma);
        imadjustStats = {'edit',low_in,high_in,low_out,high_out,gamma};
    end

    % Crop image to include COI over the lifetime of its trajectory. This is
    % the crop dimensions that will be utilized in the image sequence for this
    % COI
    crop = imcrop(raw,cropDim);
    subplot(1,2,1); imshow(crop)

    % Clear borders between COI and the background pixels (i.e. separate COI
    % from the background)
    im = imclearborder(crop);
    % Binarize the image (0 = black, 1 = white)
    im = imbinarize(im);
    % Fill in dim lit sections of the COI that turned into holes after
    % binarizing
    im = imfill(im,'hole');
    % Keep only clumps with the largest areas, which should include the COI
    % over the lifetime of its trajectory. Decide how many to keep.
    prompt11 = 'How many of the largest clump areas would you like to keep? ';
    areas2keep = input(prompt11);
    clc
    im = bwareafilt(im,areas2keep);

    % Draw boundaries around each clump (discluding holes) and visualize the
    % product of the user inputs selected above
    [B,L] = bwboundaries(im,'noholes');
    labs = label2rgb(L, @jet, [.5 .5 .5]);
    subplot(1,2,2); imshow(labs)
    hold on
    for k = 1:length(B)
        boundary = B{k};
        plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)
    end

    % Better visualization than the 'subplot' output
    % *** Do not use for analysis because the pixel grid is not true      *** %
    figure
    montage({crop, labs},'Size',[1,2]);
    title('Montage, not for analysis')

    prompt12 = ['Hopefully you have kept enough areas to identify the COI! If so, ' newline ...
        'would you like to start tracking? [1=yes, 2=no] '];
    answer3 = input(prompt12);
    clc

end

% *** SUPER IMPORTANT *** %
% *** Once satisfied with the result, the above image will be used to *** %
% *** select the coordinates of the COI and determine which of the    *** %
% *** clumps in the 'stats0' variable is the COI                      *** %

% *** In the image output window, use data tips to approximate the    *** %
% *** centroid of the COI (just click the middle of the COI) in order *** %
% *** to display the centroid coordinates. Search for the coordinates *** %
% *** in the 'stats0' structure below to acquire the initial centroid *** %
% *** to feed into the tracker *** %

% *** Each row in 'stats0' contains properties for each clump in the  *** %
% *** cropped image.                                                  *** %
% *** Open up the 'stats0' structure from the workspace and search    *** %
% *** for the Centroid coordinates that best match the coordinates    *** %
% *** displayed by the data tip in the subplot above. Write down the  *** %
% *** row #. This row # will be saved as the variable 'centRowi' and  *** %
% *** used in the following section.                                  *** %
stats0 = regionprops('struct',im,'all');
stats0.Centroid
% Manually selected clump in stats to initialize centroid tracking
prompt13 = ['Based on the estimated centroid from the data tip in the subplot,' newline 'which output in the Command Window is closest to the COI centroid ' newline '(i.g. For second ans, type 2)?' newline 'This coordinate will initialize the tracking sequence. '];
centRowi = input(prompt13);
clc
% *** End SUPER IMPORTANT :)                                          *** %

% Move the code along when ready. Make sure you press enter while in the
% Command Window
input('Press ''Enter'' to set up tracking...','s');
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Start section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TRACK THE CENTROID OF THE COI OVER THE ASSOCIATED IMAGE SPAN
% *** Automated section - little to no user input required            *** %
% This section is responsible for tracking the COI and storing its
% properties in a cell structure

close all

% Grab the centroid of the COI in order to initialize the centroid tracker.
% 'centRowi' was determined in the previous section
statsi = stats0(centRowi).Centroid;

% Parse out the X and Y component of the initial centroid coordinates in
% the image.
% Note: (x,y)_crop ~= (x,y)_raw, but (x2-x1,y2-y1) equivalent in both
% frames
centXi = statsi(1);
centYi = statsi(2);

% *** Iterations of manual input required to set the tolerances that  *** %
% *** achieve the most accurate (smoothest) trajectory                *** %
% Searching for the centroid in the next image based on the centroid of the
% COI in the previous image, within some tolerance. This works because the
% clumps don't move much image to image.
prompt14 = {'Centroid X-coord tolerance:','Centroid Y-coord tolerance:'};
% Xtol = str2double(inputdlg(prompt14(1)));
% Ytol = str2double(inputdlg(prompt14(2)));
% XYtol = [Xtol,Ytol];
XYtol = str2double(inputdlg(prompt14));
tols = {areas2keep centRowi XYtol};

% Loop through the user-selected image span associated with the COI
for i = 1:nMicro-imStart-imStop

    % Read in and 'imadjust' the first image in which the COI is detectable
    raw = imadjust(flipud(imread(imFolderInfoMicro(i+imStart-1).name)));

    % Crop the current raw image to the dimensions previously selected.
    % Same crop dimensions for all images in the span because crop should
    % contain the lifetime of the trajectory
    crop = imcrop(raw,cropDim);

    % Clear borders between COI and the background pixels (i.e. separate COI
    % from the background)
    im = imclearborder(crop);
    % Binarize the image (0 = black, 1 = white)
    im = imbinarize(im);
    % Fill in dim lit sections of the COI that turned into holes after
    % binarizing
    im = imfill(im,'hole');
    % Keep only clumps with the largest areas, which should include the COI
    % over the lifetime of its trajectory. Decide how many to keep.
    im = bwareafilt(im,areas2keep);
    % *** 'areas2keep', which was selected in the previous section    *** %
    % *** based on ONE image, may need to be increased or decreased   *** %
    % *** if the number of clumps entering the frame changes later in *** %
    % *** the image sequence. If not enough areas have been kept such *** %
    % *** that the COI is discluded from the image, then no centroid  *** %
    % *** within tolerance will exist and the loop will abort, most   *** %
    % *** likely with the following message:                          *** %
    % *** Undefined function 'stats' for input arguments of type      *** %
    % *** 'double'. Error in trackSingleClumps (line 231)             *** %
    % *** centXi = stats(i).Centroid(1);                              *** %

    % *** If too many areas have been kept and the tolerance is too   *** %
    % *** loose, then the trajectory will be spikey and unsmooth      *** %
    % *** because clumps that don't have to be included and in the    *** %
    % *** vicinity of the COI will be tracked and add noise.          *** %

    % Let's show the 'im' with the COI centroid coordinate plot overlaid
    % and update this after each iteration in the loop in order to draw the
    % trajectory over the lifetime of the COI. Also a good visual to
    % reference for fine tuning of the user-inputted parameters.
    imshow(im)
    hold on
    drawnow
    % Nothing to plot at i = 1 yet
    if i ~= 1
        % Plots the centroid of the COI at each image
        plot(centX,centY,'LineWidth',1.5)
        pause(0.1)
    end
    hold off

    % Obtain updated region properties of the COI in the current image
    statsNew = regionprops('struct',im,'all');

    % Multiple the current COI centroid coordinates by the tolerance to get
    % a fraction of the current COI centroid
    centXtol = XYtol(1)*centXi;
    centYtol = XYtol(2)*centYi;

    % 'statsNew' will contain region properties on all the clumps in the
    % image, including the COI
    for j = 1:length(statsNew)

        % Take the absolute value of the difference between the new
        % centroid coordinates and the old centroid coordinates of the COI
        centXdiff = abs(statsNew(j).Centroid(1) - centXi);
        centYdiff = abs(statsNew(j).Centroid(2) - centYi);

        % If the difference between the new, prospective centroid and the
        % old centroid of the COI is within the allowable tolerance, then
        % I've found the COI in the new image! If so, let's store its
        % properties and save its centroid coordinates in a separate file
        % to use in plotting the trajectory in the code upstairs.
        if centXdiff < centXtol && centYdiff < centYtol
            stats(i,1) = statsNew(j);
            centX(i,1) = stats(i).Centroid(1);
            centY(i,1) = stats(i).Centroid(2);
        end
        % If COI not found, don't save anything.
    end

    % Update the current centroid coordinates to be the initial ones in the
    % next iteration for the next image
    centXi = stats(i).Centroid(1);
    centYi = stats(i).Centroid(2);

end

% Output parameters and statistics to cell structure
out.cropParams = cropStats; % Dimensions of the cropped region
out.imadjustParams = imadjustStats; % Image adjustment parameters/thresholding
out.tolerances = tols; % Centroid tracking tolerance
out.initialStats = stats0; % Initial conditions of tracked clump
out.Stats = stats; % Parameters of tracked clump over entire trajectory
out.Centroid = [centX centY]; % Centroid positions of tracked clump over trajectory

disp('Tracking complete!')

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
