function [DTItracts] = ReadDiffusivity_Deterministic(DTItracts,path,nan_rows)

% ----------------- INPUT -----------------
% DTItracts: structure created in "MuscleParameters" function
% path: string of pathway directory of where the diffusivity files are located (EX: 'FCU')
% nan_rows: logical array for rows (samples) that were eliminated, created in "FindCoordinates_Deterministic" function

% ----------------- OUTPUT -----------------
% DTItracts: strcuture with the fields
%     - FA: Mx1 vector of average fractional anisotropy for each m-th fiber
%     - AD: Mx1 vector of average axial diffusivity for each m-th fiber
%     - RD: Mx1 vector of average radial diffusivity for each m-th fiber
%     - MD: Mx1 vector of average mean diffusivity for each m-th fiber

%% Load and Setup

% Load Files
DSI_coord = strcat(path,'/DSI_studio_tracts.mat');      % Specify filename
DSI_tracts = load(DSI_coord);                           % Load the DSI studio tracts file
metrics_index = DTItracts.fibindex_dti;                 % Save index of tract startpoint and endpoints

% Initialize Containers
M = length(DSI_tracts.length);                  % M = # of fibers tracked
diff_metrics = zeros(M,4);                      % Empty container for average values of each tract: each row = tract, each column = diffusivity metric (col 1 = FA, col 2 = MD, col 3 = AD, col 4 = RD)
diff_values = cell(M,4);                        % Empty container for all values for each fiber: each row = tract, each column = diffusivity metric (col 1 = FA, col 2 = MD, col 3 = AD, col 4 = RD)

% Specify filename paths
FAname = strcat(path,'/FA.txt');
MDname = strcat(path,'/MD.txt');
ADname = strcat(path,'/AD.txt');
RDname = strcat(path,'/RD.txt');

% Create string vector with filename paths
filenames = {FAname;MDname;ADname;RDname};

%% Find the Diffusivity Metrics

% Loop through each diffusivity metric
for i = 1:4     % For each of the 4 diffusivity metrics
    metrics = fscanf( fopen(filenames{i},'r') , '%f' );         % Load file for that diffusivity metric as a vector

    % Loop through each tract
    for m = 1:M
        diff_values{m,i} = metrics(metrics_index(m,1):metrics_index(m,2));      % Create Nx1 array with all the diffusivity values in that tract
        diff_metrics(m,i) = mean(diff_values{m,i});                             % Find the average of the diffusivity metric across all points in that tract     
    end
end
    
%% Eliminate samples that did not track

if ~isempty(nan_rows)                   % If there are any deleted fibers
    diff_metrics(nan_rows,:) = [];      % Eliminate
end

%% Save diffusivity values to the structure

DTItracts.FA = diff_metrics(:,1);
DTItracts.MD = diff_metrics(:,2);
DTItracts.AD = diff_metrics(:,3);
DTItracts.RD = diff_metrics(:,4);

end