function [legEndpointsStore, frameFlags, centroidCoordinates, headDirectionCoordinates] = AutomaticTrackingFunction(movName,...
    methodFlag,thresholdingFactorForCentroid, pixelAcceptanceModifierThresholding, smallComponentSizeForCentroid,...
    minDistanceThresholdFactor, thresholdingFactorForTracking, lowConfidenceQuantile,lowConfidenceCutoffDistance,targetNumLegs,...
    showMovieFlag,moviePauseTime)
% AutomaticTrackingFunction loads the movie (assumed cropped to frames of interest),
% finds fly centroid/orientation, and performs leg tracking with one of two dynamic thresholding algorithms.
%
%%%%%%% Arguments:
% movName [string]: name of movie to consider
% methodFlag [1 OR 2]: which method to use for leg detection.
    % 1 is the original distance-based, 2 is the constrained search.
    % Note that if method 2 comes up with too few regions for a frame, it will revert
    % to method 2 for that frame.
%
%
% ALL FURTHER ARGUMENTS CAN BE SET TO 'default' %
%
%
% thresholdingFactorForCentroid [double between 0 and 1 OR 'default'];  
    % should be between 0 and 1 to decrease the black/white threshold value (remember, white pixels = 1).
    % We will multiply this by the Otsu threshold level vector to take out the
    % legs/wings, giving a more consistent centroid (since the legs are moving not necessarily
    % symmetrically). However, we can't make it too small, since then it starts to take out 
    % too many parts of the body as well. Default value is .7, though a wide range works.
        % default = .7
% pixelAcceptanceModifierThresholding [double generally >= 1 (should be around 1) OR 'default']
    % pixelAcceptanceModifierThresholding*max(number of pixels in body with minimum thresholding)
    % gives the number of pixels in the fly connected component we accept. Defaults to 1.1 - if that
    % gives too much noise (random pixels that aren't thresholded out get "stuck" to fly and the
    % image is accepted) feel free to decrease this - even 1 is good since the second term is
    % max(number of pixels in body with minimum thresholding), which is typically quite a bit greater
    % than the median anyway. Can decrease below 1, but there shouldn't be a reason to. If too low,
    % everything will be rejected and the minimum thresholded image will end up getting used (less
    % accuracy on legs).
        % default = 1.1
% smallComponentSizeForCentroid [integer >= 0]
    % remove all connected components with smallComponentSize pixels or fewer (random pixels 
    % that remain even after thresholding). 5 seems to work as well as 10 or 20, but no danger
    % in setting set it high (even over a hundred, though there is no benefit to it either)
    % to remove everything that is not fly body or directly connected since fly body is huge
    % and anything "smaller" than it and disconnected should get removed)
    % HAVE to do this, otherwise inbuilt centroid jumps around to the connected components/
    % Default value is 300;
        % default = 300
% minDistanceThresholdFactor [double between 0 and 1 OR 'default']: first index where derivative of 
    % number of distances that far away from the center line flattens out sufficiently
    % (where we should start counting a "low" ("on body") distance). 10% seems to work well (deriv
    % becomes factor*deriv(t0)). Can set this higher, but then in general might find fewer than the
    % target number of legs (but those will generally be very accurate) -  this might be useful if 
    % we want high confidence in a few legs and have a human just tag e.g. the last 1 leg per frame.
    % Default value is .1.
        % default = .1
% thresholdingFactorForTracking [double >= 1 OR 'default']: increases Otsu thresholding so we 
    % can clearly see the legs, but can't be so high that we get too many extraneous regions.
    % to remove the extraneous pixels from the higher thresholding, we can use 
    % max(numBlackPixelsBody) (already calculated in function for fly not including legs) as a very 
    % safe lower bound on the  size of the region and remove the connected components that have 
    % size less than that. Default value is 1.3.
        % default = 1.3
% lowConfidenceQuantile [double between 0 and 1 OR 'default']: 
    % distance removal for 20% of "clearest candidate" frames potentially
    % does not require intervention at all. This is for marking "low" confidence frames (the vast 
    % majority of those will be correct, but just in case).
    % I found that a better alternative is to set the low confidence flag to be a min of this number 
    % and e.g. 30, since almost all body pixels are < 15 away from the center line, so
    % when the distance is 30 away from the body we are still very confident [this especially
    % happens in some of the post-amputation videos, where the .8 quintile is like 30 and we are
    % confident in all frames - the flies just tend to splay their legs more in some of the videos].
    % Note that this is not the case with dynamic thresholding, since now we have the whole leg
    % being tracked more accurately and this can be many pixels away, while some closer connected
    % components double up. But this is still a reasonable proxy.
    % See below for this "hard" suggested lower bound cutoff.
        % default = .8;
% lowConfidenceCutoffDistance [integer >= 15 recommended OR 'default']: "hard" distance bound
    % in min (we are very confident that this is a leg endpoint if the distance is greater than this
    % from the center of the body). See above.
        % default = 23;
