function [DTItracts] = FindCurvature(DTItracts,M)

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

% ----------------- OUTPUT -----------------
% DTItracts: structure with the fields 
%   - curvature: the curvature for M tracts (Mx1 vector)

%% Find curvature for all M tracts

DTItracts.curvature = zeros(M,1);   % Create container for pennation angle of each of M tracts

% Loop through 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 curvature evaluated at all t points
    kappa=zeros(1,length(t_int));           % Nx1 vector with curvature evaluated at all points of the tract
    
    % 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''

        % Find curvature, K
        kappa(i) = norm(cross(rdot,rdot2)) / norm(rdot)^3;                      % curvature = ( ||r'(t) x r''(t)|| ) / ||r'(t)||^3

    end

    % Take the average of curvature values across all t values to find the curvature, save to the structure
    DTItracts.curvature(m) = mean(kappa);

end

end