function [DTItracts] = FindInflection(DTItracts,M)

% ----------------- INPUT -----------------
% DTItracts: structure created in "LSQ_Polynomial" function
% M: the generated number of tracts

% ----------------- OUTPUT -----------------
% DTItracts: structure with the fields
%     - inflection_loc: Mx1 cell with t location of inflection points of each m-th tracts
%     - inflection_pts: Mx1 cell with number of inflection points of each of m-th tracts

%% Create empty containers

DTItracts.inflection_loc = cell(M,1);   % Create container for the locations of inflection points of each of M tracts
DTItracts.inflection_pts = zeros(M,1);  % Create container for number of inflection points of each of M fibers

% Loop through each tracts to find and count inflection points for each tract
for m = 1:M     % For each m-th tract

    % Load polynomial coefficients and t values for the m-th fiber --> x(t), y(t), z(t)
    coeff.x = DTItracts.PolyCoeff(m).x;
    coeff.y = DTItracts.PolyCoeff(m).y;
    coeff.z = DTItracts.PolyCoeff(m).z;
    coeff.t0 = DTItracts.PolyCoeff(m).t0;
    coeff.t1 = DTItracts.PolyCoeff(m).t1;

    % Sample from t with N values
    N = 1000;                               % Can be changed
    t_int = linspace(coeff.t0,coeff.t1,N);  % Create vector with N points from t0 to t1 (1xN vector)

    % Find the first derivatives for each polynomial: x'(t), y'(t), z'(t)
    dxdt = polyder(coeff.x);
    dydt = polyder(coeff.y);
    dzdt = polyder(coeff.z);

    % Find the first derivatives for each polynomial: x'(t), y'(t), z'(t)
    dx2dt2 = polyder(dxdt);
    dy2dt2 = polyder(dydt);
    dz2dt2 = polyder(dzdt);
    
    % Create container to save the inflection points evaluated at all t points
    inflection=zeros(1,length(t_int));      % 1xN vector of inflection points for each point t

    % Loop through each t-value
    for i = 1:length(t_int)     % For each of the t value

        tc = t_int(i);          % Find that t-value

        % Evaluate first derivative, r'(t)
        rdot = [polyval(dxdt,tc); polyval(dydt,tc); polyval(dzdt,tc)];          % 3x1 vector, each row is x',y',z'

        % Evaluate second derivative, r''(t)
        rdot2 = [polyval(dx2dt2,tc);polyval(dy2dt2,tc);polyval(dz2dt2,tc)];     % 3x1 vector, each row is x'',y'',z''
                
        % Determine whether r' and r'' are parallel (r' . r'' = |r'|*|r''|*cos(0 or pi))
        dot_prod = round(dot(rdot,rdot2), 2, 'significant');                    % Find r' . r''
        norm_prod = round(norm(rdot) * norm(rdot2), 2, 'significant');          % Find |r'|*|r''|
        inflection(i) = isequal(abs(dot_prod),abs(norm_prod));                  % If equal, there is an inflection point

    end
    
    % Save inflection points to the structure
    DTItracts.inflection_loc{m} = find(inflection);         % Find the t locations of each tract at which there is an inflection
    DTItracts.inflection_pts(m) = numel(find(inflection));  % Count the number of inflection points per tract

end

end