% targetNumLegs [integer between 1 and 6 OR 'default']: target number of legs to attempt to find.
    % If default is set, "preamputation" in move title sets to 6 legs, "postamputation" in movie
    % title sets to 5 legs
% showMovieFlag ['y' OR 'n' OR 'default']: show tracking "live".
    % default = 'y'
% moviePauseTime = [double >= 0 OR 'default']: how long to pause between showing movie frames "live"
    % default = .05



    
if isequal(thresholdingFactorForCentroid,'default')
    thresholdingFactorForCentroid = .7;
end

if isequal(pixelAcceptanceModifierThresholding,'default')
    pixelAcceptanceModifierThresholding = 1.1;
end

if isequal(smallComponentSizeForCentroid,'default')
    smallComponentSizeForCentroid = 300;
end

if isequal(minDistanceThresholdFactor,'default')
    minDistanceThresholdFactor = .1;
end

if isequal(thresholdingFactorForTracking,'default')
    thresholdingFactorForTracking = 1.3;
end

           
if isequal(lowConfidenceQuantile,'default')
    lowConfidenceQuantile = .8;
end
            
if isequal(lowConfidenceCutoffDistance,'default')
    lowConfidenceCutoffDistance = 23; %#ok<*NASGU>
end

    
if isequal(targetNumLegs,'default')
    if not(isempty(strfind(movName,'preamputation'))) % if file name includes "preamputation"
        targetNumLegs = 6; % find 6 legs
    elseif not(isempty(strfind(movName,'postamputation'))) % if file name includes "postamputation"
        targetNumLegs = 5; % find 5 legs
    end
end

if isequal(showMovieFlag,'default')
    showMovieFlag = 'y';
end

if isequal(moviePauseTime,'default')
    moviePauseTime = .05;
end


% Create Temporary Working Directory
workingDir = tempname;
mkdir(workingDir);
mkdir(workingDir,'images');

% Load Original Movie And Create Sorted Images From Frames
movOriginal = VideoReader(movName);
numFrames = movOriginal.NumberOfFrames; %#ok<VIDREAD>
movHeight = movOriginal.Height;
movWidth = movOriginal.Width;

for i = 1:numFrames
    img = read(movOriginal,i); %#ok<VIDREAD>
    img = img(:,:,2); % take only green channel (camera was set to process in green-scale)
    imwrite(img,fullfile(workingDir,'images',sprintf('img%d.jpg',i)));
end
    imwrite(img,fullfile(workingDir,'images',sprintf('img%d.jpg',i)));
imageNames = dir(fullfile(workingDir,'images','*.jpg'));
imageNames = {imageNames.name}';

imageStrings = regexp([imageNames{:}],'(\d*)','match');
imageNumbers = str2double(imageStrings);

[~,sortedIndices] = sort(imageNumbers);
sortedImageNames = imageNames(sortedIndices);



OtsuMethodLevelVector = zeros(1,numFrames);
for i = 1:numFrames
    img = imread(fullfile(workingDir,'images',sortedImageNames{i}));
    OtsuMethodLevelVector(i) = graythresh(img);
end

maxThreshold = median(OtsuMethodLevelVector);
minThreshold = thresholdingFactorForCentroid*maxThreshold;


% Fly Body Tracking
inbuiltCentroidXCoordinateVector = zeros(1,numFrames);
inbuiltCentroidYCoordinateVector = zeros(1,numFrames);
inbuiltOrientationAngleVector = zeros(1,numFrames);

xLineEndpointVector = zeros(2,numFrames);
yLineEndpointVector = zeros(2,numFrames);
numPixelsBodyVector = zeros(1,numFrames);

flyCentralLinePoints = cell(1,numFrames);
potentialHeadDirectionCoordinates_1 = zeros(numFrames,2);
potentialHeadDirectionCoordinates_2 = zeros(numFrames,2);

for i = 1:numFrames
    img = imread(fullfile(workingDir,'images',sortedImageNames{i}));
    imgForCentroid = im2bw(img,minThreshold); % note: black pixel = 0, white pixel = 1
    imgForCentroid = bwareaopen(~imgForCentroid,smallComponentSizeForCentroid);
    numPixelsBodyVector(i) = length(find(imgForCentroid));
    
    imageShapePropertiesForCentroid = regionprops(imgForCentroid, 'Centroid','Orientation','MajorAxisLength','MinorAxisLength');
    inbuiltCentroidXCoordinateVector(i) = imageShapePropertiesForCentroid.Centroid(1);
    inbuiltCentroidYCoordinateVector(i) = imageShapePropertiesForCentroid.Centroid(2);
    inbuiltOrientationAngle = imageShapePropertiesForCentroid.Orientation;
    inbuiltMajorAxisLength = imageShapePropertiesForCentroid.MajorAxisLength;
    inbuiltMinorAxisLength = imageShapePropertiesForCentroid.MinorAxisLength; 
    inbuiltOrientationAngleVector(i) = inbuiltOrientationAngle;    

    % define central axis and orientation
        % IMPORTANT NOTE: we need to have a "continuous" orientation (so the nominal "head" is always the head, and the nominal "tail" is always the tail, and we
        % choose which is actually the head or tail later). Since the ellipse is symmetric and Matlab only gives an angle between 90 and -90 degrees from the x-axis
        % (which flips head-tail when the fly crosses the vertical axis), we arbitrarily choose the "head" in frame 1 and then propose the two locations (head, tail) for
        % subsequent frames and choose the closest one to the previous "head" to be the current "head"
 
    headXCoordinate_proposal = inbuiltCentroidXCoordinateVector(i) - inbuiltMajorAxisLength/2 * cosd(inbuiltOrientationAngle);
    headYCoordinate_proposal = inbuiltCentroidYCoordinateVector(i) - inbuiltMajorAxisLength/2 * sind(-inbuiltOrientationAngle);
    tailXCoordinate_proposal = inbuiltCentroidXCoordinateVector(i) + inbuiltMajorAxisLength/2 * cosd(inbuiltOrientationAngle);
    tailYCoordinate_proposal = inbuiltCentroidYCoordinateVector(i) + inbuiltMajorAxisLength/2 * sind(-inbuiltOrientationAngle); 
    
    if i == 1 % choose which angle represents the "head" and which the "tail" of the ellipse, arbitrary for first frame
        xLineEndpointVector(1,i) = headXCoordinate_proposal;
        yLineEndpointVector(1,i) = headYCoordinate_proposal;
        xLineEndpointVector(2,i) = tailXCoordinate_proposal;
        yLineEndpointVector(2,i) = tailYCoordinate_proposal;
    else % i > 1
        distHeadProposalToOldHead = sqrt((xLineEndpointVector(1,i-1)-headXCoordinate_proposal)^2+(yLineEndpointVector(1,i-1)-headYCoordinate_proposal)^2);
        distTailProposalToOldHead = sqrt((xLineEndpointVector(1,i-1)-tailXCoordinate_proposal)^2+(yLineEndpointVector(1,i-1)-tailYCoordinate_proposal)^2);
        if distHeadProposalToOldHead < distTailProposalToOldHead
            % if distance between "head" proposal and last "head" position is less than  distance between tail proposal and last "head" position, accept it as new "head" position
            xLineEndpointVector(1,i) = headXCoordinate_proposal;
            yLineEndpointVector(1,i) = headYCoordinate_proposal;
            xLineEndpointVector(2,i) = tailXCoordinate_proposal;
            yLineEndpointVector(2,i) = tailYCoordinate_proposal;
        else
            xLineEndpointVector(1,i) = tailXCoordinate_proposal;
            yLineEndpointVector(1,i) = tailYCoordinate_proposal;
            xLineEndpointVector(2,i) = headXCoordinate_proposal;
            yLineEndpointVector(2,i) = headYCoordinate_proposal;
        end
    end
        
    
    
    xLinePoints = round(linspace(xLineEndpointVector(1,i),xLineEndpointVector(2,i),1000));
    yLinePoints = round(linspace(yLineEndpointVector(1,i),yLineEndpointVector(2,i),1000));
    tempLinePointsMatrix = [xLinePoints',yLinePoints'];
    tempLinePointsMatrix = unique(tempLinePointsMatrix,'rows','stable');
    xLinePoints = tempLinePointsMatrix(:,1)';
    yLinePoints = tempLinePointsMatrix(:,2)';
    linePointRemovalIndex = []; % remove points that lie outside image (if line is too long)

    for j = 1:length(xLinePoints)
        if xLinePoints(j) < 1 || yLinePoints(j) < 1 || xLinePoints(j) > movWidth || yLinePoints(j) > movHeight
            linePointRemovalIndex = [linePointRemovalIndex,j]; 
        end
    end
    
    xLinePoints(linePointRemovalIndex) = [];
    yLinePoints(linePointRemovalIndex) = [];
    flyCentralLinePoints{i} = [xLinePoints',yLinePoints'];
    
    potentialHeadDirectionCoordinates_1(i,1) = xLinePoints(1);
    potentialHeadDirectionCoordinates_1(i,2) = yLinePoints(1);
    potentialHeadDirectionCoordinates_2(i,1) = xLinePoints(end);
    potentialHeadDirectionCoordinates_2(i,2) = yLinePoints(end);

end

    % Get correct head and tail coordinate
    
centroidCoordinates = round([inbuiltCentroidXCoordinateVector',inbuiltCentroidYCoordinateVector']);
distToEnd1 = zeros(1,numFrames-1);
distToEnd2 = zeros(1,numFrames-1);
for i = 2:numFrames
    distToEnd1(i-1) = sum((centroidCoordinates(i,:) - potentialHeadDirectionCoordinates_1(i-1,:)).^2);
    distToEnd2(i-1) = sum((centroidCoordinates(i,:) - potentialHeadDirectionCoordinates_2(i-1,:)).^2);
end

if length(find(distToEnd1<distToEnd2)) > length(find(distToEnd2<distToEnd1)) % if distToEnd1 is usually smaller than distToEnd2
    headDirectionCoordinates = potentialHeadDirectionCoordinates_1; % heading towards End 1 of body axis
    tailDirectionCoordinates = potentialHeadDirectionCoordinates_2;
else
    headDirectionCoordinates = potentialHeadDirectionCoordinates_2; % heading towards End 2 of body axis
    tailDirectionCoordinates = potentialHeadDirectionCoordinates_1;
end



% Leg Tracking with Dynamic Thresholding Algorithm (two methods)
    
% Dynamic Thresholding Algorithm Overview: use minimum threshold to find "true centroid". Then, start with
    % maximum threshold and find the connected component where the centroid lies. Remove all components
    % without the centroid. If the number of pixels left is acceptable (not too high or too low, though it should never be too low by
    % definition), use that threshold for that frame and find the legs. Otherwise, drop the threshold and try again, until minimum threshold is hit.

removalBodyPixelsDistanceHighToLowVector = zeros(1,numFrames); % distance at which the tagging finishes (used for giving confidence bin to frame)
frameFlags = zeros(1,numFrames);

    % FRAME FLAGS
    % 0 --> frame might not require intervention (high confidence)
    % 1 --> frame should be looked at to verify
    % 2 --> frame REQUIRES intervention (can't find enough endpoints/found too many when going from high distance to low distance removal)

legEndpointsStore = cell(1,numFrames);
smallComponentSizeForTracking = min(numPixelsBodyVector);

for i = 1:numFrames
    img = imread(fullfile(workingDir,'images',sortedImageNames{i}));
    
    % Dynamic Thresholding
        % NOTE: Takes about 12 seconds per 60 frames. To decrease time spent, set fewer steps to try for thresholds (decrease numDynamicThresholdSteps)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    numDynamicThresholdSteps = 15; % maximum number of thresholds to try, starting from max and going to min
    thresholdAttemptVector = linspace(maxThreshold,minThreshold,numDynamicThresholdSteps);
    if maxThreshold < .7 % if maximum is too small (depends on movie), put in larger function argument
        thresholdAttemptVector = linspace(.7,minThreshold,numDynamicThresholdSteps);
    end
    acceptThresholdFlag = 0;
    thresholdStep = 0;
    while acceptThresholdFlag == 0 && thresholdStep < numDynamicThresholdSteps
        thresholdStep = thresholdStep + 1;
        tempThreshold = thresholdAttemptVector(thresholdStep);
        % set connectivity type
        if tempThreshold < numDynamicThresholdSteps/2 % for first half (highest thresholds)
            connectivityForTracking = 4; % consider 4-connected components (weakest standard connectivity) are "connected"
                % to prevent "noisy" pixels that are around from connecting to legs
        else % lower thresholds
            connectivityForTracking = 8; % stronger connectivity to prevent the leg ends from disappearing
        end

        imgForTracking = ~im2bw(img,tempThreshold); % white is signal (fly)
        imgForTracking = bwareaopen(imgForTracking,smallComponentSizeForTracking,connectivityForTracking); % decrease search by preliminary removal of small components
        tempConnectedComponents = bwconncomp(imgForTracking,connectivityForTracking);
        for k = 1:tempConnectedComponents.NumObjects
            tempPixelIdxList = tempConnectedComponents.PixelIdxList{k};
            [tempRowCoordinates, tempColumnCoordinates] = ind2sub([movHeight, movWidth],tempPixelIdxList);
            tempXYCoordinateVector = [tempColumnCoordinates,tempRowCoordinates];
            if isempty(intersect(centroidCoordinates(i,:),tempXYCoordinateVector,'rows')) % if this is the component without the fly
                imgForTracking(tempPixelIdxList) = 0; % remove it (set all those pixels to 0 = black)
            else % this is the component with the fly
                if maxThreshold >= .7
                    if length(tempPixelIdxList) <= pixelAcceptanceModifierThresholding*max(numPixelsBodyVector) % even for the higher threshold, this is a very generous bound - can play with it
                        acceptThresholdFlag = 1;
                    end
                else % if maximum is small, we give it a chance to be large, and we hardcode the number of acceptable pixels to be reasonable but not too generous (depends on data!)
                    if length(tempPixelIdxList) <= 2300
                        acceptThresholdFlag = 1;
                    end
                end
            end
        end        
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
        % Run region props again and overwrite old thing since now we have a differently thresholded
        % fly due to the dynamic thresholding algorithm
    imageShapePropertiesForTracking = regionprops(imgForTracking, 'Centroid','Orientation','MajorAxisLength','MinorAxisLength');
    inbuiltCentroidXCoordinateVector(i) = imageShapePropertiesForTracking.Centroid(1);
    inbuiltCentroidYCoordinateVector(i) = imageShapePropertiesForTracking.Centroid(2);
    inbuiltOrientationAngle = imageShapePropertiesForTracking.Orientation;
    inbuiltMinorAxisLength = imageShapePropertiesForTracking.MinorAxisLength;

    

    xCent = round(inbuiltCentroidXCoordinateVector(i)); % round to integer (pixel)
    yCent = round(inbuiltCentroidYCoordinateVector(i)); % round to integer (pixel)
    
        % where on the body is the cutoff for the region where the legs are connected
    perpendicularLine1Point = [(xCent+headDirectionCoordinates(i,1))/2,(yCent+headDirectionCoordinates(i,2))/2];
	perpendicularLine2Point = [(xCent+(xCent+perpendicularLine1Point(1))/2)/2,(yCent+(yCent+perpendicularLine1Point(2))/2)/2];
   
    centalLineSide1XCoordinate = xCent - inbuiltMajorAxisLength/2 * cosd(inbuiltOrientationAngle);
    centalLineSide1YCoordinate = yCent - inbuiltMajorAxisLength/2 * sind(-inbuiltOrientationAngle); % remember: due to the way (0,0) in image is upper left corner, we need to have y --> -y 
    centalLineSide2XCoordinate = xCent + inbuiltMajorAxisLength/2 * cosd(inbuiltOrientationAngle);
    centalLineSide2YCoordinate = yCent + inbuiltMajorAxisLength/2 * sind(-inbuiltOrientationAngle); % remember: due to the way (0,0) in image is upper left corner, we need to have y --> -y 
    
    centralLineX = round(linspace(centalLineSide1XCoordinate,centalLineSide2XCoordinate,1000));
    centralLineY = round(linspace(centalLineSide1YCoordinate,centalLineSide2YCoordinate,1000)); 

    if methodFlag == 1
        %%%% PRIMARY METHOD - cut fly into 6 connected components (if  there are a few more disconnected
        %%%% ones, that's OK too). Choose a leg in each one in the same way as above (except now we
        %%%% *have* to choose a leg in each one). We then choose NumTargetLegs from the suggested
        %%%% points.
        perpAngle = inbuiltOrientationAngle + 90;

        line1side1XCoordinate = perpendicularLine1Point(1) - inbuiltMinorAxisLength/2 * cosd(perpAngle);
        line1side1YCoordinate = perpendicularLine1Point(2) - inbuiltMinorAxisLength/2 * sind(-perpAngle);
        line1side2XCoordinate = perpendicularLine1Point(1) + inbuiltMinorAxisLength/2 * cosd(perpAngle);
        line1side2YCoordinate = perpendicularLine1Point(2) + inbuiltMinorAxisLength/2 * sind(-perpAngle);
        perpendicularLine1X = round(linspace(line1side1XCoordinate,line1side2XCoordinate,1000));
        perpendicularLine1Y = round(linspace(line1side1YCoordinate,line1side2YCoordinate,1000));

        line2side1XCoordinate = perpendicularLine2Point(1) - inbuiltMinorAxisLength/2 * cosd(perpAngle);
        line2side1YCoordinate = perpendicularLine2Point(2) - inbuiltMinorAxisLength/2 * sind(-perpAngle);
        line2side2XCoordinate = perpendicularLine2Point(1) + inbuiltMinorAxisLength/2 * cosd(perpAngle);
        line2side2YCoordinate = perpendicularLine2Point(2) + inbuiltMinorAxisLength/2 * sind(-perpAngle);
        perpendicularLine2X = round(linspace(line2side1XCoordinate,line2side2XCoordinate,1000));
        perpendicularLine2Y = round(linspace(line2side1YCoordinate,line2side2YCoordinate,1000));

        imgForTracking2 = imgForTracking;
        lengthAdditionMainLine = 10;
        lengthAdditionPerpLine1 = 10;
        lengthAdditionPerpLine2 = 10;


        rowCoordsMaskCentralLine = [yCent,centralLineY];
        colCoordsMaskCentralLine = [xCent,centralLineX];
        removalIdxCentralLine = find(centralLineX < 1 | centralLineX > movWidth | centralLineY < 1 | centralLineY > movHeight);
        rowCoordsMaskCentralLine(removalIdxCentralLine+1) = [];
        colCoordsMaskCentralLine(removalIdxCentralLine+1) = [];
        
        geoDist = bwdistgeodesic(imgForTracking2,colCoordsMaskCentralLine,rowCoordsMaskCentralLine,'cityblock');
        imgForTracking2(geoDist<1) = 0;
        connCompFirstCut = bwconncomp(imgForTracking2,4);

            % If line didn't cut completely, extend it a bit. Note: 4 connectivity here since we want to encourage splitting the region
        if connCompFirstCut.NumObjects < 2
            centalLineSide1XCoordinate = xCent - (inbuiltMajorAxisLength+lengthAdditionMainLine)/2 * cosd(inbuiltOrientationAngle);
            centalLineSide1YCoordinate = yCent - (inbuiltMajorAxisLength+lengthAdditionMainLine)/2 * sind(-inbuiltOrientationAngle);
            centalLineSide2XCoordinate = xCent + (inbuiltMajorAxisLength+lengthAdditionMainLine)/2 * cosd(inbuiltOrientationAngle);
            centalLineSide2YCoordinate = yCent + (inbuiltMajorAxisLength+lengthAdditionMainLine)/2 * sind(-inbuiltOrientationAngle);
            centralLineX = round(linspace(centalLineSide1XCoordinate,centalLineSide2XCoordinate,1000));
            centralLineY = round(linspace(centalLineSide1YCoordinate,centalLineSide2YCoordinate,1000)); 

            rowCoordsMaskCentralLine = [yCent,centralLineY];
            colCoordsMaskCentralLine = [xCent,centralLineX];
            removalIdxCentralLine = find(centralLineX < 1 | centralLineX > movWidth | centralLineY < 1 | centralLineY > movHeight);
            rowCoordsMaskCentralLine(removalIdxCentralLine+1) = [];
            colCoordsMaskCentralLine(removalIdxCentralLine+1) = [];

            geoDistSecondAttempt = bwdistgeodesic(imgForTracking2,colCoordsMaskCentralLine,rowCoordsMaskCentralLine,'cityblock');
            imgForTracking2(geoDistSecondAttempt<1) = 0;
        end




        rowCoordsMaskPerpLine1 = perpendicularLine1Y;
        colCoordsMaskPerpLine1 = perpendicularLine1X;
        removalIdxPerpLine1 = find(perpendicularLine1X < 1 | perpendicularLine1X > movWidth | perpendicularLine1Y < 1 | perpendicularLine1Y > movHeight);
        rowCoordsMaskPerpLine1(removalIdxPerpLine1) = [];
        colCoordsMaskPerpLine1(removalIdxPerpLine1) = [];
        geoDistPerpLine1 = bwdistgeodesic(imgForTracking2,colCoordsMaskPerpLine1,rowCoordsMaskPerpLine1,'cityblock');
        imgForTracking2(geoDistPerpLine1<1) = 0;

        connCompSecondCut = bwconncomp(imgForTracking2,4);
        if connCompSecondCut.NumObjects < 4;
            line1side1XCoordinate = perpendicularLine1Point(1) - (inbuiltMinorAxisLength+lengthAdditionPerpLine1)/2 * cosd(perpAngle);
            line1side1YCoordinate = perpendicularLine1Point(2) - (inbuiltMinorAxisLength+lengthAdditionPerpLine1)/2 * sind(-perpAngle);
            line1side2XCoordinate = perpendicularLine1Point(1) + (inbuiltMinorAxisLength+lengthAdditionPerpLine1)/2 * cosd(perpAngle);
            line1side2YCoordinate = perpendicularLine1Point(2) + (inbuiltMinorAxisLength+lengthAdditionPerpLine1)/2 * sind(-perpAngle);
            perpendicularLine1X = round(linspace(line1side1XCoordinate,line1side2XCoordinate,1000));
            perpendicularLine1Y = round(linspace(line1side1YCoordinate,line1side2YCoordinate,1000));

            rowCoordsMaskPerpLine1 = perpendicularLine1Y;
            colCoordsMaskPerpLine1 = perpendicularLine1X;
            removalIdxPerpLine1 = find(perpendicularLine1X < 1 | perpendicularLine1X > movWidth | perpendicularLine1Y < 1 | perpendicularLine1Y > movHeight);
            rowCoordsMaskPerpLine1(removalIdxPerpLine1) = [];
            colCoordsMaskPerpLine1(removalIdxPerpLine1) = [];
            geoDistPerpLine1SecondAttempt = bwdistgeodesic(imgForTracking2,colCoordsMaskPerpLine1,rowCoordsMaskPerpLine1,'cityblock');
            imgForTracking2(geoDistPerpLine1SecondAttempt<1) = 0;
        end

        
        
        rowCoordsMaskPerpLine2 = perpendicularLine2Y;
        colCoordsMaskPerpLine2 = perpendicularLine2X;
        removalIdxPerpLine2 = find(perpendicularLine2X < 1 | perpendicularLine2X > movWidth | perpendicularLine2Y < 1 | perpendicularLine2Y > movHeight);
        rowCoordsMaskPerpLine2(removalIdxPerpLine2) = [];
        colCoordsMaskPerpLine2(removalIdxPerpLine2) = [];
        geoDistPerpLine2 = bwdistgeodesic(imgForTracking2,colCoordsMaskPerpLine2,rowCoordsMaskPerpLine2,'cityblock');
        imgForTracking2(geoDistPerpLine2<1) = 0;

        connCompThirdCut = bwconncomp(imgForTracking2,4);
        if connCompThirdCut.NumObjects < 6
            line2side1XCoordinate = perpendicularLine2Point(1) - (inbuiltMinorAxisLength+lengthAdditionPerpLine2)/2 * cosd(perpAngle);
            line2side1YCoordinate = perpendicularLine2Point(2) - (inbuiltMinorAxisLength+lengthAdditionPerpLine2)/2 * sind(-perpAngle);  
            line2side2XCoordinate = perpendicularLine2Point(1) + (inbuiltMinorAxisLength+lengthAdditionPerpLine2)/2 * cosd(perpAngle);
            line2side2YCoordinate = perpendicularLine2Point(2) + (inbuiltMinorAxisLength+lengthAdditionPerpLine2)/2 * sind(-perpAngle);  
            perpendicularLine2X = round(linspace(line2side1XCoordinate,line2side2XCoordinate,1000));
            perpendicularLine2Y = round(linspace(line2side1YCoordinate,line2side2YCoordinate,1000));

            rowCoordsMaskPerpLine2 = perpendicularLine2Y;
            colCoordsMaskPerpLine2 = perpendicularLine2X;
            removalIdxPerpLine2 = find(perpendicularLine2X < 1 | perpendicularLine2X > movWidth | perpendicularLine2Y < 1 | perpendicularLine2Y > movHeight);
            rowCoordsMaskPerpLine2(removalIdxPerpLine2) = [];  
            colCoordsMaskPerpLine2(removalIdxPerpLine2) = [];
            geoDistPerpLine2SecondAttempt = bwdistgeodesic(imgForTracking2,colCoordsMaskPerpLine2,rowCoordsMaskPerpLine2,'cityblock');
            imgForTracking2(geoDistPerpLine2SecondAttempt<1) = 0;
        end

        AllConnectedComponents = bwconncomp(imgForTracking2,4);
        legEndpointPixelCoordinates = zeros(targetNumLegs,2);
        proposedLegEndpoints = [];
        proposedLegEndpointsDistancedFromCentralLine = [];
        if AllConnectedComponents.NumObjects >= 6 % otherwise, we will use old method
            for j = 1:AllConnectedComponents.NumObjects
                tempImg = zeros(movHeight,movWidth);
                tempImg(geoDist<2)=1; % put central line back so we can use it as a mask
                tempImg(AllConnectedComponents.PixelIdxList{j}) = 1; % consider only that connected component
                tempImg = logical(tempImg);
                pixelIdxLegs = find(tempImg==1); % all "on" pixels in image

                tempGeoDist = bwdistgeodesic(tempImg,colCoordsMaskCentralLine,rowCoordsMaskCentralLine,'cityblock');
                [~,descendingSortLegDistanceIdx] = sort(tempGeoDist(pixelIdxLegs),'descend');
                legEndpointPixel = pixelIdxLegs(descendingSortLegDistanceIdx(1));
                [proposedLegEndpoints(j,2),proposedLegEndpoints(j,1)] = ind2sub([movHeight movWidth],legEndpointPixel);
                proposedLegEndpointsDistancedFromCentralLine(j) = tempGeoDist(legEndpointPixel); %#ok<*AGROW>
            end
        
            % Choose targetNumLegs from proposed endpoints by choosing the ones that have maximal
            % distance from central line

            [~,farthestLegIndices] = sort(proposedLegEndpointsDistancedFromCentralLine,'descend');
            farthestLegIndices = farthestLegIndices(1:targetNumLegs);

            legEndpointPixelCoordinates = proposedLegEndpoints(farthestLegIndices,:); 
            legEndpointsStore{i} = legEndpointPixelCoordinates;
            frameFlags(i) = 0; % since we force correct number of choices with this method
            %%%%% END METHOD
        end
    end
    
    if methodFlag == 2 || (methodFlag == 1 && AllConnectedComponents.NumObjects < 6)

        % ALTERNATE METHOD
        frameFlags(i) = 1; % mark as less confident (we will mark it as 2 = Requires Correction if incorrect number of legs is found, below)

        rowCoordsMaskCentralLine = [yCent,centralLineY];
        colCoordsMaskCentralLine = [xCent,centralLineX];
        removalIdxCentralLine = find(centralLineX < 1 | centralLineX > movWidth | centralLineY < 1 | centralLineY > movHeight);
        rowCoordsMaskCentralLine(removalIdxCentralLine+1) = [];
        colCoordsMaskCentralLine(removalIdxCentralLine+1) = [];
        geoDist = bwdistgeodesic(imgForTracking,colCoordsMaskCentralLine,rowCoordsMaskCentralLine,'cityblock');
        removalBodyPixelsDistanceHighToLow = max(geoDist(:)); % initialize removal based on distance starting high (at max distance)


        [numPixelsDistance,~] = hist(geoDist(:),max(geoDist(:)));
        derivativeCumulativeNumPixelsDistance = diff(cumsum(numPixelsDistance)); % derivative of cumulative sum of pixels at each distance from central axis flattens out (becomes say half of initial one), because then a
            % large part of the body pixels are removed and the further changes to number of pixels at further distances are relativey small
        minBodyPixelsDistance = find(derivativeCumulativeNumPixelsDistance < minDistanceThresholdFactor*derivativeCumulativeNumPixelsDistance(1),1,'first');

        numLegsFoundHighToLow = 0; % reset
        while (removalBodyPixelsDistanceHighToLow >= minBodyPixelsDistance)
                % while the right number of legs hasn't been found in this frame AND we haven't reached the point where we include most body pixels (= removalBodyPixelsDistanceLowToHigh)
            removalBodyPixelsDistanceHighToLow = removalBodyPixelsDistanceHighToLow - 1; % decrease in first step and use <= sign below
            imgForTrackingLegsOnlyHighToLow = imgForTracking;
            imgForTrackingLegsOnlyHighToLow(find(geoDist <= removalBodyPixelsDistanceHighToLow)) = 0; %#ok<*FNDSB> % remove body, leaving only legs as connected components
            legsHighToLowConnectedComponents = bwconncomp(imgForTrackingLegsOnlyHighToLow,8);
            numLegsFoundHighToLow = legsHighToLowConnectedComponents.NumObjects; % number of connected components going high to low

            % if we find all the legs or more points are tagged
            if numLegsFoundHighToLow == targetNumLegs % if equal, break out of the loop
                break
            elseif numLegsFoundHighToLow > targetNumLegs % if greater
                removalBodyPixelsDistanceHighToLow = removalBodyPixelsDistanceHighToLow + 1; % bring removalBodyPixelsDistanceHighToLow back up by one and redo procedure
                imgForTrackingLegsOnlyHighToLow = imgForTracking;
                imgForTrackingLegsOnlyHighToLow(find(geoDist <= removalBodyPixelsDistanceHighToLow)) = 0; 
                legsHighToLowConnectedComponents = bwconncomp(imgForTrackingLegsOnlyHighToLow,8);
                numLegsFoundHighToLow = legsHighToLowConnectedComponents.NumObjects;
                break
            end
        end
        removalBodyPixelsDistanceHighToLowVector(i) = removalBodyPixelsDistanceHighToLow;

        legEndpointPixelCoordinates = zeros(numLegsFoundHighToLow,2);
        for n = 1:numLegsFoundHighToLow
            pixelIdxLegs = legsHighToLowConnectedComponents.PixelIdxList{n};
            [~, descendingSortLegDistanceIdx] = sort(geoDist(pixelIdxLegs),'descend');
            legEndpointPixel = pixelIdxLegs(descendingSortLegDistanceIdx(1));
            [legEndpointPixelCoordinates(n,2),legEndpointPixelCoordinates(n,1)] = ind2sub(size(imgForTracking),legEndpointPixel);
        end

        legEndpointsStore{i} = legEndpointPixelCoordinates;

        if numLegsFoundHighToLow ~= targetNumLegs
            frameFlags(i) = 2; % if while loop failed to find the right number of legs, flag the frame as REQUIRING intervention
        end
        % END METHOD
    end

    
    if isequal(showMovieFlag,'y')
        subplot(2,1,1)
        imshow(img)
        hold on
        scatter(legEndpointPixelCoordinates(:,1),legEndpointPixelCoordinates(:,2),'r')
        scatter(xCent,yCent,'b')
        scatter(headDirectionCoordinates(i,1),headDirectionCoordinates(i,2),'g')
        scatter(tailDirectionCoordinates(i,1),tailDirectionCoordinates(i,2),'b')
        plot(centralLineX,centralLineY,'g');
        
        titleString = ['Frame Number: ',num2str(i)];
        title(titleString)
        hold off
        subplot(2,1,2)
        imshow(imgForTracking)
        pause(moviePauseTime)
    end
    
end


% Save legEndpointsStore, frameFlags, centroid coordinates, and head coordinates (movie not saved here)
legEndpointsSaveString = strrep(movName, '.avi', '_legEndpointsStore.mat');
frameFlagsSaveString = strrep(movName,'.avi','_frameFlags.mat');
centroidCoordinatesSaveString = strrep(movName,'.avi','_centroidCoordinates.mat');
headDirectionCoordinatesSaveString = strrep(movName,'.avi','_headDirectionCoordinates.mat');

save(legEndpointsSaveString,'legEndpointsStore')
save(frameFlagsSaveString,'frameFlags')
save(centroidCoordinatesSaveString,'centroidCoordinates')
save(headDirectionCoordinatesSaveString,'headDirectionCoordinates')


end